itkParallelSparseFieldLevelSetImageFilter.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 itkParallelSparseFieldLevelSetImageFilter_h
#define itkParallelSparseFieldLevelSetImageFilter_h

#include "itkFiniteDifferenceImageFilter.h"
#include "itkSparseFieldLayer.h"
#include "itkObjectStore.h"
#include "itkNeighborhoodIterator.h"
#include "itkMultiThreaderBase.h"
#include "itkBarrier.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) const;

private:
  char                        m_Pad1[128];
  unsigned int                m_Size;
  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_ASSIGN(ParallelSparseFieldLevelSetImageFilter);

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

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

  itkNewMacro(Self);

  itkTypeMacro(ParallelSparseFieldLevelSetImageFilter, FiniteDifferenceImageFilter);

  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();

  struct ParallelSparseFieldLevelSetThreadStruct {
    ParallelSparseFieldLevelSetImageFilter *Filter;
    std::vector< TimeStepType > TimeStepList;
    std::vector< bool > ValidTimeStepList;
    TimeStepType TimeStep;
  };

  void Iterate();

  static ITK_THREAD_RETURN_FUNCTION_CALL_CONVENTION IterateThreaderCallback(void *arg);

  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 outputimage 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 NumericTraits< TimeStepType >::ZeroValue();
  }

  virtual TimeStepType ThreadedCalculateChange(ThreadIdType ThreadId);

  void ThreadedUpdateActiveLayerValues( const TimeStepType& dt,
    LayerType   *StatusUpList,
    LayerType   *StatusDownList,
    ThreadIdType ThreadId);

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

  void ClearList(ThreadIdType ThreadId, LayerPointerType ListPtr);

  void CopyInsertInterNeighborNodeTransferBufferLayers(
    ThreadIdType ThreadId,
    LayerPointerType InputList,
    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();

  virtual void ThreadedLoadBalance(ThreadIdType ThreadId);

  void WaitForAll();

  void SignalNeighborsAndWait(ThreadIdType ThreadId);

  void SignalNeighbor(unsigned int SemaphoreArrayNumber, ThreadIdType ThreadId);

  void WaitForNeighbor(unsigned int SemaphoreArrayNumber, ThreadIdType ThreadId);

  virtual void ThreadedInitializeIteration(ThreadIdType ThreadId);

  //  void WriteActivePointsToFile ();

  ThreadIdType m_NumOfThreads{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};

  typename Barrier::Pointer m_Barrier;

  struct ThreadData {
    char pad1[128];

    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;

    char m_Pad2[128];
  };

  ThreadData *m_Data;

  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