Main Page   Namespace List   Class Hierarchy   Alphabetical List   Compound List   File List   Namespace Members   Compound Members   File Members   Related Pages  

Watershed.h

Go to the documentation of this file.
00001 /***************************************************************************
00002  *cr
00003  *cr            (C) Copyright 1995-2019 The Board of Trustees of the
00004  *cr                        University of Illinois
00005  *cr                         All Rights Reserved
00006  *cr
00007  ***************************************************************************/
00008 
00009 /***************************************************************************
00010  * RCS INFORMATION:
00011  *
00012  *      $RCSfile: Watershed.h,v $
00013  *      $Author: johns $        $Locker:  $             $State: Exp $
00014  *      $Revision: 1.17 $        $Date: 2019/01/17 21:21:02 $
00015  *
00016  ***************************************************************************
00017  * DESCRIPTION:
00018  *   CUDA-accelerated Watershed image segmentation
00019  ***************************************************************************/
00020 
00021 #ifndef WATERSHED_H
00022 #define WATERSHED_H
00023 
00024 
00025 #if defined(WATERSHED_INTERNAL) || defined(ARCH_SOLARISX86_64)
00026 
00027 /* Enables verbose printing throughout the watershed algorithm */
00028 
00029 /*
00030  * Uses 18-connected 3D neighbors instead of 6
00031  * They produce very similar but not identical results
00032  * It may acutally be slightly faster to use 18 neighbors
00033  * but uses more memory per voxel.
00034  * Using 18 neighbors results in 99.9% similarity between CUDA and
00035  * CPU version instead of 100. This is not an error in the CUDA implementation.
00036  * It occurs because currently we only mark 6-connected neighbors to update
00037  * if a voxel changes, if the updating scheme is turned off
00038  * then this issue goes away.
00039  */
00040 //#define CONNECTED_18
00041 
00042 /* 
00043  * Keeps track of blocks that need to be updated and
00044  * only performs calulations for those regions
00045  */
00046 #define BLOCK_UPDATES
00047 
00048 /* 
00049  * Allocates extra memory so that we can skip some bounds
00050  * checking when doing block update bookkeeping
00051  */
00052 #define FAST_UPDATES
00053 
00054 /* Tracks how often the voxels are actually updated */
00055 //#define STATS
00056 
00057 /* Times how long the algorithm runs */
00058 #define TIMER
00059 
00060 #define FLOAT_DIFF      0.0001f
00061 #define UPDATE_SIZE     4 // Update "blocks" size
00062 #define UPDATE_SHIFT    2 // log2 of UPDATE_SIZE
00063 
00064 #define EQ_SHIFT 6
00065 #define EQ_MASK  0x3F
00066 #define LOWER_MASK 0xFFFFFFC0
00067 
00068 #define GET_N_LOWER(equal_and_lower) ((equal_and_lower) & EQ_MASK)
00069 
00070 #define GET_N_EQ(equal_and_lower) ((equal_and_lower) >> EQ_SHIFT)
00071 
00072 #define MAKE_EQUAL_AND_LOWER(equal, lower) (((equal) << EQ_SHIFT) | (lower))
00073 
00074 #define REPLACE_LOWER(equal_and_lower, lower) (((equal_and_lower) & LOWER_MASK) | (lower))
00075 
00076 // offsets from starting voxel to each neighbor
00077 // px_ny represents the neighbor at (x+1, y-1, z) position where
00078 // x, y, and z are the coordinates of the starting voxel
00079 // Given a direction, such as px_ny, px_ny_offset gives the index offset
00080 // such that starting_index + px_ny_offset = the index of the px_ny neighbor.
00081 #define px_offset    (1)
00082 #define py_offset    (width)
00083 #define pz_offset    (heightWidth)
00084 #define nx_offset    (-1)
00085 #define ny_offset    (-width)
00086 #define nz_offset    (-heightWidth)
00087 
00088 #define nx_ny_offset (-1 - width)
00089 #define nx_py_offset (-1 + width)
00090 #define px_py_offset (1 + width)
00091 #define px_ny_offset (1 - width)
00092 
00093 #define px_pz_offset (1 + heightWidth)
00094 #define nx_nz_offset (-1 - heightWidth)
00095 #define nx_pz_offset (-1 + heightWidth)
00096 #define px_nz_offset (1 - heightWidth)
00097 
00098 #define py_pz_offset (width + heightWidth)
00099 #define ny_nz_offset (-width - heightWidth)
00100 #define ny_pz_offset (-width + heightWidth)
00101 #define py_nz_offset (width - heightWidth)
00102 
00103 #define SWAP(a, b, type) {\
00104   type t = a;\
00105   a = b;\
00106   b = t;\
00107 }
00108 
00109 
00113 enum direction {
00114   px = 0x00001,
00115   py = 0x00002,
00116   pz = 0x00004,
00117   nx = 0x00008,
00118   ny = 0x00010,
00119   nz = 0x00020,
00120 
00121   nx_ny = 0x00040,
00122   nx_py = 0x00080,
00123   px_py = 0x00100,
00124   px_ny = 0x00200,
00125 
00126   px_pz = 0x00400,
00127   nx_nz = 0x00800,
00128   nx_pz = 0x01000,
00129   px_nz = 0x02000,
00130 
00131   py_pz = 0x04000,
00132   ny_nz = 0x08000,
00133   ny_pz = 0x10000,
00134   py_nz = 0x20000
00135 };
00136 
00140 enum neighbor {
00141   px_idx = 1,
00142   py_idx,
00143   pz_idx,
00144   nx_idx,
00145   ny_idx,
00146   nz_idx,
00147 
00148   nx_ny_idx,
00149   nx_py_idx,
00150   px_py_idx,
00151   px_ny_idx,
00152 
00153   px_pz_idx,
00154   nx_nz_idx,
00155   nx_pz_idx,
00156   px_nz_idx,
00157 
00158   py_pz_idx,
00159   ny_nz_idx,
00160   ny_pz_idx,
00161   py_nz_idx
00162 };
00163 
00164 #endif // WATERSHED_INTERNAL
00165 
00166 
00167 template <typename GROUP_T, typename IMAGE_T>
00168 struct state_t {
00169   GROUP_T* group;
00170   IMAGE_T* value;
00171 };
00172 
00173 template <typename IMAGE_T>
00174 struct group_t {
00175   IMAGE_T* max_value;
00176   int* max_idx;
00177   bool init;
00178 };
00179 
00180 template<typename GROUP_T, typename IMAGE_T>
00181 struct watershed_gpu_state_t {
00182   int* eq_and_lower_d;
00183   IMAGE_T* current_value_d;
00184   IMAGE_T* next_value_d;
00185   unsigned char* current_update_d;
00186   unsigned char* next_update_d;
00187   GROUP_T* segments_d;
00188   int height;
00189   int width;
00190   int depth;
00191   bool init;
00192 };
00193 
00195 template<typename GROUP_T, typename IMAGE_T>
00196 class Watershed {
00197   public:
00198 
00200     Watershed<GROUP_T, IMAGE_T>(unsigned int h, unsigned int w, unsigned int d, bool cuda=true);
00201 
00204     void getSegmentedVoxels(IMAGE_T* voxels);
00205     
00206     IMAGE_T* getRawVoxels(); 
00207 
00208     ~Watershed<GROUP_T, IMAGE_T>();
00209 
00211     void watershed(IMAGE_T* image, int imageongpu, GROUP_T* segments, bool verbose=true);
00212 
00213   private:
00214     bool use_cuda; 
00215     watershed_gpu_state_t<GROUP_T, IMAGE_T> gpu_state;
00216 
00217     int* equal_and_lower; 
00218 
00219 
00220 
00221     unsigned char* current_update; 
00222     unsigned char* next_update;  
00223 
00224     int height;                  
00225     int width;                   
00226     int depth;                   
00227     int heightWidth;             
00228     long nVoxels;                   
00229 
00230     unsigned long update_width;  
00231     unsigned int update_offset;  
00232 
00233     group_t<IMAGE_T> group_data;          
00234     long nGroups;                 
00235     int neighbor_offsets[19];    
00236 
00237     state_t<GROUP_T, IMAGE_T> current_state;    
00238     state_t<GROUP_T, IMAGE_T> next_state;       
00239 
00240     void init(IMAGE_T* image);     
00241 
00242     void init_neighbor_offsets();
00243 
00244 
00252     void gaussian3D(const float sigma, const int reps);
00253 
00254 #if defined(VMDCUDA)
00255     void watershed_gpu(GROUP_T* segments_d); 
00256 #endif
00257 
00258     void watershed_cpu(GROUP_T* segments); 
00259 
00260     unsigned int update(); 
00261 
00263     unsigned int update_blocks();
00264 
00266     unsigned int update_stats(unsigned int& numUpdates);
00267 
00268 
00271     unsigned int update_blocks_stats(unsigned int& numUpdates,
00272                                      unsigned int& numBlockUpdates); 
00273 
00276     long find_local_maxima(long x, long y, long z);
00277 };
00278 
00279 
00280 #endif

Generated on Thu Mar 28 02:44:28 2024 for VMD (current) by doxygen1.2.14 written by Dimitri van Heesch, © 1997-2002