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 ITK_TEMPLATE_EXPORT LevelSetFunction:
  public FiniteDifferenceFunction< TImageType >
{
public:
  ITK_DISALLOW_COPY_AND_ASSIGN(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 & neighborhood,
                                   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 & s, 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