[Insight-users] Active Shape Models Segmentation

[Ilaria Canova Calori]

Hi,

I am working on segmentation of AAA in CTimages.
I found really intresting the file report "Active Shape Models_itk.ppt"
in the InsightDocuments checkout. I tried to modified the
itkImagePCAShapeModelEstimatorTest.cxx to see how it handles some
training dataset but it always returns the last training image and the
eigen values are all zero.

I attach the code and the training images, any suggestion?
Thanks in advance,

Ilaria Canova Calori

--
-------------- next part --------------

// Insight classes
#include "itkImage.h"
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include "itkVector.h"
#include "vnl/vnl_matrix_fixed.h"
#include "vnl/vnl_math.h"
#include "itkImageRegionIterator.h"
#include "itkLightProcessObject.h"
#include "itkTextOutput.h"

#include "itkImagePCAShapeModelEstimator.h"

int main(int argc, char *argv[] )
{

  if( argc < 5 )
    {
    std::cerr << "Missing Parameters " << std::endl;
    std::cerr << "Usage: " << argv[0];
    std::cerr << " firstShapeModel secondShapeModel ThirdShapeModel ";
	std::cerr << " meanOutputFile principalComponentFile" <<std::endl;
    return 1;
    }

  typedef unsigned char PixelType;
  const unsigned int  Dimension = 2;
  const unsigned int NUMTRAINIMAGES = 3;
  const unsigned int NUMLARGESTPC = 1;
  typedef itk::Image<PixelType,Dimension> InputImageType; 
  typedef itk::Image<PixelType,Dimension> OutputImageType; 
  typedef itk::Image<PixelType,Dimension> MeanImageType;

  typedef itk::ImageFileReader< InputImageType >  ReaderType;
  typedef itk::ImageFileWriter< OutputImageType >  WriterType;

  typedef itk::ImageRegionIterator< InputImageType > InputImageIterator;
  typedef itk::ImageRegionIterator< OutputImageType > OutputImageIterator;

  ReaderType::Pointer reader = ReaderType::New();
  WriterType::Pointer writer = WriterType::New();

  reader->SetFileName( argv[1] );
  reader->Update();
  
  InputImageType::Pointer image1 = InputImageType::New();

  InputImageType::Pointer image2 = InputImageType::New();

  InputImageType::Pointer image3 = InputImageType::New();

  InputImageType::SizeType inputImageSize ;
  inputImageSize = reader->GetOutput()->GetRequestedRegion().GetSize();

  InputImageType::IndexType index;
  index.Fill(0);

  InputImageType::RegionType region1;
  region1.SetSize( inputImageSize );
  region1.SetIndex( index );

  //--------------------------------------------------------------------------
  // Set up Image 1 first
  //--------------------------------------------------------------------------

  image1->SetLargestPossibleRegion( region1 );
  image1->SetBufferedRegion( region1 );
  image1->Allocate();
  image1 = reader->GetOutput();

  // setup the iterators
  InputImageIterator image1It( image1, image1->GetBufferedRegion() );

  //--------------------------------------------------------------------------
  // Set up Image 2 first
  //--------------------------------------------------------------------------

  reader->SetFileName( argv[2] );
  reader->Update();
  InputImageType::RegionType region2;
  inputImageSize = reader->GetOutput()->GetRequestedRegion().GetSize();
  region2.SetSize( inputImageSize );
  region2.SetIndex( index );

  image2->SetLargestPossibleRegion( region2 );
  image2->SetBufferedRegion( region2 );
  image2->Allocate();
  image2 = reader->GetOutput();

  // setup the iterators
  InputImageIterator image2It( image2, image2->GetBufferedRegion() );

  //--------------------------------------------------------------------------
  // Set up Image 3 first
  //--------------------------------------------------------------------------

  reader->SetFileName( argv[3] );
  reader->Update();
  InputImageType::RegionType region3;
  inputImageSize = reader->GetOutput()->GetRequestedRegion().GetSize();
  region3.SetSize( inputImageSize );
  region3.SetIndex( index );

  image3->SetLargestPossibleRegion( region3 );
  image3->SetBufferedRegion( region3 );
  image3->Allocate();
  image3 = reader->GetOutput();

  // setup the iterators
  InputImageIterator image3It( image3, image3->GetBufferedRegion() );

  //----------------------------------------------------------------------
  // Test code for the Shape model estimator
  //----------------------------------------------------------------------

  //----------------------------------------------------------------------
  //Set the image model estimator
  //----------------------------------------------------------------------
  typedef itk::ImagePCAShapeModelEstimator<InputImageType, OutputImageType> 
    ImagePCAShapeModelEstimatorType;

  ImagePCAShapeModelEstimatorType::Pointer 
    applyPCAShapeEstimator = ImagePCAShapeModelEstimatorType::New();

  //----------------------------------------------------------------------
  //Set the parameters of the clusterer
  //----------------------------------------------------------------------
  applyPCAShapeEstimator->SetNumberOfTrainingImages( NUMTRAINIMAGES );
  //applyPCAShapeEstimator->SetNumberOfPrincipalComponentsRequired( NUMLARGESTPC + 1 );
  applyPCAShapeEstimator->SetNumberOfPrincipalComponentsRequired( NUMLARGESTPC );
  applyPCAShapeEstimator->SetInput(0, image1);
  applyPCAShapeEstimator->SetInput(1, image2);
  applyPCAShapeEstimator->SetInput(2, image3);

  applyPCAShapeEstimator->Update();

  //Print the eigen vectors
  vnl_vector<double> eigenValues = applyPCAShapeEstimator->GetEigenValues();
  unsigned int numEigVal =  eigenValues.size();
  std::cout << "Number of returned eign-values: " << numEigVal << std::endl;

  std::cout << "The " << 
    applyPCAShapeEstimator->GetNumberOfPrincipalComponentsRequired() << 
    " largest eigen values are:" << std::endl;

  for(unsigned int i= 0; i< vnl_math_min( numEigVal, (unsigned int)NUMLARGESTPC ); i++ )
    {
    std::cout << eigenValues[ i ] << std::ends; 
    }  
  std::cout << "" << std::endl;
  std::cout << "" << std::endl;

  
  //Print the MeanImage
  OutputImageType::Pointer outImage = applyPCAShapeEstimator->GetOutput( 0 );
  writer->SetFileName( argv[4]);
  writer->SetInput(outImage);
  writer->Update();
  OutputImageIterator outImageIt( outImage, outImage->GetBufferedRegion() );
  outImageIt.GoToBegin();

  //Print the largest two eigen vectors
  for (unsigned int j=1; j< NUMLARGESTPC + 1; j++ )
    {
    OutputImageType::Pointer outImage = applyPCAShapeEstimator->GetOutput( j );
	writer->SetFileName( argv[5] );
    OutputImageIterator outImageIt( outImage, outImage->GetBufferedRegion() );
	outImageIt.GoToBegin();
	writer->SetInput(outImage);
	writer->Update();

    std::cout << "" << std::endl;
    std::cout << "The eigen vector number: " << j << " is:" << std::endl;
    while(!outImageIt.IsAtEnd() )
      {
      std::cout << (double) (outImageIt.Get()) << ";"  << std::ends;  
      ++outImageIt; 
      } 
    std::cout << "  " << std::endl;

    }

  //Test for the eigen values for the test case precomputed using Matlab/Splus
  std::cout << "" << std::endl;
  if( (eigenValues[2] < 6 || eigenValues[2] > 6.1) || (eigenValues[1] >0.1) )
    std::cout<< "Test Passed" << std::endl;
  else 
    std::cout<< "Test failed" << std::endl;


  return 0;
}
-------------- next part --------------
A non-text attachment was scrubbed...
Name: ASM3-1.jpg
Type: image/pjpeg
Size: 852 bytes
Desc: not available
Url : http://public.kitware.com/pipermail/insight-users/attachments/20041214/aae13e1f/ASM3-1.bin
-------------- next part --------------
A non-text attachment was scrubbed...
Name: ASM3-2.jpg
Type: image/tiff
Size: 1668 bytes
Desc: not available
Url : http://public.kitware.com/pipermail/insight-users/attachments/20041214/aae13e1f/ASM3-2.tiff
-------------- next part --------------
A non-text attachment was scrubbed...
Name: ASM3-4.jpg
Type: image/tiff
Size: 1668 bytes
Desc: not available
Url : http://public.kitware.com/pipermail/insight-users/attachments/20041214/aae13e1f/ASM3-4.tiff