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;
 
  const RadiusType &
  GetRadius() const
  {
    return m_Radius;
  }
 
  const unsigned int &
  GetArrayIndex(unsigned int i) const
  {
    return m_ArrayIndex[i];
  }
 
  const OffsetType &
  GetNeighborhoodOffset(unsigned int i) const
  {
    return m_NeighborhoodOffset[i];
  }
 
  const unsigned int &
  GetSize() const
  {
    return m_Size;
  }
 
  unsigned int
  GetStride(unsigned int i)
  {
    return m_StrideTable[i];
  }
 
  ParallelSparseFieldCityBlockNeighborList();
 
  ~ParallelSparseFieldCityBlockNeighborList()
  {
    m_ArrayIndex.clear();
    m_NeighborhoodOffset.clear();
  }
 
  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 StatusImageType = Image<StatusType, Self::ImageDimension>;
 
  using LayerNodeStorageType = ObjectStore<LayerNodeType>;
 
  using OffsetType = Offset<Self::ImageDimension>;
 
  itkSetMacro(NumberOfLayers, StatusType);
  itkGetConstMacro(NumberOfLayers, StatusType);
  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];
  }
 
#ifdef ITK_USE_CONCEPT_CHECKING
  // Begin concept checking
  itkConceptMacro(OutputEqualityComparableCheck, (Concept::EqualityComparable<PixelType>));
  itkConceptMacro(DoubleConvertibleToOutputCheck, (Concept::Convertible<double, PixelType>));
  itkConceptMacro(OutputOStreamWritableCheck, (Concept::OStreamWritable<PixelType>));
  // End concept checking
#endif
 
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{};
 
  StatusType 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