[Insight-users] a question about itkImagePCAShapeModelEstimator

[yang.xu]

Hi:everyone.

  I have wrote an email yesterday but I haven't got a reply by now. I think maybe I didn't make my question clear enough to understand. So I write a email now to ask my question again and hope to find a solution.

  I am trying to do liver segmentation based on CT Image datasets. My plan is following:

        1.     use 10 liver CT datasets  to make a statistical shape model
     
        2.     find a method to segment liver based on the statistical shape model (Now I am trying the method suggested by ME.Leventon   (Insight\Examples\Segmentation\GeodesicActiveContourShapePriorLevelSetImageFilter.cxx)

But when I making   statistical shape model , I found the result produced by itkImagePCAShapeModelEstimator may not correct. The mean model result is correct but I can't get correct principal component models. In order to speed up the calculation, the following code  is only based on 5 CT liver dataset. The input of my program are 5 signed distance image datasets. And their link is http://www.siat.ac.cn/~jiafucang/ssmdata 

// ImagePCAShapeModel.cpp : Defines the entry point for the console application.
//       based on itkImagePCAShapeModelEstimatortest.cxx
#if defined(_MSC_VER)
#pragma warning(disable : 4786)
#endif

#include "itkImage.h"
#include "itkVector.h"

#include "itkImageRegionIterator.h"
#include "itkLightProcessObject.h"
#include "itkTextOutput.h"
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"

#include "itkImagePCAShapeModelEstimator.h"

//Data definitions
#define IMGWIDTH		 2
#define IMGHEIGHT		 2
#define NDIMENSION       3
// #define NUMTRAINIMAGES   13
#define NUMTRAINIMAGES   5
#define NUMLARGESTPC     2

// class to support progress feeback

//class ShowProgressObject
//{
//public:
//	ShowProgressObject(itk::LightProcessObject *o)
//	  {m_Process = o;}
//    void ShowProgress()
//	{std::cout << "Progress " << m_Process->GetProgress() <<std::endl;}
//	itk::LightProcessObject::Pointer m_Process;
//
//};

int main(int argc,char* argv[])
{
	if(argc<9)
	{
		std::cerr << "Usage: " << argv[0] << " InputImage1 InputImage2 InputImage3 InputImage4 InputImage5 OutputMean OutPutImage1 OutImage2 \n";
		return -1;
	}
	
//	itk::OutputWindow::SetInstance(itk::TextOutput::New().GetPointer());

	const unsigned int ImageDimension=3;
	typedef float InputPixelType;
	typedef float OutputPixelType;

	typedef itk::Image<InputPixelType, ImageDimension> InputImageType;
	typedef itk::Image<OutputPixelType,ImageDimension> OutputImageType;
	typedef itk::ImageFileReader<InputImageType> ReaderType;
	typedef itk::ImageFileWriter<OutputImageType> WriterType;
	typedef InputImageType::SizeType InputSizeType;


	//InputImageType::Pointer image1 = InputImageType::New();
	//InputImageType::Pointer image2 = InputImageType::New();
	//InputImageType::Pointer image3 = InputImageType::New();



	ReaderType::Pointer reader1 = ReaderType::New();
	reader1->SetFileName(argv[1]);
	reader1->Update();
	image1 = reader1->GetOutput(); 
	

	ReaderType::Pointer reader2 = ReaderType::New();
	reader2->SetFileName(argv[2]);
	reader2->Update();
	image2 = reader2->GetOutput();

	ReaderType::Pointer reader3= ReaderType::New();
	reader3->SetFileName(argv[3]);
	reader3->Update();
	image3=reader3->GetOutput();

	ReaderType::Pointer reader4= ReaderType::New();
	reader4->SetFileName(argv[4]);
	reader4->Update();

	ReaderType::Pointer reader5= ReaderType::New();
	reader5->SetFileName(argv[5]);
	reader5->Update();

	//ReaderType::Pointer reader6= ReaderType::New();
	//reader6->SetFileName(argv[6]);
	//reader6->Update();

	//ReaderType::Pointer reader7= ReaderType::New();
	//reader7->SetFileName(argv[7]);
	//reader7->Update();

	//ReaderType::Pointer reader8= ReaderType::New();
	//reader8->SetFileName(argv[8]);
	//reader8->Update();

	//ReaderType::Pointer reader9= ReaderType::New();
	//reader9->SetFileName(argv[9]);
	//reader9->Update();

	//ReaderType::Pointer reader10= ReaderType::New();
	//reader10->SetFileName(argv[10]);
	//reader10->Update();

	//ReaderType::Pointer reader11= ReaderType::New();
	//reader11->SetFileName(argv[11]);
	//reader11->Update();

	//ReaderType::Pointer reader12= ReaderType::New();
	//reader12->SetFileName(argv[12]);
	//reader12->Update();

	//ReaderType::Pointer reader13= ReaderType::New();
	//reader13->SetFileName(argv[13]);
	//reader13->Update();

	typedef itk::ImagePCAShapeModelEstimator<InputImageType, OutputImageType>
		ImagePCAShapeModelEstimatorType;
    
	ImagePCAShapeModelEstimatorType::Pointer
		applyPCAShapeEstimator = ImagePCAShapeModelEstimatorType::New();
	

	applyPCAShapeEstimator->SetNumberOfTrainingImages(NUMTRAINIMAGES);
	// applyPCAShapeEstimator->SetNumberOfPrincipalComponentsRequired(NUMLARGESTPC+1);
	applyPCAShapeEstimator->SetNumberOfPrincipalComponentsRequired(NUMLARGESTPC);
	

	applyPCAShapeEstimator->SetInput(0,reader1->GetOutput());
	applyPCAShapeEstimator->SetInput(1,reader2->GetOutput());
	applyPCAShapeEstimator->SetInput(2,reader3->GetOutput());
	applyPCAShapeEstimator->SetInput(3,reader4->GetOutput());
	applyPCAShapeEstimator->SetInput(4,reader5->GetOutput());
	//applyPCAShapeEstimator->SetInput(5,reader6->GetOutput());
	//applyPCAShapeEstimator->SetInput(6,reader7->GetOutput());
	//applyPCAShapeEstimator->SetInput(7,reader8->GetOutput());
	//applyPCAShapeEstimator->SetInput(8,reader9->GetOutput());
	//applyPCAShapeEstimator->SetInput(9,reader10->GetOutput());
	//applyPCAShapeEstimator->SetInput(10,reader11->GetOutput());
	//applyPCAShapeEstimator->SetInput(11,reader12->GetOutput());
	//applyPCAShapeEstimator->SetInput(12,reader13->GetOutput());


		applyPCAShapeEstimator->Update();
	
	
	//Test the Printfself function to increase coverage
	applyPCAShapeEstimator->Print(std::cout);

	//Exercise TypeMacro in superclass
	typedef ImagePCAShapeModelEstimatorType::Superclass GenericEstimatorType;
	std::cout<< applyPCAShapeEstimator->GenericEstimatorType::GetNameOfClass()<<std::endl;

	//Print out the number of training images and the number of principal 
	//components
	std::cout << "The number of training images are: " <<
		applyPCAShapeEstimator->GetNumberOfTrainingImages() << std::endl;

	std::cout << "The number of principal components desired are: " <<
		applyPCAShapeEstimator->GetNumberOfPrincipalComponentsRequired() << std::endl;

	//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::endl;
	}
	std::cout << "" << std::endl;
	std::cout << "" << std::endl;

	//Output the MeanImage
	WriterType::Pointer writer = WriterType::New();
	writer->SetInput(applyPCAShapeEstimator->GetOutput(0));
	writer->SetFileName(argv[6]);
	writer->Update();

	//Output the Components
	writer->SetInput(applyPCAShapeEstimator->GetOutput(1));
	writer->SetFileName(argv[7]);
	writer->Update();

	writer->SetInput(applyPCAShapeEstimator->GetOutput(2));
	writer->SetFileName(argv[8]);
	writer->Update();

///*	writer->SetInput(applyPCAShapeEstimator->GetOutput(3));
//	writer->SetFileName(argv[7]);
//	writer->Update();


	//Print the largest two eigen vectors
	typedef
		itk::ImageRegionIterator<OutputImageType > OutputImageIterator;

	for(unsigned int j = 1; j< NUMLARGESTPC + 1; j++)
	{
		OutputImageType::Pointer outImage2 = applyPCAShapeEstimator->GetOutput(j);
		OutputImageIterator outImage2It( outImage2, outImage2->GetBufferedRegion());
		outImage2It.GoToBegin();
		
		std::cout << "" << std::endl;
		std::cout << "The eigen vector number: "<< j <<" is:" << std::endl;
		while(!outImage2It.IsAtEnd())
		{
			std::cout << (double)(outImage2It.Get()) << ":" << std::endl;
			++outImage2It;
		}
		std::cout << "	"<<std::endl;
	}


	std::cout<<"PCA PASSED!"<<std::endl;



	return 0;
}

 Hoping to get some advices.Thank you very much!

      Anthony Xu
     2008-12-03