Hi Paul,
Thanks for your detailed explanation and for
including the source code and the screenshots
of the images that you are working on.
The seed points for Fast Marching must be
placed inside of the object that you are
segmenting. (not in the border, but well inside).
Note that the Fast Marching filter is used here
just to create a set of initial circles (or spheres
if in 3D) that will then be used as input level
set for the GeodesicActive contour.
You should visually verify the output of the
Fast Marching filter, before you proceed further.
Save the output of the Fast marching filter to
a file, and then load it into any application
capable of visualizing images of pixel type
"float". (for example Paraview or Slicer).
Also, to prevent the Fast Marching filter from
running until the border of the image, please
make a call such as
fastMarching->SetStoppingValue( 500 );
(with this the fast marching front will not go
futher than 500 time steps from your seed
points).
The image should look like a distance map
from your seed points.
Please do this, and let us know what you
find.
Regards,
Luis
--------------------------------------
On Wed, May 5, 2010 at 6:36 AM, Paul <pare85 at googlemail.com> wrote:
> Hello mailinglist,
>
> I've got some problems using the Geodesic Active Contours Segmentation that
> is explained in the ITK Software Guide. I hope you can help me. In Figure
> 9.21 you can see the pipeline. Creating the Edge Image works fine. The
> problems start with creating the Input Levelset. First I'm not sure, were I
> have to set the seedpoints, in the object, on the contour, outside? I'm
> working with ct images of the spine (DICOM, 3D). Here is a part of the code.
> And a Link to the distance map und the segmentation result. I'm using MITK,
> a toolkit with interaction tools, like setting seed points. Region Growing
> und Fast Marching Segmentation works fine with it.
>
> IMG1: http://250kb.de/f8Wg0ug , shows edge image (black and white) and the
> distancemap ( red points)
>
> IMG2: http://250kb.de/NZ8lIUz , shows edge image (black and white) and the
> segmentation result ( red points)
>
> (here I set the seed points along the horizontal and the verticle axis of
> the shown image)
>
>
> CODE:
>
> typedef float InternalPixelType;
> const unsigned int Dimension = 3;
>
> typedef itk::Image< InternalPixelType, Dimension > InternalImageType;
>
> typedef unsigned char OutputPixelType;
> typedef itk::Image<OutputPixelType, VImageDimension> OutputImageType;
>
> typedef unsigned char CharPixelType;
> typedef short ShortPixelType;
>
> typedef itk::Image<CharPixelType, Dimension> CharImageType;
> typedef itk::Image<ShortPixelType, Dimension> ShortImageType;
>
> typedef itk::RescaleIntensityImageFilter<itk::Image<TPixel,
> VImageDimension>, InternalImageType > FilterType;
>
> typename FilterType::Pointer filter = FilterType::New();
> filter->SetOutputMinimum( -961.00 ); // MINIMUM PIXEL VALUE OF MY DATA
> SET
> filter->SetOutputMaximum( 3106.00 ); // MAXIMUM PIXEL VALUE OF MY DATA
> SET
>
> typedef itk::BinaryThresholdImageFilter<InternalImageType, OutputImageType>
> ThresholdingFilterType;
> typename ThresholdingFilterType:: Pointer
> thresholder=ThresholdingFilterType::New();
>
> thresholder-> SetLowerThreshold(getLower()); // I'M USING THE VALUES FROM
> THE GUIDE, DEFAULT: -1000 (FROM SHAPE DETECTION SEGMENTATION)
> thresholder-> SetUpperThreshold(getUpper()); // I'M USING THE VALUES FROM
> THE GUIDE, DEFAULT: 0 (FROM SHAPE DETECTION SEGMENTATION)
>
> thresholder-> SetOutsideValue(255);
> thresholder-> SetInsideValue(0);
>
> typedef itk::CurvatureAnisotropicDiffusionImageFilter<InternalImageType,
> InternalImageType > SmoothingFilterType;
> typename SmoothingFilterType::Pointer smoothing =
> SmoothingFilterType::New();
>
> typedef
> itk::GradientMagnitudeRecursiveGaussianImageFilter<InternalImageType,
> InternalImageType > GradientFilterType;
> typedef itk::SigmoidImageFilter<InternalImageType, InternalImageType >
> SigmoidFilterType;
>
> typename GradientFilterType::Pointer gradientMagnitude =
> GradientFilterType::New();
> typename SigmoidFilterType::Pointer sigmoid = SigmoidFilterType::New();
>
> sigmoid->SetOutputMinimum( 0.0 );
> sigmoid->SetOutputMaximum( 1.0 );
>
> typedef itk::FastMarchingImageFilter< InternalImageType, InternalImageType
>> FastMarchingFilterType;
> typename FastMarchingFilterType::Pointer fastMarching =
> FastMarchingFilterType::New();
>
> typedef itk::ShapeDetectionLevelSetImageFilter< InternalImageType,
> InternalImageType > ShapeDetectionFilterType;
> typename ShapeDetectionFilterType::Pointer shapeDetection =
> ShapeDetectionFilterType::New();
>
> typedef itk::GeodesicActiveContourLevelSetImageFilter< InternalImageType,
> InternalImageType > GeodesicActiveContourFilterType;
> typename GeodesicActiveContourFilterType::Pointer geodesicActiveContour =
> GeodesicActiveContourFilterType::New();
>
> filter->SetInput( itkImage ); // RESCALING TO FLOAT AS INPUT FOR THE
> CurvatureAnisotropicDiffusionImageFilter
> smoothing->SetInput( filter->GetOutput() );
> gradientMagnitude->SetInput( smoothing->GetOutput() );
> sigmoid->SetInput( gradientMagnitude->GetOutput() );
> fastMarching->SetInput( filter->GetOutput() ); // MITK DOES NOT ACCEPT
> IMAGES WITH A DIMENSION OF ZERO, FOR THIS REASON I HAVE TO SET AN INPUT FOR
> THE FAST MARCHING
> geodesicActiveContour->SetInput( fastMarching->GetOutput() );
> geodesicActiveContour->SetFeatureImage( sigmoid->GetOutput() );
> thresholder->SetInput( geodesicActiveContour->GetOutput() );
>
>
> smoothing->SetTimeStep( 0.0624 );
> smoothing->SetNumberOfIterations( 5 );
> smoothing->SetConductanceParameter( 9.0 );
>
> gradientMagnitude->SetSigma( getSigma() ); // I'M USING THE VALUES FROM THE
> GUIDE, DEFAULT: 1.0
>
> sigmoid->SetAlpha( getAlpha() ); // I'M USING THE VALUES FROM THE GUIDE,
> DEFAULT: -0.5
> sigmoid->SetBeta( getBeta() ); // I'M USING THE VALUES FROM THE GUIDE,
> DEFAULT: 2.0
>
>
>
> typedef typename FastMarchingFilterType::NodeContainer NodeContainer;
> typedef typename FastMarchingFilterType::NodeType NodeType;
> typename NodeContainer::Pointer seeds = NodeContainer::New();
>
> itk::Index<3> seedIndex;
>
> NodeType node;
> double seedValue= -5.0; // - initialDistance
>
> seeds->Initialize();
> int i=0;
>
> mitk::PointSet::PointsContainer* points =
> m_PointSet->GetPointSet()->GetPoints();
> for ( mitk::PointSet::PointsConstIterator pointsIterator = points->Begin();
> pointsIterator != points->End();
> ++pointsIterator )
> {
> if ( !imageGeometry->IsInside( pointsIterator.Value()) )
> {
> continue;
> }
>
> // convert world coordinates to image indices
> imageGeometry->WorldToIndex( pointsIterator.Value(), seedIndex);
>
> node.SetIndex( seedIndex );
> node.SetValue( seedValue );
>
> seeds->InsertElement( i, node );
> i+=1;
> }
>
>
> fastMarching->SetTrialPoints( seeds );
>
> fastMarching->SetOutputSize(
> filter->GetOutput()->GetBufferedRegion().GetSize() );
>
> fastMarching->SetSpeedConstant( 1.0 );
>
> geodesicActiveContour->SetPropagationScaling( 2.0 );
> geodesicActiveContour->SetCurvatureScaling( 1.0 );
> geodesicActiveContour->SetAdvectionScaling( 1.0 );
>
>
> try
> {
>
> thresholder->Update();
>
> // PUTTING THE RESULT IMAGE AS NODE IN THE DATATREE OF MITK
> mitk::Image::Pointer resultImage= mitk::ImportItkImage (thresholder->
> GetOutput());
> mitk::DataTreeNode::Pointer newNode= mitk::DataTreeNode::New();
> newNode->SetData( resultImage );
>
> newNode->SetProperty("binary", mitk::BoolProperty::New(true));
> newNode->SetProperty("name", mitk::StringProperty::New("dumb
> segmentation"));
> newNode->SetProperty("color", mitk::ColorProperty::New(1.0,0.0,0.0));
> newNode->SetProperty("volumerendering", mitk::BoolProperty::New(true));
> newNode->SetProperty("layer", mitk::IntProperty::New(1));
> newNode->SetProperty("opacity", mitk::FloatProperty::New(0.5));
>
> this->GetDefaultDataStorage()->Add(newNode);
>
> mitk::RenderingManager::GetInstance()->RequestUpdateAll();
> }
> catch( itk::ExceptionObject & excep )
> {
> std::cerr << "Exception caught !" << std::endl;
> std::cerr << excep << std::endl;
> }
>
>
> }
>
>
> best regards
> _____________________________________
> Powered by www.kitware.com
>
> Visit other Kitware open-source projects at
> http://www.kitware.com/opensource/opensource.html
>
> Kitware offers ITK Training Courses, for more information visit:
> http://www.kitware.com/products/protraining.html
>
> Please keep messages on-topic and check the ITK FAQ at:
> http://www.itk.org/Wiki/ITK_FAQ
>
> Follow this link to subscribe/unsubscribe:
> http://www.itk.org/mailman/listinfo/insight-users
>