itkLevelSetFunction.h

Go to the documentation of this file.

/*=========================================================================

  Program:   Insight Segmentation & Registration Toolkit
  Module:    $RCSfile: itkLevelSetFunction.h,v $
  Language:  C++
  Date:      $Date: 2009-11-24 02:27:27 $
  Version:   $Revision: 1.26 $

  Copyright (c) Insight Software Consortium. All rights reserved.
  See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details.

     This software is distributed WITHOUT ANY WARRANTY; without even 
     the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR 
     PURPOSE.  See the above copyright notices for more information.

=========================================================================*/
#ifndef __itkLevelSetFunction_h
#define __itkLevelSetFunction_h

#include "itkFiniteDifferenceFunction.h"
#include "vnl/vnl_matrix_fixed.h"

namespace itk {

template <class TImageType>
class ITK_EXPORT LevelSetFunction
  : public FiniteDifferenceFunction<TImageType>
{
public:
  typedef LevelSetFunction                     Self;
  typedef FiniteDifferenceFunction<TImageType> Superclass;
  typedef SmartPointer<Self>                   Pointer;
  typedef SmartPointer<const Self>             ConstPointer;

  itkNewMacro(Self);

  itkTypeMacro( LevelSetFunction, FiniteDifferenceFunction );

  itkStaticConstMacro(ImageDimension, unsigned int,Superclass::ImageDimension);

  typedef double                                      TimeStepType;
  typedef typename Superclass::ImageType              ImageType;
  typedef typename Superclass::PixelType              PixelType;
  typedef                      PixelType              ScalarValueType;
  typedef typename Superclass::PixelRealType          PixelRealType;
  typedef typename Superclass::RadiusType             RadiusType;
  typedef typename Superclass::NeighborhoodType       NeighborhoodType;
  typedef typename Superclass::NeighborhoodScalesType NeighborhoodScalesType;
  typedef typename Superclass::FloatOffsetType        FloatOffsetType;

  //  typedef
  typedef FixedArray<ScalarValueType, itkGetStaticConstMacro(ImageDimension)> VectorType;

  struct GlobalDataStruct
    {
    ScalarValueType m_MaxAdvectionChange;
    ScalarValueType m_MaxPropagationChange;
    ScalarValueType m_MaxCurvatureChange;

    vnl_matrix_fixed<ScalarValueType,
                     itkGetStaticConstMacro(ImageDimension),
                     itkGetStaticConstMacro(ImageDimension)> m_dxy;
    
    ScalarValueType m_dx[itkGetStaticConstMacro(ImageDimension)];
    
    ScalarValueType m_dx_forward[itkGetStaticConstMacro(ImageDimension)];
    ScalarValueType m_dx_backward[itkGetStaticConstMacro(ImageDimension)];
    
    ScalarValueType m_GradMagSqr;
    };

  virtual VectorType AdvectionField(const NeighborhoodType &,
                                    const FloatOffsetType &, GlobalDataStruct * = 0)  const
    { return m_ZeroVectorConstant; }

  virtual ScalarValueType PropagationSpeed(
    const NeighborhoodType& ,
    const FloatOffsetType &, GlobalDataStruct * = 0 ) const
    { return NumericTraits<ScalarValueType>::Zero; }

  virtual ScalarValueType CurvatureSpeed(const NeighborhoodType &,
                                         const FloatOffsetType &, GlobalDataStruct * = 0
                                         ) const
    { return NumericTraits<ScalarValueType>::One; }

  virtual ScalarValueType LaplacianSmoothingSpeed(
    const NeighborhoodType &,
    const FloatOffsetType &, GlobalDataStruct * = 0) const
    { return NumericTraits<ScalarValueType>::One; }

  virtual void SetAdvectionWeight(const ScalarValueType a)
    { m_AdvectionWeight = a; }
  ScalarValueType GetAdvectionWeight() const
    { return m_AdvectionWeight; }
  
  virtual void SetPropagationWeight(const ScalarValueType p)
    { m_PropagationWeight = p; }
  ScalarValueType GetPropagationWeight() const
    { return m_PropagationWeight; }
  
  virtual void SetCurvatureWeight(const ScalarValueType c)
    { m_CurvatureWeight = c; }
  ScalarValueType GetCurvatureWeight() const
    { return m_CurvatureWeight; }
  
  void SetLaplacianSmoothingWeight(const ScalarValueType c)
    { m_LaplacianSmoothingWeight = c; }
  ScalarValueType GetLaplacianSmoothingWeight() const
    { return m_LaplacianSmoothingWeight; }
  
  void SetEpsilonMagnitude(const ScalarValueType e)
    { m_EpsilonMagnitude = e; }
  ScalarValueType GetEpsilonMagnitude() const
    { return m_EpsilonMagnitude; }

  virtual PixelType ComputeUpdate(const NeighborhoodType &neighborhood,
                                  void *globalData,
                                  const FloatOffsetType& = FloatOffsetType(0.0));

  virtual TimeStepType ComputeGlobalTimeStep(void *GlobalData) const;

  virtual void *GetGlobalDataPointer() const
    {
    GlobalDataStruct *ans = new GlobalDataStruct();
    ans->m_MaxAdvectionChange   = NumericTraits<ScalarValueType>::Zero;
    ans->m_MaxPropagationChange = NumericTraits<ScalarValueType>::Zero;
    ans->m_MaxCurvatureChange   = NumericTraits<ScalarValueType>::Zero;
    return ans; 
    }

  virtual void Initialize(const RadiusType &r);
  
  virtual void ReleaseGlobalDataPointer(void *GlobalData) const
    { delete (GlobalDataStruct *) GlobalData; }

  virtual ScalarValueType ComputeCurvatureTerm(const NeighborhoodType &,
                                               const FloatOffsetType &,
                                               GlobalDataStruct *gd = 0
                                               );
  virtual ScalarValueType ComputeMeanCurvature(const NeighborhoodType &,
                                               const FloatOffsetType &,
                                               GlobalDataStruct *gd = 0
                                               );

  virtual ScalarValueType ComputeMinimalCurvature(const NeighborhoodType &,
                                                  const FloatOffsetType &,
                                                  GlobalDataStruct *gd = 0
                                                  );
  
  virtual ScalarValueType Compute3DMinimalCurvature(const NeighborhoodType &,
                                                    const FloatOffsetType &,
                                                    GlobalDataStruct *gd = 0
                                                    );
  
  void SetUseMinimalCurvature( bool b )
    {
    m_UseMinimalCurvature = b;
    }
  bool GetUseMinimalCurvature() const
    {
    return m_UseMinimalCurvature;
    }
  void UseMinimalCurvatureOn()
    {
    this->SetUseMinimalCurvature(true);
    }
  void UseMinimalCurvatureOff()
    {
    this->SetUseMinimalCurvature(false);
    }

  static void SetMaximumCurvatureTimeStep(double n)
    {
    m_DT = n;
    }
  static double GetMaximumCurvatureTimeStep()
    {
    return m_DT;
    }

  static void SetMaximumPropagationTimeStep(double n)
    {
    m_WaveDT = n;
    }
  static double GetMaximumPropagationTimeStep()
    {
    return m_WaveDT;
    }

protected:
  LevelSetFunction()
    {
    m_EpsilonMagnitude = static_cast<ScalarValueType>( 1.0e-5 );
    m_AdvectionWeight = m_PropagationWeight 
      = m_CurvatureWeight = m_LaplacianSmoothingWeight 
      = NumericTraits<ScalarValueType>::Zero;
    m_UseMinimalCurvature = false;
    }
  virtual ~LevelSetFunction() {}
  void PrintSelf(std::ostream &s, Indent indent) const;
  
  static double m_WaveDT;
  static double m_DT;

  std::slice x_slice[itkGetStaticConstMacro(ImageDimension)];

  ::size_t m_Center;

  ::size_t m_xStride[itkGetStaticConstMacro(ImageDimension)];

  bool m_UseMinimalCurvature;

  static VectorType InitializeZeroVectorConstant();
  
  static VectorType m_ZeroVectorConstant;

  ScalarValueType m_EpsilonMagnitude;
  
  ScalarValueType m_AdvectionWeight;

  ScalarValueType m_PropagationWeight;

  ScalarValueType m_CurvatureWeight;

  ScalarValueType m_LaplacianSmoothingWeight;

private:
  LevelSetFunction(const Self&); //purposely not implemented
  void operator=(const Self&);   //purposely not implemented
};

} // namespace itk

// Define instantiation macro for this template.
#define ITK_TEMPLATE_LevelSetFunction(_, EXPORT, x, y) namespace itk { \
  _(1(class EXPORT LevelSetFunction< ITK_TEMPLATE_1 x >)) \
  namespace Templates { typedef LevelSetFunction< ITK_TEMPLATE_1 x > \
                                     LevelSetFunction##y; } \
  }

#if ITK_TEMPLATE_EXPLICIT
# include "Templates/itkLevelSetFunction+-.h"
#endif

#if ITK_TEMPLATE_TXX
# include "itkLevelSetFunction.txx"
#endif

#endif