KWStyle - itkKernelTransform.h

Matrix View

Description

---

1		/*=========================================================================
2
3		Program:   Insight Segmentation & Registration Toolkit
4		Module:    $RCSfile: itkKernelTransform.h.html,v $
5		Language:  C++
6		Date:      $Date: 2006/01/17 19:15:40 $
7		Version:   $Revision: 1.4 $
8
9		Copyright (c) Insight Software Consortium. All rights reserved.
10		See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details.
11
12		This software is distributed WITHOUT ANY WARRANTY; without even
13		the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
14		PURPOSE.  See the above copyright notices for more information.
15
16		=========================================================================*/
17		#ifndef __itkKernelTransform_h
18		#define __itkKernelTransform_h
19
20		#include "itkTransform.h"
21		#include "itkPoint.h"
22		#include "itkVector.h"
23		#include "itkMatrix.h"
24		#include "itkPointSet.h"
25		#include <deque>
26		#include <math.h>
27		#include "vnl/vnl_matrix_fixed.h"
28		#include "vnl/vnl_matrix.h"
29		#include "vnl/vnl_vector.h"
30		#include "vnl/vnl_vector_fixed.h"
31		#include "vnl/algo/vnl_svd.h"
32		#include "vnl/vnl_sample.h"
33
34		namespace itk
35		{
36
37		/** \class KernelTransform
38		* Intended to be a base class for elastic body spline and thin plate spline.
39		* This is implemented in as straightforward a manner as possible from the
40		* IEEE TMI paper by Davis, Khotanzad, Flamig, and Harms, Vol. 16,
41		* No. 3 June 1997. Notation closely follows their paper, so if you have it
42		* in front of you, this code will make a lot more sense.
43		*
44		* KernelTransform:
45		*  Provides support for defining source and target landmarks
46		*  Defines a number of data types used in the computations
47		*  Defines the mathematical framework used to compute all splines,
48		*    so that subclasses need only provide a kernel specific to
49		*    that spline
50		*
51		* This formulation allows the stiffness of the spline to
52		* be adjusted, allowing the spline to vary from interpolating the
53		* landmarks to approximating the landmarks.  This part of the
54		* formulation is based on the short paper by R. Sprengel, K. Rohr,
55		* H. Stiehl. "Thin-Plate Spline Approximation for Image
56		* Registration". In 18th International Conference of the IEEE
57		* Engineering in Medicine and Biology Society. 1996.
58		*
59		* \ingroup Transforms
60		*
61		*/
62		template <class TScalarType, // probably only float and double make sense here
63		unsigned int NDimensions>   // Number of dimensions
64		class ITK_EXPORT KernelTransform :
65	IND	********************public Transform<TScalarType, NDimensions,NDimensions>
66		{
67		public:
68		/** Standard class typedefs. */
69		typedef KernelTransform Self;
70	TDA	typedef Transform<TScalarType, NDimensions, NDimensions >   Superclass;
71	TDA	typedef SmartPointer<Self>        Pointer;
72	TDA	typedef SmartPointer<const Self>  ConstPointer;
73
74		/** Run-time type information (and related methods). */
75		itkTypeMacro( KernelTransform, Transform );
76
77		/** New macro for creation of through a Smart Pointer */
78		itkNewMacro( Self );
79
80		/** Dimension of the domain space. */
81		itkStaticConstMacro(SpaceDimension, unsigned int, NDimensions);
82
83		/** Scalar type. */
84		typedef typename Superclass::ScalarType  ScalarType;
85
86		/** Parameters type. */
87		typedef typename Superclass::ParametersType  ParametersType;
88
89		/** Jacobian type. */
90		typedef typename Superclass::JacobianType  JacobianType;
91
92		/** Standard coordinate point type for this class. */
93		typedef typename Superclass::InputPointType   InputPointType;
94		typedef typename Superclass::OutputPointType  OutputPointType;
95
96		/** Standard vector type for this class. */
97		typedef typename Superclass::InputVectorType   InputVectorType;
98		typedef typename Superclass::OutputVectorType  OutputVectorType;
99
100		/** PointList typedef. This type is used for maintaining lists of points,
101		* specifically, the source and target landmark lists. */
102		typedef DefaultStaticMeshTraits<TScalarType,
103		NDimensions,
104		NDimensions,
105		TScalarType,
106		TScalarType> PointSetTraitsType;
107	LEN,TDA	typedef PointSet<InputPointType, NDimensions, PointSetTraitsType> PointSetType;
108	TDA	typedef typename PointSetType::Pointer                        PointSetPointer;
109	TDA	typedef typename PointSetType::PointsContainer                PointsContainer;
110	TDA	typedef typename PointSetType::PointsContainerIterator        PointsIterator;
111	LEN,TDA	typedef typename PointSetType::PointsContainerConstIterator   PointsConstIterator;
112
113		/** VectorSet typedef. */
114		typedef itk::VectorContainer<unsigned long,InputVectorType> VectorSetType;
115	TDA	typedef typename VectorSetType::Pointer        VectorSetPointer;
116
117		/** Get the source landmarks list, which we will denote \f$ p \f$. */
118		itkGetObjectMacro( SourceLandmarks, PointSetType);
119
120		/** Set the source landmarks list. */
121		virtual void SetSourceLandmarks(PointSetType *);
122
123		/** Get the target landmarks list, which we will denote  \f$ q \f$. */
124		itkGetObjectMacro( TargetLandmarks, PointSetType);
125
126		/** Set the target landmarks list. */
127		virtual void SetTargetLandmarks(PointSetType *);
128
129		/** Get the displacements list, which we will denote \f$ d \f$,
130		* where \f$ d_i = q_i - p_i \f$. */
131		itkGetObjectMacro( Displacements, VectorSetType );
132
133		/** Compute W matrix. */
134		void ComputeWMatrix(void);
135
136		/** Compute the position of point in the new space */
137		virtual OutputPointType TransformPoint(const InputPointType& thisPoint) const;
138
139		/** 'I' (identity) matrix typedef. */
140		typedef vnl_matrix_fixed<TScalarType, NDimensions, NDimensions> IMatrixType;
141
142
143		/** Compute the Jacobian Matrix of the transformation at one point */
144		virtual const JacobianType & GetJacobian(const InputPointType  &point ) const;
145
146		/** Set the Transformation Parameters and update the internal transformation.
147	LEN	* The parameters represent the source landmarks. Each landmark point is represented
148		* by NDimensions doubles. All the landmarks are concatenated to form one flat
149		* Array<double>. */
150		virtual void SetParameters(const ParametersType &);
151
152		/** Set Transform Fixed Parameters:
153		*     To support the transform file writer this function was
154		*     added to set the target landmarks similar to the
155		*     SetParameters function setting the source landmarks
156		*/
157		virtual void SetFixedParameters(const ParametersType &);
158
159		/** Update the Parameters array from the landmarks corrdinates. */
160		virtual void UpdateParameters(void) const;
161
162		/** Get the Transformation Parameters - Gets the Source Landmarks */
163		virtual const ParametersType& GetParameters(void) const;
164
165		/** Get Transform Fixed Parameters - Gets the Target Landmarks */
166		virtual const ParametersType& GetFixedParameters(void) const;
167
168		/** Stiffness of the spline.  A stiffness of zero results in the
169		* standard interpolating spline.  A non-zero stiffness allows the
170		* spline to approximate rather than interpolate the landmarks.
171		* Stiffness values are usually rather small, typically in the range
172		* of 0.001 to 0.1. The approximating spline formulation is based on
173		* the short paper by R. Sprengel, K. Rohr, H. Stiehl. "Thin-Plate
174		* Spline Approximation for Image Registration". In 18th
175		* International Conference of the IEEE Engineering in Medicine and
176		* Biology Society. 1996.
177		*/
178		itkSetClampMacro(Stiffness, double, 0.0, NumericTraits<double>::max());
179		itkGetMacro(Stiffness, double);
180
181
182		protected:
183		KernelTransform();
184		virtual ~KernelTransform();
185		void PrintSelf(std::ostream& os, Indent indent) const;
186
187		public:
188		/** 'G' matrix typedef. */
189		typedef vnl_matrix_fixed<TScalarType, NDimensions, NDimensions> GMatrixType;
190
191		/** 'L' matrix typedef. */
192		typedef vnl_matrix<TScalarType> LMatrixType;
193
194		/** 'K' matrix typedef. */
195		typedef vnl_matrix<TScalarType> KMatrixType;
196
197		/** 'P' matrix typedef. */
198		typedef vnl_matrix<TScalarType> PMatrixType;
199
200		/** 'Y' matrix typedef. */
201		typedef vnl_matrix<TScalarType> YMatrixType;
202
203		/** 'W' matrix typedef. */
204		typedef vnl_matrix<TScalarType> WMatrixType;
205
206		/** 'D' matrix typedef. Deformation component */
207		typedef vnl_matrix<TScalarType> DMatrixType;
208
209		/** 'A' matrix typedef. Rotational part of the Affine component */
210		typedef vnl_matrix_fixed<TScalarType,NDimensions,NDimensions> AMatrixType;
211
212		/** 'B' matrix typedef. Translational part of the Affine component */
213		typedef vnl_vector_fixed<TScalarType,NDimensions> BMatrixType;
214
215		/** Row matrix typedef. */
216		typedef vnl_matrix_fixed<TScalarType, 1, NDimensions> RowMatrixType;
217
218		/** Column matrix typedef. */
219		typedef vnl_matrix_fixed<TScalarType, NDimensions, 1> ColumnMatrixType;
220
221		/** The list of source landmarks, denoted 'p'. */
222		PointSetPointer m_SourceLandmarks;
223
224		/** The list of target landmarks, denoted 'q'. */
225		PointSetPointer m_TargetLandmarks;
226
227		protected:
228		/** Compute G(x)
229		* This is essentially the kernel of the transform.
230		* By overriding this method, we can obtain (among others):
231		*    Elastic body spline
232		*    Thin plate spline
233		*    Volume spline */
234	LEN	virtual const GMatrixType & ComputeG(const InputVectorType & landmarkVector) const;
235
236		/** Compute a G(x) for a point to itself (i.e. for the block
237		* diagonal elements of the matrix K. Parameter indicates for which
238		* landmark the reflexive G is to be computed. The default
239		* implementation for the reflexive contribution is a diagonal
240		* matrix where the diagonal elements are the stiffness of the
241		* spline. */
242		virtual const GMatrixType & ComputeReflexiveG(PointsIterator) const;
243
244
245		/** Compute the contribution of the landmarks weighted by the kernel funcion
246	IND	******to the global deformation of the space  */
247	LEN	virtual void ComputeDeformationContribution( const InputPointType & inputPoint,
248	LEN	OutputPointType & result ) const;
249
250		/** Compute K matrix. */
251		void ComputeK();
252
253		/** Compute L matrix. */
254		void ComputeL();
255
256		/** Compute P matrix. */
257		void ComputeP();
258
259		/** Compute Y matrix. */
260		void ComputeY();
261
262		/** Compute displacements \f$ q_i - p_i \f$. */
263		void ComputeD();
264
265		/** Reorganize the components of W into
266	IND	****D (deformable), A (rotation part of affine)
267	IND	****and B (translational part of affine ) components.
268	IND	****\warning This method release the memory of the W Matrix  */
269		void ReorganizeW(void);
270
271		/** Stiffness parameter */
272		double m_Stiffness;
273
274		/** The list of displacements.
275		* d[i] = q[i] - p[i]; */
276		VectorSetPointer m_Displacements;
277
278		/** The L matrix. */
279		LMatrixType m_LMatrix;
280
281		/** The K matrix. */
282		KMatrixType m_KMatrix;
283
284		/** The P matrix. */
285		PMatrixType m_PMatrix;
286
287		/** The Y matrix. */
288		YMatrixType m_YMatrix;
289
290		/** The W matrix. */
291		WMatrixType m_WMatrix;
292
293		/** The Deformation matrix.
294	IND	******This is an auxiliary matrix that will hold the
295	IND	******Deformation (non-affine) part of the transform.
296	IND	******Those are the coefficients that will multiply the
297	IND	******Kernel function */
298		DMatrixType m_DMatrix;
299
300		/** Rotatinoal/Shearing part of the Affine component of the Transformation */
301		AMatrixType m_AMatrix;
302
303		/** Translational part of the Affine component of the Transformation */
304		BMatrixType m_BVector;
305
306		/** The G matrix.
307	IND	****It is made mutable because m_GMatrix was made an ivar
308	IND	****only to avoid copying the matrix at return time */
309		mutable GMatrixType m_GMatrix;
310
311		/** Has the W matrix been computed? */
312		bool m_WMatrixComputed;
313
314		/** Identity matrix. */
315		IMatrixType m_I;
316
317	IND	*private:
318		KernelTransform(const Self&); //purposely not implemented
319		void operator=(const Self&); //purposely not implemented
320
321		};
322
323		} // end namespace itk
324
325		#ifndef ITK_MANUAL_INSTANTIATION
326		#include "itkKernelTransform.txx"
327		#endif
328
329		#endif // __itkKernelTransform_h
330

---

Generated by KWStyle 1.0b on Tuesday January,17 at 02:14:17PM

© Kitware Inc.