ITK  4.1.0
Insight Segmentation and Registration Toolkit
itkParallelSparseFieldLevelSetImageFilter.h
Go to the documentation of this file.
00001 /*=========================================================================
00002  *
00003  *  Copyright Insight Software Consortium
00004  *
00005  *  Licensed under the Apache License, Version 2.0 (the "License");
00006  *  you may not use this file except in compliance with the License.
00007  *  You may obtain a copy of the License at
00008  *
00009  *         http://www.apache.org/licenses/LICENSE-2.0.txt
00010  *
00011  *  Unless required by applicable law or agreed to in writing, software
00012  *  distributed under the License is distributed on an "AS IS" BASIS,
00013  *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
00014  *  See the License for the specific language governing permissions and
00015  *  limitations under the License.
00016  *
00017  *=========================================================================*/
00018 #ifndef __itkParallelSparseFieldLevelSetImageFilter_h
00019 #define __itkParallelSparseFieldLevelSetImageFilter_h
00020 
00021 #include <vector>
00022 #include "itkFiniteDifferenceImageFilter.h"
00023 #include "itkSparseFieldLayer.h"
00024 #include "itkObjectStore.h"
00025 #include "itkNeighborhoodIterator.h"
00026 #include "itkMultiThreader.h"
00027 #include "itkBarrier.h"
00028 
00029 namespace itk
00030 {
00036 template< class TNodeIndexType >
00037 class ParallelSparseFieldLevelSetNode
00038 {
00039 public:
00040   TNodeIndexType                   m_Index;
00041   float                            m_Value;
00042   ParallelSparseFieldLevelSetNode *Next;
00043   ParallelSparseFieldLevelSetNode *Previous;
00044 };
00045 
00074 template< class TNeighborhoodType >
00075 class ParallelSparseFieldCityBlockNeighborList
00076 {
00077 public:
00078   typedef TNeighborhoodType                     NeighborhoodType;
00079   typedef typename NeighborhoodType::OffsetType OffsetType;
00080   typedef typename NeighborhoodType::RadiusType RadiusType;
00081 
00082   itkStaticConstMacro(Dimension, unsigned int, NeighborhoodType::Dimension);
00083 
00084   const RadiusType & GetRadius() const
00085   {
00086     return m_Radius;
00087   }
00088 
00089   const unsigned int & GetArrayIndex(unsigned int i) const
00090   {
00091     return m_ArrayIndex[i];
00092   }
00093 
00094   const OffsetType & GetNeighborhoodOffset(unsigned int i) const
00095   {
00096     return m_NeighborhoodOffset[i];
00097   }
00098 
00099   const unsigned int & GetSize() const
00100   {
00101     return m_Size;
00102   }
00103 
00104   unsigned int GetStride(unsigned int i)
00105   {
00106     return m_StrideTable[i];
00107   }
00108 
00109   ParallelSparseFieldCityBlockNeighborList();
00110 
00111   ~ParallelSparseFieldCityBlockNeighborList()
00112   {
00113     m_ArrayIndex.clear();
00114     m_NeighborhoodOffset.clear();
00115   }
00116 
00117   void Print(std::ostream & os) const;
00118 
00119 private:
00120   char                        m_Pad1[128];
00121   unsigned int                m_Size;
00122   RadiusType                  m_Radius;
00123   std::vector< unsigned int > m_ArrayIndex;
00124   std::vector< OffsetType >   m_NeighborhoodOffset;
00125 
00128   unsigned int m_StrideTable[Dimension];
00129   char         m_Pad2[128];
00130 };
00131 
00248 template< class TInputImage, class TOutputImage >
00249 class ITK_EXPORT ParallelSparseFieldLevelSetImageFilter:
00250   public FiniteDifferenceImageFilter< TInputImage, TOutputImage >
00251 {
00252 public:
00253 
00255   typedef ParallelSparseFieldLevelSetImageFilter                   Self;
00256   typedef FiniteDifferenceImageFilter< TInputImage, TOutputImage > Superclass;
00257   typedef SmartPointer< Self >                                     Pointer;
00258   typedef SmartPointer< const Self >                               ConstPointer;
00259 
00261   typedef typename Superclass::TimeStepType                 TimeStepType;
00262   typedef typename Superclass::FiniteDifferenceFunctionType FiniteDifferenceFunctionType;
00263   typedef typename Superclass::RadiusType                   RadiusType;
00264   typedef typename Superclass::NeighborhoodScalesType       NeighborhoodScalesType;
00265 
00267   itkNewMacro(Self);
00268 
00270   itkTypeMacro(ParallelSparseFieldLevelSetImageFilter, FiniteDifferenceImageFilter);
00271 
00273   typedef TInputImage                         InputImageType;
00274   typedef TOutputImage                        OutputImageType;
00275   typedef typename OutputImageType::IndexType IndexType;
00276 
00277   itkStaticConstMacro(ImageDimension, unsigned int, TOutputImage::ImageDimension);
00278 
00279   typedef typename OutputImageType::PixelType PixelType;
00280 
00281   typedef typename OutputImageType::RegionType ThreadRegionType;
00282 
00285   typedef typename OutputImageType::ValueType ValueType;
00286 
00288   typedef ParallelSparseFieldLevelSetNode< IndexType > LayerNodeType;
00289 
00291   typedef SparseFieldLayer< LayerNodeType > LayerType;
00292   typedef typename LayerType::Pointer       LayerPointerType;
00293 
00295   typedef std::vector< LayerPointerType > LayerListType;
00296 
00298   typedef signed char StatusType;
00299 
00302   typedef Image< StatusType, itkGetStaticConstMacro(ImageDimension) > StatusImageType;
00303 
00306   typedef ObjectStore< LayerNodeType > LayerNodeStorageType;
00307 
00308   typedef Offset< itkGetStaticConstMacro(ImageDimension) > OffsetType;
00309 
00313   itkSetMacro(NumberOfLayers, StatusType);
00314   itkGetConstMacro(NumberOfLayers, StatusType);
00316 
00318   itkSetMacro(IsoSurfaceValue, ValueType);
00319   itkGetConstMacro(IsoSurfaceValue, ValueType);
00321 
00322   LayerPointerType GetActiveListForIndex(const IndexType index)
00323   {
00324     // get the 'z' value for the index
00325     const unsigned int indexZ = index[m_SplitAxis];
00326     // get the thread in whose region the index lies
00327     const unsigned int ThreadNum = this->GetThreadNumber (indexZ);
00328 
00329     // get the active list for that thread
00330     return m_Data[ThreadNum].m_Layers[0];
00331   }
00332 
00333 #ifdef ITK_USE_CONCEPT_CHECKING
00334 
00335   itkConceptMacro( OutputEqualityComparableCheck,
00336                    ( Concept::EqualityComparable< PixelType > ) );
00337   itkConceptMacro( DoubleConvertibleToOutputCheck,
00338                    ( Concept::Convertible< double, PixelType > ) );
00339   itkConceptMacro( OutputOStreamWritableCheck,
00340                    ( Concept::OStreamWritable< PixelType > ) );
00341 
00343 #endif
00344 protected:
00345   ParallelSparseFieldLevelSetImageFilter();
00346   ~ParallelSparseFieldLevelSetImageFilter() {}
00347   virtual void PrintSelf(std::ostream & os, Indent indent) const;
00349 
00351   ParallelSparseFieldCityBlockNeighborList< NeighborhoodIterator< OutputImageType > >
00352   m_NeighborList;
00353 
00356   double m_ConstantGradientValue;
00357 
00359   static ValueType m_ValueOne;
00360 
00362   static ValueType m_ValueZero;
00363 
00366   static StatusType m_StatusActiveChangingUp;
00367 
00370   static StatusType m_StatusActiveChangingDown;
00371 
00374   static StatusType m_StatusNull;
00375 
00378   static StatusType m_StatusChanging;
00379 
00382   static StatusType m_StatusBoundaryPixel;
00383 
00388   typename OutputImageType::Pointer m_ShiftedImage;
00389 
00394   LayerListType m_Layers;
00395 
00400   StatusType m_NumberOfLayers;
00401 
00403   typename StatusImageType::Pointer m_StatusImage;
00404   typename OutputImageType::Pointer m_OutputImage;
00405 
00408   typename StatusImageType::Pointer m_StatusImageTemp;
00409   typename OutputImageType::Pointer m_OutputImageTemp;
00410 
00412   typename LayerNodeStorageType::Pointer m_LayerNodeStore;
00413 
00415   ValueType m_IsoSurfaceValue;
00416 
00420   //  ValueType m_RMSChange;
00421 
00424   virtual void GenerateData();
00425 
00430   void CopyInputToOutput();
00431 
00433   void AllocateUpdateBuffer() {}
00434 
00437   void Initialize();
00438 
00443   void ConstructActiveLayer();
00444 
00446   void InitializeActiveLayerValues();
00447 
00451   void ConstructLayer(const StatusType& from, const StatusType& to);
00452 
00454   void ProcessStatusList(LayerType *InputList, const StatusType& ChangeToStatus,
00455                          const StatusType& SearchForStatus, ThreadIdType ThreadId);
00456 
00461   void PropagateAllLayerValues();
00462 
00470   void PropagateLayerValues(const StatusType& from, const StatusType& to,
00471                             const StatusType& promote, unsigned int InOrOut);
00472 
00477   virtual void InitializeBackgroundPixels();
00478 
00482   void ThreadedAllocateData(ThreadIdType ThreadId);
00483 
00484   void ThreadedInitializeData(ThreadIdType ThreadId, const ThreadRegionType & ThreadRegion);
00485 
00495   void ComputeInitialThreadBoundaries();
00496 
00498   unsigned int GetThreadNumber(unsigned int splitAxisValue);
00499 
00501   void GetThreadRegionSplitByBoundary(ThreadIdType ThreadId, ThreadRegionType & ThreadRegion);
00502 
00505   void DeallocateData();
00506 
00508   struct ParallelSparseFieldLevelSetThreadStruct {
00509     ParallelSparseFieldLevelSetImageFilter *Filter;
00510     std::vector< TimeStepType > TimeStepList;
00511     std::vector< bool > ValidTimeStepList;
00512     TimeStepType TimeStep;
00513   };
00514 
00519   void Iterate();
00520 
00521   static ITK_THREAD_RETURN_TYPE IterateThreaderCallback(void *arg);
00522 
00527   inline virtual ValueType ThreadedCalculateUpdateValue(const ThreadIdType itkNotUsed(ThreadId),
00528                                                         const IndexType itkNotUsed(index),
00529                                                         const TimeStepType & dt,
00530                                                         const ValueType & value,
00531                                                         const ValueType & change)
00532   {
00533     return ( value + static_cast< ValueType >( dt ) * change );
00534   }
00535 
00536   // This method can be overridden in derived classes.
00537   // The pixel at 'index' is entering the active layer for thread 'ThreadId'.
00538   // The outputimage at 'index' will have the value as given by the
00539   // 'value' parameter.
00540   virtual void ThreadedProcessPixelEnteringActiveLayer( const IndexType& itkNotUsed(index),
00541                                                         const ValueType& itkNotUsed(value),
00542                                                         ThreadIdType itkNotUsed(ThreadId) );
00543 
00545   void ApplyUpdate(const TimeStepType&) {}
00546 
00549   virtual void ThreadedApplyUpdate(const TimeStepType& dt, ThreadIdType ThreadId);
00550 
00552   TimeStepType CalculateChange()
00553   {
00554     return NumericTraits< TimeStepType >::Zero;
00555   }
00556 
00559   virtual TimeStepType ThreadedCalculateChange(ThreadIdType ThreadId);
00560 
00567   void ThreadedUpdateActiveLayerValues( const TimeStepType& dt,
00568     LayerType   *StatusUpList,
00569     LayerType   *StatusDownList,
00570     ThreadIdType ThreadId);
00571 
00574   void CopyInsertList( ThreadIdType ThreadId,
00575     LayerPointerType FromListPtr,
00576     LayerPointerType ToListPtr);
00577 
00579   void ClearList(ThreadIdType ThreadId, LayerPointerType ListPtr);
00580 
00583   void CopyInsertInterNeighborNodeTransferBufferLayers(
00584     ThreadIdType ThreadId,
00585     LayerPointerType InputList,
00586     unsigned int InOrOut,
00587     unsigned int BufferLayerNumber);
00588 
00591   void ClearInterNeighborNodeTransferBufferLayers(
00592     ThreadIdType ThreadId, unsigned int InOrOut,
00593     unsigned int BufferLayerNumber);
00594 
00602   void ThreadedProcessFirstLayerStatusLists(
00603     unsigned int InputLayerNumber,
00604     unsigned int OutputLayerNumber,
00605     const StatusType& SearchForStatus,
00606     unsigned int InOrOut,
00607     unsigned int BufferLayerNumber,
00608     ThreadIdType ThreadId);
00609 
00615   void ThreadedProcessStatusList(
00616     unsigned int InputLayerNumber,
00617     unsigned int OutputLayerNumber,
00618     const StatusType& ChangeToStatus,
00619     const StatusType& SearchForStatus,
00620     unsigned int InOrOut,
00621     unsigned int BufferLayerNumber,
00622     ThreadIdType ThreadId);
00623 
00627   void ThreadedProcessOutsideList(
00628     unsigned int InputLayerNumber,
00629     const StatusType& ChangeToStatus,
00630     unsigned int InOrOut,
00631     unsigned int BufferLayerNumber,
00632     ThreadIdType ThreadId);
00633 
00635   void ThreadedPropagateLayerValues(
00636     const StatusType& from,
00637     const StatusType& to,
00638     const StatusType& promote,
00639     unsigned int InorOut,
00640     ThreadIdType ThreadId);
00641 
00644   void GetThreadRegionSplitUniformly(
00645     ThreadIdType ThreadId, ThreadRegionType & ThreadRegion);
00646 
00652   void ThreadedPostProcessOutput(const ThreadRegionType & regionToProcess);
00653 
00664   virtual void CheckLoadBalance();
00665 
00668   virtual void ThreadedLoadBalance(ThreadIdType ThreadId);
00669 
00671   void WaitForAll();
00672 
00673   void SignalNeighborsAndWait(ThreadIdType ThreadId);
00674 
00675   void SignalNeighbor(unsigned int SemaphoreArrayNumber, ThreadIdType ThreadId);
00676 
00677   void WaitForNeighbor(unsigned int SemaphoreArrayNumber, ThreadIdType ThreadId);
00678 
00681   virtual void ThreadedInitializeIteration(ThreadIdType ThreadId);
00682 
00685   //  void WriteActivePointsToFile ();
00686 
00688   ThreadIdType m_NumOfThreads;
00689 
00691   unsigned int m_SplitAxis;
00692 
00694   unsigned int m_ZSize;
00695 
00698   bool m_BoundaryChanged;
00699 
00701   unsigned int *m_Boundary;
00702 
00704   int *m_GlobalZHistogram;
00705 
00708   unsigned int *m_MapZToThreadNumber;
00709 
00712   int *m_ZCumulativeFrequency;
00713 
00715   typename Barrier::Pointer m_Barrier;
00716 
00718   struct ThreadData {
00719     char pad1[128];
00720 
00721     TimeStepType TimeStep;
00722     ThreadRegionType ThreadRegion;
00723     ValueType m_RMSChange;
00724     unsigned int m_Count;
00725 
00727     LayerListType m_Layers;
00728 
00730     LayerListType *m_LoadTransferBufferLayers;
00731 
00733     typename LayerNodeStorageType::Pointer m_LayerNodeStore;
00734 
00735     LayerPointerType UpList[2];
00736     LayerPointerType DownList[2];
00737 
00740     LayerPointerType **m_InterNeighborNodeTransferBufferLayers[2];
00741 
00744     void *globalData;
00745 
00747     int *m_ZHistogram;
00748 
00753     int m_Semaphore[2];
00754 
00755     SimpleMutexLock            m_Lock[2];
00756     ConditionVariable::Pointer m_Condition[2];
00757 
00760     unsigned int m_SemaphoreArrayNumber;
00761 
00762     char m_Pad2[128];
00763   };
00764 
00767   ThreadData *m_Data;
00768 
00771   bool m_Stop;
00772 
00777   bool m_InterpolateSurfaceLocation;
00778 private:
00779 
00780   ParallelSparseFieldLevelSetImageFilter(const Self &); // purposely not
00781                                                         // implemented
00782   void operator=(const Self &);                         // purposely not
00783 
00784   // implemented
00785 
00788   bool m_BoundsCheckingActive;
00789 };
00790 } // end namespace itk
00791 
00792 #ifndef ITK_MANUAL_INSTANTIATION
00793 #include "itkParallelSparseFieldLevelSetImageFilter.hxx"
00794 #endif
00795 
00796 #endif
00797