← Back to branch summary

~ubuntu-branches/ubuntu/precise/insighttoolkit/precise

« back to all changes in this revision

Viewing changes to Testing/Code/Review/itkOptImageToImageMetricsTest.h

  • Committer: Bazaar Package Importer
  • Author(s): Steve M. Robbins
  • Date: 2008-12-19 20:16:49 UTC
  • mfrom: (1.2.1 upstream) (4.1.1 sid)
  • Revision ID: james.westby@ubuntu.com-20081219201649-drt97guwl2ryt0cn

* New upstream version.
  - patches/nifti-versioning.patch: Remove.  Applied upstream.
  - control:
  - rules: Update version numbers, package names.

* control: Build-depend on uuid-dev (gdcm uses it).

* copyright: Update download URL.

* rules: Adhere to parallel=N in DEB_BUILD_OPTIONS by setting MAKEFLAGS.

* compat: Set to 7.
* control: Update build-dep on debhelper to version >= 7.

* CMakeCache.txt.debian: Set CMAKE_BUILD_TYPE to "RELEASE" so that we
  build with -O3 (not -O2), necessary to optimize the templated code.

Show diffs side-by-side

added added

removed removed

Lines of Context:
Testing/Code/Review/itkOptImageToImageMetricsTest.h
 
1
/*=========================================================================
 
2
 
 
3
  Program:   Insight Segmentation & Registration Toolkit
 
4
  Module:    $RCSfile: itkOptImageToImageMetricsTest.h,v $
 
5
  Language:  C++
 
6
  Date:      $Date: 2007-12-28 15:04:56 $
 
7
  Version:   $Revision: 1.4 $
 
8
 
 
9
  Copyright (c) Insight Software Consortium. All rights reserved.
 
10
  See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details.
 
11
 
 
12
     This software is distributed WITHOUT ANY WARRANTY; without even 
 
13
     the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR 
 
14
     PURPOSE.  See the above copyright notices for more information.
 
15
 
 
16
=========================================================================*/
 
17
 
 
18
#ifndef __itkOptImageToImageMetricsTest_h
 
19
#define __itkOptImageToImageMetricsTest_h
 
20
 
 
21
#include "itkTimeProbe.h"
 
22
#include "itkMersenneTwisterRandomVariateGenerator.h"
 
23
 
 
24
namespace itk {
 
25
 
 
26
template <typename FixedImageType, 
 
27
          typename MovingImageType, 
 
28
          typename InterpolatorType,
 
29
          typename TransformType,
 
30
          typename MetricType, 
 
31
          typename MetricInitializerType >
 
32
class OptImageToImageMetricsTest
 
33
{
 
34
public:
 
35
 
 
36
  OptImageToImageMetricsTest() {}
 
37
 
 
38
  int RunTest( FixedImageType* fixed, 
 
39
               MovingImageType* moving,
 
40
               InterpolatorType* interpolator,
 
41
               TransformType* transform,
 
42
               MetricType* metric,
 
43
               MetricInitializerType metricInitializer)
 
44
    {
 
45
    typedef typename MetricType::ParametersType ParametersType;
 
46
 
 
47
    std::cout << "-------------------------------------------------------------------" << std::endl;
 
48
    std::cout << "Testing" << std::endl;
 
49
    std::cout << "\tMetric       : " << metric->GetNameOfClass() << std::endl;
 
50
    std::cout << "\tInterpolator : " << interpolator->GetNameOfClass() << std::endl;
 
51
    std::cout << "\tTransform    : " << transform->GetNameOfClass() << std::endl;
 
52
    std::cout << "-------------------------------------------------------------------" << std::endl;
 
53
    std::cout << std::endl;
 
54
 
 
55
    int result = EXIT_SUCCESS;
 
56
 
 
57
    // connect the interpolator
 
58
    metric->SetInterpolator( interpolator );
 
59
 
 
60
    // connect the transform
 
61
    metric->SetTransform( transform );
 
62
 
 
63
    // connect the images to the metric
 
64
    metric->SetFixedImage( fixed );
 
65
    metric->SetMovingImage( moving );
 
66
 
 
67
    // call custom initialization for the metric
 
68
    metricInitializer.Initialize();
 
69
 
 
70
    // Always use the same seed value.
 
71
    // All instances are the same since MersenneTwisterRandomVariateGenerator
 
72
    // uses a singleton pattern.
 
73
    itk::Statistics::MersenneTwisterRandomVariateGenerator::GetInstance()->SetSeed( 42 );
 
74
 
 
75
    // initialize the metric
 
76
    // Samples are drawn here in metric->Initialize(), 
 
77
    // so we seed the random number generator
 
78
    // immediately before this call.
 
79
    metric->Initialize();
 
80
 
 
81
    // Set the transform to identity
 
82
    transform->SetIdentity();
 
83
 
 
84
    // Get the transform parameters for identity.
 
85
    ParametersType parameters = transform->GetParameters();
 
86
 
 
87
    typename MetricType::MeasureType value;
 
88
    typename MetricType::DerivativeType derivative;
 
89
 
 
90
    // Try GetValue and GetDerivative...
 
91
    value = metric->GetValue( parameters );
 
92
    metric->GetDerivative( parameters, derivative );
 
93
 
 
94
    // Make a time probe
 
95
    itk::TimeProbe timeProbe;
 
96
 
 
97
    // Walk around the parameter value at parameterIdx
 
98
    for (unsigned int parameterIdx = 0; parameterIdx < parameters.GetSize(); parameterIdx++)
 
99
      {
 
100
      std::cout << "Param[" << parameterIdx << "]\tValue\tDerivative " << std::endl;
 
101
      double startVal = parameters[parameterIdx];
 
102
      // endVal is 10% beyond startVal.
 
103
      double endVal = 1.10 * startVal;
 
104
      // If startVal is 0, endVal needs to be fixed up.
 
105
      if( fabs(endVal - 0.0) < 1e-8 )
 
106
        {
 
107
        endVal = startVal + 1.0;
 
108
        }
 
109
      double incr = (endVal - startVal) / 10.0;
 
110
 
 
111
      for( double pval = startVal; pval <= endVal; pval += incr )
 
112
        {
 
113
        parameters[parameterIdx] = pval;
 
114
 
 
115
        timeProbe.Start();
 
116
        metric->GetValueAndDerivative( parameters, value, derivative );
 
117
        timeProbe.Stop();
 
118
 
 
119
        std::cout << pval << "\t" << value << "\t" << derivative << std::endl;
 
120
        }
 
121
      }
 
122
 
 
123
    std::cout << std::endl;
 
124
    std::cout << "Mean time for GetValueAndDerivative : " << timeProbe.GetMeanTime() << std::endl;
 
125
    std::cout << std::endl;
 
126
    std::cout << "------------------------------Done---------------------------------" << std::endl;
 
127
 
 
128
    return result;
 
129
    }
 
130
 
 
131
};
 
132
 
 
133
template <typename FixedImageType, typename MovingImageType>
 
134
class MeanSquaresMetricInitializer
 
135
{
 
136
public:
 
137
  typedef itk::MeanSquaresImageToImageMetric< FixedImageType, 
 
138
                                              MovingImageType> MetricType;
 
139
 
 
140
 
 
141
  MeanSquaresMetricInitializer(MetricType* metric)
 
142
    {
 
143
    m_Metric = metric;
 
144
    }
 
145
 
 
146
  void Initialize()
 
147
    {
 
148
    // Do stuff on m_Metric
 
149
#ifdef ITK_USE_OPTIMIZED_REGISTRATION_METHODS
 
150
    m_Metric->UseAllPixelsOn();
 
151
#endif 
 
152
    }
 
153
 
 
154
protected:
 
155
  MetricType* m_Metric;
 
156
 
 
157
};
 
158
 
 
159
template <typename FixedImageType, typename MovingImageType>
 
160
class MattesMIMetricInitializer
 
161
{
 
162
public:
 
163
  typedef itk::MattesMutualInformationImageToImageMetric< FixedImageType, 
 
164
                                                          MovingImageType> MetricType;
 
165
 
 
166
 
 
167
  MattesMIMetricInitializer(MetricType* metric)
 
168
    {
 
169
    m_Metric = metric;
 
170
    }
 
171
 
 
172
  void Initialize()
 
173
    {
 
174
    // Do stuff on m_Metric
 
175
    m_Metric->SetNumberOfHistogramBins( 50 );
 
176
    m_Metric->SetNumberOfSpatialSamples( 5000 );
 
177
    }
 
178
 
 
179
protected:
 
180
  MetricType* m_Metric;
 
181
 
 
182
};
 
183
 
 
184
template <typename FixedImageType, typename MovingImageType>
 
185
class MIMetricInitializer
 
186
{
 
187
public:
 
188
  typedef itk::MutualInformationImageToImageMetric< FixedImageType, 
 
189
                                                    MovingImageType> MetricType;
 
190
 
 
191
 
 
192
  MIMetricInitializer(MetricType* metric)
 
193
    {
 
194
    m_Metric = metric;
 
195
    }
 
196
 
 
197
  void Initialize()
 
198
    {
 
199
    // Do stuff on m_Metric
 
200
    m_Metric->SetNumberOfSpatialSamples( 400 );
 
201
    }
 
202
 
 
203
protected:
 
204
  MetricType* m_Metric;
 
205
 
 
206
};
 
207
 
 
208
template < class InterpolatorType, 
 
209
            class TransformType,
 
210
            class FixedImageReaderType,
 
211
            class MovingImageReaderType >
 
212
void BasicTest( FixedImageReaderType* fixedImageReader, 
 
213
                MovingImageReaderType* movingImageReader,
 
214
                InterpolatorType* interpolator,
 
215
                TransformType* transform
 
216
                )
 
217
{
 
218
  typedef typename FixedImageReaderType::OutputImageType    FixedImageType;
 
219
  typedef typename MovingImageReaderType::OutputImageType   MovingImageType;
 
220
 
 
221
  fixedImageReader->Update();
 
222
  movingImageReader->Update();
 
223
 
 
224
  typename FixedImageType::Pointer fixed = fixedImageReader->GetOutput();
 
225
  typename MovingImageType::Pointer moving = movingImageReader->GetOutput();
 
226
 
 
227
  // Mean squares 
 
228
  typedef itk::MeanSquaresImageToImageMetric< FixedImageType, MovingImageType > MetricType;
 
229
  typedef MeanSquaresMetricInitializer< FixedImageType, MovingImageType > MetricInitializerType;
 
230
  typename MetricType::Pointer msMetric = MetricType::New();
 
231
  MeanSquaresMetricInitializer< FixedImageType, MovingImageType > msMetricInitializer( msMetric );
 
232
 
 
233
  TestAMetric( fixedImageReader, movingImageReader, interpolator, transform, msMetric.GetPointer(), msMetricInitializer );
 
234
 
 
235
  // Mattes MI
 
236
  typedef itk::MattesMutualInformationImageToImageMetric< FixedImageType, MovingImageType > MattesMetricType;
 
237
  typedef MattesMIMetricInitializer< FixedImageType, MovingImageType > MattesMetricInitializerType;
 
238
  typename MattesMetricType::Pointer mattesMetric = MattesMetricType::New();
 
239
  MattesMIMetricInitializer< FixedImageType, MovingImageType > mattesMetricInitializer( mattesMetric );
 
240
 
 
241
  TestAMetric( fixedImageReader, movingImageReader, interpolator, transform, mattesMetric.GetPointer(), mattesMetricInitializer );
 
242
 
 
243
#if 0 // OptMutualInformationImageToImageMetric will be removed from Review.
 
244
  // MI
 
245
  typedef itk::MutualInformationImageToImageMetric< FixedImageType, MovingImageType > MIMetricType;
 
246
  typedef MIMetricInitializer< FixedImageType, MovingImageType > MIMetricInitializerType;
 
247
  typename MIMetricType::Pointer miMetric = MIMetricType::New();
 
248
  MIMetricInitializer< FixedImageType, MovingImageType> miMetricInitializer( miMetric );
 
249
  
 
250
  TestAMetric( fixedImageReader, movingImageReader, interpolator, transform, miMetric.GetPointer(), miMetricInitializer );
 
251
#endif 
 
252
}
 
253
 
 
254
template <class FixedImageReaderType,
 
255
          class MovingImageReaderType,
 
256
          class InterpolatorType,
 
257
          class TransformType,
 
258
          class MetricType,
 
259
          class MetricInitializerType>
 
260
void TestAMetric(FixedImageReaderType* fixedImageReader,
 
261
                 MovingImageReaderType* movingImageReader,
 
262
                 InterpolatorType* interpolator,
 
263
                 TransformType* transform,
 
264
                 MetricType* metric,
 
265
                 MetricInitializerType metricInitializer)
 
266
{
 
267
  typedef typename FixedImageReaderType::OutputImageType    FixedImageType;
 
268
  typedef typename MovingImageReaderType::OutputImageType   MovingImageType;
 
269
 
 
270
  metric->SetFixedImageRegion( fixedImageReader->GetOutput()->GetBufferedRegion() );
 
271
 
 
272
  OptImageToImageMetricsTest<  FixedImageType,
 
273
                          MovingImageType,
 
274
                          InterpolatorType,
 
275
                          TransformType,
 
276
                          MetricType,
 
277
                          MetricInitializerType > testMetric;
 
278
 
 
279
  testMetric.RunTest( fixedImageReader->GetOutput(), movingImageReader->GetOutput(), interpolator, transform, metric, metricInitializer );
 
280
}
 
281
 
 
282
template <class FixedImageReaderType, class MovingImageReaderType>
 
283
void AffineLinearTest( FixedImageReaderType* fixedImageReader, 
 
284
                       MovingImageReaderType* movingImageReader)
 
285
{
 
286
  typedef typename MovingImageReaderType::OutputImageType MovingImageType;
 
287
  typedef itk::LinearInterpolateImageFunction< MovingImageType, double > InterpolatorType;
 
288
  typedef itk::AffineTransform<double, 2> TransformType;
 
289
 
 
290
  typename InterpolatorType::Pointer interpolator = InterpolatorType::New();
 
291
  TransformType::Pointer transform = TransformType::New();
 
292
 
 
293
  BasicTest(fixedImageReader,
 
294
            movingImageReader, 
 
295
            interpolator.GetPointer(), 
 
296
            transform.GetPointer());
 
297
 
 
298
}
 
299
 
 
300
template <class FixedImageReaderType, class MovingImageReaderType>
 
301
void RigidLinearTest( FixedImageReaderType* fixedImageReader, 
 
302
                       MovingImageReaderType* movingImageReader)
 
303
{
 
304
  typedef typename MovingImageReaderType::OutputImageType MovingImageType;
 
305
 
 
306
  typedef itk::LinearInterpolateImageFunction< MovingImageType, double > InterpolatorType;
 
307
  typedef itk::Rigid2DTransform<double> TransformType;
 
308
 
 
309
  typename InterpolatorType::Pointer interpolator = InterpolatorType::New();
 
310
  TransformType::Pointer transform = TransformType::New();
 
311
 
 
312
  BasicTest(fixedImageReader,
 
313
            movingImageReader, 
 
314
            interpolator.GetPointer(), 
 
315
            transform.GetPointer());
 
316
}
 
317
 
 
318
template <class FixedImageReaderType, class MovingImageReaderType>
 
319
void TranslationLinearTest( FixedImageReaderType* fixedImageReader, 
 
320
                            MovingImageReaderType* movingImageReader)
 
321
{
 
322
  typedef typename MovingImageReaderType::OutputImageType MovingImageType;
 
323
 
 
324
  typedef itk::LinearInterpolateImageFunction< MovingImageType, double > InterpolatorType;
 
325
  typedef itk::TranslationTransform<double, 2> TransformType;
 
326
 
 
327
  typename InterpolatorType::Pointer interpolator = InterpolatorType::New();
 
328
  TransformType::Pointer transform = TransformType::New();
 
329
 
 
330
  BasicTest(fixedImageReader,
 
331
            movingImageReader, 
 
332
            interpolator.GetPointer(), 
 
333
            transform.GetPointer());
 
334
}
 
335
 
 
336
 
 
337
template <class FixedImageReaderType, class MovingImageReaderType>
 
338
void DoDebugTest( FixedImageReaderType* fixedImageReader,
 
339
                  MovingImageReaderType* movingImageReader)
 
340
{
 
341
  typedef typename MovingImageReaderType::OutputImageType MovingImageType;
 
342
 
 
343
  typedef itk::LinearInterpolateImageFunction< MovingImageType, double > InterpolatorType;
 
344
  typedef itk::Rigid2DTransform<double> TransformType;
 
345
 
 
346
  typename InterpolatorType::Pointer interpolator = InterpolatorType::New();
 
347
  TransformType::Pointer transform = TransformType::New();
 
348
 
 
349
  typedef typename FixedImageReaderType::OutputImageType      FixedImageType;
 
350
  typedef typename MovingImageReaderType::OutputImageType     MovingImageType;
 
351
 
 
352
  fixedImageReader->Update();
 
353
  movingImageReader->Update();
 
354
 
 
355
  typename FixedImageType::Pointer fixed = fixedImageReader->GetOutput();
 
356
  typename MovingImageType::Pointer moving = movingImageReader->GetOutput();
 
357
 
 
358
  // Mean squares 
 
359
  typedef itk::MeanSquaresImageToImageMetric< FixedImageType, MovingImageType > MetricType;
 
360
  typedef MeanSquaresMetricInitializer< FixedImageType, MovingImageType > MetricInitializerType;
 
361
  typename MetricType::Pointer metric = MetricType::New();
 
362
  MeanSquaresMetricInitializer< FixedImageType, MovingImageType > metricInitializer( metric );
 
363
 
 
364
  metric->SetFixedImageRegion( fixedImageReader->GetOutput()->GetBufferedRegion() );
 
365
 
 
366
  typedef typename MetricType::ParametersType ParametersType;
 
367
 
 
368
  std::cout << "-------------------------------------------------------------------" << std::endl;
 
369
  std::cout << "Testing" << std::endl;
 
370
  std::cout << "\tMetric       : " << metric->GetNameOfClass() << std::endl;
 
371
  std::cout << "\tInterpolator : " << interpolator->GetNameOfClass() << std::endl;
 
372
  std::cout << "\tTransform    : " << transform->GetNameOfClass() << std::endl;
 
373
  std::cout << "-------------------------------------------------------------------" << std::endl;
 
374
  std::cout << std::endl;
 
375
 
 
376
  int result = EXIT_SUCCESS;
 
377
 
 
378
  // connect the interpolator
 
379
  metric->SetInterpolator( interpolator );
 
380
 
 
381
  // connect the transform
 
382
  metric->SetTransform( transform );
 
383
 
 
384
  // connect the images to the metric
 
385
  metric->SetFixedImage( fixed );
 
386
  metric->SetMovingImage( moving );
 
387
 
 
388
  // call custom initialization for the metric
 
389
  metricInitializer.Initialize();
 
390
 
 
391
  // initialize the metric
 
392
  metric->Initialize();
 
393
 
 
394
  // Set the transform to identity
 
395
  transform->SetIdentity();
 
396
 
 
397
  // Get the transform parameters for identity.
 
398
  ParametersType parameters = transform->GetParameters();
 
399
 
 
400
  typename MetricType::MeasureType value;
 
401
  typename MetricType::DerivativeType derivative;
 
402
 
 
403
  parameters[0] = 0.1;
 
404
 
 
405
  metric->GetValueAndDerivative( parameters,
 
406
                                 value, 
 
407
                                 derivative );
 
408
 
 
409
  //metric->StopDebug();
 
410
 
 
411
  // Force the test to end here so the debug file
 
412
  // ends at the right place.
 
413
  exit(EXIT_SUCCESS);
 
414
}
 
415
 
 
416
} // end namespace itk
 
417
 
 
418
#endif