Index: Algorithms/mitkExtractDirectedPlaneImageFilter.cpp =================================================================== --- Algorithms/mitkExtractDirectedPlaneImageFilter.cpp (revision 20100) +++ Algorithms/mitkExtractDirectedPlaneImageFilter.cpp (working copy) @@ -25,6 +25,7 @@ #include #include #include +#include #include "pic2vtk.h" #include "picimage.h" @@ -341,40 +342,57 @@ m_Reslicer->ReleaseDataFlagOn(); m_Reslicer->Update(); + + // Over sampling the slice data to avoid missing some pixels in the extraction + // (otherwise results in stripes in the roi) + // original spacing + const double spacingX = m_Reslicer->GetOutputSpacing()[0]; + const double spacingY = m_Reslicer->GetOutputSpacing()[1]; + const double spacingZ = m_Reslicer->GetOutputSpacing()[2]; + // resampling factors + const double factorX = 2; + const double factorY = 2; + const double factorZ = 1; + // new spacing + const double newSpacingX = spacingX / factorX; + const double newSpacingY = spacingY / factorY; + const double newSpacingZ = spacingZ / factorZ; + + // resample + vtkImageResample* resample = vtkImageResample::New(); + resample->SetInput( m_Reslicer->GetOutput() ); + resample->SetDimensionality(3); + resample->SetAxisMagnificationFactor( 0, factorX ); + resample->SetAxisMagnificationFactor( 1, factorY ); + resample->SetAxisMagnificationFactor( 2, factorZ ); + resample->Update(); + // 1. Check the result - vtkImageData* reslicedImage = m_Reslicer->GetOutput(); + vtkImageData* outputImage = resample->GetOutput(); - mitkIpPicDescriptor *pic = Pic2vtk::convert( reslicedImage ); + mitkIpPicDescriptor *pic = Pic2vtk::convert( outputImage ); - if((reslicedImage == NULL) || (reslicedImage->GetDataDimension() < 1)) + if((outputImage == NULL) || (outputImage->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( pic ); - resultImage->SetSpacing( spacingVector ); resultImage->SetPicVolume( pic ); + // Create a new geometry to reflect our resampling + // The outputImage->GetGeometry() does not contain the plane transformation (not sure why?YM) + // copy transformation of the planeGeometry + resultImage->GetGeometry()->SetIndexToWorldTransform( const_cast(m_WorldGeometry->GetIndexToWorldTransform()) ); + // copy spacing of the outputImage (because of the resampling) + Vector3D spacingVector; + FillVector3D(spacingVector, newSpacingX, newSpacingY, newSpacingZ); + resultImage->GetGeometry()->SetSpacing( spacingVector ); + 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());*/ } Index: Algorithms/mitkOverwriteSliceImageFilter.cpp =================================================================== --- Algorithms/mitkOverwriteSliceImageFilter.cpp (revision 20100) +++ Algorithms/mitkOverwriteSliceImageFilter.cpp (working copy) @@ -15,6 +15,8 @@ =========================================================================*/ +#include + #include "mitkOverwriteSliceImageFilter.h" #include "mitkImageCast.h" #include "mitkSegmentationInterpolationController.h" @@ -25,6 +27,8 @@ #include "mitkDiffImageApplier.h" #include "mitkImageTimeSelector.h" +#include "mitkContourTool.h" + #include #include @@ -36,6 +40,7 @@ m_Dimension1(1), m_CreateUndoInformation(false) { + m_bSliceAndImageOriented = true; } mitk::OverwriteSliceImageFilter::~OverwriteSliceImageFilter() @@ -80,18 +85,35 @@ break; } - if ( slice->GetDimension() < 2 || input->GetDimension() > 4 || - slice->GetDimension(0) != input->GetDimension(m_Dimension0) || - slice->GetDimension(1) != input->GetDimension(m_Dimension1) || - m_SliceIndex >= input->GetDimension(m_SliceDimension) - ) + if ( slice->GetDimension() != 2 || ( input->GetDimension() != 3 && input->GetDimension() != 4 ) ) { itkExceptionMacro("Slice and image dimensions differ or slice index is too large. Sorry, cannot work like this."); return; } + + if ( slice->GetDimension(0) != input->GetDimension(m_Dimension0) || + slice->GetDimension(1) != input->GetDimension(m_Dimension1) || + m_SliceIndex >= input->GetDimension(m_SliceDimension) + ) + { + m_bSliceAndImageOriented = false; + } + else + { + m_bSliceAndImageOriented = true; + } if ( input->GetDimension() == 4 ) { + ImageTimeSelector::Pointer timeSelector = ImageTimeSelector::New(); + timeSelector->SetInput( input ); + timeSelector->SetTimeNr( m_TimeStep ); + timeSelector->UpdateLargestPossibleRegion(); + input3D = timeSelector->GetOutput(); + } + + if ( input->GetDimension() == 4 ) + { ImageTimeSelector::Pointer timeSelector = ImageTimeSelector::New(); timeSelector->SetInput( input ); timeSelector->SetTimeNr( m_TimeStep ); @@ -116,10 +138,13 @@ if (interpolator) { interpolator->BlockModified(true); - interpolator->SetChangedSlice( m_SliceDifferenceImage, m_SliceDimension, m_SliceIndex, m_TimeStep ); + if ( m_bSliceAndImageOriented ) + { + interpolator->SetChangedSlice( m_SliceDifferenceImage, m_SliceDimension, m_SliceIndex, m_TimeStep ); + } } - if ( m_CreateUndoInformation ) + if ( m_bSliceAndImageOriented && m_CreateUndoInformation ) { // create do/undo operations (we don't execute the doOp here, because it has already been executed during calculation of the diff image ApplyDiffImageOperation* doOp = new ApplyDiffImageOperation( OpTEST, const_cast(input.GetPointer()), m_SliceDifferenceImage, m_TimeStep, m_SliceDimension, m_SliceIndex ); @@ -182,16 +207,6 @@ typedef itk::ImageRegionConstIterator< SliceImageType > InputSliceIteratorType; typedef itk::ImageRegionIterator< DiffImageType > DiffSliceIteratorType; - typename VolumeImageType::RegionType sliceInVolumeRegion; - - sliceInVolumeRegion = outputImage->GetLargestPossibleRegion(); - sliceInVolumeRegion.SetSize( m_SliceDimension, 1 ); // just one slice - sliceInVolumeRegion.SetIndex( m_SliceDimension, m_SliceIndex ); // exactly this slice, please - - OutputSliceIteratorType outputIterator( outputImage, sliceInVolumeRegion ); - outputIterator.SetFirstDirection(m_Dimension0); - outputIterator.SetSecondDirection(m_Dimension1); - InputSliceIteratorType inputIterator( inputImage, inputImage->GetLargestPossibleRegion() ); typename DiffImageType::Pointer diffImage; @@ -199,25 +214,100 @@ DiffSliceIteratorType diffIterator( diffImage, diffImage->GetLargestPossibleRegion() ); // iterate over output slice (and over input slice simultaneously) - outputIterator.GoToBegin(); inputIterator.GoToBegin(); diffIterator.GoToBegin(); - while ( !outputIterator.IsAtEnd() ) + if ( m_bSliceAndImageOriented ) { - while ( !outputIterator.IsAtEndOfSlice() ) - { - while ( !outputIterator.IsAtEndOfLine() ) - { - diffIterator.Set( static_cast(inputIterator.Get() - outputIterator.Get()) ); // oh oh, not good for bigger values - outputIterator.Set( (TPixel2) inputIterator.Get() ); - ++outputIterator; - ++inputIterator; - ++diffIterator; - } - outputIterator.NextLine(); - } - outputIterator.NextSlice(); + typename VolumeImageType::RegionType sliceInVolumeRegion; + + sliceInVolumeRegion = outputImage->GetLargestPossibleRegion(); + sliceInVolumeRegion.SetSize( m_SliceDimension, 1 ); // just one slice + sliceInVolumeRegion.SetIndex( m_SliceDimension, m_SliceIndex ); // exactly this slice, please + + // Get an iterator to the output image slice + OutputSliceIteratorType outputIterator( outputImage, sliceInVolumeRegion ); + outputIterator.SetFirstDirection(m_Dimension0); + outputIterator.SetSecondDirection(m_Dimension1); + + outputIterator.GoToBegin(); + while ( !outputIterator.IsAtEnd() ) + { + while ( !outputIterator.IsAtEndOfSlice() ) + { + while ( !outputIterator.IsAtEndOfLine() ) + { + diffIterator.Set( static_cast(inputIterator.Get() - outputIterator.Get()) ); // oh oh, not good for bigger values + if( inputIterator.Get() == mitk::paint::addPixelValue ) + { + outputIterator.Set( (TPixel2)( 1 ) ); + } + else if( inputIterator.Get() == mitk::paint::subPixelValue ) + { + outputIterator.Set( (TPixel2)( 0 ) ); + } + ++outputIterator; + ++inputIterator; + ++diffIterator; + } + outputIterator.NextLine(); + } + outputIterator.NextSlice(); + } } + else + { + while ( !inputIterator.IsAtEnd() ) + { + // Input world point + Point3D currentPointIn2D; + currentPointIn2D.Fill(0); + for (int i = 0; i < 2; ++i) currentPointIn2D[i] = inputIterator.GetIndex()[i]; + Point3D worldPointIn3D; + worldPointIn3D.Fill(0); + m_SliceGeometry3D->IndexToWorld( currentPointIn2D, worldPointIn3D ); + + // Offset the world coordinate by one pixel to compensate for + // index/world origin differences. + Point3D offsetIndex; + offsetIndex.Fill( 1 ); + Point3D offsetWorld; + offsetWorld.Fill( 0 ); + m_SliceGeometry3D->IndexToWorld( offsetIndex, offsetWorld ); + // remove origin shift + const Point3D origin = m_SliceGeometry3D->GetOrigin(); + offsetWorld[0] -= origin[0]; + offsetWorld[1] -= origin[1]; + offsetWorld[2] -= origin[2]; + // offset world coordinate + worldPointIn3D[ 0 ] += offsetWorld[0]; + worldPointIn3D[ 1 ] += offsetWorld[1]; + worldPointIn3D[ 2 ] += offsetWorld[2]; + + // Output index + typename itk::Image::IndexType outputIndex; + m_ImageGeometry3D->WorldToIndex( worldPointIn3D, outputIndex ); + + // Only access ITK image if it's inside + TPixel2 outputPixel = 0; + if ( m_ImageGeometry3D->IsIndexInside( outputIndex ) ) + { + outputPixel = outputImage->GetPixel( outputIndex ); + if( inputIterator.Get() == mitk::paint::addPixelValue ) + { + outputImage->SetPixel( outputIndex, (TPixel2)( 1 ) ); + } + else if( inputIterator.Get() == mitk::paint::subPixelValue ) + { + outputImage->SetPixel( outputIndex, (TPixel2)( 0 ) ); + } + } + + // Set difference image + diffIterator.Set( static_cast(inputIterator.Get() - outputPixel ) ); // oh oh, not good for bigger values + ++inputIterator; + ++diffIterator; + } + } } std::string mitk::OverwriteSliceImageFilter::EventDescription( unsigned int sliceDimension, unsigned int sliceIndex, unsigned int timeStep ) @@ -244,3 +334,23 @@ return s.str(); } + +const mitk::Geometry3D* mitk::OverwriteSliceImageFilter::GetSliceGeometry3D() const +{ + return m_SliceGeometry3D; +} + +void mitk::OverwriteSliceImageFilter::SetSliceGeometry3D( const mitk::Geometry3D* val ) +{ + m_SliceGeometry3D = val; +} + +const mitk::Geometry3D* mitk::OverwriteSliceImageFilter::GetImageGeometry3D() const +{ + return m_ImageGeometry3D; +} + +void mitk::OverwriteSliceImageFilter::SetImageGeometry3D( const mitk::Geometry3D* val ) +{ + m_ImageGeometry3D = val; +} Index: Algorithms/mitkOverwriteSliceImageFilter.h =================================================================== --- Algorithms/mitkOverwriteSliceImageFilter.h (revision 20100) +++ Algorithms/mitkOverwriteSliceImageFilter.h (working copy) @@ -88,7 +88,14 @@ const Image* GetSliceImage() { return m_SliceImage.GetPointer(); } const Image* GetLastDifferenceImage() { return m_SliceDifferenceImage.GetPointer(); } + + const Geometry3D* GetSliceGeometry3D() const; + void SetSliceGeometry3D(const Geometry3D* val); + const Geometry3D* GetImageGeometry3D() const; + void SetImageGeometry3D(const Geometry3D* val); + + protected: OverwriteSliceImageFilter(); // purposely hidden @@ -114,6 +121,12 @@ unsigned int m_Dimension1; bool m_CreateUndoInformation; + + //! Slice and image are oriented on the same direction + bool m_bSliceAndImageOriented; + + const Geometry3D* m_ImageGeometry3D; + const Geometry3D* m_SliceGeometry3D; }; } // namespace Index: Interactions/mitkContourTool.cpp =================================================================== --- Interactions/mitkContourTool.cpp (revision 20100) +++ Interactions/mitkContourTool.cpp (working copy) @@ -23,6 +23,7 @@ #include "mitkBaseRenderer.h" #include "mitkRenderingManager.h" //#include "mitkProperties.h" +#include "mitkImageWriter.h" mitk::ContourTool::ContourTool(int paintingPixelValue) :FeedbackContourTool("PressMoveReleaseWithCTRLInversion"), @@ -102,7 +103,7 @@ DataTreeNode* workingNode( m_ToolManager->GetWorkingData(0) ); if (!workingNode) return false; - + Image* image = dynamic_cast(workingNode->GetData()); const PlaneGeometry* planeGeometry( dynamic_cast (positionEvent->GetSender()->GetCurrentWorldGeometry2D() ) ); if ( !image || !planeGeometry ) return false; @@ -113,51 +114,41 @@ { // 2. Slice is known, now we try to get it as a 2D image and project the contour into index coordinates of this slice Image::Pointer slice = SegTool2D::GetAffectedImageSliceAs2DImage( positionEvent, image ); - if ( slice.IsNull() ) { LOG_ERROR << "Unable to extract slice." << std::endl; return false; } + + // 3. Get the contour and project it on the 2D slice + Contour* feedbackContour( FeedbackContourTool::GetFeedbackContour() ); + // arg3: true: actually no idea why this is neccessary, but it works :-( + // arg4: false: don't constrain the contour to the image's inside + Contour::Pointer projectedContour = FeedbackContourTool::ProjectContourTo2DSlice( slice, feedbackContour, true, false ); + + if (projectedContour.IsNull()) + { + LOG_ERROR << "Unable to project contour." << std::endl; + return false; + } - /* - DataTreeNode::Pointer debugNode = DataTreeNode::New(); - debugNode->SetData( slice ); - debugNode->SetProperty( "name", StringProperty::New("extracted slice") ); - debugNode->SetProperty( "color", ColorProperty::New(1.0, 0.0, 0.0) ); - debugNode->SetProperty( "layer", FloatProperty::New(100) ); - m_ToolManager->GetDataStorage()->Add( debugNode ); - */ - - Contour* feedbackContour( FeedbackContourTool::GetFeedbackContour() ); - Contour::Pointer projectedContour = FeedbackContourTool::ProjectContourTo2DSlice( slice, feedbackContour, true, false ); // true: actually no idea why this is neccessary, but it works :-( - // false: don't constrain the contour to the image's inside - /* - Contour::Pointer back = FeedbackContourTool::BackProjectContourFrom2DSlice( slice, projectedContour, true ); // true: actually no idea why this is neccessary, but it works :-( - for (unsigned int idx = 0; idx < back->GetNumberOfPoints(); ++idx) - { - Point3D before = feedbackContour->GetPoints()->ElementAt(idx); - Point3D inbtw = projectedContour->GetPoints()->ElementAt(idx); - Point3D after = back->GetPoints()->ElementAt(idx); - - LOG_DEBUG << "before " << before << " zwischen " << inbtw << " after " << after; - } - */ - - if (projectedContour.IsNull()) return false; - - FeedbackContourTool::FillContourInSlice( projectedContour, slice, m_PaintingPixelValue ); - + // 4. Fill in the contour + FeedbackContourTool::FillContourInSlice( projectedContour, slice, m_PaintingPixelValue ); + // 5. Write the modified 2D working data slice back into the image + const Geometry3D* geometry3D = positionEvent->GetSender()->GetCurrentWorldGeometry(); + Geometry3D* imageGeometry = image->GetGeometry(0); OverwriteSliceImageFilter::Pointer slicewriter = OverwriteSliceImageFilter::New(); slicewriter->SetInput( image ); slicewriter->SetCreateUndoInformation( true ); slicewriter->SetSliceImage( slice ); slicewriter->SetSliceDimension( affectedDimension ); slicewriter->SetSliceIndex( affectedSlice ); + slicewriter->SetSliceGeometry3D( slice->GetGeometry() ); + slicewriter->SetImageGeometry3D( imageGeometry ); slicewriter->SetTimeStep( positionEvent->GetSender()->GetTimeStep( image ) ); slicewriter->Update(); - + // 6. Make sure the result is drawn again --> is visible then. assert( positionEvent->GetSender()->GetRenderWindow() ); @@ -179,14 +170,14 @@ if (!FeedbackContourTool::OnInvertLogic(action, stateEvent)) return false; // Inversion only for 0 and 1 as painting values - if (m_PaintingPixelValue == 1) + if (m_PaintingPixelValue == mitk::paint::addPixelValue) { - m_PaintingPixelValue = 0; + m_PaintingPixelValue = mitk::paint::subPixelValue; FeedbackContourTool::SetFeedbackContourColor( 1.0, 0.0, 0.0 ); } - else if (m_PaintingPixelValue == 0) + else if (m_PaintingPixelValue == mitk::paint::subPixelValue) { - m_PaintingPixelValue = 1; + m_PaintingPixelValue = mitk::paint::addPixelValue; FeedbackContourTool::SetFeedbackContourColorDefault(); } Index: Interactions/mitkContourTool.h =================================================================== --- Interactions/mitkContourTool.h (revision 20100) +++ Interactions/mitkContourTool.h (working copy) @@ -26,6 +26,14 @@ class Image; +// The intermediate data and the final one do not use the same pixel values +// in order not to confuse and re-process them. +namespace paint +{ + const int addPixelValue = 3; // unsigned + const int subPixelValue = 4; // unsigned +} + /** \brief Simple contour filling tool. Index: Interactions/mitkSegTool2D.cpp =================================================================== --- Interactions/mitkSegTool2D.cpp (revision 20100) +++ Interactions/mitkSegTool2D.cpp (working copy) @@ -24,6 +24,7 @@ #include "mitkPlaneGeometry.h" #include "mitkExtractImageFilter.h" +#include "mitkExtractDirectedPlaneImageFilter.h" #define ROUND(a) ((a)>0 ? (int)((a)+0.5) : -(int)(0.5-(a))) @@ -128,7 +129,7 @@ else { affectedDimension = -1; // no idea - return false; + return true; } // determine slice number in image @@ -165,26 +166,55 @@ int affectedSlice( -1 ); 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(); + if( affectedSlice != -1 ) + { + 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(); - - return slice; - } - catch(...) - { - // not working - return NULL; - } + // here we have a single slice that can be modified + // We need to copy the image, otherwise, when the image is used + // by another filter, the ExtractImageFilter::GenerateOutputInformation( ) + // will be called and the image will be Initialized to 0 + Image::Pointer slice = Image::New( ); + slice->Initialize( extractor->GetOutput() ); + + return slice; + } + catch(...) + { + // not working + return NULL; + } + } + else + { + //const Geometry3D* geometry3D = positionEvent->GetSender()->GetCurrentWorldGeometry(); + ExtractDirectedPlaneImageFilter::Pointer extractorPlane = ExtractDirectedPlaneImageFilter::New(); + extractorPlane->SetInput( image ); + extractorPlane->SetWorldGeometry( planeGeometry ); + extractorPlane->SetTargetTimestep( timeStep ); + extractorPlane->Update( ); + + // We need to copy the image, otherwise, when the image is used + // by another filter, the ExtractImageFilter::GenerateOutputInformation( ) + // will be called and the image will be Initialized to 0 + Image::Pointer imageReoriented = Image::New(); + imageReoriented->Initialize( extractorPlane->GetOutput( ) ); + + /*Point3D origin; + origin = image->GetGeometry( )->GetOrigin( ); + imageReoriented->GetGeometry( )->SetOrigin( origin );*/ + //imageReoriented->SetGeometry( (mitk::Geometry3D*) planeGeometry ); + + return imageReoriented; + } } else {