itkLevelSetFunction.h

Go to the documentation of this file.

/*=========================================================================
 *
 *  Copyright Insight Software Consortium
 *
 *  Licensed under the Apache License, Version 2.0 (the "License");
 *  you may not use this file except in compliance with the License.
 *  You may obtain a copy of the License at
 *
 *         http://www.apache.org/licenses/LICENSE-2.0.txt
 *
 *  Unless required by applicable law or agreed to in writing, software
 *  distributed under the License is distributed on an "AS IS" BASIS,
 *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 *  See the License for the specific language governing permissions and
 *  limitations under the License.
 *
 *=========================================================================*/
#ifndef __itkLevelSetFunction_h
#define __itkLevelSetFunction_h

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

namespace itk
{
template< typename TImageType >
class 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_Center(0),
    m_UseMinimalCurvature(false),
    m_EpsilonMagnitude(static_cast< ScalarValueType >( 1.0e-5 )),
    m_AdvectionWeight(NumericTraits< ScalarValueType >::Zero),
    m_PropagationWeight(NumericTraits< ScalarValueType >::Zero),
    m_CurvatureWeight(NumericTraits< ScalarValueType >::Zero),
    m_LaplacianSmoothingWeight(NumericTraits< ScalarValueType >::Zero)
  {
  }

  virtual ~LevelSetFunction() {}
  void PrintSelf(std::ostream & s, Indent indent) const;

  static double m_WaveDT;
  static double m_DT;

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

  OffsetValueType m_Center;

  OffsetValueType 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

#ifndef ITK_MANUAL_INSTANTIATION
#include "itkLevelSetFunction.hxx"
#endif

#endif