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 itkCompositeTransform_hxx
19 #define itkCompositeTransform_hxx
20 
21 #include "itkCompositeTransform.h"
22 
23 namespace itk
24 {
25 
26 
27 template
28 <typename TParametersValueType, unsigned int NDimensions>
CompositeTransform()29 CompositeTransform<TParametersValueType, NDimensions>::CompositeTransform()
30 {
31   this->m_TransformsToOptimizeFlags.clear();
32   this->m_TransformsToOptimizeQueue.clear();
33   this->m_PreviousTransformsToOptimizeUpdateTime = 0;
34 }
35 
36 template
37 <typename TParametersValueType, unsigned int NDimensions>
38 typename CompositeTransform<TParametersValueType, NDimensions>::TransformCategoryType
39 CompositeTransform<TParametersValueType, NDimensions>
GetTransformCategory() const40 ::GetTransformCategory() const
41 {
42   // Check if linear
43   bool isLinearTransform = this->IsLinear();
44   if( isLinearTransform )
45     {
46     return Self::Linear;
47     }
48 
49   // Check if displacement field
50   bool isDisplacementFieldTransform = true;
51   for( signed long tind = static_cast<signed long>( this->GetNumberOfTransforms() ) - 1; tind >= 0; tind-- )
52     {
53     if( this->GetNthTransformToOptimize( tind ) &&
54       ( this->GetNthTransformConstPointer( tind )->GetTransformCategory() != Self::DisplacementField ) )
55       {
56       isDisplacementFieldTransform = false;
57       break;
58       }
59     }
60 
61   if( isDisplacementFieldTransform )
62     {
63     return Self::DisplacementField;
64     }
65   else
66     {
67     return Self::UnknownTransformCategory;
68     }
69 }
70 
71 
72 template
73 <typename TParametersValueType, unsigned int NDimensions>
74 typename CompositeTransform<TParametersValueType, NDimensions>
75 ::OutputPointType
76 CompositeTransform<TParametersValueType, NDimensions>
TransformPoint(const InputPointType & inputPoint) const77 ::TransformPoint( const InputPointType& inputPoint ) const
78 {
79 
80   /* Apply in reverse queue order.  */
81   typename TransformQueueType::const_iterator it( this->m_TransformQueue.end() );
82   const typename TransformQueueType::const_iterator beginit( this->m_TransformQueue.begin() );
83   OutputPointType outputPoint( inputPoint );
84   do
85     {
86     it--;
87     outputPoint = (*it)->TransformPoint( outputPoint );
88     }
89   while( it != beginit );
90   return outputPoint;
91 }
92 
93 
94 template<typename TParametersValueType, unsigned int NDimensions>
95 typename CompositeTransform<TParametersValueType, NDimensions>
96 ::OutputVectorType
97 CompositeTransform<TParametersValueType, NDimensions>
TransformVector(const InputVectorType & inputVector) const98 ::TransformVector( const InputVectorType & inputVector ) const
99 {
100   OutputVectorType outputVector( inputVector );
101 
102   typename TransformQueueType::const_iterator it;
103   /* Apply in reverse queue order.  */
104   it = this->m_TransformQueue.end();
105 
106   do
107     {
108     it--;
109     outputVector = (*it)->TransformVector( outputVector );
110     }
111   while( it != this->m_TransformQueue.begin() );
112 
113   return outputVector;
114 }
115 
116 
117 template<typename TParametersValueType, unsigned int NDimensions>
118 typename CompositeTransform<TParametersValueType, NDimensions>
119 ::OutputVectorType
120 CompositeTransform<TParametersValueType, NDimensions>
TransformVector(const InputVectorType & inputVector,const InputPointType & inputPoint) const121 ::TransformVector( const InputVectorType & inputVector, const InputPointType & inputPoint ) const
122 {
123   OutputVectorType outputVector( inputVector );
124   OutputPointType outputPoint( inputPoint );
125 
126   typename TransformQueueType::const_iterator it;
127   /* Apply in reverse queue order.  */
128   it = this->m_TransformQueue.end();
129 
130   do
131     {
132     it--;
133     outputVector = (*it)->TransformVector( outputVector, outputPoint );
134     outputPoint = (*it)->TransformPoint( outputPoint );
135     }
136   while( it != this->m_TransformQueue.begin() );
137 
138   return outputVector;
139 }
140 
141 
142 template<typename TParametersValueType, unsigned int NDimensions>
143 typename CompositeTransform<TParametersValueType, NDimensions>
144 ::OutputVnlVectorType
145 CompositeTransform<TParametersValueType, NDimensions>
TransformVector(const InputVnlVectorType & inputVector,const InputPointType & inputPoint) const146 ::TransformVector( const InputVnlVectorType & inputVector, const InputPointType & inputPoint ) const
147 {
148   OutputVnlVectorType outputVector( inputVector );
149   OutputPointType outputPoint( inputPoint );
150 
151   typename TransformQueueType::const_iterator it;
152   /* Apply in reverse queue order.  */
153   it = this->m_TransformQueue.end();
154 
155   do
156     {
157     it--;
158     outputVector = (*it)->TransformVector( outputVector, outputPoint );
159     outputPoint = (*it)->TransformPoint( outputPoint );
160     }
161   while( it != this->m_TransformQueue.begin() );
162 
163   return outputVector;
164 }
165 
166 
167 template<typename TParametersValueType, unsigned int NDimensions>
168 typename CompositeTransform<TParametersValueType, NDimensions>
169 ::OutputVnlVectorType
170 CompositeTransform<TParametersValueType, NDimensions>
TransformVector(const InputVnlVectorType & inputVector) const171 ::TransformVector( const InputVnlVectorType & inputVector) const
172 {
173   OutputVnlVectorType outputVector( inputVector );
174 
175   typename TransformQueueType::const_iterator it;
176   /* Apply in reverse queue order.  */
177   it = this->m_TransformQueue.end();
178 
179   do
180     {
181     it--;
182     outputVector = (*it)->TransformVector( outputVector );
183     }
184   while( it != this->m_TransformQueue.begin() );
185 
186   return outputVector;
187 }
188 
189 
190 template<typename TParametersValueType, unsigned int NDimensions>
191 typename CompositeTransform<TParametersValueType, NDimensions>
192 ::OutputVectorPixelType
193 CompositeTransform<TParametersValueType, NDimensions>
TransformVector(const InputVectorPixelType & inputVector) const194 ::TransformVector( const InputVectorPixelType & inputVector ) const
195 {
196   OutputVectorPixelType outputVector( inputVector );
197 
198   typename TransformQueueType::const_iterator it;
199   /* Apply in reverse queue order.  */
200   it = this->m_TransformQueue.end();
201 
202   do
203     {
204     it--;
205     outputVector = (*it)->TransformVector( outputVector );
206     }
207   while( it != this->m_TransformQueue.begin() );
208 
209   return outputVector;
210 }
211 
212 
213 template<typename TParametersValueType, unsigned int NDimensions>
214 typename CompositeTransform<TParametersValueType, NDimensions>
215 ::OutputVectorPixelType
216 CompositeTransform<TParametersValueType, NDimensions>
TransformVector(const InputVectorPixelType & inputVector,const InputPointType & inputPoint) const217 ::TransformVector( const InputVectorPixelType & inputVector, const InputPointType & inputPoint ) const
218 {
219   OutputVectorPixelType outputVector( inputVector );
220   OutputPointType outputPoint( inputPoint );
221 
222   typename TransformQueueType::const_iterator it;
223   /* Apply in reverse queue order.  */
224   it = this->m_TransformQueue.end();
225 
226   do
227     {
228     it--;
229     outputVector = (*it)->TransformVector( outputVector, outputPoint );
230     outputPoint = (*it)->TransformPoint( outputPoint );
231     }
232   while( it != this->m_TransformQueue.begin() );
233 
234   return outputVector;
235 }
236 
237 
238 template<typename TParametersValueType, unsigned int NDimensions>
239 typename CompositeTransform<TParametersValueType, NDimensions>
240 ::OutputCovariantVectorType
241 CompositeTransform<TParametersValueType, NDimensions>
TransformCovariantVector(const InputCovariantVectorType & inputVector) const242 ::TransformCovariantVector( const InputCovariantVectorType & inputVector ) const
243 {
244   OutputCovariantVectorType outputVector( inputVector );
245 
246   typename TransformQueueType::const_iterator it;
247   /* Apply in reverse queue order.  */
248   it = this->m_TransformQueue.end();
249 
250   do
251     {
252     it--;
253     outputVector = (*it)->TransformCovariantVector( outputVector );
254     }
255   while( it != this->m_TransformQueue.begin() );
256 
257   return outputVector;
258 }
259 
260 
261 template<typename TParametersValueType, unsigned int NDimensions>
262 typename CompositeTransform<TParametersValueType, NDimensions>
263 ::OutputCovariantVectorType
264 CompositeTransform<TParametersValueType, NDimensions>
TransformCovariantVector(const InputCovariantVectorType & inputVector,const InputPointType & inputPoint) const265 ::TransformCovariantVector( const InputCovariantVectorType & inputVector, const InputPointType & inputPoint ) const
266 {
267   OutputCovariantVectorType outputVector( inputVector );
268   OutputPointType outputPoint( inputPoint );
269 
270   typename TransformQueueType::const_iterator it;
271   /* Apply in reverse queue order.  */
272   it = this->m_TransformQueue.end();
273 
274   do
275     {
276     it--;
277     outputVector = (*it)->TransformCovariantVector( outputVector, outputPoint );
278     outputPoint = (*it)->TransformPoint( outputPoint );
279     }
280   while( it != this->m_TransformQueue.begin() );
281 
282   return outputVector;
283 }
284 
285 
286 template<typename TParametersValueType, unsigned int NDimensions>
287 typename CompositeTransform<TParametersValueType, NDimensions>
288 ::OutputVectorPixelType
289 CompositeTransform<TParametersValueType, NDimensions>
TransformCovariantVector(const InputVectorPixelType & inputVector) const290 ::TransformCovariantVector( const InputVectorPixelType & inputVector ) const
291 {
292   OutputVectorPixelType outputVector( inputVector );
293 
294   typename TransformQueueType::const_iterator it;
295   /* Apply in reverse queue order.  */
296   it = this->m_TransformQueue.end();
297 
298   do
299     {
300     it--;
301     outputVector = (*it)->TransformCovariantVector( outputVector );
302     }
303   while( it != this->m_TransformQueue.begin() );
304 
305   return outputVector;
306 }
307 
308 
309 template<typename TParametersValueType, unsigned int NDimensions>
310 typename CompositeTransform<TParametersValueType, NDimensions>
311 ::OutputVectorPixelType
312 CompositeTransform<TParametersValueType, NDimensions>
TransformCovariantVector(const InputVectorPixelType & inputVector,const InputPointType & inputPoint) const313 ::TransformCovariantVector( const InputVectorPixelType & inputVector, const InputPointType & inputPoint ) const
314 {
315   OutputVectorPixelType outputVector( inputVector );
316   OutputPointType outputPoint( inputPoint );
317 
318   typename TransformQueueType::const_iterator it;
319   /* Apply in reverse queue order.  */
320   it = this->m_TransformQueue.end();
321 
322   do
323     {
324     it--;
325     outputVector = (*it)->TransformCovariantVector( outputVector, outputPoint );
326     outputPoint = (*it)->TransformPoint( outputPoint );
327     }
328   while( it != this->m_TransformQueue.begin() );
329 
330   return outputVector;
331 }
332 
333 
334 template<typename TParametersValueType, unsigned int NDimensions>
335 typename CompositeTransform<TParametersValueType, NDimensions>
336 ::OutputDiffusionTensor3DType
337 CompositeTransform<TParametersValueType, NDimensions>
TransformDiffusionTensor3D(const InputDiffusionTensor3DType & inputTensor,const InputPointType & inputPoint) const338 ::TransformDiffusionTensor3D( const InputDiffusionTensor3DType & inputTensor, const InputPointType & inputPoint ) const
339 {
340   OutputDiffusionTensor3DType outputTensor( inputTensor );
341   OutputPointType outputPoint( inputPoint );
342 
343   typename TransformQueueType::const_iterator it;
344   /* Apply in reverse queue order.  */
345   it = this->m_TransformQueue.end();
346 
347   do
348     {
349     it--;
350     outputTensor = (*it)->TransformDiffusionTensor3D( outputTensor, outputPoint );
351     outputPoint = (*it)->TransformPoint( outputPoint );
352     }
353   while( it != this->m_TransformQueue.begin() );
354 
355   return outputTensor;
356 }
357 
358 
359 template<typename TParametersValueType, unsigned int NDimensions>
360 typename CompositeTransform<TParametersValueType, NDimensions>
361 ::OutputVectorPixelType
362 CompositeTransform<TParametersValueType, NDimensions>
TransformDiffusionTensor3D(const InputVectorPixelType & inputTensor,const InputPointType & inputPoint) const363 ::TransformDiffusionTensor3D( const InputVectorPixelType & inputTensor, const InputPointType & inputPoint ) const
364 {
365   OutputVectorPixelType outputTensor( inputTensor );
366   OutputPointType outputPoint( inputPoint );
367 
368   typename TransformQueueType::const_iterator it;
369   /* Apply in reverse queue order.  */
370   it = this->m_TransformQueue.end();
371 
372   do
373     {
374     it--;
375     outputTensor = (*it)->TransformDiffusionTensor3D( outputTensor, outputPoint );
376     outputPoint = (*it)->TransformPoint( outputPoint );
377     }
378   while( it != this->m_TransformQueue.begin() );
379 
380   return outputTensor;
381 }
382 
383 
384 template<typename TParametersValueType, unsigned int NDimensions>
385 typename CompositeTransform<TParametersValueType, NDimensions>
386 ::OutputDiffusionTensor3DType
387 CompositeTransform<TParametersValueType, NDimensions>
TransformDiffusionTensor3D(const InputDiffusionTensor3DType & inputTensor) const388 ::TransformDiffusionTensor3D( const InputDiffusionTensor3DType & inputTensor ) const
389 {
390   OutputDiffusionTensor3DType outputTensor( inputTensor );
391 
392   typename TransformQueueType::const_iterator it;
393   /* Apply in reverse queue order.  */
394   it = this->m_TransformQueue.end();
395 
396   do
397     {
398     it--;
399     outputTensor = (*it)->TransformDiffusionTensor3D( outputTensor );
400     }
401   while( it != this->m_TransformQueue.begin() );
402 
403   return outputTensor;
404 }
405 
406 
407 template<typename TParametersValueType, unsigned int NDimensions>
408 typename CompositeTransform<TParametersValueType, NDimensions>
409 ::OutputVectorPixelType
410 CompositeTransform<TParametersValueType, NDimensions>
TransformDiffusionTensor3D(const InputVectorPixelType & inputTensor) const411 ::TransformDiffusionTensor3D( const InputVectorPixelType & inputTensor ) const
412 {
413   OutputVectorPixelType outputTensor( inputTensor );
414 
415   typename TransformQueueType::const_iterator it;
416   /* Apply in reverse queue order.  */
417   it = this->m_TransformQueue.end();
418 
419   do
420     {
421     it--;
422     outputTensor = (*it)->TransformDiffusionTensor3D( outputTensor );
423     }
424   while( it != this->m_TransformQueue.begin() );
425 
426   return outputTensor;
427 }
428 
429 
430 template<typename TParametersValueType, unsigned int NDimensions>
431 typename CompositeTransform<TParametersValueType, NDimensions>
432 ::OutputSymmetricSecondRankTensorType
433 CompositeTransform<TParametersValueType, NDimensions>
TransformSymmetricSecondRankTensor(const InputSymmetricSecondRankTensorType & inputTensor,const InputPointType & inputPoint) const434 ::TransformSymmetricSecondRankTensor( const InputSymmetricSecondRankTensorType & inputTensor, const InputPointType & inputPoint ) const
435 {
436   OutputSymmetricSecondRankTensorType outputTensor( inputTensor );
437   OutputPointType outputPoint( inputPoint );
438 
439   typename TransformQueueType::const_iterator it;
440   /* Apply in reverse queue order.  */
441   it = this->m_TransformQueue.end();
442 
443   do
444     {
445     it--;
446     outputTensor = (*it)->TransformSymmetricSecondRankTensor( outputTensor, outputPoint );
447     outputPoint = (*it)->TransformPoint( outputPoint );
448     }
449   while( it != this->m_TransformQueue.begin() );
450 
451   return outputTensor;
452 }
453 
454 
455 template<typename TParametersValueType, unsigned int NDimensions>
456 typename CompositeTransform<TParametersValueType, NDimensions>
457 ::OutputVectorPixelType
458 CompositeTransform<TParametersValueType, NDimensions>
TransformSymmetricSecondRankTensor(const InputVectorPixelType & inputTensor,const InputPointType & inputPoint) const459 ::TransformSymmetricSecondRankTensor( const InputVectorPixelType & inputTensor, const InputPointType & inputPoint ) const
460 {
461   OutputVectorPixelType outputTensor( inputTensor );
462   OutputPointType outputPoint( inputPoint );
463 
464   typename TransformQueueType::const_iterator it;
465   /* Apply in reverse queue order.  */
466   it = this->m_TransformQueue.end();
467 
468   do
469     {
470     it--;
471     outputTensor = (*it)->TransformSymmetricSecondRankTensor( outputTensor, outputPoint );
472     outputPoint = (*it)->TransformPoint( outputPoint );
473     }
474   while( it != this->m_TransformQueue.begin() );
475 
476   return outputTensor;
477 }
478 
479 
480 template<typename TParametersValueType, unsigned int NDimensions>
481 typename CompositeTransform<TParametersValueType, NDimensions>
482 ::OutputSymmetricSecondRankTensorType
483 CompositeTransform<TParametersValueType, NDimensions>
TransformSymmetricSecondRankTensor(const InputSymmetricSecondRankTensorType & inputTensor) const484 ::TransformSymmetricSecondRankTensor( const InputSymmetricSecondRankTensorType & inputTensor ) const
485 {
486   OutputSymmetricSecondRankTensorType outputTensor( inputTensor );
487 
488   typename TransformQueueType::const_iterator it;
489   /* Apply in reverse queue order.  */
490   it = this->m_TransformQueue.end();
491 
492   do
493     {
494     it--;
495     outputTensor = (*it)->TransformSymmetricSecondRankTensor( outputTensor );
496     }
497   while( it != this->m_TransformQueue.begin() );
498 
499   return outputTensor;
500 }
501 
502 
503 template<typename TParametersValueType, unsigned int NDimensions>
504 typename CompositeTransform<TParametersValueType, NDimensions>
505 ::OutputVectorPixelType
506 CompositeTransform<TParametersValueType, NDimensions>
TransformSymmetricSecondRankTensor(const InputVectorPixelType & inputTensor) const507 ::TransformSymmetricSecondRankTensor( const InputVectorPixelType & inputTensor ) const
508 {
509   OutputVectorPixelType outputTensor( inputTensor );
510 
511   typename TransformQueueType::const_iterator it;
512   /* Apply in reverse queue order.  */
513   it = this->m_TransformQueue.end();
514 
515   do
516     {
517     it--;
518     outputTensor = (*it)->TransformSymmetricSecondRankTensor( outputTensor );
519     }
520   while( it != this->m_TransformQueue.begin() );
521 
522   return outputTensor;
523 }
524 
525 
526 template<typename TParametersValueType, unsigned int NDimensions>
527 bool
528 CompositeTransform<TParametersValueType, NDimensions>
GetInverse(Self * inverse) const529 ::GetInverse( Self *inverse ) const
530 {
531   typename TransformQueueType::const_iterator it;
532 
533   //NOTE: CompositeTransform delegagtes to
534   //      individual transform for setting FixedParameters
535   //      inverse->SetFixedParameters( this->GetFixedParameters() );
536   inverse->ClearTransformQueue();
537   for( it = this->m_TransformQueue.begin(); it != this->m_TransformQueue.end(); ++it )
538     {
539     TransformTypePointer inverseTransform = dynamic_cast<TransformType *>( ( ( *it )->GetInverseTransform() ).GetPointer() );
540     if( !inverseTransform )
541       {
542       inverse->ClearTransformQueue();
543       return false;
544       }
545     else
546       {
547       /* Push to front to reverse the transform order */
548       inverse->PushFrontTransform( inverseTransform );
549       }
550     }
551 
552   /* Copy the optimization flags */
553   inverse->m_TransformsToOptimizeFlags.clear();
554   for( auto ofit = this->m_TransformsToOptimizeFlags.begin();
555        ofit != this->m_TransformsToOptimizeFlags.end(); ofit++ )
556     {
557     inverse->m_TransformsToOptimizeFlags.push_front( *ofit );
558     }
559 
560   return true;
561 }
562 
563 
564 template<typename TParametersValueType, unsigned int NDimensions>
565 typename CompositeTransform<TParametersValueType, NDimensions>
566 ::InverseTransformBasePointer
567 CompositeTransform<TParametersValueType, NDimensions>
GetInverseTransform() const568 ::GetInverseTransform() const
569 {
570   /* This method can't be defined in Superclass because of the call to New() */
571   Pointer inverseTransform = New();
572 
573   if( this->GetInverse( inverseTransform ) )
574     {
575     return inverseTransform.GetPointer();
576     }
577   else
578     {
579     return nullptr;
580     }
581 }
582 
583 
584 template<typename TParametersValueType, unsigned int NDimensions>
585 void
586 CompositeTransform<TParametersValueType, NDimensions>
ComputeJacobianWithRespectToParameters(const InputPointType & p,JacobianType & outJacobian) const587 ::ComputeJacobianWithRespectToParameters( const InputPointType & p, JacobianType & outJacobian ) const
588 {
589   /* Returns a concatenated MxN array, holding the Jacobian of each sub
590    * transform that is selected for optimization. The order is the same
591    * as that in which they're applied, i.e. reverse order.
592    * M rows = dimensionality of the transforms
593    * N cols = total number of parameters in the selected sub transforms. */
594   outJacobian.SetSize( NDimensions, this->GetNumberOfLocalParameters() );
595   JacobianType jacobianWithRespectToPosition;
596   this->ComputeJacobianWithRespectToParametersCachedTemporaries( p, outJacobian, jacobianWithRespectToPosition );
597 }
598 
599 template<typename TParametersValueType, unsigned int NDimensions>
600 void
601 CompositeTransform<TParametersValueType, NDimensions>
ComputeJacobianWithRespectToParametersCachedTemporaries(const InputPointType & p,JacobianType & outJacobian,JacobianType & cacheJacobian) const602 ::ComputeJacobianWithRespectToParametersCachedTemporaries( const InputPointType & p, JacobianType & outJacobian, JacobianType &cacheJacobian) const
603 {
604   //NOTE: This must have been done outside of outJacobian.SetSize( NDimensions, this->GetNumberOfLocalParameters() );
605   assert( outJacobian.rows() == NDimensions &&
606           outJacobian.cols() == this->GetNumberOfLocalParameters());
607 
608   if ( this->GetNumberOfTransforms() == 1)
609     {
610     const TransformType * const transform = this->GetNthTransformConstPointer( 0 );
611     transform->ComputeJacobianWithRespectToParameters( p, outJacobian );
612     return;
613     }
614 
615   NumberOfParametersType offset = NumericTraits< NumberOfParametersType >::ZeroValue();
616 
617   OutputPointType transformedPoint( p );
618 
619   /*
620    * Composite transform $T is composed of $T0(p0,x), $T1(p1,x) and $T2(p2, x) as:
621    *
622    * T(p0, p1, p2, x)
623    * = T0(p0, T1(p1, T2(p2, x)))
624    *
625    * p0, p1, p2 are the transform parameters for transform T0, T1, T2
626    * respectively.
627    *
628    * Let p = (p0, p1, p2).
629    *  x2 = T2(p2, x).
630    *  x1 = T1(p1, x2).
631    *
632    *
633    * The following loop computes dT/dp:
634    *
635    * dT/dp
636    * = (dT/dp0, dT/dp1, dT/dp2)
637    * = ( dT0/dp0 | x1 ),
638    *   ( dT0/dT1 | x1 ) * ( dT1/dp1 | x2 ),
639    *   ( ( dT0/dT1 | x1 ) * ( dT1/dT2 | x2 ) * ( dT2/dp2 | x )
640    *
641    * In the first iteration, it computes
642    *   dT2/dp2 | x
643    *
644    * In the second iteration, it computes
645    *   dT1/dp1 | x2
646    *
647    *  and it computes
648    *   dT1/dT2 | x2, and left multiplying to  dT2/dp2 | x
649    *
650    * In the third iteration, it computes
651    *   dT0/dp0 | x1,
652    *
653    *  and it computes
654    *   dT0/dT1 | x1, and left multiplying to
655    *    ( dT1/dT2 | x2 ) * ( dT2/dp2 | x )
656    *    and ( dT1/dp1 | x2 )
657    *
658    */
659   for( signed long tind = (signed long) this->GetNumberOfTransforms() - 1;
660        tind >= 0; --tind )
661     {
662     /* Get a raw pointer for efficiency, avoiding SmartPointer register/unregister */
663     const TransformType * const transform = this->GetNthTransformConstPointer( tind );
664 
665     const NumberOfParametersType offsetLast = offset;
666 
667     if( this->GetNthTransformToOptimize( tind ) )
668       {
669       /* Copy from another matrix, element-by-element */
670       /* The matrices are row-major, so block copy is less obviously
671        * better */
672 
673       const NumberOfParametersType numberOfLocalParameters = transform->GetNumberOfLocalParameters();
674 
675       cacheJacobian.SetSize( NDimensions, numberOfLocalParameters );
676       transform->ComputeJacobianWithRespectToParameters( transformedPoint, cacheJacobian );
677       outJacobian.update( cacheJacobian, 0, offset );
678       offset += numberOfLocalParameters;
679       }
680 
681     /* The composite transform needs to compose previous jacobians
682      * (those closer to the originating point) with the current
683      * transform's jacobian.  We therefore update the previous
684      * jacobian by multiplying the current matrix jumping over the
685      * first transform. The matrix here refers to  dT/dx at the point.
686      * For example, in the affine transform, this is the affine matrix.
687      *
688      * Also, noted the multiplication contains all the affine matrix from
689      * all transforms no matter they are going to be optimized or not.
690      *
691      * Update every old term by left multiplying dTk / dT{k-1}
692      * do this before computing the transformedPoint for the next
693      * iteration.
694      */
695     if( offsetLast > 0 )
696       {
697       JacobianPositionType jacobianWithRespectToPosition;
698       transform->ComputeJacobianWithRespectToPosition(transformedPoint, jacobianWithRespectToPosition);
699 
700       // Perform the following matrix multiplication in-place:
701       // outJacobian[0:NDimensions,0:offsetLast] = jacobianWithRespectToPosition*outJacobian[0:NDimensions,0:offsetLast]
702       assert(jacobianWithRespectToPosition.rows() == NDimensions);
703       double temp[NDimensions];
704       for (unsigned int c = 0; c < offsetLast; ++c)
705         {
706         for (unsigned int r = 0; r < NDimensions; ++r)
707           {
708           temp[r] = 0.0;
709           //update_j[r][c] = 0.0;
710           for (unsigned int k = 0; k < NDimensions; ++k)
711             {
712             //update_j[r][c]
713             temp[r] += jacobianWithRespectToPosition[r][k] * outJacobian[k][c];
714             }
715           }
716         for (unsigned int r = 0; r < NDimensions; ++r)
717           {
718           outJacobian[r][c] = temp[r];
719           }
720         }
721       }
722 
723     /* Transform the point so it's ready for next transform's Jacobian */
724     transformedPoint = transform->TransformPoint( transformedPoint );
725     }
726 }
727 
728 
729 template<typename TParametersValueType, unsigned int NDimensions>
730 const typename CompositeTransform<TParametersValueType, NDimensions>::ParametersType &
731 CompositeTransform<TParametersValueType, NDimensions>
GetParameters() const732 ::GetParameters() const
733 {
734   const TransformQueueType & transforms = this->GetTransformsToOptimizeQueue();
735   if( transforms.size() == 1 )
736     {
737     // Return directly to avoid copying. Most often we'll have only a single
738     // active transform, so we'll end up here.
739     return transforms[0]->GetParameters();
740     }
741   else
742     {
743     /* Resize destructively. But if it's already this size, nothing is done so
744          * it's efficient. */
745     this->m_Parameters.SetSize( this->GetNumberOfParameters() );
746 
747     NumberOfParametersType offset = NumericTraits< NumberOfParametersType >::ZeroValue();
748 
749     auto it = transforms.end();
750 
751     do
752       {
753       it--;
754       const ParametersType & subParameters = (*it)->GetParameters();
755       /* use vnl_vector data_block() to get data ptr */
756       std::copy(subParameters.data_block(),
757                 subParameters.data_block()+subParameters.Size(),
758                 &(this->m_Parameters.data_block() )[offset]);
759       offset += subParameters.Size();
760 
761       }
762     while( it != transforms.begin() );
763     }
764 
765   return this->m_Parameters;
766 }
767 
768 
769 template<typename TParametersValueType, unsigned int NDimensions>
770 void
771 CompositeTransform<TParametersValueType, NDimensions>
SetParameters(const ParametersType & inputParameters)772 ::SetParameters(const ParametersType & inputParameters)
773 {
774   /* We do not copy inputParameters into m_Parameters,
775      * to avoid unnecessary copying. */
776 
777   /* Assumes input params are concatenation of the parameters of the
778      sub transforms currently selected for optimization, in
779      the order of the queue from begin() to end(). */
780   TransformQueueType transforms = this->GetTransformsToOptimizeQueue();
781 
782   /* Verify proper input size. */
783   if( inputParameters.Size() != this->GetNumberOfParameters() )
784     {
785     itkExceptionMacro(<< "Input parameter list size is not expected size. "
786                       << inputParameters.Size() << " instead of "
787                       << this->GetNumberOfParameters() << ".");
788     }
789 
790   if( transforms.size() == 1 )
791     {
792     /* Avoid unnecessary copying. See comments below */
793     if( &inputParameters == &this->m_Parameters )
794       {
795       transforms[0]->SetParameters( transforms[0]->GetParameters() );
796       }
797     else
798       {
799       transforms[0]->SetParameters(inputParameters);
800       }
801     }
802   else
803     {
804     NumberOfParametersType offset = NumericTraits< NumberOfParametersType >::ZeroValue();
805     auto it = transforms.end();
806 
807     do
808       {
809       it--;
810       /* If inputParams is same object as m_Parameters, we just pass
811        * each sub-transforms own m_Parameters in. This is needed to
812        * avoid unnecessary copying of parameters in the sub-transforms,
813        * while still allowing SetParameters to do any oeprations on the
814        * parameters to update member variable states. A hack. */
815       if( &inputParameters == &this->m_Parameters )
816         {
817         (*it)->SetParameters( (*it)->GetParameters() );
818         }
819       else
820         {
821         const size_t parameterSize = (*it)->GetParameters().Size();
822         (*it)->CopyInParameters(&(inputParameters.data_block() )[offset],
823                                 &(inputParameters.data_block() )[offset]+parameterSize );
824         offset += static_cast<NumberOfParametersType>(parameterSize);
825         }
826 
827       }
828     while( it != transforms.begin() );
829     }
830 }
831 
832 
833 template<typename TParametersValueType, unsigned int NDimensions>
834 const typename CompositeTransform<TParametersValueType, NDimensions>::FixedParametersType &
835 CompositeTransform<TParametersValueType, NDimensions>
GetFixedParameters() const836 ::GetFixedParameters() const
837   {
838   TransformQueueType transforms = this->GetTransformsToOptimizeQueue();
839   /* Resize destructively. But if it's already this size, nothing is done so
840    * it's efficient. */
841   this->m_FixedParameters.SetSize( this->GetNumberOfFixedParameters() );
842 
843   NumberOfParametersType offset = NumericTraits< NumberOfParametersType >::ZeroValue();
844   typename TransformQueueType::const_iterator it;
845 
846   it = transforms.end();
847 
848   do
849     {
850     it--;
851     const FixedParametersType & subFixedParameters = (*it)->GetFixedParameters();
852     /* use vnl_vector data_block() to get data ptr */
853     std::copy(subFixedParameters.data_block(),
854               subFixedParameters.data_block()+subFixedParameters.Size(),
855               &(this->m_FixedParameters.data_block() )[offset]);
856     offset += subFixedParameters.Size();
857     }
858   while( it != transforms.begin() );
859 
860   return this->m_FixedParameters;
861   }
862 
863 template<typename TParametersValueType, unsigned int NDimensions>
864 void
865 CompositeTransform<TParametersValueType, NDimensions>
SetFixedParameters(const FixedParametersType & inputParameters)866 ::SetFixedParameters(const FixedParametersType & inputParameters)
867 {
868   /* Assumes input params are concatenation of the parameters of the
869    * sub transforms currently selected for optimization. */
870   TransformQueueType transforms = this->GetTransformsToOptimizeQueue();
871 
872   NumberOfParametersType offset = NumericTraits< NumberOfParametersType >::ZeroValue();
873 
874 
875   /* Verify proper input size. */
876   if( inputParameters.Size() != this->GetNumberOfFixedParameters() )
877     {
878     itkExceptionMacro(<< "Input parameter list size is not expected size. "
879                       << inputParameters.Size() << " instead of "
880                       << this->GetNumberOfFixedParameters() << ".");
881     }
882   this->m_FixedParameters = inputParameters;
883 
884   typename TransformQueueType::const_iterator it = transforms.end();
885 
886   do
887     {
888     it--;
889     const size_t fixedParameterSize=(*it)->GetFixedParameters().Size();
890     (*it)->CopyInFixedParameters(&(this->m_FixedParameters.data_block() )[offset],
891               &(this->m_FixedParameters.data_block() )[offset]+fixedParameterSize);
892     offset += static_cast<NumberOfParametersType>(fixedParameterSize);
893     }
894   while( it != transforms.begin() );
895 }
896 
897 
898 template<typename TParametersValueType, unsigned int NDimensions>
899 typename CompositeTransform<TParametersValueType, NDimensions>::NumberOfParametersType
900 CompositeTransform<TParametersValueType, NDimensions>
GetNumberOfParameters() const901 ::GetNumberOfParameters() const
902 {
903   /* Returns to total number of params in all transforms currently
904    * set to be used for optimized.
905    * NOTE: We might want to optimize this only to store the result and
906    * only re-calc when the composite object has been modified.
907    * However, it seems that number of parameter might change for dense
908    * field transforms (deformation, bspline) during processing and
909    * we wouldn't know that in this class, so this is safest. */
910   NumberOfParametersType result = NumericTraits< NumberOfParametersType >::ZeroValue();
911 
912 
913   for( signed long tind = (signed long) this->GetNumberOfTransforms() - 1; tind >= 0; tind-- )
914     {
915     if( this->GetNthTransformToOptimize( tind ) )
916       {
917       const TransformType * transform = this->GetNthTransformConstPointer( tind );
918       result += transform->GetNumberOfParameters();
919       }
920     }
921   return result;
922 }
923 
924 
925 template<typename TParametersValueType, unsigned int NDimensions>
926 typename CompositeTransform<TParametersValueType, NDimensions>::NumberOfParametersType
927 CompositeTransform<TParametersValueType, NDimensions>
GetNumberOfLocalParameters() const928 ::GetNumberOfLocalParameters() const
929 {
930   if ( this->GetMTime() == this->m_LocalParametersUpdateTime )
931    {
932    return this->m_NumberOfLocalParameters;
933    }
934 
935   this->m_LocalParametersUpdateTime = this->GetMTime();
936 
937   /* Returns to total number of *local* params in all transforms currently
938    * set to be used for optimized.
939    * Note that unlike in GetNumberOfParameters(), we don't expect the
940    * number of local parameters to possibly change. */
941   NumberOfParametersType result = NumericTraits< NumberOfParametersType >::ZeroValue();
942 
943   for( signed long tind = (signed long) this->GetNumberOfTransforms() - 1; tind >= 0; tind-- )
944     {
945     if( this->GetNthTransformToOptimize( tind ) )
946       {
947       const TransformType * transform = this->GetNthTransformConstPointer( tind );
948       result += transform->GetNumberOfLocalParameters();
949       }
950     }
951   this->m_NumberOfLocalParameters = result;
952   return result;
953 }
954 
955 
956 template<typename TParametersValueType, unsigned int NDimensions>
957 typename CompositeTransform<TParametersValueType, NDimensions>::NumberOfParametersType
958 CompositeTransform<TParametersValueType, NDimensions>
GetNumberOfFixedParameters() const959 ::GetNumberOfFixedParameters() const
960 {
961   /* Returns to total number of params in all transforms currently
962    * set to be used for optimized.
963    * NOTE: We might want to optimize this only to store the result and
964    * only re-calc when the composite object has been modified. */
965   NumberOfParametersType result = NumericTraits< NumberOfParametersType >::ZeroValue();
966 
967   for( signed long tind = (signed long) this->GetNumberOfTransforms() - 1;
968        tind >= 0; tind-- )
969     {
970     if( this->GetNthTransformToOptimize( tind ) )
971       {
972       const TransformType * transform = this->GetNthTransformConstPointer( tind );
973       result += transform->GetFixedParameters().Size();
974       }
975     }
976   return result;
977 }
978 
979 
980 template<typename TParametersValueType, unsigned int NDimensions>
981 void
982 CompositeTransform<TParametersValueType, NDimensions>
UpdateTransformParameters(const DerivativeType & update,ScalarType factor)983 ::UpdateTransformParameters( const DerivativeType & update, ScalarType  factor )
984 {
985   /* Update parameters within the sub-transforms set to be optimized. */
986   /* NOTE: We might want to thread this over each sub-transform, if we
987    * find we're working with longer lists of sub-transforms that do
988    * not implement any threading of their own for UpdateTransformParameters.
989    * Since the plan is for a UpdateTransformParameters functor that is
990    * user-assignable, we would need a method in the
991    * functor to return whether or not it does therading. If all sub-transforms
992    * return that they don't thread, we could do each sub-transform in its
993    * own thread from here. */
994   NumberOfParametersType numberOfParameters = this->GetNumberOfParameters();
995 
996   if( update.Size() != numberOfParameters )
997     {
998     itkExceptionMacro("Parameter update size, " << update.Size() << ", must "
999                       " be same as transform parameter size, " << numberOfParameters << std::endl);
1000     }
1001 
1002   NumberOfParametersType offset = NumericTraits< NumberOfParametersType >::ZeroValue();
1003 
1004 
1005   for( signed long tind = (signed long) this->GetNumberOfTransforms() - 1;
1006        tind >= 0; tind-- )
1007     {
1008     if( this->GetNthTransformToOptimize( tind ) )
1009       {
1010       TransformType * subtransform = this->GetNthTransformModifiablePointer( tind );
1011       /* The input values are in a monolithic block, so we have to point
1012        * to the subregion corresponding to the individual subtransform.
1013        * This simply creates an Array object with data pointer, no
1014        * memory is allocated or copied.
1015        * NOTE: the use of const_cast is used to avoid a deep copy in the underlying vnl_vector
1016        * by using LetArrayManageMemory=false, and being very careful here we can
1017        * ensure that casting away consteness does not result in memory corruption. */
1018       auto * nonConstDataRefForPerformance = const_cast< typename DerivativeType::ValueType * >(
1019                                                &( (update.data_block() )[offset]) );
1020       const DerivativeType subUpdate( nonConstDataRefForPerformance,
1021                                 subtransform->GetNumberOfParameters(), false );
1022       /* This call will also call SetParameters, so don't need to call it
1023        * expliclity here. */
1024       subtransform->UpdateTransformParameters( subUpdate, factor );
1025       offset += subtransform->GetNumberOfParameters();
1026       }
1027     }
1028   this->Modified();
1029 }
1030 
1031 
1032 template<typename TParametersValueType, unsigned int NDimensions>
1033 typename CompositeTransform<TParametersValueType, NDimensions>::TransformQueueType &
1034 CompositeTransform<TParametersValueType, NDimensions>
GetTransformsToOptimizeQueue() const1035 ::GetTransformsToOptimizeQueue() const
1036 {
1037   /* Update the list of transforms to use for optimization only if
1038    the selection of transforms to optimize may have changed */
1039   if( this->GetMTime() > this->m_PreviousTransformsToOptimizeUpdateTime )
1040     {
1041     this->m_TransformsToOptimizeQueue.clear();
1042     for( size_t n = 0; n < this->m_TransformQueue.size(); n++ )
1043       {
1044       /* Return them in the same order as they're found in the main list */
1045       if( this->GetNthTransformToOptimize(static_cast<SizeValueType>( n ) ) )
1046         {
1047         this->m_TransformsToOptimizeQueue.push_back( this->GetNthTransformModifiablePointer(static_cast<SizeValueType>( n ) ) );
1048         }
1049       }
1050     this->m_PreviousTransformsToOptimizeUpdateTime = this->GetMTime();
1051     }
1052   return this->m_TransformsToOptimizeQueue;
1053 }
1054 
1055 
1056 template<typename TParametersValueType, unsigned int NDimensions>
1057 void
1058 CompositeTransform<TParametersValueType, NDimensions>
FlattenTransformQueue()1059 ::FlattenTransformQueue()
1060 {
1061   TransformQueueType             transformQueue;
1062   TransformQueueType             transformsToOptimizeQueue;
1063   TransformsToOptimizeFlagsType  transformsToOptimizeFlags;
1064 
1065   for( SizeValueType m = 0; m < this->GetNumberOfTransforms(); m++ )
1066     {
1067     auto * nestedCompositeTransform = dynamic_cast<Self *>( this->m_TransformQueue[m].GetPointer() );
1068     if( nestedCompositeTransform )
1069       {
1070       nestedCompositeTransform->FlattenTransformQueue();
1071       for( SizeValueType n = 0; n < nestedCompositeTransform->GetNumberOfTransforms(); n++ )
1072         {
1073         transformQueue.push_back( nestedCompositeTransform->GetNthTransformModifiablePointer( n ) );
1074         if( nestedCompositeTransform->GetNthTransformToOptimize( n ) )
1075           {
1076           transformsToOptimizeFlags.push_back( true );
1077           transformsToOptimizeQueue.push_back( nestedCompositeTransform->GetNthTransformModifiablePointer( n ) );
1078           }
1079         else
1080           {
1081           transformsToOptimizeFlags.push_back( false );
1082           }
1083         }
1084       }
1085     else
1086       {
1087       transformQueue.push_back( this->m_TransformQueue[m] );
1088       if( this->m_TransformsToOptimizeFlags[m] )
1089         {
1090         transformsToOptimizeFlags.push_back( true );
1091         transformsToOptimizeQueue.push_back( this->m_TransformQueue[m] );
1092         }
1093       else
1094         {
1095         transformsToOptimizeFlags.push_back( false );
1096         }
1097       }
1098     }
1099 
1100   this->m_TransformQueue = transformQueue;
1101   this->m_TransformsToOptimizeQueue = transformsToOptimizeQueue;
1102   this->m_TransformsToOptimizeFlags = transformsToOptimizeFlags;
1103 }
1104 
1105 
1106 template<typename TParametersValueType, unsigned int NDimensions>
1107 void
1108 CompositeTransform<TParametersValueType, NDimensions>
PrintSelf(std::ostream & os,Indent indent) const1109 ::PrintSelf( std::ostream& os, Indent indent ) const
1110 {
1111   Superclass::PrintSelf( os, indent );
1112 
1113   if( this->GetNumberOfTransforms() == 0 )
1114     {
1115     return;
1116     }
1117 
1118   os << indent << "TransformsToOptimizeFlags, begin() to end(): " << std::endl << indent << indent;
1119   for(  auto
1120         it = this->m_TransformsToOptimizeFlags.begin();
1121         it != this->m_TransformsToOptimizeFlags.end(); it++ )
1122     {
1123     os << *it << " ";
1124     }
1125   os << std::endl;
1126 
1127   os << indent <<  "TransformsToOptimize in queue, from begin to end:" << std::endl;
1128   typename TransformQueueType::const_iterator cit;
1129   for( cit = this->m_TransformsToOptimizeQueue.begin();
1130        cit != this->m_TransformsToOptimizeQueue.end(); ++cit )
1131     {
1132     os << indent << ">>>>>>>>>" << std::endl;
1133     (*cit)->Print( os, indent );
1134     }
1135   os << indent <<  "End of TransformsToOptimizeQueue." << std::endl << "<<<<<<<<<<" << std::endl;
1136 
1137   os << indent <<  "End of CompositeTransform." << std::endl << "<<<<<<<<<<" << std::endl;
1138 }
1139 
1140 
1141 template<typename TParametersValueType, unsigned int NDimensions>
1142 typename LightObject::Pointer
1143 CompositeTransform<TParametersValueType, NDimensions>
InternalClone() const1144 ::InternalClone() const
1145 {
1146   // This class doesn't use its superclass implementation
1147   // TODO: is it really the right behavior?
1148   // LightObject::Pointer loPtr = Superclass::InternalClone();
1149 
1150   LightObject::Pointer loPtr = CreateAnother();
1151   typename Self::Pointer clone =
1152     dynamic_cast<Self *>(loPtr.GetPointer());
1153   if(clone.IsNull())
1154     {
1155     itkExceptionMacro(<< "downcast to type " << this->GetNameOfClass() << " failed.");
1156     }
1157 
1158   auto tqIt = this->m_TransformQueue.begin();
1159 
1160   auto tfIt = this->m_TransformsToOptimizeFlags.begin();
1161 
1162   for(int i = 0; tqIt != this->m_TransformQueue.end() &&
1163         tfIt != this->m_TransformsToOptimizeFlags.end();
1164       ++tqIt, ++tfIt, ++i)
1165     {
1166     clone->AddTransform((*tqIt)->Clone().GetPointer());
1167     clone->SetNthTransformToOptimize(i,(*tfIt));
1168     }
1169   return loPtr;
1170 }
1171 
1172 } // end namespace itk
1173 
1174 #endif
1175