Page MenuHomePhabricator
Files F255

contour_reoriented-svn20100.patch
Public
Actions

Authored By
• martelli
Feb 22 2010, 12:09 PM
Size
21 KB
Referenced Files
None
Subscribers
None

contour_reoriented-svn20100.patch
View Options

Index: Algorithms/mitkExtractDirectedPlaneImageFilter.cpp
===================================================================
--- Algorithms/mitkExtractDirectedPlaneImageFilter.cpp (revision 20100)
+++ Algorithms/mitkExtractDirectedPlaneImageFilter.cpp (working copy)
@@ -25,6 +25,7 @@
#include <vtkImageData.h>
#include <vtkImageChangeInformation.h>
#include <vtkPoints.h>
+#include <vtkImageResample.h>
#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<mitk::AffineTransform3D*>(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 <fstream>
+
#include "mitkOverwriteSliceImageFilter.h"
#include "mitkImageCast.h"
#include "mitkSegmentationInterpolationController.h"
@@ -25,6 +27,8 @@
#include "mitkDiffImageApplier.h"
#include "mitkImageTimeSelector.h"
+#include "mitkContourTool.h"
+
#include <itkImageSliceIteratorWithIndex.h>
#include <itkImageRegionIterator.h>
@@ -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<Image*>(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<short signed int>(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<short signed int>(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<TPixel2,VImageDimension2>::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<short signed int>(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<Image*>(workingNode->GetData());
const PlaneGeometry* planeGeometry( dynamic_cast<const PlaneGeometry*> (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
{

File Metadata

Mime Type
text/plain
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
246
Default Alt Text
contour_reoriented-svn20100.patch (21 KB)

Event Timeline

• martelli added a comment.Feb 22 2010, 12:09 PM

Patch for contour on reoriented planes to mitk svn 20100

MITK · MIC · IMSY · DKFZ