]> Creatis software - clitk.git/blob - registration/clitkShapedBLUTSpatioTemporalDeformableTransform.h
Moved from repository clitk to clitk.private/tests_dav
[clitk.git] / registration / clitkShapedBLUTSpatioTemporalDeformableTransform.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 __clitkShapedBLUTSpatioTemporalDeformableTransform_h
19 #define __clitkShapedBLUTSpatioTemporalDeformableTransform_h
20 #include "clitkVectorBSplineResampleImageFunctionWithLUT.h"
21 #include "clitkExtractImageFilter.h"
22 #include "clitkLinearCombinationImageFilter.h"
23
24 //itk include
25 #include "itkTransform.h"
26 #include "itkImage.h"
27 #include "itkImageRegion.h"
28 #include "itkSpatialObject.h"
29 #include "itkPasteImageFilter.h"
30
31 namespace clitk
32 {
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 ShapedBLUTSpatioTemporalDeformableTransform : 
39     public itk::Transform< TCoordRep, NInputDimensions, NOutputDimensions >
40   {
41   public:
42
43     //====================================================================
44     // Typedefs
45     //====================================================================
46     typedef ShapedBLUTSpatioTemporalDeformableTransform                         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( ShapedBLUTSpatioTemporalDeformableTransform, Transform );
56
57     /** Dimension of the output space. */
58     itkStaticConstMacro(OutputDimension, unsigned int, NOutputDimensions);
59
60     /** Dimension of the input space. */
61     itkStaticConstMacro(InputDimension, unsigned int, NInputDimensions);
62
63     //JV the number of spatial dimensions
64     itkStaticConstMacro(SpaceDimension, unsigned int, NInputDimensions-1);
65
66     /** Standard scalar type for this class. */
67     typedef typename Superclass::ScalarType ScalarType;
68
69     /** Standard parameters container. */
70     typedef typename Superclass::ParametersType ParametersType;
71 #if ITK_VERSION_MAJOR >= 4
72     typedef typename Superclass::NumberOfParametersType NumberOfParametersType;
73 #endif
74
75     /** Standard Jacobian container. */
76     typedef typename Superclass::JacobianType JacobianType;
77
78     /** Standard vector type for this class. */
79     typedef itk::Vector<TCoordRep,
80                         InputDimension> InputVectorType;
81     typedef itk::Vector<TCoordRep,
82                         OutputDimension> OutputVectorType;
83   
84     /** Standard covariant vector type for this class. */
85     typedef itk::CovariantVector<TCoordRep,
86                                  InputDimension> InputCovariantVectorType;
87     typedef itk::CovariantVector<TCoordRep,
88                                  OutputDimension> OutputCovariantVectorType;
89
90     /** Standard vnl_vector type for this class. */
91     typedef vnl_vector_fixed<TCoordRep,
92                              InputDimension> InputVnlVectorType;
93     typedef vnl_vector_fixed<TCoordRep,
94                              OutputDimension> OutputVnlVectorType;
95
96     /** Standard coordinate point type for this class. */
97     typedef itk::Point<TCoordRep,
98                        InputDimension> InputPointType;
99     typedef itk::Point<TCoordRep,
100                        OutputDimension> OutputPointType;
101
102     //JV Parameters as images with OutputDimension number of components per Pixel
103     typedef typename ParametersType::ValueType                          ParametersValueType;
104     typedef typename itk::Vector<ParametersValueType, SpaceDimension>   PixelType;
105     typedef itk::Image<PixelType, InputDimension>                       CoefficientImageType;
106     typedef typename CoefficientImageType::Pointer                      CoefficientImagePointer;
107
108   
109     /** Typedefs for specifying the extend to the grid. */
110     typedef itk::ImageRegion<InputDimension>                RegionType;
111     typedef typename RegionType::IndexType                  IndexType;
112     typedef typename RegionType::SizeType                   SizeType;
113     typedef typename CoefficientImageType::SpacingType      SpacingType;
114     typedef typename CoefficientImageType::DirectionType    DirectionType;
115     typedef typename CoefficientImageType::PointType        OriginType;
116     typedef itk::ContinuousIndex<TCoordRep, InputDimension> ContinuousIndexType;
117   
118     //JV added for the BLUT interpolator
119     typedef itk::Vector<TCoordRep, InputDimension> OutputSpacingType;
120
121     //JV m_VectorInterpolator
122     typedef VectorBSplineResampleImageFunctionWithLUT
123     <CoefficientImageType, TCoordRep> VectorInterpolatorType;
124     typedef  typename  VectorInterpolatorType::CoefficientDataType CoefficientDataType;
125     typedef  typename  VectorInterpolatorType::CoefficientDataType WeightsDataType;
126
127     /** Typedef of the bulk transform. */
128     typedef itk::Transform<ScalarType, InputDimension, OutputDimension> BulkTransformType;
129     typedef BulkTransformType*                  BulkTransformPointer;
130
131     /** Typedef of the mask */
132     typedef itk::SpatialObject<  InputDimension >   MaskType;
133     typedef MaskType*    MaskPointer;
134
135     //====================================================================
136     // Set et Gets
137     //====================================================================
138     //JV  added for the BLUT interpolator
139     void SetSplineOrder(const unsigned int &  splineOrder);
140     void SetSplineOrders(const SizeType &  splineOrders);
141     itkGetMacro( SplineOrders, SizeType );
142     itkGetConstMacro( SplineOrders, SizeType );
143     void SetLUTSamplingFactor(const int &  samplingFactor);
144     void SetLUTSamplingFactors(const SizeType &  samplingFactors);
145     itkGetMacro( LUTSamplingFactors, SizeType );
146     itkGetConstMacro( LUTSamplingFactors,SizeType );
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     // The padded coefficient image
170     virtual CoefficientImagePointer GetPaddedCoefficientImage()
171     { return m_PaddedCoefficientImage; }
172     virtual const CoefficientImagePointer  GetPaddedCoefficientImage() const
173     { return m_PaddedCoefficientImage; }
174     // virtual void SetPaddedCoefficientImage(CoefficientImagePointer image);  
175
176     /** This method specifies the region over which the grid resides. */
177     virtual void SetGridRegion( const RegionType& region );
178     itkGetMacro( GridRegion, RegionType );
179     itkGetConstMacro( GridRegion, RegionType );
180
181     /** This method specifies the grid spacing or resolution. */
182     virtual void SetGridSpacing( const SpacingType& spacing );
183     itkGetMacro( GridSpacing, SpacingType );
184     itkGetConstMacro( GridSpacing, SpacingType );
185
186     /** This method specifies the grid directions . */
187     virtual void SetGridDirection( const DirectionType & spacing );
188     itkGetMacro( GridDirection, DirectionType );
189     itkGetConstMacro( GridDirection, DirectionType );
190
191     /** This method specifies the grid origin. */
192     virtual void SetGridOrigin( const OriginType& origin );
193     itkGetMacro( GridOrigin, OriginType );
194     itkGetConstMacro( GridOrigin, OriginType );
195     
196     // Set the bulk transform, real pointer
197     // itkSetConstObjectMacro( BulkTransform, BulkTransformType );
198     // itkGetConstObjectMacro( BulkTransform, BulkTransformType );
199     void SetBulkTransform(BulkTransformPointer b){m_BulkTransform=b;}
200     BulkTransformPointer GetBulkTransform(void) {return m_BulkTransform;}
201
202     //Set mask, inside transform applies, outside zero, real pointer
203     void SetMask(MaskPointer m){m_Mask=m;}
204
205     // JV the shape
206     itkSetMacro( TransformShape , unsigned int  );
207     itkGetMacro( TransformShape , unsigned int  );
208     itkGetConstMacro( TransformShape, unsigned int  );
209
210     /** Transform points by a BSpline deformable transformation. */
211     OutputPointType  TransformPoint(const InputPointType  &point ) const;
212   
213     // JV added for just the deformable part, without bulk
214     OutputPointType  DeformablyTransformPoint(const InputPointType  &point ) const;
215   
216     /** Parameter index array type. */
217     typedef itk::Array<unsigned long> ParameterIndexArrayType;
218
219     /** Transform points by a BSpline deformable transformation. 
220      * On return, weights contains the interpolation weights used to compute the 
221      * deformation and indices of the x (zeroth) dimension coefficient parameters
222      * in the support region used to compute the deformation.
223      * Parameter indices for the i-th dimension can be obtained by adding
224      * ( i * this->GetNumberOfParametersPerDimension() ) to the indices array.
225      */
226
227     // JV not implemented
228     //  virtual void TransformPoint( const InputPointType & inputPoint,
229     //                           OutputPointType & outputPoint,
230     //                           WeightsType & weights,
231     //                           ParameterIndexArrayType & indices, 
232     //                           bool & inside ) const;
233     //     virtual void DeformablyTransformPoint( const InputPointType & inputPoint,
234     //                                     OutputPointType & outputPoint,
235     //                                     WeightsType & weights,
236     //                                     ParameterIndexArrayType & indices, 
237     //                                     bool & inside ) const;
238     //     virtual void GetJacobian( const InputPointType & inputPoint,
239     //                                WeightsType & weights,
240     //                                ParameterIndexArrayType & indices
241     //                                ) const;
242    
243     /** Method to transform a vector - 
244      *  not applicable for this type of transform. */
245     virtual OutputVectorType TransformVector(const InputVectorType &) const
246     { 
247       itkExceptionMacro(<< "Method not applicable for deformable transform." );
248       return OutputVectorType(); 
249     }
250
251     /** Method to transform a vnl_vector - 
252      *  not applicable for this type of transform */
253     virtual OutputVnlVectorType TransformVector(const InputVnlVectorType &) const
254     { 
255       itkExceptionMacro(<< "Method not applicable for deformable transform. ");
256       return OutputVnlVectorType(); 
257     }
258
259     /** Method to transform a CovariantVector - 
260      *  not applicable for this type of transform */
261     virtual OutputCovariantVectorType TransformCovariantVector(
262                                                                const InputCovariantVectorType &) const
263     { 
264       itkExceptionMacro(<< "Method not applicable for deformable transfrom. ");
265       return OutputCovariantVectorType(); 
266     } 
267     
268     /** Compute the Jacobian Matrix of the transformation at one point */
269 #if ITK_VERSION_MAJOR >= 4
270     virtual void ComputeJacobianWithRespectToParameters (const InputPointType &p, JacobianType &jacobian) const;
271     virtual void ComputeJacobianWithRespectToPosition (const InputPointType &p, JacobianType &jacobian) const
272     {
273       itkExceptionMacro( "ComputeJacobianWithRespectToPosition not yet implemented for " << this->GetNameOfClass() );
274     }
275 #else
276     virtual const JacobianType& GetJacobian(const InputPointType  &point ) const;
277 #endif
278
279     /** Return the number of parameters that completely define the Transfom */
280 #if ITK_VERSION_MAJOR >= 4
281     virtual NumberOfParametersType GetNumberOfParameters(void) const;
282 #else
283     virtual unsigned int GetNumberOfParameters(void) const;
284 #endif
285
286     //JV Return the padded number of parameters
287     virtual unsigned int GetPaddedNumberOfParameters(void) const;
288
289     /** Return the number of parameters per dimension */
290     unsigned int GetNumberOfParametersPerDimension(void) const;
291
292     /** Return the region of the grid wholly within the support region */
293     itkGetConstReferenceMacro( ValidRegion, RegionType );
294
295     /** Indicates that this transform is linear. That is, given two
296      * points P and Q, and scalar coefficients a and b, then
297      *
298      *           T( a*P + b*Q ) = a * T(P) + b * T(Q)
299      */
300     virtual bool IsLinear() const { return false; }
301
302    //unsigned int GetNumberOfAffectedWeights() const;
303
304   protected:
305     /** Print contents of an BSplineSpatioTemporalDeformableTransform. */
306     void PrintSelf(std::ostream &os, itk::Indent indent) const;
307
308
309     ShapedBLUTSpatioTemporalDeformableTransform();
310     virtual ~ShapedBLUTSpatioTemporalDeformableTransform();
311
312     /** Wrap flat array into images of coefficients. */
313     void WrapAsImages();
314
315     // JV Pad/Extract the coefficient image
316     void PadCoefficientImage(void);
317     typename  CoefficientImageType::Pointer ExtractTemporalRow(const typename CoefficientImageType::Pointer& coefficientImage, unsigned int temporalIndex);
318     //void ExtractCoefficientImage(void);
319
320     /** Wrap flat array into images of coefficients. */
321     inline void WrapRegion(const RegionType& support, 
322                            RegionType& first, 
323                            RegionType& second, 
324                            RegionType& third, 
325                            std::vector<RegionType>& bc, 
326                            std::vector<double>& bcValues,
327                            std::vector<RegionType>& bc2, 
328                            std::vector<double>& bc2Values,
329                            std::vector<RegionType>& bc3, 
330                            std::vector<double>& bc3Values,
331                            unsigned int& m_InitialOffset ) const;
332
333     /** Convert an input point to a continuous index inside the BSpline grid */
334     void TransformPointToContinuousIndex( const InputPointType & point, ContinuousIndexType & index ) const;
335
336   private:
337     ShapedBLUTSpatioTemporalDeformableTransform(const Self&); //purposely not implemented
338     void operator=(const Self&); //purposely not implemented
339
340     /** The bulk transform. */
341     BulkTransformPointer  m_BulkTransform;
342     MaskPointer m_Mask;
343
344     // JV added for BLUT interpolator
345     SizeType m_SplineOrders;
346     //OutputSpacingType m_OutputSpacing;
347     SizeType m_LUTSamplingFactors;
348
349     /** Variables defining the coefficient grid extend. */
350     RegionType    m_GridRegion;
351     SpacingType   m_GridSpacing;
352     DirectionType m_GridDirection;
353     OriginType    m_GridOrigin;
354
355     // JV additional variables for the padded region
356     RegionType    m_PaddedGridRegion;
357     OriginType    m_PaddedGridOrigin;    
358   
359
360     DirectionType m_PointToIndex;
361     DirectionType m_IndexToPoint;
362   
363     RegionType    m_ValidRegion;
364
365     /** Variables defining the interpolation support region. */
366     SizeType m_Offset;
367     itk::FixedArray<bool,InputDimension>    m_SplineOrderOdd;
368     SizeType      m_SupportSize;
369     IndexType     m_ValidRegionLast;
370   
371     /** Array holding images wrapped from the flat parameters. */
372     CoefficientImagePointer   m_WrappedImage;
373
374     /** Vector image representing the B-spline coefficients 
375      *  in each dimension. */
376     CoefficientImagePointer   m_CoefficientImage;
377     CoefficientImagePointer   m_PaddedCoefficientImage;
378
379     /** Jacobian as OutputDimension number of images. */
380     typedef typename JacobianType::ValueType JacobianValueType;
381     typedef typename itk::Vector<JacobianValueType,SpaceDimension> JacobianPixelType;
382     typedef itk::Image<JacobianPixelType, OutputDimension> JacobianImageType;
383     typename JacobianImageType::Pointer m_JacobianImage[OutputDimension];
384     typename JacobianImageType::Pointer m_PaddedJacobianImage[OutputDimension];
385     typedef itk::ImageRegionIterator<JacobianImageType> IteratorType;
386
387     //JV for J calculation
388     IndexType m_NullIndex;
389     SizeType m_NullSize;
390
391     mutable RegionType    m_SupportRegion;
392     mutable IndexType    m_SupportIndex;
393     mutable RegionType    m_FirstRegion;
394     mutable RegionType    m_SecondRegion;
395     mutable RegionType    m_ThirdRegion;
396     mutable unsigned int m_ThirdSize;
397     mutable unsigned int m_InitialOffset;
398
399     mutable std::vector<IteratorType>  m_BCIterators[SpaceDimension];
400     mutable std::vector<double>   m_BCValues;
401     mutable std::vector<RegionType>  m_BCRegions;
402     mutable unsigned int m_BCSize;
403
404     mutable std::vector<IteratorType>  m_BC2Iterators[SpaceDimension];
405     mutable std::vector<double>   m_BC2Values;
406     mutable std::vector<RegionType>  m_BC2Regions;
407     mutable unsigned int m_BC2Size;
408
409     mutable std::vector<IteratorType>  m_BC3Iterators[SpaceDimension];
410     mutable std::vector<double>   m_BC3Values;
411     mutable std::vector<RegionType>  m_BC3Regions;
412     mutable unsigned int m_BC3Size;
413
414     mutable RegionType m_BCRegion;
415     mutable IteratorType m_FirstIterator[SpaceDimension];
416     mutable IteratorType m_SecondIterator[SpaceDimension];
417     mutable IteratorType m_ThirdIterator[SpaceDimension];
418     mutable ContinuousIndexType m_Index;
419
420     //JV add a padded Jacobian matrix
421     mutable JacobianType m_PaddedJacobian;
422     mutable JacobianPixelType m_ZeroVector;
423
424     /** Keep track of last support region used in computing the Jacobian
425      * for fast resetting of Jacobian to zero.
426      */
427     mutable IndexType m_LastJacobianIndex;
428    
429     /** Keep a pointer to the input parameters. */
430     const ParametersType *  m_InputParametersPointer;
431
432     /** Internal parameters buffer. */
433     ParametersType          m_InternalParametersBuffer;
434
435     //JV  the BLUT interpolator
436     typename VectorInterpolatorType::Pointer m_VectorInterpolator;
437
438     // the coefficients to apply the BC
439     std::vector<std::vector<double> > m_Weights;
440     std::vector<std::vector<double> > m_WeightRatio;
441   
442     /** Check if a continuous index is inside the valid region. */
443     bool InsideValidRegion( const ContinuousIndexType& index ) const;
444
445     // JV Shape
446     unsigned int m_TransformShape;
447
448 #if ITK_VERSION_MAJOR >= 4
449     mutable JacobianType                            m_SharedDataBSplineJacobian;
450 #endif
451
452   }; //class ShapedBLUTSpatioTemporalDeformableTransform
453
454
455 }  // namespace itk
456
457 #if ITK_TEMPLATE_TXX
458 # include "clitkShapedBLUTSpatioTemporalDeformableTransform.txx"
459 #endif
460
461
462 #endif // __clitkShapedBLUTSpatioTemporalDeformableTransform_h