1 /*=========================================================================
2 *
3 * Copyright Insight Software Consortium
4 *
5 * Licensed under the Apache License, Version 2.0 (the "License");
6 * you may not use this file except in compliance with the License.
7 * You may obtain a copy of the License at
8 *
9 * http://www.apache.org/licenses/LICENSE-2.0.txt
10 *
11 * Unless required by applicable law or agreed to in writing, software
12 * distributed under the License is distributed on an "AS IS" BASIS,
13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 * See the License for the specific language governing permissions and
15 * limitations under the License.
16 *
17 *=========================================================================*/
18 #include "itkQuasiNewtonOptimizerv4.h"
19 #include "itkMeanSquaresImageToImageMetricv4.h"
20 #include "itkRegistrationParameterScalesFromPhysicalShift.h"
21 #include "itkRegistrationParameterScalesFromIndexShift.h"
22 #include "itkRegistrationParameterScalesFromJacobian.h"
23
24 #include "itkImageRegistrationMethodImageSource.h"
25
26 #include "itkIdentityTransform.h"
27 #include "itkAffineTransform.h"
28 #include "itkTranslationTransform.h"
29 #include "itkTestingMacros.h"
30
31 /**
32 * This is a test using GradientDescentOptimizerv4 and parameter scales
33 * estimator. The scales are estimated before the first iteration by
34 * RegistrationParameterScalesFromShift. The learning rates are estimated
35 * at each iteration according to the shift of each step.
36 */
37
38 template< typename TMovingTransform >
itkQuasiNewtonOptimizerv4TestTemplated(int numberOfIterations,double shiftOfStep,std::string scalesOption,bool usePhysicalSpaceForShift=true)39 int itkQuasiNewtonOptimizerv4TestTemplated(int numberOfIterations,
40 double shiftOfStep,
41 std::string scalesOption,
42 bool usePhysicalSpaceForShift = true)
43 {
44 const unsigned int Dimension = TMovingTransform::SpaceDimension;
45 using PixelType = double;
46
47 // Fixed Image Type
48 using FixedImageType = itk::Image<PixelType,Dimension>;
49
50 // Moving Image Type
51 using MovingImageType = itk::Image<PixelType,Dimension>;
52
53 // Size Type
54 using SizeType = typename MovingImageType::SizeType;
55
56 // ImageSource
57 using ImageSourceType = typename itk::testhelper::ImageRegistrationMethodImageSource<
58 typename FixedImageType::PixelType,
59 typename MovingImageType::PixelType,
60 Dimension >;
61
62 typename FixedImageType::ConstPointer fixedImage;
63 typename MovingImageType::ConstPointer movingImage;
64 typename ImageSourceType::Pointer imageSource;
65
66 imageSource = ImageSourceType::New();
67
68 SizeType size;
69 size[0] = 100;
70 size[1] = 100;
71
72 imageSource->GenerateImages( size );
73
74 fixedImage = imageSource->GetFixedImage();
75 movingImage = imageSource->GetMovingImage();
76
77 // Transform for the moving image
78 using MovingTransformType = TMovingTransform;
79 typename MovingTransformType::Pointer movingTransform = MovingTransformType::New();
80 movingTransform->SetIdentity();
81
82 // Transform for the fixed image
83 using FixedTransformType = itk::IdentityTransform<double, Dimension>;
84 typename FixedTransformType::Pointer fixedTransform = FixedTransformType::New();
85 fixedTransform->SetIdentity();
86
87 // ParametersType for the moving transform
88 using ParametersType = typename MovingTransformType::ParametersType;
89
90 // Metric
91 using MetricType = itk::MeanSquaresImageToImageMetricv4
92 < FixedImageType, MovingImageType, FixedImageType >;
93 typename MetricType::Pointer metric = MetricType::New();
94
95 // Assign images and transforms to the metric.
96 metric->SetFixedImage( fixedImage );
97 metric->SetMovingImage( movingImage );
98 metric->SetVirtualDomainFromImage( const_cast<FixedImageType *>(fixedImage.GetPointer()) );
99
100 metric->SetFixedTransform( fixedTransform );
101 metric->SetMovingTransform( movingTransform );
102
103 // Initialize the metric to prepare for use
104 metric->Initialize();
105
106 // Optimizer
107 using OptimizerType = itk::QuasiNewtonOptimizerv4;
108 OptimizerType::Pointer optimizer = OptimizerType::New();
109
110 optimizer->SetMetric( metric );
111 optimizer->SetNumberOfIterations( numberOfIterations );
112
113 // Instantiate an Observer to report the progress of the Optimization
114 using CommandIterationType = itk::CommandIterationUpdate< OptimizerType >;
115 CommandIterationType::Pointer iterationCommand = CommandIterationType::New();
116 iterationCommand->SetOptimizer( optimizer );
117
118 // Optimizer parameter scales estimator
119 typename itk::OptimizerParameterScalesEstimator::Pointer scalesEstimator;
120
121 using PhysicalShiftScalesEstimatorType = itk::RegistrationParameterScalesFromPhysicalShift< MetricType >;
122 using IndexShiftScalesEstimatorType = itk::RegistrationParameterScalesFromIndexShift< MetricType >;
123 using JacobianScalesEstimatorType = itk::RegistrationParameterScalesFromJacobian< MetricType >;
124
125 if (scalesOption.compare("shift") == 0)
126 {
127 if( usePhysicalSpaceForShift )
128 {
129 std::cout << "Testing RegistrationParameterScalesFrom*Physical*Shift" << std::endl;
130 typename PhysicalShiftScalesEstimatorType::Pointer shiftScalesEstimator = PhysicalShiftScalesEstimatorType::New();
131 shiftScalesEstimator->SetMetric(metric);
132 shiftScalesEstimator->SetTransformForward(true); //default
133 scalesEstimator = shiftScalesEstimator;
134 }
135 else
136 {
137 std::cout << "Testing RegistrationParameterScalesFrom*Index*Shift" << std::endl;
138 typename IndexShiftScalesEstimatorType::Pointer shiftScalesEstimator = IndexShiftScalesEstimatorType::New();
139 shiftScalesEstimator->SetMetric(metric);
140 shiftScalesEstimator->SetTransformForward(true); //default
141 scalesEstimator = shiftScalesEstimator;
142 }
143 }
144 else
145 {
146 std::cout << "Testing RegistrationParameterScalesFromJacobian" << std::endl;
147 typename JacobianScalesEstimatorType::Pointer jacobianScalesEstimator
148 = JacobianScalesEstimatorType::New();
149 jacobianScalesEstimator->SetMetric(metric);
150 scalesEstimator = jacobianScalesEstimator;
151 }
152
153 optimizer->SetScalesEstimator(scalesEstimator);
154 // If SetTrustedStepScale is not called, it will use voxel spacing.
155 optimizer->SetMaximumStepSizeInPhysicalUnits(shiftOfStep);
156 optimizer->SetMaximumNewtonStepSizeInPhysicalUnits(shiftOfStep*3.0);
157
158 std::cout << "Start optimization..." << std::endl
159 << "Number of iterations: " << numberOfIterations << std::endl;
160
161 try
162 {
163 optimizer->StartOptimization();
164 }
165 catch( itk::ExceptionObject & e )
166 {
167 std::cout << "Exception thrown ! " << std::endl;
168 std::cout << "An error occurred during Optimization:" << std::endl;
169 std::cout << e.GetLocation() << std::endl;
170 std::cout << e.GetDescription() << std::endl;
171 std::cout << e.what() << std::endl;
172 return EXIT_FAILURE;
173 }
174
175 std::cout << "...finished. " << std::endl
176 << "StopCondition: " << optimizer->GetStopConditionDescription()
177 << std::endl
178 << "Metric: NumberOfValidPoints: "
179 << metric->GetNumberOfValidPoints()
180 << std::endl;
181
182 //
183 // results
184 //
185 ParametersType finalParameters = movingTransform->GetParameters();
186 ParametersType fixedParameters = movingTransform->GetFixedParameters();
187 std::cout << "Estimated scales = " << optimizer->GetScales() << std::endl;
188 std::cout << "finalParameters = " << finalParameters << std::endl;
189 std::cout << "fixedParameters = " << fixedParameters << std::endl;
190 bool pass = true;
191
192 ParametersType actualParameters = imageSource->GetActualParameters();
193 std::cout << "actualParameters = " << actualParameters << std::endl;
194 const unsigned int numbeOfParameters = actualParameters.Size();
195
196 // We know that for the Affine transform the Translation parameters are at
197 // the end of the list of parameters.
198 const unsigned int offsetOrder = finalParameters.Size()-actualParameters.Size();
199
200 constexpr double tolerance = 1.0; // equivalent to 1 pixel.
201
202 for(unsigned int i=0; i<numbeOfParameters; i++)
203 {
204 // the parameters are negated in order to get the inverse transformation.
205 // this only works for comparing translation parameters....
206 std::cout << finalParameters[i+offsetOrder] << " == " << -actualParameters[i] << std::endl;
207 if( itk::Math::abs ( finalParameters[i+offsetOrder] - (-actualParameters[i]) ) > tolerance )
208 {
209 std::cout << "Tolerance exceeded at component " << i << std::endl;
210 pass = false;
211 }
212 }
213
214 if( !pass )
215 {
216 std::cout << "Test FAILED." << std::endl;
217 return EXIT_FAILURE;
218 }
219 else
220 {
221 std::cout << "Test PASSED." << std::endl;
222 return EXIT_SUCCESS;
223 }
224 }
225
itkQuasiNewtonOptimizerv4Test(int argc,char ** const argv)226 int itkQuasiNewtonOptimizerv4Test(int argc, char ** const argv)
227 {
228 if( argc > 3 )
229 {
230 std::cerr << "Missing Parameters " << std::endl;
231 std::cerr << "Usage: " << itkNameOfTestExecutableMacro(argv);
232 std::cerr << " [numberOfIterations=50 shiftOfStep=1] ";
233 std::cerr << std::endl;
234 return EXIT_FAILURE;
235 }
236 unsigned int numberOfIterations = 50;
237 double shiftOfStep = 1.0;
238
239 if( argc >= 2 )
240 {
241 numberOfIterations = std::stoi( argv[1] );
242 }
243 if (argc >= 3)
244 {
245 shiftOfStep = std::stod( argv[2] );
246 }
247
248 constexpr unsigned int Dimension = 2;
249
250 std::cout << std::endl << "Optimizing translation transform with shift scales" << std::endl;
251 using TranslationTransformType = itk::TranslationTransform<double, Dimension>;
252 int ret1 = itkQuasiNewtonOptimizerv4TestTemplated<TranslationTransformType>(numberOfIterations, shiftOfStep, "shift");
253
254 std::cout << std::endl << "Optimizing translation transform with Jacobian scales" << std::endl;
255 using TranslationTransformType = itk::TranslationTransform<double, Dimension>;
256 int ret2 = itkQuasiNewtonOptimizerv4TestTemplated<TranslationTransformType>(numberOfIterations, shiftOfStep, "jacobian");
257
258 if ( ret1 == EXIT_SUCCESS && ret2 == EXIT_SUCCESS )
259 {
260 std::cout << std::endl << "Tests PASSED." << std::endl;
261 return EXIT_SUCCESS;
262 }
263 else
264 {
265 std::cout << std::endl << "Tests FAILED." << std::endl;
266 return EXIT_FAILURE;
267 }
268 }
269