[Insight-users] problems using geodesic active contours segmentation
Paul
pare85 at googlemail.com
Wed May 12 09:16:20 EDT 2010
Hi Luis,
Thanks for your answer. Unfortunatly I did not find my mistake. I
verified the output of the Fast Marching filter and added the line
"fastMarching->SetStoppingValue( 500 );". You will see the output, if
you follow this links.
distance map (slicer): http://250kb.de/u85R1vR
distance map (mitk): http://250kb.de/L73Yfb5
distance map 3D (mitk): http://250kb.de/pqiHgT1
Here a link to the segmentation result using shown distance map:
http://250kb.de/aGDfpuA
Maybe you have got another tip for me.
regards, Paul
Am 09.05.2010 23:39, schrieb Luis Ibanez:
> 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
>>
>>
>
More information about the Insight-users
mailing list