diff --git a/Core/Code/Algorithms/mitkExtractSliceFilter.cpp b/Core/Code/Algorithms/mitkExtractSliceFilter.cpp index 7caf170b13..c74ef1a934 100644 --- a/Core/Code/Algorithms/mitkExtractSliceFilter.cpp +++ b/Core/Code/Algorithms/mitkExtractSliceFilter.cpp @@ -1,175 +1,487 @@ /*========================================================================= 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 + mitk::ExtractSliceFilter::ExtractSliceFilter(){ m_Reslicer = vtkSmartPointer::New(); m_TimeStep = 0; - m_Reslicer->ReleaseDataFlagOn();//todo check - interpolationMode = ExtractSliceFilter::ResliceInterpolation::RESLICE_NEAREST; + m_Reslicer->ReleaseDataFlagOn();//todo check what it does + m_InterpolationMode = ExtractSliceFilter::RESLICE_NEAREST; + m_ResliceTransform = NULL; + m_InPlaneResampleExtentByGeometry = false; } mitk::ExtractSliceFilter::~ExtractSliceFilter(){ } +void mitk::ExtractSliceFilter::GenerateOutputInformation(){ + + /*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(){ + ImageToImageFilter::InputImagePointer input = const_cast< ImageToImageFilter::InputImageType* > ( this->GetInput() ); + input->SetRequestedRegionToLargestPossibleRegion(); +} -void mitk::ExtractSliceFilter::GenerateData(){ +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; } - - - if(input->GetDimension() == 2){ - //if the image is already 2D, we have nothing to resice and we can set the input image as the output - - /*Image::Pointer resultImage = ImageToImageFilter::GetOutput(); - resultImage = const_cast(input); - ImageToImageFilter::SetNthOutput( 0, resultImage ); - return;*/ - } - const TimeSlicedGeometry *inputTimeGeometry = this->GetInput()->GetTimeSlicedGeometry(); if ( ( inputTimeGeometry == NULL ) || ( inputTimeGeometry->GetTimeSteps() == 0 ) ) { itkWarningMacro(<<"Error reading input image TimeSlicedGeometry."); return; } 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================*/ - m_Reslicer->SetInput(input->GetVtkImageData(m_TimeStep)); + + Point3D origin; + Vector3D right, bottom, normal; + double widthInMM, heightInMM; + Vector2D extent; + double mmPerPixel[2]; + + + const PlaneGeometry* planeGeometry = dynamic_cast(m_WorldGeometry); + if ( planeGeometry != NULL ) + { + origin = planeGeometry->GetOrigin(); + right = planeGeometry->GetAxisVector( 0 ); // right = Extent of Image in mm (worldspace) + 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 ); + + + mmPerPixel[0] = widthInMM / extent[0]; + mmPerPixel[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 * ( mmPerPixel[0] * 0.5 ); + origin += bottom * ( mmPerPixel[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) + 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 ); + + + // 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, sliceBounds ); + + } + ////Code for curved planes + // 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); + + // localStorage->m_mmPerPixel[0] = widthInMM / extent[0]; + // localStorage->m_mmPerPixel[1] = heightInMM / extent[1]; + + // localStorage->m_Origin = abstractGeometry->GetPlane()->GetOrigin(); + + // localStorage->m_Right = abstractGeometry->GetPlane()->GetAxisVector(0); + // localStorage->m_Right.Normalize(); + + // localStorage->m_Bottom = abstractGeometry->GetPlane()->GetAxisVector(1); + // localStorage->m_Bottom.Normalize(); + + // 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 ); + // } + //} + + vtkImageChangeInformation * unitSpacingImageFilter = vtkImageChangeInformation::New() ; + unitSpacingImageFilter->ReleaseDataFlagOn(); + + unitSpacingImageFilter->SetOutputSpacing( 1.0, 1.0, 1.0 ); + unitSpacingImageFilter->SetInput( input->GetVtkImageData(m_TimeStep) ); + + m_Reslicer->SetInput(unitSpacingImageFilter->GetOutput() ); + /*m_Reslicer->SetInput(input->GetVtkImageData(m_TimeStep));*/ + /*setup the plane where vktImageReslice extracts the slice*/ //ResliceAxesOrigin is the ancor point of the plane - double origin[3]; - itk2vtk(m_WorldGeometry->GetOrigin(), origin); - m_Reslicer->SetResliceAxesOrigin(origin); + double originArray[3]; + itk2vtk(origin, originArray); + m_Reslicer->SetResliceAxesOrigin(originArray); + //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]; - Vector3D right, bottom; - right = m_WorldGeometry->GetAxisVector(0); - right.Normalize(); - vnl2vtk(right.GetVnlVector(), cosines); + vnl2vtk(right.GetVnlVector(), cosines);//x - bottom = m_WorldGeometry->GetAxisVector(1); - bottom.Normalize(); - vnl2vtk(bottom.GetVnlVector(), cosines + 3); + vnl2vtk(bottom.GetVnlVector(), cosines + 3);//y Vector3D normalVector = CrossProduct(m_WorldGeometry->GetAxisVector(0), m_WorldGeometry->GetAxisVector(1)); normalVector.Normalize(); - vnl2vtk(normalVector.GetVnlVector(), cosines + 6); + vnl2vtk(normal.GetVnlVector(), cosines + 6);//n m_Reslicer->SetResliceAxesDirectionCosines(cosines); //we only have one slice, not a volume m_Reslicer->SetOutputDimensionality(2); //set the interpolation mode for slicing - switch(this->interpolationMode){ - case ResliceInterpolation::RESLICE_NEAREST: + switch(this->m_InterpolationMode){ + case RESLICE_NEAREST: m_Reslicer->SetInterpolationModeToNearestNeighbor(); break; - case ResliceInterpolation::RESLICE_LINEAR: + case RESLICE_LINEAR: m_Reslicer->SetInterpolationModeToLinear(); break; - case ResliceInterpolation::RESLICE_CUBIC: + case RESLICE_CUBIC: m_Reslicer->SetInterpolationModeToCubic(); default: //the default interpolation used by mitk m_Reslicer->SetInterpolationModeToNearestNeighbor(); } + + + int xMin, xMax, yMin, yMax; + bool boundsInitialized = false; + if(m_WorldGeometry->GetReferenceGeometry()){ + vtkFloatingPointType sliceBounds[6]; + + for ( int i = 0; i < 6; ++i ) + { + sliceBounds[i] = 0.0; + } + this->CalculateClippedPlaneBounds( m_WorldGeometry->GetReferenceGeometry(), planeGeometry, sliceBounds ); + + // Calculate output extent (integer values) + xMin = static_cast< int >( sliceBounds[0] / mmPerPixel[0] + 0.5 ); + xMax = static_cast< int >( sliceBounds[1] / mmPerPixel[0] + 0.5 ); + yMin = static_cast< int >( sliceBounds[2] / mmPerPixel[1] + 0.5 ); + yMax = static_cast< int >( sliceBounds[3] / 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]); + yMax = static_cast< int >( extent[1]); + } + + m_Reslicer->SetOutputExtent(xMin, xMax-1, yMin, yMax-1, 0, 0 ); + + + + m_Reslicer->SetOutputOrigin( 0.0, 0.0, 0.0 ); + m_Reslicer->SetOutputSpacing( mmPerPixel[0], mmPerPixel[1], 1.0 ); + //m_Reslicer->SetOutputSpacing(m_WorldGeometry->GetSpacing()[0],m_WorldGeometry->GetSpacing()[1], 1.0); + //m_Reslicer->SetOutputSpacing( 1.0, 1.0, 1.0 ); + + //start the pipeline m_Reslicer->Modified(); m_Reslicer->Update(); /*================ #END setup vtkImageRslice properties================*/ - /*================ #BEGIN Get the slice from vtkImageReslice and convert it to mitk 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 AffineGeometryFrame3D::Pointer originalGeometryAGF = m_WorldGeometry->Clone(); Geometry2D::Pointer originalGeometry = dynamic_cast( originalGeometryAGF.GetPointer() ); originalGeometry->ChangeImageGeometryConsideringOriginOffset(true); resultImage->SetGeometry( originalGeometry ); - //resultImage->GetGeometry()->TransferItkToVtkTransform(); + resultImage->GetGeometry()->TransferItkToVtkTransform(); resultImage->GetGeometry()->Modified(); /*================ #END Get the slice from vtkImageReslice and convert it to mitk Image================*/ +} + + +bool mitk::ExtractSliceFilter::GetBounds(vtkFloatingPointType bounds[6]){ + + return this->CalculateClippedPlaneBounds(m_WorldGeometry->GetReferenceGeometry(), dynamic_cast< const PlaneGeometry * >( m_WorldGeometry ), bounds); +} + +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 index ab9d20cb8f..312ebe7b37 100644 --- a/Core/Code/Algorithms/mitkExtractSliceFilter.h +++ b/Core/Code/Algorithms/mitkExtractSliceFilter.h @@ -1,70 +1,99 @@ /*========================================================================= 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. * */ class MITK_CORE_EXPORT ExtractSliceFilter : public ImageToImageFilter { public: mitkClassMacro(ExtractSliceFilter, ImageToImageFilter); itkNewMacro(ExtractSliceFilter); + int* GetPlaneBounds(){ return this->m_Reslicer->GetOutput()->GetWholeExtent();} + void SetWorldGeometry(const Geometry2D* geometry ){ this->m_WorldGeometry = geometry; this->Modified(); } + 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. + */ + void SetResliceTransformByGeometry(const Geometry3D* transform){ this->m_ResliceTransform = transform; } + + + void SetInPlaneResampleExtentByGeometry(bool inPlaneResampleExtentByGeometry){ this->m_InPlaneResampleExtentByGeometry = inPlaneResampleExtentByGeometry; } + + bool GetBounds(vtkFloatingPointType bounds[6]); + vtkMatrix4x4* GetResliceAxes(){ return this->m_Reslicer->GetResliceAxes(); } enum ResliceInterpolation { RESLICE_NEAREST, RESLICE_LINEAR, RESLICE_CUBIC }; - void SetInterPolationMode( ExtractSliceFilter::ResliceInterpolation interpolation){ this->interpolationMode = interpolation; } + void SetInterPolationMode( ExtractSliceFilter::ResliceInterpolation interpolation){ this->m_InterpolationMode = interpolation; } protected: ExtractSliceFilter(); virtual ~ExtractSliceFilter(); virtual void GenerateData(); + virtual void GenerateOutputInformation(); + virtual void GenerateInputRequestedRegion(); + const Geometry2D* m_WorldGeometry; vtkSmartPointer m_Reslicer; unsigned int m_TimeStep; - ResliceInterpolation interpolationMode; + ResliceInterpolation m_InterpolationMode; + + const Geometry3D* m_ResliceTransform; + + bool m_InPlaneResampleExtentByGeometry; + + bool CalculateClippedPlaneBounds( const Geometry3D *boundingGeometry, const PlaneGeometry *planeGeometry, vtkFloatingPointType *bounds ); + bool LineIntersectZero( vtkPoints *points, int p1, int p2, 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 1a1d6050ff..99f605beab 100644 --- a/Core/Code/Rendering/mitkImageVtkMapper2D.cpp +++ b/Core/Code/Rendering/mitkImageVtkMapper2D.cpp @@ -1,965 +1,761 @@ /*========================================================================= 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. =========================================================================*/ //MITK #include #include #include #include #include #include #include #include #include #include #include #include //#include #include //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 ) { return; } //check if there is a valid worldGeometry const Geometry2D *worldGeometry = renderer->GetCurrentWorldGeometry2D(); if( ( worldGeometry == NULL ) || ( !worldGeometry->IsValid() ) || ( !worldGeometry->HasReferenceGeometry() )) { return; } input->Update(); - //localStorage->m_Reslicer = mitk::ExtractSliceFilter::New(); localStorage->m_Reslicer->SetInput(input); localStorage->m_Reslicer->SetWorldGeometry(worldGeometry); - localStorage->m_Reslicer->Update(); + localStorage->m_Reslicer->SetResliceTransformByGeometry( input->GetTimeSlicedGeometry()->GetGeometry3D( this->GetTimestep() ) ); - - // 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() ); + bool inPlaneResampleExtentByGeometry = false; + GetDataNode()->GetBoolProperty("in plane resample extent by geometry", inPlaneResampleExtentByGeometry, renderer); + localStorage->m_Reslicer->SetInPlaneResampleExtentByGeometry(inPlaneResampleExtentByGeometry); + + localStorage->m_Reslicer->UpdateLargestPossibleRegion(); + // Bounds information for reslicing (only reuqired if reference geometry // is present) + MITK_INFO<<"ViewDirection: "<GetSliceNavigationController()->GetViewDirection()<<"World-Origin: "<GetOrigin(); 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]; - - // 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( - worldGeometry->GetReferenceGeometry(), planeGeometry, sliceBounds ); - } - - ////Code for curved planes - // 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); - - // localStorage->m_mmPerPixel[0] = widthInMM / extent[0]; - // localStorage->m_mmPerPixel[1] = heightInMM / extent[1]; - - // localStorage->m_Origin = abstractGeometry->GetPlane()->GetOrigin(); - - // localStorage->m_Right = abstractGeometry->GetPlane()->GetAxisVector(0); - // localStorage->m_Right.Normalize(); - - // localStorage->m_Bottom = abstractGeometry->GetPlane()->GetAxisVector(1); - // localStorage->m_Bottom.Normalize(); - - // 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 ); - // } - //} + localStorage->m_Reslicer->GetBounds(sliceBounds); + localStorage->m_ReslicedImage = localStorage->m_Reslicer->GetOutput()->GetVtkImageData(); //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()); + 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)); 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)); 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")); //...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 )); 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->GetScalarValueMin(); maxValue = centralSliceImage->GetScalarValueMax(); min2ndValue = centralSliceImage->GetScalarValue2ndMin(); max2ndValue = centralSliceImage->GetScalarValue2ndMax(); } if ( minValue == maxValue ) { // centralSlice is strange, lets look at all data minValue = image->GetScalarValueMin(); maxValue = image->GetScalarValueMaxNoRecompute(); min2ndValue = image->GetScalarValue2ndMinNoRecompute(); max2ndValue = image->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().GetItkTypeId() && *(image->GetPixelType().GetItkTypeId()) == 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); int* dims = localStorage->m_ReslicedImage->GetDimensions(); //dimensions of the image int line = dims[0]; //how many pixels per line? int x = 0; //pixel index x int y = 0; //pixel index y char* currentPixel; int nn = dims[0]*dims[1]; //max pixel(n,n) //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 for (int ii = 0; ii(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 if (ii >= line && *(currentPixel-line) == 0) { //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 (ii <= nn-line && *(currentPixel+line) == 0) { //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 (ii > 1 && *(currentPixel-1) == 0) { //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 (ii < nn-1 && *(currentPixel+1) == 0) { //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); } } //reached end of line x++; if (x >= line) { x = 0; y++; } } // 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); + //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 = 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_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/Modules/ImageExtraction/Testing/mitkExtractSliceFilterTest.cpp b/Modules/ImageExtraction/Testing/mitkExtractSliceFilterTest.cpp index 3b596ea19b..ffcca9bec8 100644 --- a/Modules/ImageExtraction/Testing/mitkExtractSliceFilterTest.cpp +++ b/Modules/ImageExtraction/Testing/mitkExtractSliceFilterTest.cpp @@ -1,1002 +1,1020 @@ /*========================================================================= 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 //#define CREATE_VOLUME //#define CALC_TESTFAILURE_DEVIATION //#define SHOW_SLICE_IN_RENDER_WINDOW #define 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, "resliced Image returned"); 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;} //do think about the deviation MITK_TEST_CONDITION(std::abs(100 - testResults.percentageAreaCalcToPixel) < (std::sqrt(devArea)), "testing area"); MITK_TEST_CONDITION(std::abs(100 - testResults.percentageRadiusToPixel) < (2 * devDiameter), "testing diameter"); #ifdef 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 InputImageType::IndexType currentIndexX; currentIndexX[0] = (int)(TestvolumeSize / 2.0); currentIndexX[1] = 0; 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); - int sumpixels = 0; + 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 TestPerPixelOrthogonal() { /* 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();*/ //set the seed of the rand function srand((unsigned)time(0)); /* setup a random orthogonal plane */ int sliceindex = rand() % 32; int planeOrientation = rand() % 3; bool isFrontside = true; bool isRotated = false; //sagittal mitk::PlaneGeometry::PlaneOrientation orientation = mitk::PlaneGeometry::PlaneOrientation::Sagittal; if( planeOrientation == 0) //transversal { isFrontside = false; isRotated = true; orientation = mitk::PlaneGeometry::PlaneOrientation::Transversal; }else if( planeOrientation == 2) //coronal { orientation = mitk::PlaneGeometry::PlaneOrientation::Frontal; } mitk::PlaneGeometry::Pointer plane = mitk::PlaneGeometry::New(); - // Maybe have a look at this method.. is it reaaaally correct? + plane->InitializeStandardPlane(imageInMitk->GetGeometry(), orientation, sliceindex, isFrontside, isRotated); plane->ChangeImageGeometryConsideringOriginOffset(true); mitk::Vector3D rotationVector; rotationVector[0] = randFloat(); rotationVector[1] = randFloat(); rotationVector[2] = randFloat(); float degree = randFloat() * 180.0; //we dont need this any more, because we are only testing orthogonal planes /*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 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 = rand() % 32; unsigned char n2 = 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, pp; p[0] = testPoint2DInIndex[0]; p[1] = testPoint2DInIndex[1]; p[2] = testPoint2DInIndex[2]; slice->GetGeometry()->IndexToWorld(p, pp); MITK_INFO << pp; p[0] = testPoint3DInIndex[0]; p[1] = testPoint3DInIndex[1]; p[2] = testPoint3DInIndex[2]; imageInMitk->GetGeometry()->IndexToWorld(p, pp); MITK_INFO << pp; //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"); #ifdef DEBUG MITK_INFO << "index: " << sliceindex << " orientation: " << planeOrientation << " frontside: " << isFrontside << " rotated: " << isRotated; MITK_INFO << "\n" << "point" << testPoint3DInWorld << " is " << testPoint2DInIndex << " in 2D"; MITK_INFO << "\n" << "inside plane: " << plane->IsInside(testPoint3DInWorld) << " and volume: " << imageInMitk->GetGeometry()->IsInside(testPoint3DInWorld); MITK_INFO << "\n" << "volume idx: " << testPoint3DInIndex << " = " << valueAt3DVolume ; MITK_INFO << "\n" << "volume w: " << 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 //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(); //save in file mitk::ImageWriter::Pointer writer = mitk::ImageWriter::New(); writer->SetInput(TestVolume); std::string file;/* = "C:\\home\\schroedt\\MITK\\Modules\\ImageExtraction\\Testing\\Data\\sphere_"; file += TestvolumeSize; file += ".nrrd";*/ 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 #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 = 512.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::TestPerPixelOrthogonal(); //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(256.0, 256.0, spacing, mitk::PlaneGeometry::PlaneOrientation::Transversal, sphereCenter, false, true); + geometryTransversal->InitializeStandardPlane(planeSize, planeSize, spacing, mitk::PlaneGeometry::PlaneOrientation::Transversal, sphereCenter, false, true); geometryTransversal->ChangeImageGeometryConsideringOriginOffset(true); mitkExtractSliceFilterTestClass::TestSlice(geometryTransversal, "Testing transversal plane"); /* end transversal plane */ /* sagittal plane */ mitk::PlaneGeometry::Pointer geometrySagital = mitk::PlaneGeometry::New(); - geometrySagital->InitializeStandardPlane(256.0, 256.0, spacing, mitk::PlaneGeometry::PlaneOrientation::Sagittal, sphereCenter, true, false); + geometrySagital->InitializeStandardPlane(planeSize, planeSize, spacing, mitk::PlaneGeometry::PlaneOrientation::Sagittal, sphereCenter, true, false); geometrySagital->ChangeImageGeometryConsideringOriginOffset(true); mitkExtractSliceFilterTestClass::TestSlice(geometrySagital, "Testing sagittal plane"); /* sagittal plane */ /* sagittal shifted plane */ mitk::PlaneGeometry::Pointer geometrySagitalShifted = mitk::PlaneGeometry::New(); - geometrySagitalShifted->InitializeStandardPlane(256.0, 256.0, spacing, mitk::PlaneGeometry::PlaneOrientation::Sagittal, (sphereCenter - 14), true, false); + geometrySagitalShifted->InitializeStandardPlane(planeSize, planeSize, spacing, mitk::PlaneGeometry::PlaneOrientation::Sagittal, (sphereCenter - 14), true, false); geometrySagitalShifted->ChangeImageGeometryConsideringOriginOffset(true); mitkExtractSliceFilterTestClass::TestSlice(geometrySagitalShifted, "Testing sagittal plane shifted"); /* end sagittal shifted plane */ /* coronal plane */ mitk::PlaneGeometry::Pointer geometryCoronal = mitk::PlaneGeometry::New(); - geometryCoronal->InitializeStandardPlane(256.0, 256.0, spacing, mitk::PlaneGeometry::PlaneOrientation::Frontal, sphereCenter, true, false); + geometryCoronal->InitializeStandardPlane(planeSize, planeSize, spacing, mitk::PlaneGeometry::PlaneOrientation::Frontal, sphereCenter, true, false); geometryCoronal->ChangeImageGeometryConsideringOriginOffset(true); mitkExtractSliceFilterTestClass::TestSlice(geometryCoronal, "Testing coronal plane"); /* end coronal plane */ /* oblique plane */ mitk::PlaneGeometry::Pointer obliquePlane = mitk::PlaneGeometry::New(); - obliquePlane->InitializeStandardPlane(256.0, 256.0, spacing, mitk::PlaneGeometry::PlaneOrientation::Sagittal, sphereCenter, true, false); + obliquePlane->InitializeStandardPlane(planeSize, planeSize, spacing, mitk::PlaneGeometry::PlaneOrientation::Sagittal, sphereCenter, true, false); obliquePlane->ChangeImageGeometryConsideringOriginOffset(true); mitk::Vector3D rotationVector; rotationVector[0] = mitkExtractSliceFilterTestClass::randFloat(); rotationVector[1] = mitkExtractSliceFilterTestClass::randFloat(); rotationVector[2] = mitkExtractSliceFilterTestClass::randFloat(); float degree = mitkExtractSliceFilterTestClass::randFloat() * 180; 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(); + mitk::ItkImageFileReader::Pointer reader = mitk::ItkImageFileReader::New(); std::string filename = "C:\\home\\schroedt\\MITK\\Modules\\ImageExtraction\\Testing\\Data\\Pic3D.nrrd"; reader->SetFileName(filename); reader->Update(); - mitk::Image::Pointer pic = reader->GetOutput();*/ + mitk::Image::Pointer pic = reader->GetOutput(); mitk::ExtractSliceFilter::Pointer slicer = mitk::ExtractSliceFilter::New(); - slicer->SetInput(mitkExtractSliceFilterTestClass::TestVolume); + slicer->SetInput(pic);//mitkExtractSliceFilterTestClass::TestVolume); + + + mitk::PlaneGeometry::Pointer obliquePl = mitk::PlaneGeometry::New(); + obliquePl->InitializeStandardPlane(pic->GetGeometry()->get, mitk::PlaneGeometry::PlaneOrientation::Sagittal, pic->GetGeometry()->GetCenter(),, true, false); + obliquePl->ChangeImageGeometryConsideringOriginOffset(true); + + mitk::Vector3D rotationVector; + rotationVector[0] = mitkExtractSliceFilterTestClass::randFloat(); + rotationVector[1] = mitkExtractSliceFilterTestClass::randFloat(); + rotationVector[2] = mitkExtractSliceFilterTestClass::randFloat(); + + float degree = mitkExtractSliceFilterTestClass::randFloat() * 180; + + mitk::RotationOperation* op = new mitk::RotationOperation(mitk::OpROTATE, obliquePl->GetCenter(), rotationVector, degree); + obliquePl->ExecuteOperation(op); + delete op; slicer->SetWorldGeometry(geometryCoronal); 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()->GetVtkImageData()); 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() } \ No newline at end of file