ANIMA  4.0
animaMCMMeasureTest.cxx
Go to the documentation of this file.
6 #include <itkImageRegionIterator.h>
8 #include <itkTimeProbe.h>
9 
10 #include <fstream>
11 #include <tclap/CmdLine.h>
14 #include <animaMCMFileReader.h>
15 
16 int main(int ac, const char** av)
17 {
18  // Parsing arguments
19  TCLAP::CmdLine cmd("INRIA / IRISA - VisAGeS/Empenn Team", ' ', ANIMA_VERSION);
20 
21  // Setting up parameters
22  TCLAP::ValueArg<std::string> inArg("i","im","input image",true,"","input image",cmd);
23  TCLAP::ValueArg<std::string> bvalArg("b","bvals","B-values",true,"","b-values",cmd);
24  TCLAP::ValueArg<std::string> bvecArg("v","bvec","Gradient direction",true,"","gradient directions",cmd);
25  TCLAP::ValueArg<std::string> outArg("o","out","Output file",true,"","output file",cmd);
26 
27  try
28  {
29  cmd.parse(ac,av);
30  }
31  catch (TCLAP::ArgException& e)
32  {
33  std::cerr << "Error: " << e.error() << "for argument " << e.argId() << std::endl;
34  return EXIT_FAILURE;
35  }
36 
37  typedef anima::MCMFileReader <double,3> MCMReaderType;
38  MCMReaderType refReader;
39  refReader.SetFileName(inArg.getValue());
40  refReader.Update();
41 
42  MCMReaderType movingReader;
43  movingReader.SetFileName(inArg.getValue());
44  movingReader.Update();
45 
46  typedef anima::GradientFileReader <vnl_vector_fixed <double,3>, double> GradientReaderType;
47  GradientReaderType gradientReader;
48  gradientReader.SetGradientFileName(bvecArg.getValue());
49  gradientReader.SetBValueBaseString(bvalArg.getValue());
50  gradientReader.SetGradientIndependentNormalization(true);
51 
52  gradientReader.Update();
53 
55  typedef BasicMetricType::Pointer BasicMetricPointer;
56 
57  BasicMetricPointer testBasicMetric = BasicMetricType::New();
58  testBasicMetric->SetFixedImage(refReader.GetModelVectorImage());
59  testBasicMetric->SetMovingImage(movingReader.GetModelVectorImage());
60 
61  typedef MCMReaderType::OutputImageType ImageType;
62  ImageType::RegionType region;
63  region.SetIndex(0,68);
64  region.SetIndex(1,56);
65  region.SetIndex(2,62);
66  region.SetSize(0,5);
67  region.SetSize(1,5);
68  region.SetSize(2,5);
69  testBasicMetric->SetFixedImageRegion(region);
70 
71  typedef anima::MCMLinearInterpolateImageFunction <ImageType> InterpolateFunctionType;
72  InterpolateFunctionType::Pointer interpolator = InterpolateFunctionType::New();
73  interpolator->SetInputImage(movingReader.GetModelVectorImage());
74  interpolator->SetReferenceOutputModel(refReader.GetModelVectorImage()->GetDescriptionModel());
75  testBasicMetric->SetInterpolator(interpolator);
76 
78  typedef MCMMetricType::Pointer MCMMetricPointer;
79  MCMMetricPointer testMetric = MCMMetricType::New();
80 
81  testMetric->SetFixedImage(refReader.GetModelVectorImage());
82  testMetric->SetMovingImage(movingReader.GetModelVectorImage());
83  testMetric->SetFixedImageRegion(region);
84  // Use default small and big delta values
85  testMetric->SetGradientStrengths(gradientReader.GetGradientStrengths());
86  testMetric->SetGradientDirections(gradientReader.GetGradients());
87  testMetric->SetInterpolator(interpolator);
88 
89  typedef anima::MCMCorrelationImageToImageMetric <double, double, 3> MCMCorrelationMetricType;
90  typedef MCMCorrelationMetricType::Pointer MCMCorrelationMetricPointer;
91  MCMCorrelationMetricPointer testCorrelationMetric = MCMCorrelationMetricType::New();
92 
93  testCorrelationMetric->SetFixedImage(refReader.GetModelVectorImage());
94  testCorrelationMetric->SetMovingImage(movingReader.GetModelVectorImage());
95  testCorrelationMetric->SetFixedImageRegion(region);
96  // Use default small and big delta values
97  testCorrelationMetric->SetGradientStrengths(gradientReader.GetGradientStrengths());
98  testCorrelationMetric->SetGradientDirections(gradientReader.GetGradients());
99  testCorrelationMetric->SetInterpolator(interpolator);
100 
101  typedef anima::LogRigid3DTransform <double> TransformType;
102  TransformType::Pointer trsf = TransformType::New();
103  trsf->SetIdentity();
104  ImageType::IndexType centerIndex;
105  centerIndex[0] = 70;
106  centerIndex[1] = 58;
107  centerIndex[2] = 64;
108  TransformType::InputPointType center;
109  movingReader.GetModelVectorImage()->TransformIndexToPhysicalPoint(centerIndex,center);
110  trsf->SetCenter(center);
111  testBasicMetric->SetTransform(trsf);
112  testBasicMetric->SetModelRotation(BasicMetricType::Superclass::PPD);
113  testBasicMetric->PreComputeFixedValues();
114  testMetric->SetTransform(trsf);
115  testMetric->SetModelRotation(MCMMetricType::Superclass::PPD);
116  testMetric->PreComputeFixedValues();
117  testCorrelationMetric->SetTransform(trsf);
118  testCorrelationMetric->SetModelRotation(MCMMetricType::Superclass::PPD);
119  testCorrelationMetric->PreComputeFixedValues();
120 
121  TransformType::ParametersType parameters = trsf->GetParameters();
122  std::ofstream outFile(outArg.getValue().c_str());
123 
124  itk::TimeProbe timerL2Approx, timerL2, timerPairingOneToOne, timerPairing, timerCorrelation, timerCorrelationApprox;
125  outFile.precision(10);
126  //for (int i = -90;i <= 90;++i)
127  for (int i = -100;i <= 100;++i)
128  {
129  //double angle = M_PI * i / 180.0;
130  //parameters[2] = angle;
131  double translation = i / 10.0;
132  parameters[4] = translation;
133 
134  timerL2.Start();
135  testMetric->SetForceApproximation(false);
136  double metricValue = testMetric->GetValue(parameters);
137  timerL2.Stop();
138 
139  timerL2Approx.Start();
140  testMetric->SetForceApproximation(true);
141  double metricValueApprox = testMetric->GetValue(parameters);
142  timerL2Approx.Stop();
143 
144  timerCorrelation.Start();
145  testCorrelationMetric->SetForceApproximation(false);
146  double metricCorrelationValue = testCorrelationMetric->GetValue(parameters);
147  timerCorrelation.Stop();
148 
149  timerCorrelationApprox.Start();
150  testCorrelationMetric->SetForceApproximation(true);
151  double metricApproxCorrelationValue = testCorrelationMetric->GetValue(parameters);
152  timerCorrelationApprox.Stop();
153 
154  timerPairingOneToOne.Start();
155  testBasicMetric->SetOneToOneMapping(true);
156  double basicMetricValueOneToOne = testBasicMetric->GetValue(parameters);
157  timerPairingOneToOne.Stop();
158 
159  timerPairing.Start();
160  testBasicMetric->SetOneToOneMapping(false);
161  double basicMetricValue = testBasicMetric->GetValue(parameters);
162  timerPairing.Stop();
163 
164  outFile << translation << " " << metricValue << " " << metricValueApprox << " "
165  << basicMetricValue << " " << basicMetricValueOneToOne << " "
166  << metricCorrelationValue << " " << metricApproxCorrelationValue << std::endl;
167  }
168 
169  std::cout << "Time L2: " << timerL2.GetTotal() << std::endl;
170  std::cout << "Time L2 approx: " << timerL2Approx.GetTotal() << std::endl;
171  std::cout << "Time pairing: " << timerPairing.GetTotal() << std::endl;
172  std::cout << "Time pairing one to one: " << timerPairingOneToOne.GetTotal() << std::endl;
173  std::cout << "Time correlation: " << timerCorrelation.GetTotal() << std::endl;
174  std::cout << "Time correlation approx: " << timerCorrelationApprox.GetTotal() << std::endl;
175 
176  outFile.close();
177 
178  return EXIT_SUCCESS;
179 }
void SetGradientFileName(std::string fName)
void SetFileName(std::string fileName)
int main(int ac, const char **av)