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 <vector>
#include "itkFiniteDifferenceImageFilter.h"
#include "itkSparseFieldLayer.h"
#include "itkObjectStore.h"
#include "itkNeighborhoodIterator.h"
#include "itkMultiThreader.h"
#include "itkBarrier.h"

namespace itk
{
template< typename TNodeIndexType >
class ParallelSparseFieldLevelSetNode
{
public:
  TNodeIndexType                   m_Index;
  float                            m_Value;
  ParallelSparseFieldLevelSetNode *Next;
  ParallelSparseFieldLevelSetNode *Previous;
};

template< typename TNeighborhoodType >
class ParallelSparseFieldCityBlockNeighborList
{
public:
  typedef TNeighborhoodType                     NeighborhoodType;
  typedef typename NeighborhoodType::OffsetType OffsetType;
  typedef typename NeighborhoodType::RadiusType RadiusType;

  itkStaticConstMacro(Dimension, unsigned int, 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 ParallelSparseFieldLevelSetImageFilter:
  public FiniteDifferenceImageFilter< TInputImage, TOutputImage >
{
public:

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

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

  itkNewMacro(Self);

  itkTypeMacro(ParallelSparseFieldLevelSetImageFilter, FiniteDifferenceImageFilter);

  typedef TInputImage                         InputImageType;
  typedef TOutputImage                        OutputImageType;
  typedef typename OutputImageType::IndexType IndexType;

  itkStaticConstMacro(ImageDimension, unsigned int, TOutputImage::ImageDimension);

  typedef typename OutputImageType::PixelType PixelType;

  typedef typename OutputImageType::RegionType ThreadRegionType;

  typedef typename OutputImageType::ValueType ValueType;

  typedef ParallelSparseFieldLevelSetNode< IndexType > LayerNodeType;

  typedef SparseFieldLayer< LayerNodeType > LayerType;
  typedef typename LayerType::Pointer       LayerPointerType;

  typedef std::vector< LayerPointerType > LayerListType;

  typedef signed char StatusType;

  typedef Image< StatusType, itkGetStaticConstMacro(ImageDimension) > StatusImageType;

  typedef ObjectStore< LayerNodeType > LayerNodeStorageType;

  typedef Offset< itkGetStaticConstMacro(ImageDimension) > OffsetType;

  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() {}
  virtual void PrintSelf(std::ostream & os, Indent indent) const;

  ParallelSparseFieldCityBlockNeighborList< NeighborhoodIterator< OutputImageType > >
  m_NeighborList;

  double m_ConstantGradientValue;

  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;

  virtual void GenerateData();

  void CopyInputToOutput();

  void AllocateUpdateBuffer() {}

  void Initialize();

  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_TYPE 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&) {}

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

  TimeStepType CalculateChange()
  {
    return NumericTraits< TimeStepType >::Zero;
  }

  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;

  unsigned int m_SplitAxis;

  unsigned int m_ZSize;

  bool m_BoundaryChanged;

  unsigned int *m_Boundary;

  int *m_GlobalZHistogram;

  unsigned int *m_MapZToThreadNumber;

  int *m_ZCumulativeFrequency;

  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];

    SimpleMutexLock            m_Lock[2];
    ConditionVariable::Pointer m_Condition[2];

    unsigned int m_SemaphoreArrayNumber;

    char m_Pad2[128];
  };

  ThreadData *m_Data;

  bool m_Stop;

  bool m_InterpolateSurfaceLocation;

private:

  ParallelSparseFieldLevelSetImageFilter(const Self &); // purposely not
                                                        // implemented
  void operator=(const Self &);                         // purposely not

  // implemented

  bool m_BoundsCheckingActive;
};
} // end namespace itk

#ifndef ITK_MANUAL_INSTANTIATION
#include "itkParallelSparseFieldLevelSetImageFilter.hxx"
#endif

#endif