ANIMA  4.0
animaDenseSVFBMRegistration.cxx
Go to the documentation of this file.
2 
3 #include <tclap/CmdLine.h>
4 
5 #include <itkTimeProbe.h>
6 
7 int main(int argc, const char** argv)
8 {
9  typedef anima::PyramidalDenseSVFMatchingBridge <3> PyramidBMType;
10  typedef PyramidBMType::InputImageType InputImageType;
11  typedef itk::ImageFileReader<InputImageType> ReaderType;
12 
13  // Parsing arguments
14  TCLAP::CmdLine cmd("INRIA / IRISA - VisAGeS/Empenn Team", ' ',ANIMA_VERSION);
15 
16  // Setting up parameters
17  TCLAP::ValueArg<std::string> fixedArg("r","refimage","Fixed image",true,"","fixed image",cmd);
18  TCLAP::ValueArg<std::string> movingArg("m","movingimage","Moving image",true,"","moving image",cmd);
19  TCLAP::ValueArg<std::string> outArg("o","outputimage","Output (registered) image",true,"","output image",cmd);
20  TCLAP::ValueArg<std::string> outputTransformArg("O","outtransform","Output transformation",false,"","output transform",cmd);
21  TCLAP::ValueArg<std::string> blockMaskArg("M","mask-im","Mask image for block generation",false,"","block mask image",cmd);
22 
23  TCLAP::ValueArg<unsigned int> blockSizeArg("","bs","Block size (default: 5)",false,5,"block size",cmd);
24  TCLAP::ValueArg<unsigned int> blockSpacingArg("","sp","Block spacing (default: 2)",false,2,"block spacing",cmd);
25  TCLAP::ValueArg<double> stdevThresholdArg("s","stdev","Threshold block standard deviation (default: 5)",false,5,"block minimal standard deviation",cmd);
26  TCLAP::ValueArg<double> percentageKeptArg("k","per-kept","Percentage of blocks with the highest variance kept (default: 0.8)",false,0.8,"percentage of blocks kept",cmd);
27 
28  TCLAP::ValueArg<unsigned int> blockTransfoArg("t","in-transform","Transformation computed between blocks (0: translation, 1: rigid, 2: affine, 3: directional affine, default: 0)",false,0,"transformation between blocks",cmd);
29  TCLAP::ValueArg<unsigned int> directionArg("d","dir","Affine direction for directional transform output (default: 1 = Y axis)",false,1,"direction of directional affine",cmd);
30  TCLAP::ValueArg<unsigned int> blockMetricArg("","metric","Similarity metric between blocks (0: squared correlation coefficient, 1: correlation coefficient, 2: mean squares, default: 0)",false,0,"similarity metric",cmd);
31  TCLAP::ValueArg<unsigned int> optimizerArg("","opt","Optimizer for optimal block search (0: Exhaustive, 1: Bobyqa, default: 1)",false,1,"optimizer",cmd);
32 
33  TCLAP::ValueArg<unsigned int> maxIterationsArg("","mi","Maximum block match iterations (default: 10)",false,10,"maximum iterations",cmd);
34  TCLAP::ValueArg<double> minErrorArg("","me","Minimal distance between consecutive estimated transforms (default: 0.01)",false,0.01,"minimal distance between transforms",cmd);
35 
36  TCLAP::ValueArg<unsigned int> optimizerMaxIterationsArg("","oi","Maximum iterations for local optimizer (default: 100)",false,100,"maximum local optimizer iterations",cmd);
37 
38  TCLAP::ValueArg<double> searchRadiusArg("","sr","Search radius in pixels (exhaustive search window, rho start for bobyqa / newuoa, default: 2)",false,2,"optimizer search radius",cmd);
39  TCLAP::ValueArg<double> searchAngleRadiusArg("","sar","Search angle radius in degrees (rho start for bobyqa / newuoa, default: 5)",false,5,"optimizer search angle radius",cmd);
40  TCLAP::ValueArg<double> searchScaleRadiusArg("","scr","Search scale radius (rho start for bobyqa / newuoa, default: 0.1)",false,0.1,"optimizer search scale radius",cmd);
41  TCLAP::ValueArg<double> finalRadiusArg("","fr","Final radius (rho end for bobyqa / newuoa, default: 0.001)",false,0.001,"optimizer final radius",cmd);
42  TCLAP::ValueArg<double> searchStepArg("","st","Search step for exhaustive search (default: 1)",false,1,"exhaustive optimizer search step",cmd);
43 
44  TCLAP::ValueArg<double> translateUpperBoundArg("","tub","Upper bound on translation for bobyqa (in voxels, default: 10)",false,10,"Bobyqa translate upper bound",cmd);
45  TCLAP::ValueArg<double> angleUpperBoundArg("","aub","Upper bound on angles for bobyqa (in degrees, default: 180)",false,180,"Bobyqa angle upper bound",cmd);
46  TCLAP::ValueArg<double> scaleUpperBoundArg("","scu","Upper bound on scale for bobyqa (default: 3)",false,3,"Bobyqa scale upper bound",cmd);
47 
48  TCLAP::ValueArg<unsigned int> symmetryArg("","sym-reg","Registration symmetry type (0: asymmetric, 1: symmetric, 2: kissing, default: 0)",false,0,"symmetry type",cmd);
49 
50  TCLAP::ValueArg<unsigned int> agregatorArg("a","agregator","Transformation agregator type (0: Baloo, 1: M-smoother, default: 0)",false,0,"agregator type",cmd);
51  TCLAP::ValueArg<double> extrapolationSigmaArg("","fs","Sigma for extrapolation of local pairings (default: 3)",false,3,"extrapolation sigma",cmd);
52  TCLAP::ValueArg<double> elasticSigmaArg("","es","Sigma for elastic regularization (default: 3)",false,3,"elastic regularization sigma",cmd);
53  TCLAP::ValueArg<double> outlierSigmaArg("","os","Sigma for outlier rejection among local pairings (default: 3)",false,3,"outlier rejection sigma",cmd);
54  TCLAP::ValueArg<double> mEstimateConvergenceThresholdArg("","met","Threshold to consider m-estimator converged (default: 0.01)",false,0.01,"m-estimation convergence threshold",cmd);
55  TCLAP::ValueArg<double> neighborhoodApproximationArg("","na","Half size of the neighborhood approximation (multiplied by extrapolation sigma, default: 2.5)",false,2.5,"half size of neighborhood approximation",cmd);
56  TCLAP::ValueArg<unsigned int> bchOrderArg("b","bch-order","BCH composition order (default: 1)",false,1,"BCH order",cmd);
57  TCLAP::ValueArg<unsigned int> expOrderArg("e","exp-order","Order of field exponentiation approximation (in between 0 and 1, default: 0)",false,0,"exponentiation order",cmd);
58 
59  TCLAP::ValueArg<unsigned int> numPyramidLevelsArg("p","pyr","Number of pyramid levels (default: 3)",false,3,"number of pyramid levels",cmd);
60  TCLAP::ValueArg<unsigned int> lastPyramidLevelArg("l","last-level","Index of the last pyramid level explored (default: 0)",false,0,"last pyramid level",cmd);
61  TCLAP::ValueArg<unsigned int> numThreadsArg("T","threads","Number of execution threads (default: 0 = all cores)",false,0,"number of threads",cmd);
62 
63  try
64  {
65  cmd.parse(argc,argv);
66  }
67  catch (TCLAP::ArgException& e)
68  {
69  std::cerr << "Error: " << e.error() << "for argument " << e.argId() << std::endl;
70  return EXIT_FAILURE;
71  }
72 
73  PyramidBMType::Pointer matcher = PyramidBMType::New();
74 
75  ReaderType::Pointer tmpRead = ReaderType::New();
76  tmpRead->SetFileName(fixedArg.getValue());
77  tmpRead->Update();
78 
79  matcher->SetReferenceImage(tmpRead->GetOutput());
80 
81  tmpRead = ReaderType::New();
82  tmpRead->SetFileName(movingArg.getValue());
83  tmpRead->Update();
84 
85  matcher->SetFloatingImage(tmpRead->GetOutput());
86 
87  // Setting matcher arguments
88  matcher->SetBlockSize( blockSizeArg.getValue() );
89  matcher->SetBlockSpacing( blockSpacingArg.getValue() );
90  matcher->SetStDevThreshold( stdevThresholdArg.getValue() );
91  matcher->SetTransform( (PyramidBMType::Transform) blockTransfoArg.getValue() );
92  matcher->SetAffineDirection(directionArg.getValue());
93  matcher->SetMetric( (PyramidBMType::Metric) blockMetricArg.getValue() );
94  matcher->SetOptimizer( (PyramidBMType::Optimizer) optimizerArg.getValue() );
95  matcher->SetMaximumIterations( maxIterationsArg.getValue() );
96  matcher->SetMinimalTransformError( minErrorArg.getValue() );
97  matcher->SetFinalRadius(finalRadiusArg.getValue());
98  matcher->SetOptimizerMaximumIterations( optimizerMaxIterationsArg.getValue() );
99  matcher->SetSearchRadius( searchRadiusArg.getValue() );
100  matcher->SetSearchAngleRadius( searchAngleRadiusArg.getValue() );
101  matcher->SetSearchScaleRadius( searchScaleRadiusArg.getValue() );
102  matcher->SetStepSize( searchStepArg.getValue() );
103  matcher->SetTranslateUpperBound( translateUpperBoundArg.getValue() );
104  matcher->SetAngleUpperBound( angleUpperBoundArg.getValue() );
105  matcher->SetScaleUpperBound( scaleUpperBoundArg.getValue() );
106  matcher->SetSymmetryType( (PyramidBMType::SymmetryType) symmetryArg.getValue() );
107  matcher->SetAgregator( (PyramidBMType::Agregator) agregatorArg.getValue() );
108  matcher->SetExtrapolationSigma(extrapolationSigmaArg.getValue());
109  matcher->SetElasticSigma(elasticSigmaArg.getValue());
110  matcher->SetOutlierSigma(outlierSigmaArg.getValue());
111  matcher->SetMEstimateConvergenceThreshold(mEstimateConvergenceThresholdArg.getValue());
112  matcher->SetNeighborhoodApproximation(neighborhoodApproximationArg.getValue());
113  matcher->SetBCHCompositionOrder(bchOrderArg.getValue());
114  matcher->SetExponentiationOrder(expOrderArg.getValue());
115  matcher->SetNumberOfPyramidLevels( numPyramidLevelsArg.getValue() );
116  matcher->SetLastPyramidLevel( lastPyramidLevelArg.getValue() );
117 
118  if (blockMaskArg.getValue() != "")
119  matcher->SetBlockGenerationMask(anima::readImage<PyramidBMType::MaskImageType>(blockMaskArg.getValue()));
120 
121  if (numThreadsArg.getValue() != 0)
122  matcher->SetNumberOfWorkUnits( numThreadsArg.getValue() );
123 
124  matcher->SetPercentageKept( percentageKeptArg.getValue() );
125 
126  matcher->SetResultFile( outArg.getValue() );
127  matcher->SetOutputTransformFile( outputTransformArg.getValue() );
128 
129  itk::TimeProbe timer;
130  timer.Start();
131 
132  try
133  {
134  matcher->Update();
135  matcher->WriteOutputs();
136  }
137  catch (itk::ExceptionObject &e)
138  {
139  std::cerr << e << std::endl;
140  return EXIT_FAILURE;
141  }
142 
143  timer.Stop();
144 
145  std::cout << "Elapsed Time: " << timer.GetTotal() << timer.GetUnit() << std::endl;
146 
147  return EXIT_SUCCESS;
148 }
int main(int argc, const char **argv)