1 /*=========================================================================
2 Program: vv http://www.creatis.insa-lyon.fr/rio/vv
5 - University of LYON http://www.universite-lyon.fr/
6 - Léon Bérard cancer center http://oncora1.lyon.fnclcc.fr
7 - CREATIS CNRS laboratory http://www.creatis.insa-lyon.fr
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.
13 It is distributed under dual licence
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"
25 #include "itkTransform.h"
27 #include "itkImageRegion.h"
28 #include "itkSpatialObject.h"
29 #include "itkPasteImageFilter.h"
30 #include "itkMultiplyByConstantImageFilter.h"
36 class TCoordRep = double, // Data type for scalars, coordinate representation,vectors
37 unsigned int NInputDimensions = 3, // Number of input dimensions
38 unsigned int NOutputDimensions = 3 > // Number of output dimensions
39 class ITK_EXPORT ShapedBLUTSpatioTemporalDeformableTransform :
40 public itk::Transform< TCoordRep, NInputDimensions, NOutputDimensions >
44 //====================================================================
46 //====================================================================
47 typedef ShapedBLUTSpatioTemporalDeformableTransform Self;
48 typedef itk::Transform< TCoordRep, NInputDimensions, NOutputDimensions > Superclass;
49 typedef itk::SmartPointer<Self> Pointer;
50 typedef itk::SmartPointer<const Self> ConstPointer;
52 /** New macro for creation of through the object factory.*/
55 /** Run-time type information (and related methods). */
56 itkTypeMacro( ShapedBLUTSpatioTemporalDeformableTransform, Transform );
58 /** Dimension of the output space. */
59 itkStaticConstMacro(OutputDimension, unsigned int, NOutputDimensions);
61 /** Dimension of the input space. */
62 itkStaticConstMacro(InputDimension, unsigned int, NInputDimensions);
64 //JV the number of spatial dimensions
65 itkStaticConstMacro(SpaceDimension, unsigned int, NInputDimensions-1);
67 /** Standard scalar type for this class. */
68 typedef typename Superclass::ScalarType ScalarType;
70 /** Standard parameters container. */
71 typedef typename Superclass::ParametersType ParametersType;
73 /** Standard Jacobian container. */
74 typedef typename Superclass::JacobianType JacobianType;
76 /** Standard vector type for this class. */
77 typedef itk::Vector<TCoordRep,
78 InputDimension> InputVectorType;
79 typedef itk::Vector<TCoordRep,
80 OutputDimension> OutputVectorType;
82 /** Standard covariant vector type for this class. */
83 typedef itk::CovariantVector<TCoordRep,
84 InputDimension> InputCovariantVectorType;
85 typedef itk::CovariantVector<TCoordRep,
86 OutputDimension> OutputCovariantVectorType;
88 /** Standard vnl_vector type for this class. */
89 typedef vnl_vector_fixed<TCoordRep,
90 InputDimension> InputVnlVectorType;
91 typedef vnl_vector_fixed<TCoordRep,
92 OutputDimension> OutputVnlVectorType;
94 /** Standard coordinate point type for this class. */
95 typedef itk::Point<TCoordRep,
96 InputDimension> InputPointType;
97 typedef itk::Point<TCoordRep,
98 OutputDimension> OutputPointType;
100 //JV Parameters as images with OutputDimension number of components per Pixel
101 typedef typename ParametersType::ValueType ParametersValueType;
102 typedef typename itk::Vector<ParametersValueType, SpaceDimension> PixelType;
103 typedef itk::Image<PixelType, InputDimension> CoefficientImageType;
104 typedef typename CoefficientImageType::Pointer CoefficientImagePointer;
107 /** Typedefs for specifying the extend to the grid. */
108 typedef itk::ImageRegion<InputDimension> RegionType;
109 typedef typename RegionType::IndexType IndexType;
110 typedef typename RegionType::SizeType SizeType;
111 typedef typename CoefficientImageType::SpacingType SpacingType;
112 typedef typename CoefficientImageType::DirectionType DirectionType;
113 typedef typename CoefficientImageType::PointType OriginType;
114 typedef itk::ContinuousIndex<TCoordRep, InputDimension> ContinuousIndexType;
116 //JV added for the BLUT interpolator
117 typedef itk::Vector<TCoordRep, InputDimension> OutputSpacingType;
119 //JV m_VectorInterpolator
120 typedef VectorBSplineResampleImageFunctionWithLUT
121 <CoefficientImageType, TCoordRep> VectorInterpolatorType;
122 typedef typename VectorInterpolatorType::CoefficientDataType CoefficientDataType;
123 typedef typename VectorInterpolatorType::CoefficientDataType WeightsDataType;
125 /** Typedef of the bulk transform. */
126 typedef itk::Transform<ScalarType, InputDimension, OutputDimension> BulkTransformType;
127 typedef BulkTransformType* BulkTransformPointer;
129 /** Typedef of the mask */
130 typedef itk::SpatialObject< InputDimension > MaskType;
131 typedef MaskType* MaskPointer;
133 //====================================================================
135 //====================================================================
136 //JV added for the BLUT interpolator
137 void SetSplineOrder(const unsigned int & splineOrder);
138 void SetSplineOrders(const SizeType & splineOrders);
139 itkGetMacro( SplineOrders, SizeType );
140 itkGetConstMacro( SplineOrders, SizeType );
141 void SetLUTSamplingFactor(const int & samplingFactor);
142 void SetLUTSamplingFactors(const SizeType & samplingFactors);
143 itkGetMacro( LUTSamplingFactors, SizeType );
144 itkGetConstMacro( LUTSamplingFactors,SizeType );
146 void SetParameters(const ParametersType & parameters);
148 void SetFixedParameters(const ParametersType & parameters);
150 void SetParametersByValue(const ParametersType & parameters);
154 /** Get the Transformation Parameters. */
155 virtual const ParametersType& GetParameters(void) const;
157 /** Get the Transformation Fixed Parameters. */
158 virtual const ParametersType& GetFixedParameters(void) const;
160 // The coefficientImage
161 virtual CoefficientImagePointer GetCoefficientImage()
162 { return m_CoefficientImage; }
163 virtual const CoefficientImagePointer GetCoefficientImage() const
164 { return m_CoefficientImage; }
165 virtual void SetCoefficientImage(CoefficientImagePointer image);
167 // The padded coefficient image
168 virtual CoefficientImagePointer GetPaddedCoefficientImage()
169 { return m_PaddedCoefficientImage; }
170 virtual const CoefficientImagePointer GetPaddedCoefficientImage() const
171 { return m_PaddedCoefficientImage; }
172 // virtual void SetPaddedCoefficientImage(CoefficientImagePointer image);
174 /** This method specifies the region over which the grid resides. */
175 virtual void SetGridRegion( const RegionType& region );
176 itkGetMacro( GridRegion, RegionType );
177 itkGetConstMacro( GridRegion, RegionType );
179 /** This method specifies the grid spacing or resolution. */
180 virtual void SetGridSpacing( const SpacingType& spacing );
181 itkGetMacro( GridSpacing, SpacingType );
182 itkGetConstMacro( GridSpacing, SpacingType );
184 /** This method specifies the grid directions . */
185 virtual void SetGridDirection( const DirectionType & spacing );
186 itkGetMacro( GridDirection, DirectionType );
187 itkGetConstMacro( GridDirection, DirectionType );
189 /** This method specifies the grid origin. */
190 virtual void SetGridOrigin( const OriginType& origin );
191 itkGetMacro( GridOrigin, OriginType );
192 itkGetConstMacro( GridOrigin, OriginType );
194 // Set the bulk transform, real pointer
195 // itkSetConstObjectMacro( BulkTransform, BulkTransformType );
196 // itkGetConstObjectMacro( BulkTransform, BulkTransformType );
197 void SetBulkTransform(BulkTransformPointer b){m_BulkTransform=b;}
198 BulkTransformPointer GetBulkTransform(void) {return m_BulkTransform;}
200 //Set mask, inside transform applies, outside zero, real pointer
201 void SetMask(MaskPointer m){m_Mask=m;}
204 itkSetMacro( TransformShape , unsigned int );
205 itkGetMacro( TransformShape , unsigned int );
206 itkGetConstMacro( TransformShape, unsigned int );
208 /** Transform points by a BSpline deformable transformation. */
209 OutputPointType TransformPoint(const InputPointType &point ) const;
211 // JV added for just the deformable part, without bulk
212 OutputPointType DeformablyTransformPoint(const InputPointType &point ) const;
214 /** Parameter index array type. */
215 typedef itk::Array<unsigned long> ParameterIndexArrayType;
217 /** Transform points by a BSpline deformable transformation.
218 * On return, weights contains the interpolation weights used to compute the
219 * deformation and indices of the x (zeroth) dimension coefficient parameters
220 * in the support region used to compute the deformation.
221 * Parameter indices for the i-th dimension can be obtained by adding
222 * ( i * this->GetNumberOfParametersPerDimension() ) to the indices array.
225 // JV not implemented
226 // virtual void TransformPoint( const InputPointType & inputPoint,
227 // OutputPointType & outputPoint,
228 // WeightsType & weights,
229 // ParameterIndexArrayType & indices,
230 // bool & inside ) const;
231 // virtual void DeformablyTransformPoint( const InputPointType & inputPoint,
232 // OutputPointType & outputPoint,
233 // WeightsType & weights,
234 // ParameterIndexArrayType & indices,
235 // bool & inside ) const;
236 // virtual void GetJacobian( const InputPointType & inputPoint,
237 // WeightsType & weights,
238 // ParameterIndexArrayType & indices
241 /** Method to transform a vector -
242 * not applicable for this type of transform. */
243 virtual OutputVectorType TransformVector(const InputVectorType &) const
245 itkExceptionMacro(<< "Method not applicable for deformable transform." );
246 return OutputVectorType();
249 /** Method to transform a vnl_vector -
250 * not applicable for this type of transform */
251 virtual OutputVnlVectorType TransformVector(const InputVnlVectorType &) const
253 itkExceptionMacro(<< "Method not applicable for deformable transform. ");
254 return OutputVnlVectorType();
257 /** Method to transform a CovariantVector -
258 * not applicable for this type of transform */
259 virtual OutputCovariantVectorType TransformCovariantVector(
260 const InputCovariantVectorType &) const
262 itkExceptionMacro(<< "Method not applicable for deformable transfrom. ");
263 return OutputCovariantVectorType();
266 /** Compute the Jacobian Matrix of the transformation at one point */
267 virtual const JacobianType& GetJacobian(const InputPointType &point ) const;
269 /** Return the number of parameters that completely define the Transfom */
270 virtual unsigned int GetNumberOfParameters(void) const;
272 //JV Return the padded number of parameters
273 virtual unsigned int GetPaddedNumberOfParameters(void) const;
275 /** Return the number of parameters per dimension */
276 unsigned int GetNumberOfParametersPerDimension(void) const;
278 /** Return the region of the grid wholly within the support region */
279 itkGetConstReferenceMacro( ValidRegion, RegionType );
281 /** Indicates that this transform is linear. That is, given two
282 * points P and Q, and scalar coefficients a and b, then
284 * T( a*P + b*Q ) = a * T(P) + b * T(Q)
286 virtual bool IsLinear() const { return false; }
288 //unsigned int GetNumberOfAffectedWeights() const;
291 /** Print contents of an BSplineSpatioTemporalDeformableTransform. */
292 void PrintSelf(std::ostream &os, itk::Indent indent) const;
295 ShapedBLUTSpatioTemporalDeformableTransform();
296 virtual ~ShapedBLUTSpatioTemporalDeformableTransform();
298 /** Wrap flat array into images of coefficients. */
301 // JV Pad/Extract the coefficient image
302 void PadCoefficientImage(void);
303 typename CoefficientImageType::Pointer ExtractTemporalRow(const typename CoefficientImageType::Pointer& coefficientImage, unsigned int temporalIndex);
304 //void ExtractCoefficientImage(void);
306 /** Wrap flat array into images of coefficients. */
307 inline void WrapRegion(const RegionType& support,
311 std::vector<RegionType>& bc,
312 std::vector<double>& bcValues,
313 std::vector<RegionType>& bc2,
314 std::vector<double>& bc2Values,
315 std::vector<RegionType>& bc3,
316 std::vector<double>& bc3Values,
317 unsigned int& m_InitialOffset ) const;
319 /** Convert an input point to a continuous index inside the BSpline grid */
320 void TransformPointToContinuousIndex( const InputPointType & point, ContinuousIndexType & index ) const;
323 ShapedBLUTSpatioTemporalDeformableTransform(const Self&); //purposely not implemented
324 void operator=(const Self&); //purposely not implemented
326 /** The bulk transform. */
327 BulkTransformPointer m_BulkTransform;
330 // JV added for BLUT interpolator
331 SizeType m_SplineOrders;
332 //OutputSpacingType m_OutputSpacing;
333 SizeType m_LUTSamplingFactors;
335 /** Variables defining the coefficient grid extend. */
336 RegionType m_GridRegion;
337 SpacingType m_GridSpacing;
338 DirectionType m_GridDirection;
339 OriginType m_GridOrigin;
341 // JV additional variables for the padded region
342 RegionType m_PaddedGridRegion;
343 OriginType m_PaddedGridOrigin;
346 DirectionType m_PointToIndex;
347 DirectionType m_IndexToPoint;
349 RegionType m_ValidRegion;
351 /** Variables defining the interpolation support region. */
353 itk::FixedArray<bool,InputDimension> m_SplineOrderOdd;
354 SizeType m_SupportSize;
355 IndexType m_ValidRegionLast;
357 /** Array holding images wrapped from the flat parameters. */
358 CoefficientImagePointer m_WrappedImage;
360 /** Vector image representing the B-spline coefficients
361 * in each dimension. */
362 CoefficientImagePointer m_CoefficientImage;
363 CoefficientImagePointer m_PaddedCoefficientImage;
365 /** Jacobian as OutputDimension number of images. */
366 typedef typename JacobianType::ValueType JacobianValueType;
367 typedef typename itk::Vector<JacobianValueType,SpaceDimension> JacobianPixelType;
368 typedef itk::Image<JacobianPixelType, OutputDimension> JacobianImageType;
369 typename JacobianImageType::Pointer m_JacobianImage[OutputDimension];
370 typename JacobianImageType::Pointer m_PaddedJacobianImage[OutputDimension];
371 typedef itk::ImageRegionIterator<JacobianImageType> IteratorType;
373 //JV for J calculation
374 IndexType m_NullIndex;
377 mutable RegionType m_SupportRegion;
378 mutable IndexType m_SupportIndex;
379 mutable RegionType m_FirstRegion;
380 mutable RegionType m_SecondRegion;
381 mutable RegionType m_ThirdRegion;
382 mutable unsigned int m_ThirdSize;
383 mutable unsigned int m_InitialOffset;
385 mutable std::vector<IteratorType> m_BCIterators[SpaceDimension];
386 mutable std::vector<double> m_BCValues;
387 mutable std::vector<RegionType> m_BCRegions;
388 mutable unsigned int m_BCSize;
390 mutable std::vector<IteratorType> m_BC2Iterators[SpaceDimension];
391 mutable std::vector<double> m_BC2Values;
392 mutable std::vector<RegionType> m_BC2Regions;
393 mutable unsigned int m_BC2Size;
395 mutable std::vector<IteratorType> m_BC3Iterators[SpaceDimension];
396 mutable std::vector<double> m_BC3Values;
397 mutable std::vector<RegionType> m_BC3Regions;
398 mutable unsigned int m_BC3Size;
400 mutable RegionType m_BCRegion;
401 mutable IteratorType m_FirstIterator[SpaceDimension];
402 mutable IteratorType m_SecondIterator[SpaceDimension];
403 mutable IteratorType m_ThirdIterator[SpaceDimension];
404 mutable ContinuousIndexType m_Index;
406 //JV add a padded Jacobian matrix
407 mutable JacobianType m_PaddedJacobian;
408 mutable JacobianPixelType m_ZeroVector;
410 /** Keep track of last support region used in computing the Jacobian
411 * for fast resetting of Jacobian to zero.
413 mutable IndexType m_LastJacobianIndex;
415 /** Keep a pointer to the input parameters. */
416 const ParametersType * m_InputParametersPointer;
418 /** Internal parameters buffer. */
419 ParametersType m_InternalParametersBuffer;
421 //JV the BLUT interpolator
422 typename VectorInterpolatorType::Pointer m_VectorInterpolator;
424 // the coefficients to apply the BC
425 std::vector<std::vector<double> > m_Weights;
426 std::vector<std::vector<double> > m_WeightRatio;
428 /** Check if a continuous index is inside the valid region. */
429 bool InsideValidRegion( const ContinuousIndexType& index ) const;
432 unsigned int m_TransformShape;
435 }; //class ShapedBLUTSpatioTemporalDeformableTransform
441 # include "clitkShapedBLUTSpatioTemporalDeformableTransform.txx"
445 #endif // __clitkShapedBLUTSpatioTemporalDeformableTransform_h