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