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
 *
 *         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 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);
 
  itkTypeMacro(LevelSetFunction, FiniteDifferenceFunction);
 
  static constexpr unsigned int ImageDimension = Superclass::ImageDimension;
 
  using TimeStepType = double;
  using ImageType = typename Superclass::ImageType;
  using PixelType = typename Superclass::PixelType;
  using ScalarValueType = PixelType;
  using PixelRealType = typename Superclass::PixelRealType;
  using RadiusType = typename Superclass::RadiusType;
  using NeighborhoodType = typename Superclass::NeighborhoodType;
  using NeighborhoodScalesType = typename Superclass::NeighborhoodScalesType;
  using FloatOffsetType = 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 NumericTraits<ScalarValueType>::ZeroValue();
  }
 
  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 = NumericTraits<ScalarValueType>::ZeroValue();
    ans->m_MaxPropagationChange = NumericTraits<ScalarValueType>::ZeroValue();
    ans->m_MaxCurvatureChange = NumericTraits<ScalarValueType>::ZeroValue();
    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(NumericTraits<ScalarValueType>::ZeroValue())
    , m_PropagationWeight(NumericTraits<ScalarValueType>::ZeroValue())
    , m_CurvatureWeight(NumericTraits<ScalarValueType>::ZeroValue())
    , m_LaplacianSmoothingWeight(NumericTraits<ScalarValueType>::ZeroValue())
  {}
 
  ~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