diff --git a/Modules/ImageExtraction/mitkExtractDirectedPlaneImageFilter.cpp b/Modules/ImageExtraction/mitkExtractDirectedPlaneImageFilter.cpp index 7a910aec58..5829c201cd 100644 --- a/Modules/ImageExtraction/mitkExtractDirectedPlaneImageFilter.cpp +++ b/Modules/ImageExtraction/mitkExtractDirectedPlaneImageFilter.cpp @@ -1,509 +1,511 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "mitkExtractDirectedPlaneImageFilter.h" #include "mitkAbstractTransformGeometry.h" //#include "mitkImageMapperGL2D.h" #include #include #include #include #include "vtkMitkThickSlicesFilter.h" #include #include #include #include #include #include #include #include "pic2vtk.h" mitk::ExtractDirectedPlaneImageFilter::ExtractDirectedPlaneImageFilter() : m_WorldGeometry(NULL) { + MITK_WARN << "Class ExtractDirectedPlaneImageFilter is deprecated! Use ExtractSliceFilter instead."; + m_Reslicer = vtkImageReslice::New(); m_TargetTimestep = 0; m_InPlaneResampleExtentByGeometry = true; m_ResliceInterpolationProperty = NULL;//VtkResliceInterpolationProperty::New(); //TODO initial with value m_ThickSlicesMode = 0; m_ThickSlicesNum = 1; } mitk::ExtractDirectedPlaneImageFilter::~ExtractDirectedPlaneImageFilter() { if(m_ResliceInterpolationProperty!=NULL)m_ResliceInterpolationProperty->Delete(); m_Reslicer->Delete(); } void mitk::ExtractDirectedPlaneImageFilter::GenerateData() { // A world geometry must be set... if ( m_WorldGeometry == NULL ) { itkWarningMacro(<<"No world geometry has been set. Returning."); return; } Image *input = const_cast< ImageToImageFilter::InputImageType* >( this->GetInput() ); input->Update(); if ( input == NULL ) { itkWarningMacro(<<"No input set."); return; } const TimeSlicedGeometry *inputTimeGeometry = input->GetTimeSlicedGeometry(); if ( ( inputTimeGeometry == NULL ) || ( inputTimeGeometry->GetTimeSteps() == 0 ) ) { itkWarningMacro(<<"Error reading input image geometry."); return; } // Get the target timestep; if none is set, use the lowest given. unsigned int timestep = 0; if ( ! m_TargetTimestep ) { ScalarType time = m_WorldGeometry->GetTimeBounds()[0]; if ( time > ScalarTypeNumericTraits::NonpositiveMin() ) { timestep = inputTimeGeometry->MSToTimeStep( time ); } } else timestep = m_TargetTimestep; if ( inputTimeGeometry->IsValidTime( timestep ) == false ) { itkWarningMacro(<<"This is not a valid timestep: "<IsVolumeSet( timestep ) ) { itkWarningMacro(<<"No volume data existent at given timestep "<GetLargestPossibleRegion(); requestedRegion.SetIndex( 3, timestep ); requestedRegion.SetSize( 3, 1 ); requestedRegion.SetSize( 4, 1 ); input->SetRequestedRegion( &requestedRegion ); input->Update(); vtkImageData* inputData = input->GetVtkImageData( timestep ); if ( inputData == NULL ) { itkWarningMacro(<<"Could not extract vtk image data for given timestep"<GetSpacing( spacing ); // how big the area is in physical coordinates: widthInMM x heightInMM pixels mitk::ScalarType widthInMM, heightInMM; // where we want to sample Point3D origin; Vector3D right, bottom, normal; Vector3D rightInIndex, bottomInIndex; assert( input->GetTimeSlicedGeometry() == inputTimeGeometry ); // take transform of input image into account Geometry3D* inputGeometry = inputTimeGeometry->GetGeometry3D( timestep ); if ( inputGeometry == NULL ) { itkWarningMacro(<<"There is no Geometry3D at given timestep "<( m_WorldGeometry ) != NULL ) { const PlaneGeometry *planeGeometry = static_cast< const PlaneGeometry * >( m_WorldGeometry ); origin = planeGeometry->GetOrigin(); right = planeGeometry->GetAxisVector( 0 ); bottom = planeGeometry->GetAxisVector( 1 ); normal = planeGeometry->GetNormal(); if ( m_InPlaneResampleExtentByGeometry ) { // Resampling grid corresponds to the current world geometry. This // means that the spacing of the output 2D image depends on the // currently selected world geometry, and *not* on the image itself. extent[0] = m_WorldGeometry->GetExtent( 0 ); extent[1] = m_WorldGeometry->GetExtent( 1 ); } else { // Resampling grid corresponds to the input geometry. This means that // the spacing of the output 2D image is directly derived from the // associated input image, regardless of the currently selected world // geometry. inputGeometry->WorldToIndex( right, rightInIndex ); inputGeometry->WorldToIndex( bottom, bottomInIndex ); extent[0] = rightInIndex.GetNorm(); extent[1] = bottomInIndex.GetNorm(); } // Get the extent of the current world geometry and calculate resampling // spacing therefrom. widthInMM = m_WorldGeometry->GetExtentInMM( 0 ); heightInMM = m_WorldGeometry->GetExtentInMM( 1 ); mmPerPixel[0] = widthInMM / extent[0]; mmPerPixel[1] = heightInMM / extent[1]; right.Normalize(); bottom.Normalize(); normal.Normalize(); //origin += right * ( mmPerPixel[0] * 0.5 ); //origin += bottom * ( mmPerPixel[1] * 0.5 ); //widthInMM -= mmPerPixel[0]; //heightInMM -= mmPerPixel[1]; // Use inverse transform of the input geometry for reslicing the 3D image m_Reslicer->SetResliceTransform( inputGeometry->GetVtkTransform()->GetLinearInverse() ); // Set background level to TRANSLUCENT (see Geometry2DDataVtkMapper3D) m_Reslicer->SetBackgroundLevel( -32768 ); // Check if a reference geometry does exist (as would usually be the case for // PlaneGeometry). // Note: this is currently not strictly required, but could facilitate // correct plane clipping. if ( m_WorldGeometry->GetReferenceGeometry() ) { // Calculate the actual bounds of the transformed plane clipped by the // dataset bounding box; this is required for drawing the texture at the // correct position during 3D mapping. boundsInitialized = this->CalculateClippedPlaneBounds( m_WorldGeometry->GetReferenceGeometry(), planeGeometry, bounds ); } } // Do we have an AbstractTransformGeometry? else if ( dynamic_cast< const AbstractTransformGeometry * >( m_WorldGeometry ) ) { const mitk::AbstractTransformGeometry* abstractGeometry = dynamic_cast< const AbstractTransformGeometry * >(m_WorldGeometry); extent[0] = abstractGeometry->GetParametricExtent(0); extent[1] = abstractGeometry->GetParametricExtent(1); widthInMM = abstractGeometry->GetParametricExtentInMM(0); heightInMM = abstractGeometry->GetParametricExtentInMM(1); mmPerPixel[0] = widthInMM / extent[0]; mmPerPixel[1] = heightInMM / extent[1]; origin = abstractGeometry->GetPlane()->GetOrigin(); right = abstractGeometry->GetPlane()->GetAxisVector(0); right.Normalize(); bottom = abstractGeometry->GetPlane()->GetAxisVector(1); bottom.Normalize(); normal = abstractGeometry->GetPlane()->GetNormal(); normal.Normalize(); // Use a combination of the InputGeometry *and* the possible non-rigid // AbstractTransformGeometry for reslicing the 3D Image vtkGeneralTransform *composedResliceTransform = vtkGeneralTransform::New(); composedResliceTransform->Identity(); composedResliceTransform->Concatenate( inputGeometry->GetVtkTransform()->GetLinearInverse() ); composedResliceTransform->Concatenate( abstractGeometry->GetVtkAbstractTransform() ); m_Reslicer->SetResliceTransform( composedResliceTransform ); // Set background level to BLACK instead of translucent, to avoid // boundary artifacts (see Geometry2DDataVtkMapper3D) m_Reslicer->SetBackgroundLevel( -1023 ); composedResliceTransform->Delete(); } else { itkWarningMacro(<<"World Geometry has to be a PlaneGeometry or an AbstractTransformGeometry."); return; } // Make sure that the image to be resliced has a certain minimum size. if ( (extent[0] <= 2) && (extent[1] <= 2) ) { itkWarningMacro(<<"Image is too small to be resliced..."); return; } vtkImageChangeInformation * unitSpacingImageFilter = vtkImageChangeInformation::New() ; unitSpacingImageFilter->SetOutputSpacing( 1.0, 1.0, 1.0 ); unitSpacingImageFilter->SetInput( inputData ); m_Reslicer->SetInput( unitSpacingImageFilter->GetOutput() ); unitSpacingImageFilter->Delete(); //m_Reslicer->SetInput( inputData ); m_Reslicer->SetOutputDimensionality( 2 ); m_Reslicer->SetOutputOrigin( 0.0, 0.0, 0.0 ); Vector2D pixelsPerMM; pixelsPerMM[0] = 1.0 / mmPerPixel[0]; pixelsPerMM[1] = 1.0 / mmPerPixel[1]; //calulate the originArray and the orientations for the reslice-filter double originArray[3]; itk2vtk( origin, originArray ); m_Reslicer->SetResliceAxesOrigin( originArray ); double cosines[9]; // direction of the X-axis of the sampled result vnl2vtk( right.Get_vnl_vector(), cosines ); // direction of the Y-axis of the sampled result vnl2vtk( bottom.Get_vnl_vector(), cosines + 3 ); // normal of the plane vnl2vtk( normal.Get_vnl_vector(), cosines + 6 ); m_Reslicer->SetResliceAxesDirectionCosines( cosines ); int xMin, xMax, yMin, yMax; if ( boundsInitialized ) { xMin = static_cast< int >( bounds[0] / mmPerPixel[0] );//+ 0.5 ); xMax = static_cast< int >( bounds[1] / mmPerPixel[0] );//+ 0.5 ); yMin = static_cast< int >( bounds[2] / mmPerPixel[1] );//+ 0.5); yMax = static_cast< int >( bounds[3] / mmPerPixel[1] );//+ 0.5 ); } else { // If no reference geometry is available, we also don't know about the // maximum plane size; so the overlap is just ignored xMin = yMin = 0; xMax = static_cast< int >( extent[0] - pixelsPerMM[0] );//+ 0.5 ); yMax = static_cast< int >( extent[1] - pixelsPerMM[1] );//+ 0.5 ); } m_Reslicer->SetOutputSpacing( mmPerPixel[0], mmPerPixel[1], 1.0 ); // xMax and yMax are meant exclusive until now, whereas // SetOutputExtent wants an inclusive bound. Thus, we need // to subtract 1. m_Reslicer->SetOutputExtent( xMin, xMax-1, yMin, yMax-1, 0, 1 ); // Do the reslicing. Modified() is called to make sure that the reslicer is // executed even though the input geometry information did not change; this // is necessary when the input /em data, but not the /em geometry changes. m_Reslicer->Modified(); m_Reslicer->ReleaseDataFlagOn(); m_Reslicer->Update(); // 1. Check the result vtkImageData* reslicedImage = m_Reslicer->GetOutput(); //mitkIpPicDescriptor *pic = Pic2vtk::convert( reslicedImage ); if((reslicedImage == NULL) || (reslicedImage->GetDataDimension() < 1)) { itkWarningMacro(<<"Reslicer returned empty image"); return; } unsigned int dimensions[2]; dimensions[0] = (unsigned int)extent[0]; dimensions[1] = (unsigned int)extent[1]; Vector3D spacingVector; FillVector3D(spacingVector, mmPerPixel[0], mmPerPixel[1], 1.0); mitk::Image::Pointer resultImage = this->GetOutput(); resultImage->Initialize(input->GetPixelType(), 2, dimensions ); //resultImage->Initialize( pic ); resultImage->SetSpacing( spacingVector ); //resultImage->SetPicVolume( pic ); //mitkIpPicFree(pic); /*unsigned int dimensions[2]; dimensions[0] = (unsigned int)extent[0]; dimensions[1] = (unsigned int)extent[1]; Vector3D spacingVector; FillVector3D(spacingVector, mmPerPixel[0], mmPerPixel[1], 1.0); mitk::Image::Pointer resultImage = this->GetOutput(); resultImage->Initialize(m_Reslicer->GetOutput()); resultImage->Initialize(inputImage->GetPixelType(), 2, dimensions); resultImage->SetSpacing(spacingVector); resultImage->SetSlice(m_Reslicer->GetOutput());*/ } void mitk::ExtractDirectedPlaneImageFilter::GenerateOutputInformation() { Superclass::GenerateOutputInformation(); } bool mitk::ExtractDirectedPlaneImageFilter ::CalculateClippedPlaneBounds( const Geometry3D *boundingGeometry, const PlaneGeometry *planeGeometry, vtkFloatingPointType *bounds ) { // Clip the plane with the bounding geometry. To do so, the corner points // of the bounding box are transformed by the inverse transformation // matrix, and the transformed bounding box edges derived therefrom are // clipped with the plane z=0. The resulting min/max values are taken as // bounds for the image reslicer. const BoundingBox *boundingBox = boundingGeometry->GetBoundingBox(); BoundingBox::PointType bbMin = boundingBox->GetMinimum(); BoundingBox::PointType bbMax = boundingBox->GetMaximum(); vtkPoints *points = vtkPoints::New(); if(boundingGeometry->GetImageGeometry()) { points->InsertPoint( 0, bbMin[0]-0.5, bbMin[1]-0.5, bbMin[2]-0.5 ); points->InsertPoint( 1, bbMin[0]-0.5, bbMin[1]-0.5, bbMax[2]-0.5 ); points->InsertPoint( 2, bbMin[0]-0.5, bbMax[1]-0.5, bbMax[2]-0.5 ); points->InsertPoint( 3, bbMin[0]-0.5, bbMax[1]-0.5, bbMin[2]-0.5 ); points->InsertPoint( 4, bbMax[0]-0.5, bbMin[1]-0.5, bbMin[2]-0.5 ); points->InsertPoint( 5, bbMax[0]-0.5, bbMin[1]-0.5, bbMax[2]-0.5 ); points->InsertPoint( 6, bbMax[0]-0.5, bbMax[1]-0.5, bbMax[2]-0.5 ); points->InsertPoint( 7, bbMax[0]-0.5, bbMax[1]-0.5, bbMin[2]-0.5 ); } else { points->InsertPoint( 0, bbMin[0], bbMin[1], bbMin[2] ); points->InsertPoint( 1, bbMin[0], bbMin[1], bbMax[2] ); points->InsertPoint( 2, bbMin[0], bbMax[1], bbMax[2] ); points->InsertPoint( 3, bbMin[0], bbMax[1], bbMin[2] ); points->InsertPoint( 4, bbMax[0], bbMin[1], bbMin[2] ); points->InsertPoint( 5, bbMax[0], bbMin[1], bbMax[2] ); points->InsertPoint( 6, bbMax[0], bbMax[1], bbMax[2] ); points->InsertPoint( 7, bbMax[0], bbMax[1], bbMin[2] ); } vtkPoints *newPoints = vtkPoints::New(); vtkTransform *transform = vtkTransform::New(); transform->Identity(); transform->Concatenate( planeGeometry->GetVtkTransform()->GetLinearInverse() ); transform->Concatenate( boundingGeometry->GetVtkTransform() ); transform->TransformPoints( points, newPoints ); transform->Delete(); bounds[0] = bounds[2] = 10000000.0; bounds[1] = bounds[3] = -10000000.0; bounds[4] = bounds[5] = 0.0; this->LineIntersectZero( newPoints, 0, 1, bounds ); this->LineIntersectZero( newPoints, 1, 2, bounds ); this->LineIntersectZero( newPoints, 2, 3, bounds ); this->LineIntersectZero( newPoints, 3, 0, bounds ); this->LineIntersectZero( newPoints, 0, 4, bounds ); this->LineIntersectZero( newPoints, 1, 5, bounds ); this->LineIntersectZero( newPoints, 2, 6, bounds ); this->LineIntersectZero( newPoints, 3, 7, bounds ); this->LineIntersectZero( newPoints, 4, 5, bounds ); this->LineIntersectZero( newPoints, 5, 6, bounds ); this->LineIntersectZero( newPoints, 6, 7, bounds ); this->LineIntersectZero( newPoints, 7, 4, bounds ); // clean up vtk data points->Delete(); newPoints->Delete(); if ( (bounds[0] > 9999999.0) || (bounds[2] > 9999999.0) || (bounds[1] < -9999999.0) || (bounds[3] < -9999999.0) ) { return false; } else { // The resulting bounds must be adjusted by the plane spacing, since we // we have so far dealt with index coordinates const float *planeSpacing = planeGeometry->GetFloatSpacing(); bounds[0] *= planeSpacing[0]; bounds[1] *= planeSpacing[0]; bounds[2] *= planeSpacing[1]; bounds[3] *= planeSpacing[1]; bounds[4] *= planeSpacing[2]; bounds[5] *= planeSpacing[2]; return true; } } bool mitk::ExtractDirectedPlaneImageFilter ::LineIntersectZero( vtkPoints *points, int p1, int p2, vtkFloatingPointType *bounds ) { vtkFloatingPointType point1[3]; vtkFloatingPointType point2[3]; points->GetPoint( p1, point1 ); points->GetPoint( p2, point2 ); if ( (point1[2] * point2[2] <= 0.0) && (point1[2] != point2[2]) ) { double x, y; x = ( point1[0] * point2[2] - point1[2] * point2[0] ) / ( point2[2] - point1[2] ); y = ( point1[1] * point2[2] - point1[2] * point2[1] ) / ( point2[2] - point1[2] ); if ( x < bounds[0] ) { bounds[0] = x; } if ( x > bounds[1] ) { bounds[1] = x; } if ( y < bounds[2] ) { bounds[2] = y; } if ( y > bounds[3] ) { bounds[3] = y; } bounds[4] = bounds[5] = 0.0; return true; } return false; } diff --git a/Modules/ImageExtraction/mitkExtractDirectedPlaneImageFilter.h b/Modules/ImageExtraction/mitkExtractDirectedPlaneImageFilter.h index 8161b01819..5afd23ced0 100644 --- a/Modules/ImageExtraction/mitkExtractDirectedPlaneImageFilter.h +++ b/Modules/ImageExtraction/mitkExtractDirectedPlaneImageFilter.h @@ -1,121 +1,124 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #ifndef mitkExtractDirectedPlaneImageFilter_h_Included #define mitkExtractDirectedPlaneImageFilter_h_Included #include "ImageExtractionExports.h" #include "mitkImageToImageFilter.h" #include "vtkImageReslice.h" #include "mitkVtkResliceInterpolationProperty.h" #define setMacro(name,type) \ virtual void Set##name (type _arg) \ { \ if (this->m_##name != _arg) \ { \ this->m_##name = _arg; \ } \ } #define getMacro(name,type) \ virtual type Get##name () \ { \ return m_##name; \ } class vtkPoints; namespace mitk { - /** +/** + \deprecated This class is deprecated. Use mitk::ExtractSliceFilter instead. + \sa ExtractSliceFilter + \brief Extracts a 2D slice of arbitrary geometry from a 3D or 4D image. \sa mitkImageMapper2D \ingroup ImageToImageFilter This class takes a 3D or 4D mitk::Image as input and tries to extract one slice from it. This slice can be arbitrary oriented in space. The 2D slice is resliced by a vtk::ResliceImage filter if not perpendicular to the input image. The world geometry of the plane to be extracted image must be given as an input to the filter in order to correctly calculate world coordinates of the extracted slice. Setting a timestep from which the plane should be extracted is optional. Output will not be set if there was a problem extracting the desired slice. Last contributor: $Author: T. Schwarz$ - */ +*/ class ImageExtraction_EXPORT ExtractDirectedPlaneImageFilter : public ImageToImageFilter { public: mitkClassMacro(ExtractDirectedPlaneImageFilter, ImageToImageFilter); itkNewMacro(ExtractDirectedPlaneImageFilter); itkSetMacro( WorldGeometry, Geometry2D* ); // The Reslicer is accessible to configure the desired interpolation; // (See vtk::ImageReslice class for documentation). // Misusage is at your own risk... itkGetMacro( Reslicer, vtkImageReslice* ); // The target timestep in a 4D image from which the 2D plane is supposed // to be extracted. itkSetMacro( TargetTimestep, unsigned int ); itkGetMacro( TargetTimestep, unsigned int ); itkSetMacro( InPlaneResampleExtentByGeometry, bool ); itkGetMacro( InPlaneResampleExtentByGeometry, bool ); setMacro( ResliceInterpolationProperty, VtkResliceInterpolationProperty* ); itkGetMacro( ResliceInterpolationProperty, VtkResliceInterpolationProperty* ); setMacro( IsMapperMode, bool ); getMacro( IsMapperMode, bool ); protected: ExtractDirectedPlaneImageFilter(); // purposely hidden virtual ~ExtractDirectedPlaneImageFilter(); virtual void GenerateData(); virtual void GenerateOutputInformation(); bool CalculateClippedPlaneBounds( const Geometry3D *boundingGeometry, const PlaneGeometry *planeGeometry, vtkFloatingPointType *bounds ); bool LineIntersectZero( vtkPoints *points, int p1, int p2, vtkFloatingPointType *bounds ); const Geometry2D* m_WorldGeometry; vtkImageReslice * m_Reslicer; unsigned int m_TargetTimestep; bool m_InPlaneResampleExtentByGeometry; int m_ThickSlicesMode; int m_ThickSlicesNum; bool m_IsMapperMode; VtkResliceInterpolationProperty* m_ResliceInterpolationProperty; }; } // namespace mitk #endif // mitkExtractDirectedPlaneImageFilter_h_Included diff --git a/Modules/ImageExtraction/mitkExtractDirectedPlaneImageFilterNew.cpp b/Modules/ImageExtraction/mitkExtractDirectedPlaneImageFilterNew.cpp index ec6ce80698..d58c9cd074 100644 --- a/Modules/ImageExtraction/mitkExtractDirectedPlaneImageFilterNew.cpp +++ b/Modules/ImageExtraction/mitkExtractDirectedPlaneImageFilterNew.cpp @@ -1,296 +1,297 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "mitkExtractDirectedPlaneImageFilterNew.h" #include "mitkImageCast.h" #include "mitkImageTimeSelector.h" #include "itkImageRegionIterator.h" #include mitk::ExtractDirectedPlaneImageFilterNew::ExtractDirectedPlaneImageFilterNew() :m_CurrentWorldGeometry2D(NULL), m_ActualInputTimestep(-1) { + MITK_WARN << "Class ExtractDirectedPlaneImageFilterNew is deprecated! Use ExtractSliceFilter instead."; } mitk::ExtractDirectedPlaneImageFilterNew::~ExtractDirectedPlaneImageFilterNew() { } void mitk::ExtractDirectedPlaneImageFilterNew::GenerateData(){ mitk::Image::ConstPointer inputImage = ImageToImageFilter::GetInput(0); if ( !inputImage ) { MITK_ERROR << "mitk::ExtractDirectedPlaneImageFilterNew: No input available. Please set the input!" << std::endl; itkExceptionMacro("mitk::ExtractDirectedPlaneImageFilterNew: No input available. Please set the input!"); return; } m_ImageGeometry = inputImage->GetGeometry(); //If no timestep is set, the lowest given will be selected const mitk::TimeSlicedGeometry* inputTimeGeometry = this->GetInput()->GetTimeSlicedGeometry(); if ( m_ActualInputTimestep == -1) { ScalarType time = m_CurrentWorldGeometry2D->GetTimeBounds()[0]; if ( time > ScalarTypeNumericTraits::NonpositiveMin() ) { m_ActualInputTimestep = inputTimeGeometry->MSToTimeStep( time ); } } if ( inputImage->GetDimension() > 4 || inputImage->GetDimension() < 2) { MITK_ERROR << "mitk::ExtractDirectedPlaneImageFilterNew:GenerateData works only with 3D and 3D+t images, sorry." << std::endl; itkExceptionMacro("mitk::ExtractDirectedPlaneImageFilterNew works only with 3D and 3D+t images, sorry."); return; } else if ( inputImage->GetDimension() == 4 ) { mitk::ImageTimeSelector::Pointer timeselector = mitk::ImageTimeSelector::New(); timeselector->SetInput( inputImage ); timeselector->SetTimeNr( m_ActualInputTimestep ); timeselector->UpdateLargestPossibleRegion(); inputImage = timeselector->GetOutput(); } else if ( inputImage->GetDimension() == 2) { mitk::Image::Pointer resultImage = ImageToImageFilter::GetOutput(); resultImage = const_cast( inputImage.GetPointer() ); ImageToImageFilter::SetNthOutput( 0, resultImage); return; } if ( !m_CurrentWorldGeometry2D ) { MITK_ERROR<< "mitk::ExtractDirectedPlaneImageFilterNew::GenerateData has no CurrentWorldGeometry2D set" << std::endl; return; } AccessFixedDimensionByItk( inputImage, ItkSliceExtraction, 3 ); }//Generate Data void mitk::ExtractDirectedPlaneImageFilterNew::GenerateOutputInformation () { Superclass::GenerateOutputInformation(); } /* * The desired slice is extracted by filling the image`s corresponding pixel values in an empty 2 dimensional itk::Image * Therefor the itk image`s extent in pixel (in each direction) is doubled and its spacing (also in each direction) is divided by two * (similar to the shannon theorem). */ template void mitk::ExtractDirectedPlaneImageFilterNew::ItkSliceExtraction (itk::Image* inputImage) { typedef itk::Image InputImageType; typedef itk::Image SliceImageType; typedef itk::ImageRegionConstIterator< SliceImageType > SliceIterator; //Creating an itk::Image that represents the sampled slice typename SliceImageType::Pointer resultSlice = SliceImageType::New(); typename SliceImageType::IndexType start; start[0] = 0; start[1] = 0; Point3D origin = m_CurrentWorldGeometry2D->GetOrigin(); Vector3D right = m_CurrentWorldGeometry2D->GetAxisVector(0); Vector3D bottom = m_CurrentWorldGeometry2D->GetAxisVector(1); //Calculation the sample-spacing, i.e the half of the smallest spacing existing in the original image Vector3D newPixelSpacing = m_ImageGeometry->GetSpacing(); float minSpacing = newPixelSpacing[0]; for (unsigned int i = 1; i < newPixelSpacing.Size(); i++) { if (newPixelSpacing[i] < minSpacing ) { minSpacing = newPixelSpacing[i]; } } newPixelSpacing[0] = 0.5*minSpacing; newPixelSpacing[1] = 0.5*minSpacing; newPixelSpacing[2] = 0.5*minSpacing; float pixelSpacing[2]; pixelSpacing[0] = newPixelSpacing[0]; pixelSpacing[1] = newPixelSpacing[1]; //Calculating the size of the sampled slice typename SliceImageType::SizeType size; Vector2D extentInMM; extentInMM[0] = m_CurrentWorldGeometry2D->GetExtentInMM(0); extentInMM[1] = m_CurrentWorldGeometry2D->GetExtentInMM(1); //The maximum extent is the lenght of the diagonal of the considered plane double maxExtent = sqrt(extentInMM[0]*extentInMM[0]+extentInMM[1]*extentInMM[1]); unsigned int xTranlation = (maxExtent-extentInMM[0]); unsigned int yTranlation = (maxExtent-extentInMM[1]); size[0] = (maxExtent+xTranlation)/newPixelSpacing[0]; size[1] = (maxExtent+yTranlation)/newPixelSpacing[1]; //Creating an ImageRegion Object typename SliceImageType::RegionType region; region.SetSize( size ); region.SetIndex( start ); //Defining the image`s extent and origin by passing the region to it and allocating memory for it resultSlice->SetRegions( region ); resultSlice->SetSpacing( pixelSpacing ); resultSlice->Allocate(); /* * Here we create an new geometry so that the transformations are calculated correctly (our resulting slice has a different bounding box and spacing) * The original current worldgeometry must be cloned because we have to keep the directions of the axis vector which represents the rotation */ right.Normalize(); bottom.Normalize(); //Here we translate the origin to adapt the new geometry to the previous calculated extent origin[0] -= xTranlation*right[0]+yTranlation*bottom[0]; origin[1] -= xTranlation*right[1]+yTranlation*bottom[1]; origin[2] -= xTranlation*right[2]+yTranlation*bottom[2]; //Putting it together for the new geometry mitk::Geometry3D::Pointer newSliceGeometryTest = dynamic_cast(m_CurrentWorldGeometry2D->Clone().GetPointer()); newSliceGeometryTest->ChangeImageGeometryConsideringOriginOffset(true); //Workaround because of BUG (#6505) newSliceGeometryTest->GetIndexToWorldTransform()->SetMatrix(m_CurrentWorldGeometry2D->GetIndexToWorldTransform()->GetMatrix()); //Workaround end newSliceGeometryTest->SetOrigin(origin); ScalarType bounds[6]={0, size[0], 0, size[1], 0, 1}; newSliceGeometryTest->SetBounds(bounds); newSliceGeometryTest->SetSpacing(newPixelSpacing); newSliceGeometryTest->Modified(); //Workaround because of BUG (#6505) itk::MatrixOffsetTransformBase::MatrixType tempTransform = newSliceGeometryTest->GetIndexToWorldTransform()->GetMatrix(); //Workaround end /* * Now we iterate over the recently created slice. * For each slice - pixel we check whether there is an according * pixel in the input - image which can be set in the slice. * In this way a slice is sampled out of the input - image regrading to the given PlaneGeometry */ Point3D currentSliceIndexPointIn2D; Point3D currentImageWorldPointIn3D; typename InputImageType::IndexType inputIndex; SliceIterator sliceIterator ( resultSlice, resultSlice->GetLargestPossibleRegion() ); sliceIterator.GoToBegin(); while ( !sliceIterator.IsAtEnd() ) { /* * Here we add 0.5 to to assure that the indices are correctly transformed. * (Because of the 0.5er Bug) */ currentSliceIndexPointIn2D[0] = sliceIterator.GetIndex()[0]+0.5; currentSliceIndexPointIn2D[1] = sliceIterator.GetIndex()[1]+0.5; currentSliceIndexPointIn2D[2] = 0; newSliceGeometryTest->IndexToWorld( currentSliceIndexPointIn2D, currentImageWorldPointIn3D ); m_ImageGeometry->WorldToIndex( currentImageWorldPointIn3D, inputIndex); if ( m_ImageGeometry->IsIndexInside( inputIndex )) { resultSlice->SetPixel( sliceIterator.GetIndex(), inputImage->GetPixel(inputIndex) ); } else { resultSlice->SetPixel( sliceIterator.GetIndex(), 0); } ++sliceIterator; } Image::Pointer resultImage = ImageToImageFilter::GetOutput(); GrabItkImageMemory(resultSlice, resultImage, NULL, false); resultImage->SetClonedGeometry(newSliceGeometryTest); //Workaround because of BUG (#6505) resultImage->GetGeometry()->GetIndexToWorldTransform()->SetMatrix(tempTransform); //Workaround end } ///**TEST** May ba a little bit more efficient but doesn`t already work/ //right.Normalize(); //bottom.Normalize(); //Point3D currentImagePointIn3D = origin /*+ bottom*newPixelSpacing*/; //unsigned int columns ( 0 ); /**ENDE**/ /****TEST***/ //SliceImageType::IndexType index = sliceIterator.GetIndex(); //if ( columns == (extentInPixel[0]) ) //{ //If we are at the end of a row, then we have to go to the beginning of the next row //currentImagePointIn3D = origin; //currentImagePointIn3D += newPixelSpacing[1]*bottom*index[1]; //columns = 0; //m_ImageGeometry->WorldToIndex(currentImagePointIn3D, inputIndex); //} //else //{ //// //if ( columns != 0 ) //{ //currentImagePointIn3D += newPixelSpacing[0]*right; //} //m_ImageGeometry->WorldToIndex(currentImagePointIn3D, inputIndex); //} //if ( m_ImageGeometry->IsIndexInside( inputIndex )) //{ //resultSlice->SetPixel( sliceIterator.GetIndex(), inputImage->GetPixel(inputIndex) ); //} //else if (currentImagePointIn3D == origin) //{ //Point3D temp; //temp[0] = bottom[0]*newPixelSpacing[0]*0.5; //temp[1] = bottom[1]*newPixelSpacing[1]*0.5; //temp[2] = bottom[2]*newPixelSpacing[2]*0.5; //origin[0] += temp[0]; //origin[1] += temp[1]; //origin[2] += temp[2]; //currentImagePointIn3D = origin; //m_ImageGeometry->WorldToIndex(currentImagePointIn3D, inputIndex); //if ( m_ImageGeometry->IsIndexInside( inputIndex )) //{ //resultSlice->SetPixel( sliceIterator.GetIndex(), inputImage->GetPixel(inputIndex) ); //} //} /****TEST ENDE****/ diff --git a/Modules/ImageExtraction/mitkExtractDirectedPlaneImageFilterNew.h b/Modules/ImageExtraction/mitkExtractDirectedPlaneImageFilterNew.h index 67d0ee37ce..507abc28ea 100644 --- a/Modules/ImageExtraction/mitkExtractDirectedPlaneImageFilterNew.h +++ b/Modules/ImageExtraction/mitkExtractDirectedPlaneImageFilterNew.h @@ -1,93 +1,96 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #ifndef mitkExtractDirectedPlaneImageFilterNew_h_Included #define mitkExtractDirectedPlaneImageFilterNew_h_Included #include "ImageExtractionExports.h" #include "mitkImageToImageFilter.h" #include "itkImage.h" #include "mitkITKImageImport.h" namespace mitk { /** + \deprecated This class is deprecated. Use mitk::ExtractSliceFilter instead. + \sa ExtractSliceFilter + \brief A filter that can extract a 2D slice from a 3D or 4D image especially if the image`s axes are rotated \sa ContourTool \sa SegTool2D \sa ExtractImageFilter \sa OverwriteSliceImageFilter \sa OverwriteDirectedPlaneImageFilter \ingroup Process \ingroup Reliver There is a separate page describing the general design of QmitkInteractiveSegmentation: \ref QmitkSegmentationTechnicalPage This class takes an 3D or 4D mitk::Image as input and extracts a slice from it. If you work with a 4D image as input you have to specify the desired timestep at which the slice shall be extracted, otherwise the lowest given timestep is selected by default. The special feature of this filter is, that the planes of the input image can be rotated in any way. To assure a proper extraction you have to set the currentWorldGeometry2D with you can obtain from the BaseRenderer, respectively the positionEvent send by the renderer. The output will not be set if there was a problem with the input image $Author: fetzer $ */ class ImageExtraction_EXPORT ExtractDirectedPlaneImageFilterNew : public ImageToImageFilter { public: mitkClassMacro(ExtractDirectedPlaneImageFilterNew, ImageToImageFilter); itkNewMacro(Self); /** \brief Set macro for the current worldgeometry \a Parameter The current wordgeometry that describes the position (rotation, translation) of the plane (and therefore the slice to be extracted) in our 3D(+t) image */ itkSetMacro(CurrentWorldGeometry2D, Geometry3D* ); itkSetMacro(ImageGeometry, Geometry3D* ); /** \brief Set macro for the current timestep \a Parameter The timestep of the image from which the slice shall be extracted */ itkSetMacro(ActualInputTimestep, int); protected: ExtractDirectedPlaneImageFilterNew(); virtual ~ExtractDirectedPlaneImageFilterNew(); virtual void GenerateData(); virtual void GenerateOutputInformation(); private: const Geometry3D* m_CurrentWorldGeometry2D; const Geometry3D* m_ImageGeometry; int m_ActualInputTimestep; template void ItkSliceExtraction (itk::Image* inputImage); }; }//namespace #endif diff --git a/Modules/ImageExtraction/mitkExtractImageFilter.cpp b/Modules/ImageExtraction/mitkExtractImageFilter.cpp index 474dffb6ac..40cb6bc7ff 100644 --- a/Modules/ImageExtraction/mitkExtractImageFilter.cpp +++ b/Modules/ImageExtraction/mitkExtractImageFilter.cpp @@ -1,244 +1,245 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "mitkExtractImageFilter.h" #include "mitkImageCast.h" #include "mitkPlaneGeometry.h" #include "mitkITKImageImport.h" #include "mitkImageTimeSelector.h" #include #include mitk::ExtractImageFilter::ExtractImageFilter() :m_SliceIndex(0), m_SliceDimension(0), m_TimeStep(0) { + MITK_WARN << "Class ExtractImageFilter is deprecated! Use ExtractSliceFilter instead."; } mitk::ExtractImageFilter::~ExtractImageFilter() { } void mitk::ExtractImageFilter::GenerateData() { Image::ConstPointer input = ImageToImageFilter::GetInput(0); if ( (input->GetDimension() > 4) || (input->GetDimension() < 2) ) { MITK_ERROR << "mitk::ExtractImageFilter:GenerateData works only with 3D and 3D+t images, sorry." << std::endl; itkExceptionMacro("mitk::ExtractImageFilter works only with 3D and 3D+t images, sorry."); return; } else if (input->GetDimension() == 4) { ImageTimeSelector::Pointer timeSelector = ImageTimeSelector::New(); timeSelector->SetInput( input ); timeSelector->SetTimeNr( m_TimeStep ); timeSelector->UpdateLargestPossibleRegion(); input = timeSelector->GetOutput(); } else if (input->GetDimension() == 2) { Image::Pointer resultImage = ImageToImageFilter::GetOutput(); resultImage = const_cast(input.GetPointer()); ImageToImageFilter::SetNthOutput( 0, resultImage ); return; } if ( m_SliceDimension >= input->GetDimension() ) { MITK_ERROR << "mitk::ExtractImageFilter:GenerateData m_SliceDimension == " << m_SliceDimension << " makes no sense with an " << input->GetDimension() << "D image." << std::endl; itkExceptionMacro("This is not a sensible value for m_SliceDimension."); return; } AccessFixedDimensionByItk( input, ItkImageProcessing, 3 ); // set a nice geometry for display and point transformations Geometry3D* inputImageGeometry = ImageToImageFilter::GetInput(0)->GetGeometry(); if (!inputImageGeometry) { MITK_ERROR << "In ExtractImageFilter::ItkImageProcessing: Input image has no geometry!" << std::endl; return; } PlaneGeometry::PlaneOrientation orientation = PlaneGeometry::Transversal; switch ( m_SliceDimension ) { default: case 2: orientation = PlaneGeometry::Transversal; break; case 1: orientation = PlaneGeometry::Frontal; break; case 0: orientation = PlaneGeometry::Sagittal; break; } PlaneGeometry::Pointer planeGeometry = PlaneGeometry::New(); planeGeometry->InitializeStandardPlane( inputImageGeometry, orientation, (ScalarType)m_SliceIndex, true, false ); Image::Pointer resultImage = ImageToImageFilter::GetOutput(); planeGeometry->ChangeImageGeometryConsideringOriginOffset(true); resultImage->SetGeometry( planeGeometry ); } template void mitk::ExtractImageFilter::ItkImageProcessing( itk::Image* itkImage ) { // use the itk::ExtractImageFilter to get a 2D image typedef itk::Image< TPixel, VImageDimension > ImageType3D; typedef itk::Image< TPixel, VImageDimension-1 > ImageType2D; typename ImageType3D::RegionType inSliceRegion = itkImage->GetLargestPossibleRegion(); inSliceRegion.SetSize( m_SliceDimension, 0 ); typedef itk::ExtractImageFilter ExtractImageFilterType; typename ExtractImageFilterType::Pointer sliceExtractor = ExtractImageFilterType::New(); sliceExtractor->SetInput( itkImage ); inSliceRegion.SetIndex( m_SliceDimension, m_SliceIndex ); sliceExtractor->SetExtractionRegion( inSliceRegion ); // calculate the output sliceExtractor->UpdateLargestPossibleRegion(); typename ImageType2D::Pointer slice = sliceExtractor->GetOutput(); // re-import to MITK Image::Pointer resultImage = ImageToImageFilter::GetOutput(); GrabItkImageMemory(slice, resultImage, NULL, false); } /* * What is the input requested region that is required to produce the output * requested region? By default, the largest possible region is always * required but this is overridden in many subclasses. For instance, for an * image processing filter where an output pixel is a simple function of an * input pixel, the input requested region will be set to the output * requested region. For an image processing filter where an output pixel is * a function of the pixels in a neighborhood of an input pixel, then the * input requested region will need to be larger than the output requested * region (to avoid introducing artificial boundary conditions). This * function should never request an input region that is outside the the * input largest possible region (i.e. implementations of this method should * crop the input requested region at the boundaries of the input largest * possible region). */ void mitk::ExtractImageFilter::GenerateInputRequestedRegion() { Superclass::GenerateInputRequestedRegion(); ImageToImageFilter::InputImagePointer input = const_cast< ImageToImageFilter::InputImageType* > ( this->GetInput() ); Image::Pointer output = this->GetOutput(); if (input->GetDimension() == 2) { input->SetRequestedRegionToLargestPossibleRegion(); return; } Image::RegionType requestedRegion; requestedRegion = output->GetRequestedRegion(); requestedRegion.SetIndex(0, 0); requestedRegion.SetIndex(1, 0); requestedRegion.SetIndex(2, 0); requestedRegion.SetSize(0, input->GetDimension(0)); requestedRegion.SetSize(1, input->GetDimension(1)); requestedRegion.SetSize(2, input->GetDimension(2)); requestedRegion.SetIndex( m_SliceDimension, m_SliceIndex ); // only one slice needed requestedRegion.SetSize( m_SliceDimension, 1 ); input->SetRequestedRegion( &requestedRegion ); } /* * Generate the information decribing the output data. The default * implementation of this method will copy information from the input to the * output. A filter may override this method if its output will have different * information than its input. For instance, a filter that shrinks an image will * need to provide an implementation for this method that changes the spacing of * the pixels. Such filters should call their superclass' implementation of this * method prior to changing the information values they need (i.e. * GenerateOutputInformation() should call * Superclass::GenerateOutputInformation() prior to changing the information. */ void mitk::ExtractImageFilter::GenerateOutputInformation() { Image::Pointer output = this->GetOutput(); Image::ConstPointer input = this->GetInput(); if (input.IsNull()) return; if ( m_SliceDimension >= input->GetDimension() && input->GetDimension() != 2 ) { MITK_ERROR << "mitk::ExtractImageFilter:GenerateOutputInformation m_SliceDimension == " << m_SliceDimension << " makes no sense with an " << input->GetDimension() << "D image." << std::endl; itkExceptionMacro("This is not a sensible value for m_SliceDimension."); return; } unsigned int sliceDimension( m_SliceDimension ); if ( input->GetDimension() == 2) { sliceDimension = 2; } unsigned int tmpDimensions[2]; switch ( sliceDimension ) { default: case 2: // orientation = PlaneGeometry::Transversal; tmpDimensions[0] = input->GetDimension(0); tmpDimensions[1] = input->GetDimension(1); break; case 1: // orientation = PlaneGeometry::Frontal; tmpDimensions[0] = input->GetDimension(0); tmpDimensions[1] = input->GetDimension(2); break; case 0: // orientation = PlaneGeometry::Sagittal; tmpDimensions[0] = input->GetDimension(1); tmpDimensions[1] = input->GetDimension(2); break; } output->Initialize(input->GetPixelType(), 2, tmpDimensions, 1 /*input->GetNumberOfChannels()*/); // initialize the spacing of the output /* Vector3D spacing = input->GetSlicedGeometry()->GetSpacing(); if(input->GetDimension()>=2) spacing[2]=spacing[1]; else spacing[2] = 1.0; output->GetSlicedGeometry()->SetSpacing(spacing); */ output->SetPropertyList(input->GetPropertyList()->Clone()); } diff --git a/Modules/ImageExtraction/mitkExtractImageFilter.h b/Modules/ImageExtraction/mitkExtractImageFilter.h index ef707c168c..8196563b7f 100644 --- a/Modules/ImageExtraction/mitkExtractImageFilter.h +++ b/Modules/ImageExtraction/mitkExtractImageFilter.h @@ -1,98 +1,101 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #ifndef mitkExtractImageFilter_h_Included #define mitkExtractImageFilter_h_Included #include "mitkCommon.h" #include "ImageExtractionExports.h" #include "mitkImageToImageFilter.h" #include "itkImage.h" namespace mitk { /** + \deprecated This class is deprecated. Use mitk::ExtractSliceFilter instead. + \sa ExtractSliceFilter + \brief Extracts a 2D slice from a 3D image. \sa SegTool2D \sa OverwriteSliceImageFilter \ingroup Process \ingroup ToolManagerEtAl There is a separate page describing the general design of QmitkInteractiveSegmentation: \ref QmitkSegmentationTechnicalPage This class takes a 3D mitk::Image as input and tries to extract one slice from it. Two parameters determine which slice is extracted: the "slice dimension" is that one, which is constant for all points in the plane, e.g. transversal would mean 2. The "slice index" is the slice index in the image direction you specified with "affected dimension". Indices count from zero. Output will not be set if there was a problem extracting the desired slice. Last contributor: $Author$ */ class ImageExtraction_EXPORT ExtractImageFilter : public ImageToImageFilter { public: mitkClassMacro(ExtractImageFilter, ImageToImageFilter); itkNewMacro(ExtractImageFilter); /** \brief Which slice to extract (first one has index 0). */ itkSetMacro(SliceIndex, unsigned int); itkGetConstMacro(SliceIndex, unsigned int); /** \brief The orientation of the slice to be extracted. \a Parameter SliceDimension Number of the dimension which is constant for all pixels of the desired slice (e.g. 2 for transversal) */ itkSetMacro(SliceDimension, unsigned int); itkGetConstMacro(SliceDimension, unsigned int); /** \brief Time step of the image to be extracted. */ itkSetMacro(TimeStep, unsigned int); itkGetConstMacro(TimeStep, unsigned int); protected: ExtractImageFilter(); // purposely hidden virtual ~ExtractImageFilter(); virtual void GenerateOutputInformation(); virtual void GenerateInputRequestedRegion(); virtual void GenerateData(); template void ItkImageProcessing( itk::Image* image ); unsigned int m_SliceIndex; unsigned int m_SliceDimension; unsigned int m_TimeStep; }; } // namespace #endif