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 #include "itkTranslationTransform.h"
20 #include "itkLinearInterpolateImageFunction.h"
21 #include "itkNormalizedCorrelationPointSetToImageMetric.h"
22 #include "itkRegularStepGradientDescentOptimizer.h"
23 #include "itkPointSet.h"
24 #include "itkPointSetToImageRegistrationMethod.h"
25 #include "itkImageRegistrationMethodImageSource.h"
26 #include "itkTestingMacros.h"
27
28 #include <iostream>
29
30 /**
31 *
32 * This program tests the registration of a PointSet against an image.
33 *
34 * This test uses two 2D-Gaussians (standard deviation RegionSize/2)
35 * One is shifted by 5 pixels from the other.
36 *
37 */
38
itkPointSetToImageRegistrationTest(int,char * [])39 int itkPointSetToImageRegistrationTest( int, char* [] )
40 {
41 constexpr unsigned int ImageDimension = 2;
42
43 using PixelType = double;
44
45 using CoordinateRepresentationType = double;
46
47 using MovingImageType = itk::Image< PixelType, ImageDimension >;
48 using FixedImageType = itk::Image< PixelType, ImageDimension >;
49
50 using ImageSourceType = itk::testhelper::ImageRegistrationMethodImageSource<
51 PixelType,
52 PixelType,
53 ImageDimension >;
54
55 ImageSourceType::Pointer imageSource = ImageSourceType::New();
56
57 itk::Size< ImageDimension > size;
58 size[0] = 100;
59 size[1] = 100;
60
61 imageSource->GenerateImages( size );
62
63 // Create the two images
64 MovingImageType::ConstPointer movingImage = imageSource->GetMovingImage();
65 FixedImageType::ConstPointer fixedImage = imageSource->GetFixedImage();
66
67 // Create the point set and load it with data by sampling
68 // the fixed image.
69 //
70 using FixedPointSetType = itk::PointSet< float, ImageDimension >;
71 FixedPointSetType::Pointer fixedPointSet = FixedPointSetType::New();
72
73 constexpr unsigned int numberOfPoints = 10000;
74
75 fixedPointSet->SetPointData( FixedPointSetType::PointDataContainer::New() );
76
77 fixedPointSet->GetPoints()->Reserve( numberOfPoints );
78 fixedPointSet->GetPointData()->Reserve( numberOfPoints );
79
80 using ImageIteratorType = itk::ImageRegionConstIterator< FixedImageType >;
81
82 ImageIteratorType it( fixedImage, fixedImage->GetBufferedRegion() );
83
84 const unsigned int skip =
85 fixedImage->GetBufferedRegion().GetNumberOfPixels() / numberOfPoints;
86
87 unsigned int counter = 0;
88
89 FixedPointSetType::PointIdentifier pointId = 0;
90 FixedPointSetType::PointType point;
91
92 it.GoToBegin();
93 while( !it.IsAtEnd() )
94 {
95 if( counter == 0 )
96 {
97 fixedImage->TransformIndexToPhysicalPoint( it.GetIndex(), point );
98 fixedPointSet->SetPoint( pointId, point );
99 fixedPointSet->SetPointData( pointId, it.Get() );
100 ++pointId;
101 if( pointId == numberOfPoints )
102 {
103 break;
104 }
105 counter = skip;
106 }
107 --counter;
108 ++it;
109 }
110
111 // Set up the Metric
112 using MetricType = itk::NormalizedCorrelationPointSetToImageMetric< FixedPointSetType,
113 MovingImageType >;
114
115 using TransformBaseType = MetricType::TransformType;
116 using ParametersType = TransformBaseType::ParametersType;
117
118 MetricType::Pointer metric = MetricType::New();
119
120 // Set up the Transform
121 using TransformType = itk::TranslationTransform< CoordinateRepresentationType,
122 ImageDimension >;
123 TransformType::Pointer transform = TransformType::New();
124
125 // Set up the Interpolator
126 using InterpolatorType =
127 itk::LinearInterpolateImageFunction< MovingImageType, double >;
128 InterpolatorType::Pointer interpolator = InterpolatorType::New();
129
130 // Set up the Optimizer
131 using OptimizerType = itk::RegularStepGradientDescentOptimizer;
132 OptimizerType::Pointer optimizer = OptimizerType::New();
133
134 // Set up the Registration method
135 using RegistrationType = itk::PointSetToImageRegistrationMethod< FixedPointSetType,
136 MovingImageType >;
137
138 RegistrationType::Pointer registration = RegistrationType::New();
139
140 EXERCISE_BASIC_OBJECT_METHODS( registration, PointSetToImageRegistrationMethod,
141 ProcessObject );
142
143 using CommandIterationType = itk::CommandIterationUpdate< OptimizerType >;
144
145 // Instantiate an Observer to report the progress of the Optimization
146 CommandIterationType::Pointer iterationCommand = CommandIterationType::New();
147 iterationCommand->SetOptimizer( optimizer );
148
149 // Scale the translation components of the Transform in the Optimizer
150 OptimizerType::ScalesType scales( transform->GetNumberOfParameters() );
151 scales.Fill( 1.0 );
152
153 unsigned long numberOfIterations = 50;
154 double maximumStepLenght = 1.0; // no step will be larger than this
155 double minimumStepLenght = 0.01;
156 double gradientTolerance = 1e-6; // convergence criterion
157
158 optimizer->SetScales( scales );
159 optimizer->SetNumberOfIterations( numberOfIterations );
160 optimizer->SetMinimumStepLength( minimumStepLenght );
161 optimizer->SetMaximumStepLength( maximumStepLenght );
162 optimizer->SetGradientMagnitudeTolerance( gradientTolerance );
163 optimizer->MinimizeOn();
164
165 // Connect all the components required for the registration
166 registration->SetMetric( metric );
167 TEST_SET_GET_VALUE( metric, registration->GetMetric() );
168
169 registration->SetOptimizer( optimizer );
170 TEST_SET_GET_VALUE( optimizer, registration->GetOptimizer() );
171
172 registration->SetTransform( transform );
173 TEST_SET_GET_VALUE( transform, registration->GetTransform() );
174
175 registration->SetFixedPointSet( fixedPointSet );
176 TEST_SET_GET_VALUE( fixedPointSet, registration->GetFixedPointSet() );
177
178 registration->SetMovingImage( movingImage );
179 TEST_SET_GET_VALUE( movingImage, registration->GetMovingImage() );
180
181 registration->SetInterpolator( interpolator );
182 TEST_SET_GET_VALUE( interpolator, registration->GetInterpolator() );
183
184 // Set up transform parameters
185 ParametersType parameters( transform->GetNumberOfParameters() );
186
187 // Initialize the offset/vector part
188 for( unsigned int k = 0; k < parameters.size(); k++ )
189 {
190 parameters[k] = 0.0f;
191 }
192 transform->SetParameters( parameters );
193 registration->SetInitialTransformParameters( transform->GetParameters() );
194 TEST_SET_GET_VALUE( transform->GetParameters(),
195 registration->GetInitialTransformParameters() );
196
197
198 TRY_EXPECT_NO_EXCEPTION( registration->Update() );
199
200
201 // Print the last transform parameters to improve coverage
202 //
203 ParametersType finalParameters = registration->GetLastTransformParameters();
204
205 const unsigned int numberOfParameters = parameters.Size();
206
207 std::cout << "Last Transform Parameters: " << std::endl;
208 for( unsigned int i = 0; i < numberOfParameters; ++i )
209 {
210 std::cout << finalParameters[i] << std::endl;
211 }
212
213 return EXIT_SUCCESS;
214 }
215