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 itkANTSNeighborhoodCorrelationImageToImageMetricv4GetValueAndDerivativeThreader_hxx
19 #define itkANTSNeighborhoodCorrelationImageToImageMetricv4GetValueAndDerivativeThreader_hxx
20 
21 #include "itkANTSNeighborhoodCorrelationImageToImageMetricv4GetValueAndDerivativeThreader.h"
22 
23 namespace itk
24 {
25 
26 template < typename TDomainPartitioner, typename TImageToImageMetric, typename TNeighborhoodCorrelationMetric >
27 void
28 ANTSNeighborhoodCorrelationImageToImageMetricv4GetValueAndDerivativeThreader< TDomainPartitioner, TImageToImageMetric, TNeighborhoodCorrelationMetric >
ThreadedExecution_impl(IdentityHelper<ThreadedImageRegionPartitioner<TImageToImageMetric::VirtualImageDimension>> itkNotUsed (self),const DomainType & virtualImageSubRegion,const ThreadIdType threadId)29 ::ThreadedExecution_impl(
30     IdentityHelper<ThreadedImageRegionPartitioner<TImageToImageMetric::VirtualImageDimension> > itkNotUsed(self),
31     const DomainType& virtualImageSubRegion,
32     const ThreadIdType threadId )
33 {
34   /* Store the casted pointer to avoid dynamic casting in tight loops. */
35   this->m_ANTSAssociate = dynamic_cast< TNeighborhoodCorrelationMetric * >( this->m_Associate );
36   if( this->m_ANTSAssociate == nullptr )
37     {
38     itkExceptionMacro("Dynamic casting of associate pointer failed.");
39     }
40 
41   VirtualPointType     virtualPoint;
42   MeasureType          metricValueResult = NumericTraits< MeasureType >::ZeroValue();
43   MeasureType          metricValueSum = NumericTraits< MeasureType >::ZeroValue();
44   bool                 pointIsValid;
45   ScanIteratorType     scanIt;
46   ScanParametersType   scanParameters;
47   ScanMemType          scanMem;
48 
49   DerivativeType & localDerivativeResult = this->m_GetValueAndDerivativePerThreadVariables[threadId].LocalDerivatives;
50 
51   /* Create an iterator over the virtual sub region */
52   // this->m_ANTSAssociate->InitializeScanning( virtualImageSubRegion, scanIt, scanMem, scanParameters );
53   this->InitializeScanning( virtualImageSubRegion, scanIt, scanMem, scanParameters );
54 
55   /* Iterate over the sub region */
56   scanIt.GoToBegin();
57   while (!scanIt.IsAtEnd())
58     {
59     /* Get the virtual point */
60     this->m_ANTSAssociate->TransformVirtualIndexToPhysicalPoint( scanIt.GetIndex(), virtualPoint );
61 
62     /* Call the user method in derived classes to do the specific
63      * calculations for value and derivative. */
64     try
65       {
66       this->UpdateQueues(scanIt, scanMem, scanParameters, threadId);
67       pointIsValid = this->ComputeInformationFromQueues(scanIt, scanMem, scanParameters, threadId);
68       if( pointIsValid )
69         {
70         this->ComputeMovingTransformDerivative(scanIt, scanMem, scanParameters, localDerivativeResult, metricValueResult, threadId );
71         }
72       }
73     catch (ExceptionObject & exc)
74       {
75       //NOTE: there must be a cleaner way to do this:
76       std::string msg("Caught exception: \n");
77       msg += exc.what();
78       ExceptionObject err(__FILE__, __LINE__, msg);
79       throw err;
80       }
81 
82     /* Assign the results */
83     if ( pointIsValid )
84       {
85       this->m_GetValueAndDerivativePerThreadVariables[threadId].NumberOfValidPoints++;
86       metricValueSum -= metricValueResult;
87       /* Store the result. This depends on what type of
88        * transform is being used. */
89       if( this->GetComputeDerivative() )
90         {
91         this->StorePointDerivativeResult( scanIt.GetIndex(), threadId );
92         }
93       }
94 
95     //next index
96     ++scanIt;
97     }
98 
99   /* Store metric value result for this thread. */
100   this->m_GetValueAndDerivativePerThreadVariables[threadId].Measure = metricValueSum;
101 }
102 
103 template < typename TDomainPartitioner, typename TImageToImageMetric, typename TNeighborhoodCorrelationMetric >
104 template < typename T >
105 void
106 ANTSNeighborhoodCorrelationImageToImageMetricv4GetValueAndDerivativeThreader< TDomainPartitioner, TImageToImageMetric, TNeighborhoodCorrelationMetric >
ThreadedExecution_impl(IdentityHelper<T> itkNotUsed (self),const DomainType & domain,const ThreadIdType threadId)107 ::ThreadedExecution_impl(
108     IdentityHelper<T> itkNotUsed(self),
109     const DomainType& domain,
110     const ThreadIdType threadId )
111 {
112     /* call base method */
113     /* Store the casted pointer to avoid dynamic casting in tight loops. */
114     this->m_ANTSAssociate = dynamic_cast< TNeighborhoodCorrelationMetric * >( this->m_Associate );
115     if( this->m_ANTSAssociate == nullptr )
116       {
117       itkExceptionMacro("Dynamic casting of associate pointer failed.");
118       }
119 
120     Superclass::ThreadedExecution(domain, threadId);
121 }
122 
123 template < typename TDomainPartitioner, typename TImageToImageMetric, typename TNeighborhoodCorrelationMetric >
124 void
125 ANTSNeighborhoodCorrelationImageToImageMetricv4GetValueAndDerivativeThreader< TDomainPartitioner, TImageToImageMetric, TNeighborhoodCorrelationMetric >
UpdateQueuesAtBeginningOfLine(const ScanIteratorType & scanIt,ScanMemType & scanMem,const ScanParametersType & scanParameters,const ThreadIdType) const126 ::UpdateQueuesAtBeginningOfLine( const ScanIteratorType &scanIt, ScanMemType &scanMem, const ScanParametersType &scanParameters, const ThreadIdType ) const
127  {
128    const SizeValueType numberOfFillZero = scanParameters.numberOfFillZero;
129    const SizeValueType hoodlen = scanParameters.windowLength;
130 
131    InternalComputationValueType zero = NumericTraits<InternalComputationValueType>::ZeroValue();
132    scanMem.QsumFixed2 = SumQueueType(numberOfFillZero, zero);
133    scanMem.QsumMoving2 = SumQueueType(numberOfFillZero, zero);
134    scanMem.QsumFixed = SumQueueType(numberOfFillZero, zero);
135    scanMem.QsumMoving = SumQueueType(numberOfFillZero, zero);
136    scanMem.QsumFixedMoving = SumQueueType(numberOfFillZero, zero);
137    scanMem.Qcount = SumQueueType(numberOfFillZero, zero);
138 
139    using LocalRealType = InternalComputationValueType;
140 
141    // Now add the rest of the values from each hyperplane
142    SizeValueType diameter = 2 * scanParameters.radius[0];
143 
144    const LocalRealType localZero = NumericTraits<LocalRealType>::ZeroValue();
145    for (SizeValueType i = numberOfFillZero; i < ( diameter + NumericTraits<SizeValueType>::OneValue() ); i++)
146      {
147      LocalRealType sumFixed2      = localZero;
148      LocalRealType sumMoving2     = localZero;
149      LocalRealType sumFixed       = localZero;
150      LocalRealType sumMoving      = localZero;
151      LocalRealType sumFixedMoving = localZero;
152      LocalRealType count = localZero;
153 
154      for ( SizeValueType indct = i; indct < hoodlen; indct += ( diameter + NumericTraits<SizeValueType>::OneValue() ) )
155        {
156        typename ScanIteratorType::OffsetType internalIndex, offset;
157        bool isInBounds = scanIt.IndexInBounds( indct, internalIndex, offset );
158        if (!isInBounds)
159          {
160          // std::cout << "DEBUG: error" << std::endl;
161          continue;
162          }
163 
164        typename VirtualImageType::IndexType index = scanIt.GetIndex(indct);
165 
166        VirtualPointType        virtualPoint;
167        FixedImagePointType     mappedFixedPoint;
168        FixedImagePixelType     fixedImageValue;
169        MovingImagePointType    mappedMovingPoint;
170        MovingImagePixelType    movingImageValue;
171        bool pointIsValid;
172 
173        this->m_ANTSAssociate->TransformVirtualIndexToPhysicalPoint(index, virtualPoint);
174 
175        try
176          {
177          pointIsValid = this->m_ANTSAssociate->TransformAndEvaluateFixedPoint( virtualPoint, mappedFixedPoint, fixedImageValue );
178          if ( pointIsValid )
179            {
180            pointIsValid = this->m_ANTSAssociate->TransformAndEvaluateMovingPoint( virtualPoint, mappedMovingPoint, movingImageValue );
181            }
182          }
183        catch (ExceptionObject & exc)
184          {
185          //NOTE: there must be a cleaner way to do this:
186          std::string msg("Caught exception: \n");
187          msg += exc.what();
188          ExceptionObject err(__FILE__, __LINE__, msg);
189          throw err;
190          }
191 
192 
193        if ( pointIsValid )
194          {
195          sumFixed2 += fixedImageValue  * fixedImageValue;
196          sumMoving2 += movingImageValue * movingImageValue;
197          sumFixed += fixedImageValue;
198          sumMoving += movingImageValue;
199          sumFixedMoving += fixedImageValue * movingImageValue;
200          count += NumericTraits<LocalRealType>::OneValue();
201          }
202        }//for indct
203 
204      scanMem.QsumFixed2.push_back(sumFixed2);
205      scanMem.QsumMoving2.push_back(sumMoving2);
206      scanMem.QsumFixed.push_back(sumFixed);
207      scanMem.QsumMoving.push_back(sumMoving);
208      scanMem.QsumFixedMoving.push_back(sumFixedMoving);
209      scanMem.Qcount.push_back(count);
210      }
211  }
212 
213 template < typename TDomainPartitioner, typename TImageToImageMetric, typename TNeighborhoodCorrelationMetric >
214  void
215  ANTSNeighborhoodCorrelationImageToImageMetricv4GetValueAndDerivativeThreader< TDomainPartitioner, TImageToImageMetric, TNeighborhoodCorrelationMetric >
UpdateQueuesToNextScanWindow(const ScanIteratorType & scanIt,ScanMemType & scanMem,const ScanParametersType & scanParameters,const ThreadIdType) const216 ::UpdateQueuesToNextScanWindow( const ScanIteratorType &scanIt, ScanMemType &scanMem, const ScanParametersType &scanParameters, const ThreadIdType ) const
217 {
218  const SizeValueType hoodlen = scanParameters.windowLength;
219 
220  using LocalRealType = InternalComputationValueType;
221 
222  const LocalRealType localZero = NumericTraits<LocalRealType>::ZeroValue();
223 
224  LocalRealType sumFixed2      = localZero;
225  LocalRealType sumMoving2     = localZero;
226  LocalRealType sumFixed       = localZero;
227  LocalRealType sumMoving      = localZero;
228  LocalRealType sumFixedMoving = localZero;
229  LocalRealType count          = localZero;
230 
231  SizeValueType diameter = 2 * scanParameters.radius[0];
232 
233  for ( SizeValueType indct = diameter; indct < hoodlen; indct += (diameter + NumericTraits<SizeValueType>::OneValue()))
234    {
235    typename ScanIteratorType::OffsetType internalIndex, offset;
236    bool isInBounds = scanIt.IndexInBounds( indct, internalIndex, offset );
237 
238    if (!isInBounds)
239    {
240    continue;
241    }
242 
243    typename VirtualImageType::IndexType index = scanIt.GetIndex(indct);
244 
245    VirtualPointType virtualPoint;
246    FixedImagePointType mappedFixedPoint;
247    FixedImagePixelType fixedImageValue;
248    MovingImagePointType mappedMovingPoint;
249    MovingImagePixelType movingImageValue;
250    bool pointIsValid;
251 
252    this->m_ANTSAssociate->TransformVirtualIndexToPhysicalPoint(index, virtualPoint);
253    try
254      {
255      pointIsValid = this->m_ANTSAssociate->TransformAndEvaluateFixedPoint( virtualPoint, mappedFixedPoint, fixedImageValue );
256      if (pointIsValid)
257        {
258        pointIsValid = this->m_ANTSAssociate->TransformAndEvaluateMovingPoint( virtualPoint, mappedMovingPoint, movingImageValue );
259        }
260      }
261    catch (ExceptionObject & exc)
262      {
263      //NOTE: there must be a cleaner way to do this:
264      std::string msg("Caught exception: \n");
265      msg += exc.what();
266      ExceptionObject err(__FILE__, __LINE__, msg);
267      throw err;
268      }
269    if ( pointIsValid )
270      {
271      sumFixed2 += fixedImageValue  * fixedImageValue;
272      sumMoving2 += movingImageValue * movingImageValue;
273      sumFixed += fixedImageValue;
274      sumMoving += movingImageValue;
275      sumFixedMoving += fixedImageValue * movingImageValue;
276      count += NumericTraits<LocalRealType>::OneValue();
277      }
278    }
279    scanMem.QsumFixed2.push_back(sumFixed2);
280    scanMem.QsumMoving2.push_back(sumMoving2);
281    scanMem.QsumFixed.push_back(sumFixed);
282    scanMem.QsumMoving.push_back(sumMoving);
283    scanMem.QsumFixedMoving.push_back(sumFixedMoving);
284    scanMem.Qcount.push_back(count);
285 
286    scanMem.QsumFixed2.pop_front();
287    scanMem.QsumMoving2.pop_front();
288    scanMem.QsumFixed.pop_front();
289    scanMem.QsumMoving.pop_front();
290    scanMem.QsumFixedMoving.pop_front();
291    scanMem.Qcount.pop_front();
292 }
293 
294 template < typename TDomainPartitioner, typename TImageToImageMetric, typename TNeighborhoodCorrelationMetric >
295 void
296 ANTSNeighborhoodCorrelationImageToImageMetricv4GetValueAndDerivativeThreader< TDomainPartitioner, TImageToImageMetric, TNeighborhoodCorrelationMetric >
InitializeScanning(const ImageRegionType & scanRegion,ScanIteratorType & scanIt,ScanMemType & scanMem,ScanParametersType & scanParameters) const297 ::InitializeScanning( const ImageRegionType &scanRegion, ScanIteratorType &scanIt,
298                       ScanMemType & scanMem, ScanParametersType &scanParameters ) const
299 {
300   scanParameters.scanRegion   = scanRegion;
301   scanParameters.fixedImage   = this->m_ANTSAssociate->m_FixedImage;
302   scanParameters.movingImage  = this->m_ANTSAssociate->m_MovingImage;
303   scanParameters.virtualImage = this->m_ANTSAssociate->GetVirtualImage();
304   scanParameters.radius       = this->m_ANTSAssociate->GetRadius();
305 
306   OffsetValueType numberOfFillZero = this->m_ANTSAssociate->GetVirtualRegion().GetIndex(0)
307       - (scanRegion.GetIndex(0) - scanParameters.radius[0]);
308 
309   if (numberOfFillZero < NumericTraits<OffsetValueType>::ZeroValue())
310     {
311     numberOfFillZero = NumericTraits<OffsetValueType>::ZeroValue();
312     }
313 
314   scanParameters.numberOfFillZero = numberOfFillZero;
315 
316   scanIt = ScanIteratorType(scanParameters.radius, scanParameters.virtualImage, scanRegion);
317   scanParameters.windowLength = scanIt.Size();
318   scanParameters.scanRegionBeginIndexDim0 = scanIt.GetBeginIndex()[0];
319 
320   scanMem.fixedA = NumericTraits< QueueRealType >::ZeroValue();
321   scanMem.movingA = NumericTraits< QueueRealType >::ZeroValue();
322   scanMem.sFixedMoving = NumericTraits< QueueRealType >::ZeroValue();
323   scanMem.sFixedFixed = NumericTraits< QueueRealType >::ZeroValue();
324   scanMem.sMovingMoving = NumericTraits< QueueRealType >::ZeroValue();
325 
326   scanMem.fixedImageGradient.Fill(0.0);
327   scanMem.movingImageGradient.Fill(0.0);
328   scanMem.mappedMovingPoint.Fill(0.0);
329 }
330 
331 template < typename TDomainPartitioner, typename TImageToImageMetric, typename TNeighborhoodCorrelationMetric >
332 void
333 ANTSNeighborhoodCorrelationImageToImageMetricv4GetValueAndDerivativeThreader< TDomainPartitioner, TImageToImageMetric, TNeighborhoodCorrelationMetric >
UpdateQueues(const ScanIteratorType & scanIt,ScanMemType & scanMem,const ScanParametersType & scanParameters,const ThreadIdType threadId) const334 ::UpdateQueues( const ScanIteratorType &scanIt, ScanMemType &scanMem,
335                 const ScanParametersType &scanParameters, const ThreadIdType threadId) const
336 {
337   if (scanIt.GetIndex()[0] == scanParameters.scanRegionBeginIndexDim0 )
338     {
339     this->UpdateQueuesAtBeginningOfLine(scanIt, scanMem, scanParameters, threadId);
340     }
341   else
342     {
343     this->UpdateQueuesToNextScanWindow(scanIt, scanMem, scanParameters, threadId);
344     }
345 }
346 
347 
348 template < typename TDomainPartitioner, typename TImageToImageMetric, typename TNeighborhoodCorrelationMetric >
349 bool
350 ANTSNeighborhoodCorrelationImageToImageMetricv4GetValueAndDerivativeThreader< TDomainPartitioner, TImageToImageMetric, TNeighborhoodCorrelationMetric >
ComputeInformationFromQueues(const ScanIteratorType & scanIt,ScanMemType & scanMem,const ScanParametersType &,const ThreadIdType) const351 ::ComputeInformationFromQueues( const ScanIteratorType &scanIt, ScanMemType &scanMem, const ScanParametersType &, const ThreadIdType ) const
352 {
353  using LocalRealType = InternalComputationValueType;
354 
355  const LocalRealType localZero = NumericTraits<LocalRealType>::ZeroValue();
356 
357  LocalRealType count = localZero;
358 
359  auto itcount = scanMem.Qcount.begin();
360  while (itcount != scanMem.Qcount.end())
361    {
362    count += *itcount;
363    ++itcount;
364    }
365 
366  if (count <= localZero)
367    {
368    // no points available in the queue, perhaps out of image region
369    return false;
370    }
371 
372  // If there are values, we need to calculate the different quantities
373  LocalRealType sumFixed2      = localZero;
374  LocalRealType sumMoving2     = localZero;
375  LocalRealType sumFixed       = localZero;
376  LocalRealType sumMoving      = localZero;
377  LocalRealType sumFixedMoving = localZero;
378  auto itFixed2      = scanMem.QsumFixed2.begin();
379  auto itMoving2     = scanMem.QsumMoving2.begin();
380  auto itFixed       = scanMem.QsumFixed.begin();
381  auto itMoving      = scanMem.QsumMoving.begin();
382  auto itFixedMoving = scanMem.QsumFixedMoving.begin();
383 
384  while (itFixed2 != scanMem.QsumFixed2.end())
385    {
386    sumFixed2 += *itFixed2;
387    sumMoving2 += *itMoving2;
388    sumFixed += *itFixed;
389    sumMoving += *itMoving;
390    sumFixedMoving += *itFixedMoving;
391 
392    ++itFixed2;
393    ++itMoving2;
394    ++itFixed;
395    ++itMoving;
396    ++itFixedMoving;
397    }
398 
399  LocalRealType fixedMean  = sumFixed  / count;
400  LocalRealType movingMean = sumMoving / count;
401 
402  LocalRealType sFixedFixed   = sumFixed2 - fixedMean * sumFixed - fixedMean * sumFixed + count * fixedMean * fixedMean;
403  LocalRealType sMovingMoving = sumMoving2 - movingMean * sumMoving - movingMean * sumMoving + count * movingMean * movingMean;
404  LocalRealType sFixedMoving  = sumFixedMoving - movingMean * sumFixed - fixedMean * sumMoving + count * movingMean * fixedMean;
405 
406  typename VirtualImageType::IndexType oindex = scanIt.GetIndex();
407 
408  VirtualPointType        virtualPoint;
409  FixedImagePointType     mappedFixedPoint;
410  FixedImagePixelType     fixedImageValue;
411  FixedImageGradientType  fixedImageGradient;
412  MovingImagePointType    mappedMovingPoint;
413  MovingImagePixelType    movingImageValue;
414  MovingImageGradientType movingImageGradient;
415  bool pointIsValid;
416 
417  this->m_ANTSAssociate->TransformVirtualIndexToPhysicalPoint(oindex, virtualPoint);
418 
419  try
420    {
421    pointIsValid = this->m_ANTSAssociate->TransformAndEvaluateFixedPoint( virtualPoint, mappedFixedPoint, fixedImageValue );
422    if ( pointIsValid )
423      {
424      pointIsValid = this->m_ANTSAssociate->TransformAndEvaluateMovingPoint( virtualPoint, mappedMovingPoint, movingImageValue );
425      if( pointIsValid && this->m_ANTSAssociate->GetComputeDerivative() )
426        {
427        if( this->m_ANTSAssociate->GetGradientSourceIncludesFixed() )
428          {
429          this->m_ANTSAssociate->ComputeFixedImageGradientAtPoint( mappedFixedPoint, fixedImageGradient );
430          }
431        if( this->m_ANTSAssociate->GetGradientSourceIncludesMoving() )
432          {
433          this->m_ANTSAssociate->ComputeMovingImageGradientAtPoint( mappedMovingPoint, movingImageGradient );
434          }
435        }
436      }
437    }
438  catch (ExceptionObject & exc)
439    {
440    //NOTE: there must be a cleaner way to do this:
441    std::string msg("Caught exception: \n");
442    msg += exc.what();
443    ExceptionObject err(__FILE__, __LINE__, msg);
444    throw err;
445    }
446 
447  if ( pointIsValid )
448    {
449    scanMem.fixedA        = fixedImageValue  - fixedMean;
450    scanMem.movingA       = movingImageValue - movingMean;
451    scanMem.sFixedMoving  = sFixedMoving;
452    scanMem.sFixedFixed   = sFixedFixed;
453    scanMem.sMovingMoving = sMovingMoving;
454 
455    scanMem.fixedImageGradient  = fixedImageGradient;
456    scanMem.movingImageGradient = movingImageGradient;
457 
458    scanMem.mappedFixedPoint  = mappedFixedPoint;
459    scanMem.mappedMovingPoint = mappedMovingPoint;
460    scanMem.virtualPoint      = virtualPoint;
461    }
462 
463  return pointIsValid;
464 }
465 
466 template < typename TDomainPartitioner, typename TImageToImageMetric, typename TNeighborhoodCorrelationMetric >
467 void
468 ANTSNeighborhoodCorrelationImageToImageMetricv4GetValueAndDerivativeThreader< TDomainPartitioner, TImageToImageMetric, TNeighborhoodCorrelationMetric >
ComputeMovingTransformDerivative(const ScanIteratorType &,ScanMemType & scanMem,const ScanParametersType &,DerivativeType & deriv,MeasureType & localCC,const ThreadIdType threadId) const469 ::ComputeMovingTransformDerivative( const ScanIteratorType &, ScanMemType &scanMem, const ScanParametersType &, DerivativeType &deriv, MeasureType &localCC, const ThreadIdType threadId) const
470 {
471   MovingImageGradientType derivWRTImage;
472   localCC = NumericTraits<MeasureType>::OneValue();
473 
474   using LocalRealType = InternalComputationValueType;
475 
476   LocalRealType sFixedFixed   = scanMem.sFixedFixed;
477   LocalRealType sMovingMoving = scanMem.sMovingMoving;
478   LocalRealType sFixedMoving  = scanMem.sFixedMoving;
479   LocalRealType fixedI        = scanMem.fixedA;
480   LocalRealType movingI       = scanMem.movingA;
481 
482   LocalRealType sFixedFixed_sMovingMoving = sFixedFixed * sMovingMoving;
483 
484   if ( fabs(sFixedFixed_sMovingMoving) > NumericTraits< LocalRealType >::epsilon() )
485     {
486     localCC = sFixedMoving * sFixedMoving / (sFixedFixed_sMovingMoving);
487     }
488 
489   if( this->m_ANTSAssociate->GetComputeDerivative() )
490     {
491     const MovingImageGradientType movingImageGradient = scanMem.movingImageGradient;
492 
493     if ( ! (sFixedFixed > NumericTraits<LocalRealType>::epsilon() && sMovingMoving > NumericTraits<LocalRealType>::epsilon() ) )
494       {
495       deriv.Fill( NumericTraits<DerivativeValueType>::ZeroValue() );
496       return;
497       }
498 
499     for (ImageDimensionType qq = 0; qq < TImageToImageMetric::VirtualImageDimension; qq++)
500       {
501       derivWRTImage[qq] = 2.0 * sFixedMoving / (sFixedFixed_sMovingMoving) * (fixedI - sFixedMoving / sMovingMoving * movingI) * movingImageGradient[qq];
502       }
503 
504     /* Use a pre-allocated jacobian object for efficiency */
505     using JacobianReferenceType = JacobianType &;
506     JacobianReferenceType jacobian = this->m_GetValueAndDerivativePerThreadVariables[threadId].MovingTransformJacobian;
507     JacobianReferenceType jacobianPositional = this->m_GetValueAndDerivativePerThreadVariables[threadId].MovingTransformJacobianPositional;
508 
509     /** For dense transforms, this returns identity */
510     this->m_Associate->GetMovingTransform()->
511       ComputeJacobianWithRespectToParametersCachedTemporaries(scanMem.virtualPoint,
512                                                               jacobian,
513                                                               jacobianPositional);
514 
515     NumberOfParametersType numberOfLocalParameters = this->m_Associate->GetMovingTransform()->GetNumberOfLocalParameters();
516 
517     for (NumberOfParametersType par = 0; par < numberOfLocalParameters; par++)
518       {
519       deriv[par] = NumericTraits<DerivativeValueType>::ZeroValue();
520       for (ImageDimensionType dim = 0; dim < TImageToImageMetric::MovingImageDimension; dim++)
521         {
522         deriv[par] += derivWRTImage[dim] * jacobian(dim, par);
523         }
524       }
525     }
526 }
527 
528 /*
529  * Specific implementation for sparse threader. It reuse most of the routine from the dense threader by
530  * reinitializing the scanning at every point.
531  */
532 template < typename TDomainPartitioner, typename TImageToImageMetric, typename TNeighborhoodCorrelationMetric >
533 bool
534 ANTSNeighborhoodCorrelationImageToImageMetricv4GetValueAndDerivativeThreader< TDomainPartitioner, TImageToImageMetric, TNeighborhoodCorrelationMetric >
ProcessVirtualPoint_impl(IdentityHelper<ThreadedIndexedContainerPartitioner> itkNotUsed (self),const VirtualIndexType & virtualIndex,const VirtualPointType & itkNotUsed (virtualPoint),const ThreadIdType threadId)535 ::ProcessVirtualPoint_impl(IdentityHelper<ThreadedIndexedContainerPartitioner> itkNotUsed(self),
536     const VirtualIndexType & virtualIndex, const VirtualPointType & itkNotUsed(virtualPoint),
537     const ThreadIdType threadId )
538 {
539 
540   MeasureType          metricValueResult = NumericTraits< MeasureType >::ZeroValue();
541   bool                 pointIsValid;
542   ScanIteratorType     scanIt;
543   ScanParametersType   scanParameters;
544   ScanMemType          scanMem;
545 
546   DerivativeType & localDerivativeResult = this->m_GetValueAndDerivativePerThreadVariables[threadId].LocalDerivatives;
547 
548 
549   // convert virtualPoint to a single point region
550   ImageRegionType singlePointRegion;
551   singlePointRegion.SetIndex(virtualIndex);
552   typename ImageRegionType::SizeType singlePointSize;
553   singlePointSize.Fill(1);
554   singlePointRegion.SetSize(singlePointSize);
555 
556   // use scanning variables just for a single point region
557   // iterate over the single point and initialize the scanning variables
558   this->InitializeScanning( singlePointRegion, scanIt, scanMem, scanParameters );
559 
560   /* Iterate over the sub region of a single point*/
561   scanIt.GoToBegin();
562 
563   try
564     {
565     pointIsValid = false;
566 
567     this->UpdateQueues(scanIt, scanMem, scanParameters, threadId);
568     pointIsValid = this->ComputeInformationFromQueues(scanIt, scanMem, scanParameters, threadId);
569     if( pointIsValid )
570       {
571       this->ComputeMovingTransformDerivative(scanIt, scanMem, scanParameters, localDerivativeResult, metricValueResult, threadId );
572       }
573     }
574   catch (ExceptionObject & exc)
575     {
576     //NOTE: there must be a cleaner way to do this:
577     std::string msg("Caught exception: \n");
578     msg += exc.what();
579     ExceptionObject err(__FILE__, __LINE__, msg);
580     throw err;
581     }
582 
583   /* Assign the results */
584   if ( pointIsValid )
585     {
586     this->m_GetValueAndDerivativePerThreadVariables[threadId].NumberOfValidPoints++;
587     this->m_GetValueAndDerivativePerThreadVariables[threadId].Measure -= metricValueResult;
588     /* Store the result. This depends on what type of
589      * transform is being used. */
590     if( this->GetComputeDerivative() )
591       {
592       this->StorePointDerivativeResult( scanIt.GetIndex(), threadId );
593       }
594     }
595 
596   return pointIsValid;
597 }
598 
599 
600 } // end namespace itk
601 
602 #endif
603