[Insight-users] Problems with geodesic active contours

[Anders Almås]

Hi!
I have a problem in implementing geodesic active contours. My problem is
to segment aorta. Here is my file with typedef's and some other parameters
that might give you a feeling of what I am doing and what might be wrong:

// Number of dimensions in ImageType
  enum{ ImageDimension = 3 };
 // PixelType is the raw image pixel type
  typedef float PixelType;
  // ImageType is the raw image type
  typedef itk::Image<PixelType, ImageDimension> ImageType;
  // Derivative filter type, */
  typedef itk::GradientMagnitudeRecursiveGaussianImageFilter<ImageType>
DerivativeFilterType;
 //Rescale intensity for increase intnsityscale in images for viewing...
  typedef itk::RescaleIntensityImageFilter<ImageType, ImageType>
RescaleRangeType;
  // SigmoidFilterType for doing a sigmoid filtration
  typedef itk::SigmoidImageFilter<ImageType,ImageType> SigmoidFilterType;
  // FastMarchingFilterType for doing fastmarching filtration
  typedef itk::FastMarchingImageFilter<ImageType,ImageType>
                                                   FastMarchingFilterType;
  // GeodesicActiveContourFilterType for doing the geodesicActiveContour
filtration
  typedef itk::GeodesicActiveContourLevelSetImageFilter<ImageType,
                                ImageType,PixelType>
GeodesicActiveContourFilterType;
  // ThresholdType for doing thresholding
  typedef itk::ThresholdImageFilter<ImageType> ThresholdType;

std::cout<< "Starts preprocessing"<<std::endl;
  // Discrete CurvatureAnistropicDiffusionImageFilter

  typedef itk::CurvatureAnisotropicDiffusionImageFilter<ImageType,
ImageType> SmoothingFilterType;
  SmoothingFilterType::Pointer smoothing = SmoothingFilterType::New();
  smoothing->SetTimeStep(0.0625); // Default time step for 3D images with
itk...
  smoothing->SetIterations(1);    // Testable parameter...
  smoothing->SetConductanceParameter(1);
  smoothing->SetInput(input);
  smoothing->Update();

  RescaleRangeType::Pointer rangeSmoothing = RescaleRangeType::New();
  rangeSmoothing->SetOutputMinimum(min_intensity);
  rangeSmoothing->SetOutputMaximum(max_intensity);
  rangeSmoothing->SetInput(smoothing->GetOutput());
  rangeSmoothing->Update();

  //Setting up a gradientMagnitudeRecursiveGaussianImageFilter
  DerivativeFilterType::Pointer gradientMagnitude =
DerivativeFilterType::New();

  gradientMagnitude->SetSigma(1);
  gradientMagnitude->SetInput(smoothing->GetOutput());
  gradientMagnitude->Update();

  ThresholdType::Pointer thresholdGradient = ThresholdType::New();
  thresholdGradient->ThresholdBelow(50);
  thresholdGradient->SetOutsideValue(255);
  thresholdGradient->SetInput(rangeGradient->GetOutput());
  thresholdGradient->Update();

  //Setting up a sigmoid filter

  SigmoidFilterType::Pointer sigmoid = SigmoidFilterType::New();

  sigmoid->SetAlpha(-0.6);
  sigmoid->SetBeta(-0.05);
  sigmoid->SetInput(thresholdGradient->GetOutput());
  sigmoid->SetOutputMaximum(1.0);
  sigmoid->SetOutputMinimum(0.0);

  //Setting up a fast marching Filter

  FastMarchingFilterType::Pointer fastMarching =
FastMarchingFilterType::New();
  typedef FastMarchingFilterType::NodeType NodeType;
  typedef FastMarchingFilterType::NodeContainer NodeContainer;

  NodeType node;
  NodeType node2;
  NodeType node3;

  NodeContainer::Pointer trialpoints = NodeContainer::New();

  node.SetValue(0.0);
  node2.SetValue(0.0);
  node3.SetValue(0.0);

  ImageType::IndexType seed = {{lumenseed[0], lumenseed[1],
lumenseed[2]}};
  ImageType::IndexType seed2 = {{lumenseed2[0], lumenseed2[1],
lumenseed2[2]}};
  ImageType::IndexType seed3 = {{lumenseed3[0], lumenseed3[1],
lumenseed3[2]}};

  node.SetIndex(seed);
  node2.SetIndex(seed2);
  node3.SetIndex(seed3);

  trialpoints->InsertElement(0, node);
  trialpoints->InsertElement(1, node2);
  trialpoints->InsertElement(2, node3);

  fastMarching->SetTrialPoints(trialpoints);
  fastMarching->SetSpeedConstant(1.0);
  fastMarching->SetOutputSize(input->GetBufferedRegion().GetSize());
  //fastMarching->SetStoppingValue(10.0);
  fastMarching->Update();

   // Set GeodesicActiveContourImageFilter
  GeodesicActiveContourFilterType::Pointer geodesicActiveContour =
GeodesicActiveContourFilterType::New();

  geodesicActiveContour->SetPropagationScaling(5);
  geodesicActiveContour->SetCurvatureScaling(0.1);
  geodesicActiveContour->SetAdvectionScaling(1.0);
  geodesicActiveContour->SetMaximumRMSError( 0.02 );
  geodesicActiveContour->SetMaximumIterations( 100);
  geodesicActiveContour->SetInput(fastMarching->GetOutput());
  //geodesicActiveContour->SetFeatureImage(sigmoid->GetOutput());
  geodesicActiveContour->SetFeatureImage(thresholdGradient->GetOutput());
  geodesicActiveContour->Update();

  //Done GeodesicActiveContourImageFilter...

The rescaling intensity is used when I write the result to an image..
I get good results for my fast marching algorithm, but when I use and with
the sigmoid function.. but when I put these into GeodesicActiveContour I
get a result image with maximal pixel values (255 for me..). When I check
each slice I get 3 of 128 images with a result of an black region on a
white background(probably inverted result..)
What is then wrong here? Is it that my parameters are wrong or the code it
self?



-------------------------------------------------------------------------------
Navn:   Anders Almås
E-post: andersal at stud . ntnu . no
Tlf:    90550374
-------------------------------------------------------------------------------
























...............................................................................

Anders Almås            Min hjemmeside:

7230 Trondheim          http://stud . ntnu . no/~andersal/

...............................................................................