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_h
19 #define itkCompositeTransform_h
20 
21 #include "itkMultiTransform.h"
22 
23 #include <deque>
24 
25 namespace itk
26 {
27 
28 /** \class CompositeTransform
29  * \brief This class contains a list of transforms and concatenates them by composition.
30  *
31  * This class concatenates transforms in \b reverse \b queue order by means of composition:
32  *    \f$ T_0 o T_1 = T_0(T_1(x)) \f$
33  * Transforms are stored in a container (queue), in the following order:
34  *    \f$ T_0, T_1, ... , T_N-1 \f$
35  * Transforms are added via a single method, AddTransform(). This adds the
36  * transforms to the back of the queue. A single method for adding transforms
37  * is meant to simplify the interface and prevent errors.
38  * One use of the class is to optimize only a subset of included transforms.
39  *
40  * The sub transforms are the same dimensionality as this class.
41  *
42  * Example:
43  * A user wants to optimize two Affine transforms together, then add a
44  * Deformation Field (DF) transform, and optimize it separately.
45  * He first adds the two Affines, then runs the optimization and both Affines
46  * transforms are optimized. Next, he adds the DF transform and calls
47  * SetOnlyMostRecentTransformToOptimizeOn, which clears the optimization flags
48  * for both of the affine transforms, and leaves the flag set only for the DF
49  * transform, since it was the last transform added. Now he runs the
50  * optimization and only the DF transform is optimized, but the affines are
51  * included in the transformation during the optimization.
52  *
53  * Optimization Flags:
54  * The m_TransformsToOptimize flags hold one flag for each transform in the
55  * queue, designating if each transform is to be used for optimization. Note
56  * that all transforms in the queue are applied in TransformPoint, regardless
57  * of these flags states'. The methods GetParameters, SetParameters,
58  * ComputeJacobianWithRespectToParameters, GetTransformCategory,
59  * GetFixedParameters, and SetFixedParameters all query these
60  * flags and include only those transforms whose corresponding flag is set.
61  * Their input or output is a concatenated array of all transforms set for use
62  * in optimization. The goal is to be able to optimize multiple transforms at
63  * once, while leaving other transforms fixed. See the above example.
64  *
65  * Setting Optimization Flags:
66  * A transform's optimization flag is set when it is added to the queue, and
67  * remains set as other transforms are added. The methods
68  * SetNthTransformToOptimize* and SetAllTransformToOptimize* are used to
69  * set and clear flags arbitrarily. SetOnlyMostRecentTransformToOptimizeOn is
70  * a convenience method for setting only the most recently added transform
71  * for optimization, with the idea that this will be a common practice.
72  *
73  * Indexing:
74  * The index values used in GetNthTransform and
75  * SetNthTransformToOptimize* and SetAllTransformToOptimize* follow the
76  * order in which transforms were added. Thus, the first transform added is at
77  * index 0, the next at index 1, etc.
78  *
79  * Inverse:
80  * The inverse transform is created by retrieving the inverse from each
81  * sub transform and adding them to a composite transform in reverse order.
82  * The m_TransformsToOptimizeFlags is copied in reverse for the inverse.
83  *
84  * \ingroup ITKTransform
85  */
86 template<typename TParametersValueType=double, unsigned int NDimensions = 3>
87 class ITK_TEMPLATE_EXPORT CompositeTransform :
88   public MultiTransform<TParametersValueType, NDimensions, NDimensions>
89 {
90 public:
91   ITK_DISALLOW_COPY_AND_ASSIGN(CompositeTransform);
92 
93   /** Standard class type aliases. */
94   using Self = CompositeTransform;
95   using Superclass = MultiTransform<TParametersValueType, NDimensions, NDimensions>;
96   using Pointer = SmartPointer<Self>;
97   using ConstPointer = SmartPointer<const Self>;
98 
99   /** Run-time type information (and related methods). */
100   itkTypeMacro( CompositeTransform, Transform );
101 
102   /** New macro for creation of through a Smart Pointer */
103   itkNewMacro( Self );
104 
105   /** Sub transform type **/
106   using TransformType = typename Superclass::TransformType;
107   using TransformTypePointer = typename Superclass::TransformTypePointer;
108   /** InverseTransform type. */
109   using InverseTransformBasePointer = typename Superclass::InverseTransformBasePointer;
110   /** Scalar type. */
111   using ScalarType = typename Superclass::ScalarType;
112   /** Parameters type. */
113   using FixedParametersType = typename Superclass::FixedParametersType;
114   using FixedParametersValueType = typename Superclass::FixedParametersValueType;
115   using ParametersType = typename Superclass::ParametersType;
116   using ParametersValueType = typename Superclass::ParametersValueType;
117   /** Derivative type */
118   using DerivativeType = typename Superclass::DerivativeType;
119   /** Jacobian types. */
120   using JacobianType = typename Superclass::JacobianType;
121   using JacobianPositionType = typename Superclass::JacobianPositionType;
122   using InverseJacobianPositionType = typename Superclass::InverseJacobianPositionType;
123   /** Transform category type. */
124   using TransformCategoryType = typename Superclass::TransformCategoryType;
125   /** Standard coordinate point type for this class. */
126   using InputPointType = typename Superclass::InputPointType;
127   using OutputPointType = typename Superclass::OutputPointType;
128   /** Standard vector type for this class. */
129   using InputVectorType = typename Superclass::InputVectorType;
130   using OutputVectorType = typename Superclass::OutputVectorType;
131   /** Standard covariant vector type for this class */
132   using InputCovariantVectorType = typename Superclass::InputCovariantVectorType;
133   using OutputCovariantVectorType = typename Superclass::OutputCovariantVectorType;
134   /** Standard vnl_vector type for this class. */
135   using InputVnlVectorType = typename Superclass::InputVnlVectorType;
136   using OutputVnlVectorType = typename Superclass::OutputVnlVectorType;
137   /** Standard Vectorpixel type for this class */
138   using InputVectorPixelType = typename Superclass::InputVectorPixelType;
139   using OutputVectorPixelType = typename Superclass::OutputVectorPixelType;
140   /** Standard DiffusionTensor3D type alias for this class */
141   using InputDiffusionTensor3DType = typename Superclass::InputDiffusionTensor3DType;
142   using OutputDiffusionTensor3DType = typename Superclass::OutputDiffusionTensor3DType;
143   /** Standard SymmetricSecondRankTensor type alias for this class */
144   using InputSymmetricSecondRankTensorType = typename Superclass::InputSymmetricSecondRankTensorType;
145   using OutputSymmetricSecondRankTensorType = typename Superclass::OutputSymmetricSecondRankTensorType;
146 
147   /** Transform queue type */
148   using TransformQueueType = typename Superclass::TransformQueueType;
149 
150   /** The number of parameters defininig this transform. */
151   using NumberOfParametersType = typename Superclass::NumberOfParametersType;
152 
153   /** Optimization flags queue type */
154   using TransformsToOptimizeFlagsType = std::deque<bool>;
155 
156   /** Dimension of the domain spaces. */
157   static constexpr unsigned int InputDimension = NDimensions;
158   static constexpr unsigned int OutputDimension = NDimensions;
159 
160   /** Active Transform state manipulation */
161 
SetNthTransformToOptimize(SizeValueType i,bool state)162   virtual void SetNthTransformToOptimize( SizeValueType i, bool state )
163   {
164     this->m_TransformsToOptimizeFlags.at(i) = state;
165     this->Modified();
166   }
167 
SetNthTransformToOptimizeOn(SizeValueType i)168   virtual void SetNthTransformToOptimizeOn( SizeValueType i )
169   {
170     this->SetNthTransformToOptimize( i, true );
171   }
172 
SetNthTransformToOptimizeOff(SizeValueType i)173   virtual void SetNthTransformToOptimizeOff( SizeValueType i )
174   {
175     this->SetNthTransformToOptimize( i, false );
176   }
177 
SetAllTransformsToOptimize(bool state)178   virtual void SetAllTransformsToOptimize( bool state )
179   {
180     this->m_TransformsToOptimizeFlags.assign(
181       this->m_TransformsToOptimizeFlags.size(), state );
182     this->Modified();
183   }
184 
SetAllTransformsToOptimizeOn()185   virtual void SetAllTransformsToOptimizeOn()
186   {
187     this->SetAllTransformsToOptimize( true );
188   }
189 
SetAllTransformsToOptimizeOff()190   virtual void SetAllTransformsToOptimizeOff()
191   {
192     this->SetAllTransformsToOptimize( false );
193   }
194 
195   /* With AddTransform() as the only way to add a transform, we
196    * can have this method to easily allow user to optimize only
197    * the transform added most recenlty. */
SetOnlyMostRecentTransformToOptimizeOn()198   virtual void SetOnlyMostRecentTransformToOptimizeOn()
199   {
200     this->SetAllTransformsToOptimize( false );
201     this->SetNthTransformToOptimizeOn( this->GetNumberOfTransforms() - 1 );
202   }
203 
204   /* Get whether the Nth transform is set to be optimzied */
205   /* NOTE: ambiguous function name here - are we getting if the Nth transform
206       is set to be optimized, or the Nth of the transforms that are set to be
207       optimized? */
GetNthTransformToOptimize(SizeValueType i)208   virtual bool GetNthTransformToOptimize( SizeValueType i ) const
209   {
210     return this->m_TransformsToOptimizeFlags.at(i);
211   }
212 
213   /** Access optimize flags */
GetTransformsToOptimizeFlags()214   virtual const TransformsToOptimizeFlagsType & GetTransformsToOptimizeFlags() const
215   {
216     return this->m_TransformsToOptimizeFlags;
217   }
218 
ClearTransformQueue()219   void ClearTransformQueue() override
220   {
221     Superclass::ClearTransformQueue();
222     this->m_TransformsToOptimizeFlags.clear();
223   }
224 
225   /** Returns a boolean indicating whether it is possible or not to compute the
226    * inverse of this current Transform. If it is possible, then the inverse of
227    * the transform is returned in the inverseTransform variable passed by the user.
228    * The inverse consists of the inverse of each sub-transform, in the \b reverse order
229    * of the forward transforms. */
230   bool GetInverse( Self *inverse ) const;
231 
232   InverseTransformBasePointer GetInverseTransform() const override;
233 
234   /** Compute the position of point in the new space.
235   *
236   * Transforms are applied starting from the *back* of the
237   * queue. That is, in reverse order of which they were added, in order
238   * to work properly with ResampleFilter.
239   *
240   * Imagine a user wants to apply an Affine transform followed by a Deformation
241   * Field (DF) transform. He adds the Affine, then the DF. Because the user
242   * typically conceptualizes a transformation as being applied from the Moving
243   * image to the Fixed image, this makes intuitive sense. But since the
244   * ResampleFilter expects to transform from the Fixed image to the Moving
245   * image, the transforms are applied in reverse order of addition, i.e. from
246   * the back of the queue, and thus, DF then Affine.
247   */
248   OutputPointType TransformPoint( const InputPointType & inputPoint ) const override;
249 
250   /**  Method to transform a vector. */
251   using Superclass::TransformVector;
252   OutputVectorType TransformVector(const InputVectorType &) const override;
253 
254   OutputVnlVectorType TransformVector(const InputVnlVectorType & inputVector) const override;
255 
256   OutputVectorPixelType TransformVector(const InputVectorPixelType & inputVector ) const override;
257 
258   OutputVectorType TransformVector(const InputVectorType & inputVector,
259                                            const InputPointType & inputPoint ) const override;
260 
261   OutputVnlVectorType TransformVector(const InputVnlVectorType & inputVector,
262                                               const InputPointType & inputPoint ) const override;
263 
264   OutputVectorPixelType TransformVector(const InputVectorPixelType & inputVector,
265                                                 const InputPointType & inputPoint ) const override;
266 
267   /**  Method to transform a CovariantVector. */
268   using Superclass::TransformCovariantVector;
269   OutputCovariantVectorType TransformCovariantVector(const InputCovariantVectorType &) const override;
270 
271   OutputVectorPixelType TransformCovariantVector(const InputVectorPixelType &) const override;
272 
273   OutputCovariantVectorType TransformCovariantVector(const InputCovariantVectorType & inputVector,
274                                                              const InputPointType & inputPoint ) const override;
275 
276   OutputVectorPixelType TransformCovariantVector(const InputVectorPixelType & inputVector,
277                                                          const InputPointType & inputPoint ) const override;
278 
279   /** Method to transform a DiffusionTensor3D */
280   using Superclass::TransformDiffusionTensor3D;
281   OutputDiffusionTensor3DType TransformDiffusionTensor3D(
282     const InputDiffusionTensor3DType & inputTensor) const override;
283 
284   OutputVectorPixelType TransformDiffusionTensor3D(
285     const InputVectorPixelType & inputTensor) const override;
286 
287   OutputDiffusionTensor3DType TransformDiffusionTensor3D(
288     const InputDiffusionTensor3DType & inputTensor,
289     const InputPointType & inputPoint ) const override;
290 
291   OutputVectorPixelType TransformDiffusionTensor3D(
292     const InputVectorPixelType & inputTensor,
293     const InputPointType & inputPoint ) const override;
294 
295   /** Method to transform a SymmetricSecondRankTensor */
296   using Superclass::TransformSymmetricSecondRankTensor;
297   OutputSymmetricSecondRankTensorType TransformSymmetricSecondRankTensor(
298     const InputSymmetricSecondRankTensorType & inputTensor) const override;
299 
300   OutputVectorPixelType TransformSymmetricSecondRankTensor(
301     const InputVectorPixelType & inputTensor) const override;
302 
303   OutputSymmetricSecondRankTensorType TransformSymmetricSecondRankTensor(
304     const InputSymmetricSecondRankTensorType & inputTensor,
305     const InputPointType & inputPoint ) const override;
306 
307   OutputVectorPixelType TransformSymmetricSecondRankTensor(
308     const InputVectorPixelType & inputTensor,
309     const InputPointType & inputPoint ) const override;
310 
311   /** Special handling for composite transform. If all transforms
312    * are linear, then return category Linear. Otherwise if all
313    * transforms set to optimize are DisplacementFields, then
314    * return DisplacementField category. */
315   TransformCategoryType GetTransformCategory() const override;
316 
317   /** Get/Set Parameter functions work on the current list of transforms
318       that are set to be optimized (active) using the
319       'Set[Nth|All]TransformToOptimze' routines.
320       The parameter data from each active transform is
321       concatenated into a single ParametersType object.
322       \note The sub-transforms are read in \b reverse queue order,
323       so the returned array is ordered in the same way. That is,
324       the last sub-transform to be added is returned first in the
325       parameter array. This is the opposite of what's done in the
326       parent MultiTransform class. */
327   const ParametersType & GetParameters() const override;
328 
329   /* SetParameters only for transforms that are set to be optimized
330    * See GetParameters() for parameter ordering. */
331   void  SetParameters(const ParametersType & p) override;
332 
333   /* GetFixedParameters only for transforms that are set to be optimized
334    * See GetParameters() for parameter ordering. */
335   const FixedParametersType & GetFixedParameters() const override;
336 
337   /* SetFixedParameters only for transforms that are set to be optimized.
338    * See GetParameters() for parameter ordering. */
339   void  SetFixedParameters(const FixedParametersType & fixedParameters) override;
340 
341   /* Get total number of parameters for transforms that are set to be
342    * optimized */
343   NumberOfParametersType GetNumberOfParameters() const override;
344 
345   /* Get total number of local parameters for transforms that are set
346    * to be optimized */
347   NumberOfParametersType GetNumberOfLocalParameters() const override;
348 
349   /* Get total number of fixed parameters for transforms that are set
350    * to be optimized */
351   NumberOfParametersType GetNumberOfFixedParameters() const override;
352 
353   /** Update the transform's parameters by the values in \c update.
354    * See GetParameters() for parameter ordering. */
355   void UpdateTransformParameters( const DerivativeType & update, ScalarType factor = 1.0 ) override;
356 
357   /**
358    * Flatten the transform queue such that there are no nested composite transforms.
359    */
360   virtual void FlattenTransformQueue();
361 
362   /**
363    * Compute the Jacobian with respect to the parameters for the compositie
364    * transform using Jacobian rule. See comments in the implementation.
365    */
366   void ComputeJacobianWithRespectToParameters(const InputPointType  & p, JacobianType & j) const override;
367 
368   /**
369    * Expanded interface to Compute the Jacobian with respect to the parameters for the compositie
370    * transform using Jacobian rule. This version takes in temporary
371    * variables to avoid excessive constructions and memory allocations.
372    * NOTE: outJacobian MUST be sized correctly prior to the call;
373    * outJacobian's size should be [NDimensions, this->GetNumberOfLocalParameters() ]
374    * jacobianCache may be resized internally and will be reused between calls
375    */
376   void ComputeJacobianWithRespectToParametersCachedTemporaries( const InputPointType & p,
377                                                                 JacobianType & outJacobian,
378                                                                 JacobianType & cacheJacobian ) const override;
379 
380 protected:
381   CompositeTransform();
382   ~CompositeTransform() override = default;
383   void PrintSelf( std::ostream& os, Indent indent ) const override;
384 
385   /** Clone the current transform */
386   typename LightObject::Pointer InternalClone() const override;
387 
PushFrontTransform(TransformTypePointer t)388   void PushFrontTransform( TransformTypePointer t  ) override
389   {
390     Superclass::PushFrontTransform( t );
391     /* Add element to list of flags, and set true by default */
392     this->m_TransformsToOptimizeFlags.push_front( true );
393   }
394 
PushBackTransform(TransformTypePointer t)395   void PushBackTransform( TransformTypePointer t  ) override
396   {
397     Superclass::PushBackTransform( t );
398     /* Add element to list of flags, and set true by default */
399     this->m_TransformsToOptimizeFlags.push_back( true );
400   }
401 
PopFrontTransform()402   void PopFrontTransform() override
403   {
404     Superclass::PopFrontTransform();
405     this->m_TransformsToOptimizeFlags.pop_front();
406   }
407 
PopBackTransform()408   void PopBackTransform() override
409   {
410     Superclass::PopBackTransform();
411     this->m_TransformsToOptimizeFlags.pop_back();
412   }
413 
414   /** Get a list of transforms to optimize. Helper function. */
415   TransformQueueType & GetTransformsToOptimizeQueue() const;
416 
417   mutable TransformQueueType            m_TransformsToOptimizeQueue;
418   mutable TransformsToOptimizeFlagsType m_TransformsToOptimizeFlags;
419 
420 private:
421   mutable ModifiedTimeType m_PreviousTransformsToOptimizeUpdateTime;
422 
423 };
424 
425 } // end namespace itk
426 
427 #ifndef ITK_MANUAL_INSTANTIATION
428 #include "itkCompositeTransform.hxx"
429 #endif
430 
431 #endif // itkCompositeTransform_h
432