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 itkMeanSquareRegistrationFunction_hxx
19 #define itkMeanSquareRegistrationFunction_hxx
20
21 #include "itkMeanSquareRegistrationFunction.h"
22 #include "itkMacro.h"
23 #include "itkMath.h"
24 #include "itkMath.h"
25
26 namespace itk
27 {
28 /**
29 * Default constructor
30 */
31 template< typename TFixedImage, typename TMovingImage, typename TDisplacementField >
32 MeanSquareRegistrationFunction< TFixedImage, TMovingImage, TDisplacementField >
MeanSquareRegistrationFunction()33 ::MeanSquareRegistrationFunction()
34 {
35 RadiusType r;
36 unsigned int j;
37
38 for ( j = 0; j < ImageDimension; j++ )
39 {
40 r[j] = 0;
41 }
42 this->SetRadius(r);
43
44 this->SetEnergy(0.0);
45 m_TimeStep = 1.0;
46 m_DenominatorThreshold = 1e-9;
47 m_IntensityDifferenceThreshold = 0.001;
48 this->SetMovingImage(nullptr);
49 this->SetFixedImage(nullptr);
50 m_FixedImageGradientCalculator = GradientCalculatorType::New();
51
52 typename DefaultInterpolatorType::Pointer interp =
53 DefaultInterpolatorType::New();
54
55 m_MovingImageInterpolator = itkDynamicCastInDebugMode< InterpolatorType * >(interp.GetPointer() );
56 }
57
58 /*
59 * Standard "PrintSelf" method.
60 */
61 template< typename TFixedImage, typename TMovingImage, typename TDisplacementField >
62 void
63 MeanSquareRegistrationFunction< TFixedImage, TMovingImage, TDisplacementField >
PrintSelf(std::ostream & os,Indent indent) const64 ::PrintSelf(std::ostream & os, Indent indent) const
65 {
66 Superclass::PrintSelf(os, indent);
67 /*
68 os << indent << "MovingImageIterpolator: ";
69 os << m_MovingImageInterpolator.GetPointer() << std::endl;
70 os << indent << "FixedImageGradientCalculator: ";
71 os << m_FixedImageGradientCalculator.GetPointer() << std::endl;
72 os << indent << "DenominatorThreshold: ";
73 os << m_DenominatorThreshold << std::endl;
74 os << indent << "IntensityDifferenceThreshold: ";
75 os << m_IntensityDifferenceThreshold << std::endl;
76 */
77 }
78
79 /*
80 * Set the function state values before each iteration
81 */
82 template< typename TFixedImage, typename TMovingImage, typename TDisplacementField >
83 void
84 MeanSquareRegistrationFunction< TFixedImage, TMovingImage, TDisplacementField >
InitializeIteration()85 ::InitializeIteration()
86 {
87 if ( !this->GetMovingImage() || !this->GetFixedImage() || !m_MovingImageInterpolator )
88 {
89 itkExceptionMacro(<< "MovingImage, FixedImage and/or Interpolator not set");
90 }
91
92 // cache fixed image information
93 m_FixedImageSpacing = this->GetFixedImage()->GetSpacing();
94
95 // setup gradient calculator
96 m_FixedImageGradientCalculator->SetInputImage( this->GetFixedImage() );
97
98 // setup moving image interpolator
99 m_MovingImageInterpolator->SetInputImage( this->GetMovingImage() );
100
101 this->SetEnergy(0.0);
102 }
103
104 /**
105 * Compute update at a non boundary neighbourhood
106 */
107 template< typename TFixedImage, typename TMovingImage, typename TDisplacementField >
108 typename MeanSquareRegistrationFunction< TFixedImage, TMovingImage, TDisplacementField >
109 ::PixelType
110 MeanSquareRegistrationFunction< TFixedImage, TMovingImage, TDisplacementField >
ComputeUpdate(const NeighborhoodType & it,void * itkNotUsed (globalData),const FloatOffsetType & itkNotUsed (offset))111 ::ComputeUpdate( const NeighborhoodType & it, void *itkNotUsed(globalData),
112 const FloatOffsetType & itkNotUsed(offset) )
113 {
114 // Get fixed image related information
115 // Note: no need to check the index is within
116 // fixed image buffer. This is done by the external filter.
117 const IndexType index = it.GetIndex();
118 const auto fixedValue = (double)this->GetFixedImage()->GetPixel(index);
119 const CovariantVectorType fixedGradient = m_FixedImageGradientCalculator->EvaluateAtIndex(index);
120 double fixedGradientSquaredMagnitude = 0;
121
122 for ( unsigned int j = 0; j < ImageDimension; j++ )
123 {
124 fixedGradientSquaredMagnitude += itk::Math::sqr(fixedGradient[j]) * m_FixedImageSpacing[j];
125 }
126
127 // Get moving image related information
128 const DisplacementFieldPixelType itvec = this->GetDisplacementField()->GetPixel(index);
129 PointType mappedPoint;
130 this->GetFixedImage()->TransformIndexToPhysicalPoint(index, mappedPoint);
131 for ( unsigned int j = 0; j < ImageDimension; j++ )
132 {
133 mappedPoint[j] += itvec[j];
134 }
135 double movingValue = 0.0;
136 if ( m_MovingImageInterpolator->IsInsideBuffer(mappedPoint) )
137 {
138 movingValue = m_MovingImageInterpolator->Evaluate(mappedPoint);
139 }
140
141 // Compute update
142 const double speedValue = fixedValue - movingValue;
143 this->m_Energy += speedValue * speedValue;
144
145 const bool normalizemetric = this->GetNormalizeGradient();
146 double denominator = 1.0;
147 if ( normalizemetric )
148 {
149 denominator = speedValue * speedValue * fixedGradientSquaredMagnitude;
150 denominator = std::sqrt(denominator);
151 }
152 if ( Math::AlmostEquals( denominator, 0.0 ) )
153 {
154 denominator = 1.0;
155 }
156 PixelType update;
157 if ( itk::Math::abs(speedValue) < m_IntensityDifferenceThreshold
158 || denominator < m_DenominatorThreshold )
159 {
160 update.Fill(0.0);
161 return update;
162 }
163
164 for ( unsigned int j = 0; j < ImageDimension; j++ )
165 {
166 update[j] = speedValue * fixedGradient[j] * itk::Math::sqr(m_FixedImageSpacing[j])
167 / denominator * this->m_GradientStep;
168 }
169 return update;
170 }
171 } // end namespace itk
172
173 #endif
174