]> Creatis software - clitk.git/blob - registration/clitkBSplineDeformableTransform.h
Moved from repository clitk to clitk.private/tests_dav
[clitk.git] / registration / clitkBSplineDeformableTransform.h
1 /*=========================================================================
2   Program:   vv                     http://www.creatis.insa-lyon.fr/rio/vv
3
4   Authors belong to: 
5   - University of LYON              http://www.universite-lyon.fr/
6   - Léon Bérard cancer center       http://www.centreleonberard.fr
7   - CREATIS CNRS laboratory         http://www.creatis.insa-lyon.fr
8
9   This software is distributed WITHOUT ANY WARRANTY; without even
10   the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
11   PURPOSE.  See the copyright notices for more information.
12
13   It is distributed under dual licence
14
15   - BSD        See included LICENSE.txt file
16   - CeCILL-B   http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
17 ===========================================================================**/
18 #ifndef __clitkBSplineDeformableTransform_h
19 #define __clitkBSplineDeformableTransform_h
20 #include "clitkVectorBSplineResampleImageFunctionWithLUT.h"
21
22 //itk include
23 #include "itkTransform.h"
24 #include "itkImage.h"
25 #include "itkImageRegion.h"
26 #include "itkSpatialObject.h"
27
28 namespace clitk
29 {
30   // Forward declaration needed for friendship
31   template <class TCoordRep, unsigned int NInputDimensions, unsigned int NOutputDimensions>
32   class ITK_EXPORT MultipleBSplineDeformableTransform;
33
34   template <
35     class TCoordRep = double,               // Data type for scalars, coordinate representation,vectors
36     unsigned int NInputDimensions = 3,        // Number of input dimensions
37     unsigned int NOutputDimensions = 3 >      // Number of output dimensions
38   class ITK_EXPORT BSplineDeformableTransform : 
39     public itk::Transform< TCoordRep, NInputDimensions, NOutputDimensions >
40   {
41   public:
42
43     //====================================================================
44     // Typedefs
45     //====================================================================
46     typedef BSplineDeformableTransform                         Self;
47     typedef itk::Transform< TCoordRep, NInputDimensions, NOutputDimensions > Superclass;
48     typedef itk::SmartPointer<Self>                                 Pointer;
49     typedef itk::SmartPointer<const Self>                           ConstPointer;
50       
51     /** New macro for creation of through the object factory.*/
52     itkNewMacro( Self );
53
54     /** Run-time type information (and related methods). */
55     itkTypeMacro( BSplineDeformableTransform, Transform );
56
57     /** Dimension of the domain space. */
58     itkStaticConstMacro(OutputDimension, unsigned int, NOutputDimensions);
59
60     /** Dimension of the input model. */
61     itkStaticConstMacro(InputDimension, unsigned int, NInputDimensions);
62
63     /** Standard scalar type for this class. */
64     typedef typename Superclass::ScalarType ScalarType;
65
66     /** Standard parameters container. */
67     typedef typename Superclass::ParametersType ParametersType;
68 #if ITK_VERSION_MAJOR >= 4
69     typedef typename Superclass::NumberOfParametersType NumberOfParametersType;
70 #endif
71
72     /** Standard Jacobian container. */
73     typedef typename Superclass::JacobianType JacobianType;
74
75     /** Standard vector type for this class. */
76     typedef itk::Vector<TCoordRep,
77                         InputDimension> InputVectorType;
78     typedef itk::Vector<TCoordRep,
79                         OutputDimension> OutputVectorType;
80   
81     /** Standard covariant vector type for this class. */
82     typedef itk::CovariantVector<TCoordRep,
83                                  InputDimension> InputCovariantVectorType;
84     typedef itk::CovariantVector<TCoordRep,
85                                  OutputDimension> OutputCovariantVectorType;
86
87     /** Standard vnl_vector type for this class. */
88     typedef vnl_vector_fixed<TCoordRep,
89                              InputDimension> InputVnlVectorType;
90     typedef vnl_vector_fixed<TCoordRep,
91                              OutputDimension> OutputVnlVectorType;
92
93     /** Standard coordinate point type for this class. */
94     typedef itk::Point<TCoordRep,
95                        InputDimension> InputPointType;
96     typedef itk::Point<TCoordRep,
97                        OutputDimension> OutputPointType;
98
99     //JV Parameters as images with OutputDimension number of components per Pixel
100     typedef typename ParametersType::ValueType                          ParametersValueType;
101     typedef typename itk::Vector<ParametersValueType, OutputDimension>  PixelType;
102     typedef itk::Image<PixelType, InputDimension>                       CoefficientImageType;
103     typedef typename CoefficientImageType::Pointer                      CoefficientImagePointer;
104
105   
106     /** Typedefs for specifying the extend to the grid. */
107     typedef itk::ImageRegion<InputDimension>                RegionType;
108     typedef typename RegionType::IndexType                  IndexType;
109     typedef typename RegionType::SizeType                   SizeType;
110     typedef typename CoefficientImageType::SpacingType      SpacingType;
111     typedef typename CoefficientImageType::DirectionType    DirectionType;
112     typedef typename CoefficientImageType::PointType        OriginType;
113     typedef itk::ContinuousIndex<TCoordRep, InputDimension> ContinuousIndexType;
114   
115     //JV added for the BLUT interpolator
116     //typedef itk::Vector<TCoordRep, InputDimension> OutputSpacingType;
117
118     //JV m_VectorInterpolator
119     typedef VectorBSplineResampleImageFunctionWithLUT
120     <CoefficientImageType, TCoordRep> VectorInterpolatorType;
121     typedef  typename  VectorInterpolatorType::CoefficientDataType CoefficientDataType;
122     typedef  typename  VectorInterpolatorType::CoefficientDataType WeightsDataType;
123
124     /** Typedef of the bulk transform. */
125     typedef itk::Transform<ScalarType, InputDimension, OutputDimension> BulkTransformType;
126     typedef BulkTransformType*                  BulkTransformPointer;
127
128     /** Typedef of the mask */
129     typedef itk::SpatialObject<  InputDimension >   MaskType;
130     typedef MaskType*    MaskPointer;
131
132     //====================================================================
133     // Set et Gets
134     //====================================================================
135     //JV  added for the BLUT interpolator
136     void SetSplineOrder(const unsigned int &  splineOrder);
137     void SetSplineOrders(const SizeType &  splineOrders);
138     itkGetMacro( SplineOrders, SizeType );
139     itkGetConstMacro( SplineOrders, SizeType );
140     void SetLUTSamplingFactor(const int &  samplingFactor);
141     void SetLUTSamplingFactors(const SizeType &  samplingFactors);
142     itkGetMacro( LUTSamplingFactors, SizeType );
143     itkGetConstMacro( LUTSamplingFactors,SizeType );
144     //void SetOutputSpacing(const OutputSpacingType & outputSpacing);
145     //itkGetMacro( OutputSpacing, OutputSpacingType );
146     //itkGetConstMacro( OutputSpacing, OutputSpacingType );
147
148     void SetParameters(const ParametersType & parameters);
149   
150     void SetFixedParameters(const ParametersType & parameters);
151
152     void SetParametersByValue(const ParametersType & parameters);
153
154     void SetIdentity();
155
156     /** Get the Transformation Parameters. */
157     virtual const ParametersType& GetParameters(void) const;
158
159     /** Get the Transformation Fixed Parameters. */
160     virtual const ParametersType& GetFixedParameters(void) const;
161   
162     // The coefficientImage
163     virtual CoefficientImagePointer GetCoefficientImage()
164     { return m_CoefficientImage; }
165     virtual const CoefficientImagePointer  GetCoefficientImage() const
166     { return m_CoefficientImage; }
167     virtual void SetCoefficientImage(CoefficientImagePointer image);  
168
169     /** This method specifies the region over which the grid resides. */
170     virtual void SetGridRegion( const RegionType& region );
171     itkGetMacro( GridRegion, RegionType );
172     itkGetConstMacro( GridRegion, RegionType );
173
174     /** This method specifies the grid spacing or resolution. */
175     virtual void SetGridSpacing( const SpacingType& spacing );
176     itkGetMacro( GridSpacing, SpacingType );
177     itkGetConstMacro( GridSpacing, SpacingType );
178
179     /** This method specifies the grid directions . */
180     virtual void SetGridDirection( const DirectionType & spacing );
181     itkGetMacro( GridDirection, DirectionType );
182     itkGetConstMacro( GridDirection, DirectionType );
183
184     /** This method specifies the grid origin. */
185     virtual void SetGridOrigin( const OriginType& origin );
186     itkGetMacro( GridOrigin, OriginType );
187     itkGetConstMacro( GridOrigin, OriginType );
188     
189     // Set the bulk transform, real pointer
190     // itkSetConstObjectMacro( BulkTransform, BulkTransformType );
191     // itkGetConstObjectMacro( BulkTransform, BulkTransformType );
192     void SetBulkTransform(BulkTransformPointer b){m_BulkTransform=b;}
193     BulkTransformPointer GetBulkTransform(void) {return m_BulkTransform;}
194
195     // Set mask, inside transform applies, outside zero, real pointer
196     void SetMask(MaskPointer m){m_Mask=m;}
197     MaskPointer GetMask(void){return m_Mask;}
198     // itkSetConstObjectMacro( Mask, MaskType );
199     // itkGetConstObjectMacro( Mask, MaskType );
200
201     /** Transform points by a BSpline deformable transformation. */
202     OutputPointType  TransformPoint(const InputPointType  &point ) const;
203   
204     // JV added for just the deformable part, without bulk
205     OutputPointType  DeformablyTransformPoint(const InputPointType  &point ) const;
206   
207     /** Parameter index array type. */
208     typedef itk::Array<unsigned long> ParameterIndexArrayType;
209
210     /** Transform points by a BSpline deformable transformation. 
211      * On return, weights contains the interpolation weights used to compute the 
212      * deformation and indices of the x (zeroth) dimension coefficient parameters
213      * in the support region used to compute the deformation.
214      * Parameter indices for the i-th dimension can be obtained by adding
215      * ( i * this->GetNumberOfParametersPerDimension() ) to the indices array.
216      */
217
218     // JV not implemented
219     //  virtual void TransformPoint( const InputPointType & inputPoint,
220     //                           OutputPointType & outputPoint,
221     //                           WeightsType & weights,
222     //                           ParameterIndexArrayType & indices, 
223     //                           bool & inside ) const;
224     //     virtual void DeformablyTransformPoint( const InputPointType & inputPoint,
225     //                                     OutputPointType & outputPoint,
226     //                                     WeightsType & weights,
227     //                                     ParameterIndexArrayType & indices, 
228     //                                     bool & inside ) const;
229     //     virtual void GetJacobian( const InputPointType & inputPoint,
230     //                                WeightsType & weights,
231     //                                ParameterIndexArrayType & indices
232     //                                ) const;
233    
234     /** Method to transform a vector - 
235      *  not applicable for this type of transform. */
236     virtual OutputVectorType TransformVector(const InputVectorType &) const
237     { 
238       itkExceptionMacro(<< "Method not applicable for deformable transform." );
239       return OutputVectorType(); 
240     }
241
242     /** Method to transform a vnl_vector - 
243      *  not applicable for this type of transform */
244     virtual OutputVnlVectorType TransformVector(const InputVnlVectorType &) const
245     { 
246       itkExceptionMacro(<< "Method not applicable for deformable transform. ");
247       return OutputVnlVectorType(); 
248     }
249
250     /** Method to transform a CovariantVector - 
251      *  not applicable for this type of transform */
252     virtual OutputCovariantVectorType TransformCovariantVector(
253                                                                const InputCovariantVectorType &) const
254     { 
255       itkExceptionMacro(<< "Method not applicable for deformable transfrom. ");
256       return OutputCovariantVectorType(); 
257     } 
258     
259     /** Compute the Jacobian Matrix of the transformation at one point */
260 #if ITK_VERSION_MAJOR >= 4
261     virtual void ComputeJacobianWithRespectToParameters (const InputPointType &p, JacobianType &jacobian) const;
262     virtual void ComputeJacobianWithRespectToPosition (const InputPointType &p, JacobianType &jacobian) const
263     {
264       itkExceptionMacro( "ComputeJacobianWithRespectToPosition not yet implemented for " << this->GetNameOfClass() );
265     }
266 #else
267     virtual const JacobianType& GetJacobian(const InputPointType  &point ) const;
268 #endif
269
270     /** Return the number of parameters that completely define the Transfom */
271 #if ITK_VERSION_MAJOR >= 4
272     virtual NumberOfParametersType GetNumberOfParameters(void) const;
273 #else
274     virtual unsigned int GetNumberOfParameters(void) const;
275 #endif
276
277     /** Return the number of parameters per dimension */
278     unsigned int GetNumberOfParametersPerDimension(void) const;
279
280     /** Return the region of the grid wholly within the support region */
281     itkGetConstReferenceMacro( ValidRegion, RegionType );
282
283     /** Indicates that this transform is linear. That is, given two
284      * points P and Q, and scalar coefficients a and b, then
285      *
286      *           T( a*P + b*Q ) = a * T(P) + b * T(Q)
287      */
288     virtual bool IsLinear() const { return false; }
289
290    //unsigned int GetNumberOfAffectedWeights() const;
291
292   protected:
293     /** Print contents of an BSplineDeformableTransform. */
294     void PrintSelf(std::ostream &os, itk::Indent indent) const;
295
296
297     BSplineDeformableTransform();
298     virtual ~BSplineDeformableTransform();
299
300     /** Wrap flat array into images of coefficients. */
301     void WrapAsImages();
302
303     /** Convert an input point to a continuous index inside the BSpline grid */
304     void TransformPointToContinuousIndex( const InputPointType & point, ContinuousIndexType & index ) const;
305
306   private:
307     BSplineDeformableTransform(const Self&); //purposely not implemented
308     void operator=(const Self&); //purposely not implemented
309
310     /** The bulk transform. */
311     BulkTransformPointer  m_BulkTransform;
312     MaskPointer m_Mask;
313
314     // JV added for BLUT interpolator
315     SizeType m_SplineOrders;
316     SizeType m_LUTSamplingFactors;
317
318     /** Variables defining the coefficient grid extend. */
319     RegionType    m_GridRegion;
320     SpacingType   m_GridSpacing;
321     DirectionType m_GridDirection;
322     OriginType    m_GridOrigin;
323
324     DirectionType m_PointToIndex;
325     DirectionType m_IndexToPoint;
326   
327     RegionType    m_ValidRegion;
328
329     /** Variables defining the interpolation support region. */
330     SizeType m_Offset;
331     itk::FixedArray<bool,InputDimension>    m_SplineOrderOdd;
332     SizeType      m_SupportSize;
333     IndexType     m_ValidRegionLast;
334   
335     /** Array holding images wrapped from the flat parameters. */
336     CoefficientImagePointer   m_WrappedImage;
337
338     /** Vector image representing the B-spline coefficients 
339      *  in each dimension. */
340     CoefficientImagePointer   m_CoefficientImage;
341
342     /** Jacobian as OutputDimension number of images. */
343     typedef typename JacobianType::ValueType JacobianValueType;
344     typedef typename itk::Vector<JacobianValueType,OutputDimension> JacobianPixelType;
345     typedef itk::Image<JacobianPixelType, OutputDimension> JacobianImageType;
346     typename JacobianImageType::Pointer m_JacobianImage[OutputDimension];
347     typedef itk::ImageRegionIterator<JacobianImageType> IteratorType;
348
349     /** Keep track of last support region used in computing the Jacobian
350      * for fast resetting of Jacobian to zero.
351      */
352     //JV for J calculation
353     mutable RegionType    m_SupportRegion;
354     mutable IndexType     m_SupportIndex;
355     mutable IndexType m_LastJacobianIndex;
356     mutable IteratorType m_Iterator[OutputDimension];
357     mutable JacobianPixelType m_ZeroVector;
358     mutable ContinuousIndexType m_Index;
359     mutable bool m_NeedResetJacobian;
360
361     /** Keep a pointer to the input parameters. */
362     const ParametersType *  m_InputParametersPointer;
363
364     /** Internal parameters buffer. */
365     ParametersType          m_InternalParametersBuffer;
366
367     //JV  the BLUT interpolator
368     typename VectorInterpolatorType::Pointer m_VectorInterpolator;
369   
370     /** Check if a continuous index is inside the valid region. */
371     bool InsideValidRegion( const ContinuousIndexType& index ) const;
372
373     // VD Use external data container for JacobianImage
374     unsigned SetJacobianImageData(JacobianPixelType * jacobianDataPointer, unsigned dim);
375
376     // VD Reset Jacobian
377     void ResetJacobian() const;
378
379     // VD Add MultipleBSplineDeformableTransform as friend to facilitate wrapping
380     friend class MultipleBSplineDeformableTransform<TCoordRep, NInputDimensions, NOutputDimensions>;
381 #if ITK_VERSION_MAJOR >= 4
382     mutable JacobianType                            m_SharedDataBSplineJacobian;
383 #endif
384
385   }; //class BSplineDeformableTransform
386
387
388 }  // namespace itk
389
390 #if ITK_TEMPLATE_TXX
391 # include "clitkBSplineDeformableTransform.txx"
392 #endif
393
394
395 #endif // __clitkBSplineDeformableTransform_h