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 itkMeanSquaresImageToImageMetricv4GetValueAndDerivativeThreader_hxx
19 #define itkMeanSquaresImageToImageMetricv4GetValueAndDerivativeThreader_hxx
20
21 #include "itkMeanSquaresImageToImageMetricv4GetValueAndDerivativeThreader.h"
22 #include "itkDefaultConvertPixelTraits.h"
23
24 namespace itk
25 {
26
27 template< typename TDomainPartitioner, typename TImageToImageMetric, typename TMeanSquaresMetric >
28 bool
29 MeanSquaresImageToImageMetricv4GetValueAndDerivativeThreader< TDomainPartitioner, TImageToImageMetric, TMeanSquaresMetric >
ProcessPoint(const VirtualIndexType &,const VirtualPointType & virtualPoint,const FixedImagePointType &,const FixedImagePixelType & fixedImageValue,const FixedImageGradientType &,const MovingImagePointType &,const MovingImagePixelType & movingImageValue,const MovingImageGradientType & movingImageGradient,MeasureType & metricValueReturn,DerivativeType & localDerivativeReturn,const ThreadIdType threadId) const30 ::ProcessPoint( const VirtualIndexType &,
31 const VirtualPointType & virtualPoint,
32 const FixedImagePointType &,
33 const FixedImagePixelType & fixedImageValue,
34 const FixedImageGradientType &,
35 const MovingImagePointType & ,
36 const MovingImagePixelType & movingImageValue,
37 const MovingImageGradientType & movingImageGradient,
38 MeasureType & metricValueReturn,
39 DerivativeType & localDerivativeReturn,
40 const ThreadIdType threadId) const
41 {
42 /** Only the voxelwise contribution given the point pairs. */
43 FixedImagePixelType diff = fixedImageValue - movingImageValue;
44 const unsigned int nComponents = NumericTraits<FixedImagePixelType>::GetLength( diff );
45 metricValueReturn = NumericTraits<MeasureType>::ZeroValue();
46
47 for ( unsigned int nc = 0; nc < nComponents; nc++ )
48 {
49 MeasureType diffC = DefaultConvertPixelTraits<FixedImagePixelType>::GetNthComponent(nc, diff);
50 metricValueReturn += diffC*diffC;
51 }
52
53 if( ! this->GetComputeDerivative() )
54 {
55 return true;
56 }
57
58 /* Use a pre-allocated jacobian object for efficiency */
59 using JacobianReferenceType = typename TImageToImageMetric::JacobianType &;
60 JacobianReferenceType jacobian = this->m_GetValueAndDerivativePerThreadVariables[threadId].MovingTransformJacobian;
61 JacobianReferenceType jacobianPositional = this->m_GetValueAndDerivativePerThreadVariables[threadId].MovingTransformJacobianPositional;
62
63 /** For dense transforms, this returns identity */
64 this->m_Associate->GetMovingTransform()->
65 ComputeJacobianWithRespectToParametersCachedTemporaries(virtualPoint,
66 jacobian,
67 jacobianPositional);
68
69 for ( unsigned int par = 0; par < this->GetCachedNumberOfLocalParameters(); par++ )
70 {
71 localDerivativeReturn[par] = NumericTraits<DerivativeValueType>::ZeroValue();
72 for ( unsigned int nc = 0; nc < nComponents; nc++ )
73 {
74 MeasureType diffValue = DefaultConvertPixelTraits<FixedImagePixelType>::GetNthComponent(nc,diff);
75 for ( SizeValueType dim = 0; dim < ImageToImageMetricv4Type::MovingImageDimension; dim++ )
76 {
77 localDerivativeReturn[par] += 2.0 * diffValue * jacobian(dim, par) *
78 DefaultConvertPixelTraits<MovingImageGradientType>::GetNthComponent(
79 ImageToImageMetricv4Type::FixedImageDimension * nc + dim, movingImageGradient );
80 }
81 }
82 }
83 return true;
84 }
85
86 } // end namespace itk
87
88 #endif
89