[Insight-users] FastMarching-output is empty using GeodesicActiveContourLevelSetImageFilter for 3D data sets

Luis Ibanez luis.ibanez at kitware.com
Tue May 19 21:05:49 EDT 2009 


Hi Oliver


The output of the FastMarching filter is a time crossing map,
that can be interpreted as a distance map to the seed points.


It seems that you forgot to set the parameter: StoppingValue.

Something like:


    fastMarching->SetStoppingValue(  stoppingTime  );


Please see the example:

      Insight/Examples/Segmentation/FastMarchingImageFilter.cxx

It is unlikely that your output image is empty.

Most likely it have a very high value that is throwing out
the rescaling of intensities.


----


About your theoretical question:

The Fast marching filter is combined with the GeodesicActive Contours
for *practical* reasons. It is simply used to provide an initialization
of the level set.  Any other quick segmentation method could have been
used to initialize the LevelSet as well.




     Regards,


        Luis



----------------------
Oliver Gloger wrote:
> Dear Insight-users,
> 
> 
> I want to apply the "itkGeodesicActiveContourLevelSetImageFilter" for 3D datasets. For 2D images I had success in applying this filter!
> However, for 3D datasets (I tried .gipl and .vtk-formats by using the corresponding IO-classes) I have problems using this filter!
> 
> I store the outputs from the "itkGradientMagnitudeImageFilter" and "itkSigmoidImageFilter" and they look fine! Just how I expected, they show the 3D edges (and inverse edges after applying the sigmoid filter) from the 3D data set.
> 
> However, the Output of the "itkFastMarchingImageFilter" is empty <=> contrarily to the 2D-case. In the 2D-case there is always a signed distance image!
> In 3D there is a gipl-OutputImage, but there is nothing to see!!! (I provide the fastMarchingfilter 1 seepoint!!!)
> 
> Consequently, the final binary result of the "itkGeodesicActiveContourLevelSetImageFilter" is also empty! Does the "itkGeodesicActiveContourLevelSetImageFilter" or the "itkFastMarchingImageFilter" not function for 3D datasets? Or could anybody of you give me a hint what I made wrong?I send the code below!
> 
> I have also another theoretical question: in SETHIAN's book the fast Marching approach is used for the boundary value problem by solving the Eikonal equation. However, this approach can only expand or shrink! Hence, the geodesic active contours could also move in both directions at certain locations due to curvature or edge attraction force! So, why is it justified to use the fast Marching approach in combination with the geodesic active contours?
> 
> Thank you very much for your help!
> 
> Olli
> 
> ---
> 
> #include "itkImage.h"
> #include "itkGeodesicActiveContourLevelSetImageFilter.h"
> #include "itkCurvatureAnisotropicDiffusionImageFilter.h"
> #include "itkGradientMagnitudeRecursiveGaussianImageFilter.h"
> #include "itkSigmoidImageFilter.h"
> #include "itkFastMarchingImageFilter.h"
> #include "itkRescaleIntensityImageFilter.h"
> #include "itkBinaryThresholdImageFilter.h"
> #include "itkImageFileReader.h"
> #include "itkImageFileWriter.h"
> #include "itkGDCMImageIO.h"
> #include "itkGiplImageIO.h"
> #include "itkVTKImageIO.h"
> #include "DicomSeriesReadImageWrite2.h"
> #include "itkGDCMImageIO.h"
> #include "itkGDCMSeriesFileNames.h"
> #include "itkImageSeriesReader.h"
> #include "itkImageFileWriter.h"
> #include "DicomSeriesReadImageWrite2.h"
> using namespace std;
> 
> //reading DicomSeries and storing volume data in gipl-file
> void DicomSeriesReadImageWrite2(char* DicomDirectory, char* outputFileName, char* seriesName)
> //This function is okay => I do not send it here!
> 
> int main( int argc, char *argv[] )
> {
> if( argc < 11 )
> {
> std::cerr << "Missing Parameters " << std::endl;
> std::cerr << "Usage: " << argv[0];
> std::cerr << " seedX seedY seedZ InitialDistance";
> std::cerr << " Sigma SigmoidAlpha SigmoidBeta";
> std::cerr << " PropagationScaling CurvatureScaling AdvectionScaling"<< std::endl;
> return 1;
> }
> 
> //char *path = "C:\\InsightToolkit\\Greifswald_Dicom_Datensaetze\\Proband_711286\\watersaturated3D\\";
> // String path = "C:\\InsightToolkit\\Greifswald_Dicom_Datensaetze\\Proband_711286\\watersaturated3D\\";
> 
> char *src_dicom = "C:\\InsightToolkit\\Greifswald_Dicom_Datensaetze\\Proband_711286\\watersaturated";
> char *dst_3D = "C:\\InsightToolkit\\Greifswald_Dicom_Datensaetze\\Proband_711286\\watersaturated3D\\dreiD.gipl";
> char *serName ="";
> 
> DicomSeriesReadImageWrite2(src_dicom,dst_3D,serName);
> 
> typedef float InternalPixelType;
> const unsigned int Dimension = 3; //here all should be in 3D !!!
> typedef itk::Image< InternalPixelType, Dimension > InternalImageType;
> 
> typedef unsigned char OutputPixelType;
> typedef itk::Image< OutputPixelType, Dimension > OutputImageType;
> 
> typedef itk::BinaryThresholdImageFilter< 
> InternalImageType, 
> OutputImageType > ThresholdingFilterType;
> 
> ThresholdingFilterType::Pointer thresholder = ThresholdingFilterType::New();
> 
> thresholder->SetLowerThreshold( -1000.0 );
> thresholder->SetUpperThreshold( 0.0 );
> 
> thresholder->SetOutsideValue( 0 );
> thresholder->SetInsideValue( 255 );
> 
> typedef itk::ImageFileReader< InternalImageType > ReaderType;
> typedef itk::ImageFileWriter< OutputImageType > WriterType;
> 
> ReaderType::Pointer reader = ReaderType::New();
> WriterType::Pointer writer = WriterType::New();
> 
> 
> typedef itk::GDCMImageIO ImageIOType;
> ImageIOType::Pointer gdcmImageIO = ImageIOType::New();
> typedef itk::VTKImageIO VTKIOType;
> VTKIOType::Pointer vtkIO = VTKIOType::New();
> typedef itk::GiplImageIO GiplIOType;
> GiplIOType::Pointer giplIO = GiplIOType::New(); 
> 
> reader->SetImageIO(giplIO);
> reader->SetFileName(dst_3D);
> 
> char *name = "C:\\InsightToolkit\\Greifswald_Dicom_Datensaetze\\Proband_711286\\watersaturated3D\\volume_dataset.gipl";
> writer->SetFileName(name);
> writer->SetImageIO(giplIO);
> typedef itk::RescaleIntensityImageFilter< 
> InternalImageType, 
> OutputImageType > CastFilterType;
> 
> typedef itk::CurvatureAnisotropicDiffusionImageFilter< 
> InternalImageType, 
> InternalImageType > SmoothingFilterType;
> 
> SmoothingFilterType::Pointer smoothing = SmoothingFilterType::New();
> 
> typedef itk::GradientMagnitudeRecursiveGaussianImageFilter< 
> InternalImageType, 
> InternalImageType > GradientFilterType;
> 
> typedef itk::SigmoidImageFilter< 
> InternalImageType, 
> InternalImageType > SigmoidFilterType;
> 
> GradientFilterType::Pointer gradientMagnitude = GradientFilterType::New();
> 
> SigmoidFilterType::Pointer sigmoid = SigmoidFilterType::New();
> 
> sigmoid->SetOutputMinimum( 0.0 );
> sigmoid->SetOutputMaximum( 1.0 );
> 
> typedef itk::FastMarchingImageFilter< 
> InternalImageType, 
> InternalImageType > FastMarchingFilterType;
> 
> 
> FastMarchingFilterType::Pointer fastMarching = FastMarchingFilterType::New();
> 
> typedef itk::GeodesicActiveContourLevelSetImageFilter< InternalImageType, 
> InternalImageType > GeodesicActiveContourFilterType;
> GeodesicActiveContourFilterType::Pointer geodesicActiveContour = 
> GeodesicActiveContourFilterType::New();
> 
> const double propagationScaling = atof( argv[8] );
> 
> geodesicActiveContour->SetPropagationScaling( propagationScaling );
> geodesicActiveContour->SetCurvatureScaling( atof(argv[9]) );
> geodesicActiveContour->SetAdvectionScaling( atof(argv[10]) );
> geodesicActiveContour->SetMaximumRMSError( 0.02 );
> geodesicActiveContour->SetNumberOfIterations( 5000 );
> 
> smoothing->SetInput( reader->GetOutput() );
> gradientMagnitude->SetInput( smoothing->GetOutput() );
> sigmoid->SetInput( gradientMagnitude->GetOutput() );
> 
> geodesicActiveContour->SetInput( fastMarching->GetOutput() );
> geodesicActiveContour->SetFeatureImage( sigmoid->GetOutput() );
> 
> thresholder->SetInput( geodesicActiveContour->GetOutput() );
> writer->SetInput( thresholder->GetOutput() );
> 
> smoothing->SetTimeStep( 0.0625 );
> smoothing->SetNumberOfIterations( 10 );
> smoothing->SetConductanceParameter( 1.0 );
> 
> const double sigma = atof( argv[5] );
> gradientMagnitude->SetSigma( sigma );
> 
> const double alpha = atof( argv[6] );
> const double beta = atof( argv[7] );
> 
> sigmoid->SetAlpha( alpha );
> sigmoid->SetBeta( beta );
> 
> typedef FastMarchingFilterType::NodeContainer NodeContainer;
> typedef FastMarchingFilterType::NodeType NodeType;
> 
> NodeContainer::Pointer seeds = NodeContainer::New();
> 
> InternalImageType::IndexType seedPosition;
> 
> seedPosition[0] = atoi( argv[1] );
> seedPosition[1] = atoi( argv[2] );
> seedPosition[2] = atoi( argv[3] );
> 
> const double initialDistance = atof( argv[4] );
> 
> NodeType node;
> 
> const double seedValue = - initialDistance;
> 
> node.SetValue( seedValue );
> node.SetIndex( seedPosition );
> 
> seeds->Initialize();
> seeds->InsertElement( 0, node );
> 
> fastMarching->SetTrialPoints( seeds );
> fastMarching->SetSpeedConstant( 1.0 );
> fastMarching->DebugOn();
> 
> CastFilterType::Pointer caster1 = CastFilterType::New();
> CastFilterType::Pointer caster2 = CastFilterType::New();
> CastFilterType::Pointer caster3 = CastFilterType::New();
> CastFilterType::Pointer caster4 = CastFilterType::New();
> 
> WriterType::Pointer writer1 = WriterType::New();
> WriterType::Pointer writer2 = WriterType::New();
> WriterType::Pointer writer3 = WriterType::New();
> WriterType::Pointer writer4 = WriterType::New();
> 
> caster1->SetInput( smoothing->GetOutput() );
> writer1->SetInput( caster1->GetOutput() );
> char *name1 = "C:\\InsightToolkit\\Greifswald_Dicom_Datensaetze\\Proband_711286\\watersaturated3D\\smoothed_dataset.gipl";
> writer1->SetFileName(name1);
> caster1->SetOutputMinimum( 0 );
> caster1->SetOutputMaximum( 255 );
> writer1->SetImageIO(giplIO);
> try{
> writer1->Update();
> }
> catch( itk::ExceptionObject & excep ){
> std::cerr << "Exception caught at writer1!" << std::endl;
> std::cerr << excep << std::endl;
> }
> 
> caster2->SetInput( gradientMagnitude->GetOutput() );
> writer2->SetInput( caster2->GetOutput() );
> char *name2 = "C:\\InsightToolkit\\Greifswald_Dicom_Datensaetze\\Proband_711286\\watersaturated3D\\3DKanten.gipl";
> writer2->SetFileName(name2);
> caster2->SetOutputMinimum( 0 );
> caster2->SetOutputMaximum( 255 );
> writer2->SetImageIO(giplIO);
> try{
> writer2->Update();
> }
> catch( itk::ExceptionObject & excep ){
> std::cerr << "Exception caught at writer2!" << std::endl;
> std::cerr << excep << std::endl;
> }
> 
> caster3->SetInput( sigmoid->GetOutput() );
> writer3->SetInput( caster3->GetOutput() );
> char *name3 = "C:\\InsightToolkit\\Greifswald_Dicom_Datensaetze\\Proband_711286\\watersaturated3D\\inverse3DKanten.gipl";
> writer3->SetFileName(name3);
> caster3->SetOutputMinimum( 0 );
> caster3->SetOutputMaximum( 255 );
> writer3->SetImageIO(giplIO);
> try{
> writer3->Update();
> }
> catch( itk::ExceptionObject & excep ){
> std::cerr << "Exception caught at writer3 !" << std::endl;
> std::cerr << excep << std::endl;
> }
> 
> caster4->SetInput( fastMarching->GetOutput() );
> writer4->SetInput( caster4->GetOutput() );
> char *name4 = "C:\\InsightToolkit\\Greifswald_Dicom_Datensaetze\\Proband_711286\\watersaturated3D\\FFM.gipl";
> writer4->SetFileName(name4);
> caster4->SetOutputMinimum( 0 );
> caster4->SetOutputMaximum( 255 );
> writer4->SetImageIO(giplIO);
> try{
> writer4->Update();
> std::cout<< "FMIF: GetOutputSize: "<<fastMarching->GetOutputSize()<<std::endl;
> }
> catch( itk::ExceptionObject & excep ){
> std::cerr << "Exception caught at writer4!" << std::endl;
> std::cerr << excep << std::endl;
> }
> 
> fastMarching->SetOutputSize( 
> reader->GetOutput()->GetBufferedRegion().GetSize() );
> try
> {
> writer->Update();
> }
> catch( itk::ExceptionObject & excep )
> {
> std::cerr << "Exception caught at main writer!" << std::endl;
> std::cerr << excep << std::endl;
> }
> std::cout << std::endl;
> std::cout << "Max. no. iterations: " << geodesicActiveContour->GetNumberOfIterations() << std::endl;
> std::cout << "Max. RMS error: " << geodesicActiveContour->GetMaximumRMSError() << std::endl;
> std::cout << std::endl;
> std::cout << "No. elpased iterations: " << geodesicActiveContour->GetElapsedIterations() << std::endl;
> std::cout << "RMS change: " << geodesicActiveContour->GetRMSChange() << std::endl;
> 
> return 0;
> }
> 
> 
> I call the filter with those parameters (for instance, but I also changed the values in a small range!!!)
> sigma=1, sigmiodAlpha = -4, SigmoidBeta=0, PrpagationScaling=3, CurvatureScaling=2, AdvectionScaling=3
> 
> 
> 
> ______________________________________________________
> GRATIS für alle WEB.DE-Nutzer: Die maxdome Movie-FLAT!
> Jetzt freischalten unter http://movieflat.web.de
> 
> _____________________________________
> Powered by www.kitware.com
> 
> Visit other Kitware open-source projects at
> http://www.kitware.com/opensource/opensource.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