itkParallelSparseFieldLevelSetImageFilter.h

Go to the documentation of this file.

/*=========================================================================
  
Program:   Insight Segmentation & Registration Toolkit
Module:    $RCSfile: itkParallelSparseFieldLevelSetImageFilter.h,v $
Language:  C++
Date:      $Date: 2009-04-25 12:27:47 $
Version:   $Revision: 1.18 $

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

  itkConceptMacro(OutputEqualityComparableCheck,
                  (Concept::EqualityComparable<PixelType>));
  itkConceptMacro(DoubleConvertibleToOutputCheck,
                  (Concept::Convertible<double, PixelType>));
  itkConceptMacro(OutputOStreamWritableCheck,
                  (Concept::OStreamWritable<PixelType>));

#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(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 ThreadedCalculateUpdateValue(const unsigned int itkNotUsed(ThreadId),
                                                        const IndexType itkNotUsed(index),
                                                        const TimeStepType &dt,
                                                        const ValueType &value,
                                                        const ValueType &change)
    {
    return (value + 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),
                                                        const unsigned int itkNotUsed(ThreadId));
  
  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);

  virtual void ThreadedInitializeIteration (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 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.txx"
#endif

#endif