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 #ifndef itkImageToImageMetric_hxx
19 #define itkImageToImageMetric_hxx
20 
21 #include "itkImageToImageMetric.h"
22 #include "itkImageRandomConstIteratorWithIndex.h"
23 #include "itkMath.h"
24 
25 namespace itk
26 {
27 /**
28  * Constructor
29  */
30 template< typename TFixedImage, typename TMovingImage >
31 ImageToImageMetric< TFixedImage, TMovingImage >
ImageToImageMetric()32 ::ImageToImageMetric():
33   m_FixedImageIndexes(0),
34   m_FixedImageSamplesIntensityThreshold(0),
35   m_FixedImageSamples(0),
36   m_FixedImage(nullptr), // has to be provided by the user.
37   m_MovingImage(nullptr), // has to be provided by the user.
38   m_Transform(nullptr), // has to be provided by the user.
39   m_ThreaderTransform(nullptr), // constructed at initialization.
40   m_Interpolator(nullptr),  // metric computes gradient by default
41   m_GradientImage(nullptr),   // computed at initialization
42   m_FixedImageMask(nullptr),
43   m_MovingImageMask(nullptr),
44   m_RandomSeed(Statistics::MersenneTwisterRandomVariateGenerator::GetNextSeed()),
45   m_BSplineTransform(nullptr),
46   m_BSplineTransformWeightsArray(),
47   m_BSplineTransformIndicesArray(),
48   m_BSplinePreTransformPointsArray(0),
49   m_WithinBSplineSupportRegionArray(0),
50   m_BSplineParametersOffset(),
51   m_BSplineTransformWeights(),
52   m_BSplineTransformIndices(),
53   m_ThreaderBSplineTransformWeights(nullptr),
54   m_ThreaderBSplineTransformIndices(nullptr),
55   m_BSplineInterpolator(nullptr),
56   m_DerivativeCalculator(nullptr),
57   m_Threader(MultiThreaderType::New())
58 {
59   m_ConstSelfWrapper = new ConstantPointerWrapper(this);
60   this->m_NumberOfWorkUnits = this->m_Threader->GetNumberOfWorkUnits();
61 
62   /* if 100% backward compatible, we should include this...but...
63   typename BSplineTransformType::Pointer transformer =
64            BSplineTransformType::New();
65   this->SetTransform (transformer);
66 
67   typename BSplineInterpolatorType::Pointer interpolator =
68            BSplineInterpolatorType::New();
69   this->SetInterpolator (interpolator);
70   */
71 }
72 
73 template< typename TFixedImage, typename TMovingImage >
74 ImageToImageMetric< TFixedImage, TMovingImage >
~ImageToImageMetric()75 ::~ImageToImageMetric()
76 {
77   delete m_ConstSelfWrapper;
78 
79   delete[] m_ThreaderNumberOfMovingImageSamples;
80   m_ThreaderNumberOfMovingImageSamples = nullptr;
81 
82   delete[] m_ThreaderTransform;
83   m_ThreaderTransform = nullptr;
84 
85   delete[] this->m_ThreaderBSplineTransformWeights;
86   this->m_ThreaderBSplineTransformWeights = nullptr;
87 
88   delete[] this->m_ThreaderBSplineTransformIndices;
89   this->m_ThreaderBSplineTransformIndices = nullptr;
90 }
91 
92 /**
93  * Set the number of work units. This will be clamped by the
94  * multithreader, so we must check to see if it is accepted.
95  */
96 template< typename TFixedImage, typename TMovingImage >
97 void
98 ImageToImageMetric< TFixedImage, TMovingImage >
SetNumberOfWorkUnits(ThreadIdType numberOfThreads)99 ::SetNumberOfWorkUnits(ThreadIdType numberOfThreads)
100 {
101   m_Threader->SetNumberOfWorkUnits( numberOfThreads );
102   m_NumberOfWorkUnits = m_Threader->GetNumberOfWorkUnits();
103 }
104 
105 /**
106  * Set the parameters that define a unique transform
107  */
108 template< typename TFixedImage, typename TMovingImage >
109 void
110 ImageToImageMetric< TFixedImage, TMovingImage >
SetTransformParameters(const ParametersType & parameters) const111 ::SetTransformParameters(const ParametersType & parameters) const
112 {
113   if ( !m_Transform )
114     {
115     itkExceptionMacro(<< "Transform has not been assigned");
116     }
117   m_Transform->SetParameters(parameters);
118 
119 }
120 
121 template< typename TFixedImage, typename TMovingImage >
122 void
123 ImageToImageMetric< TFixedImage, TMovingImage >
SetNumberOfFixedImageSamples(SizeValueType numSamples)124 ::SetNumberOfFixedImageSamples(SizeValueType numSamples)
125 {
126   if ( numSamples != m_NumberOfFixedImageSamples )
127     {
128     m_NumberOfFixedImageSamples = numSamples;
129     if ( m_NumberOfFixedImageSamples != this->m_FixedImageRegion.GetNumberOfPixels() )
130       {
131       this->SetUseAllPixels(false);
132       }
133     this->Modified();
134     }
135 }
136 
137 template< typename TFixedImage, typename TMovingImage >
138 void
139 ImageToImageMetric< TFixedImage, TMovingImage >
SetFixedImageIndexes(const FixedImageIndexContainer & indexes)140 ::SetFixedImageIndexes(const FixedImageIndexContainer & indexes)
141 {
142   this->SetUseFixedImageIndexes(true);
143   m_NumberOfFixedImageSamples = indexes.size();
144   m_FixedImageIndexes.resize(m_NumberOfFixedImageSamples);
145   for ( unsigned int i = 0; i < m_NumberOfFixedImageSamples; i++ )
146     {
147     m_FixedImageIndexes[i] = indexes[i];
148     }
149 }
150 
151 template< typename TFixedImage, typename TMovingImage >
152 void
153 ImageToImageMetric< TFixedImage, TMovingImage >
SetUseFixedImageIndexes(bool useIndexes)154 ::SetUseFixedImageIndexes(bool useIndexes)
155 {
156   if ( useIndexes != m_UseFixedImageIndexes )
157     {
158     m_UseFixedImageIndexes = useIndexes;
159     if ( m_UseFixedImageIndexes )
160       {
161       this->SetUseAllPixels(false);
162       }
163     else
164       {
165       this->Modified();
166       }
167     }
168 }
169 
170 template< typename TFixedImage, typename TMovingImage >
171 void
172 ImageToImageMetric< TFixedImage, TMovingImage >
SetFixedImageSamplesIntensityThreshold(const FixedImagePixelType & thresh)173 ::SetFixedImageSamplesIntensityThreshold(const FixedImagePixelType & thresh)
174 {
175   if ( thresh != m_FixedImageSamplesIntensityThreshold )
176     {
177     m_FixedImageSamplesIntensityThreshold = thresh;
178     this->SetUseFixedImageSamplesIntensityThreshold(true);
179     this->Modified();
180     }
181 }
182 
183 template< typename TFixedImage, typename TMovingImage >
184 void
185 ImageToImageMetric< TFixedImage, TMovingImage >
SetUseFixedImageSamplesIntensityThreshold(bool useThresh)186 ::SetUseFixedImageSamplesIntensityThreshold(bool useThresh)
187 {
188   if ( useThresh != m_UseFixedImageSamplesIntensityThreshold )
189     {
190     m_UseFixedImageSamplesIntensityThreshold = useThresh;
191     if ( m_UseFixedImageSamplesIntensityThreshold )
192       {
193       this->SetUseAllPixels(false);
194       }
195     else
196       {
197       this->Modified();
198       }
199     }
200 }
201 
202 template< typename TFixedImage, typename TMovingImage >
203 void
204 ImageToImageMetric< TFixedImage, TMovingImage >
SetFixedImageRegion(const FixedImageRegionType reg)205 ::SetFixedImageRegion(const FixedImageRegionType reg)
206 {
207   if ( reg != m_FixedImageRegion )
208     {
209     m_FixedImageRegion = reg;
210     if ( this->GetUseAllPixels() )
211       {
212       this->SetNumberOfFixedImageSamples( this->m_FixedImageRegion.GetNumberOfPixels() );
213       }
214     }
215 }
216 
217 template< typename TFixedImage, typename TMovingImage >
218 void
219 ImageToImageMetric< TFixedImage, TMovingImage >
SetUseAllPixels(bool useAllPixels)220 ::SetUseAllPixels(bool useAllPixels)
221 {
222   if ( useAllPixels != m_UseAllPixels )
223     {
224     m_UseAllPixels = useAllPixels;
225     if ( m_UseAllPixels )
226       {
227       this->SetUseFixedImageSamplesIntensityThreshold(false);
228       this->SetNumberOfFixedImageSamples( this->m_FixedImageRegion.GetNumberOfPixels() );
229       this->SetUseSequentialSampling(true);
230       }
231     else
232       {
233       this->SetUseSequentialSampling(false);
234       this->Modified();
235       }
236     }
237 }
238 
239 template< typename TFixedImage, typename TMovingImage >
240 void
241 ImageToImageMetric< TFixedImage, TMovingImage >
SetUseSequentialSampling(bool useSequential)242 ::SetUseSequentialSampling(bool useSequential)
243 {
244   if ( useSequential != m_UseSequentialSampling )
245     {
246     m_UseSequentialSampling = useSequential;
247     if ( !m_UseSequentialSampling )
248       {
249       this->SetUseAllPixels(false);
250       }
251     else
252       {
253       this->Modified();
254       }
255     }
256 }
257 
258 /**
259  * Initialize
260  */
261 template< typename TFixedImage, typename TMovingImage >
262 void
263 ImageToImageMetric< TFixedImage, TMovingImage >
Initialize()264 ::Initialize()
265 {
266   if ( !m_Transform )
267     {
268     itkExceptionMacro(<< "Transform is not present");
269     }
270   m_NumberOfParameters = m_Transform->GetNumberOfParameters();
271 
272   if ( !m_Interpolator )
273     {
274     itkExceptionMacro(<< "Interpolator is not present");
275     }
276 
277   if ( !m_MovingImage )
278     {
279     itkExceptionMacro(<< "MovingImage is not present");
280     }
281 
282   if ( !m_FixedImage )
283     {
284     itkExceptionMacro(<< "FixedImage is not present");
285     }
286 
287   // If the image is provided by a source, update the source.
288   if ( m_MovingImage->GetSource() )
289     {
290     m_MovingImage->GetSource()->Update();
291     }
292 
293   // If the image is provided by a source, update the source.
294   if ( m_FixedImage->GetSource() )
295     {
296     m_FixedImage->GetSource()->Update();
297     }
298 
299   //The use of FixedImageIndexes and the use of FixedImageRegion
300   //are mutually exclusive, so they should not both be checked.
301   if ( this->m_UseFixedImageIndexes  == true )
302     {
303     if( this->m_FixedImageIndexes.empty() )
304       {
305       itkExceptionMacro(<< "FixedImageIndexes list is empty");
306       }
307     }
308   else
309     {
310     // Make sure the FixedImageRegion is within the FixedImage buffered region
311     if ( m_FixedImageRegion.GetNumberOfPixels() == 0 )
312       {
313       itkExceptionMacro(<< "FixedImageRegion is empty");
314       }
315 
316     if ( !m_FixedImageRegion.Crop( m_FixedImage->GetBufferedRegion() ) )
317       {
318       itkExceptionMacro(
319         << "FixedImageRegion does not overlap the fixed image buffered region");
320       }
321     }
322 
323   m_Interpolator->SetInputImage(m_MovingImage);
324 
325   if ( m_ComputeGradient )
326     {
327     ComputeGradient();
328     }
329 
330   // If there are any observers on the metric, call them to give the
331   // user code a chance to set parameters on the metric
332   this->InvokeEvent( InitializeEvent() );
333 }
334 
335 /**
336  * MultiThreading Initialize
337  */
338 template< typename TFixedImage, typename TMovingImage >
339 void
340 ImageToImageMetric< TFixedImage, TMovingImage >
MultiThreadingInitialize()341 ::MultiThreadingInitialize()
342 {
343   this->SetNumberOfWorkUnits(m_NumberOfWorkUnits);
344 
345   delete[] m_ThreaderNumberOfMovingImageSamples;
346   m_ThreaderNumberOfMovingImageSamples = new unsigned int[m_NumberOfWorkUnits - 1];
347 
348   // Allocate the array of transform clones to be used in every thread
349   delete[] m_ThreaderTransform;
350   m_ThreaderTransform = new TransformPointer[m_NumberOfWorkUnits - 1];
351   for ( ThreadIdType ithread = 0; ithread < m_NumberOfWorkUnits - 1; ++ithread )
352     {
353     this->m_ThreaderTransform[ithread] = this->m_Transform->Clone();
354     }
355 
356   m_FixedImageSamples.resize(m_NumberOfFixedImageSamples);
357   if ( m_UseSequentialSampling )
358     {
359     //
360     // Take all the pixels within the fixed image region)
361     // to create the sample points list.
362     //
363     SampleFullFixedImageRegion(m_FixedImageSamples);
364     }
365   else
366     {
367     if ( m_UseFixedImageIndexes )
368       {
369       //
370       //  Use the list of indexes passed to the SetFixedImageIndexes
371       //  member function .
372       //
373       SampleFixedImageIndexes(m_FixedImageSamples);
374       }
375     else
376       {
377       //
378       // Uniformly sample the fixed image (within the fixed image region)
379       // to create the sample points list.
380       //
381       SampleFixedImageRegion(m_FixedImageSamples);
382       }
383     }
384 
385   //
386   //  Check if the interpolator is of type BSplineInterpolateImageFunction.
387   //  If so, we can make use of its EvaluateDerivatives method.
388   //  Otherwise, we instantiate an external central difference
389   //  derivative calculator.
390   //
391   m_InterpolatorIsBSpline = true;
392 
393   auto * testPtr = dynamic_cast< BSplineInterpolatorType * >( this->m_Interpolator.GetPointer() );
394   if ( !testPtr )
395     {
396     m_InterpolatorIsBSpline = false;
397 
398     m_DerivativeCalculator = DerivativeFunctionType::New();
399     m_DerivativeCalculator->UseImageDirectionOn();
400 
401     m_DerivativeCalculator->SetInputImage(this->m_MovingImage);
402 
403     m_BSplineInterpolator = nullptr;
404     itkDebugMacro("Interpolator is not BSpline");
405     }
406   else
407     {
408     m_BSplineInterpolator = testPtr;
409     m_BSplineInterpolator->SetNumberOfWorkUnits(m_NumberOfWorkUnits);
410     m_BSplineInterpolator->UseImageDirectionOn();
411 
412     m_DerivativeCalculator = nullptr;
413     itkDebugMacro("Interpolator is BSpline");
414     }
415 
416   //
417   //  Check if the transform is of type BSplineTransform.
418   //
419   //  If so, several speed up features are implemented.
420   //  [1] Precomputing the results of bulk transform for each sample point.
421   //  [2] Precomputing the BSpline weights for each sample point,
422   //      to be used later to directly compute the deformation vector
423   //  [3] Precomputing the indices of the parameters within the
424   //      the support region of each sample point.
425   //
426   m_TransformIsBSpline = true;
427 
428   auto * testPtr2 = dynamic_cast< BSplineTransformType * >( this->m_Transform.GetPointer() );
429   if ( !testPtr2 )
430     {
431     m_TransformIsBSpline = false;
432     m_BSplineTransform = nullptr;
433     itkDebugMacro("Transform is not BSplineDeformable");
434     }
435   else
436     {
437     m_BSplineTransform = testPtr2;
438     m_NumBSplineWeights = m_BSplineTransform->GetNumberOfWeights();
439     itkDebugMacro("Transform is BSplineDeformable");
440     }
441 
442   if ( this->m_TransformIsBSpline )
443     {
444     // First, deallocate memory that may have been used from previous run of the Metric
445     this->m_BSplineTransformWeightsArray.SetSize(1, 1);
446     this->m_BSplineTransformIndicesArray.SetSize(1, 1);
447     this->m_BSplinePreTransformPointsArray.resize(1);
448     this->m_WithinBSplineSupportRegionArray.resize(1);
449     this->m_BSplineTransformWeights.SetSize(1);
450     this->m_BSplineTransformIndices.SetSize(1);
451 
452     delete[] this->m_ThreaderBSplineTransformWeights;
453     this->m_ThreaderBSplineTransformWeights = nullptr;
454 
455     delete[] this->m_ThreaderBSplineTransformIndices;
456     this->m_ThreaderBSplineTransformIndices = nullptr;
457 
458     if ( this->m_UseCachingOfBSplineWeights )
459       {
460       m_BSplineTransformWeightsArray.SetSize(
461         m_NumberOfFixedImageSamples, m_NumBSplineWeights);
462       m_BSplineTransformIndicesArray.SetSize(
463         m_NumberOfFixedImageSamples, m_NumBSplineWeights);
464       m_BSplinePreTransformPointsArray.resize(m_NumberOfFixedImageSamples);
465       m_WithinBSplineSupportRegionArray.resize(m_NumberOfFixedImageSamples);
466 
467       this->PreComputeTransformValues();
468       }
469     else
470       {
471       this->m_BSplineTransformWeights.SetSize(this->m_NumBSplineWeights);
472       this->m_BSplineTransformIndices.SetSize(this->m_NumBSplineWeights);
473 
474       this->m_ThreaderBSplineTransformWeights = new BSplineTransformWeightsType[m_NumberOfWorkUnits - 1];
475       this->m_ThreaderBSplineTransformIndices = new BSplineTransformIndexArrayType[m_NumberOfWorkUnits - 1];
476 
477       for ( ThreadIdType ithread = 0; ithread < m_NumberOfWorkUnits - 1; ++ithread )
478         {
479         this->m_ThreaderBSplineTransformWeights[ithread].SetSize(this->m_NumBSplineWeights);
480         this->m_ThreaderBSplineTransformIndices[ithread].SetSize(this->m_NumBSplineWeights);
481         }
482       }
483 
484     for ( unsigned int j = 0; j < FixedImageDimension; j++ )
485       {
486       this->m_BSplineParametersOffset[j] = j * this->m_BSplineTransform->GetNumberOfParametersPerDimension();
487       }
488     }
489 }
490 
491 /**
492  * Use the indexes that have been passed to the metric
493  */
494 template< typename TFixedImage, typename TMovingImage >
495 void
496 ImageToImageMetric< TFixedImage, TMovingImage >
SampleFixedImageIndexes(FixedImageSampleContainer & samples) const497 ::SampleFixedImageIndexes(FixedImageSampleContainer & samples) const
498 {
499   typename FixedImageSampleContainer::iterator iter;
500 
501   auto len = static_cast<SizeValueType>( m_FixedImageIndexes.size() );
502   if ( len != m_NumberOfFixedImageSamples
503        || samples.size() != m_NumberOfFixedImageSamples )
504     {
505     throw ExceptionObject(__FILE__, __LINE__,
506                           "Index list size does not match desired number of samples");
507     }
508 
509   iter = samples.begin();
510   for ( SizeValueType i = 0; i < len; i++ )
511     {
512     // Get sampled index
513     FixedImageIndexType index = m_FixedImageIndexes[i];
514     // Translate index to point
515     m_FixedImage->TransformIndexToPhysicalPoint(index, ( *iter ).point);
516 
517     // Get sampled fixed image value
518     ( *iter ).value = m_FixedImage->GetPixel(index);
519     ( *iter ).valueIndex = 0;
520 
521     ++iter;
522     }
523 }
524 
525 /**
526  * Sample the fixed image using a random walk
527  */
528 template< typename TFixedImage, typename TMovingImage >
529 void
530 ImageToImageMetric< TFixedImage, TMovingImage >
SampleFixedImageRegion(FixedImageSampleContainer & samples) const531 ::SampleFixedImageRegion(FixedImageSampleContainer & samples) const
532 {
533   if ( samples.size() != m_NumberOfFixedImageSamples )
534     {
535     throw ExceptionObject(__FILE__, __LINE__,
536                           "Sample size does not match desired number of samples");
537     }
538 
539   // Set up a random interator within the user specified fixed image region.
540   using RandomIterator = ImageRandomConstIteratorWithIndex< FixedImageType >;
541   RandomIterator randIter( m_FixedImage, GetFixedImageRegion() );
542   randIter.ReinitializeSeed(Statistics::MersenneTwisterRandomVariateGenerator::GetInstance()->GetSeed());
543   if (m_ReseedIterator)
544     {
545     randIter.ReinitializeSeed();
546     }
547   else
548     {
549     randIter.ReinitializeSeed(m_RandomSeed++);
550     }
551   typename FixedImageSampleContainer::iterator iter;
552   typename FixedImageSampleContainer::const_iterator end = samples.end();
553 
554   if ( m_FixedImageMask.IsNotNull()
555        || m_UseFixedImageSamplesIntensityThreshold )
556     {
557     InputPointType inputPoint;
558 
559     iter = samples.begin();
560     SizeValueType samplesFound = 0;
561     randIter.SetNumberOfSamples(m_NumberOfFixedImageSamples * 1000);
562     randIter.GoToBegin();
563     while ( iter != end )
564       {
565       if ( randIter.IsAtEnd() )
566         {
567         // Must be a small mask since after many random samples we don't
568         // have enough to fill the desired number.   So, we will replicate
569         // the samples we've found so far to fill-in the desired number
570         // of samples
571         SizeValueType count = 0;
572         while ( iter != end )
573           {
574           ( *iter ).point = samples[count].point;
575           ( *iter ).value = samples[count].value;
576           ( *iter ).valueIndex = 0;
577           ++count;
578           if ( count >= samplesFound )
579             {
580             count = 0;
581             }
582           ++iter;
583           }
584         break;
585         }
586 
587       // Get sampled index
588       FixedImageIndexType index = randIter.GetIndex();
589       // Check if the Index is inside the mask, translate index to point
590       m_FixedImage->TransformIndexToPhysicalPoint(index, inputPoint);
591 
592       if ( m_FixedImageMask.IsNotNull() )
593         {
594         double val;
595         if ( m_FixedImageMask->ValueAtInWorldSpace(inputPoint, val) )
596           {
597           if ( Math::AlmostEquals( val, 0.0 ) )
598             {
599             ++randIter; // jump to another random position
600             continue;
601             }
602           }
603         else
604           {
605           ++randIter; // jump to another random position
606           continue;
607           }
608         }
609 
610       if ( m_UseFixedImageSamplesIntensityThreshold
611            && randIter.Get() < m_FixedImageSamplesIntensityThreshold )
612         {
613         ++randIter;
614         continue;
615         }
616 
617       // Translate index to point
618       ( *iter ).point = inputPoint;
619       // Get sampled fixed image value
620       ( *iter ).value = randIter.Get();
621       ( *iter ).valueIndex = 0;
622 
623       ++samplesFound;
624       ++randIter;
625       ++iter;
626       }
627     }
628   else
629     {
630     randIter.SetNumberOfSamples(m_NumberOfFixedImageSamples);
631     randIter.GoToBegin();
632     for ( iter = samples.begin(); iter != end; ++iter )
633       {
634       // Get sampled index
635       FixedImageIndexType index = randIter.GetIndex();
636       // Translate index to point
637       m_FixedImage->TransformIndexToPhysicalPoint(index,
638                                                   ( *iter ).point);
639       // Get sampled fixed image value
640       ( *iter ).value = randIter.Get();
641       ( *iter ).valueIndex = 0;
642 
643       // Jump to random position
644       ++randIter;
645       }
646     }
647 }
648 
649 /**
650  * Sample the fixed image domain using all pixels in the Fixed image region
651  */
652 template< typename TFixedImage, typename TMovingImage >
653 void
654 ImageToImageMetric< TFixedImage, TMovingImage >
SampleFullFixedImageRegion(FixedImageSampleContainer & samples) const655 ::SampleFullFixedImageRegion(FixedImageSampleContainer & samples) const
656 {
657   if ( samples.size() != m_NumberOfFixedImageSamples )
658     {
659     throw ExceptionObject(__FILE__, __LINE__,
660                           "Sample size does not match desired number of samples");
661     }
662 
663   // Set up a region interator within the user specified fixed image region.
664   using RegionIterator = ImageRegionConstIteratorWithIndex< FixedImageType >;
665   RegionIterator regionIter( m_FixedImage, GetFixedImageRegion() );
666 
667   regionIter.GoToBegin();
668 
669   typename FixedImageSampleContainer::iterator iter;
670   typename FixedImageSampleContainer::const_iterator end = samples.end();
671 
672   if ( m_FixedImageMask.IsNotNull()
673        || m_UseFixedImageSamplesIntensityThreshold )
674     {
675     InputPointType inputPoint;
676 
677     // repeat until we get enough samples to fill the array
678     iter = samples.begin();
679     while ( iter != end )
680       {
681       // Get sampled index
682       FixedImageIndexType index = regionIter.GetIndex();
683       // Check if the Index is inside the mask, translate index to point
684       m_FixedImage->TransformIndexToPhysicalPoint(index, inputPoint);
685 
686       if ( m_FixedImageMask.IsNotNull() )
687         {
688         // If not inside the mask, ignore the point
689         if ( !m_FixedImageMask->IsInsideInWorldSpace(inputPoint) )
690           {
691           ++regionIter; // jump to next pixel
692           if ( regionIter.IsAtEnd() )
693             {
694             regionIter.GoToBegin();
695             }
696           continue;
697           }
698         }
699 
700       if ( m_UseFixedImageSamplesIntensityThreshold
701            && regionIter.Get() < m_FixedImageSamplesIntensityThreshold )
702         {
703         ++regionIter; // jump to next pixel
704         if ( regionIter.IsAtEnd() )
705           {
706           regionIter.GoToBegin();
707           }
708         continue;
709         }
710 
711       // Translate index to point
712       ( *iter ).point = inputPoint;
713       // Get sampled fixed image value
714       ( *iter ).value = regionIter.Get();
715       ( *iter ).valueIndex = 0;
716 
717       ++regionIter;
718       if ( regionIter.IsAtEnd() )
719         {
720         regionIter.GoToBegin();
721         }
722       ++iter;
723       }
724     }
725   else // not restricting sample throwing to a mask
726     {
727     for ( iter = samples.begin(); iter != end; ++iter )
728       {
729       // Get sampled index
730       FixedImageIndexType index = regionIter.GetIndex();
731 
732       // Translate index to point
733       m_FixedImage->TransformIndexToPhysicalPoint(index,
734                                                   ( *iter ).point);
735       // Get sampled fixed image value
736       ( *iter ).value = regionIter.Get();
737       ( *iter ).valueIndex = 0;
738 
739       ++regionIter;
740       if ( regionIter.IsAtEnd() )
741         {
742         regionIter.GoToBegin();
743         }
744       }
745     }
746 }
747 
748 /**
749  * Compute the gradient image and assign it to m_GradientImage.
750  */
751 template< typename TFixedImage, typename TMovingImage >
752 void
753 ImageToImageMetric< TFixedImage, TMovingImage >
ComputeGradient()754 ::ComputeGradient()
755 {
756   GradientImageFilterPointer gradientFilter = GradientImageFilterType::New();
757 
758   gradientFilter->SetInput(m_MovingImage);
759 
760   const typename MovingImageType::SpacingType & spacing = m_MovingImage
761                                                           ->GetSpacing();
762   double maximumSpacing = 0.0;
763   for ( unsigned int i = 0; i < MovingImageDimension; i++ )
764     {
765     if ( spacing[i] > maximumSpacing )
766       {
767       maximumSpacing = spacing[i];
768       }
769     }
770   gradientFilter->SetSigma(maximumSpacing);
771   gradientFilter->SetNormalizeAcrossScale(true);
772   gradientFilter->SetNumberOfWorkUnits(m_NumberOfWorkUnits);
773   gradientFilter->SetUseImageDirection(true);
774   gradientFilter->Update();
775 
776   m_GradientImage = gradientFilter->GetOutput();
777 }
778 
779 // Method to reinitialize the seed of the random number generator
780 template< typename TFixedImage, typename TMovingImage  >
781 void
782 ImageToImageMetric< TFixedImage, TMovingImage >
ReinitializeSeed()783 ::ReinitializeSeed()
784 {
785   m_ReseedIterator = true;
786 }
787 
788 // Method to reinitialize the seed of the random number generator
789 template< typename TFixedImage, typename TMovingImage  >
790 void
791 ImageToImageMetric< TFixedImage, TMovingImage >
ReinitializeSeed(int seed)792 ::ReinitializeSeed(int seed)
793 {
794   m_ReseedIterator = false;
795   m_RandomSeed = seed;
796 }
797 
798 /**
799  * Cache pre-transformed points, weights and indices.
800  */
801 template< typename TFixedImage, typename TMovingImage >
802 void
803 ImageToImageMetric< TFixedImage, TMovingImage >
PreComputeTransformValues()804 ::PreComputeTransformValues()
805 {
806   // Note: This code is specific to the b-spline deformable transform.
807 
808   // Unfortunately, the BSplineTransform stores a
809   // pointer to parameters passed to SetParameters(). Since
810   // we're creating a dummy set of parameters below on the
811   // stack, this can cause a crash if the transform's
812   // parameters are not later reset with a more properly
813   // scoped set of parameters. In addition, we're overwriting
814   // any previously set parameters. In order to be kinder,
815   // we'll save a pointer to the current set of parameters
816   // and restore them after we're done.
817 
818   // Note the address operator.
819   // const TransformParametersType* previousParameters = &
820   // m_Transform->GetParameters();
821 
822   // Create all zero dummy transform parameters
823   ParametersType dummyParameters(m_NumberOfParameters);
824 
825   dummyParameters.Fill(0.0);
826   m_Transform->SetParameters(dummyParameters);
827 
828   // Cycle through each sampled fixed image point
829   BSplineTransformWeightsType    weights(m_NumBSplineWeights);
830   BSplineTransformIndexArrayType indices(m_NumBSplineWeights);
831   bool                           valid;
832   MovingImagePointType           mappedPoint;
833 
834   // Declare iterators for iteration over the sample container
835   typename FixedImageSampleContainer::const_iterator fiter;
836   typename FixedImageSampleContainer::const_iterator fend =
837     m_FixedImageSamples.end();
838   SizeValueType counter = 0;
839 
840   for ( fiter = m_FixedImageSamples.begin(); fiter != fend; ++fiter, counter++ )
841     {
842     m_BSplineTransform->TransformPoint(m_FixedImageSamples[counter].point,
843                                        mappedPoint, weights, indices, valid);
844 
845     for ( SizeValueType k = 0; k < m_NumBSplineWeights; k++ )
846       {
847       m_BSplineTransformWeightsArray[counter][k] = weights[k];
848       m_BSplineTransformIndicesArray[counter][k] = indices[k];
849       }
850 
851     m_BSplinePreTransformPointsArray[counter]      = mappedPoint;
852     m_WithinBSplineSupportRegionArray[counter]     = valid;
853     }
854 
855   // Restore the previous parameters.
856   // m_Transform->SetParameters( *previousParameters );
857 }
858 
859 /**
860  * Transform a point from FixedImage domain to MovingImage domain.
861  * This function also checks if mapped point is within support region.
862  */
863 template< typename TFixedImage, typename TMovingImage >
864 void
865 ImageToImageMetric< TFixedImage, TMovingImage >
TransformPoint(unsigned int sampleNumber,MovingImagePointType & mappedPoint,bool & sampleOk,double & movingImageValue,ThreadIdType threadId) const866 ::TransformPoint(unsigned int sampleNumber,
867                  MovingImagePointType & mappedPoint,
868                  bool & sampleOk,
869                  double & movingImageValue,
870                  ThreadIdType threadId) const
871 {
872   sampleOk = true;
873   TransformType *transform;
874 
875   if ( threadId > 0 )
876     {
877     transform = this->m_ThreaderTransform[threadId - 1];
878     }
879   else
880     {
881     transform = this->m_Transform;
882     }
883 
884   if ( !m_TransformIsBSpline )
885     {
886     // Use generic transform to compute mapped position
887     mappedPoint = transform->TransformPoint(m_FixedImageSamples[sampleNumber].point);
888     sampleOk = true;
889     }
890   else
891     {
892     if ( this->m_UseCachingOfBSplineWeights )
893       {
894       sampleOk = m_WithinBSplineSupportRegionArray[sampleNumber];
895 
896       if ( sampleOk )
897         {
898         // If the transform is BSplineDeformable, we can use the precomputed
899         // weights and indices to obtained the mapped position
900         const WeightsValueType *weights =
901           m_BSplineTransformWeightsArray[sampleNumber];
902         const IndexValueType *indices =
903           m_BSplineTransformIndicesArray[sampleNumber];
904 
905         for ( unsigned int j = 0; j < FixedImageDimension; j++ )
906           {
907           mappedPoint[j] = m_BSplinePreTransformPointsArray[sampleNumber][j];
908           }
909 
910         const ParametersType &LocalParameters = m_Transform->GetParameters();
911         for ( unsigned int k = 0; k < m_NumBSplineWeights; k++ )
912           {
913           for ( unsigned int j = 0; j < FixedImageDimension; j++ )
914             {
915             mappedPoint[j] += weights[k] * LocalParameters[indices[k]
916                                                         + m_BSplineParametersOffset[j]];
917             }
918           }
919         }
920       }
921     else
922       {
923       BSplineTransformWeightsType *   weightsHelper;
924       BSplineTransformIndexArrayType *indicesHelper;
925 
926       if ( threadId > 0 )
927         {
928         weightsHelper = &( this->m_ThreaderBSplineTransformWeights[threadId - 1] );
929         indicesHelper = &( this->m_ThreaderBSplineTransformIndices[threadId - 1] );
930         }
931       else
932         {
933         weightsHelper = &( this->m_BSplineTransformWeights );
934         indicesHelper = &( this->m_BSplineTransformIndices );
935         }
936 
937       // If not caching values, we invoke the Transform to recompute the
938       // mapping of the point.
939       this->m_BSplineTransform->TransformPoint(
940         this->m_FixedImageSamples[sampleNumber].point,
941         mappedPoint, *weightsHelper, *indicesHelper, sampleOk);
942       }
943     }
944 
945   if ( sampleOk )
946     {
947     // If user provided a mask over the Moving image
948     if ( m_MovingImageMask )
949       {
950       // Check if mapped point is within the support region of the moving image
951       // mask
952       sampleOk = sampleOk &&
953         m_MovingImageMask->IsInsideInWorldSpace(mappedPoint);
954       }
955 
956     if ( m_InterpolatorIsBSpline )
957       {
958       // Check if mapped point inside image buffer
959       sampleOk = sampleOk && m_BSplineInterpolator->IsInsideBuffer(mappedPoint);
960       if ( sampleOk )
961         {
962         movingImageValue = m_BSplineInterpolator->Evaluate(mappedPoint, threadId);
963         }
964       }
965     else
966       {
967       // Check if mapped point inside image buffer
968       sampleOk = sampleOk && m_Interpolator->IsInsideBuffer(mappedPoint);
969       if ( sampleOk )
970         {
971         movingImageValue = m_Interpolator->Evaluate(mappedPoint);
972         }
973       }
974     }
975 }
976 
977 /**
978  * Transform a point from FixedImage domain to MovingImage domain.
979  * This function also checks if mapped point is within support region.
980  */
981 template< typename TFixedImage, typename TMovingImage >
982 void
983 ImageToImageMetric< TFixedImage, TMovingImage >
TransformPointWithDerivatives(unsigned int sampleNumber,MovingImagePointType & mappedPoint,bool & sampleOk,double & movingImageValue,ImageDerivativesType & movingImageGradient,ThreadIdType threadId) const984 ::TransformPointWithDerivatives(unsigned int sampleNumber,
985                                 MovingImagePointType & mappedPoint,
986                                 bool & sampleOk,
987                                 double & movingImageValue,
988                                 ImageDerivativesType & movingImageGradient,
989                                 ThreadIdType threadId) const
990 {
991   TransformType *transform;
992 
993   sampleOk = true;
994 
995   if ( threadId > 0 )
996     {
997     transform = this->m_ThreaderTransform[threadId - 1];
998     }
999   else
1000     {
1001     transform = this->m_Transform;
1002     }
1003 
1004   if ( !m_TransformIsBSpline )
1005     {
1006     // Use generic transform to compute mapped position
1007     mappedPoint = transform->TransformPoint(m_FixedImageSamples[sampleNumber].point);
1008     sampleOk = true;
1009     }
1010   else
1011     {
1012     if ( this->m_UseCachingOfBSplineWeights )
1013       {
1014       sampleOk = m_WithinBSplineSupportRegionArray[sampleNumber];
1015 
1016       if ( sampleOk )
1017         {
1018         // If the transform is BSplineDeformable, we can use the precomputed
1019         // weights and indices to obtained the mapped position
1020         const WeightsValueType *weights =
1021           m_BSplineTransformWeightsArray[sampleNumber];
1022         const IndexValueType *indices =
1023           m_BSplineTransformIndicesArray[sampleNumber];
1024 
1025         const ParametersType &Local_Parameters=this->m_Transform->GetParameters();
1026         for ( unsigned int j = 0; j < FixedImageDimension; j++ )
1027           {
1028           mappedPoint[j] = m_BSplinePreTransformPointsArray[sampleNumber][j];
1029           }
1030 
1031         for ( unsigned int k = 0; k < m_NumBSplineWeights; k++ )
1032           {
1033           for ( unsigned int j = 0; j < FixedImageDimension; j++ )
1034             {
1035             mappedPoint[j] += weights[k] * Local_Parameters[indices[k]
1036                                                         + m_BSplineParametersOffset[j]];
1037             }
1038           }
1039         }
1040       }
1041     else
1042       {
1043       BSplineTransformWeightsType *   weightsHelper;
1044       BSplineTransformIndexArrayType *indicesHelper;
1045 
1046       if ( threadId > 0 )
1047         {
1048         weightsHelper = &( this->m_ThreaderBSplineTransformWeights[threadId - 1] );
1049         indicesHelper = &( this->m_ThreaderBSplineTransformIndices[threadId - 1] );
1050         }
1051       else
1052         {
1053         weightsHelper = &( this->m_BSplineTransformWeights );
1054         indicesHelper = &( this->m_BSplineTransformIndices );
1055         }
1056 
1057       // If not caching values, we invoke the Transform to recompute the
1058       // mapping of the point.
1059       this->m_BSplineTransform->TransformPoint(
1060         this->m_FixedImageSamples[sampleNumber].point,
1061         mappedPoint, *weightsHelper, *indicesHelper, sampleOk);
1062       }
1063     }
1064 
1065   if ( sampleOk )
1066     {
1067     // If user provided a mask over the Moving image
1068     if ( m_MovingImageMask )
1069       {
1070       // Check if mapped point is within the support region of the moving image
1071       //   mask
1072       // We assume the mask and movingImage are in a common "world" space,
1073       //   so we use the mappedPoint which is in the moving image space.
1074       sampleOk = sampleOk &&
1075         m_MovingImageMask->IsInsideInWorldSpace(mappedPoint);
1076       }
1077 
1078     if ( m_InterpolatorIsBSpline )
1079       {
1080       // Check if mapped point inside image buffer
1081       sampleOk = sampleOk && m_BSplineInterpolator->IsInsideBuffer(mappedPoint);
1082       if ( sampleOk )
1083         {
1084         this->m_BSplineInterpolator->EvaluateValueAndDerivative(mappedPoint,
1085                                                                 movingImageValue,
1086                                                                 movingImageGradient,
1087                                                                 threadId);
1088         }
1089       }
1090     else
1091       {
1092       // Check if mapped point inside image buffer
1093       sampleOk = sampleOk && m_Interpolator->IsInsideBuffer(mappedPoint);
1094       if ( sampleOk )
1095         {
1096         this->ComputeImageDerivatives(mappedPoint, movingImageGradient, threadId);
1097         movingImageValue = this->m_Interpolator->Evaluate(mappedPoint);
1098         }
1099       }
1100     }
1101 }
1102 
1103 /**
1104  * Compute image derivatives using a central difference function
1105  * if we are not using a BSplineInterpolator, which includes
1106  * derivatives.
1107  */
1108 template< typename TFixedImage, typename TMovingImage >
1109 void
1110 ImageToImageMetric< TFixedImage, TMovingImage >
ComputeImageDerivatives(const MovingImagePointType & mappedPoint,ImageDerivativesType & gradient,ThreadIdType threadId) const1111 ::ComputeImageDerivatives(const MovingImagePointType & mappedPoint,
1112                           ImageDerivativesType & gradient,
1113                           ThreadIdType threadId) const
1114 {
1115   if ( m_InterpolatorIsBSpline )
1116     {
1117     // Computed moving image gradient using derivative BSpline kernel.
1118     gradient = m_BSplineInterpolator->EvaluateDerivative(mappedPoint,
1119                                                          threadId);
1120     }
1121   else
1122     {
1123     if ( m_ComputeGradient )
1124       {
1125       ContinuousIndex< double, MovingImageDimension > tempIndex;
1126       m_MovingImage->TransformPhysicalPointToContinuousIndex(mappedPoint,
1127                                                              tempIndex);
1128       MovingImageIndexType mappedIndex;
1129       mappedIndex.CopyWithRound(tempIndex);
1130       gradient = m_GradientImage->GetPixel(mappedIndex);
1131       }
1132     else
1133       {
1134       // if not using the gradient image
1135       gradient = m_DerivativeCalculator->Evaluate(mappedPoint);
1136       }
1137     }
1138 }
1139 
1140 template< typename TFixedImage, typename TMovingImage  >
1141 void
1142 ImageToImageMetric< TFixedImage, TMovingImage >
GetValueMultiThreadedInitiate() const1143 ::GetValueMultiThreadedInitiate() const
1144 {
1145   this->SynchronizeTransforms();
1146 
1147   m_Threader->SetSingleMethod( GetValueMultiThreaded, static_cast<void *>( m_ConstSelfWrapper ) );
1148   m_Threader->SingleMethodExecute();
1149 
1150   for ( ThreadIdType threadId = 0; threadId < m_NumberOfWorkUnits - 1; threadId++ )
1151     {
1152     this->m_NumberOfPixelsCounted += m_ThreaderNumberOfMovingImageSamples[threadId];
1153     }
1154 }
1155 
1156 template< typename TFixedImage, typename TMovingImage  >
1157 void
1158 ImageToImageMetric< TFixedImage, TMovingImage >
GetValueMultiThreadedPostProcessInitiate() const1159 ::GetValueMultiThreadedPostProcessInitiate() const
1160 {
1161   m_Threader->SetSingleMethod( GetValueMultiThreadedPostProcess, static_cast<void *>( m_ConstSelfWrapper ) );
1162   m_Threader->SingleMethodExecute();
1163 }
1164 
1165 /**
1166  * Get the match Measure
1167  */
1168 template< typename TFixedImage, typename TMovingImage  >
1169 ITK_THREAD_RETURN_FUNCTION_CALL_CONVENTION
1170 ImageToImageMetric< TFixedImage, TMovingImage >
GetValueMultiThreaded(void * workunitInfoAsVoid)1171 ::GetValueMultiThreaded(void *workunitInfoAsVoid)
1172 {
1173   const MultiThreaderWorkUnitInfoImageToImageMetricWrapper helper(workunitInfoAsVoid);
1174   helper.GetConstImageToImageMetricPointer()->GetValueThread(helper.GetThreadId());
1175 
1176   return ITK_THREAD_RETURN_DEFAULT_VALUE;
1177 }
1178 
1179 /**
1180  * Get the match Measure
1181  */
1182 template< typename TFixedImage, typename TMovingImage  >
1183 ITK_THREAD_RETURN_FUNCTION_CALL_CONVENTION
1184 ImageToImageMetric< TFixedImage, TMovingImage >
GetValueMultiThreadedPostProcess(void * workunitInfoAsVoid)1185 ::GetValueMultiThreadedPostProcess(void *workunitInfoAsVoid)
1186 {
1187   const MultiThreaderWorkUnitInfoImageToImageMetricWrapper helper(workunitInfoAsVoid);
1188   helper.GetConstImageToImageMetricPointer()->GetValueThreadPostProcess(helper.GetThreadId(),false);
1189 
1190   return ITK_THREAD_RETURN_DEFAULT_VALUE;
1191 }
1192 
1193 template< typename TFixedImage, typename TMovingImage  >
1194 void
1195 ImageToImageMetric< TFixedImage, TMovingImage >
GetValueThread(ThreadIdType threadId) const1196 ::GetValueThread(ThreadIdType threadId) const
1197 {
1198   // Figure out how many samples to process
1199   int chunkSize = m_NumberOfFixedImageSamples / m_NumberOfWorkUnits;
1200 
1201   // Skip to this thread's samples to process
1202   unsigned int fixedImageSample = threadId * chunkSize;
1203 
1204   if ( threadId == m_NumberOfWorkUnits - 1 )
1205     {
1206     chunkSize = m_NumberOfFixedImageSamples
1207                 - ( ( m_NumberOfWorkUnits - 1 )
1208                     * chunkSize );
1209     }
1210 
1211 
1212   if ( m_WithinThreadPreProcess )
1213     {
1214     this->GetValueThreadPreProcess(threadId, true);
1215     }
1216 
1217   // Process the samples
1218   int numSamples = 0;
1219   for ( int count = 0; count < chunkSize; ++count, ++fixedImageSample )
1220     {
1221     MovingImagePointType mappedPoint;
1222     bool                 sampleOk;
1223     double               movingImageValue;
1224     // Get moving image value
1225     this->TransformPoint(fixedImageSample, mappedPoint, sampleOk, movingImageValue,
1226                          threadId);
1227 
1228     if ( sampleOk )
1229       {
1230       // CALL USER FUNCTION
1231       if ( GetValueThreadProcessSample(threadId, fixedImageSample,
1232                                        mappedPoint, movingImageValue) )
1233         {
1234         ++numSamples;
1235         }
1236       }
1237     }
1238 
1239   if ( threadId > 0 )
1240     {
1241     m_ThreaderNumberOfMovingImageSamples[threadId - 1] = numSamples;
1242     }
1243   else
1244     {
1245     m_NumberOfPixelsCounted = numSamples;
1246     }
1247 
1248   if ( m_WithinThreadPostProcess )
1249     {
1250     this->GetValueThreadPostProcess(threadId, true);
1251     }
1252 }
1253 
1254 template< typename TFixedImage, typename TMovingImage  >
1255 void
1256 ImageToImageMetric< TFixedImage, TMovingImage >
GetValueAndDerivativeMultiThreadedInitiate() const1257 ::GetValueAndDerivativeMultiThreadedInitiate() const
1258 {
1259   this->SynchronizeTransforms();
1260 
1261   m_Threader->SetSingleMethod( GetValueAndDerivativeMultiThreaded, static_cast<void *>( m_ConstSelfWrapper ) );
1262   m_Threader->SingleMethodExecute();
1263 
1264   for ( ThreadIdType threadId = 0; threadId < m_NumberOfWorkUnits - 1; threadId++ )
1265     {
1266     this->m_NumberOfPixelsCounted += m_ThreaderNumberOfMovingImageSamples[threadId];
1267     }
1268 }
1269 
1270 template< typename TFixedImage, typename TMovingImage  >
1271 void
1272 ImageToImageMetric< TFixedImage, TMovingImage >
GetValueAndDerivativeMultiThreadedPostProcessInitiate() const1273 ::GetValueAndDerivativeMultiThreadedPostProcessInitiate() const
1274 {
1275   m_Threader->SetSingleMethod( GetValueAndDerivativeMultiThreadedPostProcess, static_cast<void *>(m_ConstSelfWrapper ) );
1276   m_Threader->SingleMethodExecute();
1277 }
1278 
1279 /**
1280  * Get the match Measure
1281  */
1282 template< typename TFixedImage, typename TMovingImage  >
1283 ITK_THREAD_RETURN_FUNCTION_CALL_CONVENTION
1284 ImageToImageMetric< TFixedImage, TMovingImage >
GetValueAndDerivativeMultiThreaded(void * workunitInfoAsVoid)1285 ::GetValueAndDerivativeMultiThreaded(void *workunitInfoAsVoid)
1286 {
1287   const MultiThreaderWorkUnitInfoImageToImageMetricWrapper helper(workunitInfoAsVoid);
1288   helper.GetConstImageToImageMetricPointer()->GetValueAndDerivativeThread(helper.GetThreadId());
1289 
1290   return ITK_THREAD_RETURN_DEFAULT_VALUE;
1291 }
1292 
1293 /**
1294  * Get the match Measure
1295  */
1296 template< typename TFixedImage, typename TMovingImage  >
1297 ITK_THREAD_RETURN_FUNCTION_CALL_CONVENTION
1298 ImageToImageMetric< TFixedImage, TMovingImage >
GetValueAndDerivativeMultiThreadedPostProcess(void * workunitInfoAsVoid)1299 ::GetValueAndDerivativeMultiThreadedPostProcess(void *workunitInfoAsVoid)
1300 {
1301   const MultiThreaderWorkUnitInfoImageToImageMetricWrapper helper(workunitInfoAsVoid);
1302   helper.GetConstImageToImageMetricPointer()->GetValueAndDerivativeThreadPostProcess(helper.GetThreadId(),false);
1303 
1304   return ITK_THREAD_RETURN_DEFAULT_VALUE;
1305 }
1306 
1307 template< typename TFixedImage, typename TMovingImage  >
1308 void
1309 ImageToImageMetric< TFixedImage, TMovingImage >
GetValueAndDerivativeThread(ThreadIdType threadId) const1310 ::GetValueAndDerivativeThread(ThreadIdType threadId) const
1311 {
1312   // Figure out how many samples to process
1313   int chunkSize = m_NumberOfFixedImageSamples / m_NumberOfWorkUnits;
1314 
1315   // Skip to this thread's samples to process
1316   unsigned int fixedImageSample = threadId * chunkSize;
1317 
1318   if ( threadId == m_NumberOfWorkUnits - 1 )
1319     {
1320     chunkSize = m_NumberOfFixedImageSamples
1321                 - ( ( m_NumberOfWorkUnits - 1 )
1322                     * chunkSize );
1323     }
1324 
1325   int numSamples = 0;
1326 
1327   if ( m_WithinThreadPreProcess )
1328     {
1329     this->GetValueAndDerivativeThreadPreProcess(threadId, true);
1330     }
1331 
1332   // Process the samples
1333   MovingImagePointType mappedPoint;
1334   bool                 sampleOk;
1335   double               movingImageValue;
1336   ImageDerivativesType movingImageGradientValue;
1337   for ( int count = 0; count < chunkSize; ++count, ++fixedImageSample )
1338     {
1339     // Get moving image value
1340     TransformPointWithDerivatives(fixedImageSample, mappedPoint, sampleOk,
1341                                   movingImageValue, movingImageGradientValue,
1342                                   threadId);
1343 
1344     if ( sampleOk )
1345       {
1346       // CALL USER FUNCTION
1347       if ( this->GetValueAndDerivativeThreadProcessSample(threadId,
1348                                                           fixedImageSample,
1349                                                           mappedPoint,
1350                                                           movingImageValue,
1351                                                           movingImageGradientValue) )
1352         {
1353         ++numSamples;
1354         }
1355       }
1356     }
1357 
1358   if ( threadId > 0 )
1359     {
1360     m_ThreaderNumberOfMovingImageSamples[threadId - 1] = numSamples;
1361     }
1362   else
1363     {
1364     m_NumberOfPixelsCounted = numSamples;
1365     }
1366 
1367   if ( m_WithinThreadPostProcess )
1368     {
1369     this->GetValueAndDerivativeThreadPostProcess(threadId, true);
1370     }
1371 }
1372 
1373 /**
1374  * PrintSelf
1375  */
1376 template< typename TFixedImage, typename TMovingImage >
1377 void
1378 ImageToImageMetric< TFixedImage, TMovingImage >
PrintSelf(std::ostream & os,Indent indent) const1379 ::PrintSelf(std::ostream & os, Indent indent) const
1380 {
1381   Superclass::PrintSelf(os, indent);
1382 
1383   os << indent << "NumberOfFixedImageSamples: ";
1384   os << m_NumberOfFixedImageSamples << std::endl;
1385 
1386   os << indent << "FixedImageSamplesIntensityThreshold: "
1387      << static_cast< typename NumericTraits< FixedImagePixelType >::PrintType >( m_FixedImageSamplesIntensityThreshold )
1388      << std::endl;
1389 
1390   os << indent << "UseFixedImageSamplesIntensityThreshold: ";
1391   os << m_UseFixedImageSamplesIntensityThreshold << std::endl;
1392 
1393   if ( m_UseFixedImageIndexes )
1394     {
1395     os << indent << "Use Fixed Image Indexes: True" << std::endl;
1396     os << indent << "Number of Fixed Image Indexes = "
1397        << m_FixedImageIndexes.size() << std::endl;
1398     }
1399   else
1400     {
1401     os << indent << "Use Fixed Image Indexes: False" << std::endl;
1402     }
1403 
1404   if ( m_UseSequentialSampling )
1405     {
1406     os << indent << "Use Sequential Sampling: True" << std::endl;
1407     }
1408   else
1409     {
1410     os << indent << "Use Sequential Sampling: False" << std::endl;
1411     }
1412 
1413   os << indent << "UseAllPixels: ";
1414   os << m_UseAllPixels << std::endl;
1415 
1416   os << indent << "ReseedIterator: " << m_ReseedIterator << std::endl;
1417   os << indent << "RandomSeed: " << m_RandomSeed << std::endl;
1418 
1419   os << indent << "Threader: " << m_Threader << std::endl;
1420   os << indent << "Number of Work units: " << m_NumberOfWorkUnits << std::endl;
1421   os << indent << "ThreaderParameter: " << std::endl;
1422   os << indent << "ThreaderNumberOfMovingImageSamples: " << std::endl;
1423   if ( m_ThreaderNumberOfMovingImageSamples )
1424     {
1425     for ( ThreadIdType i = 0; i < m_NumberOfWorkUnits - 1; i++ )
1426       {
1427       os << "  Thread[" << i << "]= " << (unsigned int)m_ThreaderNumberOfMovingImageSamples[i] << std::endl;
1428       }
1429     }
1430 
1431   os << indent << "ComputeGradient: "
1432      << static_cast< typename NumericTraits< bool >::PrintType >( m_ComputeGradient )
1433      << std::endl;
1434   os << indent << "Moving Image: " << m_MovingImage.GetPointer()  << std::endl;
1435   os << indent << "Fixed  Image: " << m_FixedImage.GetPointer()   << std::endl;
1436   os << indent << "Gradient Image: " << m_GradientImage.GetPointer()
1437      << std::endl;
1438   os << indent << "Transform:    " << m_Transform.GetPointer()    << std::endl;
1439   os << indent << "Interpolator: " << m_Interpolator.GetPointer() << std::endl;
1440   os << indent << "FixedImageRegion: " << m_FixedImageRegion << std::endl;
1441   os << indent << "Moving Image Mask: " << m_MovingImageMask.GetPointer()
1442      << std::endl;
1443   os << indent << "Fixed Image Mask: " << m_FixedImageMask.GetPointer()
1444      << std::endl;
1445   os << indent << "Number of Moving Image Samples: " << m_NumberOfPixelsCounted
1446      << std::endl;
1447 
1448   os << indent << "UseCachingOfBSplineWeights: ";
1449   os << this->m_UseCachingOfBSplineWeights << std::endl;
1450 }
1451 
1452 /** This method can be const because we are not altering the m_ThreaderTransform
1453  *  pointer. We are altering the object that m_ThreaderTransform[idx] points at.
1454  *  This is allowed under C++ const rules.
1455  */
1456 template< typename TFixedImage, typename TMovingImage >
1457 void
1458 ImageToImageMetric< TFixedImage, TMovingImage >
SynchronizeTransforms() const1459 ::SynchronizeTransforms() const
1460 {
1461   for ( ThreadIdType threadId = 0; threadId < m_NumberOfWorkUnits - 1; threadId++ )
1462     {
1463     /** Set the fixed parameters first. Some transforms have parameters which depend on
1464         the values of the fixed parameters. For instance, the BSplineTransform
1465         checks the grid size (part of the fixed parameters) before setting the parameters. */
1466     this->m_ThreaderTransform[threadId]->SetFixedParameters( this->m_Transform->GetFixedParameters() );
1467     this->m_ThreaderTransform[threadId]->SetParameters( this->m_Transform->GetParameters() );
1468     }
1469 }
1470 } // end namespace itk
1471 
1472 #endif
1473