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 using a DistanceMap in order to increase the performance.
23 // There is of course a trade-off between the time needed for computing the
24 // DistanceMap and the time saved by its repeated use during the
25 // iterative computation of the point-to-point distances. It is then necessary
26 // in practice to ponder both factors.
27 //
28 // \doxygen{EuclideanDistancePointMetric}.
29 //
30 // Software Guide : EndLatex
31
32 // Software Guide : BeginLatex
33 //
34 // The first step is to include the relevant headers.
35 //
36 // Software Guide : EndLatex
37
38 // Software Guide : BeginCodeSnippet
39 #include "itkTranslationTransform.h"
40 #include "itkEuclideanDistancePointMetric.h"
41 #include "itkLevenbergMarquardtOptimizer.h"
42 #include "itkPointSetToPointSetRegistrationMethod.h"
43 #include "itkDanielssonDistanceMapImageFilter.h"
44 #include "itkPointSetToImageFilter.h"
45 #include <iostream>
46 #include <fstream>
47 // Software Guide : EndCodeSnippet
48
main(int argc,char * argv[])49 int main(int argc, char * argv[] )
50 {
51
52 if( argc < 3 )
53 {
54 std::cerr << "Arguments Missing. " << std::endl;
55 std::cerr <<
56 "Usage: IterativeClosestPoint3 fixedPointsFile movingPointsFile "
57 << std::endl;
58 return EXIT_FAILURE;
59 }
60
61 constexpr unsigned int Dimension = 2;
62
63 // Software Guide : BeginLatex
64 //
65 // Next, define the necessary types for the fixed and moving point sets.
66 //
67 // Software Guide : EndLatex
68
69 // Software Guide : BeginCodeSnippet
70 using PointSetType = itk::PointSet< float, Dimension >;
71
72 PointSetType::Pointer fixedPointSet = PointSetType::New();
73 PointSetType::Pointer movingPointSet = PointSetType::New();
74
75 using PointType = PointSetType::PointType;
76
77 using PointsContainer = PointSetType::PointsContainer;
78
79 PointsContainer::Pointer fixedPointContainer = PointsContainer::New();
80 PointsContainer::Pointer movingPointContainer = PointsContainer::New();
81
82 PointType fixedPoint;
83 PointType movingPoint;
84 // Software Guide : EndCodeSnippet
85
86 // Read the file containing coordinates of fixed points.
87 std::ifstream fixedFile;
88 fixedFile.open( argv[1] );
89 if( fixedFile.fail() )
90 {
91 std::cerr << "Error opening points file with name : " << std::endl;
92 std::cerr << argv[1] << std::endl;
93 return EXIT_FAILURE;
94 }
95
96 unsigned int pointId = 0;
97 fixedFile >> fixedPoint;
98 while( !fixedFile.eof() )
99 {
100 fixedPointContainer->InsertElement( pointId, fixedPoint );
101 fixedFile >> fixedPoint;
102 pointId++;
103 }
104 fixedPointSet->SetPoints( fixedPointContainer );
105 std::cout << "Number of fixed Points = "
106 << fixedPointSet->GetNumberOfPoints() << std::endl;
107
108 // Read the file containing coordinates of moving points.
109 std::ifstream movingFile;
110 movingFile.open( argv[2] );
111 if( movingFile.fail() )
112 {
113 std::cerr << "Error opening points file with name : " << std::endl;
114 std::cerr << argv[2] << std::endl;
115 return EXIT_FAILURE;
116 }
117
118 pointId = 0;
119 movingFile >> movingPoint;
120 while( !movingFile.eof() )
121 {
122 movingPointContainer->InsertElement( pointId, movingPoint );
123 movingFile >> movingPoint;
124 pointId++;
125 }
126 movingPointSet->SetPoints( movingPointContainer );
127 std::cout << "Number of moving Points = "
128 << movingPointSet->GetNumberOfPoints() << std::endl;
129
130 // Software Guide : BeginLatex
131 //
132 // Setup the metric, transform, optimizers and registration in a manner
133 // similar to the previous two examples.
134 //
135 // Software Guide : EndLatex
136
137 using MetricType = itk::EuclideanDistancePointMetric<
138 PointSetType,
139 PointSetType>;
140
141 MetricType::Pointer metric = MetricType::New();
142
143 //-----------------------------------------------------------
144 // Set up a Transform
145 //-----------------------------------------------------------
146
147 using TransformType = itk::TranslationTransform< double, Dimension >;
148
149 TransformType::Pointer transform = TransformType::New();
150
151 // Optimizer Type
152 using OptimizerType = itk::LevenbergMarquardtOptimizer;
153
154 OptimizerType::Pointer optimizer = OptimizerType::New();
155 optimizer->SetUseCostFunctionGradient(false);
156
157 // Registration Method
158 using RegistrationType = itk::PointSetToPointSetRegistrationMethod<
159 PointSetType,
160 PointSetType >;
161
162 RegistrationType::Pointer registration = RegistrationType::New();
163
164 // Scale the translation components of the Transform in the Optimizer
165 OptimizerType::ScalesType scales( transform->GetNumberOfParameters() );
166 scales.Fill( 0.01 );
167
168 constexpr unsigned long numberOfIterations = 100;
169 const double gradientTolerance = 1e-5; // convergence criterion
170 const double valueTolerance = 1e-5; // convergence criterion
171 const double epsilonFunction = 1e-6; // convergence criterion
172
173 optimizer->SetScales( scales );
174 optimizer->SetNumberOfIterations( numberOfIterations );
175 optimizer->SetValueTolerance( valueTolerance );
176 optimizer->SetGradientTolerance( gradientTolerance );
177 optimizer->SetEpsilonFunction( epsilonFunction );
178
179 // Start from an Identity transform (in a normal case, the user
180 // can probably provide a better guess than the identity...
181 transform->SetIdentity();
182
183 registration->SetInitialTransformParameters( transform->GetParameters() );
184
185 //------------------------------------------------------
186 // Connect all the components required for Registration
187 //------------------------------------------------------
188 registration->SetMetric( metric );
189 registration->SetOptimizer( optimizer );
190 registration->SetTransform( transform );
191 registration->SetFixedPointSet( fixedPointSet );
192 registration->SetMovingPointSet( movingPointSet );
193
194 //------------------------------------------------------
195 // Prepare the Distance Map in order to accelerate
196 // distance computations.
197 //------------------------------------------------------
198 //
199 // First map the Fixed Points into a binary image.
200 // This is needed because the DanielssonDistance
201 // filter expects an image as input.
202 //
203 //-------------------------------------------------
204 // Software Guide : BeginLatex
205 //
206 // In the preparation of the distance map, we first need to map the fixed
207 // points into a binary image.
208 //
209 // Software Guide : EndLatex
210
211 // Software Guide : BeginCodeSnippet
212 using BinaryImageType = itk::Image< unsigned char, Dimension >;
213
214 using PointsToImageFilterType = itk::PointSetToImageFilter<
215 PointSetType,
216 BinaryImageType>;
217
218 PointsToImageFilterType::Pointer
219 pointsToImageFilter = PointsToImageFilterType::New();
220
221 pointsToImageFilter->SetInput( fixedPointSet );
222
223 BinaryImageType::SpacingType spacing;
224 spacing.Fill( 1.0 );
225
226 BinaryImageType::PointType origin;
227 origin.Fill( 0.0 );
228 // Software Guide : EndCodeSnippet
229
230 // Software Guide : BeginLatex
231 //
232 // Continue to prepare the distance map, in order to accelerate the distance
233 // computations.
234 //
235 // Software Guide : EndLatex
236
237 // Software Guide : BeginCodeSnippet
238 pointsToImageFilter->SetSpacing( spacing );
239 pointsToImageFilter->SetOrigin( origin );
240 pointsToImageFilter->Update();
241 BinaryImageType::Pointer binaryImage = pointsToImageFilter->GetOutput();
242
243 using DistanceImageType = itk::Image< unsigned short, Dimension >;
244 using DistanceFilterType = itk::DanielssonDistanceMapImageFilter<
245 BinaryImageType, DistanceImageType>;
246
247 DistanceFilterType::Pointer distanceFilter = DistanceFilterType::New();
248 distanceFilter->SetInput( binaryImage );
249 distanceFilter->Update();
250 metric->SetDistanceMap( distanceFilter->GetOutput() );
251 // Software Guide : EndCodeSnippet
252
253 try
254 {
255 registration->Update();
256 }
257 catch( itk::ExceptionObject & e )
258 {
259 std::cerr << e << std::endl;
260 return EXIT_FAILURE;
261 }
262
263 std::cout << "Solution = " << transform->GetParameters() << std::endl;
264
265 return EXIT_SUCCESS;
266 }
267