itkRegionBasedLevelSetFunction.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 itkRegionBasedLevelSetFunction_h
#define itkRegionBasedLevelSetFunction_h
 
#include "itkFiniteDifferenceFunction.h"
#include "itkRegularizedHeavisideStepFunction.h"
#include "vnl/vnl_matrix_fixed.h"
 
namespace itk
{
template <typename TInput,   // LevelSetImageType
          typename TFeature, // FeatureImageType
          typename TSharedData>
class ITK_TEMPLATE_EXPORT RegionBasedLevelSetFunction : public FiniteDifferenceFunction<TInput>
{
public:
  ITK_DISALLOW_COPY_AND_MOVE(RegionBasedLevelSetFunction);

  using Self = RegionBasedLevelSetFunction;
  using Superclass = FiniteDifferenceFunction<TInput>;
  using Pointer = SmartPointer<Self>;
  using ConstPointer = SmartPointer<const Self>;
 
  static constexpr unsigned int ImageDimension = Superclass::ImageDimension;
 
  // itkNewMacro() is not provided since this is an abstract class.

  itkOverrideGetNameOfClassMacro(RegionBasedLevelSetFunction);

  using TimeStepType = double;
  using typename Superclass::ImageType;
  using typename Superclass::PixelType;
  using ScalarValueType = PixelType;
  using typename Superclass::RadiusType;
  using typename Superclass::NeighborhoodType;
  using typename Superclass::NeighborhoodScalesType;
  using typename Superclass::FloatOffsetType;
  using VectorType = FixedArray<ScalarValueType, Self::ImageDimension>;
 
  /* This structure is derived from LevelSetFunction and stores intermediate
  values for computing time step sizes */
  struct GlobalDataStruct
  {
    GlobalDataStruct()
    {
      ScalarValueType null_value{};
 
      m_MaxCurvatureChange = null_value;
      m_MaxAdvectionChange = null_value;
      m_MaxGlobalChange = null_value;
    }
 
    ~GlobalDataStruct() = default;
 
    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;
    ScalarValueType m_GradMag;
 
    ScalarValueType m_MaxCurvatureChange;
    ScalarValueType m_MaxAdvectionChange;
    ScalarValueType m_MaxGlobalChange;
  };
 
  using InputImageType = TInput;
  using InputImageConstPointer = typename InputImageType::ConstPointer;
  using InputImagePointer = typename InputImageType::Pointer;
  using InputPixelType = typename InputImageType::PixelType;
  using InputIndexType = typename InputImageType::IndexType;
  using InputIndexValueType = typename InputImageType::IndexValueType;
  using InputSizeType = typename InputImageType::SizeType;
  using InputSizeValueType = typename InputImageType::SizeValueType;
  using InputRegionType = typename InputImageType::RegionType;
  using InputPointType = typename InputImageType::PointType;
 
  using FeatureImageType = TFeature;
  using FeatureImageConstPointer = typename FeatureImageType::ConstPointer;
  using FeaturePixelType = typename FeatureImageType::PixelType;
  using FeatureIndexType = typename FeatureImageType::IndexType;
  using FeatureSpacingType = typename FeatureImageType::SpacingType;
  using FeatureOffsetType = typename FeatureImageType::OffsetType;
 
  using SharedDataType = TSharedData;
  using SharedDataPointer = typename SharedDataType::Pointer;
 
  using HeavisideFunctionType = HeavisideStepFunctionBase<InputPixelType, InputPixelType>;
  using HeavisideFunctionConstPointer = typename HeavisideFunctionType::ConstPointer;
 
  void
  SetDomainFunction(const HeavisideFunctionType * f)
  {
    this->m_DomainFunction = f;
  }
 
  virtual void
  Initialize(const RadiusType & r)
  {
    this->SetRadius(r);
 
    // Dummy neighborhood.
    NeighborhoodType it;
    it.SetRadius(r);
 
    // Find the center index of the neighborhood.
    m_Center = it.Size() / 2;
 
    // Get the stride length for each axis.
    for (unsigned int i = 0; i < ImageDimension; ++i)
    {
      m_xStride[i] = it.GetStride(i);
    }
  }
 
#if !defined(ITK_WRAPPING_PARSER)
  void
  SetSharedData(SharedDataPointer sharedDataIn)
  {
    this->m_SharedData = sharedDataIn;
  }
#endif
 
  void
  UpdateSharedData(bool forceUpdate);
 
  void *
  GetGlobalDataPointer() const override
  {
    return new GlobalDataStruct;
  }
 
  TimeStepType
  ComputeGlobalTimeStep(void * GlobalData) const override;

  PixelType
  ComputeUpdate(const NeighborhoodType & neighborhood,
                void *                   globalData,
                const FloatOffsetType & = FloatOffsetType(0.0)) override;
 
  void
  SetInitialImage(InputImageType * f)
  {
    m_InitialImage = f;
  }
 
  virtual const FeatureImageType *
  GetFeatureImage() const
  {
    return m_FeatureImage.GetPointer();
  }
  virtual void
  SetFeatureImage(const FeatureImageType * f)
  {
    m_FeatureImage = f;
 
    FeatureSpacingType spacing = m_FeatureImage->GetSpacing();
    for (unsigned int i = 0; i < ImageDimension; ++i)
    {
      this->m_InvSpacing[i] = 1 / spacing[i];
    }
  }

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

  void
  SetAreaWeight(const ScalarValueType & nu)
  {
    this->m_AreaWeight = nu;
  }
  ScalarValueType
  GetAreaWeight() const
  {
    return this->m_AreaWeight;
  }

  void
  SetLambda1(const ScalarValueType & lambda1)
  {
    this->m_Lambda1 = lambda1;
  }
  ScalarValueType
  GetLambda1() const
  {
    return this->m_Lambda1;
  }

  void
  SetLambda2(const ScalarValueType & lambda2)
  {
    this->m_Lambda2 = lambda2;
  }
  ScalarValueType
  GetLambda2() const
  {
    return this->m_Lambda2;
  }

  void
  SetOverlapPenaltyWeight(const ScalarValueType & gamma)
  {
    this->m_OverlapPenaltyWeight = gamma;
  }
  ScalarValueType
  GetOverlapPenaltyWeight() const
  {
    return this->m_OverlapPenaltyWeight;
  }

  virtual void
  SetCurvatureWeight(const ScalarValueType c)
  {
    m_CurvatureWeight = c;
  }
  ScalarValueType
  GetCurvatureWeight() const
  {
    return m_CurvatureWeight;
  }
 
  void
  SetAdvectionWeight(const ScalarValueType & iA)
  {
    this->m_AdvectionWeight = iA;
  }
  ScalarValueType
  GetAdvectionWeight() const
  {
    return this->m_AdvectionWeight;
  }

  void
  SetReinitializationSmoothingWeight(const ScalarValueType c)
  {
    m_ReinitializationSmoothingWeight = c;
  }
  ScalarValueType
  GetReinitializationSmoothingWeight() const
  {
    return m_ReinitializationSmoothingWeight;
  }

  void
  SetVolumeMatchingWeight(const ScalarValueType & tau)
  {
    this->m_VolumeMatchingWeight = tau;
  }
  ScalarValueType
  GetVolumeMatchingWeight() const
  {
    return this->m_VolumeMatchingWeight;
  }

  void
  SetVolume(const ScalarValueType & volume)
  {
    this->m_Volume = volume;
  }
  ScalarValueType
  GetVolume() const
  {
    return this->m_Volume;
  }

  void
  SetFunctionId(const unsigned int iFid)
  {
    this->m_FunctionId = iFid;
  }
 
  void
  ReleaseGlobalDataPointer(void * GlobalData) const override
  {
    delete (GlobalDataStruct *)GlobalData;
  }
 
  virtual ScalarValueType
  ComputeCurvature(const NeighborhoodType &, const FloatOffsetType &, GlobalDataStruct * gd);

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

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

  virtual void
  CalculateAdvectionImage()
  {}
 
protected:
  RegionBasedLevelSetFunction();
  ~RegionBasedLevelSetFunction() override = default;

  InputImageConstPointer m_InitialImage{};

  FeatureImageConstPointer m_FeatureImage{};
 
  SharedDataPointer m_SharedData{};
 
  HeavisideFunctionConstPointer m_DomainFunction{};

  ScalarValueType m_AreaWeight{};

  ScalarValueType m_Lambda1{};

  ScalarValueType m_Lambda2{};

  ScalarValueType m_OverlapPenaltyWeight{};

  ScalarValueType m_VolumeMatchingWeight{};

  ScalarValueType m_Volume{};

  ScalarValueType m_CurvatureWeight{};
 
  ScalarValueType m_AdvectionWeight{};

  ScalarValueType m_ReinitializationSmoothingWeight{};
 
  unsigned int m_FunctionId{};
 
  std::slice      x_slice[Self::ImageDimension];
  OffsetValueType m_Center{};
  OffsetValueType m_xStride[Self::ImageDimension]{};
  double          m_InvSpacing[Self::ImageDimension]{};
 
  static double m_WaveDT;
  static double m_DT;
 
  void
  ComputeHImage();

  ScalarValueType
  ComputeGlobalTerm(const ScalarValueType & imagePixel, const InputIndexType & inputIndex);

  virtual ScalarValueType
  ComputeInternalTerm(const FeaturePixelType & iValue, const FeatureIndexType & iIdx) = 0;

  virtual ScalarValueType
  ComputeExternalTerm(const FeaturePixelType & iValue, const FeatureIndexType & iIdx) = 0;

  virtual ScalarValueType
  ComputeOverlapParameters(const FeatureIndexType & featIndex, ScalarValueType & pr) = 0;

  ScalarValueType
  ComputeVolumeRegularizationTerm();


  ScalarValueType
  ComputeLaplacian(GlobalDataStruct * gd);

  void
  ComputeHessian(const NeighborhoodType & it, GlobalDataStruct * globalData);

  virtual void
  ComputeParameters() = 0;

  virtual void
  UpdateSharedDataParameters() = 0;
 
  bool m_UpdateC{};

  static VectorType
  InitializeZeroVectorConstant();

  static VectorType m_ZeroVectorConstant;
};
} // end namespace itk
 
#ifndef ITK_MANUAL_INSTANTIATION
#  include "itkRegionBasedLevelSetFunction.hxx"
#endif
 
#endif