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