itkParallelSparseFieldLevelSetImageFilter.h

Go to the documentation of this file.

/*=========================================================================
  
Program:   Insight Segmentation & Registration Toolkit
Module:    $RCSfile: itkParallelSparseFieldLevelSetImageFilter.h,v $
Language:  C++
Date:      $Date: 2003/09/10 14:28:54 $
Version:   $Revision: 1.7 $

Copyright (c) Insight Software Consortium. All rights reserved.
See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details.

This software is distributed WITHOUT ANY WARRANTY; without even 
the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR 
PURPOSE.  See the above copyright notices for more information.

=========================================================================*/
#ifndef __itkParallelSparseFieldLevelSetImageFilter_h_
#define __itkParallelSparseFieldLevelSetImageFilter_h_

#include <vector>
#include "itkFiniteDifferenceImageFilter.h"
#include "itkSparseFieldLayer.h"
#include "itkObjectStore.h"
#include "itkNeighborhoodIterator.h"
#include "itkConstNeighborhoodIterator.h"
#include "itkMultiThreader.h"
#include "itkBarrier.h"
#include "itkSemaphore.h"

namespace itk
{

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

template <class 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 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 pad2[128];
};

template <class TInputImage, class TOutputImage>
class ITK_EXPORT 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;
  
  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);
  itkGetMacro(NumberOfLayers, StatusType);

  itkSetMacro(IsoSurfaceValue, ValueType);
  itkGetMacro(IsoSurfaceValue, ValueType);

  itkGetMacro(RMSChange, ValueType);

protected:
  ParallelSparseFieldLevelSetImageFilter();
  ~ParallelSparseFieldLevelSetImageFilter() {}
  virtual void PrintSelf(std::ostream& os, Indent indent) const;

  ParallelSparseFieldCityBlockNeighborList < NeighborhoodIterator<OutputImageType> >
    m_NeighborList;
  
  static 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(StatusType from, StatusType to);
  
  void ProcessStatusList(LayerType *InputList, StatusType ChangeToStatus,
                         StatusType SearchForStatus, unsigned int ThreadId);
  
  void PropagateAllLayerValues();
  
  void PropagateLayerValues(StatusType from, StatusType to, StatusType promote,
                            unsigned int InOrOut);
  
  virtual void InitializeBackgroundPixels();
  
  void ThreadedAllocateData(unsigned int ThreadId);
  void ThreadedInitializeData(unsigned int ThreadId, const ThreadRegionType & ThreadRegion);
  
  void ComputeInitialThreadBoundaries();
  
  unsigned int GetThreadNumber(unsigned int splitAxisValue);
  
  void GetThreadRegionSplitByBoundary(unsigned int ThreadId, ThreadRegionType& ThreadRegion);
  
  void DeallocateData();

  struct ParallelSparseFieldLevelSetThreadStruct
  {
    ParallelSparseFieldLevelSetImageFilter* Filter;
    TimeStepType* TimeStepList;
    bool* ValidTimeStepList;
    TimeStepType TimeStep;
  };
  
  void Iterate();
  static ITK_THREAD_RETURN_TYPE IterateThreaderCallback(void * arg);
  
  inline virtual ValueType CalculateUpdateValue(const IndexType,
                                                const TimeStepType &dt,
                                                const ValueType &value,
                                                const ValueType &change)
    {
      return (value + dt * change);
    }
  
  void ApplyUpdate(TimeStepType) {}

  virtual void ThreadedApplyUpdate(TimeStepType dt, unsigned int ThreadId);
  
  TimeStepType CalculateChange()
  {
    return NumericTraits<TimeStepType>::Zero;
  }
  
  virtual TimeStepType ThreadedCalculateChange(unsigned int ThreadId);

  void ThreadedUpdateActiveLayerValues(TimeStepType dt, LayerType *StatusUpList,
                                       LayerType *StatusDownList, unsigned int ThreadId);
  
  void CopyInsertList(unsigned int ThreadId, LayerPointerType FromListPtr,
                      LayerPointerType ToListPtr);

  void ClearList(unsigned int ThreadId, LayerPointerType ListPtr);

  void CopyInsertInterNeighborNodeTransferBufferLayers(unsigned int ThreadId,
                                                       LayerPointerType InputList,
                                                       unsigned int InOrOut,
                                                       unsigned int BufferLayerNumber);
  
  void ClearInterNeighborNodeTransferBufferLayers(unsigned int ThreadId, unsigned int InOrOut,
                                                  unsigned int BufferLayerNumber);
  
  void ThreadedProcessFirstLayerStatusLists(unsigned int InputLayerNumber,
                                            unsigned int OutputLayerNumber,
                                            StatusType SearchForStatus,
                                            unsigned int InOrOut,
                                            unsigned int BufferLayerNumber, unsigned int ThreadId);
  
  void ThreadedProcessStatusList(unsigned int InputLayerNumber, unsigned int OutputLayerNumber,
                                 StatusType ChangeToStatus, StatusType SearchForStatus,
                                 unsigned int InOrOut,
                                 unsigned int BufferLayerNumber, unsigned int ThreadId);
  
  void ThreadedProcessOutsideList(unsigned int InputLayerNumber, StatusType ChangeToStatus,
                                  unsigned int InOrOut, unsigned int BufferLayerNumber, unsigned int ThreadId);
  
  void ThreadedPropagateLayerValues(StatusType from, StatusType to, StatusType promote,
                                    unsigned int InorOut, unsigned int ThreadId);
  
  void GetThreadRegionSplitUniformly(unsigned int ThreadId, ThreadRegionType& ThreadRegion);
  
  void ThreadedPostProcessOutput(const ThreadRegionType & regionToProcess);
  
  virtual void CheckLoadBalance();
  
  virtual void ThreadedLoadBalance(unsigned int ThreadId);
  
  void WaitForAll();
  void SignalNeighborsAndWait (unsigned int ThreadId);
  void SignalNeighbor  (unsigned int SemaphoreArrayNumber, unsigned int ThreadId);
  void WaitForNeighbor (unsigned int SemaphoreArrayNumber, unsigned int ThreadId);
  
  //  void WriteActivePointsToFile ();
  
  unsigned int 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;
    
    typename Semaphore::Pointer m_Semaphore[2];
    
    unsigned int m_SemaphoreArrayNumber;
    
    char 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.txx"
#endif

#endif

---