00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017 #ifndef __itkParallelSparseFieldLevelSetImageFilter_h
00018 #define __itkParallelSparseFieldLevelSetImageFilter_h
00019
00020 #include <vector>
00021 #include "itkFiniteDifferenceImageFilter.h"
00022 #include "itkSparseFieldLayer.h"
00023 #include "itkObjectStore.h"
00024 #include "itkNeighborhoodIterator.h"
00025 #include "itkConstNeighborhoodIterator.h"
00026 #include "itkMultiThreader.h"
00027 #include "itkBarrier.h"
00028 #include "itkSemaphore.h"
00029
00030 namespace itk
00031 {
00032
00037 template <class TNodeIndexType>
00038 class ParallelSparseFieldLevelSetNode
00039 {
00040 public:
00041 TNodeIndexType m_Index;
00042 float m_Value;
00043 ParallelSparseFieldLevelSetNode *Next;
00044 ParallelSparseFieldLevelSetNode *Previous;
00045 };
00046
00073 template <class TNeighborhoodType>
00074 class ParallelSparseFieldCityBlockNeighborList
00075 {
00076 public:
00077 typedef TNeighborhoodType NeighborhoodType;
00078 typedef typename NeighborhoodType::OffsetType OffsetType;
00079 typedef typename NeighborhoodType::RadiusType RadiusType;
00080
00081 itkStaticConstMacro(Dimension, unsigned int, NeighborhoodType::Dimension);
00082
00083 const RadiusType &GetRadius() const
00084 {
00085 return m_Radius;
00086 }
00087
00088 const unsigned int &GetArrayIndex(unsigned int i) const
00089 {
00090 return m_ArrayIndex[i];
00091 }
00092
00093 const OffsetType &GetNeighborhoodOffset(unsigned int i) const
00094 {
00095 return m_NeighborhoodOffset[i];
00096 }
00097
00098 const unsigned int &GetSize() const
00099 {
00100 return m_Size;
00101 }
00102
00103 unsigned int GetStride(unsigned int i)
00104 {
00105 return m_StrideTable[i];
00106 }
00107
00108 ParallelSparseFieldCityBlockNeighborList();
00109
00110 ~ParallelSparseFieldCityBlockNeighborList()
00111 {
00112 m_ArrayIndex.clear();
00113 m_NeighborhoodOffset.clear();
00114 }
00115
00116 void Print(std::ostream &os) const;
00117
00118 private:
00119 char m_Pad1[128];
00120 unsigned int m_Size;
00121 RadiusType m_Radius;
00122 std::vector<unsigned int> m_ArrayIndex;
00123 std::vector<OffsetType> m_NeighborhoodOffset;
00124
00127 unsigned int m_StrideTable[Dimension];
00128 char m_Pad2[128];
00129 };
00130
00246 template <class TInputImage, class TOutputImage>
00247 class ITK_EXPORT ParallelSparseFieldLevelSetImageFilter :
00248 public FiniteDifferenceImageFilter<TInputImage, TOutputImage>
00249 {
00250 public:
00251
00253 typedef ParallelSparseFieldLevelSetImageFilter Self;
00254 typedef FiniteDifferenceImageFilter<TInputImage, TOutputImage> Superclass;
00255 typedef SmartPointer<Self> Pointer;
00256 typedef SmartPointer<const Self> ConstPointer;
00257
00259 typedef typename Superclass::TimeStepType TimeStepType;
00260 typedef typename Superclass::FiniteDifferenceFunctionType FiniteDifferenceFunctionType;
00261 typedef typename Superclass::RadiusType RadiusType;
00262 typedef typename Superclass::NeighborhoodScalesType NeighborhoodScalesType;
00263
00265 itkNewMacro(Self);
00266
00268 itkTypeMacro(ParallelSparseFieldLevelSetImageFilter, FiniteDifferenceImageFilter);
00269
00271 typedef TInputImage InputImageType;
00272 typedef TOutputImage OutputImageType;
00273 typedef typename OutputImageType::IndexType IndexType;
00274
00275 itkStaticConstMacro(ImageDimension, unsigned int, TOutputImage::ImageDimension);
00276
00277 typedef typename OutputImageType::PixelType PixelType;
00278
00279 typedef typename OutputImageType::RegionType ThreadRegionType;
00280
00283 typedef typename OutputImageType::ValueType ValueType;
00284
00286 typedef ParallelSparseFieldLevelSetNode<IndexType> LayerNodeType;
00287
00289 typedef SparseFieldLayer<LayerNodeType> LayerType;
00290 typedef typename LayerType::Pointer LayerPointerType;
00291
00293 typedef std::vector<LayerPointerType> LayerListType;
00294
00296 typedef signed char StatusType;
00297
00300 typedef Image<StatusType, itkGetStaticConstMacro(ImageDimension)> StatusImageType;
00301
00304 typedef ObjectStore<LayerNodeType> LayerNodeStorageType;
00305
00306 typedef Offset<itkGetStaticConstMacro(ImageDimension)> OffsetType;
00307
00311 itkSetMacro(NumberOfLayers, StatusType);
00312 itkGetMacro(NumberOfLayers, StatusType);
00314
00316 itkSetMacro(IsoSurfaceValue, ValueType);
00317 itkGetMacro(IsoSurfaceValue, ValueType);
00319
00320 LayerPointerType GetActiveListForIndex (const IndexType index)
00321 {
00322
00323 const unsigned int indexZ= index[m_SplitAxis];
00324
00325 const unsigned int ThreadNum= this->GetThreadNumber (indexZ);
00326
00327 return m_Data[ThreadNum].m_Layers[0];
00328 }
00329
00330 #ifdef ITK_USE_CONCEPT_CHECKING
00331
00332 itkConceptMacro(OutputEqualityComparableCheck,
00333 (Concept::EqualityComparable<PixelType>));
00334 itkConceptMacro(DoubleConvertibleToOutputCheck,
00335 (Concept::Convertible<double, PixelType>));
00336 itkConceptMacro(OutputOStreamWritableCheck,
00337 (Concept::OStreamWritable<PixelType>));
00338
00340 #endif
00341
00342 protected:
00343 ParallelSparseFieldLevelSetImageFilter();
00344 ~ParallelSparseFieldLevelSetImageFilter() {}
00345 virtual void PrintSelf(std::ostream& os, Indent indent) const;
00346
00348 ParallelSparseFieldCityBlockNeighborList < NeighborhoodIterator<OutputImageType> >
00349 m_NeighborList;
00350
00353 double m_ConstantGradientValue;
00354
00356 static ValueType m_ValueOne;
00357
00359 static ValueType m_ValueZero;
00360
00363 static StatusType m_StatusActiveChangingUp;
00364
00367 static StatusType m_StatusActiveChangingDown;
00368
00371 static StatusType m_StatusNull;
00372
00375 static StatusType m_StatusChanging;
00376
00379 static StatusType m_StatusBoundaryPixel;
00380
00385 typename OutputImageType::Pointer m_ShiftedImage;
00386
00391 LayerListType m_Layers;
00392
00397 StatusType m_NumberOfLayers;
00398
00400 typename StatusImageType::Pointer m_StatusImage;
00401 typename OutputImageType::Pointer m_OutputImage;
00402
00404 typename StatusImageType::Pointer m_StatusImageTemp;
00405 typename OutputImageType::Pointer m_OutputImageTemp;
00406
00408 typename LayerNodeStorageType::Pointer m_LayerNodeStore;
00409
00411 ValueType m_IsoSurfaceValue;
00412
00416
00417
00420 virtual void GenerateData();
00421
00426 void CopyInputToOutput();
00427
00429 void AllocateUpdateBuffer() {}
00430
00433 void Initialize();
00434
00439 void ConstructActiveLayer();
00440
00442 void InitializeActiveLayerValues();
00443
00447 void ConstructLayer(StatusType from, StatusType to);
00448
00450 void ProcessStatusList(LayerType *InputList, StatusType ChangeToStatus,
00451 StatusType SearchForStatus, unsigned int ThreadId);
00452
00457 void PropagateAllLayerValues();
00458
00466 void PropagateLayerValues(StatusType from, StatusType to, StatusType promote,
00467 unsigned int InOrOut);
00468
00473 virtual void InitializeBackgroundPixels();
00474
00478 void ThreadedAllocateData(unsigned int ThreadId);
00479 void ThreadedInitializeData(unsigned int ThreadId, const ThreadRegionType & ThreadRegion);
00481
00491 void ComputeInitialThreadBoundaries();
00492
00494 unsigned int GetThreadNumber(unsigned int splitAxisValue);
00495
00497 void GetThreadRegionSplitByBoundary(unsigned int ThreadId, ThreadRegionType& ThreadRegion);
00498
00501 void DeallocateData();
00502
00504 struct ParallelSparseFieldLevelSetThreadStruct
00505 {
00506 ParallelSparseFieldLevelSetImageFilter* Filter;
00507 TimeStepType* TimeStepList;
00508 bool* ValidTimeStepList;
00509 TimeStepType TimeStep;
00510 };
00511
00516 void Iterate();
00517 static ITK_THREAD_RETURN_TYPE IterateThreaderCallback(void * arg);
00519
00524 inline virtual ValueType ThreadedCalculateUpdateValue(const unsigned int itkNotUsed(ThreadId),
00525 const IndexType itkNotUsed(index),
00526 const TimeStepType &dt,
00527 const ValueType &value,
00528 const ValueType &change)
00529 {
00530 return (value + dt * change);
00531 }
00532
00533
00534
00535
00536
00537 virtual void ThreadedProcessPixelEnteringActiveLayer (const IndexType itkNotUsed(index),
00538 const ValueType itkNotUsed(value),
00539 const unsigned int itkNotUsed(ThreadId));
00540
00542 void ApplyUpdate(TimeStepType) {}
00543
00546 virtual void ThreadedApplyUpdate(TimeStepType dt, unsigned int ThreadId);
00547
00549 TimeStepType CalculateChange()
00550 {
00551 return NumericTraits<TimeStepType>::Zero;
00552 }
00553
00556 virtual TimeStepType ThreadedCalculateChange(unsigned int ThreadId);
00557
00564 void ThreadedUpdateActiveLayerValues(TimeStepType dt, LayerType *StatusUpList,
00565 LayerType *StatusDownList, unsigned int ThreadId);
00566
00568 void CopyInsertList(unsigned int ThreadId, LayerPointerType FromListPtr,
00569 LayerPointerType ToListPtr);
00570
00572 void ClearList(unsigned int ThreadId, LayerPointerType ListPtr);
00573
00576 void CopyInsertInterNeighborNodeTransferBufferLayers(unsigned int ThreadId,
00577 LayerPointerType InputList,
00578 unsigned int InOrOut,
00579 unsigned int BufferLayerNumber);
00580
00583 void ClearInterNeighborNodeTransferBufferLayers(unsigned int ThreadId, unsigned int InOrOut,
00584 unsigned int BufferLayerNumber);
00585
00593 void ThreadedProcessFirstLayerStatusLists(unsigned int InputLayerNumber,
00594 unsigned int OutputLayerNumber,
00595 StatusType SearchForStatus,
00596 unsigned int InOrOut,
00597 unsigned int BufferLayerNumber, unsigned int ThreadId);
00598
00604 void ThreadedProcessStatusList(unsigned int InputLayerNumber, unsigned int OutputLayerNumber,
00605 StatusType ChangeToStatus, StatusType SearchForStatus,
00606 unsigned int InOrOut,
00607 unsigned int BufferLayerNumber, unsigned int ThreadId);
00608
00612 void ThreadedProcessOutsideList(unsigned int InputLayerNumber, StatusType ChangeToStatus,
00613 unsigned int InOrOut, unsigned int BufferLayerNumber, unsigned int ThreadId);
00614
00616 void ThreadedPropagateLayerValues(StatusType from, StatusType to, StatusType promote,
00617 unsigned int InorOut, unsigned int ThreadId);
00618
00621 void GetThreadRegionSplitUniformly(unsigned int ThreadId, ThreadRegionType& ThreadRegion);
00622
00628 void ThreadedPostProcessOutput(const ThreadRegionType & regionToProcess);
00629
00640 virtual void CheckLoadBalance();
00641
00644 virtual void ThreadedLoadBalance(unsigned int ThreadId);
00645
00647 void WaitForAll();
00648 void SignalNeighborsAndWait (unsigned int ThreadId);
00649 void SignalNeighbor (unsigned int SemaphoreArrayNumber, unsigned int ThreadId);
00650 void WaitForNeighbor (unsigned int SemaphoreArrayNumber, unsigned int ThreadId);
00652
00655 virtual void ThreadedInitializeIteration (unsigned int ThreadId);
00656
00659
00660
00662 unsigned int m_NumOfThreads;
00663
00665 unsigned int m_SplitAxis;
00666
00668 unsigned int m_ZSize;
00669
00672 bool m_BoundaryChanged;
00673
00675 unsigned int * m_Boundary;
00676
00678 int * m_GlobalZHistogram;
00679
00681 unsigned int * m_MapZToThreadNumber;
00682
00685 int * m_ZCumulativeFrequency;
00686
00688 typename Barrier::Pointer m_Barrier;
00689
00691 struct ThreadData
00692 {
00693 char pad1 [128];
00694
00695 TimeStepType TimeStep;
00696 ThreadRegionType ThreadRegion;
00697 ValueType m_RMSChange;
00698 unsigned int m_Count;
00699
00701 LayerListType m_Layers;
00702
00704 LayerListType * m_LoadTransferBufferLayers;
00705
00707 typename LayerNodeStorageType::Pointer m_LayerNodeStore;
00708
00709 LayerPointerType UpList[2];
00710 LayerPointerType DownList[2];
00711
00714 LayerPointerType** m_InterNeighborNodeTransferBufferLayers[2];
00715
00718 void * globalData;
00719
00721 int * m_ZHistogram;
00722
00726 typename Semaphore::Pointer m_Semaphore[2];
00727
00729 unsigned int m_SemaphoreArrayNumber;
00730
00731 char m_Pad2 [128];
00732 };
00733
00735 ThreadData *m_Data;
00736
00739 bool m_Stop;
00740
00745 bool m_InterpolateSurfaceLocation;
00746
00747 private:
00748
00749 ParallelSparseFieldLevelSetImageFilter(const Self&);
00750 void operator=(const Self&);
00751
00754 bool m_BoundsCheckingActive;
00755 };
00756
00757 }
00758
00759 #ifndef ITK_MANUAL_INSTANTIATION
00760 #include "itkParallelSparseFieldLevelSetImageFilter.txx"
00761 #endif
00762
00763 #endif
00764