[Insight-users] problems using geodesic active contours segmentation
Luis Ibanez
luis.ibanez at kitware.com
Tue May 25 15:23:16 EDT 2010
Hi Paul,
Thanks for posting screenshots of your
intermediate images.
The don't quite have the expected appearance.
For example:
http://250kb.de/u85R1vR
should have had a constant value of 500 on the
outside, and a gradient ramp that goes darker
towards the original seed points.
...
Maybe I'm just misinterpreting your image.
Are you placing seed points inside every cell ?
If so, then you should look for the isovalue of
the FastMarching output for which the
spheres around the seed points have sizes
closer to the ones of the cells in the image.
This is the image that is passed as an initial
level set.
Once you have a good initial level set, you
should fine tune the speed image. Ideally this
speed image should look "almost" like a
segmentation of the objects (or at least of their
boundaries).
Regards,
Luis
------------------------------------------------------------
On Wed, May 12, 2010 at 9:16 AM, Paul <pare85 at googlemail.com> wrote:
> 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