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 // Software Guide : BeginLatex
20 //
21 // This example illustrates how to perform Iterative Closest Point (ICP)
22 // registration in ITK. The main class featured in this section is the
23 // \doxygen{EuclideanDistancePointMetric}.
24 //
25 // Software Guide : EndLatex
26
27 // Software Guide : BeginLatex
28 //
29 // The first step is to include the relevant headers.
30 //
31 // Software Guide : EndLatex
32
33 // Software Guide : BeginCodeSnippet
34 #include "itkTranslationTransform.h"
35 #include "itkEuclideanDistancePointMetric.h"
36 #include "itkLevenbergMarquardtOptimizer.h"
37 #include "itkPointSetToPointSetRegistrationMethod.h"
38 // Software Guide : EndCodeSnippet
39
40 #include <iostream>
41 #include <fstream>
42
43 class CommandIterationUpdate : public itk::Command
44 {
45 public:
46 using Self = CommandIterationUpdate;
47 using Superclass = itk::Command;
48 using Pointer = itk::SmartPointer<Self>;
49 itkNewMacro( Self );
50
51 protected:
52 CommandIterationUpdate() = default;
53
54 public:
55 using OptimizerType = itk::LevenbergMarquardtOptimizer;
56 using OptimizerPointer = const OptimizerType *;
57
Execute(itk::Object * caller,const itk::EventObject & event)58 void Execute(itk::Object *caller, const itk::EventObject & event) override
59 {
60 Execute( (const itk::Object *)caller, event);
61 }
62
Execute(const itk::Object * object,const itk::EventObject & event)63 void Execute(const itk::Object * object, const itk::EventObject & event) override
64 {
65 auto optimizer = dynamic_cast< OptimizerPointer >( object );
66 if( optimizer == nullptr )
67 {
68 itkExceptionMacro( "Could not cast optimizer." );
69 }
70
71 if( ! itk::IterationEvent().CheckEvent( &event ) )
72 {
73 return;
74 }
75
76 std::cout << "Value = " << optimizer->GetCachedValue() << std::endl;
77 std::cout << "Position = " << optimizer->GetCachedCurrentPosition();
78 std::cout << std::endl << std::endl;
79 }
80 };
81
82
main(int argc,char * argv[])83 int main(int argc, char * argv[] )
84 {
85
86 if( argc < 3 )
87 {
88 std::cerr << "Arguments Missing. " << std::endl;
89 std::cerr
90 << "Usage: IterativeClosestPoint1 fixedPointsFile movingPointsFile "
91 << std::endl;
92 return EXIT_FAILURE;
93 }
94
95 // Software Guide : BeginLatex
96 //
97 // Next, define the necessary types for the fixed and moving pointsets and
98 // point containers.
99 //
100 // Software Guide : EndLatex
101
102 // Software Guide : BeginCodeSnippet
103 constexpr unsigned int Dimension = 2;
104
105 using PointSetType = itk::PointSet< float, Dimension >;
106
107 PointSetType::Pointer fixedPointSet = PointSetType::New();
108 PointSetType::Pointer movingPointSet = PointSetType::New();
109
110 using PointType = PointSetType::PointType;
111
112 using PointsContainer = PointSetType::PointsContainer;
113
114 PointsContainer::Pointer fixedPointContainer = PointsContainer::New();
115 PointsContainer::Pointer movingPointContainer = PointsContainer::New();
116
117 PointType fixedPoint;
118 PointType movingPoint;
119 // Software Guide : EndCodeSnippet
120
121 // Read the file containing coordinates of fixed points.
122 std::ifstream fixedFile;
123 fixedFile.open( argv[1] );
124 if( fixedFile.fail() )
125 {
126 std::cerr << "Error opening points file with name : " << std::endl;
127 std::cerr << argv[1] << std::endl;
128 return EXIT_FAILURE;
129 }
130
131 unsigned int pointId = 0;
132 fixedFile >> fixedPoint;
133 while( !fixedFile.eof() )
134 {
135 fixedPointContainer->InsertElement( pointId, fixedPoint );
136 fixedFile >> fixedPoint;
137 pointId++;
138 }
139 fixedPointSet->SetPoints( fixedPointContainer );
140 std::cout <<
141 "Number of fixed Points = " <<
142 fixedPointSet->GetNumberOfPoints() << std::endl;
143
144 // Read the file containing coordinates of moving points.
145 std::ifstream movingFile;
146 movingFile.open( argv[2] );
147 if( movingFile.fail() )
148 {
149 std::cerr << "Error opening points file with name : " << std::endl;
150 std::cerr << argv[2] << std::endl;
151 return EXIT_FAILURE;
152 }
153
154 pointId = 0;
155 movingFile >> movingPoint;
156 while( !movingFile.eof() )
157 {
158 movingPointContainer->InsertElement( pointId, movingPoint );
159 movingFile >> movingPoint;
160 pointId++;
161 }
162 movingPointSet->SetPoints( movingPointContainer );
163 std::cout << "Number of moving Points = "
164 << movingPointSet->GetNumberOfPoints() << std::endl;
165
166 // Software Guide : BeginLatex
167 //
168 // After the points are read in from files, set up the metric type.
169 //
170 // Software Guide : EndLatex
171
172 // Software Guide : BeginCodeSnippet
173 using MetricType = itk::EuclideanDistancePointMetric<
174 PointSetType, PointSetType>;
175
176 MetricType::Pointer metric = MetricType::New();
177 // Software Guide : EndCodeSnippet
178
179 // Software Guide : BeginLatex
180 //
181 // Now, setup the transform, optimizers, and registration method using the
182 // point set types defined earlier.
183 //
184 // Software Guide : EndLatex
185
186 // Software Guide : BeginCodeSnippet
187 using TransformType = itk::TranslationTransform< double, Dimension >;
188
189 TransformType::Pointer transform = TransformType::New();
190
191
192 // Optimizer Type
193 using OptimizerType = itk::LevenbergMarquardtOptimizer;
194
195 OptimizerType::Pointer optimizer = OptimizerType::New();
196 optimizer->SetUseCostFunctionGradient(false);
197
198 // Registration Method
199 using RegistrationType = itk::PointSetToPointSetRegistrationMethod<
200 PointSetType, PointSetType >;
201
202 RegistrationType::Pointer registration = RegistrationType::New();
203
204 // Scale the translation components of the Transform in the Optimizer
205 OptimizerType::ScalesType scales( transform->GetNumberOfParameters() );
206 scales.Fill( 0.01 );
207 // Software Guide : EndCodeSnippet
208
209 // Software Guide : BeginLatex
210 //
211 // Next we setup the convergence criteria, and other properties required
212 // by the optimizer.
213 //
214 // Software Guide : EndLatex
215
216 // Software Guide : BeginCodeSnippet
217 unsigned long numberOfIterations = 100;
218 double gradientTolerance = 1e-5; // convergence criterion
219 double valueTolerance = 1e-5; // convergence criterion
220 double epsilonFunction = 1e-6; // convergence criterion
221
222
223 optimizer->SetScales( scales );
224 optimizer->SetNumberOfIterations( numberOfIterations );
225 optimizer->SetValueTolerance( valueTolerance );
226 optimizer->SetGradientTolerance( gradientTolerance );
227 optimizer->SetEpsilonFunction( epsilonFunction );
228 // Software Guide : EndCodeSnippet
229
230 // Software Guide : BeginLatex
231 //
232 // In this case we start from an identity transform, but in reality the user
233 // will usually be able to provide a better guess than this.
234 //
235 // Software Guide : EndLatex
236
237 // Software Guide : BeginCodeSnippet
238 transform->SetIdentity();
239 // Software Guide : EndCodeSnippet
240
241 registration->SetInitialTransformParameters( transform->GetParameters() );
242
243 // Software Guide : BeginLatex
244 //
245 // Finally, connect all the components required for the registration, and an
246 // observer.
247 //
248 // Software Guide : EndLatex
249
250 // Software Guide : BeginCodeSnippet
251 registration->SetMetric( metric );
252 registration->SetOptimizer( optimizer );
253 registration->SetTransform( transform );
254 registration->SetFixedPointSet( fixedPointSet );
255 registration->SetMovingPointSet( movingPointSet );
256
257 // Connect an observer
258 CommandIterationUpdate::Pointer observer = CommandIterationUpdate::New();
259 optimizer->AddObserver( itk::IterationEvent(), observer );
260 // Software Guide : EndCodeSnippet
261
262 try
263 {
264 registration->Update();
265 }
266 catch( itk::ExceptionObject & e )
267 {
268 std::cerr << e << std::endl;
269 return EXIT_FAILURE;
270 }
271
272 std::cout << "Solution = " << transform->GetParameters() << std::endl;
273 return EXIT_SUCCESS;
274 }
275