KWStyle - itkLevelSetFunction.h

Matrix View

Description

---

1		/*=========================================================================
2
3		Program:   Insight Segmentation & Registration Toolkit
4		Module:    $RCSfile: itkLevelSetFunction.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	DEF	#ifndef __itkLevelSetFunction_h_
18	DEF	#define __itkLevelSetFunction_h_
19
20		#include "itkFiniteDifferenceFunction.h"
21		#include "vnl/vnl_matrix_fixed.h"
22
23		namespace itk {
24
25		/** \class LevelSetFunction
26	IND	*** \brief The LevelSetFunction class is a generic function object which can be
27		* used to create a level set method filter when combined with an appropriate
28		* finite difference image filter.  (See FiniteDifferenceImageFilter.)
29		*
30		* LevelSetFunction implements a generic level set function.  This function is
31		* an expanded form of the basic equation developed in [1].
32		*
33		* \f$\phi_{t} + \alpha
34		* \stackrel{\rightharpoonup}{A}(\mathbf{x})\cdot\nabla\phi + \beta
35	LEN	* P(\mathbf{x})\mid\nabla\phi\mid = \gamma Z(\mathbf{x})\kappa\mid\nabla\phi\mid\f$
36		*
37		* where \f$ \stackrel{\rightharpoonup}{A} \f$ is an advection term, \f$ P \f$
38		* is a propagation (growth) term, and \f$ Z \f$ is a spatial modifier term for
39		* the mean curvature \f$ \kappa \f$.  \f$ \alpha \f$, \f$ \beta \f$, and \f$
40		* \gamma \f$ are all scalar constants.
41		*
42		* Terms in the equation above are supplied through virtual methods, which must
43		* be subclassed to complete an implementation.  Terms can be eliminated from
44		* the equation by setting the corresponding constants to zero. A wide variety
45		* of level set methods can be implemented by subclassing this basic equation.
46		*
47		* In ITK, the usual sign convention is that the INSIDE of a surface contains
48		* NEGATIVE values and the OUTSIDE of the surface contains POSITIVE values.
49		*
50		* \warning You MUST call Initialize() in the constructor of subclasses of this
51		* object to set it up properly to do level-set Calculations.  The argument
52		* that you pass Initialize is the radius of the neighborhood needed to perform
53		* the calculations.  If your subclass does not do any additional neighborhood
54		* processing, then the default radius should be 1 in each direction.
55		*
56		* \par REFERENCES
57		* \par
58		* [1] Sethian, J.A. Level Set Methods. Cambridge University Press. 1996.
59		*
60		* \ingroup FiniteDifferenceFunctions
61		* \ingroup Functions
62		*/
63		template <class TImageType>
64		class ITK_EXPORT LevelSetFunction
65	IND	**: public FiniteDifferenceFunction<TImageType>
66		{
67		public:
68		/** Standard class typedefs. */
69		typedef LevelSetFunction Self;
70	TDA	typedef FiniteDifferenceFunction<TImageType> Superclass;
71	TDA	typedef SmartPointer<Self> Pointer;
72	TDA	typedef SmartPointer<const Self> ConstPointer;
73
74		/** Method for creation through the object factory. */
75		itkNewMacro(Self);
76
77		/** Run-time type information (and related methods) */
78		itkTypeMacro( LevelSetFunction, FiniteDifferenceFunction );
79
80		/** Extract some parameters from the superclass. */
81		itkStaticConstMacro(ImageDimension, unsigned int,Superclass::ImageDimension);
82
83		/** Convenient typedefs. */
84		typedef double TimeStepType;
85	TDA	typedef typename Superclass::ImageType  ImageType;
86	TDA	typedef typename Superclass::PixelType  PixelType;
87	TDA	typedef                      PixelType  ScalarValueType;
88	TDA	typedef typename Superclass::RadiusType RadiusType;
89	TDA	typedef typename Superclass::NeighborhoodType NeighborhoodType;
90	TDA	typedef typename Superclass::FloatOffsetType FloatOffsetType;
91
92		/** The vector type that will be used in the calculations. */
93		//  typedef
94	LEN	typedef FixedArray<ScalarValueType, itkGetStaticConstMacro(ImageDimension)> VectorType;
95
96		/** A global data type for this class of equations.  Used to store
97		* values that are needed in calculating the time step and other intermediate
98		* products such as derivatives that may be used by virtual functions called
99		* from ComputeUpdate.  Caching these values here allows the ComputeUpdate
100		* function to be const and thread safe.*/
101		struct GlobalDataStruct
102	IND	**{
103		ScalarValueType m_MaxAdvectionChange;
104		ScalarValueType m_MaxPropagationChange;
105		ScalarValueType m_MaxCurvatureChange;
106
107		/** Hessian matrix */
108		vnl_matrix_fixed<ScalarValueType,
109		itkGetStaticConstMacro(ImageDimension),
110		itkGetStaticConstMacro(ImageDimension)> m_dxy;
111
112		/** Array of first derivatives*/
113		ScalarValueType m_dx[itkGetStaticConstMacro(ImageDimension)];
114
115		ScalarValueType m_dx_forward[itkGetStaticConstMacro(ImageDimension)];
116		ScalarValueType m_dx_backward[itkGetStaticConstMacro(ImageDimension)];
117
118		ScalarValueType m_GradMagSqr;
119	IND	**};
120
121		/** Advection field.  Default implementation returns a vector of zeros. */
122		virtual VectorType AdvectionField(const NeighborhoodType &,
123	LEN	const FloatOffsetType &, GlobalDataStruct * = 0)  const
124		{ return m_ZeroVectorConstant; }
125
126		/** Propagation speed.  This term controls surface expansion/contraction.
127		*  Default implementation returns zero. */
128		virtual ScalarValueType PropagationSpeed(
129		const NeighborhoodType& ,
130		const FloatOffsetType &, GlobalDataStruct * = 0 ) const
131		{ return NumericTraits<ScalarValueType>::Zero; }
132
133		/** Curvature speed.  Can be used to spatially modify the effects of
134	IND	******curvature . The default implementation returns one. */
135		virtual ScalarValueType CurvatureSpeed(const NeighborhoodType &,
136	LEN	const FloatOffsetType &, GlobalDataStruct * = 0
137	IND	*****************************************) const
138		{ return NumericTraits<ScalarValueType>::One; }
139
140	IND	****/** Laplacian smoothing speed.  Can be used to spatially modify the
141	IND	******effects of laplacian smoothing of the level set function */
142		virtual ScalarValueType LaplacianSmoothingSpeed(
143		const NeighborhoodType &,
144		const FloatOffsetType &, GlobalDataStruct * = 0) const
145		{ return NumericTraits<ScalarValueType>::One; }
146
147		/** Alpha.  Scales all advection term values.*/
148		virtual void SetAdvectionWeight(const ScalarValueType a)
149		{ m_AdvectionWeight = a; }
150		ScalarValueType GetAdvectionWeight() const
151		{ return m_AdvectionWeight; }
152
153		/** Beta.  Scales all propagation term values. */
154		virtual void SetPropagationWeight(const ScalarValueType p)
155		{ m_PropagationWeight = p; }
156		ScalarValueType GetPropagationWeight() const
157		{ return m_PropagationWeight; }
158
159		/** Gamma. Scales all curvature weight values */
160		virtual void SetCurvatureWeight(const ScalarValueType c)
161		{ m_CurvatureWeight = c; }
162		ScalarValueType GetCurvatureWeight() const
163		{ return m_CurvatureWeight; }
164
165		/** Weight of the laplacian smoothing term */
166		void SetLaplacianSmoothingWeight(const ScalarValueType c)
167		{ m_LaplacianSmoothingWeight = c; }
168		ScalarValueType GetLaplacianSmoothingWeight() const
169		{ return m_LaplacianSmoothingWeight; }
170
171		/** Epsilon. */
172		void SetEpsilonMagnitude(const ScalarValueType e)
173		{ m_EpsilonMagnitude = e; }
174		ScalarValueType GetEpsilonMagnitude() const
175		{ return m_EpsilonMagnitude; }
176
177		/** Compute the equation value. */
178		virtual PixelType ComputeUpdate(const NeighborhoodType &neighborhood,
179		void *globalData,
180	LEN	const FloatOffsetType& = FloatOffsetType(0.0));
181
182	IND	*/** Computes the time step for an update given a global data structure.
183		* The data used in the computation may take different forms depending on
184		* the nature of the equations.  This global data cannot be kept in the
185		* instance of the equation object itself since the equation object must
186		* remain stateless for thread safety.  The global data is therefore managed
187		* for each thread by the finite difference solver filters. */
188		virtual TimeStepType ComputeGlobalTimeStep(void *GlobalData) const;
189
190		/** Returns a pointer to a global data structure that is passed to this
191		* object from the solver at each calculation.  The idea is that the solver
192		* holds the state of any global values needed to calculate the time step,
193		* while the equation object performs the actual calculations.  The global
194		* data should also be initialized in this method.  Global data can be used
195		* for caching any values used or reused by the FunctionObject.  Each thread
196		* should receive its own global data struct. */
197		virtual void *GetGlobalDataPointer() const
198		{
199	IND	******GlobalDataStruct *ans = new GlobalDataStruct();
200	IND	******ans->m_MaxAdvectionChange   = NumericTraits<ScalarValueType>::Zero;
201	IND	******ans->m_MaxPropagationChange = NumericTraits<ScalarValueType>::Zero;
202	IND	******ans->m_MaxCurvatureChange   = NumericTraits<ScalarValueType>::Zero;
203	IND	******return ans;
204		}
205
206		/** This method creates the appropriate member variable operators for the
207		* level-set calculations.  The argument to this function is a the radius
208		* necessary for performing the level-set calculations. */
209		virtual void Initialize(const RadiusType &r);
210
211		/** When the finite difference solver filter has finished using a global
212		* data pointer, it passes it to this method, which frees the memory.
213		* The solver cannot free the memory because it does not know the type
214		* to which the pointer points. */
215		virtual void ReleaseGlobalDataPointer(void *GlobalData) const
216		{ delete (GlobalDataStruct *) GlobalData; }
217
218		/**  */
219		virtual ScalarValueType ComputeCurvatureTerm(const NeighborhoodType &,
220		const FloatOffsetType &,
221		GlobalDataStruct *gd = 0
222	IND	***********************************************);
223		virtual ScalarValueType ComputeMeanCurvature(const NeighborhoodType &,
224		const FloatOffsetType &,
225		GlobalDataStruct *gd = 0
226	IND	***********************************************);
227
228		virtual ScalarValueType ComputeMinimalCurvature(const NeighborhoodType &,
229		const FloatOffsetType &,
230		GlobalDataStruct *gd = 0
231	IND	**************************************************);
232
233		virtual ScalarValueType Compute3DMinimalCurvature(const NeighborhoodType &,
234		const FloatOffsetType &,
235		GlobalDataStruct *gd = 0
236	IND	****************************************************);
237
238		/** */
239		void SetUseMinimalCurvature( bool b )
240	IND	**{
241		m_UseMinimalCurvature = b;
242	IND	**}
243		bool GetUseMinimalCurvature() const
244	IND	**{
245		return m_UseMinimalCurvature;
246	IND	**}
247		void UseMinimalCurvatureOn()
248	IND	**{
249		this->SetUseMinimalCurvature(true);
250	IND	**}
251		void UseMinimalCurvatureOff()
252	IND	**{
253		this->SetUseMinimalCurvature(false);
254	IND	**}
255
256	LEN	/** Set/Get the maximum constraint for the curvature term factor in the time step
257	IND	******calculation.  Changing this value from the default is not recommended or
258	IND	******necessary, but can be used to speed up the surface evolution at the risk
259	IND	******of creating an unstable solution.*/
260		static void SetMaximumCurvatureTimeStep(double n)
261	IND	**{
262		m_DT = n;
263	IND	**}
264		static double GetMaximumCurvatureTimeStep()
265	IND	**{
266		return m_DT;
267	IND	**}
268
269	LEN	/** Set/Get the maximum constraint for the scalar/vector term factor of the time step
270	IND	******calculation.  Changing this value from the default is not recommended or
271	IND	******necessary, but can be used to speed up the surface evolution at the risk
272	IND	******of creating an unstable solution.*/
273		static void SetMaximumPropagationTimeStep(double n)
274	IND	**{
275		m_WaveDT = n;
276	IND	**}
277		static double GetMaximumPropagationTimeStep()
278	IND	**{
279		return m_WaveDT;
280	IND	**}
281
282		protected:
283		LevelSetFunction()
284	IND	**{
285		m_EpsilonMagnitude = 1.0e-5;
286		m_AdvectionWeight = m_PropagationWeight
287	IND	******= m_CurvatureWeight = m_LaplacianSmoothingWeight
288	IND	******= NumericTraits<ScalarValueType>::Zero;
289		m_UseMinimalCurvature = false;
290	IND	**}
291		virtual ~LevelSetFunction() {}
292		void PrintSelf(std::ostream &s, Indent indent) const;
293
294		/** Constants used in the time step calculation. */
295		static double m_WaveDT;
296		static double m_DT;
297
298		/** Slices for the ND neighborhood. */
299		std::slice x_slice[itkGetStaticConstMacro(ImageDimension)];
300
301		/** The offset of the center pixel in the neighborhood. */
302		::size_t m_Center;
303
304		/** Stride length along the y-dimension. */
305		::size_t m_xStride[itkGetStaticConstMacro(ImageDimension)];
306
307		bool m_UseMinimalCurvature;
308
309		/** This method's only purpose is to initialize the zero vector
310		* constant. */
311		static VectorType InitializeZeroVectorConstant();
312
313		/** Zero vector constant. */
314		static VectorType m_ZeroVectorConstant;
315
316		/** Epsilon magnitude controls the lower limit for gradient magnitude. */
317		ScalarValueType m_EpsilonMagnitude;
318
319		/** Alpha. */
320		ScalarValueType m_AdvectionWeight;
321
322		/** Beta. */
323		ScalarValueType m_PropagationWeight;
324
325		/** Gamma. */
326		ScalarValueType m_CurvatureWeight;
327
328		/** Laplacean smoothing term */
329		ScalarValueType m_LaplacianSmoothingWeight;
330
331		private:
332		LevelSetFunction(const Self&); //purposely not implemented
333		void operator=(const Self&);   //purposely not implemented
334		};
335
336		} // namespace itk
337
338		#ifndef ITK_MANUAL_INSTANTIATION
339		#include "itkLevelSetFunction.txx"
340		#endif
341
342		#endif
343
344	EOF

---

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

© Kitware Inc.