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