[Insight-users] a question about itkImagePCAShapeModelEstimator
Luis Ibanez
luis.ibanez at kitware.com
Sat Dec 6 19:03:54 EST 2008
Hi Yang
Thanks for your detailed description of the problem and
for making your input data available.
An important question:
How do you conclude that the resulting
principal components are incorrect ?
Please note that this class takes all the image pixels
as independent components of a large vector, and performs
PCA on that collection of vectors (one vector per image).
There is therefore the assumption that the images are somehow
normalized, both in intensities and in the scaling and orientation
of the shapes contained inside of the images.
Please let us know about the question above,
Thanks
Luis
----------------
yang.xu wrote:
> 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
>
>
>
>
>
>
>
>
> _______________________________________________
> Insight-users mailing list
> Insight-users at itk.org
> http://www.itk.org/mailman/listinfo/insight-users
>
More information about the Insight-users
mailing list