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
19 //
20 // This example illustrates how to do registration with a 2D Rigid Transform
21 // and with the Normalized Mutual Information metric.
22 //
23
24
25 #include "itkImageRegistrationMethod.h"
26
27 #include "itkCenteredRigid2DTransform.h"
28 #include "itkCenteredTransformInitializer.h"
29
30 #include "itkNormalizedMutualInformationHistogramImageToImageMetric.h"
31 #include "itkOnePlusOneEvolutionaryOptimizer.h"
32 #include "itkNormalVariateGenerator.h"
33
34 #include "itkImageFileReader.h"
35 #include "itkImageFileWriter.h"
36
37 #include "itkResampleImageFilter.h"
38 #include "itkCastImageFilter.h"
39
40
41 // The following section of code implements a Command observer
42 // used to monitor the evolution of the registration process.
43 //
44 #include "itkCommand.h"
45 class CommandIterationUpdate : public itk::Command
46 {
47 public:
48 using Self = CommandIterationUpdate;
49 using Superclass = itk::Command;
50 using Pointer = itk::SmartPointer<Self>;
51 itkNewMacro( Self );
52
53 protected:
CommandIterationUpdate()54 CommandIterationUpdate() {m_LastMetricValue = 0;}
55
56 public:
57 using OptimizerType = itk::OnePlusOneEvolutionaryOptimizer;
58 using OptimizerPointer = const OptimizerType *;
59
Execute(itk::Object * caller,const itk::EventObject & event)60 void Execute(itk::Object *caller,
61 const itk::EventObject & event) override
62 {
63 Execute( (const itk::Object *)caller, event);
64 }
65
Execute(const itk::Object * object,const itk::EventObject & event)66 void Execute(const itk::Object * object,
67 const itk::EventObject & event) override
68 {
69 auto optimizer = static_cast< OptimizerPointer >( object );
70 if( ! itk::IterationEvent().CheckEvent( &event ) )
71 {
72 return;
73 }
74 double currentValue = optimizer->GetValue();
75 // Only print out when the Metric value changes
76 if( std::fabs( m_LastMetricValue - currentValue ) > 1e-7 )
77 {
78 std::cout << optimizer->GetCurrentIteration() << " ";
79 std::cout << currentValue << " ";
80 std::cout << optimizer->GetFrobeniusNorm() << " ";
81 std::cout << optimizer->GetCurrentPosition() << std::endl;
82 m_LastMetricValue = currentValue;
83 }
84 }
85
86 private:
87 double m_LastMetricValue;
88 };
89
90
91 #include "itkTestDriverIncludeRequiredIOFactories.h"
main(int argc,char * argv[])92 int main( int argc, char *argv[] )
93 {
94 RegisterRequiredFactories();
95 if( argc < 4 )
96 {
97 std::cerr << "Missing Parameters " << std::endl;
98 std::cerr << "Usage: " << argv[0];
99 std::cerr << " fixedImageFile movingImageFile ";
100 std::cerr << "outputImagefile [numberOfHistogramBins] ";
101 std::cerr << "[initialRadius] [epsilon]" << std::endl;
102 std::cerr << "[initialAngle(radians)] [initialTx] [initialTy]"
103 << std::endl;
104 return EXIT_FAILURE;
105 }
106
107 constexpr unsigned int Dimension = 2;
108 using PixelType = unsigned char;
109
110 using FixedImageType = itk::Image< PixelType, Dimension >;
111 using MovingImageType = itk::Image< PixelType, Dimension >;
112
113 using TransformType = itk::CenteredRigid2DTransform< double >;
114
115 using OptimizerType = itk::OnePlusOneEvolutionaryOptimizer;
116 using InterpolatorType = itk::LinearInterpolateImageFunction<
117 MovingImageType,
118 double >;
119 using RegistrationType = itk::ImageRegistrationMethod<
120 FixedImageType,
121 MovingImageType >;
122
123
124 using MetricType =
125 itk::NormalizedMutualInformationHistogramImageToImageMetric<
126 FixedImageType,
127 MovingImageType >;
128
129 TransformType::Pointer transform = TransformType::New();
130 OptimizerType::Pointer optimizer = OptimizerType::New();
131 InterpolatorType::Pointer interpolator = InterpolatorType::New();
132 RegistrationType::Pointer registration = RegistrationType::New();
133
134 registration->SetOptimizer( optimizer );
135 registration->SetTransform( transform );
136 registration->SetInterpolator( interpolator );
137
138
139 MetricType::Pointer metric = MetricType::New();
140 registration->SetMetric( metric );
141
142
143 unsigned int numberOfHistogramBins = 32;
144 if( argc > 4 )
145 {
146 numberOfHistogramBins = std::stoi( argv[4] );
147 std::cout << "Using " << numberOfHistogramBins << " Histogram bins"
148 << std::endl;
149 }
150
151 MetricType::HistogramType::SizeType histogramSize;
152
153 histogramSize.SetSize(2);
154
155 histogramSize[0] = numberOfHistogramBins;
156 histogramSize[1] = numberOfHistogramBins;
157 metric->SetHistogramSize( histogramSize );
158
159 const unsigned int numberOfParameters = transform->GetNumberOfParameters();
160 using ScalesType = MetricType::ScalesType;
161 ScalesType scales( numberOfParameters );
162
163 scales.Fill( 1.0 );
164
165 metric->SetDerivativeStepLengthScales(scales);
166
167 using FixedImageReaderType = itk::ImageFileReader< FixedImageType >;
168 using MovingImageReaderType = itk::ImageFileReader< MovingImageType >;
169
170 FixedImageReaderType::Pointer fixedImageReader = FixedImageReaderType::New();
171 MovingImageReaderType::Pointer movingImageReader
172 = MovingImageReaderType::New();
173
174 fixedImageReader->SetFileName( argv[1] );
175 movingImageReader->SetFileName( argv[2] );
176
177 registration->SetFixedImage( fixedImageReader->GetOutput() );
178 registration->SetMovingImage( movingImageReader->GetOutput() );
179 fixedImageReader->Update();
180 movingImageReader->Update();
181
182 FixedImageType::ConstPointer fixedImage = fixedImageReader->GetOutput();
183 registration->SetFixedImageRegion( fixedImage->GetBufferedRegion() );
184
185 using TransformInitializerType = itk::CenteredTransformInitializer<
186 TransformType, FixedImageType,
187 MovingImageType >;
188 TransformInitializerType::Pointer initializer
189 = TransformInitializerType::New();
190 initializer->SetTransform( transform );
191 initializer->SetFixedImage( fixedImageReader->GetOutput() );
192 initializer->SetMovingImage( movingImageReader->GetOutput() );
193 initializer->GeometryOn();
194 initializer->InitializeTransform();
195
196 double initialAngle = 0.0;
197 if( argc > 7 )
198 {
199 initialAngle = std::stod( argv[7] );
200 }
201 transform->SetAngle( initialAngle );
202 TransformType::OutputVectorType initialTranslation
203 = transform->GetTranslation();
204 if( argc > 9 )
205 {
206 initialTranslation[0] += std::stod( argv[8] );
207 initialTranslation[1] += std::stod( argv[9] );
208 }
209 transform->SetTranslation( initialTranslation );
210
211 using ParametersType = RegistrationType::ParametersType;
212 ParametersType initialParameters = transform->GetParameters();
213 registration->SetInitialTransformParameters( initialParameters );
214 std::cout << "Initial transform parameters = ";
215 std::cout << initialParameters << std::endl;
216
217 using OptimizerScalesType = OptimizerType::ScalesType;
218 OptimizerScalesType optimizerScales( transform->GetNumberOfParameters() );
219
220 FixedImageType::RegionType region = fixedImage->GetLargestPossibleRegion();
221 FixedImageType::SizeType size = region.GetSize();
222 FixedImageType::SpacingType spacing = fixedImage->GetSpacing();
223
224 optimizerScales[0] = 1.0 / 0.1; // make angle move slowly
225 optimizerScales[1] = 10000.0; // prevent the center from moving
226 optimizerScales[2] = 10000.0; // prevent the center from moving
227 optimizerScales[3] = 1.0 / ( 0.1 * size[0] * spacing[0] );
228 optimizerScales[4] = 1.0 / ( 0.1 * size[1] * spacing[1] );
229 std::cout << "optimizerScales = " << optimizerScales << std::endl;
230 optimizer->SetScales( optimizerScales );
231
232 using GeneratorType = itk::Statistics::NormalVariateGenerator;
233 GeneratorType::Pointer generator = GeneratorType::New();
234 generator->Initialize(12345);
235 optimizer->MaximizeOn();
236 optimizer->SetNormalVariateGenerator( generator );
237
238 double initialRadius = 0.05;
239 if( argc > 5 )
240 {
241 initialRadius = std::stod( argv[5] );
242 std::cout << "Using initial radius = " << initialRadius << std::endl;
243 }
244 optimizer->Initialize( initialRadius );
245 double epsilon = 0.001;
246 if( argc > 6 )
247 {
248 epsilon = std::stod( argv[6] );
249 std::cout << "Using epsilon = " << epsilon << std::endl;
250 }
251 optimizer->SetEpsilon( epsilon );
252 optimizer->SetMaximumIteration( 2000 );
253
254 // Create the Command observer and register it with the optimizer.
255 //
256 CommandIterationUpdate::Pointer observer = CommandIterationUpdate::New();
257 optimizer->AddObserver( itk::IterationEvent(), observer );
258 try
259 {
260 registration->Update();
261 std::cout << "Optimizer stop condition: "
262 << registration->GetOptimizer()->GetStopConditionDescription()
263 << std::endl;
264 }
265 catch( itk::ExceptionObject & err )
266 {
267 std::cout << "ExceptionObject caught !" << std::endl;
268 std::cout << err << std::endl;
269 return EXIT_FAILURE;
270 }
271
272 using ParametersType = RegistrationType::ParametersType;
273 ParametersType finalParameters = registration->GetLastTransformParameters();
274 const double finalAngle = finalParameters[0];
275 const double finalRotationCenterX = finalParameters[1];
276 const double finalRotationCenterY = finalParameters[2];
277 const double finalTranslationX = finalParameters[3];
278 const double finalTranslationY = finalParameters[4];
279
280 const unsigned int numberOfIterations = optimizer->GetCurrentIteration();
281 const double bestValue = optimizer->GetValue();
282
283 // Print out results
284 const double finalAngleInDegrees = finalAngle * 180.0 / itk::Math::pi;
285 std::cout << " Result = " << std::endl;
286 std::cout << " Angle (radians) " << finalAngle << std::endl;
287 std::cout << " Angle (degrees) " << finalAngleInDegrees << std::endl;
288 std::cout << " Center X = " << finalRotationCenterX << std::endl;
289 std::cout << " Center Y = " << finalRotationCenterY << std::endl;
290 std::cout << " Translation X = " << finalTranslationX << std::endl;
291 std::cout << " Translation Y = " << finalTranslationY << std::endl;
292 std::cout << " Iterations = " << numberOfIterations << std::endl;
293 std::cout << " Metric value = " << bestValue << std::endl;
294
295 using ResampleFilterType = itk::ResampleImageFilter<
296 MovingImageType, FixedImageType >;
297 TransformType::Pointer finalTransform = TransformType::New();
298 finalTransform->SetParameters( finalParameters );
299 finalTransform->SetFixedParameters( transform->GetFixedParameters() );
300
301 ResampleFilterType::Pointer resample = ResampleFilterType::New();
302 resample->SetTransform( finalTransform );
303 resample->SetInput( movingImageReader->GetOutput() );
304 resample->SetSize( fixedImage->GetLargestPossibleRegion().GetSize() );
305 resample->SetOutputOrigin( fixedImage->GetOrigin() );
306 resample->SetOutputSpacing( fixedImage->GetSpacing() );
307 resample->SetOutputDirection( fixedImage->GetDirection() );
308 resample->SetDefaultPixelValue( 100 );
309
310 using OutputImageType = itk::Image< PixelType, Dimension >;
311 using WriterType = itk::ImageFileWriter< OutputImageType >;
312
313 WriterType::Pointer writer = WriterType::New();
314 writer->SetFileName( argv[3] );
315 writer->SetInput( resample->GetOutput() );
316 writer->Update();
317 return EXIT_SUCCESS;
318 }
319