[Insight-users] problems about itkImagePCAShapeModelEstimator

Tue Dec 2 06:43:14 EST 2008

I'm trying to make a statistical shape model to do image segmentation. Now I have 10 liver CT datasets and I have tansformed them to meshes and found the correspondence between the landmarks between them. But when I try to make SSM(statistical shape model ), I met a problem.
    I use itkImagePCAShapeModelEstimator to do PCA. First, I transformed three liver dataset
    by SignedMaurerDistanceMapImageFilter. Then I use the itkImagePCAShapeModelEstimator to do PCA in order to get a mean image and some Principal Component images. Unfortunately, I can't get principal componet images. Following is my code .
    #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   3
    #define NUMLARGESTPC     3
- Ignored:
    // 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<8)
    {
    std::cerr << "Usage: " << argv[0] << " InputImage1 InputImage2 InputImage3 OutputMean OutPutImage1 OutImage2 OutImage3\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,image1);
    applyPCAShapeEstimator->SetInput(1,image2);
    applyPCAShapeEstimator->SetInput(2,image3);
    //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[4]);
    writer->Update();
    //Output the Components
    writer->SetInput(applyPCAShapeEstimator->GetOutput(1));
    writer->SetFileName(argv[5]);
    writer->Update();
    writer->SetInput(applyPCAShapeEstimator->GetOutput(2));
    writer->SetFileName(argv[6]);
    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;
    }

    2008-12-02

- Done.
------------------------------------------------------------
From:  "yang.xu" <yang.xu at siat.ac.cn>
To:  "insight-users-request" <insight-users-request at itk.org>
Subject:  Problem about itkImagePCAShapeModelEstimator
Date:  Tue2 Dec 2008 19:04:38 +0800

hi:
I'm trying to make a statistical shape model to do image segmentation. Now I have 10 liver CT datasets and I have tansformed them to meshes and found the correspondence between the landmarks between them. But when I try to make SSM(statistical shape model ), I met a problem.
I use itkImagePCAShapeModelEstimator to do PCA. First, I transformed three liver dataset
by SignedMaurerDistanceMapImageFilter. Then I use the itkImagePCAShapeModelEstimator to do PCA in order to get a mean image and some sigma images. Unfortunately, I can't get correct sigma images. Following is my code .
#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   3
#define NUMLARGESTPC     3
// 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<8)
{
std::cerr << "Usage: " << argv[0] << " InputImage1 InputImage2 InputImage3 OutputMean OutPutImage1 OutImage2 OutImage3\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,image1);
applyPCAShapeEstimator->SetInput(1,image2);
applyPCAShapeEstimator->SetInput(2,image3);
//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[4]);
writer->Update();
//Output the Components
writer->SetInput(applyPCAShapeEstimator->GetOutput(1));
writer->SetFileName(argv[5]);
writer->Update();
writer->SetInput(applyPCAShapeEstimator->GetOutput(2));
writer->SetFileName(argv[6]);
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;
}
2008-12-02

                 Anthony Xu  

2008-12-02