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://www.centreleonberard.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"
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 >
43 //====================================================================
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;
51 /** New macro for creation of through the object factory.*/
54 /** Run-time type information (and related methods). */
55 itkTypeMacro( ShapedBLUTSpatioTemporalDeformableTransform, Transform );
57 /** Dimension of the output space. */
58 itkStaticConstMacro(OutputDimension, unsigned int, NOutputDimensions);
60 /** Dimension of the input space. */
61 itkStaticConstMacro(InputDimension, unsigned int, NInputDimensions);
63 //JV the number of spatial dimensions
64 itkStaticConstMacro(SpaceDimension, unsigned int, NInputDimensions-1);
66 /** Standard scalar type for this class. */
67 typedef typename Superclass::ScalarType ScalarType;
69 /** Standard parameters container. */
70 typedef typename Superclass::ParametersType ParametersType;
71 typedef typename Superclass::NumberOfParametersType NumberOfParametersType;
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 void ComputeJacobianWithRespectToParameters (const InputPointType &p, JacobianType &jacobian) const;
268 virtual void ComputeJacobianWithRespectToPosition (const InputPointType &p, JacobianType &jacobian) const
270 itkExceptionMacro( "ComputeJacobianWithRespectToPosition not yet implemented for " << this->GetNameOfClass() );
273 /** Return the number of parameters that completely define the Transfom */
274 virtual NumberOfParametersType GetNumberOfParameters(void) const;
276 //JV Return the padded number of parameters
277 virtual unsigned int GetPaddedNumberOfParameters(void) const;
279 /** Return the number of parameters per dimension */
280 unsigned int GetNumberOfParametersPerDimension(void) const;
282 /** Return the region of the grid wholly within the support region */
283 itkGetConstReferenceMacro( ValidRegion, RegionType );
285 /** Indicates that this transform is linear. That is, given two
286 * points P and Q, and scalar coefficients a and b, then
288 * T( a*P + b*Q ) = a * T(P) + b * T(Q)
290 virtual bool IsLinear() const { return false; }
292 //unsigned int GetNumberOfAffectedWeights() const;
295 /** Print contents of an BSplineSpatioTemporalDeformableTransform. */
296 void PrintSelf(std::ostream &os, itk::Indent indent) const;
299 ShapedBLUTSpatioTemporalDeformableTransform();
300 virtual ~ShapedBLUTSpatioTemporalDeformableTransform();
302 /** Wrap flat array into images of coefficients. */
305 // JV Pad/Extract the coefficient image
306 void PadCoefficientImage(void);
307 typename CoefficientImageType::Pointer ExtractTemporalRow(const typename CoefficientImageType::Pointer& coefficientImage, unsigned int temporalIndex);
308 //void ExtractCoefficientImage(void);
310 /** Wrap flat array into images of coefficients. */
311 inline void WrapRegion(const RegionType& support,
315 std::vector<RegionType>& bc,
316 std::vector<double>& bcValues,
317 std::vector<RegionType>& bc2,
318 std::vector<double>& bc2Values,
319 std::vector<RegionType>& bc3,
320 std::vector<double>& bc3Values,
321 unsigned int& m_InitialOffset ) const;
323 /** Convert an input point to a continuous index inside the BSpline grid */
324 void TransformPointToContinuousIndex( const InputPointType & point, ContinuousIndexType & index ) const;
327 ShapedBLUTSpatioTemporalDeformableTransform(const Self&); //purposely not implemented
328 void operator=(const Self&); //purposely not implemented
330 /** The bulk transform. */
331 BulkTransformPointer m_BulkTransform;
334 // JV added for BLUT interpolator
335 SizeType m_SplineOrders;
336 //OutputSpacingType m_OutputSpacing;
337 SizeType m_LUTSamplingFactors;
339 /** Variables defining the coefficient grid extend. */
340 RegionType m_GridRegion;
341 SpacingType m_GridSpacing;
342 DirectionType m_GridDirection;
343 OriginType m_GridOrigin;
345 // JV additional variables for the padded region
346 RegionType m_PaddedGridRegion;
347 OriginType m_PaddedGridOrigin;
350 DirectionType m_PointToIndex;
351 DirectionType m_IndexToPoint;
353 RegionType m_ValidRegion;
355 /** Variables defining the interpolation support region. */
357 itk::FixedArray<bool,InputDimension> m_SplineOrderOdd;
358 SizeType m_SupportSize;
359 IndexType m_ValidRegionLast;
361 /** Array holding images wrapped from the flat parameters. */
362 CoefficientImagePointer m_WrappedImage;
364 /** Vector image representing the B-spline coefficients
365 * in each dimension. */
366 CoefficientImagePointer m_CoefficientImage;
367 CoefficientImagePointer m_PaddedCoefficientImage;
369 /** Jacobian as OutputDimension number of images. */
370 typedef typename JacobianType::ValueType JacobianValueType;
371 typedef typename itk::Vector<JacobianValueType,SpaceDimension> JacobianPixelType;
372 typedef itk::Image<JacobianPixelType, OutputDimension> JacobianImageType;
373 typename JacobianImageType::Pointer m_JacobianImage[OutputDimension];
374 typename JacobianImageType::Pointer m_PaddedJacobianImage[OutputDimension];
375 typedef itk::ImageRegionIterator<JacobianImageType> IteratorType;
377 //JV for J calculation
378 IndexType m_NullIndex;
381 mutable RegionType m_SupportRegion;
382 mutable IndexType m_SupportIndex;
383 mutable RegionType m_FirstRegion;
384 mutable RegionType m_SecondRegion;
385 mutable RegionType m_ThirdRegion;
386 mutable unsigned int m_ThirdSize;
387 mutable unsigned int m_InitialOffset;
389 mutable std::vector<IteratorType> m_BCIterators[SpaceDimension];
390 mutable std::vector<double> m_BCValues;
391 mutable std::vector<RegionType> m_BCRegions;
392 mutable unsigned int m_BCSize;
394 mutable std::vector<IteratorType> m_BC2Iterators[SpaceDimension];
395 mutable std::vector<double> m_BC2Values;
396 mutable std::vector<RegionType> m_BC2Regions;
397 mutable unsigned int m_BC2Size;
399 mutable std::vector<IteratorType> m_BC3Iterators[SpaceDimension];
400 mutable std::vector<double> m_BC3Values;
401 mutable std::vector<RegionType> m_BC3Regions;
402 mutable unsigned int m_BC3Size;
404 mutable RegionType m_BCRegion;
405 mutable IteratorType m_FirstIterator[SpaceDimension];
406 mutable IteratorType m_SecondIterator[SpaceDimension];
407 mutable IteratorType m_ThirdIterator[SpaceDimension];
408 mutable ContinuousIndexType m_Index;
410 //JV add a padded Jacobian matrix
411 mutable JacobianType m_PaddedJacobian;
412 mutable JacobianPixelType m_ZeroVector;
414 /** Keep track of last support region used in computing the Jacobian
415 * for fast resetting of Jacobian to zero.
417 mutable IndexType m_LastJacobianIndex;
419 /** Keep a pointer to the input parameters. */
420 const ParametersType * m_InputParametersPointer;
422 /** Internal parameters buffer. */
423 ParametersType m_InternalParametersBuffer;
425 //JV the BLUT interpolator
426 typename VectorInterpolatorType::Pointer m_VectorInterpolator;
428 // the coefficients to apply the BC
429 std::vector<std::vector<double> > m_Weights;
430 std::vector<std::vector<double> > m_WeightRatio;
432 /** Check if a continuous index is inside the valid region. */
433 bool InsideValidRegion( const ContinuousIndexType& index ) const;
436 unsigned int m_TransformShape;
438 mutable JacobianType m_SharedDataBSplineJacobian;
440 }; //class ShapedBLUTSpatioTemporalDeformableTransform
445 #ifndef ITK_MANUAL_INSTANTIATION
446 # include "clitkShapedBLUTSpatioTemporalDeformableTransform.txx"
450 #endif // __clitkShapedBLUTSpatioTemporalDeformableTransform_h