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 itkObjectToObjectMultiMetricv4_hxx
19 #define itkObjectToObjectMultiMetricv4_hxx
20
21 #include "itkObjectToObjectMultiMetricv4.h"
22 #include "itkCompositeTransform.h"
23
24 namespace itk
25 {
26
27 /** Constructor */
28 template<unsigned int TFixedDimension, unsigned int TMovingDimension, typename TVirtualImage, typename TInternalComputationValueType>
29 ObjectToObjectMultiMetricv4<TFixedDimension, TMovingDimension, TVirtualImage, TInternalComputationValueType>
ObjectToObjectMultiMetricv4()30 ::ObjectToObjectMultiMetricv4()
31 {
32 this->m_MetricQueue.clear();
33
34 //We want the moving transform to be nullptr by default
35 this->m_MovingTransform = nullptr;
36 }
37
38 /** Add a metric to the queue. */
39 template<unsigned int TFixedDimension, unsigned int TMovingDimension, typename TVirtualImage, typename TInternalComputationValueType>
40 void
41 ObjectToObjectMultiMetricv4<TFixedDimension, TMovingDimension, TVirtualImage, TInternalComputationValueType>
AddMetric(MetricType * metric)42 ::AddMetric (MetricType* metric)
43 {
44 this->m_MetricQueue.push_back(metric);
45 }
46
47 /** Clear the queue */
48 template<unsigned int TFixedDimension, unsigned int TMovingDimension, typename TVirtualImage, typename TInternalComputationValueType>
49 void
50 ObjectToObjectMultiMetricv4<TFixedDimension, TMovingDimension, TVirtualImage, TInternalComputationValueType>
ClearMetricQueue()51 ::ClearMetricQueue()
52 {
53 this->m_MetricQueue.clear();
54 }
55
56 /** Get the number of metrics */
57 template<unsigned int TFixedDimension, unsigned int TMovingDimension, typename TVirtualImage, typename TInternalComputationValueType>
58 itk::SizeValueType
59 ObjectToObjectMultiMetricv4<TFixedDimension, TMovingDimension, TVirtualImage, TInternalComputationValueType>
GetNumberOfMetrics() const60 ::GetNumberOfMetrics() const
61 {
62 return static_cast<itk::SizeValueType>( this->m_MetricQueue.size() );
63 }
64
65 template<unsigned int TFixedDimension, unsigned int TMovingDimension, typename TVirtualImage, typename TInternalComputationValueType>
66 void
67 ObjectToObjectMultiMetricv4<TFixedDimension, TMovingDimension, TVirtualImage, TInternalComputationValueType>
SetMovingTransform(MovingTransformType * transform)68 ::SetMovingTransform( MovingTransformType * transform )
69 {
70 if( this->GetNumberOfMetrics() == 0 )
71 {
72 itkExceptionMacro("No metrics are assigned. Cannot assign transform.");
73 }
74 Superclass::SetMovingTransform( transform );
75 for (SizeValueType j = 0; j < this->GetNumberOfMetrics(); j++)
76 {
77 this->m_MetricQueue[j]->SetMovingTransform( transform );
78 }
79 }
80
81 template<unsigned int TFixedDimension, unsigned int TMovingDimension, typename TVirtualImage, typename TInternalComputationValueType>
82 void
83 ObjectToObjectMultiMetricv4<TFixedDimension, TMovingDimension, TVirtualImage, TInternalComputationValueType>
SetFixedTransform(FixedTransformType * transform)84 ::SetFixedTransform( FixedTransformType * transform )
85 {
86 if( this->GetNumberOfMetrics() == 0 )
87 {
88 itkExceptionMacro("No metrics are assigned. Cannot assign transform.");
89 }
90 Superclass::SetFixedTransform( transform );
91 for (SizeValueType j = 0; j < this->GetNumberOfMetrics(); j++)
92 {
93 this->m_MetricQueue[j]->SetFixedTransform( transform );
94 }
95 }
96
97 template<unsigned int TFixedDimension, unsigned int TMovingDimension, typename TVirtualImage, typename TInternalComputationValueType>
98 void
99 ObjectToObjectMultiMetricv4<TFixedDimension, TMovingDimension, TVirtualImage, TInternalComputationValueType>
Initialize()100 ::Initialize()
101 {
102 if( this->GetNumberOfMetrics() == 0 )
103 {
104 itkExceptionMacro("No metrics are assigned. Cannot evaluate.");
105 }
106
107 /* Verify derivative weights and initialize if appropriate */
108 if( this->m_MetricWeights.Size() > 0 )
109 {
110 if( this->m_MetricWeights.Size() != this->GetNumberOfMetrics() )
111 {
112 itkExceptionMacro("The derivative weights are not of the proper size. "
113 "Number of metrics: " << this->GetNumberOfMetrics() << ", "
114 "Number of weights: " << this->m_MetricWeights.Size() );
115 }
116 /* normalize the weights */
117 WeightValueType sum = NumericTraits<WeightValueType>::ZeroValue();
118 for (SizeValueType j = 0; j < this->GetNumberOfMetrics(); j++)
119 {
120 sum += this->m_MetricWeights[j];
121 }
122 if( sum <= NumericTraits<WeightValueType>::epsilon() )
123 {
124 itkExceptionMacro("The derivative weights are too small: " << this->m_MetricWeights );
125 }
126 for (SizeValueType j = 0; j < this->GetNumberOfMetrics(); j++)
127 {
128 this->m_MetricWeights[j] = this->m_MetricWeights[j] / sum;
129 }
130 }
131 else
132 {
133 /* Initialize to defaults */
134 this->m_MetricWeights.SetSize( this->GetNumberOfMetrics() );
135 this->m_MetricWeights.Fill( NumericTraits<WeightValueType>::OneValue() / static_cast<WeightValueType>(this->GetNumberOfMetrics()) );
136 }
137
138 /* resize */
139 this->m_MetricValueArray.SetSize( this->GetNumberOfMetrics() );
140
141 /* Verify the same transform is in all metrics. */
142 const MovingTransformType * firstTransform = nullptr;
143 for (SizeValueType j = 0; j < this->GetNumberOfMetrics(); j++)
144 {
145 const MovingTransformType * transform = this->m_MetricQueue[j]->GetMovingTransform();
146 //Check if it's a composite. If so, there must be only one transform set to be
147 // optimized, and it must be the same as in other metrics.
148 using CompositeType = CompositeTransform<typename MovingTransformType::ScalarType, TFixedDimension>;
149 const auto * composite = dynamic_cast<const CompositeType*>(transform);
150 if( composite != nullptr )
151 {
152 SizeValueType count = 0;
153 for( size_t n = 0; n < composite->GetNumberOfTransforms(); n++ )
154 {
155 if( composite->GetNthTransformToOptimize( static_cast<SizeValueType>( n ) ) )
156 {
157 count++;
158 transform = composite->GetNthTransformConstPointer(static_cast<SizeValueType>( n ) );
159 }
160 }
161 if( count != 1 )
162 {
163 itkExceptionMacro("Expected exactly one transform set to be optimized within the composite transform. Error with metric " << j << ".");
164 }
165 }
166
167 if( j == 0 )
168 {
169 firstTransform = transform;
170 }
171 else
172 {
173 if( transform != firstTransform )
174 {
175 itkExceptionMacro("One or more component metrics have different active transforms. "
176 "Each metric must be using the same transform object. For CompositeTransform, "
177 "there must be only one transform set to optimize, and it must be the same "
178 "as other metric transforms." );
179 }
180 }
181 }
182
183 /* Assign local pointers to common transforms */
184 if( this->m_MovingTransform.IsNull() )
185 {
186 Superclass::SetMovingTransform( const_cast<MovingTransformType*>(this->m_MetricQueue[0]->GetMovingTransform()) );
187 }
188 if( this->m_FixedTransform.IsNull() )
189 {
190 Superclass::SetFixedTransform( const_cast<MovingTransformType*>(this->m_MetricQueue[0]->GetFixedTransform()) );
191 }
192
193 /* Initialize individual metrics. */
194 for (SizeValueType j = 0; j < this->GetNumberOfMetrics(); j++)
195 {
196 try
197 {
198 this->m_MetricQueue[j]->Initialize();
199 }
200 catch(ExceptionObject &exc)
201 {
202 std::string msg("Caught exception initializing metric: \n");
203 msg += exc.what();
204 ExceptionObject err(__FILE__, __LINE__, msg);
205 throw err;
206 }
207 }
208
209 /* Get the first valid virtual domain and assign
210 * it to this metric as a common virtual domain,
211 * for direct use by calling classes. */
212 for (SizeValueType j = 0; j < this->GetNumberOfMetrics(); j++)
213 {
214 if( this->m_MetricQueue[j]->GetVirtualImage() )
215 {
216 this->SetVirtualDomainFromImage( this->m_MetricQueue[j]->GetVirtualImage() );
217 break;
218 }
219 }
220
221 /* Do this after we've setup local copy of virtual domain */
222 Superclass::Initialize();
223 }
224
225 template<unsigned int TFixedDimension, unsigned int TMovingDimension, typename TVirtualImage, typename TInternalComputationValueType>
226 typename ObjectToObjectMultiMetricv4<TFixedDimension, TMovingDimension, TVirtualImage, TInternalComputationValueType>::MeasureType
227 ObjectToObjectMultiMetricv4<TFixedDimension, TMovingDimension, TVirtualImage, TInternalComputationValueType>
GetValue() const228 ::GetValue() const
229 {
230 for (SizeValueType j = 0; j < this->GetNumberOfMetrics(); j++)
231 {
232 this->m_MetricValueArray[j] = this->m_MetricQueue[j]->GetValue();
233 }
234
235 MeasureType firstValue = this->m_MetricValueArray[0];
236 this->m_Value = firstValue;
237 return firstValue;
238 }
239
240 template<unsigned int TFixedDimension, unsigned int TMovingDimension, typename TVirtualImage, typename TInternalComputationValueType>
241 void
242 ObjectToObjectMultiMetricv4<TFixedDimension, TMovingDimension, TVirtualImage, TInternalComputationValueType>
GetDerivative(DerivativeType & derivativeResult) const243 ::GetDerivative(DerivativeType & derivativeResult) const
244 {
245 MeasureType firstValue;
246 this->GetValueAndDerivative( firstValue, derivativeResult );
247 }
248
249 template<unsigned int TFixedDimension, unsigned int TMovingDimension, typename TVirtualImage, typename TInternalComputationValueType>
250 void
251 ObjectToObjectMultiMetricv4<TFixedDimension, TMovingDimension, TVirtualImage, TInternalComputationValueType>
GetValueAndDerivative(MeasureType & firstValue,DerivativeType & derivativeResult) const252 ::GetValueAndDerivative(MeasureType & firstValue, DerivativeType & derivativeResult) const
253 {
254 if ( derivativeResult.GetSize() != this->GetNumberOfParameters() )
255 {
256 derivativeResult.SetSize( this->GetNumberOfParameters() );
257 }
258 derivativeResult.Fill( NumericTraits<DerivativeValueType>::ZeroValue() );
259
260 DerivativeType metricDerivative;
261 MeasureType metricValue = NumericTraits<MeasureType>::ZeroValue();
262
263 // Loop over metrics
264 DerivativeValueType totalMagnitude = NumericTraits<DerivativeValueType>::ZeroValue();
265 for (SizeValueType j = 0; j < this->GetNumberOfMetrics(); j++)
266 {
267 this->m_MetricQueue[j]->GetValueAndDerivative( metricValue, metricDerivative);
268 this->m_MetricValueArray[j] = metricValue;
269
270 DerivativeValueType magnitude = metricDerivative.magnitude();
271 DerivativeValueType weightOverMagnitude = NumericTraits<DerivativeValueType>::ZeroValue();
272 totalMagnitude += magnitude;
273
274 if( magnitude > NumericTraits<DerivativeValueType>::epsilon() )
275 {
276 weightOverMagnitude = this->m_MetricWeights[j] / magnitude;
277 }
278 // derivative = \sum_j w_j * (dM_j / ||dM_j||)
279 for( NumberOfParametersType p = 0; p < this->GetNumberOfParameters(); p++ )
280 {
281 // roll our own loop to avoid temporary variable that could be large when using displacement fields.
282 derivativeResult[p] += ( metricDerivative[p] * weightOverMagnitude );
283 }
284 }
285
286 // Scale by totalMagnitude to prevent what amounts to implicit step estimation from magnitude scaling.
287 // This keeps the behavior of this metric the same as a regular metric, with respect to derivative
288 // magnitudes.
289 totalMagnitude /= this->GetNumberOfMetrics();
290 for( NumberOfParametersType p = 0; p < this->GetNumberOfParameters(); p++ )
291 {
292 derivativeResult[p] *= totalMagnitude;
293 }
294
295 firstValue = this->m_MetricValueArray[0];
296 this->m_Value = firstValue;
297 }
298
299 template<unsigned int TFixedDimension, unsigned int TMovingDimension, typename TVirtualImage, typename TInternalComputationValueType>
300 typename ObjectToObjectMultiMetricv4<TFixedDimension, TMovingDimension, TVirtualImage, TInternalComputationValueType>::MetricValueArrayType
301 ObjectToObjectMultiMetricv4<TFixedDimension, TMovingDimension, TVirtualImage, TInternalComputationValueType>
GetValueArray() const302 ::GetValueArray() const
303 {
304 return this->m_MetricValueArray;
305 }
306
307 template<unsigned int TFixedDimension, unsigned int TMovingDimension, typename TVirtualImage, typename TInternalComputationValueType>
308 typename ObjectToObjectMultiMetricv4<TFixedDimension, TMovingDimension, TVirtualImage, TInternalComputationValueType>::MeasureType
309 ObjectToObjectMultiMetricv4<TFixedDimension, TMovingDimension, TVirtualImage, TInternalComputationValueType>
GetWeightedValue() const310 ::GetWeightedValue() const
311 {
312 MeasureType value = NumericTraits<MeasureType>::ZeroValue();
313
314 for (SizeValueType j = 0; j < this->GetNumberOfMetrics(); j++)
315 {
316 // value = sum_j w_j * M_j
317 value += this->m_MetricValueArray[j] * this->m_MetricWeights[j];
318 }
319 return value;
320 }
321
322 template<unsigned int TFixedDimension, unsigned int TMovingDimension, typename TVirtualImage, typename TInternalComputationValueType>
323 const typename ObjectToObjectMultiMetricv4<TFixedDimension, TMovingDimension, TVirtualImage, TInternalComputationValueType>::MetricQueueType &
324 ObjectToObjectMultiMetricv4<TFixedDimension, TMovingDimension, TVirtualImage, TInternalComputationValueType>
GetMetricQueue() const325 ::GetMetricQueue() const
326 {
327 return this->m_MetricQueue;
328 }
329
330
331 template<unsigned int TFixedDimension, unsigned int TMovingDimension, typename TVirtualImage, typename TInternalComputationValueType>
332 bool
333 ObjectToObjectMultiMetricv4<TFixedDimension, TMovingDimension, TVirtualImage, TInternalComputationValueType>
SupportsArbitraryVirtualDomainSamples() const334 ::SupportsArbitraryVirtualDomainSamples() const
335 {
336 for (SizeValueType j = 0; j < this->GetNumberOfMetrics(); j++)
337 {
338 if( ! this->m_MetricQueue[j]->SupportsArbitraryVirtualDomainSamples() )
339 {
340 return false;
341 }
342 }
343 return true;
344 }
345
346 template<unsigned int TFixedDimension, unsigned int TMovingDimension, typename TVirtualImage, typename TInternalComputationValueType>
347 void
348 ObjectToObjectMultiMetricv4<TFixedDimension, TMovingDimension, TVirtualImage, TInternalComputationValueType>
PrintSelf(std::ostream & os,Indent indent) const349 ::PrintSelf(std::ostream & os, Indent indent) const
350 {
351 os << indent << "Weights of metric derivatives: " << this->m_MetricWeights << std::endl;
352 os << indent << "The multivariate contains the following metrics: " << std::endl << std::endl;
353 for (SizeValueType i = 0; i < this->GetNumberOfMetrics(); i++)
354 {
355 os << indent << "~~~ Metric " << i << " ~~~" << std::endl;
356 this->m_MetricQueue[i]->Print(os, indent.GetNextIndent() );
357 }
358 }
359
360 } // end namespace itk
361
362 #endif //itkObjectToObjectMultiMetricv4_hxx
363