|
GEL
2
GEL is a library for Geometry and Linear Algebra
|
/Users/jab/Documents/Teaching/02585/GEL2_and_demos/GEL/src/Geometry/GridAlgorithm.h
00001 #ifndef __GEOMETRY_GRIDALGORITHM_H 00002 #define __GEOMETRY_GRIDALGORITHM_H 00003 00004 /* 00005 Functions: 00006 00007 void for_each_voxel(Grid& g, F& f); 00008 void for_each_voxel(Grid& g, const Vec3i& p0, const Vec3i& p7, F& f); 00009 void for_each_voxel_ordered(Grid& g, const Vec3i& p0, const Vec3i& p7, F& f); 00010 void for_each_voxel_ordered(Grid& g, F&); 00011 void for_each_voxel_const(const Grid& g, 00012 const Vec3i& p0, const Vec3i& p7, F& f); 00013 void for_each_voxel_const(const Grid& g, F& f); 00014 void for_each_voxel_ordered_const(Grid& g, 00015 const Vec3i& p0, const Vec3i& p7, 00016 F& f); 00017 void for_each_voxel_ordered_const(Grid& g, F& f); 00018 00019 Purpose: 00020 ---------- 00021 00022 Visit all voxels (or voxels in a region og interest) in grid. A function 00023 is invoked on each voxel. 00024 00025 The "ordered" functions traverse the grid in a systematic way: A 00026 slice at a time and for each slice one row at a time. The other 00027 functions make no guarantees about how the volume (roi) is traversed. 00028 00029 Types: 00030 ---------- 00031 00032 Grid - a grid type, either HGrid<T,CellT> or RGrid<T> 00033 00034 F - functor type a funtion or class with the function call operator 00035 overloaded. The functor should look either like this 00036 00037 void fun(const Vec3i&, T& x) 00038 00039 or this 00040 00041 void fun(const Vec3i&, const T& x) 00042 00043 in the case of the const functions. 00044 00045 Arguments: 00046 ---------- 00047 00048 g - The voxel grid, we wish to traverse. 00049 p0 - (xmin, ymin, zmin) coordinates of the window we wish to 00050 traverse. 00051 p7 - (xmax, ymax, zmax) coordinates of the window we wish to 00052 traverse. 00053 f - functor. 00054 00055 */ 00056 00057 #include <iostream> 00058 #include "RGrid.h" 00059 #include "HGrid.h" 00060 00061 namespace Geometry 00062 { 00063 template<class T, class F> 00064 void _for_each_voxel(T* data, 00065 int x_dim, int xy_dim, 00066 const CGLA::Vec3i& p0, 00067 const CGLA::Vec3i& p7, 00068 F& functor, 00069 const CGLA::Vec3i& offset = CGLA::Vec3i(0)) 00070 { 00071 const int Amin = p0[2]*xy_dim; 00072 const int Amax = p7[2]*xy_dim; 00073 int Bmin = Amin +p0[1]*x_dim; 00074 int Bmax = Amin +p7[1]*x_dim; 00075 CGLA::Vec3i p0o = p0+offset; 00076 CGLA::Vec3i p(p0o); 00077 for(int A=Amin; A<Amax; A+=xy_dim, ++p[2]) 00078 { 00079 p[1] = p0o[1]; 00080 for(int B = Bmin; B<Bmax; B+=x_dim, ++p[1]) 00081 { 00082 p[0] = p0o[0]; 00083 int Cmin = B+p0[0]; 00084 int Cmax = B+p7[0]; 00085 for(int C=Cmin; C<Cmax; ++C, ++p[0]) 00086 functor(p, data[C]); 00087 } 00088 Bmin += xy_dim; 00089 Bmax += xy_dim; 00090 } 00091 } 00092 00093 template<class T, class F> 00094 void _for_each_voxel(T* data, 00095 const CGLA::Vec3i& dims, 00096 F& functor, 00097 const CGLA::Vec3i& offset = CGLA::Vec3i(0)) 00098 { 00099 int l=0; 00100 CGLA::Vec3i p(offset); 00101 CGLA::Vec3i p7(offset); 00102 p7 += dims; 00103 00104 for(; p[2]<p7[2]; ++p[2]) 00105 for(p[1]=offset[1]; p[1]<p7[1]; ++p[1]) 00106 for(p[0]=offset[0]; p[0]<p7[0]; ++p[0]) 00107 functor(p, data[l++]); 00108 } 00109 00110 00116 template<class T, class F> 00117 void for_each_voxel(RGrid<T>& grid, 00118 const CGLA::Vec3i& p0, 00119 const CGLA::Vec3i& p7, 00120 F& functor) 00121 { 00122 _for_each_voxel(grid.get(), grid.get_x_dim(), grid.get_xy_dim(), 00123 CGLA::v_max(p0, CGLA::Vec3i(0)), 00124 CGLA::v_min(p7, grid.get_dims()), 00125 functor); 00126 } 00127 00132 template<class T, class F> 00133 void for_each_voxel(RGrid<T>& grid, F& functor) 00134 { 00135 _for_each_voxel(grid.get(), grid.get_dims(), functor); 00136 } 00137 00138 00146 template<class T, class F> 00147 void for_each_voxel_ordered(RGrid<T>& grid, 00148 const CGLA::Vec3i& p0, 00149 const CGLA::Vec3i& p7, 00150 F& functor) 00151 { 00152 _for_each_voxel(grid.get(), grid.get_x_dim(), grid.get_xy_dim(), 00153 CGLA::v_max(p0, CGLA::Vec3i(0)), 00154 CGLA::v_min(p7, grid.get_dims()), 00155 functor); 00156 } 00157 00158 template<class T, class F> 00159 void for_each_voxel_ordered(RGrid<T>& grid, F& functor) 00160 { 00161 _for_each_voxel(grid.get(), grid.get_dims(), functor); 00162 } 00163 00164 00165 template<class T, class CellT, class F> 00166 void for_each_cell(HGrid<T,CellT>& grid, 00167 const CGLA::Vec3i& p0, 00168 const CGLA::Vec3i& p7, 00169 F& functor) 00170 { 00171 CGLA::Vec3i p0t = p0/grid.get_bottom_dim(); 00172 CGLA::Vec3i p7t = CGLA::v_min(p7/grid.get_bottom_dim()+ 00173 CGLA::Vec3i(1), 00174 grid.get_top_dims()); 00175 for(CGLA::Vec3i pt(p0t); pt[2]<p7t[2]; ++pt[2]) 00176 for(pt[1]=p0t[1]; pt[1]<p7t[1]; ++pt[1]) 00177 for(pt[0]=p0t[0]; pt[0]<p7t[0]; ++pt[0]) 00178 functor(pt*CellT::get_dim(), grid.get_cell(pt)); 00179 } 00180 00181 template<class T, class CellT, class F> 00182 void for_each_cell(HGrid<T,CellT>& grid, 00183 F& functor) 00184 { 00185 CGLA::Vec3i p0t; 00186 CGLA::Vec3i p7t = grid.get_dims(); 00187 const int inc = CellT::get_dim(); 00188 int l=0; 00189 for(CGLA::Vec3i pt(p0t); pt[2]<p7t[2]; pt[2]+=inc) 00190 for(pt[1]=p0t[1]; pt[1]<p7t[1]; pt[1]+=inc) 00191 for(pt[0]=p0t[0]; pt[0]<p7t[0]; pt[0]+=inc) 00192 functor(pt, grid.get_cell(l++)); 00193 } 00194 00195 00196 template<class CellT, class F> 00197 class _HGridCellFunctor 00198 { 00199 const CGLA::Vec3i p0; 00200 const CGLA::Vec3i p7; 00201 F& functor; 00202 00203 public: 00204 _HGridCellFunctor(const CGLA::Vec3i _p0, 00205 const CGLA::Vec3i _p7, 00206 F& _functor): p0(_p0), p7(_p7), functor(_functor) {} 00207 00208 void operator()(const CGLA::Vec3i& offset, 00209 CellT& cell) 00210 { 00211 CGLA::Vec3i p0c = CGLA::v_max(p0-offset, CGLA::Vec3i(0)); 00212 CGLA::Vec3i p7c = CGLA::v_min(p7-offset, CGLA::Vec3i(CellT::get_dim())); 00213 00214 if(cell.is_coalesced()) 00215 cell.split(); 00216 00217 _for_each_voxel(cell.get(), 00218 CellT::get_dim(), 00219 CGLA::sqr(CellT::get_dim()), 00220 p0c, p7c, functor, offset); 00221 } 00222 }; 00223 00224 00225 template<class T, class CellT, class F> 00226 void for_each_voxel(HGrid<T,CellT>& grid, 00227 const CGLA::Vec3i& _p0, 00228 const CGLA::Vec3i& _p7, 00229 F& functor) 00230 { 00231 CGLA::Vec3i p0 = CGLA::v_max(_p0, CGLA::Vec3i(0)); 00232 CGLA::Vec3i p7 = CGLA::v_min(_p7, grid.get_dims()); 00233 _HGridCellFunctor<CellT,F> cell_functor(p0, p7, functor); 00234 for_each_cell(grid, p0, p7, cell_functor); 00235 } 00236 00237 template<class T, class CellT, class F> 00238 void for_each_voxel(HGrid<T,CellT>& grid, F& functor) 00239 { 00240 _HGridCellFunctor<CellT,F> cell_functor(CGLA::Vec3i(0), 00241 grid.get_dims(), functor); 00242 for_each_cell(grid, cell_functor); 00243 } 00244 00245 template<class T, class CellT, class F> 00246 void for_each_voxel_ordered(HGrid<T,CellT>& grid, 00247 const CGLA::Vec3i& _p0, 00248 const CGLA::Vec3i& _p7, 00249 F& functor) 00250 { 00251 CGLA::Vec3i p0 = CGLA::v_max(_p0, CGLA::Vec3i(0)); 00252 CGLA::Vec3i p7 = CGLA::v_min(_p7, grid.get_dims()); 00253 for(int k=p0[2];k<p7[2];++k) 00254 for(int j=p0[1];j<p7[1];++j) 00255 for(int i=p0[0];i<p7[0];++i) 00256 { 00257 CGLA::Vec3i p(i,j,k); 00258 float val = grid[p]; 00259 float nval = val; 00260 functor(p, nval); 00261 if(nval != val) 00262 grid.store(p, nval); 00263 } 00264 } 00265 00266 template<class T, class CellT, class F> 00267 void for_each_voxel_ordered(HGrid<T,CellT>& grid, F& functor) 00268 { 00269 for_each_voxel_ordered(grid, CGLA::Vec3i(0), grid.get_dims(), functor); 00270 } 00271 00272 00273 template<typename T> 00274 class _AssignFun 00275 { 00276 T val; 00277 public: 00278 _AssignFun(const T& _val): val(_val) {} 00279 void operator()(const CGLA::Vec3i& pi, T& vox_val) 00280 { 00281 vox_val = val; 00282 } 00283 }; 00284 00285 template<class G> 00286 void clear_region(G& grid, const typename G::DataType& value) 00287 { 00288 _AssignFun<typename G::DataType> afun(value); 00289 for_each_voxel(grid, afun); 00290 } 00291 00292 template<class G> 00293 void clear_region(G& grid, 00294 const CGLA::Vec3i& p0, 00295 const CGLA::Vec3i& p7, 00296 const typename G::DataType& value) 00297 { 00298 _AssignFun<typename G::DataType>afun(value) ; 00299 for_each_voxel(grid, p0, p7, afun); 00300 } 00301 00302 00303 //---------------------------------------------------------------------- 00304 // const versions. 00305 00311 template<class T, class F> 00312 void for_each_voxel_const(const RGrid<T>& grid, 00313 const CGLA::Vec3i& p0, 00314 const CGLA::Vec3i& p7, 00315 F& functor) 00316 { 00317 _for_each_voxel(grid.get(), grid.get_x_dim(), grid.get_xy_dim(), 00318 CGLA::v_max(p0, CGLA::Vec3i(0)), 00319 CGLA::v_min(p7, grid.get_dims()), 00320 functor); 00321 } 00322 00327 template<class T, class F> 00328 void for_each_voxel_const(const RGrid<T>& grid, F& functor) 00329 { 00330 _for_each_voxel(grid.get(), grid.get_dims(), functor); 00331 } 00332 00333 00341 template<class T, class F> 00342 void for_each_voxel_ordered_const(const RGrid<T>& grid, 00343 const CGLA::Vec3i& p0, 00344 const CGLA::Vec3i& p7, 00345 F& functor) 00346 { 00347 _for_each_voxel(grid.get(), grid.get_x_dim(), grid.get_xy_dim(), 00348 CGLA::v_max(p0, CGLA::Vec3i(0)), 00349 CGLA::v_min(p7, grid.get_dims()), 00350 functor); 00351 } 00352 00353 template<class T, class F> 00354 void for_each_voxel_ordered_const(const RGrid<T>& grid, F& functor) 00355 { 00356 _for_each_voxel(grid.get(), grid.get_dims(), functor); 00357 } 00358 00359 00360 template<class T, class CellT, class F> 00361 void for_each_cell_const(const HGrid<T,CellT>& grid, 00362 const CGLA::Vec3i& p0, 00363 const CGLA::Vec3i& p7, 00364 F& functor) 00365 { 00366 CGLA::Vec3i p0t = p0/grid.get_bottom_dim(); 00367 CGLA::Vec3i p7t = CGLA::v_min(p7/grid.get_bottom_dim()+ 00368 CGLA::Vec3i(1), 00369 grid.get_top_dims()); 00370 for(CGLA::Vec3i pt(p0t); pt[2]<p7t[2]; ++pt[2]) 00371 for(pt[1]=p0t[1]; pt[1]<p7t[1]; ++pt[1]) 00372 for(pt[0]=p0t[0]; pt[0]<p7t[0]; ++pt[0]) 00373 functor(pt*CellT::get_dim(), grid.get_cell(pt)); 00374 } 00375 00376 template<class T, class CellT, class F> 00377 void for_each_cell_const(const HGrid<T,CellT>& grid, 00378 F& functor) 00379 { 00380 CGLA::Vec3i p0t; 00381 CGLA::Vec3i p7t = grid.get_dims(); 00382 const int inc = CellT::get_dim(); 00383 int l=0; 00384 for(CGLA::Vec3i pt(p0t); pt[2]<p7t[2]; pt[2]+=inc) 00385 for(pt[1]=p0t[1]; pt[1]<p7t[1]; pt[1]+=inc) 00386 for(pt[0]=p0t[0]; pt[0]<p7t[0]; pt[0]+=inc) 00387 functor(pt, grid.get_cell(l++)); 00388 } 00389 00390 00391 template<class CellT, class F> 00392 class _HGridCellFunctorConst 00393 { 00394 const CGLA::Vec3i p0; 00395 const CGLA::Vec3i p7; 00396 F& functor; 00397 00398 public: 00399 _HGridCellFunctorConst(const CGLA::Vec3i _p0, 00400 const CGLA::Vec3i _p7, 00401 F& _functor): p0(_p0), p7(_p7), functor(_functor) {} 00402 00403 void operator()(const CGLA::Vec3i& offset, 00404 const CellT& cell) 00405 { 00406 CGLA::Vec3i p0c = CGLA::v_max(p0-offset, CGLA::Vec3i(0)); 00407 CGLA::Vec3i p7c = CGLA::v_min(p7-offset, CGLA::Vec3i(CellT::get_dim())); 00408 00409 if(cell.is_coalesced()) 00410 { 00411 typename CellT::DataType val = *cell.get(); 00412 for(CGLA::Vec3i p(p0c); p[2]<p7c[2]; ++p[2]) 00413 for(p[1]=p0c[1]; p[1]<p7c[1]; ++p[1]) 00414 for(p[0]=p0c[0]; p[0]<p7c[0]; ++p[0]) 00415 functor(p+offset, val); 00416 } 00417 _for_each_voxel(cell.get(), 00418 CellT::get_dim(), 00419 CGLA::sqr(CellT::get_dim()), 00420 p0c, p7c, functor, offset); 00421 } 00422 }; 00423 00424 00425 template<class T, class CellT, class F> 00426 void for_each_voxel_const(const HGrid<T,CellT>& grid, 00427 const CGLA::Vec3i& _p0, 00428 const CGLA::Vec3i& _p7, 00429 F& functor) 00430 { 00431 CGLA::Vec3i p0 = CGLA::v_max(_p0, CGLA::Vec3i(0)); 00432 CGLA::Vec3i p7 = CGLA::v_min(_p7, grid.get_dims()); 00433 _HGridCellFunctorConst<CellT,F> cell_functor(p0, p7, functor); 00434 for_each_cell_const(grid, p0, p7, cell_functor); 00435 } 00436 00437 template<class T, class CellT, class F> 00438 void for_each_voxel_const(const HGrid<T,CellT>& grid, F& functor) 00439 { 00440 _HGridCellFunctorConst<CellT,F> cell_functor(CGLA::Vec3i(0), 00441 grid.get_dims(), functor); 00442 for_each_cell_const(grid, cell_functor); 00443 } 00444 00445 template<class T, class CellT, class F> 00446 void for_each_voxel_ordered_const(const HGrid<T,CellT>& grid, 00447 const CGLA::Vec3i& _p0, 00448 const CGLA::Vec3i& _p7, 00449 F& functor) 00450 { 00451 CGLA::Vec3i p0 = CGLA::v_max(_p0, CGLA::Vec3i(0)); 00452 CGLA::Vec3i p7 = CGLA::v_min(_p7, grid.get_dims()); 00453 for(int k=p0[2];k<p7[2];++k) 00454 for(int j=p0[1];j<p7[1];++j) 00455 for(int i=p0[0];i<p7[0];++i) 00456 { 00457 CGLA::Vec3i p(i,j,k); 00458 functor(p, grid[p]); 00459 } 00460 } 00461 00462 template<class T, class CellT, class F> 00463 void for_each_voxel_ordered_const(const HGrid<T,CellT>& grid, F& functor) 00464 { 00465 for_each_voxel_ordered_const(grid, 00466 CGLA::Vec3i(0), 00467 grid.get_dims(), 00468 functor); 00469 } 00470 00471 00472 } 00473 00474 #endif
