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 : BeginCommandLineArgs
20 // INPUTS: RatLungSlice1.mha
21 // INPUTS: RatLungSlice2.mha
22 // ARGUMENTS: DeformableRegistration2Output.mha
23 // ARGUMENTS: DeformableRegistration2Field.mha
24 // Software Guide : EndCommandLineArgs
25
26
27 #include "itkImageFileReader.h"
28 #include "itkImageFileWriter.h"
29
30 #include "itkImageRegionIterator.h"
31
32 // Software Guide : BeginLatex
33 //
34 // This example demonstrates how to use the ``demons'' algorithm to deformably
35 // register two images. The first step is to include the header files.
36 //
37 // Software Guide : EndLatex
38
39 // Software Guide : BeginCodeSnippet
40 #include "itkDemonsRegistrationFilter.h"
41 #include "itkHistogramMatchingImageFilter.h"
42 #include "itkCastImageFilter.h"
43 #include "itkResampleImageFilter.h"
44 #include "itkDisplacementFieldTransform.h"
45 // Software Guide : EndCodeSnippet
46
47
48 // The following section of code implements a Command observer
49 // that will monitor the evolution of the registration process.
50 //
51 class CommandIterationUpdate : public itk::Command
52 {
53 public:
54 using Self = CommandIterationUpdate;
55 using Superclass = itk::Command;
56 using Pointer = itk::SmartPointer<CommandIterationUpdate>;
57 itkNewMacro( CommandIterationUpdate );
58 protected:
59 CommandIterationUpdate() = default;
60
61 using InternalImageType = itk::Image< float, 2 >;
62 using VectorPixelType = itk::Vector< float, 2 >;
63 using DisplacementFieldType = itk::Image< VectorPixelType, 2 >;
64
65 using RegistrationFilterType = itk::DemonsRegistrationFilter<
66 InternalImageType,
67 InternalImageType,
68 DisplacementFieldType>;
69
70 public:
71
Execute(itk::Object * caller,const itk::EventObject & event)72 void Execute(itk::Object *caller, const itk::EventObject & event) override
73 {
74 Execute( (const itk::Object *)caller, event);
75 }
76
Execute(const itk::Object * object,const itk::EventObject & event)77 void Execute(const itk::Object * object, const itk::EventObject & event) override
78 {
79 const auto * filter = static_cast< const RegistrationFilterType * >( object );
80 if( !(itk::IterationEvent().CheckEvent( &event )) )
81 {
82 return;
83 }
84 std::cout << filter->GetMetric() << std::endl;
85 }
86 };
87
88
main(int argc,char * argv[])89 int main( int argc, char *argv[] )
90 {
91 if( argc < 4 )
92 {
93 std::cerr << "Missing Parameters " << std::endl;
94 std::cerr << "Usage: " << argv[0];
95 std::cerr << " fixedImageFile movingImageFile ";
96 std::cerr << " outputImageFile " << std::endl;
97 std::cerr << " [outputDisplacementFieldFile] " << std::endl;
98 return EXIT_FAILURE;
99 }
100
101 // Software Guide : BeginLatex
102 //
103 // Second, we declare the types of the images.
104 //
105 // Software Guide : EndLatex
106
107 // Software Guide : BeginCodeSnippet
108 constexpr unsigned int Dimension = 2;
109 using PixelType = unsigned short;
110
111 using FixedImageType = itk::Image< PixelType, Dimension >;
112 using MovingImageType = itk::Image< PixelType, Dimension >;
113 // Software Guide : EndCodeSnippet
114
115 // Set up the file readers
116 using FixedImageReaderType = itk::ImageFileReader< FixedImageType >;
117 using MovingImageReaderType = itk::ImageFileReader< MovingImageType >;
118
119 FixedImageReaderType::Pointer fixedImageReader = FixedImageReaderType::New();
120 MovingImageReaderType::Pointer movingImageReader = MovingImageReaderType::New();
121
122 fixedImageReader->SetFileName( argv[1] );
123 movingImageReader->SetFileName( argv[2] );
124
125
126 // Software Guide : BeginLatex
127 //
128 // Image file readers are set up in a similar fashion to previous examples.
129 // To support the re-mapping of the moving image intensity, we declare an
130 // internal image type with a floating point pixel type and cast the input
131 // images to the internal image type.
132 //
133 // Software Guide : EndLatex
134
135 // Software Guide : BeginCodeSnippet
136 using InternalPixelType = float;
137 using InternalImageType = itk::Image< InternalPixelType, Dimension >;
138 using FixedImageCasterType = itk::CastImageFilter< FixedImageType,
139 InternalImageType >;
140 using MovingImageCasterType = itk::CastImageFilter< MovingImageType,
141 InternalImageType >;
142
143 FixedImageCasterType::Pointer fixedImageCaster = FixedImageCasterType::New();
144 MovingImageCasterType::Pointer movingImageCaster
145 = MovingImageCasterType::New();
146
147 fixedImageCaster->SetInput( fixedImageReader->GetOutput() );
148 movingImageCaster->SetInput( movingImageReader->GetOutput() );
149 // Software Guide : EndCodeSnippet
150
151 // Software Guide : BeginLatex
152 //
153 // The demons algorithm relies on the assumption that pixels representing the
154 // same homologous point on an object have the same intensity on both the
155 // fixed and moving images to be registered. In this example, we will
156 // preprocess the moving image to match the intensity between the images
157 // using the \doxygen{HistogramMatchingImageFilter}.
158 //
159 // \index{itk::HistogramMatchingImageFilter}
160 //
161 // The basic idea is to match the histograms of the two images at a
162 // user-specified number of quantile values. For robustness, the histograms
163 // are matched so that the background pixels are excluded from both
164 // histograms. For MR images, a simple procedure is to exclude all gray
165 // values that are smaller than the mean gray value of the image.
166 //
167 // Software Guide : EndLatex
168
169 // Software Guide : BeginCodeSnippet
170 using MatchingFilterType = itk::HistogramMatchingImageFilter<
171 InternalImageType,
172 InternalImageType >;
173 MatchingFilterType::Pointer matcher = MatchingFilterType::New();
174 // Software Guide : EndCodeSnippet
175
176
177 // Software Guide : BeginLatex
178 //
179 // For this example, we set the moving image as the source or input image and
180 // the fixed image as the reference image.
181 //
182 // \index{itk::HistogramMatchingImageFilter!SetInput()}
183 // \index{itk::HistogramMatchingImageFilter!SetSourceImage()}
184 // \index{itk::HistogramMatchingImageFilter!SetReferenceImage()}
185 //
186 // Software Guide : EndLatex
187
188 // Software Guide : BeginCodeSnippet
189 matcher->SetInput( movingImageCaster->GetOutput() );
190 matcher->SetReferenceImage( fixedImageCaster->GetOutput() );
191 // Software Guide : EndCodeSnippet
192
193
194 // Software Guide : BeginLatex
195 //
196 // We then select the number of bins to represent the histograms and the
197 // number of points or quantile values where the histogram is to be
198 // matched.
199 //
200 // \index{itk::HistogramMatchingImageFilter!SetNumberOfHistogramLevels()}
201 // \index{itk::HistogramMatchingImageFilter!SetNumberOfMatchPoints()}
202 //
203 // Software Guide : EndLatex
204
205 // Software Guide : BeginCodeSnippet
206 matcher->SetNumberOfHistogramLevels( 1024 );
207 matcher->SetNumberOfMatchPoints( 7 );
208 // Software Guide : EndCodeSnippet
209
210
211 // Software Guide : BeginLatex
212 //
213 // Simple background extraction is done by thresholding at the mean
214 // intensity.
215 //
216 // \index{itk::HistogramMatchingImageFilter!ThresholdAtMeanIntensityOn()}
217 //
218 // Software Guide : EndLatex
219
220 // Software Guide : BeginCodeSnippet
221 matcher->ThresholdAtMeanIntensityOn();
222 // Software Guide : EndCodeSnippet
223
224
225 // Software Guide : BeginLatex
226 //
227 // In the \doxygen{DemonsRegistrationFilter}, the deformation field is
228 // represented as an image whose pixels are floating point vectors.
229 //
230 // \index{itk::DemonsRegistrationFilter}
231 //
232 // Software Guide : EndLatex
233
234 // Software Guide : BeginCodeSnippet
235 using VectorPixelType = itk::Vector< float, Dimension >;
236 using DisplacementFieldType = itk::Image< VectorPixelType, Dimension >;
237 using RegistrationFilterType = itk::DemonsRegistrationFilter<
238 InternalImageType,
239 InternalImageType,
240 DisplacementFieldType>;
241 RegistrationFilterType::Pointer filter = RegistrationFilterType::New();
242 // Software Guide : EndCodeSnippet
243
244
245 // Create the Command observer and register it with the registration filter.
246 //
247 CommandIterationUpdate::Pointer observer = CommandIterationUpdate::New();
248 filter->AddObserver( itk::IterationEvent(), observer );
249
250
251 // Software Guide : BeginLatex
252 //
253 // The input fixed image is simply the output of the fixed image casting
254 // filter. The input moving image is the output of the histogram matching
255 // filter.
256 //
257 // \index{itk::DemonsRegistrationFilter!SetFixedImage()}
258 // \index{itk::DemonsRegistrationFilter!SetMovingImage()}
259 //
260 // Software Guide : EndLatex
261
262 // Software Guide : BeginCodeSnippet
263 filter->SetFixedImage( fixedImageCaster->GetOutput() );
264 filter->SetMovingImage( matcher->GetOutput() );
265 // Software Guide : EndCodeSnippet
266
267
268 // Software Guide : BeginLatex
269 //
270 // The demons registration filter has two parameters: the number of
271 // iterations to be performed and the standard deviation of the Gaussian
272 // smoothing kernel to be applied to the deformation field after each
273 // iteration.
274 // \index{itk::DemonsRegistrationFilter!SetNumberOfIterations()}
275 // \index{itk::DemonsRegistrationFilter!SetStandardDeviations()}
276 //
277 // Software Guide : EndLatex
278
279 // Software Guide : BeginCodeSnippet
280 filter->SetNumberOfIterations( 50 );
281 filter->SetStandardDeviations( 1.0 );
282 // Software Guide : EndCodeSnippet
283
284
285 // Software Guide : BeginLatex
286 //
287 // The registration algorithm is triggered by updating the filter. The
288 // filter output is the computed deformation field.
289 //
290 // Software Guide : EndLatex
291
292 // Software Guide : BeginCodeSnippet
293 filter->Update();
294 // Software Guide : EndCodeSnippet
295
296
297 // Software Guide : BeginLatex
298 //
299 // The \doxygen{ResampleImageFilter} can be used to warp the moving image with
300 // the output deformation field. The default interpolator of the \doxygen{ResampleImageFilter},
301 // is used but specification of the output image spacing and origin
302 // are required.
303 //
304 // \index{itk::ResampleImageFilter}
305 // \index{itk::ResampleImageFilter!SetInput()}
306 // \index{itk::ResampleImageFilter!SetInterpolator()}
307 // \index{itk::ResampleImageFilter!SetOutputSpacing()}
308 // \index{itk::ResampleImageFilter!SetOutputOrigin()}
309 //
310 // Software Guide : EndLatex
311
312 // Software Guide : BeginCodeSnippet
313 using OutputPixelType = unsigned char;
314 using OutputImageType = itk::Image< OutputPixelType, Dimension >;
315
316 using InterpolatorPrecisionType = double;
317 using WarperType =
318 itk::ResampleImageFilter< MovingImageType, OutputImageType,
319 InterpolatorPrecisionType, float >;
320 WarperType::Pointer warper = WarperType::New();
321
322 warper->SetInput( movingImageReader->GetOutput() );
323 warper->UseReferenceImageOn();
324 warper->SetReferenceImage( fixedImageReader->GetOutput() );
325
326 // Software Guide : EndCodeSnippet
327
328
329 // Software Guide : BeginLatex
330 //
331 // The \code{ResampleImageFilter} requires a transform, so a
332 // \doxygen{DisplacementFieldTransform} must be constructed then set
333 // as the transform of the \code{ResampleImageFilter}. The resulting
334 // warped or resampled image is written to file as per previous
335 // examples.
336 //
337 // Software Guide : EndLatex
338
339 // Software Guide : BeginCodeSnippet
340
341 using DisplacementFieldTransformType =
342 itk::DisplacementFieldTransform<InternalPixelType, Dimension>;
343 auto displacementTransform = DisplacementFieldTransformType::New();
344 displacementTransform->SetDisplacementField( filter->GetOutput() );
345
346 warper->SetTransform( displacementTransform );
347 // Software Guide : EndCodeSnippet
348
349
350 // Write warped image out to file
351 using WriterType = itk::ImageFileWriter< OutputImageType >;
352
353 WriterType::Pointer writer = WriterType::New();
354
355 writer->SetFileName( argv[3] );
356 writer->SetInput( warper->GetOutput() );
357 writer->Update();
358
359
360 // Software Guide : BeginLatex
361 //
362 // Let's execute this example using the rat lung data from the previous example.
363 // The associated data files can be found in \code{Examples/Data}:
364 //
365 // \begin{itemize}
366 // \item \code{RatLungSlice1.mha}
367 // \item \code{RatLungSlice2.mha}
368 // \end{itemize}
369 //
370 // \begin{figure} \center
371 // \includegraphics[width=0.44\textwidth]{DeformableRegistration2CheckerboardBefore}
372 // \includegraphics[width=0.44\textwidth]{DeformableRegistration2CheckerboardAfter}
373 // \itkcaption[Demon's deformable registration output]{Checkerboard comparisons
374 // before and after demons-based deformable registration.}
375 // \label{fig:DeformableRegistration2Output}
376 // \end{figure}
377 //
378 // The result of the demons-based deformable registration is presented in
379 // Figure \ref{fig:DeformableRegistration2Output}. The checkerboard
380 // comparison shows that the algorithm was able to recover the misalignment
381 // due to expiration.
382 //
383 // Software Guide : EndLatex
384
385
386 // Software Guide : BeginLatex
387 //
388 // It may be also desirable to write the deformation field as an image of
389 // vectors. This can be done with the following code.
390 //
391 // Software Guide : EndLatex
392
393 if( argc > 4 ) // if a fourth line argument has been provided...
394 {
395
396 // Software Guide : BeginCodeSnippet
397 using FieldWriterType = itk::ImageFileWriter< DisplacementFieldType >;
398 FieldWriterType::Pointer fieldWriter = FieldWriterType::New();
399 fieldWriter->SetFileName( argv[4] );
400 fieldWriter->SetInput( filter->GetOutput() );
401
402 fieldWriter->Update();
403 // Software Guide : EndCodeSnippet
404
405 // Software Guide : BeginLatex
406 //
407 // Note that the file format used for writing the deformation field must be
408 // capable of representing multiple components per pixel. This is the case
409 // for the MetaImage and VTK file formats for example.
410 //
411 // Software Guide : EndLatex
412
413 }
414
415
416 if( argc > 5 ) // if a fifth line argument has been provided...
417 {
418
419 using VectorImage2DType = DisplacementFieldType;
420 using Vector2DType = DisplacementFieldType::PixelType;
421
422 VectorImage2DType::ConstPointer vectorImage2D = filter->GetOutput();
423
424 VectorImage2DType::RegionType region2D = vectorImage2D->GetBufferedRegion();
425 VectorImage2DType::IndexType index2D = region2D.GetIndex();
426 VectorImage2DType::SizeType size2D = region2D.GetSize();
427
428
429 using Vector3DType = itk::Vector< float, 3 >;
430 using VectorImage3DType = itk::Image< Vector3DType, 3 >;
431
432 using VectorImage3DWriterType = itk::ImageFileWriter< VectorImage3DType >;
433
434 VectorImage3DWriterType::Pointer writer3D = VectorImage3DWriterType::New();
435
436 VectorImage3DType::Pointer vectorImage3D = VectorImage3DType::New();
437
438 VectorImage3DType::RegionType region3D;
439 VectorImage3DType::IndexType index3D;
440 VectorImage3DType::SizeType size3D;
441
442 index3D[0] = index2D[0];
443 index3D[1] = index2D[1];
444 index3D[2] = 0;
445
446 size3D[0] = size2D[0];
447 size3D[1] = size2D[1];
448 size3D[2] = 1;
449
450 region3D.SetSize( size3D );
451 region3D.SetIndex( index3D );
452
453 vectorImage3D->SetRegions( region3D );
454 vectorImage3D->Allocate();
455
456 using Iterator2DType = itk::ImageRegionConstIterator< VectorImage2DType >;
457
458 using Iterator3DType = itk::ImageRegionIterator< VectorImage3DType >;
459
460 Iterator2DType it2( vectorImage2D, region2D );
461 Iterator3DType it3( vectorImage3D, region3D );
462
463 it2.GoToBegin();
464 it3.GoToBegin();
465
466 Vector2DType vector2D;
467 Vector3DType vector3D;
468
469 vector3D[2] = 0; // set Z component to zero.
470
471 while( !it2.IsAtEnd() )
472 {
473 vector2D = it2.Get();
474 vector3D[0] = vector2D[0];
475 vector3D[1] = vector2D[1];
476 it3.Set( vector3D );
477 ++it2;
478 ++it3;
479 }
480
481
482 writer3D->SetInput( vectorImage3D );
483
484 writer3D->SetFileName( argv[5] );
485
486 try
487 {
488 writer3D->Update();
489 }
490 catch( itk::ExceptionObject & excp )
491 {
492 std::cerr << excp << std::endl;
493 return EXIT_FAILURE;
494 }
495
496
497 }
498
499 return EXIT_SUCCESS;
500 }
501