<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
<html>
<head>
<meta http-equiv="content-type" content="text/html; charset=ISO-8859-1">
</head>
<body bgcolor="#ffffff" text="#3333ff">
Hi all,<br>
<br>
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. <br>
<br>
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.<br>
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. <br>
<br>
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 ?<br>
<br>
<br>
#include "itkImage.h"<br>
#include "itkGeodesicActiveContourLevelSetImageFilter.h"<br>
#include "itkCurvatureAnisotropicDiffusionImageFilter.h"<br>
#include "itkGradientMagnitudeRecursiveGaussianImageFilter.h"<br>
#include "itkSigmoidImageFilter.h"<br>
//#include "itkFastMarchingImageFilter.h"<br>
#include "itkRescaleIntensityImageFilter.h"<br>
#include "itkBinaryThresholdImageFilter.h"<br>
#include "itkImageFileReader.h"<br>
#include "itkImageFileWriter.h"<br>
<br>
#include "itkDanielssonDistanceMapImageFilter.h" //Include for
distance map<br>
<br>
int main( int argc, char *argv[] )<br>
{<br>
if( argc < 10 )<br>
{<br>
std::cerr << "Missing Parameters " << std::endl;<br>
std::cerr << "Usage: " << argv[0];<br>
std::cerr << " inputImage initialLevelSet outputImage";<br>
std::cerr << " Sigma SigmoidAlpha SigmoidBeta";<br>
std::cerr << " CurvatureScaling AdvectionTerm
PropagationScaling" << std::endl;<br>
return 1;<br>
}<br>
<br>
std::cout << "Debug Point 1" << std::endl;<br>
<br>
typedef float InternalPixelType;<br>
const unsigned int Dimension = 2;<br>
typedef itk::Image< InternalPixelType, Dimension >
InternalImageType;<br>
typedef unsigned char OutputPixelType;<br>
typedef itk::Image< OutputPixelType, Dimension >
OutputImageType;<br>
typedef itk::BinaryThresholdImageFilter< <br>
InternalImageType, <br>
OutputImageType >
ThresholdingFilterType;<br>
<br>
ThresholdingFilterType::Pointer thresholder =
ThresholdingFilterType::New();<br>
<br>
thresholder->SetLowerThreshold( -1000.0 );<br>
thresholder->SetUpperThreshold( 0.0 );<br>
<br>
thresholder->SetOutsideValue( 0 );<br>
thresholder->SetInsideValue( 255 );<br>
<br>
typedef itk::DanielssonDistanceMapImageFilter<<br>
InternalImageType, InternalImageType >
DistanceMapFilterType;<br>
DistanceMapFilterType::Pointer distanceMap =
DistanceMapFilterType::New(); //Instantiating the DistanceMap
filter<br>
<br>
<br>
std::cout << "Debug Point 3" << std::endl; <br>
typedef itk::ImageFileReader< InternalImageType >
ReaderType;<br>
typedef itk::ImageFileWriter< OutputImageType >
WriterType;<br>
<br>
ReaderType::Pointer reader = ReaderType::New(); //reads image<br>
ReaderType::Pointer reader2 = ReaderType::New(); //reads
thresholded initial level-set //reading the inital binary
level-set image<br>
WriterType::Pointer writer = WriterType::New(); <br>
<br>
reader->SetFileName( argv[1] );<br>
reader2->SetFileName( argv[2] );<br>
writer->SetFileName( argv[3] );<br>
typedef itk::RescaleIntensityImageFilter< <br>
InternalImageType, <br>
OutputImageType >
CastFilterType;<br>
typedef itk::CurvatureAnisotropicDiffusionImageFilter< <br>
InternalImageType, <br>
InternalImageType >
SmoothingFilterType;<br>
<br>
SmoothingFilterType::Pointer smoothing =
SmoothingFilterType::New();<br>
typedef itk::GradientMagnitudeRecursiveGaussianImageFilter< <br>
InternalImageType, <br>
InternalImageType >
GradientFilterType;<br>
typedef itk::SigmoidImageFilter<<br>
InternalImageType, <br>
InternalImageType >
SigmoidFilterType;<br>
<br>
GradientFilterType::Pointer gradientMagnitude =
GradientFilterType::New();<br>
<br>
SigmoidFilterType::Pointer sigmoid = SigmoidFilterType::New();<br>
sigmoid->SetOutputMinimum( 0.0 );<br>
sigmoid->SetOutputMaximum( 1.0 );<br>
<br>
typedef itk::GeodesicActiveContourLevelSetImageFilter<
InternalImageType, <br>
InternalImageType >
GeodesicActiveContourFilterType;<br>
GeodesicActiveContourFilterType::Pointer geodesicActiveContour = <br>
GeodesicActiveContourFilterType::New(); <br>
<br>
const double propagationScaling = atof( argv[9] );<br>
const double curvatureScaling = atof( argv[7] );<br>
const double advectionScaling = atof( argv[8] );<br>
geodesicActiveContour->SetPropagationScaling(
propagationScaling );<br>
geodesicActiveContour->SetCurvatureScaling( curvatureScaling );<br>
geodesicActiveContour->SetAdvectionScaling( advectionScaling );<br>
std::cout << "Debug Point 7" << std::endl;<br>
geodesicActiveContour->SetMaximumRMSError( 0.02 );<br>
geodesicActiveContour->SetNumberOfIterations(800);<br>
<br>
smoothing->SetInput( reader->GetOutput() );<br>
gradientMagnitude->SetInput( smoothing->GetOutput() );<br>
sigmoid->SetInput( gradientMagnitude->GetOutput() );<br>
<br>
distanceMap->SetInput( reader2->GetOutput() );<br>
geodesicActiveContour->SetInput( distanceMap->GetOutput()
); //our distance Mapped initial level-set<br>
geodesicActiveContour->SetFeatureImage( sigmoid->GetOutput()
);<br>
<br>
thresholder->SetInput( geodesicActiveContour->GetOutput() );<br>
writer->SetInput( thresholder->GetOutput() );<br>
<br>
smoothing->SetTimeStep( 0.125 );<br>
smoothing->SetNumberOfIterations( 5 );<br>
smoothing->SetConductanceParameter( 9.0 );<br>
std::cout << "Debug Point 9" << std::endl;<br>
const double sigma = atof( argv[4] );<br>
gradientMagnitude->SetSigma( sigma );<br>
const double alpha = atof( argv[5] );<br>
const double beta = atof( argv[6] );<br>
sigmoid->SetAlpha( alpha );<br>
sigmoid->SetBeta( beta );<br>
std::cout << "Debug Point 9.1" << std::endl;<br>
<br>
CastFilterType::Pointer caster1 = CastFilterType::New();<br>
CastFilterType::Pointer caster2 = CastFilterType::New();<br>
CastFilterType::Pointer caster3 = CastFilterType::New();<br>
CastFilterType::Pointer caster4 = CastFilterType::New();<br>
<br>
WriterType::Pointer writer1 = WriterType::New();<br>
WriterType::Pointer writer2 = WriterType::New();<br>
WriterType::Pointer writer3 = WriterType::New();<br>
WriterType::Pointer writer4 = WriterType::New();<br>
std::cout << "Debug Point 9.15" << std::endl;<br>
<br>
caster1->SetInput( smoothing->GetOutput() );<br>
writer1->SetInput( caster1->GetOutput() );<br>
writer1->SetFileName("myGeodesicActiveContourImageFilterOutput1.png");<br>
caster1->SetOutputMinimum( 0 );<br>
caster1->SetOutputMaximum( 255 );<br>
writer1->Update();<br>
std::cout << "Debug Point 9.2" << std::endl;<br>
<br>
caster2->SetInput( gradientMagnitude->GetOutput() );<br>
writer2->SetInput( caster2->GetOutput() );<br>
writer2->SetFileName("myGeodesicActiveContourImageFilterOutput2.png");<br>
caster2->SetOutputMinimum( 0 );<br>
caster2->SetOutputMaximum( 255 );<br>
writer2->Update();<br>
std::cout << "Debug Point 9.3" << std::endl;<br>
<br>
caster3->SetInput( sigmoid->GetOutput() );<br>
writer3->SetInput( caster3->GetOutput() );<br>
writer3->SetFileName("myGeodesicActiveContourImageFilterOutput3.png");<br>
caster3->SetOutputMinimum( 0 );<br>
caster3->SetOutputMaximum( 255 );<br>
writer3->Update();<br>
<br>
//caster4->SetInput( distanceMap->GetOutput() );<br>
caster4->SetInput( geodesicActiveContour->GetInput());<br>
writer4->SetInput( caster4->GetOutput() );<br>
writer4->SetFileName("myGeodesicActiveContourImageFilterOutput4.png");<br>
caster4->SetOutputMinimum( 0 );<br>
caster4->SetOutputMaximum( 255 );<br>
//writer4->Update();<br>
<br>
std::cout << "Debug Point 10" << std::endl;<br>
<br>
try<br>
{<br>
writer->Update();<br>
}<br>
catch( itk::ExceptionObject & excep )<br>
{<br>
std::cerr << "Exception caught !" << std::endl;<br>
std::cerr << excep << std::endl;<br>
}<br>
<br>
<br>
// Print out some useful information <br>
std::cout << std::endl;<br>
std::cout << "Max. no. iterations: " <<
geodesicActiveContour->GetNumberOfIterations() <<
std::endl;<br>
std::cout << "Max. RMS error: " <<
geodesicActiveContour->GetMaximumRMSError() << std::endl;<br>
std::cout << std::endl;<br>
std::cout << "No. elpased iterations: " <<
geodesicActiveContour->GetElapsedIterations() << std::endl;<br>
std::cout << "RMS change: " <<
geodesicActiveContour->GetRMSChange() << std::endl;<br>
<br>
writer4->Update();<br>
return 0;<br>
}
</body>
</html>