itkMultiphaseSparseFiniteDifferenceImageFilter.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 itkMultiphaseSparseFiniteDifferenceImageFilter_h
#define itkMultiphaseSparseFiniteDifferenceImageFilter_h
 
#include "itkMultiphaseFiniteDifferenceImageFilter.h"
#include "itkZeroCrossingImageFilter.h"
#include "itkImageRegionConstIterator.h"
#include "itkShiftScaleImageFilter.h"
#include "itkNeighborhoodAlgorithm.h"
#include "itkSparseFieldLevelSetImageFilter.h"
 
#include <vector>
 
namespace itk
{
template <typename TInputImage,
          typename TFeatureImage,
          typename TOutputImage,
          typename TFunction,
          typename TIdCell = unsigned int>
class ITK_TEMPLATE_EXPORT MultiphaseSparseFiniteDifferenceImageFilter
  : public MultiphaseFiniteDifferenceImageFilter<TInputImage, TFeatureImage, TOutputImage, TFunction, TIdCell>
{
public:
  ITK_DISALLOW_COPY_AND_MOVE(MultiphaseSparseFiniteDifferenceImageFilter);

  using Self = MultiphaseSparseFiniteDifferenceImageFilter;
  using Superclass =
    MultiphaseFiniteDifferenceImageFilter<TInputImage, TFeatureImage, TOutputImage, TFunction, TIdCell>;
  using Pointer = SmartPointer<Self>;
  using ConstPointer = SmartPointer<const Self>;

  itkNewMacro(Self);

  itkOverrideGetNameOfClassMacro(MultiphaseSparseFiniteDifferenceImageFilter);
 
  static constexpr unsigned int ImageDimension = Superclass::ImageDimension;

  using typename Superclass::TimeStepType;

  using typename Superclass::InputImageType;
  using typename Superclass::InputImagePointer;
  using typename Superclass::InputRegionType;
  using typename Superclass::InputSizeType;
  using typename Superclass::InputSizeValueType;
  using typename Superclass::InputIndexType;
  using typename Superclass::InputIndexValueType;
  using typename Superclass::InputPixelType;
  using typename Superclass::InputPointType;
  using typename Superclass::InputSpacingType;
  using typename Superclass::InputOffsetValueType;
 
  using typename Superclass::FeatureImageType;
  using typename Superclass::FeatureSizeType;
  using typename Superclass::FeatureImagePointer;
  using typename Superclass::FeatureRegionType;
  using typename Superclass::FeatureSpacingType;
  using typename Superclass::FeaturePointType;
 
  using typename Superclass::OutputImageType;
  using typename Superclass::OutputImagePointer;
  using typename Superclass::OutputRegionType;
  using typename Superclass::OutputSizeType;
  using typename Superclass::OutputIndexType;
  using typename Superclass::OutputIndexValueType;
  using typename Superclass::OutputPixelType;
 
  using ValueType = typename InputImageType::ValueType;
  using typename Superclass::IdCellType;
 
  using typename Superclass::FiniteDifferenceFunctionType;
  using typename Superclass::FiniteDifferenceFunctionPointer;
  using FiniteDifferenceFunctionFloatOffsetType = typename FiniteDifferenceFunctionType::FloatOffsetType;

  using LayerNodeType = SparseFieldLevelSetNode<OutputIndexType>;

  using LayerType = SparseFieldLayer<LayerNodeType>;
  using LayerPointerType = typename LayerType::Pointer;
  using LayerIterator = typename LayerType::Iterator;
  using LayerConstIterator = typename LayerType::ConstIterator;

  using LayerListType = std::vector<LayerPointerType>;
  using LayerListIterator = typename LayerListType::iterator;
  using LayerListConstIterator = typename LayerListType::const_iterator;

  using StatusType = signed char;

  using StatusImageType = Image<StatusType, Self::ImageDimension>;
  using StatusImagePointer = typename StatusImageType::Pointer;
 
  using ZeroCrossingFilterType = ZeroCrossingImageFilter<InputImageType, InputImageType>;
  using ZeroCrossingFilterPointer = typename ZeroCrossingFilterType::Pointer;
 
  using BFCType = NeighborhoodAlgorithm::ImageBoundaryFacesCalculator<StatusImageType>;

  using LayerNodeStorageType = ObjectStore<LayerNodeType>;
  using LayerNodeStoragePointer = typename LayerNodeStorageType::Pointer;

  using UpdateBufferType = std::vector<ValueType>;
  using UpdateBufferConstIterator = typename UpdateBufferType::const_iterator;
 
  using NeighborListType = SparseFieldCityBlockNeighborList<NeighborhoodIterator<OutputImageType>>;
  using OffsetType = typename NeighborListType::OffsetType;

  itkSetMacro(NumberOfLayers, unsigned int);
  itkGetConstMacro(NumberOfLayers, unsigned int);

  itkSetMacro(IsoSurfaceValue, ValueType);
  itkGetConstMacro(IsoSurfaceValue, ValueType);

  itkSetMacro(InterpolateSurfaceLocation, bool);
  itkGetConstMacro(InterpolateSurfaceLocation, bool);
  itkBooleanMacro(InterpolateSurfaceLocation);
 
  void
  SetFunctionCount(const IdCellType & n)
  {
    this->Superclass::SetFunctionCount(n);
 
    m_SparseData.resize(this->m_FunctionCount, nullptr);
 
    for (IdCellType i = 0; i < this->m_FunctionCount; ++i)
    {
      m_SparseData[i] = new SparseDataStruct(i);
    }
  }
 
  itkConceptMacro(OutputEqualityComparableCheck, (Concept::EqualityComparable<typename TOutputImage::PixelType>));
  itkConceptMacro(DoubleConvertibleToOutputCheck, (Concept::Convertible<double, typename TOutputImage::PixelType>));
  itkConceptMacro(OutputOStreamWritableCheck, (Concept::OStreamWritable<typename TOutputImage::PixelType>));
 
protected:
  MultiphaseSparseFiniteDifferenceImageFilter();
  ~MultiphaseSparseFiniteDifferenceImageFilter() override
  {
    while (!m_SparseData.empty())
    {
      delete m_SparseData.back();
      m_SparseData.pop_back();
    }
  }
 
  void
  PrintSelf(std::ostream & os, Indent indent) const override;
 
  // This data structure is created for each phase
  struct SparseDataStruct
  {
    SparseDataStruct(const IdCellType & index)
    {
      m_LayerNodeStore = LayerNodeStorageType::New();
      m_LayerNodeStore->SetGrowthStrategyToExponential();
      m_Index = index;
    }

    LayerListType m_Layers;

    StatusImagePointer m_StatusImage;

    LayerNodeStoragePointer m_LayerNodeStore;

    UpdateBufferType m_UpdateBuffer;
 
    IdCellType m_Index;
  };

  NeighborListType m_NeighborList{};

  std::vector<ValueType> m_PixelDistance{};

  inline virtual ValueType
  CalculateUpdateValue(const OutputIndexType & itkNotUsed(idx),
                       const TimeStepType &    dt,
                       const ValueType &       value,
                       const ValueType &       change)
  {
    return (value + dt * change);
  }

  void
  PostProcessOutput() override;

  virtual void
  InitializeBackgroundPixels();

  void
  Initialize() override;

  void
  CopyInputToOutput() override;

  void
  AllocateUpdateBuffer() override
  {}

  void
  ApplyUpdate(TimeStepType dt) override;

  TimeStepType
  CalculateChange() override;

  void
  ConstructLayer(SparseDataStruct * sparsePtr, StatusType from, StatusType to);

  void
  ConstructActiveLayer();

  void
  InitializeActiveLayerValues();

  void
  InitializeBackgroundConstants();

  void
  PropagateLayerValues(SparseDataStruct * sparsePtr, StatusType from, StatusType to, StatusType promote, int InOrOut);

  void
  PropagateAllLayerValues();
 
  void
  PropagateFunctionLayerValues(unsigned int functionIndex);

  void
  UpdateActiveLayerValues(TimeStepType dt, LayerType * StatusUpList, LayerType * StatusDownList);

  void
  ProcessStatusList(LayerType * InputList,
                    LayerType * OutputList,
                    StatusType  ChangeToStatus,
                    StatusType  SearchForStatus);

  void
  ProcessOutsideList(LayerType * OutsideList, StatusType ChangeToStatus);
 
  void
  InitializeIteration() override;
 
  virtual void
  UpdatePixel(unsigned int                           itkNotUsed(functionIndex),
              unsigned int                           itkNotUsed(idx),
              NeighborhoodIterator<InputImageType> & itkNotUsed(iterator),
              ValueType &                            itkNotUsed(newValue),
              bool &                                 itkNotUsed(status))
  {}
 
  itkGetConstMacro(ValueZero, ValueType);
  itkGetConstMacro(ValueOne, ValueType);

  static double m_ConstantGradientValue;

  static const ValueType m_ValueOne;

  static const ValueType m_ValueZero;

  static const StatusType m_StatusChanging;

  static const StatusType m_StatusActiveChangingUp;

  static const StatusType m_StatusActiveChangingDown;

  static const StatusType m_StatusBoundaryPixel;

  static const StatusType m_StatusNull;
 
  std::vector<SparseDataStruct *> m_SparseData{};

  unsigned int m_NumberOfLayers{};

  ValueType m_IsoSurfaceValue{};

  ValueType m_BackgroundValue{};

  bool m_InterpolateSurfaceLocation{};
 
private:
  unsigned int m_CurrentFunctionIndex{};
 
  double       m_RMSSum{ 0. };
  unsigned int m_RMSCounter{ 0 };

  bool m_BoundsCheckingActive{};
};
} // end namespace itk
 
#ifndef ITK_MANUAL_INSTANTIATION
#  include "itkMultiphaseSparseFiniteDifferenceImageFilter.hxx"
#endif
 
#endif