diff --git a/Core/Code/Algorithms/mitkExtractSliceFilter.cpp b/Core/Code/Algorithms/mitkExtractSliceFilter.cpp new file mode 100644 index 0000000000..1f0effc099 --- /dev/null +++ b/Core/Code/Algorithms/mitkExtractSliceFilter.cpp @@ -0,0 +1,593 @@ +/*========================================================================= + +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; + + vtkFloatingPointType sliceBounds[6]; + if(m_WorldGeometry->GetReferenceGeometry()){ + + + for ( int i = 0; i < 6; ++i ) + { + sliceBounds[i] = 0.0; + } + 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() ); + + originalGeometry->GetIndexToWorldTransform()->SetMatrix(m_WorldGeometry->GetIndexToWorldTransform()->GetMatrix()); + + //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(); + + + /*At this point we have to adjust the geometry because the origin isn't correct. + The wrong origin is related to the rotation of the current world geometry plane. + This causes errors on transfering world to index coordinates. We just shift the + origin in each direction about the amount of the expanding (needed while rotating + the plane). + */ + Vector3D axis0 = originalGeometry->GetAxisVector(0); + Vector3D axis1 = originalGeometry->GetAxisVector(1); + axis0.Normalize(); + axis1.Normalize(); + + + //adapt the origin. Note that for orthogonal planes the minima are '0' and thus the origin stays the same. + sliceOrigin += (axis0 * (xMin * m_OutPutSpacing[0])) + (axis1 * (yMin * m_OutPutSpacing[1])); + + originalGeometry->SetOrigin(sliceOrigin); + + originalGeometry->Modified(); + + + resultImage->SetGeometry( originalGeometry ); + + + /*the bounds as well as the extent of the worldGeometry are not adapted correctly during crosshair rotation. + This is only a quick fix and has to be evaluated. + The new bounds are set via the max values of the calcuted slice extent. It will look like [ 0, x, 0, y, 0, 1]. + */ + mitk::BoundingBox::BoundsArrayType boundsCopy; + boundsCopy[0] = boundsCopy[2] = boundsCopy[4] = 0; + boundsCopy[5] = 1; + boundsCopy[1] = xMax - xMin; + boundsCopy[3] = yMax - yMin; + resultImage->GetGeometry()->SetBounds(boundsCopy); + + /*================ #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->GetClippedPlaneBounds(m_WorldGeometry->GetReferenceGeometry(), dynamic_cast< const PlaneGeometry * >( m_WorldGeometry ), bounds); +} + + +bool mitk::ExtractSliceFilter::GetClippedPlaneBounds( const Geometry3D *boundingGeometry, + const PlaneGeometry *planeGeometry, vtkFloatingPointType *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 diff --git a/Core/Code/Algorithms/mitkExtractSliceFilter.h b/Core/Code/Algorithms/mitkExtractSliceFilter.h new file mode 100644 index 0000000000..7645d8432a --- /dev/null +++ b/Core/Code/Algorithms/mitkExtractSliceFilter.h @@ -0,0 +1,185 @@ +/*========================================================================= + +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. + +=========================================================================*/ + +#ifndef mitkExtractSliceFilter_h_Included +#define mitkExtractSliceFilter_h_Included + + +#include "MitkExports.h" +#include "mitkImageToImageFilter.h" +#include + +#include +#include +#include +#include +#include +#include + +namespace mitk +{ + /** + \brief ExtractSliceFilter extracts a 2D abitrary oriented slice from a 3D volume. + + The filter can reslice in all orthogonal planes such as sagittal, coronal and transversal, + and is also able to reslice a abitrary oriented oblique plane. + Curved planes are specified via an AbstractTransformGeometry as the input worldgeometry. + + The convinient workflow is: + 1. Set an image as input. + 2. Set the worldGeometry2D. This defines a grid where the slice is being extracted + 3. And then start the pipeline. + + There are a few more properties that can be set to modify the behavior of the slicing. + The properties are: + - interpolation mode either Nearestneighbor, Linear or Cubic. + - a transform this is a convinient way to adapt the reslice axis for the case + that the image is transformed e.g. rotated. + - time step the time step in a timesliced volume. + - resample by geometry wether the resampling grid corresponds to the specs of the + worldgeometry or is directly derived from the input image + + By default the properties are set to: + - interpolation mode Nearestneighbor. + - a transform NULL (No transform is set). + - time step 0. + - resample by geometry false (Corresponds to input image). + */ + class MITK_CORE_EXPORT ExtractSliceFilter : public ImageToImageFilter + { + public: + + mitkClassMacro(ExtractSliceFilter, ImageToImageFilter); + itkNewMacro(ExtractSliceFilter); + + mitkNewMacro1Param(Self, vtkImageReslice*); + + /** \brief Set the axis where to reslice at.*/ + void SetWorldGeometry(const Geometry2D* geometry ){ this->m_WorldGeometry = geometry; } + + /** \brief Set the time step in the 4D volume */ + void SetTimeStep( unsigned int timestep){ this->m_TimeStep = timestep; } + unsigned int GetTimeStep(){ return this->m_TimeStep; } + + /** \brief Set a transform for the reslice axes. + * This transform is needed if the image volume itself is transformed. (Effects the reslice axis) + */ + void SetResliceTransformByGeometry(const Geometry3D* transform){ this->m_ResliceTransform = transform; } + + /** \brief Resampling grid corresponds to: false->image true->worldgeometry*/ + void SetInPlaneResampleExtentByGeometry(bool inPlaneResampleExtentByGeometry){ this->m_InPlaneResampleExtentByGeometry = inPlaneResampleExtentByGeometry; } + + /** \brief Sets the output dimension of the slice*/ + void SetOutputDimensionality(unsigned int dimension){ this->m_OutputDimension = dimension; } + + /** \brief Set the spacing in z direction manually. + * Required if the outputDimension is > 2. + */ + void SetOutputSpacingZDirection(double zSpacing){ this->m_ZSpacing = zSpacing; } + + /** \brief Set the extent in pixel for direction z manualy. + Required if the output dimension is > 2. + */ + void SetOutputExtentZDirection(int zMin, int zMax) { this->m_ZMin = zMin; this->m_ZMax = zMax; } + + /** \brief Get the bounding box of the slice [xMin, xMax, yMin, yMax, zMin, zMax] + * The method uses the input of the filter to calculate the bounds. + * It is recommended to use + * GetClippedPlaneBounds(const Geometry3D*, const PlaneGeometry*, vtkFloatingPointType*) + * if you are not sure about the input. + */ + bool GetClippedPlaneBounds(double bounds[6]); + + /** \brief Get the bounding box of the slice [xMin, xMax, yMin, yMax, zMin, zMax]*/ + bool GetClippedPlaneBounds( const Geometry3D *boundingGeometry, + const PlaneGeometry *planeGeometry, vtkFloatingPointType *bounds ); + + /** \brief Get the spacing of the slice. returns mitk::ScalarType[2] */ + mitk::ScalarType* GetOutputSpacing(); + + /** \brief Get Output as vtkImageData. + * Note: + * SetVtkOutputRequest(true) has to be called at least once before + * GetVtkOutput(). Otherwise the output is empty for the first update step. + */ + vtkImageData* GetVtkOutput(){ m_VtkOutputRequested = true; return m_Reslicer->GetOutput(); } + + /** Set VtkOutPutRequest to suppress the convertion of the image. + * It is suggested to use this with GetVtkOutput(). + * Note: + * SetVtkOutputRequest(true) has to be called at least once before + * GetVtkOutput(). Otherwise the output is empty for the first update step. + */ + void SetVtkOutputRequest(bool isRequested){ m_VtkOutputRequested = isRequested; } + + /** \brief Get the reslices axis matrix. + * Note: the axis are recalculated when calling SetResliceTransformByGeometry. + */ + vtkMatrix4x4* GetResliceAxes(){ + return this->m_Reslicer->GetResliceAxes(); + } + + enum ResliceInterpolation { RESLICE_NEAREST, RESLICE_LINEAR, RESLICE_CUBIC }; + + void SetInterpolationMode( ExtractSliceFilter::ResliceInterpolation interpolation){ this->m_InterpolationMode = interpolation; } + + protected: + ExtractSliceFilter(vtkImageReslice* reslicer = NULL); + virtual ~ExtractSliceFilter(); + + virtual void GenerateData(); + virtual void GenerateOutputInformation(); + virtual void GenerateInputRequestedRegion(); + + + const Geometry2D* m_WorldGeometry; + vtkSmartPointer m_Reslicer; + + unsigned int m_TimeStep; + + unsigned int m_OutputDimension; + + double m_ZSpacing; + + int m_ZMin; + + int m_ZMax; + + ResliceInterpolation m_InterpolationMode; + + const Geometry3D* m_ResliceTransform; + + bool m_InPlaneResampleExtentByGeometry;//Resampling grid corresponds to: false->image true->worldgeometry + + mitk::ScalarType* m_OutPutSpacing; + + + bool m_VtkOutputRequested; + + /** \brief Internal helper method for intersection testing used only in CalculateClippedPlaneBounds() */ + bool LineIntersectZero( vtkPoints *points, int p1, int p2, + vtkFloatingPointType *bounds ); + + /** \brief Calculate the bounding box of the resliced image. This is necessary for + * arbitrarily rotated planes in an image volume. A rotated plane (e.g. in swivel mode) + * will have a new bounding box, which needs to be calculated. */ + bool CalculateClippedPlaneBounds( const Geometry3D *boundingGeometry, + const PlaneGeometry *planeGeometry, vtkFloatingPointType *bounds ); + }; +} + +#endif // mitkExtractSliceFilter_h_Included \ No newline at end of file diff --git a/Core/Code/Rendering/mitkImageVtkMapper2D.cpp b/Core/Code/Rendering/mitkImageVtkMapper2D.cpp index 809d5e42d1..83a04c1562 100644 --- a/Core/Code/Rendering/mitkImageVtkMapper2D.cpp +++ b/Core/Code/Rendering/mitkImageVtkMapper2D.cpp @@ -1,1238 +1,959 @@ /*=================================================================== 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. ===================================================================*/ //MITK #include #include #include #include #include #include #include #include #include #include #include #include //#include #include #include "mitkImageStatisticsHolder.h" //MITK Rendering #include "mitkImageVtkMapper2D.h" #include "vtkMitkThickSlicesFilter.h" #include "vtkMitkApplyLevelWindowToRGBFilter.h" //VTK #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include //ITK #include mitk::ImageVtkMapper2D::ImageVtkMapper2D() { } mitk::ImageVtkMapper2D::~ImageVtkMapper2D() { //The 3D RW Mapper (Geometry2DDataVtkMapper3D) is listening to this event, //in order to delete the images from the 3D RW. this->InvokeEvent( itk::DeleteEvent() ); } //set the two points defining the textured plane according to the dimension and spacing void mitk::ImageVtkMapper2D::GeneratePlane(mitk::BaseRenderer* renderer, vtkFloatingPointType planeBounds[6]) { LocalStorage *localStorage = m_LSH.GetLocalStorage(renderer); + float depth = this->CalculateLayerDepth(renderer); //Set the origin to (xMin; yMin; depth) of the plane. This is necessary for obtaining the correct //plane size in crosshair rotation and swivel mode. localStorage->m_Plane->SetOrigin(planeBounds[0], planeBounds[2], depth); //These two points define the axes of the plane in combination with the origin. //Point 1 is the x-axis and point 2 the y-axis. //Each plane is transformed according to the view (transversal, coronal and saggital) afterwards. - localStorage->m_Plane->SetPoint1(planeBounds[1], planeBounds[2], depth); //P1: (xMax, yMin, depth) + localStorage->m_Plane->SetPoint1(planeBounds[1] , planeBounds[2], depth); //P1: (xMax, yMin, depth) localStorage->m_Plane->SetPoint2(planeBounds[0], planeBounds[3], depth); //P2: (xMin, yMax, depth) } float mitk::ImageVtkMapper2D::CalculateLayerDepth(mitk::BaseRenderer* renderer) { //get the clipping range to check how deep into z direction we can render images double maxRange = renderer->GetVtkRenderer()->GetActiveCamera()->GetClippingRange()[1]; //Due to a VTK bug, we cannot use the whole clipping range. /100 is empirically determined float depth = -maxRange*0.01; // divide by 100 int layer = 0; GetDataNode()->GetIntProperty( "layer", layer, renderer); //add the layer property for each image to render images with a higher layer on top of the others depth += layer*10; //*10: keep some room for each image (e.g. for QBalls in between) if(depth > 0.0f) { depth = 0.0f; MITK_WARN << "Layer value exceeds clipping range. Set to minimum instead."; } return depth; } const mitk::Image* mitk::ImageVtkMapper2D::GetInput( void ) { return static_cast< const mitk::Image * >( this->GetData() ); } vtkProp* mitk::ImageVtkMapper2D::GetVtkProp(mitk::BaseRenderer* renderer) { //return the actor corresponding to the renderer return m_LSH.GetLocalStorage(renderer)->m_Actor; } void mitk::ImageVtkMapper2D::MitkRenderOverlay(BaseRenderer* renderer) { if ( this->IsVisible(renderer)==false ) return; if ( this->GetVtkProp(renderer)->GetVisibility() ) { this->GetVtkProp(renderer)->RenderOverlay(renderer->GetVtkRenderer()); } } void mitk::ImageVtkMapper2D::MitkRenderOpaqueGeometry(BaseRenderer* renderer) { if ( this->IsVisible( renderer )==false ) return; if ( this->GetVtkProp(renderer)->GetVisibility() ) { this->GetVtkProp(renderer)->RenderOpaqueGeometry( renderer->GetVtkRenderer() ); } } void mitk::ImageVtkMapper2D::MitkRenderTranslucentGeometry(BaseRenderer* renderer) { if ( this->IsVisible(renderer)==false ) return; if ( this->GetVtkProp(renderer)->GetVisibility() ) { this->GetVtkProp(renderer)->RenderTranslucentPolygonalGeometry(renderer->GetVtkRenderer()); } } void mitk::ImageVtkMapper2D::MitkRenderVolumetricGeometry(BaseRenderer* renderer) { if(IsVisible(renderer)==false) return; if ( GetVtkProp(renderer)->GetVisibility() ) { this->GetVtkProp(renderer)->RenderVolumetricGeometry(renderer->GetVtkRenderer()); } } void mitk::ImageVtkMapper2D::GenerateDataForRenderer( mitk::BaseRenderer *renderer ) { LocalStorage *localStorage = m_LSH.GetLocalStorage(renderer); mitk::Image *input = const_cast< mitk::Image * >( this->GetInput() ); if ( input == NULL || input->IsInitialized() == false ) { return; } //check if there is a valid worldGeometry const Geometry2D *worldGeometry = renderer->GetCurrentWorldGeometry2D(); if( ( worldGeometry == NULL ) || ( !worldGeometry->IsValid() ) || ( !worldGeometry->HasReferenceGeometry() )) { return; } - // check if there is something to display - if ( !input->IsVolumeSet( this->GetTimestep() ) ) return; - input->Update(); - vtkImageData* inputData = input->GetVtkImageData( this->GetTimestep() ); - if ( inputData == NULL ) - { - return; - } - - // how big the area is in physical coordinates: widthInMM x heightInMM pixels - mitk::ScalarType widthInMM, heightInMM; - - // take transform of input image into account - const TimeSlicedGeometry *inputTimeGeometry = input->GetTimeSlicedGeometry(); - const Geometry3D* inputGeometry = inputTimeGeometry->GetGeometry3D( this->GetTimestep() ); - - // Bounds information for reslicing (only reuqired if reference geometry - // is present) - vtkFloatingPointType sliceBounds[6]; - bool boundsInitialized = false; - for ( int i = 0; i < 6; ++i ) - { - sliceBounds[i] = 0.0; - } - - //Extent (in pixels) of the image - Vector2D extent; - - // Do we have a simple PlaneGeometry? - // This is the "regular" case (e.g. slicing through an image axis-parallel or even oblique) - const PlaneGeometry *planeGeometry = dynamic_cast< const PlaneGeometry * >( worldGeometry ); - if ( planeGeometry != NULL ) - { - localStorage->m_Origin = planeGeometry->GetOrigin(); - localStorage->m_Right = planeGeometry->GetAxisVector( 0 ); // right = Extent of Image in mm (worldspace) - localStorage->m_Bottom = planeGeometry->GetAxisVector( 1 ); - localStorage->m_Normal = planeGeometry->GetNormal(); - - bool inPlaneResampleExtentByGeometry = false; - GetDataNode()->GetBoolProperty("in plane resample extent by geometry", inPlaneResampleExtentByGeometry, renderer); - - if ( 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] = worldGeometry->GetExtent( 0 ); - extent[1] = 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; - inputGeometry->WorldToIndex( localStorage->m_Right, rightInIndex ); - inputGeometry->WorldToIndex( localStorage->m_Bottom, bottomInIndex ); - extent[0] = rightInIndex.GetNorm(); - extent[1] = bottomInIndex.GetNorm(); - } - - // Get the extent of the current world geometry and calculate resampling - // spacing therefrom. - widthInMM = worldGeometry->GetExtentInMM( 0 ); - heightInMM = worldGeometry->GetExtentInMM( 1 ); - - localStorage->m_mmPerPixel[0] = widthInMM / extent[0]; - localStorage->m_mmPerPixel[1] = heightInMM / extent[1]; - - localStorage->m_Right.Normalize(); - localStorage->m_Bottom.Normalize(); - localStorage->m_Normal.Normalize(); - - //transform the origin to corner based coordinates, because VTK is corner based. - localStorage->m_Origin += localStorage->m_Right * ( localStorage->m_mmPerPixel[0] * 0.5 ); - localStorage->m_Origin += localStorage->m_Bottom * ( localStorage->m_mmPerPixel[1] * 0.5 ); - - // Use inverse transform of the input geometry for reslicing the 3D image - localStorage->m_Reslicer->SetResliceTransform( - inputGeometry->GetVtkTransform()->GetLinearInverse() ); - - // Set background level to TRANSLUCENT (see Geometry2DDataVtkMapper3D) - localStorage->m_Reslicer->SetBackgroundLevel( -32768 ); - - // 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. - if (worldGeometry->HasReferenceGeometry()) - { - boundsInitialized = this->CalculateClippedPlaneBounds( - worldGeometry->GetReferenceGeometry(), planeGeometry, sliceBounds ); - } - } - else - { - // Do we have an AbstractTransformGeometry? - // This is the case for AbstractTransformGeometry's (e.g. a thin-plate-spline transform) - const mitk::AbstractTransformGeometry* abstractGeometry = - dynamic_cast< const AbstractTransformGeometry * >(worldGeometry); - - if(abstractGeometry != NULL) - { - - extent[0] = abstractGeometry->GetParametricExtent(0); - extent[1] = abstractGeometry->GetParametricExtent(1); - widthInMM = abstractGeometry->GetParametricExtentInMM(0); - heightInMM = abstractGeometry->GetParametricExtentInMM(1); + //set main input for ExtractSliceFilter + localStorage->m_Reslicer->SetInput(input); + localStorage->m_Reslicer->SetWorldGeometry(worldGeometry); + localStorage->m_Reslicer->SetTimeStep( this->GetTimestep() ); - localStorage->m_mmPerPixel[0] = widthInMM / extent[0]; - localStorage->m_mmPerPixel[1] = heightInMM / extent[1]; - localStorage->m_Origin = abstractGeometry->GetPlane()->GetOrigin(); + //set the transformation of the image to adapt reslice axis + localStorage->m_Reslicer->SetResliceTransformByGeometry( input->GetTimeSlicedGeometry()->GetGeometry3D( this->GetTimestep() ) ); - localStorage->m_Right = abstractGeometry->GetPlane()->GetAxisVector(0); - localStorage->m_Right.Normalize(); - localStorage->m_Bottom = abstractGeometry->GetPlane()->GetAxisVector(1); - localStorage->m_Bottom.Normalize(); + //is the geometry of the slice based on the input image or the worldgeometry? + bool inPlaneResampleExtentByGeometry = false; + GetDataNode()->GetBoolProperty("in plane resample extent by geometry", inPlaneResampleExtentByGeometry, renderer); + localStorage->m_Reslicer->SetInPlaneResampleExtentByGeometry(inPlaneResampleExtentByGeometry); - localStorage->m_Normal = abstractGeometry->GetPlane()->GetNormal(); - localStorage->m_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() - ); - - localStorage->m_Reslicer->SetResliceTransform( composedResliceTransform ); - composedResliceTransform->UnRegister( NULL ); // decrease RC - - // Set background level to BLACK instead of translucent, to avoid - // boundary artifacts (see Geometry2DDataVtkMapper3D) - localStorage->m_Reslicer->SetBackgroundLevel( -1023 ); - } - else - { - //no geometry => we can't reslice - return; - } - } - - // Make sure that the image to display has a certain minimum size. - if ( (extent[0] <= 2) && (extent[1] <= 2) ) - { - return; - } - - //### begin set reslice interpolation // Initialize the interpolation mode for resampling; switch to nearest // neighbor if the input image is too small. if ( (input->GetDimension() >= 3) && (input->GetDimension(2) > 1) ) { VtkResliceInterpolationProperty *resliceInterpolationProperty; this->GetDataNode()->GetProperty( - resliceInterpolationProperty, "reslice interpolation" ); + resliceInterpolationProperty, "reslice interpolation" ); int interpolationMode = VTK_RESLICE_NEAREST; if ( resliceInterpolationProperty != NULL ) { interpolationMode = resliceInterpolationProperty->GetInterpolation(); } switch ( interpolationMode ) { case VTK_RESLICE_NEAREST: - localStorage->m_Reslicer->SetInterpolationModeToNearestNeighbor(); + localStorage->m_Reslicer->SetInterpolationMode(ExtractSliceFilter::RESLICE_NEAREST); break; case VTK_RESLICE_LINEAR: - localStorage->m_Reslicer->SetInterpolationModeToLinear(); + localStorage->m_Reslicer->SetInterpolationMode(ExtractSliceFilter::RESLICE_LINEAR); break; case VTK_RESLICE_CUBIC: - localStorage->m_Reslicer->SetInterpolationModeToCubic(); + localStorage->m_Reslicer->SetInterpolationMode(ExtractSliceFilter::RESLICE_CUBIC); break; } } else { - localStorage->m_Reslicer->SetInterpolationModeToNearestNeighbor(); + localStorage->m_Reslicer->SetInterpolationMode(ExtractSliceFilter::RESLICE_NEAREST); } - //### end set reslice interpolation + + //set the vtk output property to true, makes sure that no unneeded mitk image convertion + //is done. + localStorage->m_Reslicer->SetVtkOutputRequest(true); + //Thickslicing int thickSlicesMode = 0; int thickSlicesNum = 1; - // Thick slices parameters - if( inputData->GetNumberOfScalarComponents() == 1 ) // for now only single component are allowed + if( input->GetPixelType().GetNumberOfComponents() == 1 ) // for now only single component are allowed { DataNode *dn=renderer->GetCurrentWorldGeometry2DNode(); if(dn) { ResliceMethodProperty *resliceMethodEnumProperty=0; if( dn->GetProperty( resliceMethodEnumProperty, "reslice.thickslices" ) && resliceMethodEnumProperty ) thickSlicesMode = resliceMethodEnumProperty->GetValueAsId(); IntProperty *intProperty=0; if( dn->GetProperty( intProperty, "reslice.thickslices.num" ) && intProperty ) { thickSlicesNum = intProperty->GetValue(); if(thickSlicesNum < 1) thickSlicesNum=1; if(thickSlicesNum > 10) thickSlicesNum=10; } } else { MITK_WARN << "no associated widget plane data tree node found"; } } - localStorage->m_UnitSpacingImageFilter->SetInput( inputData ); - localStorage->m_Reslicer->SetInput( localStorage->m_UnitSpacingImageFilter->GetOutput() ); - - //number of pixels per mm in x- and y-direction of the resampled - Vector2D pixelsPerMM; - pixelsPerMM[0] = 1.0 / localStorage->m_mmPerPixel[0]; - pixelsPerMM[1] = 1.0 / localStorage->m_mmPerPixel[1]; - - //calulate the originArray and the orientations for the reslice-filter - double originArray[3]; - itk2vtk( localStorage->m_Origin, originArray ); - - localStorage->m_Reslicer->SetResliceAxesOrigin( originArray ); - - double cosines[9]; - // direction of the X-axis of the sampled result - vnl2vtk( localStorage->m_Right.Get_vnl_vector(), cosines ); - // direction of the Y-axis of the sampled result - vnl2vtk( localStorage->m_Bottom.Get_vnl_vector(), cosines + 3 );//fill next 3 elements - // normal of the plane - vnl2vtk( localStorage->m_Normal.Get_vnl_vector(), cosines + 6 );//fill the last 3 elements - - localStorage->m_Reslicer->SetResliceAxesDirectionCosines( cosines ); - - int xMin, xMax, yMin, yMax; - if ( boundsInitialized ) - { - // Calculate output extent (integer values) - xMin = static_cast< int >( sliceBounds[0] / localStorage->m_mmPerPixel[0] + 0.5 ); - xMax = static_cast< int >( sliceBounds[1] / localStorage->m_mmPerPixel[0] + 0.5 ); - yMin = static_cast< int >( sliceBounds[2] / localStorage->m_mmPerPixel[1] + 0.5 ); - yMax = static_cast< int >( sliceBounds[3] / localStorage->m_mmPerPixel[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] - - pixelsPerMM[0] + 0.5); - yMax = static_cast< int >( extent[1] - - pixelsPerMM[1] + 0.5); - } - - // Disallow huge dimensions - if ( (xMax-xMin) * (yMax-yMin) > 4096*4096 ) - { - return; - } - // Calculate dataset spacing in plane z direction (NOT spacing of current - // world geometry) - double dataZSpacing = 1.0; + if(thickSlicesMode > 0) + { + double dataZSpacing = 1.0; + + Vector3D normInIndex, normal; + const PlaneGeometry *planeGeometry = dynamic_cast< const PlaneGeometry * >( worldGeometry ); + if ( planeGeometry != NULL ){ + normal = planeGeometry->GetNormal(); + }else{ + const mitk::AbstractTransformGeometry* abstractGeometry = dynamic_cast< const AbstractTransformGeometry * >(worldGeometry); + if(abstractGeometry != NULL) + normal = abstractGeometry->GetPlane()->GetNormal(); + else + return; //no fitting geometry set + } + normal.Normalize(); - Vector3D normInIndex; - inputGeometry->WorldToIndex( localStorage->m_Normal, normInIndex ); + input->GetTimeSlicedGeometry()->GetGeometry3D( this->GetTimestep() )->WorldToIndex( normal, normInIndex ); - if(thickSlicesMode > 0) - { dataZSpacing = 1.0 / normInIndex.GetNorm(); + localStorage->m_Reslicer->SetOutputDimensionality( 3 ); - localStorage->m_Reslicer->SetOutputExtent( xMin, xMax-1, yMin, yMax-1, -thickSlicesNum, 0+thickSlicesNum ); - } - else - { - localStorage->m_Reslicer->SetOutputDimensionality( 2 ); - localStorage->m_Reslicer->SetOutputExtent( xMin, xMax-1, yMin, yMax-1, 0, 0 ); - } + localStorage->m_Reslicer->SetOutputSpacingZDirection(dataZSpacing); + localStorage->m_Reslicer->SetOutputExtentZDirection( -thickSlicesNum, 0+thickSlicesNum ); - localStorage->m_Reslicer->SetOutputOrigin( 0.0, 0.0, 0.0 ); - localStorage->m_Reslicer->SetOutputSpacing( localStorage->m_mmPerPixel[0], localStorage->m_mmPerPixel[1], dataZSpacing ); - // xMax and yMax are meant exclusive until now, whereas - // SetOutputExtent wants an inclusive bound. Thus, we need - // to subtract 1. - - //Make sure the updateExtent is equal to the wholeExtent of vtkImageReslice, - //else the updateExtent is one step ahead of the wholeExtent and vtk errors - //are thrown. - localStorage->m_Reslicer->UpdateWholeExtent(); - - // 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. - if(thickSlicesMode>0) - { + // 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. localStorage->m_TSFilter->SetThickSliceMode( thickSlicesMode-1 ); - localStorage->m_TSFilter->SetInput( localStorage->m_Reslicer->GetOutput() ); + localStorage->m_TSFilter->SetInput( localStorage->m_Reslicer->GetVtkOutput() ); + + //vtkFilter=>mitkFilter=>vtkFilter update mechanism will fail without calling manually + localStorage->m_Reslicer->Modified(); + localStorage->m_Reslicer->Update(); + localStorage->m_TSFilter->Modified(); localStorage->m_TSFilter->Update(); localStorage->m_ReslicedImage = localStorage->m_TSFilter->GetOutput(); } else { + //this is needed when thick mode was enable bevore. These variable have to be reset to default values + localStorage->m_Reslicer->SetOutputDimensionality( 2 ); + localStorage->m_Reslicer->SetOutputSpacingZDirection(1.0); + localStorage->m_Reslicer->SetOutputExtentZDirection( 0, 0 ); + + localStorage->m_Reslicer->Modified(); - localStorage->m_Reslicer->Update(); - localStorage->m_ReslicedImage = localStorage->m_Reslicer->GetOutput(); + //start the pipeline with updating the largest possible, needed if the geometry of the input has changed + localStorage->m_Reslicer->UpdateLargestPossibleRegion(); + localStorage->m_ReslicedImage = localStorage->m_Reslicer->GetVtkOutput(); } - if((localStorage->m_ReslicedImage == NULL) || (localStorage->m_ReslicedImage->GetDataDimension() < 1)) + + // Bounds information for reslicing (only reuqired if reference geometry + // is present) + //this used for generating a vtkPLaneSource with the right size + vtkFloatingPointType sliceBounds[6]; + for ( int i = 0; i < 6; ++i ) { - MITK_WARN << "reslicer returned empty image"; - return; + sliceBounds[i] = 0.0; } + localStorage->m_Reslicer->GetClippedPlaneBounds(sliceBounds); + + //get the spacing of the slice + localStorage->m_mmPerPixel = localStorage->m_Reslicer->GetOutputSpacing(); + //get the number of scalar components to distinguish between different image types int numberOfComponents = localStorage->m_ReslicedImage->GetNumberOfScalarComponents(); //get the binary property bool binary = false; bool binaryOutline = false; this->GetDataNode()->GetBoolProperty( "binary", binary, renderer ); if(binary) //binary image { this->GetDataNode()->GetBoolProperty( "outline binary", binaryOutline, renderer ); if(binaryOutline) //contour rendering { if ( this->GetInput()->GetPixelType().GetBpe() <= 8 ) { //generate ontours/outlines localStorage->m_OutlinePolyData = CreateOutlinePolyData(renderer); float binaryOutlineWidth(1.0); if (this->GetDataNode()->GetFloatProperty( "outline width", binaryOutlineWidth, renderer )) { localStorage->m_Actor->GetProperty()->SetLineWidth(binaryOutlineWidth); } } else { binaryOutline = false; this->ApplyLookuptable(renderer); MITK_WARN << "Type of all binary images should be (un)signed char. Outline does not work on other pixel types!"; } } else //standard binary image { if(numberOfComponents != 1) { MITK_ERROR << "Rendering Error: Binary Images with more then 1 component are not supported!"; } } this->ApplyLookuptable(renderer); //Interpret the values as binary values localStorage->m_Texture->MapColorScalarsThroughLookupTableOn(); } else if( numberOfComponents == 1 ) //gray images { //Interpret the values as gray values localStorage->m_Texture->MapColorScalarsThroughLookupTableOn(); this->ApplyLookuptable(renderer); } else if ( (numberOfComponents == 3) || (numberOfComponents == 4) ) //RBG(A) images { //Interpret the RGB(A) images values correctly localStorage->m_Texture->MapColorScalarsThroughLookupTableOff(); this->ApplyLookuptable(renderer); this->ApplyRBGALevelWindow(renderer); } else { MITK_ERROR << "2D Reindering Error: Unknown number of components!!! Please report to rendering task force or check your data!"; } this->ApplyColor( renderer ); this->ApplyOpacity( renderer ); this->TransformActor( renderer ); //connect mapper with the data if(binary && binaryOutline) //connect the mapper with the polyData which contains the lines { //We need the contour for the binary oultine property as actor localStorage->m_Mapper->SetInput(localStorage->m_OutlinePolyData); localStorage->m_Actor->SetTexture(NULL); //no texture for contours } else { //Connect the mapper with the input texture. This is the standard case. //setup the textured plane this->GeneratePlane( renderer, sliceBounds ); //set the plane as input for the mapper localStorage->m_Mapper->SetInputConnection(localStorage->m_Plane->GetOutputPort()); //set the texture for the actor localStorage->m_Actor->SetTexture(localStorage->m_Texture); } // We have been modified => save this for next Update() localStorage->m_LastUpdateTime.Modified(); } void mitk::ImageVtkMapper2D::ApplyColor( mitk::BaseRenderer* renderer ) { LocalStorage *localStorage = this->GetLocalStorage( renderer ); // check for interpolation properties bool textureInterpolation = false; GetDataNode()->GetBoolProperty( "texture interpolation", textureInterpolation, renderer ); //set the interpolation modus according to the property localStorage->m_Texture->SetInterpolate(textureInterpolation); bool useColor = true; this->GetDataNode()->GetBoolProperty( "use color", useColor, renderer ); if( useColor ) { float rgb[3]= { 1.0f, 1.0f, 1.0f }; // check for color prop and use it for rendering if it exists // binary image hovering & binary image selection bool hover = false; bool selected = false; GetDataNode()->GetBoolProperty("binaryimage.ishovering", hover, renderer); GetDataNode()->GetBoolProperty("selected", selected, renderer); if(hover && !selected) { mitk::ColorProperty::Pointer colorprop = dynamic_cast(GetDataNode()->GetProperty - ("binaryimage.hoveringcolor", renderer)); + ("binaryimage.hoveringcolor", renderer)); if(colorprop.IsNotNull()) { memcpy(rgb, colorprop->GetColor().GetDataPointer(), 3*sizeof(float)); } else { GetColor( rgb, renderer ); } } if(selected) { mitk::ColorProperty::Pointer colorprop = dynamic_cast(GetDataNode()->GetProperty - ("binaryimage.selectedcolor", renderer)); + ("binaryimage.selectedcolor", renderer)); if(colorprop.IsNotNull()) { memcpy(rgb, colorprop->GetColor().GetDataPointer(), 3*sizeof(float)); } else { GetColor( rgb, renderer ); } } if(!hover && !selected) { GetColor( rgb, renderer ); } double rgbConv[3] = {(double)rgb[0], (double)rgb[1], (double)rgb[2]}; //conversion to double for VTK localStorage->m_Actor->GetProperty()->SetColor(rgbConv); } else { //If the user defines a lut, we dont want to use the color and take white instead. localStorage->m_Actor->GetProperty()->SetColor(1.0, 1.0, 1.0); } } void mitk::ImageVtkMapper2D::ApplyOpacity( mitk::BaseRenderer* renderer ) { LocalStorage* localStorage = this->GetLocalStorage( renderer ); float opacity = 1.0f; // check for opacity prop and use it for rendering if it exists GetOpacity( opacity, renderer ); //set the opacity according to the properties localStorage->m_Actor->GetProperty()->SetOpacity(opacity); } void mitk::ImageVtkMapper2D::ApplyLookuptable( mitk::BaseRenderer* renderer ) { bool binary = false; bool CTFcanBeApplied = false; this->GetDataNode()->GetBoolProperty( "binary", binary, renderer ); LocalStorage* localStorage = this->GetLocalStorage(renderer); //default lookuptable localStorage->m_Texture->SetLookupTable( localStorage->m_LookupTable ); if(binary) { //default lookuptable for binary images localStorage->m_Texture->GetLookupTable()->SetRange(0.0, 1.0); } else { bool useColor = true; this->GetDataNode()->GetBoolProperty( "use color", useColor, renderer ); if((!useColor)) { //BEGIN PROPERTY user-defined lut //currently we do not allow a lookuptable if it is a binary image // If lookup table use is requested... mitk::LookupTableProperty::Pointer LookupTableProp; LookupTableProp = dynamic_cast - (this->GetDataNode()->GetProperty("LookupTable")); + (this->GetDataNode()->GetProperty("LookupTable")); //...check if there is a lookuptable provided by the user if ( LookupTableProp.IsNotNull() ) { // If lookup table use is requested and supplied by the user: // only update the lut, when the properties have changed... if( LookupTableProp->GetLookupTable()->GetMTime() <= this->GetDataNode()->GetPropertyList()->GetMTime() ) - { + { LookupTableProp->GetLookupTable()->ChangeOpacityForAll( LookupTableProp->GetLookupTable()->GetVtkLookupTable()->GetAlpha()*localStorage->m_Actor->GetProperty()->GetOpacity() ); LookupTableProp->GetLookupTable()->ChangeOpacity(0, 0.0); } //we use the user-defined lookuptable localStorage->m_Texture->SetLookupTable( LookupTableProp->GetLookupTable()->GetVtkLookupTable() ); } else { CTFcanBeApplied = true; } }//END PROPERTY user-defined lut LevelWindow levelWindow; this->GetLevelWindow( levelWindow, renderer ); //set up the lookuptable with the level window range localStorage->m_Texture->GetLookupTable()->SetRange( levelWindow.GetLowerWindowBound(), levelWindow.GetUpperWindowBound() ); } //the color function can be applied if the user does not want to use color //and does not provide a lookuptable if(CTFcanBeApplied) { ApplyColorTransferFunction(renderer); } localStorage->m_Texture->SetInput( localStorage->m_ReslicedImage ); } void mitk::ImageVtkMapper2D::ApplyColorTransferFunction(mitk::BaseRenderer* renderer) { mitk::TransferFunctionProperty::Pointer transferFunctionProperty = - dynamic_cast(this->GetDataNode()->GetProperty("Image Rendering.Transfer Function",renderer )); + dynamic_cast(this->GetDataNode()->GetProperty("Image Rendering.Transfer Function",renderer )); LocalStorage* localStorage = m_LSH.GetLocalStorage(renderer); if(transferFunctionProperty.IsNotNull()) { localStorage->m_Texture->SetLookupTable(transferFunctionProperty->GetValue()->GetColorTransferFunction()); } else { MITK_WARN << "Neither a lookuptable nor a transfer function is set and use color is off."; } } -bool mitk::ImageVtkMapper2D::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; -} - -bool mitk::ImageVtkMapper2D::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 mitk::BoundingBox *boundingBox = boundingGeometry->GetBoundingBox(); - - mitk::BoundingBox::PointType bbMin = boundingBox->GetMinimum(); - mitk::BoundingBox::PointType bbMax = boundingBox->GetMaximum(); - - vtkSmartPointer points = vtkSmartPointer::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] ); - } - - vtkSmartPointer newPoints = vtkSmartPointer::New(); - - vtkSmartPointer transform = vtkSmartPointer::New(); - transform->Identity(); - transform->Concatenate( planeGeometry->GetVtkTransform()->GetLinearInverse() ); - - transform->Concatenate( boundingGeometry->GetVtkTransform() ); - - transform->TransformPoints( points, newPoints ); - 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 ); - - 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; - } -} void mitk::ImageVtkMapper2D::ApplyRBGALevelWindow( mitk::BaseRenderer* renderer ) { LocalStorage* localStorage = this->GetLocalStorage( renderer ); //pass the LuT to the RBG filter localStorage->m_LevelWindowToRGBFilterObject->SetLookupTable(localStorage->m_Texture->GetLookupTable()); mitk::LevelWindow opacLevelWindow; if( this->GetLevelWindow( opacLevelWindow, renderer, "opaclevelwindow" ) ) {//pass the opaque level window to the filter localStorage->m_LevelWindowToRGBFilterObject->SetMinOpacity(opacLevelWindow.GetLowerWindowBound()); localStorage->m_LevelWindowToRGBFilterObject->SetMaxOpacity(opacLevelWindow.GetUpperWindowBound()); } else {//no opaque level window localStorage->m_LevelWindowToRGBFilterObject->SetMinOpacity(0.0); localStorage->m_LevelWindowToRGBFilterObject->SetMaxOpacity(255.0); } localStorage->m_LevelWindowToRGBFilterObject->SetInput(localStorage->m_ReslicedImage); //connect the texture with the output of the RGB filter localStorage->m_Texture->SetInputConnection(localStorage->m_LevelWindowToRGBFilterObject->GetOutputPort()); } void mitk::ImageVtkMapper2D::Update(mitk::BaseRenderer* renderer) { if ( !this->IsVisible( renderer ) ) { return; } mitk::Image* data = const_cast( this->GetInput() ); if ( data == NULL ) { return; } // Calculate time step of the input data for the specified renderer (integer value) this->CalculateTimeStep( renderer ); // Check if time step is valid const TimeSlicedGeometry *dataTimeGeometry = data->GetTimeSlicedGeometry(); if ( ( dataTimeGeometry == NULL ) || ( dataTimeGeometry->GetTimeSteps() == 0 ) || ( !dataTimeGeometry->IsValidTime( this->GetTimestep() ) ) ) - { + { return; } const DataNode *node = this->GetDataNode(); data->UpdateOutputInformation(); LocalStorage *localStorage = m_LSH.GetLocalStorage(renderer); //check if something important has changed and we need to rerender if ( (localStorage->m_LastUpdateTime < node->GetMTime()) //was the node modified? || (localStorage->m_LastUpdateTime < data->GetPipelineMTime()) //Was the data modified? || (localStorage->m_LastUpdateTime < renderer->GetCurrentWorldGeometry2DUpdateTime()) //was the geometry modified? || (localStorage->m_LastUpdateTime < renderer->GetCurrentWorldGeometry2D()->GetMTime()) || (localStorage->m_LastUpdateTime < node->GetPropertyList()->GetMTime()) //was a property modified? || (localStorage->m_LastUpdateTime < node->GetPropertyList(renderer)->GetMTime()) ) - { + { this->GenerateDataForRenderer( renderer ); } // since we have checked that nothing important has changed, we can set // m_LastUpdateTime to the current time localStorage->m_LastUpdateTime.Modified(); } void mitk::ImageVtkMapper2D::SetDefaultProperties(mitk::DataNode* node, mitk::BaseRenderer* renderer, bool overwrite) { mitk::Image::Pointer image = dynamic_cast(node->GetData()); // Properties common for both images and segmentations node->AddProperty( "use color", mitk::BoolProperty::New( true ), renderer, overwrite ); node->AddProperty( "depthOffset", mitk::FloatProperty::New( 0.0 ), renderer, overwrite ); node->AddProperty( "outline binary", mitk::BoolProperty::New( false ), renderer, overwrite ); node->AddProperty( "outline width", mitk::FloatProperty::New( 1.0 ), renderer, overwrite ); if(image->IsRotated()) node->AddProperty( "reslice interpolation", mitk::VtkResliceInterpolationProperty::New(VTK_RESLICE_CUBIC) ); else node->AddProperty( "reslice interpolation", mitk::VtkResliceInterpolationProperty::New() ); node->AddProperty( "texture interpolation", mitk::BoolProperty::New( mitk::DataNodeFactory::m_TextureInterpolationActive ) ); // set to user configurable default value (see global options) node->AddProperty( "in plane resample extent by geometry", mitk::BoolProperty::New( false ) ); node->AddProperty( "bounding box", mitk::BoolProperty::New( false ) ); std::string modality; if ( node->GetStringProperty( "dicom.series.Modality", modality ) ) { // modality provided by DICOM or other reader if ( modality == "PT") // NOT a typo, PT is the abbreviation for PET used in DICOM { node->SetProperty( "use color", mitk::BoolProperty::New( false ), renderer ); node->SetProperty( "opacity", mitk::FloatProperty::New( 0.5 ), renderer ); } } bool isBinaryImage(false); if ( ! node->GetBoolProperty("binary", isBinaryImage) ) { // ok, property is not set, use heuristic to determine if this // is a binary image mitk::Image::Pointer centralSliceImage; ScalarType minValue = 0.0; ScalarType maxValue = 0.0; ScalarType min2ndValue = 0.0; ScalarType max2ndValue = 0.0; mitk::ImageSliceSelector::Pointer sliceSelector = mitk::ImageSliceSelector::New(); sliceSelector->SetInput(image); sliceSelector->SetSliceNr(image->GetDimension(2)/2); sliceSelector->SetTimeNr(image->GetDimension(3)/2); sliceSelector->SetChannelNr(image->GetDimension(4)/2); sliceSelector->Update(); centralSliceImage = sliceSelector->GetOutput(); if ( centralSliceImage.IsNotNull() && centralSliceImage->IsInitialized() ) { minValue = centralSliceImage->GetStatistics()->GetScalarValueMin(); maxValue = centralSliceImage->GetStatistics()->GetScalarValueMax(); min2ndValue = centralSliceImage->GetStatistics()->GetScalarValue2ndMin(); max2ndValue = centralSliceImage->GetStatistics()->GetScalarValue2ndMax(); } if ( minValue == maxValue ) { // centralSlice is strange, lets look at all data minValue = image->GetStatistics()->GetScalarValueMin(); maxValue = image->GetStatistics()->GetScalarValueMaxNoRecompute(); min2ndValue = image->GetStatistics()->GetScalarValue2ndMinNoRecompute(); max2ndValue = image->GetStatistics()->GetScalarValue2ndMaxNoRecompute(); } isBinaryImage = ( maxValue == min2ndValue && minValue == max2ndValue ); } // some more properties specific for a binary... if (isBinaryImage) { node->AddProperty( "opacity", mitk::FloatProperty::New(0.3f), renderer, overwrite ); node->AddProperty( "color", ColorProperty::New(1.0,0.0,0.0), renderer, overwrite ); node->AddProperty( "binaryimage.selectedcolor", ColorProperty::New(1.0,0.0,0.0), renderer, overwrite ); node->AddProperty( "binaryimage.selectedannotationcolor", ColorProperty::New(1.0,0.0,0.0), renderer, overwrite ); node->AddProperty( "binaryimage.hoveringcolor", ColorProperty::New(1.0,0.0,0.0), renderer, overwrite ); node->AddProperty( "binaryimage.hoveringannotationcolor", ColorProperty::New(1.0,0.0,0.0), renderer, overwrite ); node->AddProperty( "binary", mitk::BoolProperty::New( true ), renderer, overwrite ); node->AddProperty("layer", mitk::IntProperty::New(10), renderer, overwrite); } else //...or image type object { node->AddProperty( "opacity", mitk::FloatProperty::New(1.0f), renderer, overwrite ); node->AddProperty( "color", ColorProperty::New(1.0,1.0,1.0), renderer, overwrite ); node->AddProperty( "binary", mitk::BoolProperty::New( false ), renderer, overwrite ); node->AddProperty("layer", mitk::IntProperty::New(0), renderer, overwrite); } if(image.IsNotNull() && image->IsInitialized()) { if((overwrite) || (node->GetProperty("levelwindow", renderer)==NULL)) { /* initialize level/window from DICOM tags */ std::string sLevel; std::string sWindow; if ( node->GetStringProperty( "dicom.voilut.WindowCenter", sLevel ) && node->GetStringProperty( "dicom.voilut.WindowWidth", sWindow ) ) - { + { float level = atof( sLevel.c_str() ); float window = atof( sWindow.c_str() ); mitk::LevelWindow contrast; std::string sSmallestPixelValueInSeries; std::string sLargestPixelValueInSeries; if ( node->GetStringProperty( "dicom.series.SmallestPixelValueInSeries", sSmallestPixelValueInSeries ) && node->GetStringProperty( "dicom.series.LargestPixelValueInSeries", sLargestPixelValueInSeries ) ) - { + { float smallestPixelValueInSeries = atof( sSmallestPixelValueInSeries.c_str() ); float largestPixelValueInSeries = atof( sLargestPixelValueInSeries.c_str() ); contrast.SetRangeMinMax( smallestPixelValueInSeries-1, largestPixelValueInSeries+1 ); // why not a little buffer? // might remedy some l/w widget challenges } else { contrast.SetAuto( static_cast(node->GetData()), false, true ); // we need this as a fallback } contrast.SetLevelWindow( level, window); node->SetProperty( "levelwindow", LevelWindowProperty::New( contrast ), renderer ); } } if(((overwrite) || (node->GetProperty("opaclevelwindow", renderer)==NULL)) && image->GetPixelType().GetPixelTypeId() == typeid(itk::RGBAPixel)) - { + { mitk::LevelWindow opaclevwin; opaclevwin.SetRangeMinMax(0,255); opaclevwin.SetWindowBounds(0,255); mitk::LevelWindowProperty::Pointer prop = mitk::LevelWindowProperty::New(opaclevwin); node->SetProperty( "opaclevelwindow", prop, renderer ); } if((overwrite) || (node->GetProperty("LookupTable", renderer)==NULL)) { // add a default rainbow lookup table for color mapping mitk::LookupTable::Pointer mitkLut = mitk::LookupTable::New(); vtkLookupTable* vtkLut = mitkLut->GetVtkLookupTable(); vtkLut->SetHueRange(0.6667, 0.0); vtkLut->SetTableRange(0.0, 20.0); vtkLut->Build(); mitk::LookupTableProperty::Pointer mitkLutProp = mitk::LookupTableProperty::New(); mitkLutProp->SetLookupTable(mitkLut); node->SetProperty( "LookupTable", mitkLutProp ); } } Superclass::SetDefaultProperties(node, renderer, overwrite); } mitk::ImageVtkMapper2D::LocalStorage* mitk::ImageVtkMapper2D::GetLocalStorage(mitk::BaseRenderer* renderer) { return m_LSH.GetLocalStorage(renderer); } vtkSmartPointer mitk::ImageVtkMapper2D::CreateOutlinePolyData(mitk::BaseRenderer* renderer ){ LocalStorage* localStorage = this->GetLocalStorage(renderer); //get the min and max index values of each direction int* extent = localStorage->m_ReslicedImage->GetExtent(); int xMin = extent[0]; int xMax = extent[1]; int yMin = extent[2]; int yMax = extent[3]; int* dims = localStorage->m_ReslicedImage->GetDimensions(); //dimensions of the image int line = dims[0]; //how many pixels per line? int x = xMin; //pixel index x int y = yMin; //pixel index y char* currentPixel; - + //get the depth for each contour float depth = CalculateLayerDepth(renderer); vtkSmartPointer points = vtkSmartPointer::New(); //the points to draw vtkSmartPointer lines = vtkSmartPointer::New(); //the lines to connect the points while (y <= yMax) { currentPixel = static_cast(localStorage->m_ReslicedImage->GetScalarPointer(x, y, 0)); //if the current pixel value is set to something if ((currentPixel) && (*currentPixel != 0)) - { + { //check in which direction a line is necessary - //a line is added if the neighbor of the current pixel has the value 0 - //and if the pixel is located at the edge of the image + //a line is added if the neighbor of the current pixel has the value 0 + //and if the pixel is located at the edge of the image - //if vvvvv not the first line vvvvv + //if vvvvv not the first line vvvvv if (y > yMin && *(currentPixel-line) == 0) - { //x direction - bottom edge of the pixel + { //x direction - bottom edge of the pixel //add the 2 points vtkIdType p1 = points->InsertNextPoint(x*localStorage->m_mmPerPixel[0], y*localStorage->m_mmPerPixel[1], depth); vtkIdType p2 = points->InsertNextPoint((x+1)*localStorage->m_mmPerPixel[0], y*localStorage->m_mmPerPixel[1], depth); //add the line between both points lines->InsertNextCell(2); lines->InsertCellPoint(p1); lines->InsertCellPoint(p2); } - //if vvvvv not the last line vvvvv + //if vvvvv not the last line vvvvv if (y < yMax && *(currentPixel+line) == 0) - { //x direction - top edge of the pixel + { //x direction - top edge of the pixel vtkIdType p1 = points->InsertNextPoint(x*localStorage->m_mmPerPixel[0], (y+1)*localStorage->m_mmPerPixel[1], depth); vtkIdType p2 = points->InsertNextPoint((x+1)*localStorage->m_mmPerPixel[0], (y+1)*localStorage->m_mmPerPixel[1], depth); lines->InsertNextCell(2); lines->InsertCellPoint(p1); lines->InsertCellPoint(p2); } - //if vvvvv not the first pixel vvvvv + //if vvvvv not the first pixel vvvvv if ( (x > xMin || y > yMin) && *(currentPixel-1) == 0) - { //y direction - left edge of the pixel + { //y direction - left edge of the pixel vtkIdType p1 = points->InsertNextPoint(x*localStorage->m_mmPerPixel[0], y*localStorage->m_mmPerPixel[1], depth); vtkIdType p2 = points->InsertNextPoint(x*localStorage->m_mmPerPixel[0], (y+1)*localStorage->m_mmPerPixel[1], depth); lines->InsertNextCell(2); lines->InsertCellPoint(p1); lines->InsertCellPoint(p2); } - //if vvvvv not the last pixel vvvvv + //if vvvvv not the last pixel vvvvv if ( (y < yMax || (x < xMax) ) && *(currentPixel+1) == 0) - { //y direction - right edge of the pixel + { //y direction - right edge of the pixel vtkIdType p1 = points->InsertNextPoint((x+1)*localStorage->m_mmPerPixel[0], y*localStorage->m_mmPerPixel[1], depth); vtkIdType p2 = points->InsertNextPoint((x+1)*localStorage->m_mmPerPixel[0], (y+1)*localStorage->m_mmPerPixel[1], depth); lines->InsertNextCell(2); lines->InsertCellPoint(p1); lines->InsertCellPoint(p2); } - /* now consider pixels at the edge of the image */ + /* now consider pixels at the edge of the image */ - //if vvvvv left edge of image vvvvv + //if vvvvv left edge of image vvvvv if (x == xMin) - { //draw left edge of the pixel + { //draw left edge of the pixel vtkIdType p1 = points->InsertNextPoint(x*localStorage->m_mmPerPixel[0], y*localStorage->m_mmPerPixel[1], depth); vtkIdType p2 = points->InsertNextPoint(x*localStorage->m_mmPerPixel[0], (y+1)*localStorage->m_mmPerPixel[1], depth); lines->InsertNextCell(2); lines->InsertCellPoint(p1); lines->InsertCellPoint(p2); } - //if vvvvv right edge of image vvvvv + //if vvvvv right edge of image vvvvv if (x == xMax) - { //draw right edge of the pixel + { //draw right edge of the pixel vtkIdType p1 = points->InsertNextPoint((x+1)*localStorage->m_mmPerPixel[0], y*localStorage->m_mmPerPixel[1], depth); vtkIdType p2 = points->InsertNextPoint((x+1)*localStorage->m_mmPerPixel[0], (y+1)*localStorage->m_mmPerPixel[1], depth); lines->InsertNextCell(2); lines->InsertCellPoint(p1); lines->InsertCellPoint(p2); } - //if vvvvv bottom edge of image vvvvv + //if vvvvv bottom edge of image vvvvv if (y == yMin) - { //draw bottom edge of the pixel + { //draw bottom edge of the pixel vtkIdType p1 = points->InsertNextPoint(x*localStorage->m_mmPerPixel[0], y*localStorage->m_mmPerPixel[1], depth); vtkIdType p2 = points->InsertNextPoint((x+1)*localStorage->m_mmPerPixel[0], y*localStorage->m_mmPerPixel[1], depth); lines->InsertNextCell(2); lines->InsertCellPoint(p1); lines->InsertCellPoint(p2); } - //if vvvvv top edge of image vvvvv + //if vvvvv top edge of image vvvvv if (y == yMax) - { //draw top edge of the pixel + { //draw top edge of the pixel vtkIdType p1 = points->InsertNextPoint(x*localStorage->m_mmPerPixel[0], (y+1)*localStorage->m_mmPerPixel[1], depth); vtkIdType p2 = points->InsertNextPoint((x+1)*localStorage->m_mmPerPixel[0], (y+1)*localStorage->m_mmPerPixel[1], depth); lines->InsertNextCell(2); lines->InsertCellPoint(p1); lines->InsertCellPoint(p2); } }//end if currentpixel is set x++; if (x > xMax) - { //reached end of line + { //reached end of line x = xMin; y++; } }//end of while // Create a polydata to store everything in vtkSmartPointer polyData = vtkSmartPointer::New(); // Add the points to the dataset polyData->SetPoints(points); // Add the lines to the dataset polyData->SetLines(lines); return polyData; } void mitk::ImageVtkMapper2D::TransformActor(mitk::BaseRenderer* renderer) { LocalStorage *localStorage = m_LSH.GetLocalStorage(renderer); //get the transformation matrix of the reslicer in order to render the slice as transversal, coronal or saggital vtkSmartPointer trans = vtkSmartPointer::New(); vtkSmartPointer matrix = localStorage->m_Reslicer->GetResliceAxes(); trans->SetMatrix(matrix); //transform the plane/contour (the actual actor) to the corresponding view (transversal, coronal or saggital) localStorage->m_Actor->SetUserTransform(trans); //transform the origin to center based coordinates, because MITK is center based. localStorage->m_Actor->SetPosition( -0.5*localStorage->m_mmPerPixel[0], -0.5*localStorage->m_mmPerPixel[1], 0.0); } mitk::ImageVtkMapper2D::LocalStorage::LocalStorage() { //Do as much actions as possible in here to avoid double executions. m_Plane = vtkSmartPointer::New(); m_Texture = vtkSmartPointer::New(); m_LookupTable = vtkSmartPointer::New(); m_Mapper = vtkSmartPointer::New(); m_Actor = vtkSmartPointer::New(); - m_Reslicer = vtkSmartPointer::New(); + m_Reslicer = mitk::ExtractSliceFilter::New(); m_TSFilter = vtkSmartPointer::New(); - m_UnitSpacingImageFilter = vtkSmartPointer::New(); m_OutlinePolyData = vtkSmartPointer::New(); m_ReslicedImage = vtkSmartPointer::New(); //the following actions are always the same and thus can be performed //in the constructor for each image (i.e. the image-corresponding local storage) m_TSFilter->ReleaseDataFlagOn(); - m_Reslicer->ReleaseDataFlagOn(); - - m_UnitSpacingImageFilter->SetOutputSpacing( 1.0, 1.0, 1.0 ); //built a default lookuptable m_LookupTable->SetRampToLinear(); m_LookupTable->SetSaturationRange( 0.0, 0.0 ); m_LookupTable->SetHueRange( 0.0, 0.0 ); m_LookupTable->SetValueRange( 0.0, 1.0 ); m_LookupTable->Build(); //map all black values to transparent m_LookupTable->SetTableValue(0, 0.0, 0.0, 0.0, 0.0); //do not repeat the texture (the image) m_Texture->RepeatOff(); //set the mapper for the actor m_Actor->SetMapper(m_Mapper); //filter for RGB(A) images m_LevelWindowToRGBFilterObject = new vtkMitkApplyLevelWindowToRGBFilter(); } diff --git a/Core/Code/Rendering/mitkImageVtkMapper2D.h b/Core/Code/Rendering/mitkImageVtkMapper2D.h index 2c10a08a21..c3ac358103 100644 --- a/Core/Code/Rendering/mitkImageVtkMapper2D.h +++ b/Core/Code/Rendering/mitkImageVtkMapper2D.h @@ -1,282 +1,265 @@ /*=================================================================== 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 MITKIMAGEVTKMAPPER2D_H_HEADER_INCLUDED_C10E906E #define MITKIMAGEVTKMAPPER2D_H_HEADER_INCLUDED_C10E906E //MITK #include //MITK Rendering #include "mitkBaseRenderer.h" #include "mitkVtkMapper2D.h" +#include "mitkExtractSliceFilter.h" //VTK #include class vtkActor; class vtkPolyDataMapper; class vtkPlaneSource; class vtkImageData; class vtkLookupTable; class vtkImageReslice; class vtkImageChangeInformation; class vtkPoints; class vtkMitkThickSlicesFilter; class vtkPolyData; class vtkMitkApplyLevelWindowToRGBFilter; namespace mitk { /** \brief Mapper to resample and display 2D slices of a 3D image. * * The following image gives a brief overview of the mapping and the involved parts. * * \image html imageVtkMapper2Darchitecture.png * * First, the image is resliced by means of vtkImageReslice. The volume image * serves as input to the mapper in addition to spatial placement of the slice and a few other * properties such as thick slices. This code was already present in the old version * (mitkImageMapperGL2D). * * Next, the obtained slice (m_ReslicedImage) is used to create a texture * (m_Texture) and a plane onto which the texture is rendered (m_Plane). For * mapping purposes, a vtkPolyDataMapper (m_Mapper) is utilized. Orthographic * projection is applied to create the effect of a 2D image. The mapper and the * texture are assigned to the actor (m_Actor) which is passed to the VTK rendering * pipeline via the method GetVtkProp(). * * In order to transform the textured plane to the correct position in space, the * same transformation as used for reslicing is applied to both the camera and the * vtkActor. All important steps are explained in more detail below. The resulting * 2D image (by reslicing the underlying 3D input image appropriately) can either * be directly rendered in a 2D view or just be calculated to be used later by another * rendering entity, e.g. in texture mapping in a 3D view. * * Properties that can be set for images and influence the imageMapper2D are: * * - \b "opacity": (FloatProperty) Opacity of the image * - \b "color": (ColorProperty) Color of the image * - \b "use color": (BoolProperty) Use the color of the image or not * - \b "binary": (BoolProperty) is the image a binary image or not * - \b "outline binary": (BoolProperty) show outline of the image or not * - \b "texture interpolation": (BoolProperty) texture interpolation of the image * - \b "reslice interpolation": (VtkResliceInterpolationProperty) reslice interpolation of the image * - \b "in plane resample extent by geometry": (BoolProperty) Do it or not * - \b "bounding box": (BoolProperty) Is the Bounding Box of the image shown or not * - \b "layer": (IntProperty) Layer of the image * - \b "volume annotation color": (ColorProperty) color of the volume annotation, TODO has to be reimplemented * - \b "volume annotation unit": (StringProperty) annotation unit as string (does not implicit convert the unit!) unit is ml or cm3, TODO has to be reimplemented * The default properties are: * - \b "opacity", mitk::FloatProperty::New(0.3f), renderer, overwrite ) * - \b "color", ColorProperty::New(1.0,0.0,0.0), renderer, overwrite ) * - \b "use color", mitk::BoolProperty::New( true ), renderer, overwrite ) * - \b "binary", mitk::BoolProperty::New( true ), renderer, overwrite ) * - \b "outline binary", mitk::BoolProperty::New( false ), renderer, overwrite ) * - \b "texture interpolation", mitk::BoolProperty::New( mitk::DataNodeFactory::m_TextureInterpolationActive ) ) * - \b "reslice interpolation", mitk::VtkResliceInterpolationProperty::New() ) * - \b "in plane resample extent by geometry", mitk::BoolProperty::New( false ) ) * - \b "bounding box", mitk::BoolProperty::New( false ) ) * - \b "layer", mitk::IntProperty::New(10), renderer, overwrite) * If the modality-property is set for an image, the mapper uses modality-specific default properties, * e.g. color maps, if they are defined. * \ingroup Mapper */ class MITK_CORE_EXPORT ImageVtkMapper2D : public VtkMapper2D { public: /** Standard class typedefs. */ mitkClassMacro( ImageVtkMapper2D,VtkMapper2D ); /** Method for creation through the object factory. */ itkNewMacro(Self); /** \brief Get the Image to map */ const mitk::Image *GetInput(void); /** \brief Checks whether this mapper needs to update itself and generate * data. */ virtual void Update(mitk::BaseRenderer * renderer); //### methods of MITK-VTK rendering pipeline virtual vtkProp* GetVtkProp(mitk::BaseRenderer* renderer); virtual void MitkRenderOverlay(BaseRenderer* renderer); virtual void MitkRenderOpaqueGeometry(BaseRenderer* renderer); virtual void MitkRenderTranslucentGeometry(BaseRenderer* renderer); virtual void MitkRenderVolumetricGeometry(BaseRenderer* renderer); //### end of methods of MITK-VTK rendering pipeline /** \brief Internal class holding the mapper, actor, etc. for each of the 3 2D render windows */ /** * To render transveral, coronal, and sagittal, the mapper is called three times. * For performance reasons, the corresponding data for each view is saved in the * internal helper class LocalStorage. This allows rendering n views with just * 1 mitkMapper using n vtkMapper. * */ class MITK_CORE_EXPORT LocalStorage : public mitk::Mapper::BaseLocalStorage { public: /** \brief Actor of a 2D render window. */ vtkSmartPointer m_Actor; /** \brief Mapper of a 2D render window. */ vtkSmartPointer m_Mapper; /** \brief Current slice of a 2D render window. */ vtkSmartPointer m_ReslicedImage; /** \brief Plane on which the slice is rendered as texture. */ vtkSmartPointer m_Plane; /** \brief The texture which is used to render the current slice. */ vtkSmartPointer m_Texture; /** \brief The lookuptable for colors and level window */ vtkSmartPointer m_LookupTable; /** \brief The actual reslicer (one per renderer) */ - vtkSmartPointer m_Reslicer; - /** \brief Thickslices post filtering. */ - vtkSmartPointer m_TSFilter; - /** \brief Using unit spacing for resampling makes life easier TODO improve docu ...*/ - vtkSmartPointer m_UnitSpacingImageFilter; + mitk::ExtractSliceFilter::Pointer m_Reslicer; + /** \brief Filter for thick slices */ + vtkSmartPointer m_TSFilter; /** \brief PolyData object containg all lines/points needed for outlining the contour. This container is used to save a computed contour for the next rendering execution. For instance, if you zoom or pann, there is no need to recompute the contour. */ vtkSmartPointer m_OutlinePolyData; + /** \brief Timestamp of last update of stored data. */ itk::TimeStamp m_LastUpdateTime; - /** \brief Origin of the 2D geometry. */ - mitk::Point3D m_Origin; - /** \brief Bottom end point of the y-axis of the 2D geometry. */ - mitk::Vector3D m_Bottom; - /** \brief Right end point of the x-axis of the 2D geometry. */ - mitk::Vector3D m_Right; - /** \brief Normal of the 2D geometry. */ - mitk::Vector3D m_Normal; /** \brief mmPerPixel relation between pixel and mm. (World spacing).*/ - mitk::ScalarType m_mmPerPixel[2]; + mitk::ScalarType* m_mmPerPixel; /** \brief This filter is used to apply the level window to RBG(A) images. */ vtkMitkApplyLevelWindowToRGBFilter* m_LevelWindowToRGBFilterObject; /** \brief Default constructor of the local storage. */ LocalStorage(); /** \brief Default deconstructor of the local storage. */ ~LocalStorage() { } }; /** \brief The LocalStorageHandler holds all (three) LocalStorages for the three 2D render windows. */ mitk::Mapper::LocalStorageHandler m_LSH; /** \brief Get the LocalStorage corresponding to the current renderer. */ LocalStorage* GetLocalStorage(mitk::BaseRenderer* renderer); /** \brief Set the default properties for general image rendering. */ static void SetDefaultProperties(mitk::DataNode* node, mitk::BaseRenderer* renderer = NULL, bool overwrite = false); protected: /** \brief Transforms the actor to the actual position in 3D. * \param renderer The current renderer corresponding to the render window. */ void TransformActor(mitk::BaseRenderer* renderer); /** \brief Generates a plane according to the size of the resliced image in milimeters. * * \image html texturedPlane.png * * In VTK a vtkPlaneSource is defined through three points. The origin and two * points defining the axes of the plane (see VTK documentation). The origin is * set to (xMin; yMin; Z), where xMin and yMin are the minimal bounds of the * resliced image in space. Z is relevant for blending and the layer property. * The center of the plane (C) is also the center of the view plane (cf. the image above). * * \note For the standard MITK view with three 2D render windows showing three * different slices, three such planes are generated. All these planes are generated * in the XY-plane (even if they depict a YZ-slice of the volume). * */ void GeneratePlane(mitk::BaseRenderer* renderer, vtkFloatingPointType planeBounds[6]); /** \brief Generates a vtkPolyData object containing the outline of a given binary slice. \param binarySlice - The binary image slice. (Volumes are not supported.) \param mmPerPixel - Spacing of the binary image slice. Hence it's 2D, only in x/y-direction. \note This code has been taken from the deprecated library iil. */ vtkSmartPointer CreateOutlinePolyData(mitk::BaseRenderer* renderer); /** Default constructor */ ImageVtkMapper2D(); /** Default deconstructor */ virtual ~ImageVtkMapper2D(); /** \brief Does the actual resampling, without rendering the image yet. * All the data is generated inside this method. The vtkProp (or Actor) * is filled with content (i.e. the resliced image). * * After generation, a 4x4 transformation matrix(t) of the current slice is obtained * from the vtkResliceImage object via GetReslicesAxis(). This matrix is * applied to each textured plane (actor->SetUserTransform(t)) to transform everything * to the actual 3D position (cf. the following image). * * \image html cameraPositioning3D.png * */ virtual void GenerateDataForRenderer(mitk::BaseRenderer *renderer); - /** \brief Internal helper method for intersection testing used only in CalculateClippedPlaneBounds() */ - bool LineIntersectZero( vtkPoints *points, int p1, int p2, - vtkFloatingPointType *bounds ); - - /** \brief Calculate the bounding box of the resliced image. This is necessary for - arbitrarily rotated planes in an image volume. A rotated plane (e.g. in swivel mode) - will have a new bounding box, which needs to be calculated. */ - bool CalculateClippedPlaneBounds( const Geometry3D *boundingGeometry, - const PlaneGeometry *planeGeometry, vtkFloatingPointType *bounds ); /** \brief This method uses the vtkCamera clipping range and the layer property * to calcualte the depth of the object (e.g. image or contour). The depth is used * to keep the correct order for the final VTK rendering.*/ float CalculateLayerDepth(mitk::BaseRenderer* renderer); /** \brief This method applies a level window on RBG(A) images. * It should only be called for internally for RGB(A) images. */ void ApplyRBGALevelWindow( mitk::BaseRenderer* renderer ); /** \brief This method applies (or modifies) the lookuptable for all types of images. */ void ApplyLookuptable( mitk::BaseRenderer* renderer ); /** \brief This method applies a color transfer function, if no LookuptableProperty is set. Internally, a vtkColorTransferFunction is used. This is usefull for coloring continous images (e.g. float) */ void ApplyColorTransferFunction(mitk::BaseRenderer* renderer); /** \brief Set the color of the image/polydata */ void ApplyColor( mitk::BaseRenderer* renderer ); /** \brief Set the opacity of the actor. */ void ApplyOpacity( mitk::BaseRenderer* renderer ); }; } // namespace mitk #endif /* MITKIMAGEVTKMAPPER2D_H_HEADER_INCLUDED_C10E906E */ diff --git a/Core/Code/Testing/files.cmake b/Core/Code/Testing/files.cmake index da2c1ec1ae..8f24b6f0bc 100644 --- a/Core/Code/Testing/files.cmake +++ b/Core/Code/Testing/files.cmake @@ -1,112 +1,113 @@ # tests with no extra command line parameter set(MODULE_TESTS mitkAccessByItkTest.cpp mitkCoreObjectFactoryTest.cpp mitkMaterialTest.cpp mitkActionTest.cpp mitkEnumerationPropertyTest.cpp mitkEventTest.cpp mitkFocusManagerTest.cpp mitkGenericPropertyTest.cpp mitkGeometry3DTest.cpp mitkGeometryDataToSurfaceFilterTest.cpp mitkGlobalInteractionTest.cpp mitkImageDataItemTest.cpp #mitkImageMapper2DTest.cpp mitkImageGeneratorTest.cpp mitkBaseDataTest.cpp #mitkImageToItkTest.cpp mitkInstantiateAccessFunctionTest.cpp mitkInteractorTest.cpp mitkITKThreadingTest.cpp # mitkLevelWindowManagerTest.cpp mitkLevelWindowTest.cpp mitkMessageTest.cpp #mitkPipelineSmartPointerCorrectnessTest.cpp mitkPixelTypeTest.cpp mitkPlaneGeometryTest.cpp mitkPointSetFileIOTest.cpp mitkPointSetTest.cpp mitkPointSetWriterTest.cpp mitkPointSetReaderTest.cpp mitkPointSetInteractorTest.cpp mitkPropertyTest.cpp mitkPropertyListTest.cpp #mitkRegistrationBaseTest.cpp #mitkSegmentationInterpolationTest.cpp mitkSlicedGeometry3DTest.cpp mitkSliceNavigationControllerTest.cpp mitkStateMachineTest.cpp mitkStateTest.cpp mitkSurfaceTest.cpp mitkSurfaceToSurfaceFilterTest.cpp mitkTimeSlicedGeometryTest.cpp mitkTransitionTest.cpp mitkUndoControllerTest.cpp mitkVtkWidgetRenderingTest.cpp mitkVerboseLimitedLinearUndoTest.cpp mitkWeakPointerTest.cpp mitkTransferFunctionTest.cpp #mitkAbstractTransformGeometryTest.cpp mitkStepperTest.cpp itkTotalVariationDenoisingImageFilterTest.cpp mitkRenderingManagerTest.cpp vtkMitkThickSlicesFilterTest.cpp mitkNodePredicateSourceTest.cpp mitkVectorTest.cpp mitkClippedSurfaceBoundsCalculatorTest.cpp #QmitkRenderingTestHelper.cpp mitkExceptionTest.cpp + mitkExtractSliceFilterTest.cpp ) # test with image filename as an extra command line parameter set(MODULE_IMAGE_TESTS mitkPlanePositionManagerTest.cpp mitkSurfaceVtkWriterTest.cpp #mitkImageSliceSelectorTest.cpp mitkImageTimeSelectorTest.cpp # mitkVtkPropRendererTest.cpp mitkDataNodeFactoryTest.cpp #mitkSTLFileReaderTest.cpp ) # list of images for which the tests are run set(MODULE_TESTIMAGES # Pic-Factory no more available in Core, test images now in .nrrd format US4DCyl.nrrd Pic3D.nrrd Pic2DplusT.nrrd BallBinary30x30x30.nrrd binary.stl ball.stl ) set(MODULE_CUSTOM_TESTS #mitkLabeledImageToSurfaceFilterTest.cpp #mitkExternalToolsTest.cpp mitkDataStorageTest.cpp mitkDataNodeTest.cpp mitkDicomSeriesReaderTest.cpp mitkDICOMLocaleTest.cpp mitkEventMapperTest.cpp mitkNodeDependentPointSetInteractorTest.cpp mitkStateMachineFactoryTest.cpp mitkPointSetLocaleTest.cpp mitkImageTest.cpp mitkImageWriterTest.cpp mitkImageVtkMapper2DTest.cpp mitkImageVtkMapper2DLevelWindowTest.cpp mitkImageVtkMapper2DOpacityTest.cpp mitkImageVtkMapper2DColorTest.cpp ) # Create an artificial module initializing class for # the usServiceListenerTest.cpp usFunctionGenerateModuleInit(testdriver_init_file NAME ${MODULE_NAME}TestDriver DEPENDS "Mitk" VERSION "0.1.0" EXECUTABLE ) set(TEST_CPP_FILES ${testdriver_init_file} mitkRenderingTestHelper.cpp) diff --git a/Core/Code/Testing/mitkExtractSliceFilterTest.cpp b/Core/Code/Testing/mitkExtractSliceFilterTest.cpp new file mode 100644 index 0000000000..0d6067d787 --- /dev/null +++ b/Core/Code/Testing/mitkExtractSliceFilterTest.cpp @@ -0,0 +1,1161 @@ +/*========================================================================= + +Program: Medical Imaging & Interaction Toolkit +Language: C++ +Date: $Date$ +Version: $Revision: 7837 $ + +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 +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include + + +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + + +//use this to create the test volume on the fly +#define CREATE_VOLUME + +//use this to save the created volume +//#define SAVE_VOLUME + +//use this to calculate the error from the sphere mathematical model to our pixel based one +//#define CALC_TESTFAILURE_DEVIATION + +//use this to render an oblique slice through a specified image +//#define SHOW_SLICE_IN_RENDER_WINDOW + +//use this to have infos printed in mbilog +//#define EXTRACTOR_DEBUG + + +/*these are the deviations calculated by the function CalcTestFailureDeviation (see for details)*/ +#define Testfailure_Deviation_Mean_128 0.853842 + +#define Testfailure_Deviation_Volume_128 0.145184 +#define Testfailure_Deviation_Diameter_128 1.5625 + + + +#define Testfailure_Deviation_Mean_256 0.397693 + +#define Testfailure_Deviation_Volume_256 0.0141357 +#define Testfailure_Deviation_Diameter_256 0.78125 + + + +#define Testfailure_Deviation_Mean_512 0.205277 + +#define Testfailure_Deviation_Volume_512 0.01993 +#define Testfailure_Deviation_Diameter_512 0.390625 + + + +class mitkExtractSliceFilterTestClass{ + +public: + + + static void TestSlice(mitk::PlaneGeometry* planeGeometry, std::string testname) + { + TestPlane = planeGeometry; + TestName = testname; + + float centerCoordValue = TestvolumeSize / 2.0; + float center[3] = {centerCoordValue, centerCoordValue, centerCoordValue}; + mitk::Point3D centerIndex(center); + + double radius = TestvolumeSize / 4.0; + if(TestPlane->Distance(centerIndex) >= radius ) return;//outside sphere + + + + //feed ExtractSliceFilter + mitk::ExtractSliceFilter::Pointer slicer = mitk::ExtractSliceFilter::New(); + + slicer->SetInput(TestVolume); + slicer->SetWorldGeometry(TestPlane); + slicer->Update(); + + MITK_TEST_CONDITION_REQUIRED(slicer->GetOutput() != NULL, "Extractor returned a slice"); + + mitk::Image::Pointer reslicedImage = slicer->GetOutput(); + + AccessFixedDimensionByItk(reslicedImage, TestSphereRadiusByItk, 2); + AccessFixedDimensionByItk(reslicedImage, TestSphereAreaByItk, 2); + + double devArea, devDiameter; + if(TestvolumeSize == 128.0){ devArea = Testfailure_Deviation_Volume_128; devDiameter = Testfailure_Deviation_Diameter_128; } + else if(TestvolumeSize == 256.0){devArea = Testfailure_Deviation_Volume_256; devDiameter = Testfailure_Deviation_Diameter_256;} + else if (TestvolumeSize == 512.0){devArea = Testfailure_Deviation_Volume_512; devDiameter = Testfailure_Deviation_Diameter_512;} + else{devArea = Testfailure_Deviation_Volume_128; devDiameter = Testfailure_Deviation_Diameter_128;} + + + std::string areatestName = TestName.append(" area"); + std::string diametertestName = TestName.append(" testing diameter"); + + //TODO think about the deviation, 1% makes no sense at all + MITK_TEST_CONDITION(std::abs(100 - testResults.percentageAreaCalcToPixel) < 1, areatestName ); + MITK_TEST_CONDITION(std::abs(100 - testResults.percentageRadiusToPixel) < 1, diametertestName ); + + #ifdef EXTRACTOR_DEBUG + MITK_INFO << TestName << " >>> " << "planeDistanceToSphereCenter: " << testResults.planeDistanceToSphereCenter; + MITK_INFO << "area in pixels: " << testResults.areaInPixel << " <-> area in mm: " << testResults.areaCalculated << " = " << testResults.percentageAreaCalcToPixel << "%"; + + MITK_INFO << "calculated diameter: " << testResults.diameterCalculated << " <-> diameter in mm: " << testResults.diameterInMM << " <-> diameter in pixel: " << testResults.diameterInPixel << " = " << testResults.percentageRadiusToPixel << "%"; + #endif + } + + + /* + * get the radius of the slice of a sphere based on pixel distance from edge to edge of the circle. + */ + template + static void TestSphereRadiusByItk (itk::Image* inputImage) + { + typedef itk::Image InputImageType; + + + //set the index to the middle of the image's edge at x and y axis + typename InputImageType::IndexType currentIndexX; + currentIndexX[0] = (int)(TestvolumeSize / 2.0); + currentIndexX[1] = 0; + + typename InputImageType::IndexType currentIndexY; + currentIndexY[0] = 0; + currentIndexY[1] = (int)(TestvolumeSize / 2.0); + + + //remember the last pixel value + double lastValueX = inputImage->GetPixel(currentIndexX); + double lastValueY = inputImage->GetPixel(currentIndexY); + + + //storage for the index marks + std::vector indicesX; + std::vector indicesY; + + + /*Get four indices on the edge of the circle*/ + while(currentIndexX[1] < TestvolumeSize && currentIndexX[0] < TestvolumeSize) + { + //move x direction + currentIndexX[1] += 1; + + //move y direction + currentIndexY[0] += 1; + + if(inputImage->GetPixel(currentIndexX) > lastValueX) + { + //mark the current index + typename InputImageType::IndexType markIndex; + markIndex[0] = currentIndexX[0]; + markIndex[1] = currentIndexX[1]; + + + indicesX.push_back(markIndex); + } + else if( inputImage->GetPixel(currentIndexX) < lastValueX) + { + //mark the current index + typename InputImageType::IndexType markIndex; + markIndex[0] = currentIndexX[0]; + markIndex[1] = currentIndexX[1] - 1;//value inside the sphere + + + indicesX.push_back(markIndex); + } + + if(inputImage->GetPixel(currentIndexY) > lastValueY) + { + //mark the current index + typename InputImageType::IndexType markIndex; + markIndex[0] = currentIndexY[0]; + markIndex[1] = currentIndexY[1]; + + + indicesY.push_back(markIndex); + } + else if( inputImage->GetPixel(currentIndexY) < lastValueY) + { + //mark the current index + typename InputImageType::IndexType markIndex; + markIndex[0] = currentIndexY[0]; + markIndex[1] = currentIndexY[1] - 1;//value inside the sphere + + + indicesY.push_back(markIndex); + } + + + //found both marks? + if(indicesX.size() == 2 && indicesY.size() == 2) break; + + //the new 'last' values + lastValueX = inputImage->GetPixel(currentIndexX); + lastValueY = inputImage->GetPixel(currentIndexY); + } + + /* + *If we are here we found the four marks on the edge of the circle. + *For the case our plane is rotated and shifted, we have to calculate the center of the circle, + *else the center is the intersection of both straight lines between the marks. + *When we have the center, the diameter of the circle will be checked to the reference value(math!). + */ + + //each distance from the first mark of each direction to the center of the straight line between the marks + double distanceToCenterX = std::abs(indicesX[0][1] - indicesX[1][1]) / 2.0; + double distanceToCenterY = std::abs(indicesY[0][0] - indicesY[1][0]) / 2.0; + + + //the center of the straight lines + typename InputImageType::IndexType centerX, centerY; + + centerX[0] = indicesX[0][0]; + centerX[1] = indicesX[0][1] + distanceToCenterX; + //TODO think about implicit cast to int. this is not the real center of the image, which could be between two pixels + + //centerY[0] = indicesY[0][0] + distanceToCenterY; + //centerY[1] = inidcesY[0][1]; + + typename InputImageType::IndexType currentIndex(centerX); + lastValueX = inputImage->GetPixel(currentIndex); + + long sumpixels = 0; + + std::vector diameterIndices; + + //move up + while(currentIndex[1] < TestvolumeSize) + { + currentIndex[1] += 1; + + if( inputImage->GetPixel(currentIndex) != lastValueX) + { + typename InputImageType::IndexType markIndex; + markIndex[0] = currentIndex[0]; + markIndex[1] = currentIndex[1] - 1; + + diameterIndices.push_back(markIndex); + break; + } + + sumpixels++; + lastValueX = inputImage->GetPixel(currentIndex); + } + + currentIndex[1] -= sumpixels; //move back to center to go in the other direction + lastValueX = inputImage->GetPixel(currentIndex); + + //move down + while(currentIndex[1] >= 0) + { + currentIndex[1] -= 1; + + if( inputImage->GetPixel(currentIndex) != lastValueX) + { + typename InputImageType::IndexType markIndex; + markIndex[0] = currentIndex[0]; + markIndex[1] = currentIndex[1];//outside sphere because we want to calculate the distance from edge to edge + + diameterIndices.push_back(markIndex); + break; + } + + sumpixels++; + lastValueX = inputImage->GetPixel(currentIndex); + } + + + + /* + *Now sumpixels should be the apromximate diameter of the circle. This is checked with the calculated diameter from the plane transformation(math). + */ + mitk::Point3D volumeCenter; + volumeCenter[0] = volumeCenter[1] = volumeCenter[2] = TestvolumeSize / 2.0; + + + double planeDistanceToSphereCenter = TestPlane->Distance(volumeCenter); + + + double sphereRadius = TestvolumeSize/4.0; + + //calculate the radius of the circle cut from the sphere by the plane + double diameter = 2.0 * std::sqrt(std::pow(sphereRadius, 2) - std::pow( planeDistanceToSphereCenter , 2)); + + double percentageRadiusToPixel = 100 / diameter * sumpixels; + + + + /* + *calculate the radius in mm by the both marks of the center line by using the world coordinates + */ + //get the points as 3D coordinates + mitk::Vector3D diameterPointRight, diameterPointLeft; + + diameterPointRight[2] = diameterPointLeft[2] = 0.0; + + diameterPointLeft[0] = diameterIndices[0][0]; + diameterPointLeft[1] = diameterIndices[0][1]; + + diameterPointRight[0] = diameterIndices[1][0]; + diameterPointRight[1] = diameterIndices[1][1]; + + //transform to worldcoordinates + TestVolume->GetGeometry()->IndexToWorld(diameterPointLeft, diameterPointLeft); + TestVolume->GetGeometry()->IndexToWorld(diameterPointRight, diameterPointRight); + + //euklidian distance + double diameterInMM = ( (diameterPointLeft * -1.0) + diameterPointRight).GetNorm(); + + + testResults.diameterInMM = diameterInMM; + testResults.diameterCalculated = diameter; + testResults.diameterInPixel = sumpixels; + testResults.percentageRadiusToPixel = percentageRadiusToPixel; + testResults.planeDistanceToSphereCenter = planeDistanceToSphereCenter; + } + + + + + + /*brute force the area pixel by pixel*/ + template + static void TestSphereAreaByItk (itk::Image* inputImage) + { + typedef itk::Image InputImageType; + typedef itk::ImageRegionConstIterator< InputImageType > ImageIterator; + + + ImageIterator imageIterator( inputImage, inputImage->GetLargestPossibleRegion() ); + imageIterator.GoToBegin(); + + int sumPixelsInArea = 0; + + while( !imageIterator.IsAtEnd() ) + { + if(inputImage->GetPixel(imageIterator.GetIndex()) == pixelValueSet) sumPixelsInArea++; + ++imageIterator; + } + + mitk::Point3D volumeCenter; + volumeCenter[0] = volumeCenter[1] = volumeCenter[2] = TestvolumeSize / 2.0; + + + double planeDistanceToSphereCenter = TestPlane->Distance(volumeCenter); + + double sphereRadius = TestvolumeSize/4.0; + + //calculate the radius of the circle cut from the sphere by the plane + double radius = std::sqrt(std::pow(sphereRadius, 2) - std::pow( planeDistanceToSphereCenter , 2)); + + double areaInMM = 3.14159265358979 * std::pow(radius, 2); + + + testResults.areaCalculated = areaInMM; + testResults.areaInPixel = sumPixelsInArea; + testResults.percentageAreaCalcToPixel = 100 / areaInMM * sumPixelsInArea; + } + + + + /* + * random a voxel. define plane through this voxel. reslice at the plane. compare the pixel vaues of the voxel + * in the volume with the pixel value in the resliced image. + * there are some indice shifting problems which causes the test to fail for oblique planes. seems like the chosen + * worldcoordinate is not corrresponding to the index in the 2D image. and so the pixel values are not the same as + * expected. + */ + static void PixelvalueBasedTest() + { + /* setup itk image */ + typedef itk::Image ImageType; + + typedef itk::ImageRegionConstIterator< ImageType > ImageIterator; + + ImageType::Pointer image = ImageType::New(); + + ImageType::IndexType start; + start[0] = start[1] = start[2] = 0; + + ImageType::SizeType size; + size[0] = size[1] = size[2] = 32; + + ImageType::RegionType imgRegion; + imgRegion.SetSize(size); + imgRegion.SetIndex(start); + + image->SetRegions(imgRegion); + image->SetSpacing(1.0); + image->Allocate(); + + + ImageIterator imageIterator( image, image->GetLargestPossibleRegion() ); + imageIterator.GoToBegin(); + + + unsigned short pixelValue = 0; + + //fill the image with distinct values + while ( !imageIterator.IsAtEnd() ) + { + image->SetPixel(imageIterator.GetIndex(), pixelValue); + ++imageIterator; + ++pixelValue; + } + /* end setup itk image */ + + + + mitk::Image::Pointer imageInMitk; + CastToMitkImage(image, imageInMitk); + + + + /*mitk::ImageWriter::Pointer writer = mitk::ImageWriter::New(); + writer->SetInput(imageInMitk); + std::string file = "C:\\Users\\schroedt\\Desktop\\cube.nrrd"; + writer->SetFileName(file); + writer->Update();*/ + + PixelvalueBasedTestByPlane(imageInMitk, mitk::PlaneGeometry::Frontal); + PixelvalueBasedTestByPlane(imageInMitk, mitk::PlaneGeometry::Sagittal); + PixelvalueBasedTestByPlane(imageInMitk, mitk::PlaneGeometry::Transversal); + + } + + static void PixelvalueBasedTestByPlane(mitk::Image* imageInMitk, mitk::PlaneGeometry::PlaneOrientation orientation){ + + typedef itk::Image ImageType; + + //set the seed of the rand function + srand((unsigned)time(0)); + + /* setup a random orthogonal plane */ + int sliceindex = 17;//rand() % 32; + bool isFrontside = true; + bool isRotated = false; + + + if( orientation == mitk::PlaneGeometry::Transversal) + { + /*isFrontside = false; + isRotated = true;*/ + } + + + + mitk::PlaneGeometry::Pointer plane = mitk::PlaneGeometry::New(); + + + plane->InitializeStandardPlane(imageInMitk->GetGeometry(), orientation, sliceindex, isFrontside, isRotated); + + mitk::Point3D origin = plane->GetOrigin(); + mitk::Vector3D normal; + normal = plane->GetNormal(); + + + normal.Normalize(); + + origin += normal * 0.5;//pixelspacing is 1, so half the spacing is 0.5 + + plane->SetOrigin(origin); + + //we dont need this any more, because we are only testing orthogonal planes + /*mitk::Vector3D rotationVector; + rotationVector[0] = randFloat(); + rotationVector[1] = randFloat(); + rotationVector[2] = randFloat(); + + + float degree = randFloat() * 180.0; + + mitk::RotationOperation* op = new mitk::RotationOperation(mitk::OpROTATE, plane->GetCenter(), rotationVector, degree); + plane->ExecuteOperation(op); + delete op;*/ + + /* end setup plane */ + + + /* define a point in the 3D volume. + * add the two axis vectors of the plane (each multiplied with a + * random number) to the origin. now the two random numbers + * become our index coordinates in the 2D image, because the + * length of the axis vectors is 1. + */ + mitk::Point3D planeOrigin = plane->GetOrigin(); + mitk::Vector3D axis0, axis1; + axis0 = plane->GetAxisVector(0); + axis1 = plane->GetAxisVector(1); + axis0.Normalize(); + axis1.Normalize(); + + + unsigned char n1 = 7;// rand() % 32; + unsigned char n2 = 13;// rand() % 32; + + + mitk::Point3D testPoint3DInWorld; + testPoint3DInWorld = planeOrigin + (axis0 * n1) + (axis1 * n2); + + //get the index of the point in the 3D volume + ImageType::IndexType testPoint3DInIndex; + imageInMitk->GetGeometry()->WorldToIndex(testPoint3DInWorld, testPoint3DInIndex); + + mitk::Index3D testPoint2DInIndex; + + /* end define a point in the 3D volume.*/ + + + //do reslicing at the plane + mitk::ExtractSliceFilter::Pointer slicer = mitk::ExtractSliceFilter::New(); + slicer->SetInput(imageInMitk); + slicer->SetWorldGeometry(plane); + + slicer->Update(); + + mitk::Image::Pointer slice = slicer->GetOutput(); + + // Get TestPoiont3D as Index in Slice + slice->GetGeometry()->WorldToIndex(testPoint3DInWorld,testPoint2DInIndex); + + + mitk::Point3D p, sliceIndexToWorld, imageIndexToWorld; + p[0] = testPoint2DInIndex[0]; + p[1] = testPoint2DInIndex[1]; + p[2] = testPoint2DInIndex[2]; + slice->GetGeometry()->IndexToWorld(p, sliceIndexToWorld); + + p[0] = testPoint3DInIndex[0]; + p[1] = testPoint3DInIndex[1]; + p[2] = testPoint3DInIndex[2]; + imageInMitk->GetGeometry()->IndexToWorld(p, imageIndexToWorld); + + + + //compare the pixelvalues of the defined point in the 3D volume with the value of the resliced image + unsigned short valueAt3DVolume = imageInMitk->GetPixelValueByIndex(testPoint3DInIndex);//image->GetPixel(testPoint3DInIndex); + unsigned short valueAt3DVolumeByWorld = imageInMitk->GetPixelValueByWorldCoordinate(testPoint3DInWorld); + unsigned short valueAtSlice = slice->GetPixelValueByIndex(testPoint2DInIndex); + + //valueAt3DVolume == valueAtSlice is not always working. because of rounding errors + //indices are shifted + MITK_TEST_CONDITION(valueAt3DVolume == valueAtSlice, "comparing pixelvalues for orthogonal plane"); + + + vtkSmartPointer imageInVtk = vtkSmartPointer::New(); + imageInVtk = imageInMitk->GetVtkImageData(); + vtkSmartPointer sliceInVtk = vtkSmartPointer::New(); + sliceInVtk = slice->GetVtkImageData(); + + double PixelvalueByMitkOutput = sliceInVtk->GetScalarComponentAsDouble(n1, n2, 0, 0); + double valueVTKinImage = imageInVtk->GetScalarComponentAsDouble(testPoint3DInIndex[0], testPoint3DInIndex[1], testPoint3DInIndex[2], 0); + + + /* Test that everything is working equally if vtkoutput is used instead of the default output + * from mitk ImageToImageFilter + */ + mitk::ExtractSliceFilter::Pointer slicerWithVtkOutput = mitk::ExtractSliceFilter::New(); + slicerWithVtkOutput->SetInput(imageInMitk); + slicerWithVtkOutput->SetWorldGeometry(plane); + slicerWithVtkOutput->SetVtkOutputRequest(true); + + slicerWithVtkOutput->Update(); + vtkSmartPointer vtkImageByVtkOutput = vtkSmartPointer::New(); + vtkImageByVtkOutput = slicerWithVtkOutput->GetVtkOutput(); + double PixelvalueByVtkOutput = vtkImageByVtkOutput->GetScalarComponentAsDouble(n1, n2, 0, 0); + + MITK_TEST_CONDITION(PixelvalueByMitkOutput == PixelvalueByVtkOutput, "testing convertion of image output vtk->mitk by reslicer"); + + + +/*================ mbilog outputs ===========================*/ +#ifdef EXTRACTOR_DEBUG + MITK_INFO << "\n" << "TESTINFO index: " << sliceindex << " orientation: " << orientation << " frontside: " << isFrontside << " rotated: " << isRotated; + MITK_INFO << "\n" << "slice index to world: " << sliceIndexToWorld; + MITK_INFO << "\n" << "image index to world: " << imageIndexToWorld; + + MITK_INFO << "\n" << "vtk: slice: " << PixelvalueByMitkOutput << ", image: "<< valueVTKinImage; + + MITK_INFO << "\n" << "testPoint3D InWorld" << testPoint3DInWorld << " is " << testPoint2DInIndex << " in 2D"; + MITK_INFO << "\n" << "randoms: " << ((int)n1) << ", " << ((int)n2); + MITK_INFO << "\n" << "point is inside plane: " << plane->IsInside(testPoint3DInWorld) << " and volume: " << imageInMitk->GetGeometry()->IsInside(testPoint3DInWorld); + + MITK_INFO << "\n" << "volume idx: " << testPoint3DInIndex << " = " << valueAt3DVolume ; + MITK_INFO << "\n" << "volume world: " << testPoint3DInWorld << " = " << valueAt3DVolumeByWorld ; + MITK_INFO << "\n" << "slice idx: " << testPoint2DInIndex << " = " << valueAtSlice ; + + mitk::Index3D curr; + curr[0] = curr[1] = curr[2] = 0; + + for( int i = 0; i < 32 ; ++i){ + for( int j = 0; j < 32; ++j){ + ++curr[1]; + if(slice->GetPixelValueByIndex(curr) == valueAt3DVolume){ + MITK_INFO << "\n" << valueAt3DVolume << " MATCHED mitk " << curr; + } + } + curr[1] = 0; + ++curr[0]; + } + + typedef itk::Image Image2DType; + + Image2DType::Pointer img = Image2DType::New(); + CastToItkImage(slice, img); + typedef itk::ImageRegionConstIterator< Image2DType > Iterator2D; + + Iterator2D iter(img, img->GetLargestPossibleRegion()); + iter.GoToBegin(); + while( !iter.IsAtEnd() ){ + + if(img->GetPixel(iter.GetIndex()) == valueAt3DVolume) MITK_INFO << "\n" << valueAt3DVolume << " MATCHED itk " << iter.GetIndex(); + + ++iter; + } +#endif //EXTRACTOR_DEBUG + } + + + /* random a float value */ + static float randFloat(){ return (((float)rand()+1.0) / ((float)RAND_MAX + 1.0)) + (((float)rand()+1.0) / ((float)RAND_MAX + 1.0)) / ((float)RAND_MAX + 1.0);} + + + + + /* create a sphere with the size of the given testVolumeSize*/ + static void InitializeTestVolume() + { + + #ifdef CREATE_VOLUME + + + //do sphere creation + ItkVolumeGeneration(); + + #ifdef SAVE_VOLUME + //save in file + mitk::ImageWriter::Pointer writer = mitk::ImageWriter::New(); + writer->SetInput(TestVolume); + + + std::string file; + + std::ostringstream filename; + filename << "C:\\home\\schroedt\\MITK\\Modules\\ImageExtraction\\Testing\\Data\\sphere_"; + filename << TestvolumeSize; + filename << ".nrrd"; + + file = filename.str(); + + writer->SetFileName(file); + writer->Update(); + #endif//SAVE_VOLUME + + #endif + + + #ifndef CREATE_VOLUME //read from file + + mitk::StandardFileLocations::Pointer locator = mitk::StandardFileLocations::GetInstance(); + + std::string filename = locator->FindFile("sphere_512.nrrd.mhd", "Modules/ImageExtraction/Testing/Data"); + + mitk::ItkImageFileReader::Pointer reader = mitk::ItkImageFileReader::New(); + reader->SetFileName(filename); + + reader->Update(); + TestVolume = reader->GetOutput(); + + #endif + + #ifdef CALC_TESTFAILURE_DEVIATION + //get the TestFailureDeviation in % + AccessFixedDimensionByItk(TestVolume, CalcTestFailureDeviation, 3); + #endif + } + + + //the test result of the sphere reslice + struct SliceProperties{ + double planeDistanceToSphereCenter; + double diameterInMM; + double diameterInPixel; + double diameterCalculated; + double percentageRadiusToPixel; + double areaCalculated; + double areaInPixel; + double percentageAreaCalcToPixel; + }; + + + static mitk::Image::Pointer TestVolume; + static double TestvolumeSize; + static mitk::PlaneGeometry::Pointer TestPlane; + static std::string TestName; + static unsigned char pixelValueSet; + static SliceProperties testResults; + static double TestFailureDeviation; + + + +private: + + /* + * Generate a sphere with a radius of TestvolumeSize / 4.0 + */ + static void ItkVolumeGeneration () + { + typedef itk::Image TestVolumeType; + + typedef itk::ImageRegionConstIterator< TestVolumeType > ImageIterator; + + TestVolumeType::Pointer sphereImage = TestVolumeType::New(); + + TestVolumeType::IndexType start; + start[0] = start[1] = start[2] = 0; + + TestVolumeType::SizeType size; + size[0] = size[1] = size[2] = TestvolumeSize; + + TestVolumeType::RegionType imgRegion; + imgRegion.SetSize(size); + imgRegion.SetIndex(start); + + sphereImage->SetRegions(imgRegion); + sphereImage->SetSpacing(1.0); + sphereImage->Allocate(); + + sphereImage->FillBuffer(0); + + + mitk::Vector3D center; + center[0] = center[1] = center[2] = TestvolumeSize / 2.0; + + + double radius = TestvolumeSize / 4.0; + + double pixelValue = pixelValueSet; + + double distanceToCenter = 0.0; + + + ImageIterator imageIterator( sphereImage, sphereImage->GetLargestPossibleRegion() ); + imageIterator.GoToBegin(); + + + mitk::Vector3D currentVoxelInIndex; + + while ( !imageIterator.IsAtEnd() ) + { + currentVoxelInIndex[0] = imageIterator.GetIndex()[0]; + currentVoxelInIndex[1] = imageIterator.GetIndex()[1]; + currentVoxelInIndex[2] = imageIterator.GetIndex()[2]; + + distanceToCenter = (center + ( currentVoxelInIndex * -1.0 )).GetNorm(); + + //if distance to center is smaller then the radius of the sphere + if( distanceToCenter < radius) + { + sphereImage->SetPixel(imageIterator.GetIndex(), pixelValue); + } + + ++imageIterator; + } + + CastToMitkImage(sphereImage, TestVolume); + } + + + + + + /* calculate the devation of the voxel object to the mathematical sphere object. + * this is use to make a statement about the accuracy of the resliced image, eg. the circle's diameter or area. + */ + template + static void CalcTestFailureDeviation (itk::Image* inputImage) + { + typedef itk::Image InputImageType; + typedef itk::ImageRegionConstIterator< InputImageType > ImageIterator; + + ImageIterator iterator(inputImage, inputImage->GetLargestPossibleRegion()); + iterator.GoToBegin(); + + int volumeInPixel = 0; + + while( !iterator.IsAtEnd() ) + { + if(inputImage->GetPixel(iterator.GetIndex()) == pixelValueSet) volumeInPixel++; + ++iterator; + } + + double diameter = TestvolumeSize / 2.0; + double volumeCalculated = (1.0 / 6.0) * 3.14159265358979 * std::pow(diameter, 3); + + + double volumeDeviation = std::abs( 100 - (100 / volumeCalculated * volumeInPixel) ); + + + typename InputImageType::IndexType index; + index[0] = index[1] = TestvolumeSize / 2.0; + index[2] = 0; + + int sumpixels = 0; + while (index[2] < TestvolumeSize ) + { + if(inputImage->GetPixel(index) == pixelValueSet) sumpixels++; + index[2] += 1; + } + + + double diameterDeviation = std::abs( 100 - (100 / diameter * sumpixels) ); + #ifdef DEBUG + MITK_INFO << "volume deviation: " << volumeDeviation << " diameter deviation:" << diameterDeviation; + #endif + mitkExtractSliceFilterTestClass::TestFailureDeviation = (volumeDeviation + diameterDeviation) / 2.0; + } + + +}; +/*================ #END class ================*/ + + + + + +/*================#BEGIN Instanciation of members ================*/ +mitk::Image::Pointer mitkExtractSliceFilterTestClass::TestVolume = mitk::Image::New(); +double mitkExtractSliceFilterTestClass::TestvolumeSize = 256.0; +mitk::PlaneGeometry::Pointer mitkExtractSliceFilterTestClass::TestPlane = mitk::PlaneGeometry::New(); +std::string mitkExtractSliceFilterTestClass::TestName = ""; +unsigned char mitkExtractSliceFilterTestClass::pixelValueSet = 255; +mitkExtractSliceFilterTestClass::SliceProperties mitkExtractSliceFilterTestClass::testResults = {-1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0}; +double mitkExtractSliceFilterTestClass::TestFailureDeviation = 0.0; +/*================ #END Instanciation of members ================*/ + + + + + +/*================ #BEGIN test main ================*/ +int mitkExtractSliceFilterTest(int argc, char* argv[]) +{ + + MITK_TEST_BEGIN("mitkExtractSliceFilterTest") + + + //pixelvalue based testing + mitkExtractSliceFilterTestClass::PixelvalueBasedTest(); + + //initialize sphere test volume + mitkExtractSliceFilterTestClass::InitializeTestVolume(); + + + mitk::Vector3D spacing = mitkExtractSliceFilterTestClass::TestVolume->GetGeometry()->GetSpacing(); + + + //the center of the sphere = center of image + double sphereCenter = mitkExtractSliceFilterTestClass::TestvolumeSize / 2.0; + + double planeSize = mitkExtractSliceFilterTestClass::TestvolumeSize; + + + + /* transversal plane */ + mitk::PlaneGeometry::Pointer geometryTransversal = mitk::PlaneGeometry::New(); + geometryTransversal->InitializeStandardPlane(planeSize, planeSize, spacing, mitk::PlaneGeometry::Transversal, sphereCenter, false, true); + geometryTransversal->ChangeImageGeometryConsideringOriginOffset(true); + + mitk::Point3D origin = geometryTransversal->GetOrigin(); + mitk::Vector3D normal; + normal = geometryTransversal->GetNormal(); + + + normal.Normalize(); + + origin += normal * 0.5;//pixelspacing is 1, so half the spacing is 0.5 + + //geometryTransversal->SetOrigin(origin); + + mitkExtractSliceFilterTestClass::TestSlice(geometryTransversal, "Testing transversal plane"); + /* end transversal plane */ + + + + /* sagittal plane */ + mitk::PlaneGeometry::Pointer geometrySagital = mitk::PlaneGeometry::New(); + geometrySagital->InitializeStandardPlane(planeSize, planeSize, spacing, mitk::PlaneGeometry::Sagittal, sphereCenter, true, false); + geometrySagital->ChangeImageGeometryConsideringOriginOffset(true); + + origin = geometrySagital->GetOrigin(); + normal = geometrySagital->GetNormal(); + + + normal.Normalize(); + + origin += normal * 0.5;//pixelspacing is 1, so half the spacing is 0.5 + + //geometrySagital->SetOrigin(origin); + + mitkExtractSliceFilterTestClass::TestSlice(geometrySagital, "Testing sagittal plane"); + /* sagittal plane */ + + + + /* sagittal shifted plane */ + mitk::PlaneGeometry::Pointer geometrySagitalShifted = mitk::PlaneGeometry::New(); + geometrySagitalShifted->InitializeStandardPlane(planeSize, planeSize, spacing, mitk::PlaneGeometry::Sagittal, (sphereCenter - 14), true, false); + geometrySagitalShifted->ChangeImageGeometryConsideringOriginOffset(true); + + origin = geometrySagitalShifted->GetOrigin(); + normal = geometrySagitalShifted->GetNormal(); + + + normal.Normalize(); + + origin += normal * 0.5;//pixelspacing is 1, so half the spacing is 0.5 + + //geometrySagitalShifted->SetOrigin(origin); + + mitkExtractSliceFilterTestClass::TestSlice(geometrySagitalShifted, "Testing sagittal plane shifted"); + /* end sagittal shifted plane */ + + + + /* coronal plane */ + mitk::PlaneGeometry::Pointer geometryCoronal = mitk::PlaneGeometry::New(); + geometryCoronal->InitializeStandardPlane(planeSize, planeSize, spacing, mitk::PlaneGeometry::Frontal, sphereCenter, true, false); + geometryCoronal->ChangeImageGeometryConsideringOriginOffset(true); + + origin = geometryCoronal->GetOrigin(); + normal = geometryCoronal->GetNormal(); + + + normal.Normalize(); + + origin += normal * 0.5;//pixelspacing is 1, so half the spacing is 0.5 + + //geometryCoronal->SetOrigin(origin); + + mitkExtractSliceFilterTestClass::TestSlice(geometryCoronal, "Testing coronal plane"); + /* end coronal plane */ + + + + /* oblique plane */ + mitk::PlaneGeometry::Pointer obliquePlane = mitk::PlaneGeometry::New(); + obliquePlane->InitializeStandardPlane(planeSize, planeSize, spacing, mitk::PlaneGeometry::Sagittal, sphereCenter, true, false); + obliquePlane->ChangeImageGeometryConsideringOriginOffset(true); + + origin = obliquePlane->GetOrigin(); + normal = obliquePlane->GetNormal(); + + + normal.Normalize(); + + origin += normal * 0.5;//pixelspacing is 1, so half the spacing is 0.5 + + //obliquePlane->SetOrigin(origin); + + mitk::Vector3D rotationVector; + rotationVector[0] = 0.2; + rotationVector[1] = 0.4; + rotationVector[2] = 0.62; + + float degree = 37.0; + + mitk::RotationOperation* op = new mitk::RotationOperation(mitk::OpROTATE, obliquePlane->GetCenter(), rotationVector, degree); + obliquePlane->ExecuteOperation(op); + delete op; + + mitkExtractSliceFilterTestClass::TestSlice(obliquePlane, "Testing oblique plane"); + /* end oblique plane */ + + + + #ifdef SHOW_SLICE_IN_RENDER_WINDOW + /*================ #BEGIN vtk render code ================*/ + + //set reslicer for renderwindow + mitk::ItkImageFileReader::Pointer reader = mitk::ItkImageFileReader::New(); + + std::string filename = "C:\\home\\Pics\\Pic3D.nrrd"; + reader->SetFileName(filename); + + reader->Update(); + + mitk::Image::Pointer pic = reader->GetOutput(); + vtkSmartPointer slicer = vtkSmartPointer::New(); + + slicer->SetInput(pic->GetVtkImageData()); + + + mitk::PlaneGeometry::Pointer obliquePl = mitk::PlaneGeometry::New(); + obliquePl->InitializeStandardPlane(pic->GetGeometry(), mitk::PlaneGeometry::Sagittal, pic->GetGeometry()->GetCenter()[0], true, false); + obliquePl->ChangeImageGeometryConsideringOriginOffset(true); + + mitk::Point3D origin2 = obliquePl->GetOrigin(); + mitk::Vector3D n; + n = obliquePl->GetNormal(); + + + n.Normalize(); + + origin2 += n * 0.5;//pixelspacing is 1, so half the spacing is 0.5 + + obliquePl->SetOrigin(origin2); + + mitk::Vector3D rotation; + rotation[0] = 0.534307; + rotation[1] = 0.000439605; + rotation[2] = 0.423017; + MITK_INFO << rotation; + + float rotateDegree = 70; + + mitk::RotationOperation* operation = new mitk::RotationOperation(mitk::OpROTATE, obliquePl->GetCenter(), rotationVector, degree); + obliquePl->ExecuteOperation(operation); + delete operation; + + + double origin[3]; + origin[0] = obliquePl->GetOrigin()[0]; + origin[1] = obliquePl->GetOrigin()[1]; + origin[2] = obliquePl->GetOrigin()[2]; + slicer->SetResliceAxesOrigin(origin); + + + mitk::Vector3D right, bottom, normal; + right = obliquePl->GetAxisVector( 0 ); + bottom = obliquePl->GetAxisVector( 1 ); + normal = obliquePl->GetNormal(); + + right.Normalize(); + bottom.Normalize(); + normal.Normalize(); + + double cosines[9]; + + mitk::vnl2vtk(right.GetVnlVector(), cosines);//x + + mitk::vnl2vtk(bottom.GetVnlVector(), cosines + 3);//y + + mitk::vnl2vtk(normal.GetVnlVector(), cosines + 6);//n + + slicer->SetResliceAxesDirectionCosines(cosines); + + slicer->SetOutputDimensionality(2); + slicer->Update(); + + + //set vtk renderwindow + vtkSmartPointer vtkPlane = vtkSmartPointer::New(); + vtkPlane->SetOrigin(0.0, 0.0, 0.0); + + //These two points define the axes of the plane in combination with the origin. + //Point 1 is the x-axis and point 2 the y-axis. + //Each plane is transformed according to the view (transversal, coronal and saggital) afterwards. + vtkPlane->SetPoint1(1.0, 0.0, 0.0); //P1: (xMax, yMin, depth) + vtkPlane->SetPoint2(0.0, 1.0, 0.0); //P2: (xMin, yMax, depth) + //these are not the correct values for all slices, only a square plane by now + + vtkSmartPointer imageMapper = vtkSmartPointer::New(); + imageMapper->SetInputConnection(vtkPlane->GetOutputPort()); + + + vtkSmartPointer lookupTable = vtkSmartPointer::New(); + + //built a default lookuptable + lookupTable->SetRampToLinear(); + lookupTable->SetSaturationRange( 0.0, 0.0 ); + lookupTable->SetHueRange( 0.0, 0.0 ); + lookupTable->SetValueRange( 0.0, 1.0 ); + lookupTable->Build(); + //map all black values to transparent + lookupTable->SetTableValue(0, 0.0, 0.0, 0.0, 0.0); + lookupTable->SetRange(-255.0, 255.0); + //lookupTable->SetRange(-1022.0, 1184.0);//pic3D range + + vtkSmartPointer texture = vtkSmartPointer::New(); + + + texture->SetInput(slicer->GetOutput()); + + texture->SetLookupTable(lookupTable); + + texture->SetMapColorScalarsThroughLookupTable(true); + + vtkSmartPointer imageActor = vtkSmartPointer::New(); + imageActor->SetMapper(imageMapper); + imageActor->SetTexture(texture); + + + // Setup renderers + vtkSmartPointer renderer = vtkSmartPointer::New(); + renderer->AddActor(imageActor); + + // Setup render window + vtkSmartPointer renderWindow = vtkSmartPointer::New(); + renderWindow->AddRenderer(renderer); + + // Setup render window interactor + vtkSmartPointer renderWindowInteractor = vtkSmartPointer::New(); + vtkSmartPointer style = vtkSmartPointer::New(); + renderWindowInteractor->SetInteractorStyle(style); + + // Render and start interaction + renderWindowInteractor->SetRenderWindow(renderWindow); + //renderer->AddViewProp(imageActor); + + + renderWindow->Render(); + renderWindowInteractor->Start(); + // always end with this! + /*================ #END vtk render code ================*/ + #endif //SHOW_SLICE_IN_RENDER_WINDOW + + + MITK_TEST_END() +} diff --git a/Core/Code/files.cmake b/Core/Code/files.cmake index 7982bf5875..d2d4b34679 100644 --- a/Core/Code/files.cmake +++ b/Core/Code/files.cmake @@ -1,300 +1,302 @@ set(H_FILES Algorithms/itkImportMitkImageContainer.h Algorithms/itkImportMitkImageContainer.txx Algorithms/itkLocalVariationImageFilter.h Algorithms/itkLocalVariationImageFilter.txx Algorithms/itkMITKScalarImageToHistogramGenerator.h Algorithms/itkMITKScalarImageToHistogramGenerator.txx Algorithms/itkTotalVariationDenoisingImageFilter.h Algorithms/itkTotalVariationDenoisingImageFilter.txx Algorithms/itkTotalVariationSingleIterationImageFilter.h Algorithms/itkTotalVariationSingleIterationImageFilter.txx Algorithms/mitkBilateralFilter.h Algorithms/mitkBilateralFilter.cpp Algorithms/mitkImageAccessByItk.h Algorithms/mitkImageCast.h Algorithms/mitkImageToItk.h Algorithms/mitkImageToItk.txx Algorithms/mitkInstantiateAccessFunctions.h Algorithms/mitkITKImageImport.h Algorithms/mitkITKImageImport.txx Algorithms/mitkPixelTypeList.h # Preprocessor macros taken from Boost Algorithms/mitkPPArithmeticDec.h Algorithms/mitkPPArgCount.h Algorithms/mitkPPCat.h Algorithms/mitkPPConfig.h Algorithms/mitkPPControlExprIIf.h Algorithms/mitkPPControlIf.h Algorithms/mitkPPControlIIf.h Algorithms/mitkPPDebugError.h Algorithms/mitkPPDetailAutoRec.h Algorithms/mitkPPDetailDMCAutoRec.h Algorithms/mitkPPExpand.h Algorithms/mitkPPFacilitiesEmpty.h Algorithms/mitkPPFacilitiesExpand.h Algorithms/mitkPPLogicalBool.h Algorithms/mitkPPRepetitionDetailDMCFor.h Algorithms/mitkPPRepetitionDetailEDGFor.h Algorithms/mitkPPRepetitionDetailFor.h Algorithms/mitkPPRepetitionDetailMSVCFor.h Algorithms/mitkPPRepetitionFor.h Algorithms/mitkPPSeqElem.h Algorithms/mitkPPSeqForEach.h Algorithms/mitkPPSeqForEachProduct.h Algorithms/mitkPPSeq.h Algorithms/mitkPPSeqEnum.h Algorithms/mitkPPSeqSize.h Algorithms/mitkPPSeqToTuple.h Algorithms/mitkPPStringize.h Algorithms/mitkPPTupleEat.h Algorithms/mitkPPTupleElem.h Algorithms/mitkPPTupleRem.h Algorithms/mitkClippedSurfaceBoundsCalculator.h - + Algorithms/mitkExtractSliceFilter.h + Interfaces/mitkIDataNodeReader.h IO/mitkPixelTypeTraits.h DataManagement/mitkCommon.h Interactions/mitkEventMapperAddOn.h Common/mitkExceptionMacro.h ) set(CPP_FILES Algorithms/mitkBaseDataSource.cpp Algorithms/mitkBaseProcess.cpp Algorithms/mitkCoreObjectFactoryBase.cpp Algorithms/mitkCoreObjectFactory.cpp Algorithms/mitkDataNodeFactory.cpp Algorithms/mitkDataNodeSource.cpp Algorithms/mitkGeometry2DDataToSurfaceFilter.cpp Algorithms/mitkHistogramGenerator.cpp Algorithms/mitkImageCaster.cpp Algorithms/mitkImageCastPart1.cpp Algorithms/mitkImageCastPart2.cpp Algorithms/mitkImageCastPart3.cpp Algorithms/mitkImageCastPart4.cpp Algorithms/mitkImageChannelSelector.cpp Algorithms/mitkImageSliceSelector.cpp Algorithms/mitkImageSource.cpp Algorithms/mitkImageTimeSelector.cpp Algorithms/mitkImageToImageFilter.cpp Algorithms/mitkPointSetSource.cpp Algorithms/mitkPointSetToPointSetFilter.cpp Algorithms/mitkRGBToRGBACastImageFilter.cpp Algorithms/mitkSubImageSelector.cpp Algorithms/mitkSurfaceSource.cpp Algorithms/mitkSurfaceToSurfaceFilter.cpp Algorithms/mitkUIDGenerator.cpp Algorithms/mitkVolumeCalculator.cpp Algorithms/mitkClippedSurfaceBoundsCalculator.cpp + Algorithms/mitkExtractSliceFilter.cpp Controllers/mitkBaseController.cpp Controllers/mitkCallbackFromGUIThread.cpp Controllers/mitkCameraController.cpp Controllers/mitkCameraRotationController.cpp Controllers/mitkCoreActivator.cpp Controllers/mitkFocusManager.cpp Controllers/mitkLimitedLinearUndo.cpp Controllers/mitkOperationEvent.cpp Controllers/mitkPlanePositionManager.cpp Controllers/mitkProgressBar.cpp Controllers/mitkRenderingManager.cpp Controllers/mitkSliceNavigationController.cpp Controllers/mitkSlicesCoordinator.cpp Controllers/mitkSlicesRotator.cpp Controllers/mitkSlicesSwiveller.cpp Controllers/mitkStatusBar.cpp Controllers/mitkStepper.cpp Controllers/mitkTestManager.cpp Controllers/mitkUndoController.cpp Controllers/mitkVerboseLimitedLinearUndo.cpp Controllers/mitkVtkInteractorCameraController.cpp Controllers/mitkVtkLayerController.cpp DataManagement/mitkAbstractTransformGeometry.cpp DataManagement/mitkAnnotationProperty.cpp DataManagement/mitkApplicationCursor.cpp DataManagement/mitkBaseData.cpp DataManagement/mitkBaseProperty.cpp DataManagement/mitkClippingProperty.cpp DataManagement/mitkChannelDescriptor.cpp DataManagement/mitkColorProperty.cpp DataManagement/mitkDataStorage.cpp #DataManagement/mitkDataTree.cpp DataManagement/mitkDataNode.cpp #DataManagement/mitkDataTreeStorage.cpp DataManagement/mitkDisplayGeometry.cpp DataManagement/mitkEnumerationProperty.cpp DataManagement/mitkGeometry2D.cpp DataManagement/mitkGeometry2DData.cpp DataManagement/mitkGeometry3D.cpp DataManagement/mitkGeometryData.cpp DataManagement/mitkGroupTagProperty.cpp DataManagement/mitkImage.cpp DataManagement/mitkImageDataItem.cpp DataManagement/mitkImageDescriptor.cpp DataManagement/mitkImageStatisticsHolder.cpp DataManagement/mitkLandmarkBasedCurvedGeometry.cpp DataManagement/mitkLandmarkProjectorBasedCurvedGeometry.cpp DataManagement/mitkLandmarkProjector.cpp DataManagement/mitkLevelWindow.cpp DataManagement/mitkLevelWindowManager.cpp DataManagement/mitkLevelWindowPreset.cpp DataManagement/mitkLevelWindowProperty.cpp DataManagement/mitkLookupTable.cpp DataManagement/mitkLookupTables.cpp # specializations of GenericLookupTable DataManagement/mitkMemoryUtilities.cpp DataManagement/mitkModalityProperty.cpp DataManagement/mitkModeOperation.cpp DataManagement/mitkNodePredicateAnd.cpp DataManagement/mitkNodePredicateBase.cpp DataManagement/mitkNodePredicateCompositeBase.cpp DataManagement/mitkNodePredicateData.cpp DataManagement/mitkNodePredicateDataType.cpp DataManagement/mitkNodePredicateDimension.cpp DataManagement/mitkNodePredicateFirstLevel.cpp DataManagement/mitkNodePredicateNot.cpp DataManagement/mitkNodePredicateOr.cpp DataManagement/mitkNodePredicateProperty.cpp DataManagement/mitkNodePredicateSource.cpp DataManagement/mitkPlaneOrientationProperty.cpp DataManagement/mitkPlaneGeometry.cpp DataManagement/mitkPlaneOperation.cpp DataManagement/mitkPointOperation.cpp DataManagement/mitkPointSet.cpp DataManagement/mitkProperties.cpp DataManagement/mitkPropertyList.cpp DataManagement/mitkRestorePlanePositionOperation.cpp DataManagement/mitkRotationOperation.cpp DataManagement/mitkSlicedData.cpp DataManagement/mitkSlicedGeometry3D.cpp DataManagement/mitkSmartPointerProperty.cpp DataManagement/mitkStandaloneDataStorage.cpp DataManagement/mitkStateTransitionOperation.cpp DataManagement/mitkStringProperty.cpp DataManagement/mitkSurface.cpp DataManagement/mitkSurfaceOperation.cpp DataManagement/mitkThinPlateSplineCurvedGeometry.cpp DataManagement/mitkTimeSlicedGeometry.cpp DataManagement/mitkTransferFunction.cpp DataManagement/mitkTransferFunctionProperty.cpp DataManagement/mitkTransferFunctionInitializer.cpp DataManagement/mitkVector.cpp DataManagement/mitkVtkInterpolationProperty.cpp DataManagement/mitkVtkRepresentationProperty.cpp DataManagement/mitkVtkResliceInterpolationProperty.cpp DataManagement/mitkVtkScalarModeProperty.cpp DataManagement/mitkVtkVolumeRenderingProperty.cpp DataManagement/mitkWeakPointerProperty.cpp DataManagement/mitkShaderProperty.cpp DataManagement/mitkResliceMethodProperty.cpp DataManagement/mitkMaterial.cpp Interactions/mitkAction.cpp Interactions/mitkAffineInteractor.cpp Interactions/mitkCoordinateSupplier.cpp Interactions/mitkDisplayCoordinateOperation.cpp Interactions/mitkDisplayInteractor.cpp Interactions/mitkDisplayPositionEvent.cpp Interactions/mitkDisplayVectorInteractor.cpp Interactions/mitkDisplayVectorInteractorLevelWindow.cpp Interactions/mitkDisplayVectorInteractorScroll.cpp Interactions/mitkEvent.cpp Interactions/mitkEventDescription.cpp Interactions/mitkEventMapper.cpp Interactions/mitkGlobalInteraction.cpp Interactions/mitkInteractor.cpp Interactions/mitkMouseModeSwitcher.cpp Interactions/mitkMouseMovePointSetInteractor.cpp Interactions/mitkMoveSurfaceInteractor.cpp Interactions/mitkNodeDepententPointSetInteractor.cpp Interactions/mitkPointSetInteractor.cpp Interactions/mitkPositionEvent.cpp Interactions/mitkPositionTracker.cpp Interactions/mitkState.cpp Interactions/mitkStateEvent.cpp Interactions/mitkStateMachine.cpp Interactions/mitkStateMachineFactory.cpp Interactions/mitkTransition.cpp Interactions/mitkWheelEvent.cpp Interactions/mitkKeyEvent.cpp Interactions/mitkVtkEventAdapter.cpp Interactions/mitkVtkInteractorStyle.cxx Interactions/mitkCrosshairPositionEvent.cpp IO/mitkBaseDataIOFactory.cpp IO/mitkCoreDataNodeReader.cpp IO/mitkDicomSeriesReader.cpp IO/mitkFileReader.cpp IO/mitkFileSeriesReader.cpp IO/mitkFileWriter.cpp #IO/mitkIpPicGet.c IO/mitkImageGenerator.cpp IO/mitkImageWriter.cpp IO/mitkImageWriterFactory.cpp IO/mitkItkImageFileIOFactory.cpp IO/mitkItkImageFileReader.cpp IO/mitkItkPictureWrite.cpp IO/mitkIOUtil.cpp IO/mitkLookupTableProperty.cpp IO/mitkOperation.cpp #IO/mitkPicFileIOFactory.cpp #IO/mitkPicFileReader.cpp #IO/mitkPicFileWriter.cpp #IO/mitkPicHelper.cpp #IO/mitkPicVolumeTimeSeriesIOFactory.cpp #IO/mitkPicVolumeTimeSeriesReader.cpp IO/mitkPixelType.cpp IO/mitkPointSetIOFactory.cpp IO/mitkPointSetReader.cpp IO/mitkPointSetWriter.cpp IO/mitkPointSetWriterFactory.cpp IO/mitkRawImageFileReader.cpp IO/mitkStandardFileLocations.cpp IO/mitkSTLFileIOFactory.cpp IO/mitkSTLFileReader.cpp IO/mitkSurfaceVtkWriter.cpp IO/mitkSurfaceVtkWriterFactory.cpp IO/mitkVtiFileIOFactory.cpp IO/mitkVtiFileReader.cpp IO/mitkVtkImageIOFactory.cpp IO/mitkVtkImageReader.cpp IO/mitkVtkSurfaceIOFactory.cpp IO/mitkVtkSurfaceReader.cpp IO/vtkPointSetXMLParser.cpp IO/mitkLog.cpp Rendering/mitkBaseRenderer.cpp Rendering/mitkVtkMapper2D.cpp Rendering/mitkVtkMapper3D.cpp Rendering/mitkRenderWindowFrame.cpp Rendering/mitkGeometry2DDataMapper2D.cpp Rendering/mitkGeometry2DDataVtkMapper3D.cpp Rendering/mitkGLMapper2D.cpp Rendering/mitkGradientBackground.cpp Rendering/mitkManufacturerLogo.cpp Rendering/mitkMapper2D.cpp Rendering/mitkMapper3D.cpp Rendering/mitkMapper.cpp Rendering/mitkPointSetGLMapper2D.cpp Rendering/mitkPointSetVtkMapper3D.cpp Rendering/mitkPolyDataGLMapper2D.cpp Rendering/mitkSurfaceGLMapper2D.cpp Rendering/mitkSurfaceVtkMapper3D.cpp Rendering/mitkVolumeDataVtkMapper3D.cpp Rendering/mitkVtkPropRenderer.cpp Rendering/mitkVtkWidgetRendering.cpp Rendering/vtkMitkRectangleProp.cpp Rendering/vtkMitkRenderProp.cpp Rendering/mitkVtkEventProvider.cpp Rendering/mitkRenderWindow.cpp Rendering/mitkRenderWindowBase.cpp Rendering/mitkShaderRepository.cpp Rendering/mitkImageVtkMapper2D.cpp Rendering/vtkMitkThickSlicesFilter.cpp Rendering/vtkMitkApplyLevelWindowToRGBFilter.cpp Common/mitkException.cpp ) list(APPEND CPP_FILES ${CppMicroServices_SOURCES}) diff --git a/Modules/ImageExtraction/Testing/files.cmake b/Modules/ImageExtraction/Testing/files.cmake index 6282bcbc93..c983818824 100644 --- a/Modules/ImageExtraction/Testing/files.cmake +++ b/Modules/ImageExtraction/Testing/files.cmake @@ -1,12 +1,13 @@ set(MODULE_IMAGE_TESTS - mitkExtractImageFilterTest.cpp + + ) set(MODULE_TESTIMAGES US4DCyl.pic.gz Pic3D.pic.gz Pic2DplusT.pic.gz BallBinary30x30x30.pic.gz Png2D-bw.png binary.stl ball.stl ) diff --git a/Modules/ImageExtraction/Testing/mitkExtractDirectedPlaneImageFilterTest.cpp b/Modules/ImageExtraction/Testing/mitkExtractDirectedPlaneImageFilterTest.cpp new file mode 100644 index 0000000000..93846b8217 --- /dev/null +++ b/Modules/ImageExtraction/Testing/mitkExtractDirectedPlaneImageFilterTest.cpp @@ -0,0 +1,297 @@ +/*========================================================================= + +Program: Medical Imaging & Interaction Toolkit +Language: C++ +Date: $Date$ +Version: $Revision: 7837 $ + +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 "mitkExtractDirectedPlaneImageFilter.h" +//#include "mitkStandardFileLocations.h" +// +//#include +//#include +//#include +//#include +//#include +// +//#include "mitkTestingMacros.h" +// +//#include +// +// +//class ExtractionTesting{ +// +//public: +// +// struct Testcase +// { +// int number; +// std::string name; +// std::string imageFilename; +// std::string referenceImageFilename; +// bool success; +// mitk::Geometry2D::Pointer (*GetPlane) (void); +// }; +// +// static void DoTesting(Testcase &testcase) +// { +// mitk::Image::Pointer image = GetImageToTest(testcase.imageFilename); +// if ( image.IsNull){ +// testcase.success = false; +// return; +// } +// +// /*mitk::Image::Pointer referenceImage = GetImageToTest(testcase.referenceImageFilename); +// if ( referenceImage.IsNull){ +// testcase.success = false; +// return; +// } +// +// mitk::Geometry2D::Pointer directedGeometry2D = testcase.GetPlane(); +// if(directedGeometry2D.IsNull){ +// testcase.success = false;*/ +// } +// +// //put testing here +// //TODO vtkIMageREslice setup +// //vtkSmartPointer colorImage = image->GetVtkImageData(); +// +// vtkSmartPointer imageMapper = vtkSmartPointer::New(); +// imageMapper->SetInput(colorImage); +// +// +// vtkSmartPointer imageActor = vtkSmartPointer::New(); +// imageActor->SetMapper(imageMapper); +// //imageActor->SetPosition(20, 20); +// +// // Setup renderers +// vtkSmartPointer renderer = vtkSmartPointer::New(); +// +// // Setup render window +// vtkSmartPointer renderWindow = vtkSmartPointer::New(); +// renderWindow->AddRenderer(renderer); +// +// // Setup render window interactor +// vtkSmartPointer renderWindowInteractor = vtkSmartPointer::New(); +// vtkSmartPointer style = vtkSmartPointer::New(); +// renderWindowInteractor->SetInteractorStyle(style); +// +// // Render and start interaction +// renderWindowInteractor->SetRenderWindow(renderWindow); +// //renderer->AddViewProp(imageActor); +// renderer->AddActor(imageActor); +// +// renderWindow->Render(); +// renderWindowInteractor->Start(); +// } +// +// +// static std::vector InitializeTestCases() +// { +// int testcounter = 0; +// std::vector tests= +// +// //#BEGIN setup TestCases +// +// { +// { +// ++testcounter, +// "TestCoronal", +// "image.nrrd", +// "coronalReference.nrrd", +// false, +// &TestCoronal +// }, +// { +// ++testcounter, +// "TestSagital", +// "image.nrrd", +// "sagitalReference.nrrd", +// false, +// &TestSagital +// }, +// { +// ++testcounter, +// "TestCoronal", +// "image.nrrd", +// "coronalReference.nrrd", +// false, +// &TestCoronal +// }, +// { +// ++testcounter, +// "Test_u_Rotation", +// "image.nrrd", +// "uRotationReference.nrrd", +// false, +// &Test_u_Rotation +// }, +// { +// ++testcounter, +// "Test_v_Rotation", +// "image.nrrd", +// "vRotationReference.nrrd", +// false, +// &Test_v_Rotation +// }, +// { +// ++testcounter, +// "TestTwoDirectionalRation", +// "image.nrrd", +// "twoDirectionalRationReference.nrrd", +// false, +// &TestTwoDirectionalRotation +// }, +// { +// ++testcounter, +// "Test4D", +// "image.nrrd", +// "twoDirectionalRationReference.nrrd", +// false, +// &Test4D +// }, +// { +// ++testcounter, +// "Test2D", +// "coronalReference.nrrd", +// "coronalReference.nrrd", +// false, +// &Test2D +// }, +// { +// ++testcounter, +// "Test2D", +// NULL, +// NULL, +// false, +// &Test1D +// } +// +// }; +// +// //#END setup TestCases +// +// return tests; +// } +// +//protected: +// +// static mitk::Image::Pointer GetImageToTest(std::string filename){ +// //retrieve referenceImage +// +//// mitk::StandardFileLocations::Pointer locator = mitk::StandardFileLocations::GetInstance(); +//// +//// const std::string filepath = locator->FindFile(filename, "Modules/MitkExt/Testing/Data"); +//// +//// if (filepath.empty()) +//// { +//// return NULL; +//// } +//// +//////TODO read imge from file +//// itk::FilenamesContainer file; +//// file.push_back( filename ); +// mitk::ItkImageFileReader::Pointer reader = mitk::ItkImageFileReader::New(); +// reader->SetFileName("C:\home\Pics\Pic3D.nrrd"); +// +// reader->Update(); +// +// mitk::Image::Pointer image = reader->GetOutput(); +// +// return image; +// } +// +// +// static mitk::Geometry2D::Pointer TestSagital() +// { +// +// return NULL; +// } +// +// static mitk::Geometry2D::Pointer TestCoronal() +// { +//return NULL; +// } +// +// static mitk::Geometry2D::Pointer TestTransversal() +// { +//return NULL; +// } +// +// static mitk::Geometry2D::Pointer Test_u_Rotation() +// { +//return NULL; +// } +// +// static mitk::Geometry2D::Pointer Test_v_Rotation() +// { +//return NULL; +// } +// +// static mitk::Geometry2D::Pointer TestTwoDirectionalRotation() +// { +//return NULL; +// } +// +// static mitk::Geometry2D::Pointer Test4DImage() +// {return NULL; +// +// } +// +// static mitk::Geometry2D::Pointer Test2DImage() +// { +//return NULL; +// } +// +// static mitk::Geometry2D::Pointer Test1DImage() +// { +//return NULL; +// } +// +//}; +// +// +//** +// * Tests for the class "mitkExtractDirectedPlaneImageFilter". +// * +// * argc and argv are the command line parameters which were passed to +// * the ADD_TEST command in the CMakeLists.txt file. For the automatic +// * tests, argv is either empty for the simple tests or contains the filename +// * of a test image for the image tests (see CMakeLists.txt). +// */ +//int mitkExtractDirectedPlaneImageFilterTest(int /* argc */, char* /*argv*/[]) +//{ +// // always start with this! +// MITK_TEST_BEGIN("mitkExtractDirectedPlaneImageFilter") +// +// +// mitk::ExtractDirectedPlaneImageFilter::Pointer extractor = mitk::ExtractDirectedPlaneImageFilter::New(); +// MITK_TEST_CONDITION_REQUIRED(extractor.IsNotNull(),"Testing instantiation") +// +// +// std::vector allTests = ExtractionTesting::InitializeTestCases(); +// +// for( int i = 0; i < allTests.size(); i++);{ +// +// ExtractionTesting::Testcase testCase = allTest[i]; +// +// ExtractionTesting::DoTesting(testCase); +// +// MITK_TEST_CONDITION(testCase.success, "Testcase #'" << testCase.number << " " << testCase.name); +// } +// +// // always end with this! +// MITK_TEST_END() +//} \ No newline at end of file diff --git a/Modules/ImageExtraction/files.cmake b/Modules/ImageExtraction/files.cmake index 86ec6e7d9c..ebb2d57d62 100644 --- a/Modules/ImageExtraction/files.cmake +++ b/Modules/ImageExtraction/files.cmake @@ -1,11 +1,12 @@ set(CPP_FILES mitkExtractDirectedPlaneImageFilter.cpp mitkExtractDirectedPlaneImageFilterNew.cpp mitkExtractImageFilter.cpp + #mitkExtractSliceFilter.cpp ) diff --git a/Modules/ImageExtraction/mitkExtractDirectedPlaneImageFilter.cpp b/Modules/ImageExtraction/mitkExtractDirectedPlaneImageFilter.cpp index 6fe33c2802..7a910aec58 100644 --- a/Modules/ImageExtraction/mitkExtractDirectedPlaneImageFilter.cpp +++ b/Modules/ImageExtraction/mitkExtractDirectedPlaneImageFilter.cpp @@ -1,498 +1,509 @@ /*=================================================================== 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) { m_Reslicer = vtkImageReslice::New(); + m_TargetTimestep = 0; - m_InPlaneResampleExtentByGeometry = false; + m_InPlaneResampleExtentByGeometry = true; + m_ResliceInterpolationProperty = NULL;//VtkResliceInterpolationProperty::New(); //TODO initial with value + m_ThickSlicesMode = 0; + m_ThickSlicesNum = 1; } mitk::ExtractDirectedPlaneImageFilter::~ExtractDirectedPlaneImageFilter() { - m_WorldGeometry = NULL; + 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(); + Superclass::GenerateOutputInformation(); } bool mitk::ExtractDirectedPlaneImageFilter ::CalculateClippedPlaneBounds( const Geometry3D *boundingGeometry, - const PlaneGeometry *planeGeometry, vtkFloatingPointType *bounds ) + 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 *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 2d51c72c9d..8161b01819 100644 --- a/Modules/ImageExtraction/mitkExtractDirectedPlaneImageFilter.h +++ b/Modules/ImageExtraction/mitkExtractDirectedPlaneImageFilter.h @@ -1,94 +1,121 @@ /*=================================================================== 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 { /** \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; + 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/Segmentation/Interactions/mitkSegTool2D.cpp b/Modules/Segmentation/Interactions/mitkSegTool2D.cpp index c3e41bca80..f5c4e58c22 100644 --- a/Modules/Segmentation/Interactions/mitkSegTool2D.cpp +++ b/Modules/Segmentation/Interactions/mitkSegTool2D.cpp @@ -1,402 +1,402 @@ /*=================================================================== 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 "mitkSegTool2D.h" #include "mitkToolManager.h" #include "mitkDataStorage.h" #include "mitkBaseRenderer.h" #include "mitkPlaneGeometry.h" #include "mitkExtractImageFilter.h" #include "mitkExtractDirectedPlaneImageFilter.h" //Include of the new ImageExtractor #include "mitkExtractDirectedPlaneImageFilterNew.h" #include "mitkPlanarCircle.h" #include "mitkOverwriteSliceImageFilter.h" #include "mitkOverwriteDirectedPlaneImageFilter.h" #include "mitkGetModuleContext.h" //Includes for 3DSurfaceInterpolation #include "mitkImageToContourFilter.h" #include "mitkSurfaceInterpolationController.h" #define ROUND(a) ((a)>0 ? (int)((a)+0.5) : -(int)(0.5-(a))) mitk::SegTool2D::SegTool2D(const char* type) :Tool(type), m_LastEventSender(NULL), m_LastEventSlice(0), m_Contourmarkername ("Position"), m_ShowMarkerNodes (true), m_3DInterpolationEnabled (true) { // great magic numbers CONNECT_ACTION( 80, OnMousePressed ); CONNECT_ACTION( 90, OnMouseMoved ); CONNECT_ACTION( 42, OnMouseReleased ); CONNECT_ACTION( 49014, OnInvertLogic ); } mitk::SegTool2D::~SegTool2D() { } bool mitk::SegTool2D::OnMousePressed (Action*, const StateEvent* stateEvent) { const PositionEvent* positionEvent = dynamic_cast(stateEvent->GetEvent()); if (!positionEvent) return false; if ( positionEvent->GetSender()->GetMapperID() != BaseRenderer::Standard2D ) return false; // we don't want anything but 2D m_LastEventSender = positionEvent->GetSender(); m_LastEventSlice = m_LastEventSender->GetSlice(); return true; } bool mitk::SegTool2D::OnMouseMoved (Action*, const StateEvent* stateEvent) { const PositionEvent* positionEvent = dynamic_cast(stateEvent->GetEvent()); if (!positionEvent) return false; if ( m_LastEventSender != positionEvent->GetSender() ) return false; if ( m_LastEventSlice != m_LastEventSender->GetSlice() ) return false; return true; } bool mitk::SegTool2D::OnMouseReleased(Action*, const StateEvent* stateEvent) { const PositionEvent* positionEvent = dynamic_cast(stateEvent->GetEvent()); if (!positionEvent) return false; if ( m_LastEventSender != positionEvent->GetSender() ) return false; if ( m_LastEventSlice != m_LastEventSender->GetSlice() ) return false; return true; } bool mitk::SegTool2D::OnInvertLogic(Action*, const StateEvent* stateEvent) { const PositionEvent* positionEvent = dynamic_cast(stateEvent->GetEvent()); if (!positionEvent) return false; if ( m_LastEventSender != positionEvent->GetSender() ) return false; if ( m_LastEventSlice != m_LastEventSender->GetSlice() ) return false; return true; } bool mitk::SegTool2D::DetermineAffectedImageSlice( const Image* image, const PlaneGeometry* plane, int& affectedDimension, int& affectedSlice ) { assert(image); assert(plane); // compare normal of plane to the three axis vectors of the image Vector3D normal = plane->GetNormal(); Vector3D imageNormal0 = image->GetSlicedGeometry()->GetAxisVector(0); Vector3D imageNormal1 = image->GetSlicedGeometry()->GetAxisVector(1); Vector3D imageNormal2 = image->GetSlicedGeometry()->GetAxisVector(2); normal.Normalize(); imageNormal0.Normalize(); imageNormal1.Normalize(); imageNormal2.Normalize(); imageNormal0.Set_vnl_vector( vnl_cross_3d(normal.Get_vnl_vector(),imageNormal0.Get_vnl_vector()) ); imageNormal1.Set_vnl_vector( vnl_cross_3d(normal.Get_vnl_vector(),imageNormal1.Get_vnl_vector()) ); imageNormal2.Set_vnl_vector( vnl_cross_3d(normal.Get_vnl_vector(),imageNormal2.Get_vnl_vector()) ); double eps( 0.00001 ); // transversal if ( imageNormal2.GetNorm() <= eps ) { affectedDimension = 2; } // sagittal else if ( imageNormal1.GetNorm() <= eps ) { affectedDimension = 1; } // frontal else if ( imageNormal0.GetNorm() <= eps ) { affectedDimension = 0; } else { affectedDimension = -1; // no idea return false; } // determine slice number in image Geometry3D* imageGeometry = image->GetGeometry(0); Point3D testPoint = imageGeometry->GetCenter(); Point3D projectedPoint; plane->Project( testPoint, projectedPoint ); Point3D indexPoint; imageGeometry->WorldToIndex( projectedPoint, indexPoint ); affectedSlice = ROUND( indexPoint[affectedDimension] ); MITK_DEBUG << "indexPoint " << indexPoint << " affectedDimension " << affectedDimension << " affectedSlice " << affectedSlice; // check if this index is still within the image if ( affectedSlice < 0 || affectedSlice >= static_cast(image->GetDimension(affectedDimension)) ) return false; return true; } mitk::Image::Pointer mitk::SegTool2D::GetAffectedImageSliceAs2DImage(const PositionEvent* positionEvent, const Image* image) { if (!positionEvent) return NULL; assert( positionEvent->GetSender() ); // sure, right? unsigned int timeStep = positionEvent->GetSender()->GetTimeStep( image ); // get the timestep of the visible part (time-wise) of the image // first, we determine, which slice is affected const PlaneGeometry* planeGeometry( dynamic_cast (positionEvent->GetSender()->GetCurrentWorldGeometry2D() ) ); if ( !image || !planeGeometry ) return NULL; int affectedDimension( -1 ); int affectedSlice( -1 ); - DetermineAffectedImageSlice( image, planeGeometry, affectedDimension, affectedSlice ); + //DetermineAffectedImageSlice( image, planeGeometry, affectedDimension, affectedSlice ); if ( DetermineAffectedImageSlice( image, planeGeometry, affectedDimension, affectedSlice ) ) { try { // now we extract the correct slice from the volume, resulting in a 2D image ExtractImageFilter::Pointer extractor= ExtractImageFilter::New(); extractor->SetInput( image ); extractor->SetSliceDimension( affectedDimension ); extractor->SetSliceIndex( affectedSlice ); extractor->SetTimeStep( timeStep ); extractor->Update(); // here we have a single slice that can be modified Image::Pointer slice = extractor->GetOutput(); //Workaround because of bug #7079 Point3D origin = slice->GetGeometry()->GetOrigin(); int affectedDimension(-1); if(positionEvent->GetSender()->GetRenderWindow() == mitk::BaseRenderer::GetRenderWindowByName("stdmulti.widget1")) { affectedDimension = 2; } if(positionEvent->GetSender()->GetRenderWindow() == mitk::BaseRenderer::GetRenderWindowByName("stdmulti.widget2")) { affectedDimension = 0; } if(positionEvent->GetSender()->GetRenderWindow() == mitk::BaseRenderer::GetRenderWindowByName("stdmulti.widget3")) { affectedDimension = 1; } if (affectedDimension != -1) { origin[affectedDimension] = planeGeometry->GetOrigin()[affectedDimension]; slice->GetGeometry()->SetOrigin(origin); } //Workaround end return slice; } catch(...) { // not working return NULL; } } else { ExtractDirectedPlaneImageFilterNew::Pointer newExtractor = ExtractDirectedPlaneImageFilterNew::New(); newExtractor->SetInput( image ); newExtractor->SetActualInputTimestep( timeStep ); newExtractor->SetCurrentWorldGeometry2D( planeGeometry ); newExtractor->Update(); Image::Pointer slice = newExtractor->GetOutput(); return slice; } } mitk::Image::Pointer mitk::SegTool2D::GetAffectedWorkingSlice(const PositionEvent* positionEvent) { DataNode* workingNode( m_ToolManager->GetWorkingData(0) ); if ( !workingNode ) return NULL; Image* workingImage = dynamic_cast(workingNode->GetData()); if ( !workingImage ) return NULL; return GetAffectedImageSliceAs2DImage( positionEvent, workingImage ); } mitk::Image::Pointer mitk::SegTool2D::GetAffectedReferenceSlice(const PositionEvent* positionEvent) { DataNode* referenceNode( m_ToolManager->GetReferenceData(0) ); if ( !referenceNode ) return NULL; Image* referenceImage = dynamic_cast(referenceNode->GetData()); if ( !referenceImage ) return NULL; return GetAffectedImageSliceAs2DImage( positionEvent, referenceImage ); } void mitk::SegTool2D::WriteBackSegmentationResult (const PositionEvent* positionEvent, Image* slice) { const PlaneGeometry* planeGeometry( dynamic_cast (positionEvent->GetSender()->GetCurrentWorldGeometry2D() ) ); DataNode* workingNode( m_ToolManager->GetWorkingData(0) ); Image* image = dynamic_cast(workingNode->GetData()); int affectedDimension( -1 ); int affectedSlice( -1 ); DetermineAffectedImageSlice( image, planeGeometry, affectedDimension, affectedSlice ); if (affectedDimension != -1) { OverwriteSliceImageFilter::Pointer slicewriter = OverwriteSliceImageFilter::New(); slicewriter->SetInput( image ); slicewriter->SetCreateUndoInformation( true ); slicewriter->SetSliceImage( slice ); slicewriter->SetSliceDimension( affectedDimension ); slicewriter->SetSliceIndex( affectedSlice ); slicewriter->SetTimeStep( positionEvent->GetSender()->GetTimeStep( image ) ); slicewriter->Update(); } else { OverwriteDirectedPlaneImageFilter::Pointer slicewriter = OverwriteDirectedPlaneImageFilter::New(); slicewriter->SetInput( image ); slicewriter->SetCreateUndoInformation( false ); slicewriter->SetSliceImage( slice ); slicewriter->SetPlaneGeometry3D( slice->GetGeometry() ); slicewriter->SetTimeStep( positionEvent->GetSender()->GetTimeStep( image ) ); slicewriter->Update(); } if ( m_3DInterpolationEnabled ) { slice->DisconnectPipeline(); ImageToContourFilter::Pointer contourExtractor = ImageToContourFilter::New(); contourExtractor->SetInput(slice); contourExtractor->Update(); mitk::Surface::Pointer contour = contourExtractor->GetOutput(); if (contour->GetVtkPolyData()->GetNumberOfPoints() > 0 ) { unsigned int pos = this->AddContourmarker(positionEvent); mitk::ServiceReference serviceRef = mitk::GetModuleContext()->GetServiceReference(); PlanePositionManagerService* service = dynamic_cast(mitk::GetModuleContext()->GetService(serviceRef)); mitk::SurfaceInterpolationController::GetInstance()->AddNewContour( contour, service->GetPlanePosition(pos)); contour->DisconnectPipeline(); } } } void mitk::SegTool2D::SetShowMarkerNodes(bool status) { m_ShowMarkerNodes = status; } void mitk::SegTool2D::Enable3DInterpolation(bool status) { m_3DInterpolationEnabled = status; } unsigned int mitk::SegTool2D::AddContourmarker ( const PositionEvent* positionEvent ) { const mitk::Geometry2D* plane = dynamic_cast (dynamic_cast< const mitk::SlicedGeometry3D*>( positionEvent->GetSender()->GetSliceNavigationController()->GetCurrentGeometry3D())->GetGeometry2D(0)); mitk::ServiceReference serviceRef = mitk::GetModuleContext()->GetServiceReference(); PlanePositionManagerService* service = dynamic_cast(mitk::GetModuleContext()->GetService(serviceRef)); unsigned int size = service->GetNumberOfPlanePositions(); unsigned int id = service->AddNewPlanePosition(plane, positionEvent->GetSender()->GetSliceNavigationController()->GetSlice()->GetPos()); mitk::PlanarCircle::Pointer contourMarker = mitk::PlanarCircle::New(); contourMarker->SetGeometry2D( const_cast(plane)); std::stringstream markerStream; mitk::DataNode* workingNode (m_ToolManager->GetWorkingData(0)); markerStream << m_Contourmarkername ; markerStream << " "; markerStream << id+1; DataNode::Pointer rotatedContourNode = DataNode::New(); rotatedContourNode->SetData(contourMarker); rotatedContourNode->SetProperty( "name", StringProperty::New(markerStream.str()) ); rotatedContourNode->SetProperty( "isContourMarker", BoolProperty::New(true)); rotatedContourNode->SetBoolProperty( "PlanarFigureInitializedWindow", true, positionEvent->GetSender() ); rotatedContourNode->SetProperty( "includeInBoundingBox", BoolProperty::New(false)); rotatedContourNode->SetProperty( "helper object", mitk::BoolProperty::New(!m_ShowMarkerNodes)); if (plane) { if ( id == size ) { m_ToolManager->GetDataStorage()->Add(rotatedContourNode, workingNode); } else { mitk::NodePredicateProperty::Pointer isMarker = mitk::NodePredicateProperty::New("isContourMarker", mitk::BoolProperty::New(true)); mitk::DataStorage::SetOfObjects::ConstPointer markers = m_ToolManager->GetDataStorage()->GetDerivations(workingNode,isMarker); for ( mitk::DataStorage::SetOfObjects::const_iterator iter = markers->begin(); iter != markers->end(); ++iter) { std::string nodeName = (*iter)->GetName(); unsigned int t = nodeName.find_last_of(" "); unsigned int markerId = atof(nodeName.substr(t+1).c_str())-1; if(id == markerId) { return id; } } m_ToolManager->GetDataStorage()->Add(rotatedContourNode, workingNode); } } return id; } void mitk::SegTool2D::InteractiveSegmentationBugMessage( const std::string& message ) { MITK_ERROR << "********************************************************************************" << std::endl << " " << message << std::endl << "********************************************************************************" << std::endl << " " << std::endl << " If your image is rotated or the 2D views don't really contain the patient image, try to press the button next to the image selection. " << std::endl << " " << std::endl << " Please file a BUG REPORT: " << std::endl << " http://bugs.mitk.org" << std::endl << " Contain the following information:" << std::endl << " - What image were you working on?" << std::endl << " - Which region of the image?" << std::endl << " - Which tool did you use?" << std::endl << " - What did you do?" << std::endl << " - What happened (not)? What did you expect?" << std::endl; }