Subversion Repositories gelsvn

Rev

Go to most recent revision | Details | Last modification | View Log | RSS feed

Rev Author Line No. Line
61 jab 1
#ifndef __GRIDALGORITHM_H
2
#define __GRIDALGORITHM_H
3
 
4
/*
5
Functions:
6
 
7
  void for_each_voxel(Grid& g, F& f);
8
  void for_each_voxel(Grid& g,	const Vec3i& p0,	const Vec3i& p7, F& f);
9
  void for_each_voxel_ordered(Grid& g, const Vec3i& p0, const Vec3i& p7, F& f);
10
  void for_each_voxel_ordered(Grid& g,	F&);
11
  void for_each_voxel_const(const Grid& g, 
12
														const Vec3i& p0, const Vec3i& p7, F& f);
13
  void for_each_voxel_const(const Grid& g,	F& f);
14
  void for_each_voxel_ordered_const(Grid& g, 
15
																		const Vec3i& p0, const Vec3i& p7, 
16
																		F& f);
17
  void for_each_voxel_ordered_const(Grid& g, F& f);
18
 
19
Purpose:  
20
----------
21
 
22
Visit all voxels (or voxels in a region og interest) in grid. A function 
23
is invoked on each voxel. 
24
 
25
The "ordered" functions traverse the grid in a systematic way: A
26
slice at a time and for each slice one row at a time. The other
27
functions make no guarantees about how the volume (roi) is traversed. 
28
 
29
Types: 
30
----------
31
 
32
Grid - a grid type, either HGrid<T,CellT> or RGrid<T>
33
 
34
F - functor type a funtion or class with the function call operator 
35
overloaded. The functor should look either like this
36
 
37
void fun(const Vec3i&, T& x)
38
 
39
or this
40
 
41
void fun(const Vec3i&, const T& x)
42
 
43
in the case of the const functions.
44
 
45
Arguments:
46
----------
47
 
48
g  - The voxel grid, we wish to traverse.
49
p0 - (xmin, ymin, zmin) coordinates of the window we wish to
50
     traverse.
51
p7 - (xmax, ymax, zmax) coordinates of the window we wish to
52
     traverse.
53
f  - functor.
54
 
55
*/
56
 
57
#include <iostream>
58
#include "RGrid.h"
59
#include "HGrid.h"
60
 
61
namespace Geometry 
62
{
63
	template<class T, class F>
64
	void _for_each_voxel(T* data, 
65
											 int x_dim, int xy_dim,
66
											 const CGLA::Vec3i& p0, 
67
											 const CGLA::Vec3i& p7,
68
											 F& functor,
69
											 const CGLA::Vec3i& offset = CGLA::Vec3i(0))
70
	{
71
		const int Amin = p0[2]*xy_dim;
72
		const int Amax = p7[2]*xy_dim;
73
		int Bmin = Amin +p0[1]*x_dim;
74
		int Bmax = Amin +p7[1]*x_dim; 
75
		CGLA::Vec3i p0o = p0+offset;
76
		CGLA::Vec3i p(p0o);
77
		for(int A=Amin; A<Amax; A+=xy_dim, ++p[2])
78
			{
79
				p[1] = p0o[1];
80
				for(int B = Bmin; B<Bmax; B+=x_dim, ++p[1])
81
					{
82
						p[0] = p0o[0];
83
						int Cmin = B+p0[0];
84
						int Cmax = B+p7[0];
85
						for(int C=Cmin; C<Cmax; ++C, ++p[0])
86
							functor(p, data[C]);
87
					}
88
				Bmin += xy_dim;
89
				Bmax += xy_dim;
90
			}
91
	}
92
 
93
	template<class T, class F>
94
	void _for_each_voxel(T* data, 
95
											 const CGLA::Vec3i& dims,
96
											 F& functor,
97
											 const CGLA::Vec3i& offset = CGLA::Vec3i(0))
98
	{
99
		int l=0;
100
		CGLA::Vec3i p(offset);
101
		CGLA::Vec3i p7(offset);
102
		p7 += dims;
103
 
104
		for(;                   p[2]<p7[2]; ++p[2])
105
			for(p[1]=offset[1];   p[1]<p7[1]; ++p[1])
106
				for(p[0]=offset[0]; p[0]<p7[0]; ++p[0])
107
					functor(p, data[l++]);
108
	}
109
 
110
 
111
	/** Loop over all voxels in a sub-region (slice) of an
112
			RGrid and invoke a functor on each voxel. 
113
			The grid is the first argument, the slice is
114
			specified by the two subsequent args, and the functor 
115
			is the last argument. */
116
	template<class T, class F>
117
	void for_each_voxel(RGrid<T>& grid,
118
											const CGLA::Vec3i& p0, 
119
											const CGLA::Vec3i& p7,
120
											F& functor)
121
	{
122
		_for_each_voxel(grid.get(), grid.get_x_dim(), grid.get_xy_dim(),
123
										CGLA::v_max(p0, CGLA::Vec3i(0)),
124
										CGLA::v_min(p7, grid.get_dims()), 
125
										functor);
126
	}
127
 
128
	/** Loop over all voxels in an entire	RGrid. 
129
			Grid is the first argument, and a functor is 
130
			the second.	For each voxel, an operation 
131
			specified by the functor is performed. */
132
	template<class T, class F>
133
	void for_each_voxel(RGrid<T>& grid,	F& functor)
134
	{
135
		_for_each_voxel(grid.get(), grid.get_dims(), functor);
136
	}
137
 
138
 
139
	/** For each voxel (ordered). The idea of ordered traversal is 
140
			that we traverse the volume in a systematic fashion as opposed
141
			to traversing simply according to the memory layout of the volume
142
			data structure. This is important e.g. if we want to save the 
143
			volume in raw format. 
144
			For an RGrid, there is no difference though.
145
	*/ 
146
	template<class T, class F>
147
	void for_each_voxel_ordered(RGrid<T>& grid,
148
															const CGLA::Vec3i& p0, 
149
															const CGLA::Vec3i& p7,
150
															F& functor)
151
	{
152
		_for_each_voxel(grid.get(), grid.get_x_dim(), grid.get_xy_dim(),
153
										CGLA::v_max(p0, CGLA::Vec3i(0)),
154
										CGLA::v_min(p7, grid.get_dims()), 
155
										functor);
156
	}
157
 
158
	template<class T, class F>
159
	void for_each_voxel_ordered(RGrid<T>& grid,	F& functor)
160
	{
161
		_for_each_voxel(grid.get(), grid.get_dims(), functor);
162
	}
163
 
164
 
165
	template<class T, class CellT, class F>
166
	void for_each_cell(HGrid<T,CellT>& grid,
167
										 const CGLA::Vec3i& p0, 
168
										 const CGLA::Vec3i& p7,
169
										 F& functor)
170
	{
171
		CGLA::Vec3i p0t = p0/grid.get_bottom_dim();
172
		CGLA::Vec3i p7t = CGLA::v_min(p7/grid.get_bottom_dim()+
173
																	CGLA::Vec3i(1),
174
																	grid.get_top_dims());
175
		for(CGLA::Vec3i pt(p0t); pt[2]<p7t[2]; ++pt[2])
176
			for(pt[1]=p0t[1];      pt[1]<p7t[1]; ++pt[1])
177
				for(pt[0]=p0t[0];    pt[0]<p7t[0]; ++pt[0])
178
					functor(pt*CellT::get_dim(), grid.get_cell(pt));
179
	}
180
 
181
	template<class T, class CellT, class F>
182
	void for_each_cell(HGrid<T,CellT>& grid,
183
										 F& functor)
184
	{
185
		CGLA::Vec3i p0t;
186
		CGLA::Vec3i p7t =	grid.get_dims();
187
		const int inc = CellT::get_dim();
188
		int l=0;
189
		for(CGLA::Vec3i pt(p0t); pt[2]<p7t[2]; pt[2]+=inc)
190
			for(pt[1]=p0t[1];      pt[1]<p7t[1]; pt[1]+=inc)
191
				for(pt[0]=p0t[0];    pt[0]<p7t[0]; pt[0]+=inc)
192
					functor(pt, grid.get_cell(l++));
193
	}
194
 
195
 
196
	template<class CellT, class F>
197
	class _HGridCellFunctor
198
	{
199
		const CGLA::Vec3i p0;
200
		const CGLA::Vec3i p7;
201
		F& functor;
202
 
203
	public:
204
		_HGridCellFunctor(const CGLA::Vec3i _p0,
205
											const CGLA::Vec3i _p7,
206
											F& _functor): p0(_p0), p7(_p7), functor(_functor) {}
207
 
208
		void operator()(const CGLA::Vec3i& offset, 
209
										CellT& cell)
210
		{
211
			CGLA::Vec3i p0c = CGLA::v_max(p0-offset, CGLA::Vec3i(0));
212
			CGLA::Vec3i p7c = CGLA::v_min(p7-offset, CGLA::Vec3i(CellT::get_dim()));
213
 
214
			if(cell.is_coalesced())
215
				cell.split();
216
 
217
			_for_each_voxel(cell.get(), 
218
											CellT::get_dim(), 
219
											CGLA::sqr(CellT::get_dim()),
220
											p0c, p7c,	functor, offset);
221
		}
222
	};
223
 
224
 
225
	template<class T, class CellT, class F>
226
	void for_each_voxel(HGrid<T,CellT>& grid,
227
											const CGLA::Vec3i& _p0, 
228
											const CGLA::Vec3i& _p7,
229
											F& functor)
230
	{
231
		CGLA::Vec3i p0 = CGLA::v_max(_p0, CGLA::Vec3i(0));
232
		CGLA::Vec3i p7 = CGLA::v_min(_p7, grid.get_dims());
233
		_HGridCellFunctor<CellT,F> cell_functor(p0, p7, functor);
234
		for_each_cell(grid, p0, p7, cell_functor);
235
	}
236
 
237
	template<class T, class CellT, class F>
238
	void for_each_voxel(HGrid<T,CellT>& grid,	F& functor)
239
	{
240
		_HGridCellFunctor<CellT,F> cell_functor(CGLA::Vec3i(0), 
241
																						grid.get_dims(), functor);
242
		for_each_cell(grid, cell_functor);
243
	}
244
 
245
	template<class T, class CellT, class F>
246
	void for_each_voxel_ordered(HGrid<T,CellT>& grid,
247
															const CGLA::Vec3i& _p0, 
248
															const CGLA::Vec3i& _p7,
249
															F& functor)
250
	{
251
		CGLA::Vec3i p0 = CGLA::v_max(_p0, CGLA::Vec3i(0));
252
		CGLA::Vec3i p7 = CGLA::v_min(_p7, grid.get_dims());
253
		for(int k=p0[2];k<p7[2];++k)
254
			for(int j=p0[1];j<p7[1];++j)
255
				for(int i=p0[0];i<p7[0];++i)
256
					{
257
						CGLA::Vec3i p(i,j,k);
258
						float val = grid[p];
259
						float nval = val;
260
						functor(p, nval);
261
						if(nval != val)
262
							grid.store(p, nval);
263
					}
264
	}
265
 
266
	template<class T, class CellT, class F>
267
	void for_each_voxel_ordered(HGrid<T,CellT>& grid,	F& functor)
268
	{
269
		for_each_voxel_ordered(grid, CGLA::Vec3i(0), grid.get_dims(), functor);
270
	}
271
 
272
 
273
	template<typename T>
274
	class _AssignFun
275
	{
276
		T val;
277
	public:
278
		_AssignFun(const T& _val): val(_val) {}
279
		void operator()(const CGLA::Vec3i& pi, T& vox_val)
280
		{
281
			vox_val = val;
282
		}
283
	};
284
 
285
 	template<class G>
286
	void clear_region(G& grid, const typename G::DataType& value)
287
	{
288
		_AssignFun<typename G::DataType> afun(value);
289
		for_each_voxel(grid, afun);
290
	}
291
 
292
	template<class G>
293
	void clear_region(G& grid, 
294
										const CGLA::Vec3i& p0,
295
										const CGLA::Vec3i& p7,
296
										const typename G::DataType& value)
297
	{
298
		_AssignFun<typename G::DataType>afun(value) ;
299
		for_each_voxel(grid, p0, p7, afun);
300
	}
301
 
302
 
303
	//----------------------------------------------------------------------
304
	// const versions.
305
 
306
	/** Loop over all voxels in a sub-region (slice) of an
307
			RGrid and invoke a functor on each voxel. 
308
			The grid is the first argument, the slice is
309
			specified by the two subsequent args, and the functor 
310
			is the last argument. */
311
	template<class T, class F>
312
	void for_each_voxel_const(const RGrid<T>& grid,
313
														const CGLA::Vec3i& p0, 
314
														const CGLA::Vec3i& p7,
315
														F& functor)
316
	{
317
		_for_each_voxel(grid.get(), grid.get_x_dim(), grid.get_xy_dim(),
318
										CGLA::v_max(p0, CGLA::Vec3i(0)),
319
										CGLA::v_min(p7, grid.get_dims()), 
320
										functor);
321
	}
322
 
323
	/** Loop over all voxels in an entire	RGrid. 
324
			Grid is the first argument, and a functor is 
325
			the second.	For each voxel, an operation 
326
			specified by the functor is performed. */
327
	template<class T, class F>
328
	void for_each_voxel_const(const RGrid<T>& grid,	F& functor)
329
	{
330
		_for_each_voxel(grid.get(), grid.get_dims(), functor);
331
	}
332
 
333
 
334
	/** For each voxel (ordered). The idea of ordered traversal is 
335
			that we traverse the volume in a systematic fashion as opposed
336
			to traversing simply according to the memory layout of the volume
337
			data structure. This is important e.g. if we want to save the 
338
			volume in raw format. 
339
			For an RGrid, there is no difference though.
340
	*/ 
341
	template<class T, class F>
342
	void for_each_voxel_ordered_const(const RGrid<T>& grid,
343
																		const CGLA::Vec3i& p0, 
344
																		const CGLA::Vec3i& p7,
345
																		F& functor)
346
	{
347
		_for_each_voxel(grid.get(), grid.get_x_dim(), grid.get_xy_dim(),
348
										CGLA::v_max(p0, CGLA::Vec3i(0)),
349
										CGLA::v_min(p7, grid.get_dims()), 
350
										functor);
351
	}
352
 
353
	template<class T, class F>
354
	void for_each_voxel_ordered_const(const RGrid<T>& grid,	F& functor)
355
	{
356
		_for_each_voxel(grid.get(), grid.get_dims(), functor);
357
	}
358
 
359
 
360
	template<class T, class CellT, class F>
361
	void for_each_cell_const(const HGrid<T,CellT>& grid,
362
													 const CGLA::Vec3i& p0, 
363
													 const CGLA::Vec3i& p7,
364
													 F& functor)
365
	{
366
		CGLA::Vec3i p0t = p0/grid.get_bottom_dim();
367
		CGLA::Vec3i p7t = CGLA::v_min(p7/grid.get_bottom_dim()+
368
																	CGLA::Vec3i(1),
369
																	grid.get_top_dims());
370
		for(CGLA::Vec3i pt(p0t); pt[2]<p7t[2]; ++pt[2])
371
			for(pt[1]=p0t[1];      pt[1]<p7t[1]; ++pt[1])
372
				for(pt[0]=p0t[0];    pt[0]<p7t[0]; ++pt[0])
373
					functor(pt*CellT::get_dim(), grid.get_cell(pt));
374
	}
375
 
376
	template<class T, class CellT, class F>
377
	void for_each_cell_const(const HGrid<T,CellT>& grid,
378
													 F& functor)
379
	{
380
		CGLA::Vec3i p0t;
381
		CGLA::Vec3i p7t =	grid.get_dims();
382
		const int inc = CellT::get_dim();
383
		int l=0;
384
		for(CGLA::Vec3i pt(p0t); pt[2]<p7t[2]; pt[2]+=inc)
385
			for(pt[1]=p0t[1];      pt[1]<p7t[1]; pt[1]+=inc)
386
				for(pt[0]=p0t[0];    pt[0]<p7t[0]; pt[0]+=inc)
387
					functor(pt, grid.get_cell(l++));
388
	}
389
 
390
 
391
	template<class CellT, class F>
392
	class _HGridCellFunctorConst
393
	{
394
		const CGLA::Vec3i p0;
395
		const CGLA::Vec3i p7;
396
		F& functor;
397
 
398
	public:
399
		_HGridCellFunctorConst(const CGLA::Vec3i _p0,
400
													 const CGLA::Vec3i _p7,
401
													 F& _functor): p0(_p0), p7(_p7), functor(_functor) {}
402
 
403
		void operator()(const CGLA::Vec3i& offset, 
404
										const CellT& cell)
405
		{
406
			CGLA::Vec3i p0c = CGLA::v_max(p0-offset, CGLA::Vec3i(0));
407
			CGLA::Vec3i p7c = CGLA::v_min(p7-offset, CGLA::Vec3i(CellT::get_dim()));
408
 
409
			if(cell.is_coalesced())
410
				{
411
					typename CellT::DataType val = *cell.get();
412
					for(CGLA::Vec3i p(p0c); p[2]<p7c[2]; ++p[2])
413
						for(p[1]=p0c[1];      p[1]<p7c[1]; ++p[1])
414
							for(p[0]=p0c[0];    p[0]<p7c[0]; ++p[0])
415
								functor(p+offset, val);
416
				}
417
			_for_each_voxel(cell.get(), 
418
											CellT::get_dim(), 
419
											CGLA::sqr(CellT::get_dim()),
420
											p0c, p7c,	functor, offset);
421
		}
422
	};
423
 
424
 
425
	template<class T, class CellT, class F>
426
	void for_each_voxel_const(const HGrid<T,CellT>& grid,
427
														const CGLA::Vec3i& _p0, 
428
														const CGLA::Vec3i& _p7,
429
														F& functor)
430
	{
431
		CGLA::Vec3i p0 = CGLA::v_max(_p0, CGLA::Vec3i(0));
432
		CGLA::Vec3i p7 = CGLA::v_min(_p7, grid.get_dims());
433
		_HGridCellFunctorConst<CellT,F> cell_functor(p0, p7, functor);
434
		for_each_cell_const(grid, p0, p7, cell_functor);
435
	}
436
 
437
	template<class T, class CellT, class F>
438
	void for_each_voxel_const(const HGrid<T,CellT>& grid,	F& functor)
439
	{
440
		_HGridCellFunctorConst<CellT,F> cell_functor(CGLA::Vec3i(0), 
441
																								 grid.get_dims(), functor);
442
		for_each_cell_const(grid, cell_functor);
443
	}
444
 
445
	template<class T, class CellT, class F>
446
	void for_each_voxel_ordered_const(const HGrid<T,CellT>& grid,
447
																		const CGLA::Vec3i& _p0, 
448
																		const CGLA::Vec3i& _p7,
449
																		F& functor)
450
	{
451
		CGLA::Vec3i p0 = CGLA::v_max(_p0, CGLA::Vec3i(0));
452
		CGLA::Vec3i p7 = CGLA::v_min(_p7, grid.get_dims());
453
		for(int k=p0[2];k<p7[2];++k)
454
			for(int j=p0[1];j<p7[1];++j)
455
				for(int i=p0[0];i<p7[0];++i)
456
					{
457
						CGLA::Vec3i p(i,j,k);
458
						functor(p, grid[p]);
459
					}
460
	}
461
 
462
	template<class T, class CellT, class F>
463
	void for_each_voxel_ordered_const(const HGrid<T,CellT>& grid,	F& functor)
464
	{
465
		for_each_voxel_ordered_const(grid, 
466
																 CGLA::Vec3i(0), 
467
																 grid.get_dims(), 
468
																 functor);
469
	}
470
 
471
 
472
}
473
 
474
#endif