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;
72 /** Standard Jacobian container. */
73 typedef typename Superclass::JacobianType JacobianType;
75 /** Standard vector type for this class. */
76 typedef itk::Vector<TCoordRep,
77 InputDimension> InputVectorType;
78 typedef itk::Vector<TCoordRep,
79 OutputDimension> OutputVectorType;
81 /** Standard covariant vector type for this class. */
82 typedef itk::CovariantVector<TCoordRep,
83 InputDimension> InputCovariantVectorType;
84 typedef itk::CovariantVector<TCoordRep,
85 OutputDimension> OutputCovariantVectorType;
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;
93 /** Standard coordinate point type for this class. */
94 typedef itk::Point<TCoordRep,
95 InputDimension> InputPointType;
96 typedef itk::Point<TCoordRep,
97 OutputDimension> OutputPointType;
99 //JV Parameters as images with OutputDimension number of components per Pixel
100 typedef typename ParametersType::ValueType ParametersValueType;
101 typedef typename itk::Vector<ParametersValueType, SpaceDimension> PixelType;
102 typedef itk::Image<PixelType, InputDimension> CoefficientImageType;
103 typedef typename CoefficientImageType::Pointer CoefficientImagePointer;
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;
115 //JV added for the BLUT interpolator
116 typedef itk::Vector<TCoordRep, InputDimension> OutputSpacingType;
118 //JV m_VectorInterpolator
119 typedef VectorBSplineResampleImageFunctionWithLUT
120 <CoefficientImageType, TCoordRep> VectorInterpolatorType;
121 typedef typename VectorInterpolatorType::CoefficientDataType CoefficientDataType;
122 typedef typename VectorInterpolatorType::CoefficientDataType WeightsDataType;
124 /** Typedef of the bulk transform. */
125 typedef itk::Transform<ScalarType, InputDimension, OutputDimension> BulkTransformType;
126 typedef BulkTransformType* BulkTransformPointer;
128 /** Typedef of the mask */
129 typedef itk::SpatialObject< InputDimension > MaskType;
130 typedef MaskType* MaskPointer;
132 //====================================================================
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 );
145 void SetParameters(const ParametersType & parameters);
147 void SetFixedParameters(const ParametersType & parameters);
149 void SetParametersByValue(const ParametersType & parameters);
153 /** Get the Transformation Parameters. */
154 virtual const ParametersType& GetParameters(void) const;
156 /** Get the Transformation Fixed Parameters. */
157 virtual const ParametersType& GetFixedParameters(void) const;
159 // The coefficientImage
160 virtual CoefficientImagePointer GetCoefficientImage()
161 { return m_CoefficientImage; }
162 virtual const CoefficientImagePointer GetCoefficientImage() const
163 { return m_CoefficientImage; }
164 virtual void SetCoefficientImage(CoefficientImagePointer image);
166 // The padded coefficient image
167 virtual CoefficientImagePointer GetPaddedCoefficientImage()
168 { return m_PaddedCoefficientImage; }
169 virtual const CoefficientImagePointer GetPaddedCoefficientImage() const
170 { return m_PaddedCoefficientImage; }
171 // virtual void SetPaddedCoefficientImage(CoefficientImagePointer image);
173 /** This method specifies the region over which the grid resides. */
174 virtual void SetGridRegion( const RegionType& region );
175 itkGetMacro( GridRegion, RegionType );
176 itkGetConstMacro( GridRegion, RegionType );
178 /** This method specifies the grid spacing or resolution. */
179 virtual void SetGridSpacing( const SpacingType& spacing );
180 itkGetMacro( GridSpacing, SpacingType );
181 itkGetConstMacro( GridSpacing, SpacingType );
183 /** This method specifies the grid directions . */
184 virtual void SetGridDirection( const DirectionType & spacing );
185 itkGetMacro( GridDirection, DirectionType );
186 itkGetConstMacro( GridDirection, DirectionType );
188 /** This method specifies the grid origin. */
189 virtual void SetGridOrigin( const OriginType& origin );
190 itkGetMacro( GridOrigin, OriginType );
191 itkGetConstMacro( GridOrigin, OriginType );
193 // Set the bulk transform, real pointer
194 // itkSetConstObjectMacro( BulkTransform, BulkTransformType );
195 // itkGetConstObjectMacro( BulkTransform, BulkTransformType );
196 void SetBulkTransform(BulkTransformPointer b){m_BulkTransform=b;}
197 BulkTransformPointer GetBulkTransform(void) {return m_BulkTransform;}
199 //Set mask, inside transform applies, outside zero, real pointer
200 void SetMask(MaskPointer m){m_Mask=m;}
203 itkSetMacro( TransformShape , unsigned int );
204 itkGetMacro( TransformShape , unsigned int );
205 itkGetConstMacro( TransformShape, unsigned int );
207 /** Transform points by a BSpline deformable transformation. */
208 OutputPointType TransformPoint(const InputPointType &point ) const;
210 // JV added for just the deformable part, without bulk
211 OutputPointType DeformablyTransformPoint(const InputPointType &point ) const;
213 /** Parameter index array type. */
214 typedef itk::Array<unsigned long> ParameterIndexArrayType;
216 /** Transform points by a BSpline deformable transformation.
217 * On return, weights contains the interpolation weights used to compute the
218 * deformation and indices of the x (zeroth) dimension coefficient parameters
219 * in the support region used to compute the deformation.
220 * Parameter indices for the i-th dimension can be obtained by adding
221 * ( i * this->GetNumberOfParametersPerDimension() ) to the indices array.
224 // JV not implemented
225 // virtual void TransformPoint( const InputPointType & inputPoint,
226 // OutputPointType & outputPoint,
227 // WeightsType & weights,
228 // ParameterIndexArrayType & indices,
229 // bool & inside ) const;
230 // virtual void DeformablyTransformPoint( const InputPointType & inputPoint,
231 // OutputPointType & outputPoint,
232 // WeightsType & weights,
233 // ParameterIndexArrayType & indices,
234 // bool & inside ) const;
235 // virtual void GetJacobian( const InputPointType & inputPoint,
236 // WeightsType & weights,
237 // ParameterIndexArrayType & indices
240 /** Method to transform a vector -
241 * not applicable for this type of transform. */
242 virtual OutputVectorType TransformVector(const InputVectorType &) const
244 itkExceptionMacro(<< "Method not applicable for deformable transform." );
245 return OutputVectorType();
248 /** Method to transform a vnl_vector -
249 * not applicable for this type of transform */
250 virtual OutputVnlVectorType TransformVector(const InputVnlVectorType &) const
252 itkExceptionMacro(<< "Method not applicable for deformable transform. ");
253 return OutputVnlVectorType();
256 /** Method to transform a CovariantVector -
257 * not applicable for this type of transform */
258 virtual OutputCovariantVectorType TransformCovariantVector(
259 const InputCovariantVectorType &) const
261 itkExceptionMacro(<< "Method not applicable for deformable transfrom. ");
262 return OutputCovariantVectorType();
265 /** Compute the Jacobian Matrix of the transformation at one point */
266 virtual const JacobianType& GetJacobian(const InputPointType &point ) const;
268 /** Return the number of parameters that completely define the Transfom */
269 virtual unsigned int GetNumberOfParameters(void) const;
271 //JV Return the padded number of parameters
272 virtual unsigned int GetPaddedNumberOfParameters(void) const;
274 /** Return the number of parameters per dimension */
275 unsigned int GetNumberOfParametersPerDimension(void) const;
277 /** Return the region of the grid wholly within the support region */
278 itkGetConstReferenceMacro( ValidRegion, RegionType );
280 /** Indicates that this transform is linear. That is, given two
281 * points P and Q, and scalar coefficients a and b, then
283 * T( a*P + b*Q ) = a * T(P) + b * T(Q)
285 virtual bool IsLinear() const { return false; }
287 //unsigned int GetNumberOfAffectedWeights() const;
290 /** Print contents of an BSplineSpatioTemporalDeformableTransform. */
291 void PrintSelf(std::ostream &os, itk::Indent indent) const;
294 ShapedBLUTSpatioTemporalDeformableTransform();
295 virtual ~ShapedBLUTSpatioTemporalDeformableTransform();
297 /** Wrap flat array into images of coefficients. */
300 // JV Pad/Extract the coefficient image
301 void PadCoefficientImage(void);
302 typename CoefficientImageType::Pointer ExtractTemporalRow(const typename CoefficientImageType::Pointer& coefficientImage, unsigned int temporalIndex);
303 //void ExtractCoefficientImage(void);
305 /** Wrap flat array into images of coefficients. */
306 inline void WrapRegion(const RegionType& support,
310 std::vector<RegionType>& bc,
311 std::vector<double>& bcValues,
312 std::vector<RegionType>& bc2,
313 std::vector<double>& bc2Values,
314 std::vector<RegionType>& bc3,
315 std::vector<double>& bc3Values,
316 unsigned int& m_InitialOffset ) const;
318 /** Convert an input point to a continuous index inside the BSpline grid */
319 void TransformPointToContinuousIndex( const InputPointType & point, ContinuousIndexType & index ) const;
322 ShapedBLUTSpatioTemporalDeformableTransform(const Self&); //purposely not implemented
323 void operator=(const Self&); //purposely not implemented
325 /** The bulk transform. */
326 BulkTransformPointer m_BulkTransform;
329 // JV added for BLUT interpolator
330 SizeType m_SplineOrders;
331 //OutputSpacingType m_OutputSpacing;
332 SizeType m_LUTSamplingFactors;
334 /** Variables defining the coefficient grid extend. */
335 RegionType m_GridRegion;
336 SpacingType m_GridSpacing;
337 DirectionType m_GridDirection;
338 OriginType m_GridOrigin;
340 // JV additional variables for the padded region
341 RegionType m_PaddedGridRegion;
342 OriginType m_PaddedGridOrigin;
345 DirectionType m_PointToIndex;
346 DirectionType m_IndexToPoint;
348 RegionType m_ValidRegion;
350 /** Variables defining the interpolation support region. */
352 itk::FixedArray<bool,InputDimension> m_SplineOrderOdd;
353 SizeType m_SupportSize;
354 IndexType m_ValidRegionLast;
356 /** Array holding images wrapped from the flat parameters. */
357 CoefficientImagePointer m_WrappedImage;
359 /** Vector image representing the B-spline coefficients
360 * in each dimension. */
361 CoefficientImagePointer m_CoefficientImage;
362 CoefficientImagePointer m_PaddedCoefficientImage;
364 /** Jacobian as OutputDimension number of images. */
365 typedef typename JacobianType::ValueType JacobianValueType;
366 typedef typename itk::Vector<JacobianValueType,SpaceDimension> JacobianPixelType;
367 typedef itk::Image<JacobianPixelType, OutputDimension> JacobianImageType;
368 typename JacobianImageType::Pointer m_JacobianImage[OutputDimension];
369 typename JacobianImageType::Pointer m_PaddedJacobianImage[OutputDimension];
370 typedef itk::ImageRegionIterator<JacobianImageType> IteratorType;
372 //JV for J calculation
373 IndexType m_NullIndex;
376 mutable RegionType m_SupportRegion;
377 mutable IndexType m_SupportIndex;
378 mutable RegionType m_FirstRegion;
379 mutable RegionType m_SecondRegion;
380 mutable RegionType m_ThirdRegion;
381 mutable unsigned int m_ThirdSize;
382 mutable unsigned int m_InitialOffset;
384 mutable std::vector<IteratorType> m_BCIterators[SpaceDimension];
385 mutable std::vector<double> m_BCValues;
386 mutable std::vector<RegionType> m_BCRegions;
387 mutable unsigned int m_BCSize;
389 mutable std::vector<IteratorType> m_BC2Iterators[SpaceDimension];
390 mutable std::vector<double> m_BC2Values;
391 mutable std::vector<RegionType> m_BC2Regions;
392 mutable unsigned int m_BC2Size;
394 mutable std::vector<IteratorType> m_BC3Iterators[SpaceDimension];
395 mutable std::vector<double> m_BC3Values;
396 mutable std::vector<RegionType> m_BC3Regions;
397 mutable unsigned int m_BC3Size;
399 mutable RegionType m_BCRegion;
400 mutable IteratorType m_FirstIterator[SpaceDimension];
401 mutable IteratorType m_SecondIterator[SpaceDimension];
402 mutable IteratorType m_ThirdIterator[SpaceDimension];
403 mutable ContinuousIndexType m_Index;
405 //JV add a padded Jacobian matrix
406 mutable JacobianType m_PaddedJacobian;
407 mutable JacobianPixelType m_ZeroVector;
409 /** Keep track of last support region used in computing the Jacobian
410 * for fast resetting of Jacobian to zero.
412 mutable IndexType m_LastJacobianIndex;
414 /** Keep a pointer to the input parameters. */
415 const ParametersType * m_InputParametersPointer;
417 /** Internal parameters buffer. */
418 ParametersType m_InternalParametersBuffer;
420 //JV the BLUT interpolator
421 typename VectorInterpolatorType::Pointer m_VectorInterpolator;
423 // the coefficients to apply the BC
424 std::vector<std::vector<double> > m_Weights;
425 std::vector<std::vector<double> > m_WeightRatio;
427 /** Check if a continuous index is inside the valid region. */
428 bool InsideValidRegion( const ContinuousIndexType& index ) const;
431 unsigned int m_TransformShape;
434 }; //class ShapedBLUTSpatioTemporalDeformableTransform
440 # include "clitkShapedBLUTSpatioTemporalDeformableTransform.txx"
444 #endif // __clitkShapedBLUTSpatioTemporalDeformableTransform_h