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