ITK/Examples/WishList/Segmentation/EstimatePCAModel

/* Author: Juan Cardelino <juan.cardelino@gmail.com>

• /

1. if defined(_MSC_VER)
2. pragma warning ( disable : 4786 )
3. endif

1. ifdef __BORLANDC__
2. define ITK_LEAN_AND_MEAN
3. endif

1. include "itkImage.h"

1. include "itkGeodesicActiveContourShapePriorLevelSetImageFilter.h"
2. include "itkChangeInformationImageFilter.h"
3. include "itkBoundedReciprocalImageFilter.h"
4. include "itkImagePCAShapeModelEstimator.h" // Software Guide : EndCodeSnippet

1. include "itkPCAShapeSignedDistanceFunction.h"
2. include "itkEuler2DTransform.h"
3. include "itkShapePriorMAPCostFunction.h"
4. include "itkOnePlusOneEvolutionaryOptimizer.h"
5. include "itkNormalVariateGenerator.h"
6. include "vnl/vnl_sample.h"
7. include "itkNumericSeriesFileNames.h"

1. include "itkCurvatureAnisotropicDiffusionImageFilter.h"
2. include "itkGradientMagnitudeRecursiveGaussianImageFilter.h"
3. include "itkFastMarchingImageFilter.h"
4. include "itkRescaleIntensityImageFilter.h"
5. include "itkBinaryThresholdImageFilter.h"
6. include "itkImageFileReader.h"
7. include "itkImageFileWriter.h"
8. include "itkSpatialFunctionImageEvaluatorFilter.h"

1. include "itkCommand.h"
2. include "itkMultiplyByConstantImageFilter.h"

using namespace std;

int main( int argc, char *argv[] ) { if( argc < 5 ) { std::cerr << "Missing Parameters " << std::endl; std::cerr << "Usage: " << argv[0]; std::cerr << " nbTrain trainFilePattern"; std::cerr << " nbModes modeFilePattern"; std::cerr << std::endl; return 1; }

for(int i=0; i<argc;i++) { cout<<"id: "<<i<<" arg: "<<argv[i]<<endl; } const unsigned int Dimension = 2; typedef float my_PixelType; typedef itk::Image< my_PixelType, Dimension > ImageType;

typedef itk::ImageFileReader< ImageType > ReaderType; typedef itk::ImageFileWriter< ImageType > WriterType; typedef itk::MultiplyByConstantImageFilter < ImageType , double ,ImageType > ScaleType;

int nb_train=atoi(argv[1]);

itk::NumericSeriesFileNames::Pointer fileNamesCreator = itk::NumericSeriesFileNames::New(); std::vector<ImageType::Pointer> trainingImages( nb_train );

fileNamesCreator->SetStartIndex( 0 ); fileNamesCreator->SetEndIndex( nb_train - 1 ); fileNamesCreator->SetSeriesFormat( argv[2] ); const std::vector<std::string> & shapeModeFileNames = fileNamesCreator->GetFileNames();

for ( unsigned int k = 0; k < nb_train; k++ ) { ReaderType::Pointer reader = ReaderType::New(); reader->SetFileName( shapeModeFileNames[k].c_str() ); reader->Update(); trainingImages[k] = reader->GetOutput(); }

typedef itk::ImagePCAShapeModelEstimator<ImageType, ImageType > my_Estimatortype; my_Estimatortype::Pointer filter = my_Estimatortype::New(); filter->SetNumberOfTrainingImages(nb_train); filter->SetNumberOfPrincipalComponentsRequired(2);

for ( unsigned int k = 0; k < nb_train; k++ ) {

filter->SetInput(k, trainingImages[k] ); }

int nb_modes=atoi(argv[3]);

itk::NumericSeriesFileNames::Pointer fileNamesOutCreator = itk::NumericSeriesFileNames::New();

fileNamesOutCreator->SetStartIndex( 0 ); fileNamesOutCreator->SetEndIndex( nb_modes-1 ); fileNamesOutCreator->SetSeriesFormat( argv[4] ); const std::vector<std::string> & outFileNames = fileNamesOutCreator->GetFileNames();

ScaleType::Pointer scaler = ScaleType::New();

filter->Update(); my_Estimatortype::VectorOfDoubleType v=filter->GetEigenValues(); double sv_mean=sqrt(v[0]);

for ( unsigned int k = 0; k < nb_modes; k++ ) { double sv=sqrt(v[k]); double sv_n=sv/sv_mean; //double sv_n=sv; std::cout<<"writing: "<<outFileNames[k]<<std::endl;

std::cout<<"svd["<<k<<"]: "<<sv<<" norm: "<<sv_n<<std::endl; WriterType::Pointer writer = WriterType::New(); writer->SetFileName( outFileNames[k].c_str() ); scaler->SetInput(filter->GetOutput(k)); scaler->SetConstant(sv_n); writer->SetInput( scaler->GetOutput() ); writer->Update(); }

return 0; }