[Insight-users] Help with GeodesicActiveContourImageFilter --> Modification to example code does not seem to work
Vimal Singh
vimal at mail.utexas.edu
Tue Feb 1 16:41:46 EST 2011
Hmmm.... I get it. Distance Map is just Euclidean Dist. Similar to
bwdist of matlab.
I just observed there is SignedDistanceMap as well... I think that
should do the trick.
thanks
Vimal Singh
Graduate Student
ECE Dept.
Univ of Texas, Austin
On 2/1/2011 3:22 PM, Kishore Mosaliganti wrote:
> Hi Vimal,
>
> The levelset function is supposed to be initialized in such a way that
> there is a negative and positive part in the range of the function.
> The fast marching filter allows one to specify negative seeds and
> therefore it generates a 0 levelset which is a circle.
>
> The output of the Danielsson Distance map is a fully positive function.
>
> To make it work, subtract some value from its output before setting it
> as input to the levelset filter.
>
> Kishore
>
> On Tue, Feb 1, 2011 at 4:15 PM, Vimal Singh<vimal at mail.utexas.edu> wrote:
>> Hi all,
>>
>> I'm a new User to ITK and have been learning it for past few days. I
>> managed to run the example codes and observe the same results available in
>> software guide document.
>>
>> Currently I'm trying to develop my GeodesicActiveContourImageFilter, which
>> is similar to the one provided in the examples/segmentation folder but does
>> not use FastMarchingImageFilter for generation of initial Level Set. Instead
>> relies on user input for this initial level-set.
>> I've modified the example code accordingly and observed that the DistanceMap
>> obtained using DanielssonDistanceMapImageFilter matches exactly the
>> distance-map of initial level set generated by Example code. However my
>> final segmentation result is not the same, despite the parameters for all
>> other filters in two codes have been kept the same. My complete code is
>> below.
>>
>> Can anyone give me hints as to what I might be missing. I've removed the
>> FastMarchingImageFilter all together. Any kind off help would be highly
>> appreciated !!! Or cite me an easier way to debug things ?
>>
>>
>> #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 "itkDanielssonDistanceMapImageFilter.h" //Include for distance map
>>
>> int main( int argc, char *argv[] )
>> {
>> if( argc< 10 )
>> {
>> std::cerr<< "Missing Parameters "<< std::endl;
>> std::cerr<< "Usage: "<< argv[0];
>> std::cerr<< " inputImage initialLevelSet outputImage";
>> std::cerr<< " Sigma SigmoidAlpha SigmoidBeta";
>> std::cerr<< " CurvatureScaling AdvectionTerm PropagationScaling"<<
>> std::endl;
>> return 1;
>> }
>>
>> std::cout<< "Debug Point 1"<< std::endl;
>>
>> typedef float InternalPixelType;
>> const unsigned int Dimension = 2;
>> 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::DanielssonDistanceMapImageFilter<
>> InternalImageType, InternalImageType>
>> DistanceMapFilterType;
>> DistanceMapFilterType::Pointer distanceMap =
>> DistanceMapFilterType::New(); //Instantiating the DistanceMap filter
>>
>>
>> std::cout<< "Debug Point 3"<< std::endl;
>> typedef itk::ImageFileReader< InternalImageType> ReaderType;
>> typedef itk::ImageFileWriter< OutputImageType> WriterType;
>>
>> ReaderType::Pointer reader = ReaderType::New(); //reads image
>> ReaderType::Pointer reader2 = ReaderType::New(); //reads thresholded
>> initial level-set //reading the inital binary level-set image
>> WriterType::Pointer writer = WriterType::New();
>>
>> reader->SetFileName( argv[1] );
>> reader2->SetFileName( argv[2] );
>> writer->SetFileName( argv[3] );
>> 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::GeodesicActiveContourLevelSetImageFilter< InternalImageType,
>> InternalImageType> GeodesicActiveContourFilterType;
>> GeodesicActiveContourFilterType::Pointer geodesicActiveContour =
>> GeodesicActiveContourFilterType::New();
>>
>> const double propagationScaling = atof( argv[9] );
>> const double curvatureScaling = atof( argv[7] );
>> const double advectionScaling = atof( argv[8] );
>> geodesicActiveContour->SetPropagationScaling( propagationScaling );
>> geodesicActiveContour->SetCurvatureScaling( curvatureScaling );
>> geodesicActiveContour->SetAdvectionScaling( advectionScaling );
>> std::cout<< "Debug Point 7"<< std::endl;
>> geodesicActiveContour->SetMaximumRMSError( 0.02 );
>> geodesicActiveContour->SetNumberOfIterations(800);
>>
>> smoothing->SetInput( reader->GetOutput() );
>> gradientMagnitude->SetInput( smoothing->GetOutput() );
>> sigmoid->SetInput( gradientMagnitude->GetOutput() );
>>
>> distanceMap->SetInput( reader2->GetOutput() );
>> geodesicActiveContour->SetInput( distanceMap->GetOutput() ); //our
>> distance Mapped initial level-set
>> geodesicActiveContour->SetFeatureImage( sigmoid->GetOutput() );
>>
>> thresholder->SetInput( geodesicActiveContour->GetOutput() );
>> writer->SetInput( thresholder->GetOutput() );
>>
>> smoothing->SetTimeStep( 0.125 );
>> smoothing->SetNumberOfIterations( 5 );
>> smoothing->SetConductanceParameter( 9.0 );
>> std::cout<< "Debug Point 9"<< std::endl;
>> const double sigma = atof( argv[4] );
>> gradientMagnitude->SetSigma( sigma );
>> const double alpha = atof( argv[5] );
>> const double beta = atof( argv[6] );
>> sigmoid->SetAlpha( alpha );
>> sigmoid->SetBeta( beta );
>> std::cout<< "Debug Point 9.1"<< std::endl;
>>
>> 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();
>> std::cout<< "Debug Point 9.15"<< std::endl;
>>
>> caster1->SetInput( smoothing->GetOutput() );
>> writer1->SetInput( caster1->GetOutput() );
>> writer1->SetFileName("myGeodesicActiveContourImageFilterOutput1.png");
>> caster1->SetOutputMinimum( 0 );
>> caster1->SetOutputMaximum( 255 );
>> writer1->Update();
>> std::cout<< "Debug Point 9.2"<< std::endl;
>>
>> caster2->SetInput( gradientMagnitude->GetOutput() );
>> writer2->SetInput( caster2->GetOutput() );
>> writer2->SetFileName("myGeodesicActiveContourImageFilterOutput2.png");
>> caster2->SetOutputMinimum( 0 );
>> caster2->SetOutputMaximum( 255 );
>> writer2->Update();
>> std::cout<< "Debug Point 9.3"<< std::endl;
>>
>> caster3->SetInput( sigmoid->GetOutput() );
>> writer3->SetInput( caster3->GetOutput() );
>> writer3->SetFileName("myGeodesicActiveContourImageFilterOutput3.png");
>> caster3->SetOutputMinimum( 0 );
>> caster3->SetOutputMaximum( 255 );
>> writer3->Update();
>>
>> //caster4->SetInput( distanceMap->GetOutput() );
>> caster4->SetInput( geodesicActiveContour->GetInput());
>> writer4->SetInput( caster4->GetOutput() );
>> writer4->SetFileName("myGeodesicActiveContourImageFilterOutput4.png");
>> caster4->SetOutputMinimum( 0 );
>> caster4->SetOutputMaximum( 255 );
>> //writer4->Update();
>>
>> std::cout<< "Debug Point 10"<< std::endl;
>>
>> try
>> {
>> writer->Update();
>> }
>> catch( itk::ExceptionObject& excep )
>> {
>> std::cerr<< "Exception caught !"<< std::endl;
>> std::cerr<< excep<< std::endl;
>> }
>>
>>
>> // Print out some useful information
>> 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;
>>
>> writer4->Update();
>> return 0;
>> }
>> _____________________________________
>> 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
>>
>>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.itk.org/pipermail/insight-users/attachments/20110201/ae956e28/attachment.htm>
More information about the Insight-users
mailing list