itkParallelSparseFieldLevelSetImageFilter.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 itkParallelSparseFieldLevelSetImageFilter_h
#define itkParallelSparseFieldLevelSetImageFilter_h
 
#include "itkBooleanStdVector.h"
#include "itkFiniteDifferenceImageFilter.h"
#include "itkSparseFieldLayer.h"
#include "itkObjectStore.h"
#include "itkNeighborhoodIterator.h"
#include "itkMultiThreaderBase.h"
#include <condition_variable>
#include <vector>
 
namespace itk
{
template <typename TNodeIndexType>
class ITK_TEMPLATE_EXPORT ParallelSparseFieldLevelSetNode
{
public:
  TNodeIndexType                    m_Index;
  float                             m_Value;
  ParallelSparseFieldLevelSetNode * Next;
  ParallelSparseFieldLevelSetNode * Previous;
};

template <typename TNeighborhoodType>
class ITK_TEMPLATE_EXPORT ParallelSparseFieldCityBlockNeighborList
{
public:
  using NeighborhoodType = TNeighborhoodType;
  using OffsetType = typename NeighborhoodType::OffsetType;
  using RadiusType = typename NeighborhoodType::RadiusType;
 
  static constexpr unsigned int Dimension = NeighborhoodType::Dimension;
 
  [[nodiscard]] const RadiusType &
  GetRadius() const
  {
    return m_Radius;
  }
 
  [[nodiscard]] const unsigned int &
  GetArrayIndex(unsigned int i) const
  {
    return m_ArrayIndex[i];
  }
 
  [[nodiscard]] const OffsetType &
  GetNeighborhoodOffset(unsigned int i) const
  {
    return m_NeighborhoodOffset[i];
  }
 
  [[nodiscard]] const unsigned int &
  GetSize() const
  {
    return m_Size;
  }
 
  unsigned int
  GetStride(unsigned int i)
  {
    return m_StrideTable[i];
  }
 
  ParallelSparseFieldCityBlockNeighborList();
 
  ~ParallelSparseFieldCityBlockNeighborList() = default;
 
  void
  Print(std::ostream & os, Indent indent) const;
 
private:
  char                      m_Pad1[128]{};
  unsigned int              m_Size{ 2 * Dimension };
  RadiusType                m_Radius{};
  std::vector<unsigned int> m_ArrayIndex{};
  std::vector<OffsetType>   m_NeighborhoodOffset{};

  unsigned int m_StrideTable[Dimension]{};
  char         m_Pad2[128]{};
};

template <typename TInputImage, typename TOutputImage>
class ITK_TEMPLATE_EXPORT ParallelSparseFieldLevelSetImageFilter
  : public FiniteDifferenceImageFilter<TInputImage, TOutputImage>
{
public:
  ITK_DISALLOW_COPY_AND_MOVE(ParallelSparseFieldLevelSetImageFilter);

  using Self = ParallelSparseFieldLevelSetImageFilter;
  using Superclass = FiniteDifferenceImageFilter<TInputImage, TOutputImage>;
  using Pointer = SmartPointer<Self>;
  using ConstPointer = SmartPointer<const Self>;

  using typename Superclass::TimeStepType;
  using typename Superclass::FiniteDifferenceFunctionType;
  using typename Superclass::RadiusType;
  using typename Superclass::NeighborhoodScalesType;

  itkNewMacro(Self);

  itkOverrideGetNameOfClassMacro(ParallelSparseFieldLevelSetImageFilter);

  using InputImageType = TInputImage;
  using OutputImageType = TOutputImage;
  using IndexType = typename OutputImageType::IndexType;
 
  static constexpr unsigned int ImageDimension = TOutputImage::ImageDimension;
 
  using PixelType = typename OutputImageType::PixelType;
 
  using ThreadRegionType = typename OutputImageType::RegionType;

  using ValueType = typename OutputImageType::ValueType;

  using LayerNodeType = ParallelSparseFieldLevelSetNode<IndexType>;

  using LayerType = SparseFieldLayer<LayerNodeType>;
  using LayerPointerType = typename LayerType::Pointer;

  using LayerListType = std::vector<LayerPointerType>;

  using StatusType = signed char;

  using LayerCountType = StatusType; /* purposefully the same as StatusType
                                        for backward compatibility. This value
                                        is compared with up/down layer values
                                        that can be signed.  Max layer count of
                                        127 is enough. */

  using StatusImageType = Image<StatusType, Self::ImageDimension>;

  using LayerNodeStorageType = ObjectStore<LayerNodeType>;
 
  using OffsetType = Offset<Self::ImageDimension>;

  itkSetMacro(NumberOfLayers, LayerCountType);
  itkGetConstMacro(NumberOfLayers, LayerCountType);
  itkSetMacro(IsoSurfaceValue, ValueType);
  itkGetConstMacro(IsoSurfaceValue, ValueType);
  LayerPointerType
  GetActiveListForIndex(const IndexType index)
  {
    // get the 'z' value for the index
    const unsigned int indexZ = index[m_SplitAxis];
    // get the thread in whose region the index lies
    const unsigned int ThreadNum = this->GetThreadNumber(indexZ);
 
    // get the active list for that thread
    return m_Data[ThreadNum].m_Layers[0];
  }
 
  itkConceptMacro(OutputEqualityComparableCheck, (Concept::EqualityComparable<PixelType>));
  itkConceptMacro(DoubleConvertibleToOutputCheck, (Concept::Convertible<double, PixelType>));
  itkConceptMacro(OutputOStreamWritableCheck, (Concept::OStreamWritable<PixelType>));
 
protected:
  ParallelSparseFieldLevelSetImageFilter();
  ~ParallelSparseFieldLevelSetImageFilter() override = default;
  void
  PrintSelf(std::ostream & os, Indent indent) const override;

  ParallelSparseFieldCityBlockNeighborList<NeighborhoodIterator<OutputImageType>> m_NeighborList{};

  double m_ConstantGradientValue{ 1.0 };

  static ValueType m_ValueOne;

  static ValueType m_ValueZero;

  static StatusType m_StatusActiveChangingUp;

  static StatusType m_StatusActiveChangingDown;

  static StatusType m_StatusNull;

  static StatusType m_StatusChanging;

  static StatusType m_StatusBoundaryPixel;

  typename OutputImageType::Pointer m_ShiftedImage{};

  LayerListType m_Layers{};

  LayerCountType m_NumberOfLayers{};

  typename StatusImageType::Pointer m_StatusImage{};
  typename OutputImageType::Pointer m_OutputImage{};

  typename StatusImageType::Pointer m_StatusImageTemp{};
  typename OutputImageType::Pointer m_OutputImageTemp{};

  typename LayerNodeStorageType::Pointer m_LayerNodeStore{};

  ValueType m_IsoSurfaceValue{};

  //  ValueType m_RMSChange;

  void
  GenerateData() override;

  void
  CopyInputToOutput() override;

  void
  AllocateUpdateBuffer() override
  {}

  void
  Initialize() override;

  void
  ConstructActiveLayer();

  void
  InitializeActiveLayerValues();

  void
  ConstructLayer(const StatusType & from, const StatusType & to);

  void
  ProcessStatusList(LayerType *        InputList,
                    const StatusType & ChangeToStatus,
                    const StatusType & SearchForStatus,
                    ThreadIdType       ThreadId);

  void
  PropagateAllLayerValues();

  void
  PropagateLayerValues(const StatusType & from,
                       const StatusType & to,
                       const StatusType & promote,
                       unsigned int       InOrOut);

  virtual void
  InitializeBackgroundPixels();

  void
  ThreadedAllocateData(ThreadIdType ThreadId);
 
  void
  ThreadedInitializeData(ThreadIdType ThreadId, const ThreadRegionType & ThreadRegion);

  void
  ComputeInitialThreadBoundaries();

  unsigned int
  GetThreadNumber(unsigned int splitAxisValue);

  void
  GetThreadRegionSplitByBoundary(ThreadIdType ThreadId, ThreadRegionType & ThreadRegion);

  void
  DeallocateData();

  void
  Iterate();

  inline virtual ValueType
  ThreadedCalculateUpdateValue(const ThreadIdType   itkNotUsed(ThreadId),
                               const IndexType      itkNotUsed(index),
                               const TimeStepType & dt,
                               const ValueType &    value,
                               const ValueType &    change)
  {
    return value + static_cast<ValueType>(dt) * change;
  }
 
  // This method can be overridden in derived classes.
  // The pixel at 'index' is entering the active layer for thread 'ThreadId'.
  // The output image at 'index' will have the value as given by the
  // 'value' parameter.
  virtual void
  ThreadedProcessPixelEnteringActiveLayer(const IndexType & itkNotUsed(index),
                                          const ValueType & itkNotUsed(value),
                                          ThreadIdType      itkNotUsed(ThreadId));

  void
  ApplyUpdate(const TimeStepType &) override
  {}

  virtual void
  ThreadedApplyUpdate(const TimeStepType & dt, ThreadIdType ThreadId);

  TimeStepType
  CalculateChange() override
  {
    return TimeStepType{};
  }

  virtual TimeStepType
  ThreadedCalculateChange(ThreadIdType ThreadId);

  void
  ThreadedUpdateActiveLayerValues(const TimeStepType & dt,
                                  LayerType *          UpList,
                                  LayerType *          DownList,
                                  ThreadIdType         ThreadId);

  void
  CopyInsertList(ThreadIdType ThreadId, LayerPointerType FromListPtr, LayerPointerType ToListPtr);

  void
  ClearList(ThreadIdType ThreadId, LayerPointerType ListPtr);

  void
  CopyInsertInterNeighborNodeTransferBufferLayers(ThreadIdType     ThreadId,
                                                  LayerPointerType List,
                                                  unsigned int     InOrOut,
                                                  unsigned int     BufferLayerNumber);

  void
  ClearInterNeighborNodeTransferBufferLayers(ThreadIdType ThreadId,
                                             unsigned int InOrOut,
                                             unsigned int BufferLayerNumber);

  void
  ThreadedProcessFirstLayerStatusLists(unsigned int       InputLayerNumber,
                                       unsigned int       OutputLayerNumber,
                                       const StatusType & SearchForStatus,
                                       unsigned int       InOrOut,
                                       unsigned int       BufferLayerNumber,
                                       ThreadIdType       ThreadId);

  void
  ThreadedProcessStatusList(unsigned int       InputLayerNumber,
                            unsigned int       OutputLayerNumber,
                            const StatusType & ChangeToStatus,
                            const StatusType & SearchForStatus,
                            unsigned int       InOrOut,
                            unsigned int       BufferLayerNumber,
                            ThreadIdType       ThreadId);

  void
  ThreadedProcessOutsideList(unsigned int       InputLayerNumber,
                             const StatusType & ChangeToStatus,
                             unsigned int       InOrOut,
                             unsigned int       BufferLayerNumber,
                             ThreadIdType       ThreadId);

  void
  ThreadedPropagateLayerValues(const StatusType & from,
                               const StatusType & to,
                               const StatusType & promote,
                               unsigned int       InOrOut,
                               ThreadIdType       ThreadId);

  void
  GetThreadRegionSplitUniformly(ThreadIdType ThreadId, ThreadRegionType & ThreadRegion);

  void
  ThreadedPostProcessOutput(const ThreadRegionType & regionToProcess);

  virtual void
  CheckLoadBalance();

  void
  ThreadedLoadBalance1(ThreadIdType ThreadId);
  void
  ThreadedLoadBalance2(ThreadIdType ThreadId);
  void
  SignalNeighborsAndWait(ThreadIdType ThreadId);
 
  void
  SignalNeighbor(unsigned int SemaphoreArrayNumber, ThreadIdType ThreadId);
 
  void
  WaitForNeighbor(unsigned int SemaphoreArrayNumber, ThreadIdType ThreadId);

  virtual void
  ThreadedInitializeIteration(ThreadIdType ThreadId);

  std::vector<TimeStepType> m_TimeStepList{};
  BooleanStdVectorType      m_ValidTimeStepList{};
  TimeStepType              m_TimeStep{};

  ThreadIdType m_NumOfWorkUnits{ 0 };

  unsigned int m_SplitAxis{ 0 };

  unsigned int m_ZSize{ 0 };

  bool m_BoundaryChanged{ false };

  unsigned int * m_Boundary{ nullptr };

  int * m_GlobalZHistogram{ nullptr };

  unsigned int * m_MapZToThreadNumber{ nullptr };

  int * m_ZCumulativeFrequency{ nullptr };

  struct ThreadDataUnaligned
  {
    TimeStepType     TimeStep;
    ThreadRegionType ThreadRegion;
    ValueType        m_RMSChange;
    unsigned int     m_Count;

    LayerListType m_Layers;

    LayerListType * m_LoadTransferBufferLayers;

    typename LayerNodeStorageType::Pointer m_LayerNodeStore;
 
    LayerPointerType UpList[2];
    LayerPointerType DownList[2];

    LayerPointerType ** m_InterNeighborNodeTransferBufferLayers[2];

    void * globalData;

    int * m_ZHistogram;

    int m_Semaphore[2];
 
    std::mutex              m_Lock[2];
    std::condition_variable m_Condition[2];

    unsigned int m_SemaphoreArrayNumber;
  };
 
  itkPadStruct(ITK_CACHE_LINE_ALIGNMENT, ThreadDataUnaligned, ThreadDataPadded);
  itkAlignedTypedef(ITK_CACHE_LINE_ALIGNMENT, ThreadDataPadded, ThreadData);

  ThreadData * m_Data{ nullptr };

  bool m_Stop{ false };

  bool m_InterpolateSurfaceLocation{ true };
 
private:
  bool m_BoundsCheckingActive{ false };
};
} // end namespace itk
 
#ifndef ITK_MANUAL_INSTANTIATION
#  include "itkParallelSparseFieldLevelSetImageFilter.hxx"
#endif
 
#endif