diff --git a/Core/Code/Algorithms/mitkExtractSliceFilter.cpp b/Core/Code/Algorithms/mitkExtractSliceFilter.cpp index 1eefe28de0..a86080ecbd 100644 --- a/Core/Code/Algorithms/mitkExtractSliceFilter.cpp +++ b/Core/Code/Algorithms/mitkExtractSliceFilter.cpp @@ -1,563 +1,565 @@ /*========================================================================= Program: Medical Imaging & Interaction Toolkit Language: C++ Date: $Date$ Version: $Revision: $ Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. See MITKCopyright.txt or http://www.mitk.org/copyright.html for details. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the above copyright notices for more information. =========================================================================*/ #include "mitkExtractSliceFilter.h" #include #include #include #include #include #include mitk::ExtractSliceFilter::ExtractSliceFilter(vtkImageReslice* reslicer ){ if(reslicer == NULL){ m_Reslicer = vtkSmartPointer::New(); } else { m_Reslicer = reslicer; } m_TimeStep = 0; m_Reslicer->ReleaseDataFlagOn(); m_InterpolationMode = ExtractSliceFilter::RESLICE_NEAREST; m_ResliceTransform = NULL; m_InPlaneResampleExtentByGeometry = false; m_OutPutSpacing = new mitk::ScalarType[2]; m_OutputDimension = 2; m_ZSpacing = 1.0; m_ZMin = 0; m_ZMax = 0; m_VtkOutputRequested = false; } mitk::ExtractSliceFilter::~ExtractSliceFilter(){ m_ResliceTransform = NULL; m_WorldGeometry = NULL; delete [] m_OutPutSpacing; } void mitk::ExtractSliceFilter::GenerateOutputInformation(){ //TODO try figure out how to set the specs of the slice before it is actually extracted /*Image::Pointer output = this->GetOutput(); Image::ConstPointer input = this->GetInput(); if (input.IsNull()) return; unsigned int dimensions[2]; dimensions[0] = m_WorldGeometry->GetExtent(0); dimensions[1] = m_WorldGeometry->GetExtent(1); output->Initialize(input->GetPixelType(), 2, dimensions, 1);*/ } void mitk::ExtractSliceFilter::GenerateInputRequestedRegion(){ //As we want all pixel information fo the image in our plane, the requested region //is set to the largest possible region in the image. //This is needed because an oblique plane has a larger extent then the image //and the in pipeline it is checked via PropagateResquestedRegion(). But the //extent of the slice is actually fitting because it is oblique within the image. ImageToImageFilter::InputImagePointer input = const_cast< ImageToImageFilter::InputImageType* > ( this->GetInput() ); input->SetRequestedRegionToLargestPossibleRegion(); } mitk::ScalarType* mitk::ExtractSliceFilter::GetOutputSpacing(){ return m_OutPutSpacing; } void mitk::ExtractSliceFilter::GenerateData(){ mitk::Image *input = const_cast< mitk::Image * >( this->GetInput() ); if (!input) { MITK_ERROR << "mitk::ExtractSliceFilter: No input image available. Please set the input!" << std::endl; itkExceptionMacro("mitk::ExtractSliceFilter: No input image available. Please set the input!"); return; } if(!m_WorldGeometry) { MITK_ERROR << "mitk::ExtractSliceFilter: No Geometry for reslicing available." << std::endl; itkExceptionMacro("mitk::ExtractSliceFilter: No Geometry for reslicing available."); return; } const TimeSlicedGeometry *inputTimeGeometry = this->GetInput()->GetTimeSlicedGeometry(); if ( ( inputTimeGeometry == NULL ) || ( inputTimeGeometry->GetTimeSteps() == 0 ) ) { itkWarningMacro(<<"Error reading input image TimeSlicedGeometry."); return; } // is it a valid timeStep? if ( inputTimeGeometry->IsValidTime( m_TimeStep ) == false ) { itkWarningMacro(<<"This is not a valid timestep: "<< m_TimeStep ); return; } // check if there is something to display. if ( ! input->IsVolumeSet( m_TimeStep ) ) { itkWarningMacro(<<"No volume data existent at given timestep "<< m_TimeStep ); return; } /*================#BEGIN setup vtkImageRslice properties================*/ Point3D origin; Vector3D right, bottom, normal; double widthInMM, heightInMM; Vector2D extent; const PlaneGeometry* planeGeometry = dynamic_cast(m_WorldGeometry); if ( planeGeometry != NULL ) { //if the worldGeomatry is a PlaneGeometry everthing is straight forward 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. Vector3D rightInIndex, bottomInIndex; inputTimeGeometry->GetGeometry3D( m_TimeStep )->WorldToIndex( right, rightInIndex ); inputTimeGeometry->GetGeometry3D( m_TimeStep )->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 ); m_OutPutSpacing[0] = widthInMM / extent[0]; m_OutPutSpacing[1] = heightInMM / extent[1]; right.Normalize(); bottom.Normalize(); normal.Normalize(); /* * Transform the origin to center based coordinates. * Note: * This is needed besause vtk's origin is center based too (!!!) ( see 'The VTK book' page 88 ) * and the worldGeometry surrouding the image is no imageGeometry. So the worldGeometry * has its origin at the corner of the voxel and needs to be transformed. */ origin += right * ( m_OutPutSpacing[0] * 0.5 ); origin += bottom * ( m_OutPutSpacing[1] * 0.5 ); //set the tranform for reslicing. // Use inverse transform of the input geometry for reslicing the 3D image. // This is needed if the image volume already transformed if(m_ResliceTransform != NULL) m_Reslicer->SetResliceTransform(m_ResliceTransform->GetVtkTransform()->GetLinearInverse()); // Set background level to TRANSLUCENT (see Geometry2DDataVtkMapper3D), // else the background of the image turns out gray m_Reslicer->SetBackgroundLevel( -32768 ); } else{ //Code for curved planes, mostly taken 1:1 from imageVtkMapper2D and not tested yet. // Do we have an AbstractTransformGeometry? // This is the case for AbstractTransformGeometry's (e.g. a ThinPlateSplineCurvedGeometry ) const mitk::AbstractTransformGeometry* abstractGeometry = dynamic_cast< const AbstractTransformGeometry * >(m_WorldGeometry); if(abstractGeometry != NULL) { m_ResliceTransform = abstractGeometry; extent[0] = abstractGeometry->GetParametricExtent(0); extent[1] = abstractGeometry->GetParametricExtent(1); widthInMM = abstractGeometry->GetParametricExtentInMM(0); heightInMM = abstractGeometry->GetParametricExtentInMM(1); m_OutPutSpacing[0] = widthInMM / extent[0]; m_OutPutSpacing[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 vtkSmartPointer composedResliceTransform = vtkGeneralTransform::New(); composedResliceTransform->Identity(); composedResliceTransform->Concatenate( inputTimeGeometry->GetGeometry3D( m_TimeStep )->GetVtkTransform()->GetLinearInverse() ); composedResliceTransform->Concatenate( abstractGeometry->GetVtkAbstractTransform() ); m_Reslicer->SetResliceTransform( composedResliceTransform ); composedResliceTransform->UnRegister( NULL ); // decrease RC // Set background level to BLACK instead of translucent, to avoid // boundary artifacts (see Geometry2DDataVtkMapper3D) m_Reslicer->SetBackgroundLevel( -1023 ); } else { itkExceptionMacro("mitk::ExtractSliceFilter: No fitting geometry for reslice axis!"); return; } } if(m_ResliceTransform != NULL){ //if the resliceTransform is set the reslice axis are recalculated. //Thus the geometry information is not fitting. Therefor a unitSpacingFilter //is used to set up a global spacing of 1 and compensate the transform. vtkSmartPointer unitSpacingImageFilter = vtkSmartPointer::New() ; unitSpacingImageFilter->ReleaseDataFlagOn(); unitSpacingImageFilter->SetOutputSpacing( 1.0, 1.0, 1.0 ); unitSpacingImageFilter->SetInput( input->GetVtkImageData(m_TimeStep) ); m_Reslicer->SetInput(unitSpacingImageFilter->GetOutput() ); } else { //if no tranform is set the image can be used directly m_Reslicer->SetInput(input->GetVtkImageData(m_TimeStep)); } /*setup the plane where vktImageReslice extracts the slice*/ //ResliceAxesOrigin is the ancor point of the plane double originInVtk[3]; itk2vtk(origin, originInVtk); m_Reslicer->SetResliceAxesOrigin(originInVtk); //the cosines define the plane: x and y are the direction vectors, n is the planes normal //this specifies a matrix 3x3 // x1 y1 n1 // x2 y2 n2 // x3 y3 n3 double cosines[9]; vnl2vtk(right.GetVnlVector(), cosines);//x vnl2vtk(bottom.GetVnlVector(), cosines + 3);//y vnl2vtk(normal.GetVnlVector(), cosines + 6);//n m_Reslicer->SetResliceAxesDirectionCosines(cosines); //we only have one slice, not a volume m_Reslicer->SetOutputDimensionality(m_OutputDimension); //set the interpolation mode for slicing switch(this->m_InterpolationMode){ case RESLICE_NEAREST: m_Reslicer->SetInterpolationModeToNearestNeighbor(); break; case RESLICE_LINEAR: m_Reslicer->SetInterpolationModeToLinear(); break; case RESLICE_CUBIC: m_Reslicer->SetInterpolationModeToCubic(); default: //the default interpolation used by mitk m_Reslicer->SetInterpolationModeToNearestNeighbor(); } /*========== BEGIN setup extent of the slice ==========*/ int xMin, xMax, yMin, yMax; bool boundsInitialized = false; if(m_WorldGeometry->GetReferenceGeometry()){ vtkFloatingPointType sliceBounds[6]; for ( int i = 0; i < 6; ++i ) { sliceBounds[i] = 0.0; } - this->CalculateClippedPlaneBounds( m_WorldGeometry->GetReferenceGeometry(), planeGeometry, sliceBounds ); + this->GetClippedPlaneBounds( m_WorldGeometry->GetReferenceGeometry(), planeGeometry, sliceBounds ); // Calculate output extent (integer values) xMin = static_cast< int >( sliceBounds[0] / m_OutPutSpacing[0] + 0.5 ); xMax = static_cast< int >( sliceBounds[1] / m_OutPutSpacing[0] + 0.5 ); yMin = static_cast< int >( sliceBounds[2] / m_OutPutSpacing[1] + 0.5 ); yMax = static_cast< int >( sliceBounds[3] / m_OutPutSpacing[1] + 0.5 ); }else { // If no reference geometry is available, we also don't know about the // maximum plane size; xMin = yMin = 0; xMax = static_cast< int >( extent[0]); yMax = static_cast< int >( extent[1]); } m_Reslicer->SetOutputExtent(xMin, xMax-1, yMin, yMax-1, m_ZMin, m_ZMax ); /*========== END setup extent of the slice ==========*/ m_Reslicer->SetOutputOrigin( 0.0, 0.0, 0.0 ); m_Reslicer->SetOutputSpacing( m_OutPutSpacing[0], m_OutPutSpacing[1], m_ZSpacing ); //TODO check the following lines, they are responsible wether vtk error outputs appear or not m_Reslicer->UpdateWholeExtent(); //this produces a bad allocation error for 2D images //m_Reslicer->GetOutput()->UpdateInformation(); //m_Reslicer->GetOutput()->SetUpdateExtentToWholeExtent(); //start the pipeline m_Reslicer->Update(); /*================ #END setup vtkImageRslice properties================*/ if(m_VtkOutputRequested){ return; //no converting to mitk //no mitk geometry will be set, as the output is vtkImageData only!!! } else { /*================ #BEGIN Get the slice from vtkImageReslice and convert it to mit::Image================*/ vtkImageData* reslicedImage; reslicedImage = m_Reslicer->GetOutput(); if(!reslicedImage) { itkWarningMacro(<<"Reslicer returned empty image"); return; } mitk::Image::Pointer resultImage = this->GetOutput(); //initialize resultimage with the specs of the vtkImageData object returned from vtkImageReslice resultImage->Initialize(reslicedImage); //transfer the voxel data resultImage->SetVolume(reslicedImage->GetScalarPointer()); //set the geometry from current worldgeometry for the reusultimage //this is needed that the image has the correct mitk geometry //the originalGeometry is the Geometry of the result slice AffineGeometryFrame3D::Pointer originalGeometryAGF = m_WorldGeometry->Clone(); Geometry2D::Pointer originalGeometry = dynamic_cast( originalGeometryAGF.GetPointer() ); //the origin of the worldGeometry is transformed to center based coordinates to be an imageGeometry Point3D sliceOrigin = originalGeometry->GetOrigin(); sliceOrigin += right * ( m_OutPutSpacing[0] * 0.5 ); sliceOrigin += bottom * ( m_OutPutSpacing[1] * 0.5 ); //a worldGeometry is no imageGeometry, thus it is manually set to true originalGeometry->ImageGeometryOn(); originalGeometry->SetOrigin(sliceOrigin); resultImage->SetGeometry( originalGeometry ); //resultImage->GetGeometry()->TransferItkToVtkTransform(); resultImage->GetGeometry()->Modified(); /*================ #END Get the slice from vtkImageReslice and convert it to mitk Image================*/ } } bool mitk::ExtractSliceFilter::GetClippedPlaneBounds(vtkFloatingPointType bounds[6]){ if(!m_WorldGeometry || !this->GetInput()) return false; - return this->CalculateClippedPlaneBounds(m_WorldGeometry->GetReferenceGeometry(), dynamic_cast< const PlaneGeometry * >( m_WorldGeometry ), bounds); + return this->GetClippedPlaneBounds(m_WorldGeometry->GetReferenceGeometry(), dynamic_cast< const PlaneGeometry * >( m_WorldGeometry ), bounds); } bool mitk::ExtractSliceFilter::GetClippedPlaneBounds( const Geometry3D *boundingGeometry, const PlaneGeometry *planeGeometry, vtkFloatingPointType *bounds ) { - return this->CalculateClippedPlaneBounds(boundingGeometry, planeGeometry, bounds); + bool b = this->CalculateClippedPlaneBounds(boundingGeometry, planeGeometry, bounds); + + return b; } bool mitk::ExtractSliceFilter ::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(); BoundingBox::PointType bbCenter = boundingBox->GetCenter(); 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::ExtractSliceFilter ::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; } \ No newline at end of file