forked from InsightSoftwareConsortium/ITK
-
Notifications
You must be signed in to change notification settings - Fork 0
/
ModelToImageRegistration1.cxx
935 lines (780 loc) · 31.8 KB
/
ModelToImageRegistration1.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
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
/*=========================================================================
*
* 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.
*
*=========================================================================*/
// Software Guide : BeginLatex
//
// This example illustrates the use of the \doxygen{SpatialObject} as a
// component of the registration framework in order to perform model based
// registration. The current example creates a geometrical model composed of
// several ellipses. Then, it uses the model to produce a synthetic binary
// image of the ellipses. Next, it introduces perturbations on the position
// and shape of the model, and finally it uses the perturbed version as the
// input to a registration problem. A metric is defined to evaluate the
// fitness between the geometric model and the image.
//
// Software Guide : EndLatex
// Software Guide : BeginLatex
//
// Let's look first at the classes required to support
// SpatialObject. In this example we use the
// \doxygen{EllipseSpatialObject} as the basic shape components and we use the
// \doxygen{GroupSpatialObject} to group them together as a representation of
// a more complex shape. Their respective headers are included below.
//
// \index{itk::EllipseSpatialObject!header}
// \index{itk::GroupSpatialObject!header}
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
#include "itkEllipseSpatialObject.h"
#include "itkGroupSpatialObject.h"
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// In order to generate the initial synthetic image of the ellipses, we use
// the \doxygen{SpatialObjectToImageFilter} that tests---for every pixel in
// the image---whether the pixel (and hence the spatial object) is
// \emph{inside} or \emph{outside} the geometric model.
//
// \index{itk::Spatial\-Object\-To\-Image\-Filter!header}
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
#include "itkSpatialObjectToImageFilter.h"
// Software Guide : EndCodeSnippet
#include "itkImageToSpatialObjectRegistrationMethod.h"
// Software Guide : BeginLatex
//
// A metric is defined to evaluate the fitness between the
// SpatialObject and the Image. The base class for this
// type of metric is the \doxygen{ImageToSpatialObjectMetric}, whose header is
// included below.
//
// \index{itk::Image\-To\-Spatial\-Object\-Metric!header}
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// As in previous registration problems, we have to evaluate the image
// intensity in non-grid positions. The
// \doxygen{LinearInterpolateImageFunction} is used here for this purpose.
//
// \index{itk::Linear\-Interpolate\-Image\-Function!header}
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
#include "itkLinearInterpolateImageFunction.h"
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// The SpatialObject is mapped from its own space into the image
// space by using a \doxygen{Transform}. In this
// example, we use the \doxygen{Euler2DTransform}.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
#include "itkEuler2DTransform.h"
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// Registration is fundamentally an optimization problem. Here we include
// the optimizer used to search the parameter space and identify the best
// transformation that will map the shape model on top of the image. The
// optimizer used in this example is the
// \doxygen{OnePlusOneEvolutionaryOptimizer} that implements an
// \href{http://www.aic.nrl.navy.mil/galist/}{evolutionary algorithm}.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
#include "itkOnePlusOneEvolutionaryOptimizer.h"
// Software Guide : EndCodeSnippet
#include "itkDiscreteGaussianImageFilter.h"
#include "itkNormalVariateGenerator.h"
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include "itkCastImageFilter.h"
#include "itkRescaleIntensityImageFilter.h"
// Software Guide : BeginLatex
//
// As in previous registration examples, it is important to
// track the evolution of the optimizer as it progresses through the parameter
// space. This is done by using the Command/Observer paradigm. The
// following lines of code implement the \doxygen{Command} observer that
// monitors the progress of the registration. The code is quite
// similar to what we have used in previous registration examples.
//
// \index{Model to Image Registration!Observer}
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
#include "itkCommand.h"
template < class TOptimizer >
class IterationCallback : public itk::Command
{
public:
typedef IterationCallback Self;
typedef itk::Command Superclass;
typedef itk::SmartPointer<Self> Pointer;
typedef itk::SmartPointer<const Self> ConstPointer;
itkTypeMacro( IterationCallback, Superclass );
itkNewMacro( Self );
/** Type defining the optimizer. */
typedef TOptimizer OptimizerType;
/** Method to specify the optimizer. */
void SetOptimizer( OptimizerType * optimizer )
{
m_Optimizer = optimizer;
m_Optimizer->AddObserver( itk::IterationEvent(), this );
}
/** Execute method will print data at each iteration */
void Execute(itk::Object *caller, const itk::EventObject & event)
{
Execute( (const itk::Object *)caller, event);
}
void Execute(const itk::Object *, const itk::EventObject & event)
{
if( typeid( event ) == typeid( itk::StartEvent ) )
{
std::cout << std::endl << "Position Value";
std::cout << std::endl << std::endl;
}
else if( typeid( event ) == typeid( itk::IterationEvent ) )
{
std::cout << m_Optimizer->GetCurrentIteration() << " ";
std::cout << m_Optimizer->GetValue() << " ";
std::cout << m_Optimizer->GetCurrentPosition() << std::endl;
}
else if( typeid( event ) == typeid( itk::EndEvent ) )
{
std::cout << std::endl << std::endl;
std::cout << "After " << m_Optimizer->GetCurrentIteration();
std::cout << " iterations " << std::endl;
std::cout << "Solution is = " << m_Optimizer->GetCurrentPosition();
std::cout << std::endl;
}
}
// Software Guide : EndCodeSnippet
protected:
IterationCallback() {};
itk::WeakPointer<OptimizerType> m_Optimizer;
};
// Software Guide : BeginLatex
//
// This command will be invoked at every iteration of the optimizer and will
// print out the current combination of transform parameters.
//
// Software Guide : EndLatex
// Software Guide : BeginLatex
//
// Consider now the most critical component of this new registration
// approach: the metric. This component evaluates the match between the
// SpatialObject and the Image. The
// smoothness and regularity of the metric determine the difficulty of the
// task assigned to the optimizer. In this case, we use a very robust
// optimizer that should be able to find its way even in the most
// discontinuous cost functions. The metric to be implemented should derive
// from the ImageToSpatialObjectMetric class.
//
// The following code implements a simple metric that computes the sum of
// the pixels that are inside the spatial object. In fact, the metric
// maximum is obtained when the model and the image are aligned. The metric
// is templated over the type of the SpatialObject and the type of
// the Image.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
template <typename TFixedImage, typename TMovingSpatialObject>
class SimpleImageToSpatialObjectMetric :
public itk::ImageToSpatialObjectMetric<TFixedImage,TMovingSpatialObject>
{
// Software Guide : EndCodeSnippet
public:
/** Standard class typedefs. */
typedef SimpleImageToSpatialObjectMetric Self;
typedef itk::ImageToSpatialObjectMetric<TFixedImage,TMovingSpatialObject>
Superclass;
typedef itk::SmartPointer<Self> Pointer;
typedef itk::SmartPointer<const Self> ConstPointer;
typedef itk::Point<double,2> PointType;
typedef std::list<PointType> PointListType;
typedef TMovingSpatialObject MovingSpatialObjectType;
typedef typename Superclass::ParametersType ParametersType;
typedef typename Superclass::DerivativeType DerivativeType;
typedef typename Superclass::MeasureType MeasureType;
/** Method for creation through the object factory. */
itkNewMacro(Self);
/** Run-time type information (and related methods). */
itkTypeMacro(SimpleImageToSpatialObjectMetric, ImageToSpatialObjectMetric);
itkStaticConstMacro( ParametricSpaceDimension, unsigned int, 3 );
/** Specify the moving spatial object. */
void SetMovingSpatialObject( const MovingSpatialObjectType * object)
{
if(!this->m_FixedImage)
{
std::cout << "Please set the image before the moving spatial object" << std::endl;
return;
}
this->m_MovingSpatialObject = object;
m_PointList.clear();
typedef itk::ImageRegionConstIteratorWithIndex<TFixedImage> myIteratorType;
myIteratorType it(this->m_FixedImage,this->m_FixedImage->GetBufferedRegion());
itk::Point<double,2> point;
while( !it.IsAtEnd() )
{
this->m_FixedImage->TransformIndexToPhysicalPoint( it.GetIndex(), point );
if(this->m_MovingSpatialObject->IsInside(point,99999))
{
m_PointList.push_back( point );
}
++it;
}
std::cout << "Number of points in the metric = " << static_cast<unsigned long>( m_PointList.size() ) << std::endl;
}
/** Get the Derivatives of the Match Measure */
void GetDerivative( const ParametersType &, DerivativeType & ) const
{
return;
}
// Software Guide : BeginLatex
//
// The fundamental operation of the metric is its \code{GetValue()} method.
// It is in this method that the fitness value is computed. In our current
// example, the fitness is computed over the points of the
// SpatialObject. For each point, its coordinates are mapped
// through the transform into image space. The resulting point is used
// to evaluate the image and the resulting value is accumulated in a sum.
// Since we are not allowing scale changes, the optimal value of the sum
// will result when all the SpatialObject points are mapped on
// the white regions of the image. Note that the argument for the
// \code{GetValue()} method is the array of parameters of the transform.
//
// \index{Image\-To\-Spatial\-Object\-Metric!GetValue()}
//
// Software Guide : EndLatex
/** Get the value for SingleValue optimizers. */
// Software Guide : BeginCodeSnippet
MeasureType GetValue( const ParametersType & parameters ) const
{
double value;
this->m_Transform->SetParameters( parameters );
value = 0;
for(PointListType::const_iterator it = m_PointList.begin();
it != m_PointList.end(); ++it)
{
PointType transformedPoint = this->m_Transform->TransformPoint(*it);
if( this->m_Interpolator->IsInsideBuffer( transformedPoint ) )
{
value += this->m_Interpolator->Evaluate( transformedPoint );
}
}
return value;
}
// Software Guide : EndCodeSnippet
/** Get Value and Derivatives for MultipleValuedOptimizers */
void GetValueAndDerivative( const ParametersType & parameters,
MeasureType & Value, DerivativeType & Derivative ) const
{
Value = this->GetValue(parameters);
this->GetDerivative(parameters,Derivative);
}
private:
PointListType m_PointList;
};
// Software Guide : BeginLatex
//
// Having defined all the registration components we are ready to put the
// pieces together and implement the registration process.
//
// Software Guide : EndLatex
int main( int argc, char *argv[] )
{
if( argc > 1 )
{
std::cerr << "Too many parameters " << std::endl;
std::cerr << "Usage: " << argv[0] << std::endl;
}
// Software Guide : BeginLatex
//
// First we instantiate the GroupSpatialObject and
// EllipseSpatialObject. These two objects are parameterized by
// the dimension of the space. In our current example a $2D$ instantiation
// is created.
//
// \index{Group\-Spatial\-Object!Instantiation}
// \index{Ellipse\-Spatial\-Object!Instantiation}
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
typedef itk::GroupSpatialObject< 2 > GroupType;
typedef itk::EllipseSpatialObject< 2 > EllipseType;
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// The image is instantiated in the following lines using the pixel
// type and the space dimension. This image uses a \code{float} pixel
// type since we plan to blur it in order to increase the capture radius of
// the optimizer. Images of real pixel type behave better under blurring
// than those of integer pixel type.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
typedef itk::Image< float, 2 > ImageType;
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// Here is where the fun begins! In the following lines we create the
// EllipseSpatialObjects using their \code{New()} methods, and
// assigning the results to SmartPointers. These lines will create
// three ellipses.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
EllipseType::Pointer ellipse1 = EllipseType::New();
EllipseType::Pointer ellipse2 = EllipseType::New();
EllipseType::Pointer ellipse3 = EllipseType::New();
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// Every class deriving from SpatialObject has particular
// parameters enabling the user to tailor its shape. In the case of the
// EllipseSpatialObject, \code{SetRadius()} is used to
// define the ellipse size. An additional \code{SetRadius(Array)} method
// allows the user to define the ellipse axes independently.
//
// \index{itk::EllipseSpatialObject!SetRadius()}
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
ellipse1->SetRadius( 10.0 );
ellipse2->SetRadius( 10.0 );
ellipse3->SetRadius( 10.0 );
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// The ellipses are created centered in space by default. We use the
// following lines of code to arrange the ellipses in a triangle.
// The spatial transform intrinsically associated with the object is
// accessed by the \code{GetTransform()} method. This transform can define
// a translation in space with the \code{SetOffset()} method. We take
// advantage of this feature to place the ellipses at particular
// points in space.
//
// Software Guide : EndLatex
// Place each ellipse at the right position to form a triangle
// Software Guide : BeginCodeSnippet
EllipseType::TransformType::OffsetType offset;
offset[ 0 ] = 100.0;
offset[ 1 ] = 40.0;
ellipse1->GetObjectToParentTransform()->SetOffset(offset);
ellipse1->ComputeObjectToWorldTransform();
offset[ 0 ] = 40.0;
offset[ 1 ] = 150.0;
ellipse2->GetObjectToParentTransform()->SetOffset(offset);
ellipse2->ComputeObjectToWorldTransform();
offset[ 0 ] = 150.0;
offset[ 1 ] = 150.0;
ellipse3->GetObjectToParentTransform()->SetOffset(offset);
ellipse3->ComputeObjectToWorldTransform();
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// Note that after a change has been made in the transform, the
// SpatialObject invokes the method
// \code{ComputeGlobalTransform()} in order to update its global
// transform. The reason for doing this is that SpatialObjects
// can be arranged in hierarchies. It is then possible to change the
// position of a set of spatial objects by moving the parent of the group.
//
// Software Guide : EndLatex
// Software Guide : BeginLatex
//
// Now we add the three EllipseSpatialObjects to a
// GroupSpatialObject that will be subsequently passed on to the
// registration method. The GroupSpatialObject facilitates the
// management of the three ellipses as a higher level structure
// representing a complex shape. Groups can be nested any number of levels
// in order to represent shapes with higher detail.
//
// \index{itk::GroupSpatialObject!New()}
// \index{itk::GroupSpatialObject!Pointer}
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
GroupType::Pointer group = GroupType::New();
group->AddSpatialObject( ellipse1 );
group->AddSpatialObject( ellipse2 );
group->AddSpatialObject( ellipse3 );
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// Having the geometric model ready, we proceed to generate the binary
// image representing the imprint of the space occupied by the ellipses.
// The SpatialObjectToImageFilter is used to that end. Note that
// this filter is instantiated over the spatial object used and the image
// type to be generated.
//
//
// \index{itk::Spatial\-Object\-To\-Image\-Filter!Instantiation}
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
typedef itk::SpatialObjectToImageFilter< GroupType, ImageType >
SpatialObjectToImageFilterType;
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// With the defined type, we construct a filter using the \code{New()}
// method. The newly created filter is assigned to a SmartPointer.
//
// \index{itk::SpatialObjectToImageFilter!New()}
// \index{itk::SpatialObjectToImageFilter!Pointer}
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
SpatialObjectToImageFilterType::Pointer imageFilter =
SpatialObjectToImageFilterType::New();
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// The GroupSpatialObject is passed as input to the filter.
//
// \index{itk::SpatialObjectToImageFilter!SetInput()}
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
imageFilter->SetInput( group );
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// The \doxygen{SpatialObjectToImageFilter} acts as a resampling filter.
// Therefore it requires the user to define the size of the desired output
// image. This is specified with the \code{SetSize()} method.
//
// \index{itk::SpatialObjectToImageFilter!SetSize()}
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
ImageType::SizeType size;
size[ 0 ] = 200;
size[ 1 ] = 200;
imageFilter->SetSize( size );
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// Finally we trigger the execution of the filter by calling the
// \code{Update()} method.
//
// \index{itk::SpatialObjectToImageFilter!Update()}
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
imageFilter->Update();
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// In order to obtain a smoother metric, we blur the image using a
// \doxygen{DiscreteGaussianImageFilter}. This extends the capture radius
// of the metric and produce a more continuous cost function to
// optimize. The following lines instantiate the Gaussian filter and
// create one object of this type using the \code{New()} method.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
typedef itk::DiscreteGaussianImageFilter< ImageType, ImageType >
GaussianFilterType;
GaussianFilterType::Pointer gaussianFilter = GaussianFilterType::New();
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// The output of the SpatialObjectToImageFilter is connected as
// input to the DiscreteGaussianImageFilter.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
gaussianFilter->SetInput( imageFilter->GetOutput() );
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// The variance of the filter is defined as a large value in order to
// increase the capture radius. Finally the execution of the filter is
// triggered using the \code{Update()} method.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
const double variance = 20;
gaussianFilter->SetVariance(variance);
gaussianFilter->Update();
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// Below we instantiate the type of the
// \doxygen{ImageToSpatialObjectRegistrationMethod} method and instantiate a
// registration object with the \code{New()} method. Note that the
// registration type is templated over the Image and the
// SpatialObject types. The spatial object in this case is the
// group of spatial objects.
//
// \index{itk::Image\-To\-Spatial\-Object\-Registration\-Method!Instantiation}
// \index{itk::Image\-To\-Spatial\-Object\-Registration\-Method!New()}
// \index{itk::Image\-To\-Spatial\-Object\-Registration\-Method!Pointer}
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
typedef itk::ImageToSpatialObjectRegistrationMethod< ImageType, GroupType >
RegistrationType;
RegistrationType::Pointer registration = RegistrationType::New();
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// Now we instantiate the metric that is templated over
// the image type and the spatial object type. As usual, the \code{New()}
// method is used to create an object.
//
// \index{itk::Image\-To\-Spatial\-Object\-Metric!Instantiation}
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
typedef SimpleImageToSpatialObjectMetric< ImageType, GroupType > MetricType;
MetricType::Pointer metric = MetricType::New();
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// An interpolator will be needed to evaluate the image at non-grid
// positions. Here we instantiate a linear interpolator type.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
typedef itk::LinearInterpolateImageFunction< ImageType, double >
InterpolatorType;
InterpolatorType::Pointer interpolator = InterpolatorType::New();
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// The following lines instantiate the evolutionary optimizer.
//
// \index{itk::One\-Plus\-One\-Evolutionary\-Optimizer!Instantiation}
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
typedef itk::OnePlusOneEvolutionaryOptimizer OptimizerType;
OptimizerType::Pointer optimizer = OptimizerType::New();
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// Next, we instantiate the transform class. In this case we use the
// Euler2DTransform that implements a rigid transform in $2D$
// space.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
typedef itk::Euler2DTransform<> TransformType;
TransformType::Pointer transform = TransformType::New();
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// Evolutionary algorithms are based on testing random variations
// of parameters. In order to support the computation of random values,
// ITK provides a family of random number generators. In this example, we
// use the \doxygen{NormalVariateGenerator} which generates values with a
// normal distribution.
//
// \index{itk::NormalVariateGenerator!New()}
// \index{itk::NormalVariateGenerator!Pointer}
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
itk::Statistics::NormalVariateGenerator::Pointer generator
= itk::Statistics::NormalVariateGenerator::New();
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// The random number generator must be initialized with a seed.
//
// \index{itk::NormalVariateGenerator!Initialize()}
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
generator->Initialize(12345);
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// The OnePlusOneEvolutionaryOptimizer is initialized by
// specifying the random number generator, the number of samples for the
// initial population and the maximum number of iterations.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
optimizer->SetNormalVariateGenerator( generator );
optimizer->Initialize( 10 );
optimizer->SetMaximumIteration( 400 );
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// As in previous registration examples, we take care to normalize the
// dynamic range of the different transform parameters. In particular, the
// we must compensate for the ranges of the angle and translations of the Euler2DTransform.
// In order to achieve this goal, we provide an array
// of scales to the optimizer.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
TransformType::ParametersType parametersScale;
parametersScale.set_size(3);
parametersScale[0] = 1000; // angle scale
for( unsigned int i=1; i<3; i++ )
{
parametersScale[i] = 2; // offset scale
}
optimizer->SetScales( parametersScale );
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// Here we instantiate the Command object that will act as an
// observer of the registration method and print out parameters at each
// iteration. Earlier, we defined this command as a class templated over the
// optimizer type. Once it is created with the \code{New()} method, we
// connect the optimizer to the command.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
typedef IterationCallback< OptimizerType > IterationCallbackType;
IterationCallbackType::Pointer callback = IterationCallbackType::New();
callback->SetOptimizer( optimizer );
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// All the components are plugged into the
// ImageToSpatialObjectRegistrationMethod object. The typical
// \code{Set()} methods are used here. Note the use of the
// \code{SetMovingSpatialObject()} method for connecting the spatial object.
// We provide the blurred version of the original synthetic binary
// image as the input image.
//
// \index{itk::Image\-To\-Spatial\-Object\-Registration\-Method!SetFixedImage()}
// \index{itk::Image\-To\-Spatial\-Object\-Registration\-Method!SetMovingSpatialObject()}
// \index{itk::Image\-To\-Spatial\-Object\-Registration\-Method!SetTransform()}
// \index{itk::Image\-To\-Spatial\-Object\-Registration\-Method!SetInterpolator()}
// \index{itk::Image\-To\-Spatial\-Object\-Registration\-Method!SetOptimizer()}
// \index{itk::Image\-To\-Spatial\-Object\-Registration\-Method!SetMetric()}
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
registration->SetFixedImage( gaussianFilter->GetOutput() );
registration->SetMovingSpatialObject( group );
registration->SetTransform( transform );
registration->SetInterpolator( interpolator );
registration->SetOptimizer( optimizer );
registration->SetMetric( metric );
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// The initial set of transform parameters is passed to the registration
// method using the \code{SetInitialTransformParameters()} method. Note that
// since our original model is already registered with the synthetic image,
// we introduce an artificial mis-registration in order to initialize
// the optimization at some point away from the optimal value.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
TransformType::ParametersType initialParameters(
transform->GetNumberOfParameters() );
initialParameters[0] = 0.2; // Angle
initialParameters[1] = 7.0; // Offset X
initialParameters[2] = 6.0; // Offset Y
registration->SetInitialTransformParameters(initialParameters);
// Software Guide : EndCodeSnippet
std::cout << "Initial Parameters : " << initialParameters << std::endl;
// Software Guide : BeginLatex
//
// Due to the character of the metric used to evaluate the fitness
// between the spatial object and the image, we must tell the optimizer that
// we are interested in finding the maximum value of the metric. Some
// metrics associate low numeric values with good matching, while others associate
// high numeric values with good matching. The \code{MaximizeOn()} and
// \code{MaximizeOff()} methods allow the user to deal with both types of
// metrics.
//
// \index{itk::Optimizer!MaximizeOn()}
// \index{itk::Optimizer!MaximizeOff()}
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
optimizer->MaximizeOn();
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// Finally, we trigger the execution of the registration process with the
// \code{Update()} method. We place this call in a
// \code{try/catch} block in case any exception is thrown during the
// process.
//
// \index{itk::Image\-To\-Spatial\-Object\-Registration\-Method!Update()}
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
try
{
registration->Update();
std::cout << "Optimizer stop condition: "
<< registration->GetOptimizer()->GetStopConditionDescription()
<< std::endl;
}
catch( itk::ExceptionObject & exp )
{
std::cerr << "Exception caught ! " << std::endl;
std::cerr << exp << std::endl;
}
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// The set of transform parameters resulting from the registration can be
// recovered with the \code{GetLastTransformParameters()} method. This
// method returns the array of transform parameters that should be
// interpreted according to the implementation of each transform. In our
// current example, the Euler2DTransform has three parameters:
// the rotation angle, the translation in $x$ and the translation in $y$.
//
// \index{itk::Image\-To\-Spatial\-Object\-Registration\-Method!Update()}
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
RegistrationType::ParametersType finalParameters
= registration->GetLastTransformParameters();
std::cout << "Final Solution is : " << finalParameters << std::endl;
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// \begin{figure}
// \center
// \includegraphics[height=0.44\textwidth]{ModelToImageRegistrationTraceAngle}
// \includegraphics[height=0.44\textwidth]{ModelToImageRegistrationTraceTranslations}
// \itkcaption[SpatialObject to Image Registration results]{Plots of the
// angle and translation parameters for a registration process between an
// spatial object and an image.}
// \label{fig:ModelToImageRegistrationPlots}
// \end{figure}
//
// The results are presented in
// Figure~\ref{fig:ModelToImageRegistrationPlots}. The left side shows the
// evolution of the angle parameter as a function of iteration
// numbers, while the right side shows the $(x,y)$ translation.
//
// Software Guide : EndLatex
return EXIT_SUCCESS;
}