<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
<html>
  <head>
    <meta content="text/html; charset=ISO-8859-1"
      http-equiv="Content-Type">
  </head>
  <body bgcolor="#ffffff" text="#3333ff">
    Hmmm.... I get it. Distance Map is just Euclidean Dist. Similar to
    bwdist of matlab. <br>
    <br>
    I just observed there is SignedDistanceMap as well... I think that
    should do the trick.<br>
    <pre class="moz-signature" cols="72">thanks
Vimal Singh
Graduate Student
ECE Dept.
Univ of Texas, Austin</pre>
    <br>
    On 2/1/2011 3:22 PM, Kishore Mosaliganti wrote:
    <blockquote
      cite="mid:AANLkTinrrJ94-t73y9riY4SjfDyKbmWFiJUbydeQN-R-@mail.gmail.com"
      type="cite">
      <pre wrap="">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 <a class="moz-txt-link-rfc2396E" href="mailto:vimal@mail.utexas.edu">&lt;vimal@mail.utexas.edu&gt;</a> wrote:
</pre>
      <blockquote type="cite">
        <pre wrap="">Hi all,

&nbsp;&nbsp;&nbsp; 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"&nbsp; //Include for distance map

int main( int argc, char *argv[] )
{
&nbsp; if( argc &lt; 10 )
&nbsp;&nbsp;&nbsp; {
&nbsp;&nbsp;&nbsp; std::cerr &lt;&lt; "Missing Parameters " &lt;&lt; std::endl;
&nbsp;&nbsp;&nbsp; std::cerr &lt;&lt; "Usage: " &lt;&lt; argv[0];
&nbsp;&nbsp;&nbsp; std::cerr &lt;&lt; " inputImage&nbsp; initialLevelSet outputImage";
&nbsp;&nbsp;&nbsp; std::cerr &lt;&lt; " Sigma SigmoidAlpha SigmoidBeta";
&nbsp;&nbsp;&nbsp; std::cerr &lt;&lt; " CurvatureScaling AdvectionTerm PropagationScaling"&nbsp; &lt;&lt;
std::endl;
&nbsp;&nbsp;&nbsp; return 1;
&nbsp;&nbsp;&nbsp; }

&nbsp; std::cout &lt;&lt; "Debug Point 1" &lt;&lt; std::endl;

&nbsp; typedef&nbsp;&nbsp; float&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; InternalPixelType;
&nbsp; const&nbsp;&nbsp;&nbsp;&nbsp; unsigned int&nbsp;&nbsp;&nbsp; Dimension = 2;
&nbsp; typedef itk::Image&lt; InternalPixelType, Dimension &gt;&nbsp; InternalImageType;
&nbsp; typedef unsigned char&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; OutputPixelType;
&nbsp; typedef itk::Image&lt; OutputPixelType, Dimension &gt; OutputImageType;
&nbsp; typedef itk::BinaryThresholdImageFilter&lt;
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; InternalImageType,
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; OutputImageType&nbsp;&nbsp;&nbsp; &gt;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; ThresholdingFilterType;

&nbsp; ThresholdingFilterType::Pointer thresholder =
ThresholdingFilterType::New();

&nbsp; thresholder-&gt;SetLowerThreshold( -1000.0 );
&nbsp; thresholder-&gt;SetUpperThreshold(&nbsp;&nbsp;&nbsp;&nbsp; 0.0 );

&nbsp; thresholder-&gt;SetOutsideValue(&nbsp; 0&nbsp; );
&nbsp; thresholder-&gt;SetInsideValue(&nbsp; 255 );

&nbsp; typedef itk::DanielssonDistanceMapImageFilter&lt;
&nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp; InternalImageType, InternalImageType &gt;
DistanceMapFilterType;
&nbsp; DistanceMapFilterType::Pointer distanceMap =
DistanceMapFilterType::New();&nbsp;&nbsp; //Instantiating the DistanceMap filter


&nbsp; std::cout &lt;&lt; "Debug Point 3" &lt;&lt; std::endl;
&nbsp; typedef&nbsp; itk::ImageFileReader&lt; InternalImageType &gt; ReaderType;
&nbsp; typedef&nbsp; itk::ImageFileWriter&lt;&nbsp; OutputImageType&nbsp; &gt; WriterType;

&nbsp; ReaderType::Pointer reader = ReaderType::New();&nbsp; //reads image
&nbsp; ReaderType::Pointer reader2 = ReaderType::New(); //reads thresholded
initial level-set&nbsp;&nbsp; //reading the inital binary level-set image
&nbsp; WriterType::Pointer writer = WriterType::New();

&nbsp; reader-&gt;SetFileName( argv[1] );
&nbsp; reader2-&gt;SetFileName( argv[2] );
&nbsp; writer-&gt;SetFileName( argv[3] );
&nbsp; typedef itk::RescaleIntensityImageFilter&lt;
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; InternalImageType,
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; OutputImageType &gt;&nbsp;&nbsp; CastFilterType;
&nbsp; typedef&nbsp;&nbsp; itk::CurvatureAnisotropicDiffusionImageFilter&lt;
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; InternalImageType,
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; InternalImageType &gt;&nbsp; SmoothingFilterType;

&nbsp; SmoothingFilterType::Pointer smoothing = SmoothingFilterType::New();
&nbsp; typedef&nbsp;&nbsp; itk::GradientMagnitudeRecursiveGaussianImageFilter&lt;
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; InternalImageType,
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; InternalImageType &gt;&nbsp; GradientFilterType;
&nbsp; typedef&nbsp;&nbsp; itk::SigmoidImageFilter&lt;
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; InternalImageType,
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; InternalImageType &gt;&nbsp; SigmoidFilterType;

&nbsp; GradientFilterType::Pointer&nbsp; gradientMagnitude =
GradientFilterType::New();

&nbsp; SigmoidFilterType::Pointer sigmoid = SigmoidFilterType::New();
&nbsp; sigmoid-&gt;SetOutputMinimum(&nbsp; 0.0&nbsp; );
&nbsp; sigmoid-&gt;SetOutputMaximum(&nbsp; 1.0&nbsp; );

&nbsp; typedef&nbsp; itk::GeodesicActiveContourLevelSetImageFilter&lt; InternalImageType,
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; InternalImageType &gt;&nbsp;&nbsp;&nbsp; GeodesicActiveContourFilterType;
&nbsp; GeodesicActiveContourFilterType::Pointer geodesicActiveContour =
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; GeodesicActiveContourFilterType::New();

&nbsp; const double propagationScaling = atof( argv[9] );
&nbsp; const double curvatureScaling = atof( argv[7] );
&nbsp; const double advectionScaling = atof( argv[8] );
&nbsp; geodesicActiveContour-&gt;SetPropagationScaling( propagationScaling );
&nbsp; geodesicActiveContour-&gt;SetCurvatureScaling( curvatureScaling );
&nbsp; geodesicActiveContour-&gt;SetAdvectionScaling( advectionScaling );
&nbsp; std::cout &lt;&lt; "Debug Point 7" &lt;&lt; std::endl;
&nbsp; geodesicActiveContour-&gt;SetMaximumRMSError( 0.02 );
&nbsp; geodesicActiveContour-&gt;SetNumberOfIterations(800);

&nbsp; smoothing-&gt;SetInput( reader-&gt;GetOutput() );
&nbsp; gradientMagnitude-&gt;SetInput( smoothing-&gt;GetOutput() );
&nbsp; sigmoid-&gt;SetInput( gradientMagnitude-&gt;GetOutput() );

&nbsp; distanceMap-&gt;SetInput( reader2-&gt;GetOutput() );
&nbsp; geodesicActiveContour-&gt;SetInput(&nbsp; distanceMap-&gt;GetOutput() );&nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp; //our
distance Mapped initial level-set
&nbsp; geodesicActiveContour-&gt;SetFeatureImage( sigmoid-&gt;GetOutput() );

&nbsp; thresholder-&gt;SetInput( geodesicActiveContour-&gt;GetOutput() );
&nbsp; writer-&gt;SetInput( thresholder-&gt;GetOutput() );

&nbsp; smoothing-&gt;SetTimeStep( 0.125 );
&nbsp; smoothing-&gt;SetNumberOfIterations(&nbsp; 5 );
&nbsp; smoothing-&gt;SetConductanceParameter( 9.0 );
&nbsp; std::cout &lt;&lt; "Debug Point 9" &lt;&lt; std::endl;
&nbsp; const double sigma = atof( argv[4] );
&nbsp; gradientMagnitude-&gt;SetSigma(&nbsp; sigma&nbsp; );
&nbsp; const double alpha =&nbsp; atof( argv[5] );
&nbsp; const double beta&nbsp; =&nbsp; atof( argv[6] );
&nbsp; sigmoid-&gt;SetAlpha( alpha );
&nbsp; sigmoid-&gt;SetBeta(&nbsp; beta&nbsp; );
&nbsp; std::cout &lt;&lt; "Debug Point 9.1" &lt;&lt; std::endl;

&nbsp; CastFilterType::Pointer caster1 = CastFilterType::New();
&nbsp; CastFilterType::Pointer caster2 = CastFilterType::New();
&nbsp; CastFilterType::Pointer caster3 = CastFilterType::New();
&nbsp; CastFilterType::Pointer caster4 = CastFilterType::New();

&nbsp; WriterType::Pointer writer1 = WriterType::New();
&nbsp; WriterType::Pointer writer2 = WriterType::New();
&nbsp; WriterType::Pointer writer3 = WriterType::New();
&nbsp; WriterType::Pointer writer4 = WriterType::New();
&nbsp; std::cout &lt;&lt; "Debug Point 9.15" &lt;&lt; std::endl;

&nbsp; caster1-&gt;SetInput( smoothing-&gt;GetOutput() );
&nbsp; writer1-&gt;SetInput( caster1-&gt;GetOutput() );
&nbsp; writer1-&gt;SetFileName("myGeodesicActiveContourImageFilterOutput1.png");
&nbsp; caster1-&gt;SetOutputMinimum(&nbsp;&nbsp; 0 );
&nbsp; caster1-&gt;SetOutputMaximum( 255 );
&nbsp; writer1-&gt;Update();
&nbsp; std::cout &lt;&lt; "Debug Point 9.2" &lt;&lt; std::endl;

&nbsp; caster2-&gt;SetInput( gradientMagnitude-&gt;GetOutput() );
&nbsp; writer2-&gt;SetInput( caster2-&gt;GetOutput() );
&nbsp; writer2-&gt;SetFileName("myGeodesicActiveContourImageFilterOutput2.png");
&nbsp; caster2-&gt;SetOutputMinimum(&nbsp;&nbsp; 0 );
&nbsp; caster2-&gt;SetOutputMaximum( 255 );
&nbsp; writer2-&gt;Update();
&nbsp; std::cout &lt;&lt; "Debug Point 9.3" &lt;&lt; std::endl;

&nbsp; caster3-&gt;SetInput( sigmoid-&gt;GetOutput() );
&nbsp; writer3-&gt;SetInput( caster3-&gt;GetOutput() );
&nbsp; writer3-&gt;SetFileName("myGeodesicActiveContourImageFilterOutput3.png");
&nbsp; caster3-&gt;SetOutputMinimum(&nbsp;&nbsp; 0 );
&nbsp; caster3-&gt;SetOutputMaximum( 255 );
&nbsp; writer3-&gt;Update();

&nbsp; //caster4-&gt;SetInput( distanceMap-&gt;GetOutput() );
&nbsp; caster4-&gt;SetInput( geodesicActiveContour-&gt;GetInput());
&nbsp; writer4-&gt;SetInput( caster4-&gt;GetOutput() );
&nbsp; writer4-&gt;SetFileName("myGeodesicActiveContourImageFilterOutput4.png");
&nbsp; caster4-&gt;SetOutputMinimum(&nbsp;&nbsp; 0 );
&nbsp; caster4-&gt;SetOutputMaximum( 255 );
&nbsp; //writer4-&gt;Update();

&nbsp; std::cout &lt;&lt; "Debug Point 10" &lt;&lt; std::endl;

&nbsp; try
&nbsp;&nbsp;&nbsp; {
&nbsp;&nbsp;&nbsp; writer-&gt;Update();
&nbsp;&nbsp;&nbsp; }
&nbsp; catch( itk::ExceptionObject &amp; excep )
&nbsp;&nbsp;&nbsp; {
&nbsp;&nbsp;&nbsp; std::cerr &lt;&lt; "Exception caught !" &lt;&lt; std::endl;
&nbsp;&nbsp;&nbsp; std::cerr &lt;&lt; excep &lt;&lt; std::endl;
&nbsp;&nbsp;&nbsp; }


&nbsp; // Print out some useful information
&nbsp; std::cout &lt;&lt; std::endl;
&nbsp; std::cout &lt;&lt; "Max. no. iterations: " &lt;&lt;
geodesicActiveContour-&gt;GetNumberOfIterations() &lt;&lt; std::endl;
&nbsp; std::cout &lt;&lt; "Max. RMS error: " &lt;&lt;
geodesicActiveContour-&gt;GetMaximumRMSError() &lt;&lt; std::endl;
&nbsp; std::cout &lt;&lt; std::endl;
&nbsp; std::cout &lt;&lt; "No. elpased iterations: " &lt;&lt;
geodesicActiveContour-&gt;GetElapsedIterations() &lt;&lt; std::endl;
&nbsp; std::cout &lt;&lt; "RMS change: " &lt;&lt; geodesicActiveContour-&gt;GetRMSChange() &lt;&lt;
std::endl;

&nbsp; writer4-&gt;Update();
&nbsp; return 0;
}
_____________________________________
Powered by <a class="moz-txt-link-abbreviated" href="http://www.kitware.com">www.kitware.com</a>

Visit other Kitware open-source projects at
<a class="moz-txt-link-freetext" href="http://www.kitware.com/opensource/opensource.html">http://www.kitware.com/opensource/opensource.html</a>

Kitware offers ITK Training Courses, for more information visit:
<a class="moz-txt-link-freetext" href="http://www.kitware.com/products/protraining.html">http://www.kitware.com/products/protraining.html</a>

Please keep messages on-topic and check the ITK FAQ at:
<a class="moz-txt-link-freetext" href="http://www.itk.org/Wiki/ITK_FAQ">http://www.itk.org/Wiki/ITK_FAQ</a>

Follow this link to subscribe/unsubscribe:
<a class="moz-txt-link-freetext" href="http://www.itk.org/mailman/listinfo/insight-users">http://www.itk.org/mailman/listinfo/insight-users</a>


</pre>
      </blockquote>
    </blockquote>
  </body>
</html>