ITK  4.2.0
Insight Segmentation and Registration Toolkit
itkRegionBasedLevelSetFunction.h
Go to the documentation of this file.
1 /*=========================================================================
2  *
3  * Copyright Insight Software Consortium
4  *
5  * Licensed under the Apache License, Version 2.0 (the "License");
6  * you may not use this file except in compliance with the License.
7  * You may obtain a copy of the License at
8  *
9  * http://www.apache.org/licenses/LICENSE-2.0.txt
10  *
11  * Unless required by applicable law or agreed to in writing, software
12  * distributed under the License is distributed on an "AS IS" BASIS,
13  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14  * See the License for the specific language governing permissions and
15  * limitations under the License.
16  *
17  *=========================================================================*/
18 #ifndef __itkRegionBasedLevelSetFunction_h
19 #define __itkRegionBasedLevelSetFunction_h
20 
23 #include "vnl/vnl_matrix_fixed.h"
24 
25 namespace itk
26 {
64 template< class TInput, // LevelSetImageType
65  class TFeature, // FeatureImageType
66  class TSharedData >
67 class ITK_EXPORT RegionBasedLevelSetFunction:public
68  FiniteDifferenceFunction< TInput >
69 {
70 public:
76 
77  itkStaticConstMacro(ImageDimension, unsigned int, Superclass::ImageDimension);
78 
79  // itkNewMacro() is purposely not provided since this is an abstract class.
80 
83 
85  typedef double TimeStepType;
86  typedef typename Superclass::ImageType ImageType;
87  typedef typename Superclass::PixelType PixelType;
89  typedef typename Superclass::RadiusType RadiusType;
90  typedef typename Superclass::NeighborhoodType NeighborhoodType;
91  typedef typename Superclass::NeighborhoodScalesType NeighborhoodScalesType;
92  typedef typename Superclass::FloatOffsetType FloatOffsetType;
95 
96  /* This structure is derived from LevelSetFunction and stores intermediate
97  values for computing time step sizes */
100  {
102 
103  m_MaxCurvatureChange = null_value;
104  m_MaxAdvectionChange = null_value;
105  m_MaxGlobalChange = null_value;
106  }
107 
109 
110  vnl_matrix_fixed< ScalarValueType,
111  itkGetStaticConstMacro(ImageDimension),
112  itkGetStaticConstMacro(ImageDimension) > m_dxy;
113 
114  ScalarValueType m_dx[itkGetStaticConstMacro(ImageDimension)];
115 
116  ScalarValueType m_dx_forward[itkGetStaticConstMacro(ImageDimension)];
117  ScalarValueType m_dx_backward[itkGetStaticConstMacro(ImageDimension)];
118 
121 
125  };
126 
127  typedef TInput InputImageType;
128  typedef typename InputImageType::ConstPointer InputImageConstPointer;
129  typedef typename InputImageType::Pointer InputImagePointer;
130  typedef typename InputImageType::PixelType InputPixelType;
131  typedef typename InputImageType::IndexType InputIndexType;
133  typedef typename InputImageType::SizeType InputSizeType;
135  typedef typename InputImageType::RegionType InputRegionType;
136  typedef typename InputImageType::PointType InputPointType;
137 
138  typedef TFeature FeatureImageType;
139  typedef typename FeatureImageType::ConstPointer FeatureImageConstPointer;
140  typedef typename FeatureImageType::PixelType FeaturePixelType;
141  typedef typename FeatureImageType::IndexType FeatureIndexType;
142  typedef typename FeatureImageType::SpacingType FeatureSpacingType;
143  typedef typename FeatureImageType::OffsetType FeatureOffsetType;
144 
145  typedef TSharedData SharedDataType;
146  typedef typename SharedDataType::Pointer SharedDataPointer;
147 
150 
151  void SetDomainFunction(const HeavisideFunctionType *f)
152  {
153  this->m_DomainFunction = f;
154  }
155 
156  virtual void Initialize(const RadiusType & r)
157  {
158  this->SetRadius(r);
159 
160  // Dummy neighborhood.
161  NeighborhoodType it;
162  it.SetRadius(r);
163 
164  // Find the center index of the neighborhood.
165  m_Center = it.Size() / 2;
166 
167  // Get the stride length for each axis.
168  for ( unsigned int i = 0; i < ImageDimension; i++ )
169  {
170  m_xStride[i] = it.GetStride(i);
171  }
172  }
173 
174 #if !defined( CABLE_CONFIGURATION )
175  void SetSharedData(SharedDataPointer sharedDataIn)
176  {
177  this->m_SharedData = sharedDataIn;
178  }
179 #endif
180 
181  void UpdateSharedData(bool forceUpdate);
182 
183  void * GetGlobalDataPointer() const
184  {
185  return new GlobalDataStruct;
186  }
187 
188  TimeStepType ComputeGlobalTimeStep(void *GlobalData) const;
189 
191  virtual PixelType ComputeUpdate( const NeighborhoodType & neighborhood,
192  void *globalData, const FloatOffsetType & = FloatOffsetType(0.0) );
193 
194  void SetInitialImage(InputImageType *f)
195  {
196  m_InitialImage = f;
197  }
198 
199  virtual const FeatureImageType * GetFeatureImage() const
200  { return m_FeatureImage.GetPointer(); }
201  virtual void SetFeatureImage(const FeatureImageType *f)
202  {
203  m_FeatureImage = f;
204 
205  FeatureSpacingType spacing = m_FeatureImage->GetSpacing();
206  for ( unsigned int i = 0; i < ImageDimension; i++ )
207  {
208  this->m_InvSpacing[i] = 1 / spacing[i];
209  }
210  }
211 
213  virtual VectorType AdvectionField(const NeighborhoodType &,
214  const FloatOffsetType &, GlobalDataStruct * = 0) const
215  { return this->m_ZeroVectorConstant; }
216 
218  void SetAreaWeight(const ScalarValueType & nu)
219  { this->m_AreaWeight = nu; }
220  ScalarValueType GetAreaWeight() const
221  { return this->m_AreaWeight; }
223 
225  void SetLambda1(const ScalarValueType & lambda1)
226  { this->m_Lambda1 = lambda1; }
227  ScalarValueType GetLambda1() const
228  { return this->m_Lambda1; }
230 
232  void SetLambda2(const ScalarValueType & lambda2)
233  { this->m_Lambda2 = lambda2; }
234  ScalarValueType GetLambda2() const
235  { return this->m_Lambda2; }
237 
239  void SetOverlapPenaltyWeight(const ScalarValueType & gamma)
240  { this->m_OverlapPenaltyWeight = gamma; }
241  ScalarValueType GetOverlapPenaltyWeight() const
242  { return this->m_OverlapPenaltyWeight; }
244 
246  virtual void SetCurvatureWeight(const ScalarValueType c)
247  { m_CurvatureWeight = c; }
248  ScalarValueType GetCurvatureWeight() const
249  { return m_CurvatureWeight; }
251 
252  void SetAdvectionWeight(const ScalarValueType & iA)
253  { this->m_AdvectionWeight = iA; }
254  ScalarValueType GetAdvectionWeight() const
255  { return this->m_AdvectionWeight; }
256 
258  void SetReinitializationSmoothingWeight(const ScalarValueType c)
259  { m_ReinitializationSmoothingWeight = c; }
260  ScalarValueType GetReinitializationSmoothingWeight() const
261  { return m_ReinitializationSmoothingWeight; }
263 
265  void SetVolumeMatchingWeight(const ScalarValueType & tau)
266  { this->m_VolumeMatchingWeight = tau; }
267  ScalarValueType GetVolumeMatchingWeight() const
268  { return this->m_VolumeMatchingWeight; }
270 
272  void SetVolume(const ScalarValueType & volume)
273  { this->m_Volume = volume; }
274  ScalarValueType GetVolume() const
275  { return this->m_Volume; }
277 
279  void SetFunctionId(const unsigned int & iFid)
280  { this->m_FunctionId = iFid; }
281 
282  virtual void ReleaseGlobalDataPointer(void *GlobalData) const
283  { delete (GlobalDataStruct *)GlobalData; }
284 
285  virtual ScalarValueType ComputeCurvature(const NeighborhoodType &,
286  const FloatOffsetType &, GlobalDataStruct *gd);
287 
290  virtual ScalarValueType LaplacianSmoothingSpeed(
291  const NeighborhoodType &,
292  const FloatOffsetType &, GlobalDataStruct * = 0) const
294 
297  virtual ScalarValueType CurvatureSpeed(const NeighborhoodType &,
298  const FloatOffsetType &, GlobalDataStruct * = 0
299  ) const
301 
306  virtual void CalculateAdvectionImage() {}
307 protected:
308 
311 
314 
317 
319 
321 
324 
327 
330 
333 
336 
339 
342 
344 
347 
348  unsigned int m_FunctionId;
349 
350  std::slice x_slice[itkGetStaticConstMacro(ImageDimension)];
352  OffsetValueType m_xStride[itkGetStaticConstMacro(ImageDimension)];
353  double m_InvSpacing[itkGetStaticConstMacro(ImageDimension)];
354 
355  static double m_WaveDT;
356  static double m_DT;
357 
358  void ComputeHImage();
359 
362  ScalarValueType ComputeGlobalTerm(
363  const ScalarValueType & imagePixel,
364  const InputIndexType & inputIndex);
365 
370  virtual ScalarValueType ComputeInternalTerm(const FeaturePixelType & iValue,
371  const FeatureIndexType & iIdx) = 0;
372 
376  virtual ScalarValueType ComputeExternalTerm(const FeaturePixelType & iValue,
377  const FeatureIndexType & iIdx) = 0;
378 
383  virtual ScalarValueType ComputeOverlapParameters(const FeatureIndexType & featIndex,
384  ScalarValueType & pr) = 0;
385 
390  ScalarValueType ComputeVolumeRegularizationTerm();
391 
405  ScalarValueType ComputeLaplacian(GlobalDataStruct *gd);
406 
408  void ComputeHessian(const NeighborhoodType & it,
409  GlobalDataStruct *globalData);
410 
412  virtual void ComputeParameters() = 0;
413 
416  virtual void UpdateSharedDataParameters() = 0;
417 
418  bool m_UpdateC;
419 
422  static VectorType InitializeZeroVectorConstant();
423 
426 private:
427  RegionBasedLevelSetFunction(const Self &); //purposely not implemented
428  void operator=(const Self &); //purposely not implemented
429 };
430 } // end namespace itk
432 
433 #ifndef ITK_MANUAL_INSTANTIATION
434 #include "itkRegionBasedLevelSetFunction.hxx"
435 #endif
436 
437 #endif
438