[Insight-users] GeodesicActiveContour...Filter for 3D images again!
Oliver Gloger
olivergloger at web.de
Mon Jul 20 07:21:32 EDT 2009
Dear Luis, dear Insight-users,
I just repeat my reply from 25th May, because I could not solve my problem with your hints.
thank you very much for your fast answer,Luis! You recommended me to provide a "stoppingValue" and to test the rescalings ! Hence, I followed your advice, which did not solve my problems!
My problem is not the separate application of the ITK-FastMarchingImageFilter, because I want to apply the FastMarchingImageFilter as Input for the ITK-GeodesicActiveContourLevelSetImageFilter (see Fig. 9.21, page 583 )
Hence, I apply this filter as follows:
=> "(geodesicActiveContour->SetInput(fastMarching->getOutput()))"
In the ITK-examples code a "stoppingValue" is used for the FastMarchingImageFilter but not for the FastMarchingImageFilter as Input for the GeodesicActiveContourLevelSetImageFilter ! However, I provided a "stoppingValue" also for this application ("the FastMarchingImageFilter as Input for the GeodesicActiveContourLevelSetImageFilter"), but there is nothing to see! (I use the ITK-Snap-Tool for Visualization!) The ouput is "black"! My problem is that the final 3d-segmentation dataset (.gipl-format) is empty due to an empty fastMarching-Output! (With "empty" I mean, that the there is nothing to see, when I examine this output in the ITK-Snap tool ! By the way: how can I export the gipl-3D-dataset and examine the values for this output???)
Why is it possible that I get nice results for 2D but not for 3D when I use the example Code for the GeodesicActiveContourLevelSetImageFilter-Filter! I changed the example-code from 2d into 3d (Dimension=3) and provided GIPLImageIO-Filters for the readers and writers! Nothing more!!!! What could be wrong???
(I applied the example-code for the "itk::FastMarchingImageFilter.cxx" and this "separate" application of the FastMarchingImageFilter yields a valid segmentation result. Also for 3D! ). I am really wondering what could be wrong, but I do'nt habe any idea! Do you understand my problem better now?
my best regards,
Oliver
If anyone of the insight-users had a similar problem, you can compare it with my code (see below) !
----------------------------------------------------------------------------------------------------------------------------------------------
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
>
--
Dipl.-Ing., Dipl.-Inform.
Oliver Gloger
Technische Universität Berlin
Fachgebiet Computer Vision & Remote Sensing
Franklinstraße 28/29
Sekretariat FR 3-1
D-10587 Berlin
Tel.: 0049 30 314- 73109
______________________________________________________
GRATIS für alle WEB.DE-Nutzer: Die maxdome Movie-FLAT!
Jetzt freischalten unter http://movieflat.web.de
More information about the Insight-users
mailing list