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 "itkFEMRegistrationFilter.h"
20 #include "itkImageFileWriter.h"
21 #include "itkTestingMacros.h"
22
23 #include <fstream>
24
25
26 // Typedefs used for registration
27 constexpr unsigned int ImageDimension = 2;
28
29 using InputImagePixelType = unsigned char;
30 using DeformationFieldPixelType = float;
31
32 using InputImageType = itk::Image< InputImagePixelType, ImageDimension >;
33 using DeformationFieldVectorType = itk::Vector< DeformationFieldPixelType, ImageDimension >;
34 using DeformationFieldImageType = itk::Image< DeformationFieldVectorType, ImageDimension >;
35
36 using ElementType = itk::fem::Element2DC0LinearQuadrilateralMembrane;
37
38 // Template function to fill in an image with a value.
39 template< typename TImage >
FillImage(TImage * image,typename TImage::PixelType value)40 void FillImage( TImage * image, typename TImage::PixelType value )
41 {
42 using Iterator = itk::ImageRegionIteratorWithIndex< TImage >;
43 Iterator it( image, image->GetBufferedRegion() );
44
45 for( it.GoToBegin(); !it.IsAtEnd(); ++it )
46 {
47 it.Set( value );
48 }
49 }
50
51 // Template function to fill in an image with a circle.
52 template< typename TImage >
FillWithCircle(TImage * image,double * center,double radius,typename TImage::PixelType foregnd,typename TImage::PixelType backgnd)53 void FillWithCircle( TImage * image, double * center, double radius,
54 typename TImage::PixelType foregnd, typename TImage::PixelType backgnd )
55 {
56 using Iterator = itk::ImageRegionIteratorWithIndex< TImage >;
57 Iterator it( image, image->GetBufferedRegion() );
58
59 typename TImage::IndexType index;
60 double r2 = itk::Math::sqr( radius );
61 for( it.GoToBegin(); !it.IsAtEnd(); ++it )
62 {
63 index = it.GetIndex();
64 double distance = 0;
65 for( unsigned int j = 0; j < TImage::ImageDimension; ++j )
66 {
67 distance += itk::Math::sqr( (double) index[j] - center[j] );
68 }
69 if( distance <= r2 )
70 {
71 it.Set( foregnd );
72 }
73 else
74 {
75 it.Set( backgnd );
76 }
77 }
78 }
79
80
itkFEMRegistrationFilter2DTest(int argc,char * argv[])81 int itkFEMRegistrationFilter2DTest( int argc, char *argv[] )
82 {
83 using IndexType = InputImageType::IndexType;
84 using SizeType = InputImageType::SizeType;
85 using RegionType = InputImageType::RegionType;
86
87
88 // Generate input images and initial deformation field
89
90 InputImageType::SizeValueType sizeArray[ImageDimension];
91 for(auto & i : sizeArray)
92 {
93 i = 32;
94 }
95
96 SizeType size;
97 size.SetSize( sizeArray );
98
99 IndexType index;
100 index.Fill( 0 );
101
102 RegionType region;
103 region.SetSize( size );
104 region.SetIndex( index );
105
106 InputImageType::Pointer movingImage = InputImageType::New();
107 InputImageType::Pointer fixedImage = InputImageType::New();
108
109 DeformationFieldImageType::Pointer initField = DeformationFieldImageType::New();
110
111 movingImage->SetLargestPossibleRegion( region );
112 movingImage->SetBufferedRegion( region );
113 movingImage->Allocate();
114
115 fixedImage->SetLargestPossibleRegion( region );
116 fixedImage->SetBufferedRegion( region );
117 fixedImage->Allocate();
118
119 initField->SetLargestPossibleRegion( region );
120 initField->SetBufferedRegion( region );
121 initField->Allocate();
122
123 double center[ImageDimension];
124 double radius;
125 InputImagePixelType fgnd = 250;
126 InputImagePixelType bgnd = 15;
127
128 // Set the circle center
129 for(double & i : center)
130 {
131 i = 16;
132 }
133
134 // Fill fixed image with a circle
135 radius = 8;
136 FillWithCircle< InputImageType >( fixedImage, center, radius, fgnd, bgnd );
137
138 // Fill moving image with a circle
139 radius = 5;
140 FillWithCircle< InputImageType >( movingImage, center, radius, fgnd, bgnd );
141
142 // Fill initial deformation with zero vectors
143 DeformationFieldVectorType zeroVec;
144 zeroVec.Fill( 0.0 );
145 FillImage< DeformationFieldImageType >( initField, zeroVec );
146
147
148 using FEMObjectType = itk::fem::FEMObject< ImageDimension >;
149 using RegistrationType = itk::fem::FEMRegistrationFilter< InputImageType,
150 InputImageType,
151 FEMObjectType >;
152
153 // Run registration and warp moving
154 for( unsigned int met = 0; met < 4; ++met )
155 {
156 RegistrationType::Pointer registrator = RegistrationType::New();
157
158 EXERCISE_BASIC_OBJECT_METHODS( registrator, FEMRegistrationFilter,
159 ImageToImageFilter );
160
161 registrator->SetFixedImage( fixedImage );
162 TEST_SET_GET_VALUE( fixedImage, registrator->GetFixedImage() );
163
164 registrator->SetMovingImage( movingImage );
165 TEST_SET_GET_VALUE( movingImage, registrator->GetMovingImage() );
166
167 unsigned int maxLevel = 1;
168 registrator->SetMaxLevel( maxLevel );
169 TEST_SET_GET_VALUE( maxLevel, registrator->GetMaxLevel() );
170
171 bool useNormalizedGradient = true;
172 TEST_SET_GET_BOOLEAN( registrator, UseNormalizedGradient, useNormalizedGradient );
173
174 registrator->ChooseMetric( met );
175
176 unsigned int numberOfMaxIterations = 5;
177 registrator->SetMaximumIterations( numberOfMaxIterations, 0 );
178
179 RegistrationType::Float elasticity = 10;
180 registrator->SetElasticity( elasticity, 0 );
181 TEST_SET_GET_VALUE( elasticity, registrator->GetElasticity() );
182
183 RegistrationType::Float rho = 1;
184 registrator->SetRho( rho, 0 );
185 //TEST_SET_GET_VALUE( rho, registrator->GetRho() );
186
187 RegistrationType::Float gamma = 1.;
188 registrator->SetGamma( gamma, 0 );
189 //TEST_SET_GET_VALUE( gamma, registrator->GetGamma() );
190
191 RegistrationType::Float alpha = 1.;
192 registrator->SetAlpha( alpha );
193 TEST_SET_GET_VALUE( alpha, registrator->GetAlpha() );
194
195 registrator->SetMeshPixelsPerElementAtEachResolution( 4, 0 );
196
197 unsigned int widthOfMetricRegion;
198 if( met == 0 || met == 3 )
199 {
200 widthOfMetricRegion = 0;
201 }
202 else
203 {
204 widthOfMetricRegion = 1;
205 }
206 registrator->SetWidthOfMetricRegion( widthOfMetricRegion, 0 );
207 TEST_SET_GET_VALUE( widthOfMetricRegion,
208 registrator->GetWidthOfMetricRegion() );
209
210 registrator->SetNumberOfIntegrationPoints( 2, 0 );
211
212 RegistrationType::Float timeStep = 1.;
213 registrator->SetTimeStep( timeStep );
214 TEST_SET_GET_VALUE( timeStep, registrator->GetTimeStep() );
215
216 unsigned int doLineSearchOnImageEnergy;
217 unsigned int employRegridding;
218 if( met == 0 )
219 {
220 doLineSearchOnImageEnergy = 2;
221 employRegridding = true;
222 }
223 else
224 {
225 doLineSearchOnImageEnergy = 0;
226 employRegridding = false;
227 }
228
229 registrator->SetDoLineSearchOnImageEnergy( doLineSearchOnImageEnergy );
230 TEST_SET_GET_VALUE( doLineSearchOnImageEnergy,
231 registrator->GetDoLineSearchOnImageEnergy() );
232
233 registrator->SetEmployRegridding( employRegridding );
234 TEST_SET_GET_VALUE( employRegridding, registrator->GetEmployRegridding() );
235
236 bool useLandmarks = false;
237 TEST_SET_GET_BOOLEAN( registrator, UseLandmarks, useLandmarks );
238
239 bool useMassMatrix = true;
240 TEST_SET_GET_BOOLEAN( registrator, UseMassMatrix, useMassMatrix );
241
242 RegistrationType::Float energyReductionFactor = 0.0;
243 registrator->SetEnergyReductionFactor( energyReductionFactor );
244 TEST_SET_GET_VALUE( energyReductionFactor,
245 registrator->GetEnergyReductionFactor() );
246
247 unsigned int lineSearchMaximumIterations = 100;
248 registrator->SetLineSearchMaximumIterations( lineSearchMaximumIterations );
249 TEST_SET_GET_VALUE( lineSearchMaximumIterations,
250 registrator->GetLineSearchMaximumIterations() );
251
252 bool createMeshFromImage = true;
253 TEST_SET_GET_BOOLEAN( registrator, CreateMeshFromImage, createMeshFromImage );
254
255 double standardDeviation = 0.5;
256 registrator->SetStandardDeviations( standardDeviation );
257 //TEST_SET_GET_VALUE( standardDeviations, registrator->GetStandardDeviations() );
258
259 standardDeviation = 1.0;
260 RegistrationType::StandardDeviationsType standardDeviations;
261 standardDeviations.Fill( standardDeviation );
262 registrator->SetStandardDeviations( standardDeviations );
263 //TEST_SET_GET_VALUE( standardDeviations, registrator->GetStandardDeviations() );
264
265 unsigned int maximumKernelWidth = 30;
266 registrator->SetMaximumKernelWidth( maximumKernelWidth );
267 TEST_SET_GET_VALUE( maximumKernelWidth, registrator->GetMaximumKernelWidth() );
268
269 double maximumError = 0.1;
270 registrator->SetMaximumError( maximumError );
271 TEST_SET_GET_VALUE( maximumError, registrator->GetMaximumError() );
272
273
274 itk::fem::MaterialLinearElasticity::Pointer material =
275 itk::fem::MaterialLinearElasticity::New();
276 material->SetGlobalNumber( 0 );
277 material->SetYoungsModulus( registrator->GetElasticity() );
278 material->SetCrossSectionalArea( 1.0 );
279 material->SetThickness( 1.0 );
280 material->SetMomentOfInertia( 1.0 );
281 material->SetPoissonsRatio( 0. ); // DON'T CHOOSE 1.0!!
282 material->SetDensityHeatProduct( 1.0 );
283
284 // Create the element type
285 ElementType::Pointer element1 = ElementType::New();
286 element1->SetMaterial( dynamic_cast< itk::fem::MaterialLinearElasticity * >( &*material ) );
287 registrator->SetElement( &*element1 );
288 registrator->SetMaterial( material );
289
290 // Register the images
291 try
292 {
293 registrator->RunRegistration();
294 }
295 catch( ::itk::ExceptionObject & err )
296 {
297 std::cerr << "ITK exception detected: " << err;
298 std::cout << "Test failed!" << std::endl;
299 return EXIT_FAILURE;
300 }
301 catch( ... )
302 {
303 // fixme - changes to femparray cause it to fail : old version works
304 std::cout << "Caught an exception: " << std::endl;
305 return EXIT_FAILURE;
306 // std::cout << err << std::endl;
307 // throw err;
308 }
309
310 if( argc == 2 )
311 {
312 std::string outFileName = argv[1];
313 std::stringstream ss;
314 ss << met;
315 outFileName += ss.str();
316 outFileName += ".mhd";
317 using ImageWriterType = itk::ImageFileWriter< RegistrationType::FieldType >;
318 ImageWriterType::Pointer writer = ImageWriterType::New();
319 writer->SetFileName( outFileName );
320 writer->SetInput( registrator->GetDisplacementField() );
321 writer->Update();
322 }
323
324 if( argc == 3 )
325 {
326 std::string outFileName = argv[2];
327 std::stringstream ss;
328 ss << met;
329 outFileName += ss.str();
330 outFileName += ".mhd";
331 using ImageWriterType = itk::ImageFileWriter< InputImageType >;
332 ImageWriterType::Pointer writer = ImageWriterType::New();
333 writer->SetFileName( outFileName );
334 writer->SetInput( registrator->GetWarpedImage() );
335 writer->Update();
336 }
337 }
338
339 /*
340 // get warped reference image
341 // ---------------------------------------------------------
342 std::cout << "Compare warped moving and fixed." << std::endl;
343
344 // compare the warp and fixed images
345 itk::ImageRegionIterator<ImageType> fixedIter( fixed,
346 fixed->GetBufferedRegion() );
347 itk::ImageRegionIterator<ImageType> warpedIter( registrator->GetWarpedImage(),
348 fixed->GetBufferedRegion() );
349
350 unsigned int numPixelsDifferent = 0;
351 while( !fixedIter.IsAtEnd() )
352 {
353 if( fixedIter.Get() != warpedIter.Get() )
354 {
355 numPixelsDifferent++;
356 }
357 ++fixedIter;
358 ++warpedIter;
359 }
360
361 std::cout << "Number of pixels different: " << numPixelsDifferent;
362 std::cout << std::endl;
363
364 if( numPixelsDifferent > 400 )
365 {
366 std::cout << "Test failed - too many pixels different." << std::endl;
367 return EXIT_FAILURE;
368 }
369 std::cout << "Test passed" << std::endl;
370 */
371
372 std::cout << "Test finished." << std::endl;
373 return EXIT_SUCCESS;
374 }
375