itkLevelSetFunction.h

Go to the documentation of this file.

/*=========================================================================
 *
 *  Copyright NumFOCUS
 *
 *  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
 *
 *         https://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 ITK_TEMPLATE_EXPORT LevelSetFunction : public FiniteDifferenceFunction<TImageType>
{
public:
  ITK_DISALLOW_COPY_AND_MOVE(LevelSetFunction);

  using Self = LevelSetFunction;
  using Superclass = FiniteDifferenceFunction<TImageType>;
  using Pointer = SmartPointer<Self>;
  using ConstPointer = SmartPointer<const Self>;

  itkNewMacro(Self);

  itkOverrideGetNameOfClassMacro(LevelSetFunction);

  static constexpr unsigned int ImageDimension = Superclass::ImageDimension;

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

  //  typedef
  using VectorType = FixedArray<ScalarValueType, Self::ImageDimension>;

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

    vnl_matrix_fixed<ScalarValueType, Self::ImageDimension, Self::ImageDimension> m_dxy;

    ScalarValueType m_dx[Self::ImageDimension];
 
    ScalarValueType m_dx_forward[Self::ImageDimension];
    ScalarValueType m_dx_backward[Self::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 ScalarValueType{};
  }

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

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

  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;
  }


  PixelType
  ComputeUpdate(const NeighborhoodType & it,
                void *                   globalData,
                const FloatOffsetType & = FloatOffsetType(0.0)) override;

  TimeStepType
  ComputeGlobalTimeStep(void * GlobalData) const override;

  void *
  GetGlobalDataPointer() const override
  {
    auto * ans = new GlobalDataStruct();
 
    ans->m_MaxAdvectionChange = ScalarValueType{};
    ans->m_MaxPropagationChange = ScalarValueType{};
    ans->m_MaxCurvatureChange = ScalarValueType{};
    return ans;
  }

  virtual void
  Initialize(const RadiusType & r);

  void
  ReleaseGlobalDataPointer(void * GlobalData) const override
  {
    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(ScalarValueType{})
    , m_PropagationWeight(ScalarValueType{})
    , m_CurvatureWeight(ScalarValueType{})
    , m_LaplacianSmoothingWeight(ScalarValueType{})
  {}
 
  ~LevelSetFunction() override = default;
  void
  PrintSelf(std::ostream & os, Indent indent) const override;

  static double m_WaveDT;
  static double m_DT;

  std::slice x_slice[Self::ImageDimension];

  OffsetValueType m_Center{ 0 };

  OffsetValueType m_xStride[Self::ImageDimension]{};
 
  bool m_UseMinimalCurvature{ false };

  static VectorType
  InitializeZeroVectorConstant();

  static VectorType m_ZeroVectorConstant;

  ScalarValueType m_EpsilonMagnitude{};

  ScalarValueType m_AdvectionWeight{};

  ScalarValueType m_PropagationWeight{};

  ScalarValueType m_CurvatureWeight{};

  ScalarValueType m_LaplacianSmoothingWeight{};
};
} // namespace itk
 
#ifndef ITK_MANUAL_INSTANTIATION
#  include "itkLevelSetFunction.hxx"
#endif
 
#endif