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