WebSVN – gelsvn – Diff – /branches/cpp11-devel/src/Geometry/ThreeDDDA.h


/* ----------------------------------------------------------------------- *
 * This file is part of GEL, http://www.imm.dtu.dk/GEL
 * Copyright (C) the authors and DTU Informatics
 * For license and list of authors, see ../../doc/intro.pdf
 * ----------------------------------------------------------------------- */
 
/**
 * @file ThreeDDDA.h
 * @brief a class for traversing cells in a regular grid pierced by a ray.
 */
#ifndef __GEOMETRY_THREEDDDA_H
#define __GEOMETRY_THREEDDDA_H
 
#ifdef _MSC_VER 
typedef __int64 Integer64Bits;
#else
typedef int64_t Integer64Bits;
#endif
 
 
#include "../CGLA/Vec3f.h"
#include "../CGLA/Vec3i.h"
 
namespace Geometry
{
 
 
	/** \brief A ThreeDDDA is a class for traversing a grid of cells. 
 
      We wish to 
			enumerate all cells pierced by the ray going from a point p0 to a 
			point p1. It is dangerous to use floating point arithmetic since rounding
			errors may accumulate entailing that the 3ddda misses the cell containing
			the end point of the ray.
		
			Daniel Cohen-Or devised a technique for exact ray traversal using only 
			integer arithmetic. Using integer arithmetic, the computations are
			exact, but the ray must start and terminate exactly at grid points.
			I.e. the ray must start and terminate in cell corners.
 
			To overcome this issue, I have implemented a scheme where the ray
			starts and terminates on a grid that is finer than the cell grid.
			This means that we have fast and exact computations, and the ray can
			travel between two almost arbitrary points.
 
			A central problem with this scheme was that - in order to simplify
			things during the iterative traversal - the ray direction vector
			is always mirrored into the first octant (+x +y +z).  If the fine
			grid used for the ray start and end points is a superset of the
			grid points of the cell grid, this mirroring is almost impossible
			to implement correctly. This is due to the fact that the first
			cell on the right side of the y axis is called 0,y whereas the
			first cell on the left is called -1,y. This lack of symmetry makes
			things very difficult, but by ensuring that the endpoints of a ray
			can never lie on an a corner, edge, or face of the cell grid, we
			can make the problems go away.
 
			Hence, in this implementation, the fine grid is obtained by dividing 
			the cell into NxNxN subcells, but the grid points	of the fine grid lie 
			at the CENTRES of these subcells.  
		
	*/
 
	class ThreeDDDA
	{
		/** Resolution of the grid. This number indicates how many subcells
				(along each axis) a cell is divided into. */
		const int PRECISION;
 
		/// Sign indicates which octant contains the direction vector.
		const CGLA::Vec3i sgn;
 
		/// Fine grid position where the ray begins
		const CGLA::Vec3i p0_int;
 
		// Fine grid position where the ray ends.
		const CGLA::Vec3i p1_int;
 
		/// The direction of the ray 
		const CGLA::Vec3i dir;
 
		/// The cell containing the initial point on the ray.
		const CGLA::Vec3i first_cell;
 
		/// The cell containing the last point on the ray.
		const CGLA::Vec3i last_cell;
 
		/** Number of steps. We can compute the number of steps that the
				ray will take in advance. */
		const int no_steps; 
 
		/** The direction vector, premultiplied by the step length.
				This constant is used to update the equations governing the
				traversal. */
		const CGLA::Vec3i dir_pre_mul;
 
		/** Discriminators. These values are delta[i]*dir[j] where delta[i]
				is the distance from the origin of the ray to the current position
				along the i axis multiplied by two). dir[j] is the direction vector
				coordinate j. */
	
		Integer64Bits disc[3][3];
 
		/** The current cell containing the ray. This value along with the 
				discriminators mentioned above represent the state of a ThreeDDDA. */
		CGLA::Vec3i cell;
 
 
		/** Return the position on the fine grid of a vector v. 
				We simply multiply by the precision and round down. Note that 
				it is necessary to use floor since simply clamping the value
				will always round toward zero.
		*/
		const CGLA::Vec3i grid_pos(const CGLA::Vec3f& v) const;
	
		/// Compute the cell containing a given fine grid position.
		CGLA::Vec3i get_cell(const CGLA::Vec3i& v) const;
	
		/// Mirror a fine grid position
		CGLA::Vec3i mirror(const CGLA::Vec3i& v) const;
 
 
	public:
 
		/** Create a 3DDDA based on the two end-points of the ray. An optional
				last argument indicates the resolution of the grid. */
		ThreeDDDA(const CGLA::Vec3f& p0, const CGLA::Vec3f& p1, 
							int _PRECISION = 0x100);
	
		/// Return the total number of steps the ray will traverse.
		int get_no_steps() const
		{
			return no_steps;
		}
 
		/// Yields the current cell containing the ray.
		const CGLA::Vec3i& get_current_cell() const
		{
			return cell;
		}
 
		/// Get the initial cell containing the ray.
		const CGLA::Vec3i& get_first_cell() const
		{
			return first_cell;
		}
 
		/// Get the last cell containing the ray.
		const CGLA::Vec3i& get_last_cell() const
		{
			return last_cell;
		}
 
		/// Get the initial point containing the ray.
		const CGLA::Vec3f get_first_point() const
		{
			return (CGLA::Vec3f(p0_int)+CGLA::Vec3f(0.5))/PRECISION;
		}
 
 
		/// Get the last cell containing the ray.
		const CGLA::Vec3f get_last_point() const
		{
			return (CGLA::Vec3f(p1_int)+CGLA::Vec3f(0.5))/PRECISION;
		}
	
		/// Returns true if we have reached the last cell.
		bool at_end() const
		{
			return cell == last_cell;
		}
 
		/// Take a step along the ray.
		void operator++()
		{
			int step_dir;
			if(disc[1][0] > disc[0][1])
				if(disc[2][0] > disc[0][2])
					step_dir = 0;
				else
					step_dir = 2;
			else
				if(disc[2][1] > disc[1][2])
					step_dir = 1;
				else
					step_dir = 2;
				
			const int k1 = (step_dir + 1) % 3;
			const int k2 = (step_dir + 2) % 3;
			cell[step_dir] += sgn[step_dir];
			disc[step_dir][k1] += dir_pre_mul[k1];
			disc[step_dir][k2] += dir_pre_mul[k2];
		}
 
		/// Get the initial point containing the ray.
			const CGLA::Vec3f step();
 
 
	};
}
 
#endif
 

Subversion Repositories gelsvn

(root)/branches/cpp11-devel/src/Geometry/ThreeDDDA.h – Rev 601 → 630

Rev 601		Rev 630
1	`/* ----------------------------------------------------------------------- *`	1	`/* ----------------------------------------------------------------------- *`
2	`* This file is part of GEL, http://www.imm.dtu.dk/GEL`	2	`* This file is part of GEL, http://www.imm.dtu.dk/GEL`
3	`* Copyright (C) the authors and DTU Informatics`	3	`* Copyright (C) the authors and DTU Informatics`
4	`* For license and list of authors, see ../../doc/intro.pdf`	4	`* For license and list of authors, see ../../doc/intro.pdf`
5	`* ----------------------------------------------------------------------- */`	5	`* ----------------------------------------------------------------------- */`
6		6
7	`/**`	7	`/**`
8	`* @file ThreeDDDA.h`	8	`* @file ThreeDDDA.h`
9	`* @brief a class for traversing cells in a regular grid pierced by a ray.`	9	`* @brief a class for traversing cells in a regular grid pierced by a ray.`
10	`*/`	10	`*/`
11	`#ifndef __GEOMETRY_THREEDDDA_H`	11	`#ifndef __GEOMETRY_THREEDDDA_H`
12	`#define __GEOMETRY_THREEDDDA_H`	12	`#define __GEOMETRY_THREEDDDA_H`
13		13
14	`#ifdef _MSC_VER`	14	`#ifdef _MSC_VER`
15	`typedef __int64 Integer64Bits;`	15	`typedef __int64 Integer64Bits;`
16	`#else`	16	`#else`
17	`typedef int64_t Integer64Bits;`	17	`typedef int64_t Integer64Bits;`
18	`#endif`	18	`#endif`
19		19
20		20
21	`#include "../CGLA/Vec3f.h"`	21	`#include "../CGLA/Vec3f.h"`
22	`#include "../CGLA/Vec3i.h"`	22	`#include "../CGLA/Vec3i.h"`
23		23
24	`namespace Geometry`	24	`namespace Geometry`
25	`{`	25	`{`
26		26
27		27
28	`/** \brief A ThreeDDDA is a class for traversing a grid of cells.`	28	`/** \brief A ThreeDDDA is a class for traversing a grid of cells.`
29		29
30	`We wish to`	30	`We wish to`
31	`enumerate all cells pierced by the ray going from a point p0 to a`	31	`enumerate all cells pierced by the ray going from a point p0 to a`
32	`point p1. It is dangerous to use floating point arithmetic since rounding`	32	`point p1. It is dangerous to use floating point arithmetic since rounding`
33	`errors may accumulate entailing that the 3ddda misses the cell containing`	33	`errors may accumulate entailing that the 3ddda misses the cell containing`
34	`the end point of the ray.`	34	`the end point of the ray.`
35		35
36	`Daniel Cohen-Or devised a technique for exact ray traversal using only`	36	`Daniel Cohen-Or devised a technique for exact ray traversal using only`
37	`integer arithmetic. Using integer arithmetic, the computations are`	37	`integer arithmetic. Using integer arithmetic, the computations are`
38	`exact, but the ray must start and terminate exactly at grid points.`	38	`exact, but the ray must start and terminate exactly at grid points.`
39	`I.e. the ray must start and terminate in cell corners.`	39	`I.e. the ray must start and terminate in cell corners.`
40		40
41	`To overcome this issue, I have implemented a scheme where the ray`	41	`To overcome this issue, I have implemented a scheme where the ray`
42	`starts and terminates on a grid that is finer than the cell grid.`	42	`starts and terminates on a grid that is finer than the cell grid.`
43	`This means that we have fast and exact computations, and the ray can`	43	`This means that we have fast and exact computations, and the ray can`
44	`travel between two almost arbitrary points.`	44	`travel between two almost arbitrary points.`
45		45
46	`A central problem with this scheme was that - in order to simplify`	46	`A central problem with this scheme was that - in order to simplify`
47	`things during the iterative traversal - the ray direction vector`	47	`things during the iterative traversal - the ray direction vector`
48	`is always mirrored into the first octant (+x +y +z). If the fine`	48	`is always mirrored into the first octant (+x +y +z). If the fine`
49	`grid used for the ray start and end points is a superset of the`	49	`grid used for the ray start and end points is a superset of the`
50	`grid points of the cell grid, this mirroring is almost impossible`	50	`grid points of the cell grid, this mirroring is almost impossible`
51	`to implement correctly. This is due to the fact that the first`	51	`to implement correctly. This is due to the fact that the first`
52	`cell on the right side of the y axis is called 0,y whereas the`	52	`cell on the right side of the y axis is called 0,y whereas the`
53	`first cell on the left is called -1,y. This lack of symmetry makes`	53	`first cell on the left is called -1,y. This lack of symmetry makes`
54	`things very difficult, but by ensuring that the endpoints of a ray`	54	`things very difficult, but by ensuring that the endpoints of a ray`
55	`can never lie on an a corner, edge, or face of the cell grid, we`	55	`can never lie on an a corner, edge, or face of the cell grid, we`
56	`can make the problems go away.`	56	`can make the problems go away.`
57		57
58	`Hence, in this implementation, the fine grid is obtained by dividing`	58	`Hence, in this implementation, the fine grid is obtained by dividing`
59	`the cell into NxNxN subcells, but the grid points of the fine grid lie`	59	`the cell into NxNxN subcells, but the grid points of the fine grid lie`
60	`at the CENTRES of these subcells.`	60	`at the CENTRES of these subcells.`
61		61
62	`*/`	62	`*/`
63		63
64	`class ThreeDDDA`	64	`class ThreeDDDA`
65	`{`	65	`{`
66	`/** Resolution of the grid. This number indicates how many subcells`	66	`/** Resolution of the grid. This number indicates how many subcells`
67	`(along each axis) a cell is divided into. */`	67	`(along each axis) a cell is divided into. */`
68	`const int PRECISION;`	68	`const int PRECISION;`
69		69
70	`/// Sign indicates which octant contains the direction vector.`	70	`/// Sign indicates which octant contains the direction vector.`
71	`const CGLA::Vec3i sgn;`	71	`const CGLA::Vec3i sgn;`
72		72
73	`/// Fine grid position where the ray begins`	73	`/// Fine grid position where the ray begins`
74	`const CGLA::Vec3i p0_int;`	74	`const CGLA::Vec3i p0_int;`
75		75
76	`// Fine grid position where the ray ends.`	76	`// Fine grid position where the ray ends.`
77	`const CGLA::Vec3i p1_int;`	77	`const CGLA::Vec3i p1_int;`
78		78
79	`/// The direction of the ray`	79	`/// The direction of the ray`
80	`const CGLA::Vec3i dir;`	80	`const CGLA::Vec3i dir;`
81		81
82	`/// The cell containing the initial point on the ray.`	82	`/// The cell containing the initial point on the ray.`
83	`const CGLA::Vec3i first_cell;`	83	`const CGLA::Vec3i first_cell;`
84		84
85	`/// The cell containing the last point on the ray.`	85	`/// The cell containing the last point on the ray.`
86	`const CGLA::Vec3i last_cell;`	86	`const CGLA::Vec3i last_cell;`
87		87
88	`/** Number of steps. We can compute the number of steps that the`	88	`/** Number of steps. We can compute the number of steps that the`
89	`ray will take in advance. */`	89	`ray will take in advance. */`
90	`const int no_steps;`	90	`const int no_steps;`
91		91
92	`/** The direction vector, premultiplied by the step length.`	92	`/** The direction vector, premultiplied by the step length.`
93	`This constant is used to update the equations governing the`	93	`This constant is used to update the equations governing the`
94	`traversal. */`	94	`traversal. */`
95	`const CGLA::Vec3i dir_pre_mul;`	95	`const CGLA::Vec3i dir_pre_mul;`
96		96
97	`/** Discriminators. These values are delta[i]*dir[j] where delta[i]`	97	`/** Discriminators. These values are delta[i]*dir[j] where delta[i]`
98	`is the distance from the origin of the ray to the current position`	98	`is the distance from the origin of the ray to the current position`
99	`along the i axis multiplied by two). dir[j] is the direction vector`	99	`along the i axis multiplied by two). dir[j] is the direction vector`
100	`coordinate j. */`	100	`coordinate j. */`
101		101
102	`Integer64Bits disc[3][3];`	102	`Integer64Bits disc[3][3];`
103		103
104	`/** The current cell containing the ray. This value along with the`	104	`/** The current cell containing the ray. This value along with the`
105	`discriminators mentioned above represent the state of a ThreeDDDA. */`	105	`discriminators mentioned above represent the state of a ThreeDDDA. */`
106	`CGLA::Vec3i cell;`	106	`CGLA::Vec3i cell;`
107		107
108		108
109	`/** Return the position on the fine grid of a vector v.`	109	`/** Return the position on the fine grid of a vector v.`
110	`We simply multiply by the precision and round down. Note that`	110	`We simply multiply by the precision and round down. Note that`
111	`it is necessary to use floor since simply clamping the value`	111	`it is necessary to use floor since simply clamping the value`
112	`will always round toward zero.`	112	`will always round toward zero.`
113	`*/`	113	`*/`
114	`const CGLA::Vec3i grid_pos(const CGLA::Vec3f& v) const;`	114	`const CGLA::Vec3i grid_pos(const CGLA::Vec3f& v) const;`
115		115
116	`/// Compute the cell containing a given fine grid position.`	116	`/// Compute the cell containing a given fine grid position.`
117	`CGLA::Vec3i get_cell(const CGLA::Vec3i& v) const;`	117	`CGLA::Vec3i get_cell(const CGLA::Vec3i& v) const;`
118		118
119	`/// Mirror a fine grid position`	119	`/// Mirror a fine grid position`
120	`CGLA::Vec3i mirror(const CGLA::Vec3i& v) const;`	120	`CGLA::Vec3i mirror(const CGLA::Vec3i& v) const;`
121		121
122		122
123	`public:`	123	`public:`
124		124
125	`/** Create a 3DDDA based on the two end-points of the ray. An optional`	125	`/** Create a 3DDDA based on the two end-points of the ray. An optional`
126	`last argument indicates the resolution of the grid. */`	126	`last argument indicates the resolution of the grid. */`
127	`ThreeDDDA(const CGLA::Vec3f& p0, const CGLA::Vec3f& p1,`	127	`ThreeDDDA(const CGLA::Vec3f& p0, const CGLA::Vec3f& p1,`
128	`int _PRECISION = 0x100);`	128	`int _PRECISION = 0x100);`
129		129
130	`/// Return the total number of steps the ray will traverse.`	130	`/// Return the total number of steps the ray will traverse.`
131	`int get_no_steps() const`	131	`int get_no_steps() const`
132	`{`	132	`{`
133	`return no_steps;`	133	`return no_steps;`
134	`}`	134	`}`
135		135
136	`/// Yields the current cell containing the ray.`	136	`/// Yields the current cell containing the ray.`
137	`const CGLA::Vec3i& get_current_cell() const`	137	`const CGLA::Vec3i& get_current_cell() const`
138	`{`	138	`{`
139	`return cell;`	139	`return cell;`
140	`}`	140	`}`
141		141
142	`/// Get the initial cell containing the ray.`	142	`/// Get the initial cell containing the ray.`
143	`const CGLA::Vec3i& get_first_cell() const`	143	`const CGLA::Vec3i& get_first_cell() const`
144	`{`	144	`{`
145	`return first_cell;`	145	`return first_cell;`
146	`}`	146	`}`
147		147
148	`/// Get the last cell containing the ray.`	148	`/// Get the last cell containing the ray.`
149	`const CGLA::Vec3i& get_last_cell() const`	149	`const CGLA::Vec3i& get_last_cell() const`
150	`{`	150	`{`
151	`return last_cell;`	151	`return last_cell;`
152	`}`	152	`}`
153		153
154	`/// Get the initial point containing the ray.`	154	`/// Get the initial point containing the ray.`
155	`const CGLA::Vec3f get_first_point() const`	155	`const CGLA::Vec3f get_first_point() const`
156	`{`	156	`{`
157	`return (CGLA::Vec3f(p0_int)+CGLA::Vec3f(0.5))/PRECISION;`	157	`return (CGLA::Vec3f(p0_int)+CGLA::Vec3f(0.5))/PRECISION;`
158	`}`	158	`}`
159		159
160		160
161	`/// Get the last cell containing the ray.`	161	`/// Get the last cell containing the ray.`
162	`const CGLA::Vec3f get_last_point() const`	162	`const CGLA::Vec3f get_last_point() const`
163	`{`	163	`{`
164	`return (CGLA::Vec3f(p1_int)+CGLA::Vec3f(0.5))/PRECISION;`	164	`return (CGLA::Vec3f(p1_int)+CGLA::Vec3f(0.5))/PRECISION;`
165	`}`	165	`}`
166		166
167	`/// Returns true if we have reached the last cell.`	167	`/// Returns true if we have reached the last cell.`
168	`bool at_end() const`	168	`bool at_end() const`
169	`{`	169	`{`
170	`return cell == last_cell;`	170	`return cell == last_cell;`
171	`}`	171	`}`
172		172
173	`/// Take a step along the ray.`	173	`/// Take a step along the ray.`
174	`void operator++()`	174	`void operator++()`
175	`{`	175	`{`
176	`int step_dir;`	176	`int step_dir;`
177	`if(disc[1][0] > disc[0][1])`	177	`if(disc[1][0] > disc[0][1])`
178	`if(disc[2][0] > disc[0][2])`	178	`if(disc[2][0] > disc[0][2])`
179	`step_dir = 0;`	179	`step_dir = 0;`
180	`else`	180	`else`
181	`step_dir = 2;`	181	`step_dir = 2;`
182	`else`	182	`else`
183	`if(disc[2][1] > disc[1][2])`	183	`if(disc[2][1] > disc[1][2])`
184	`step_dir = 1;`	184	`step_dir = 1;`
185	`else`	185	`else`
186	`step_dir = 2;`	186	`step_dir = 2;`
187		187
188	`const int k1 = (step_dir + 1) % 3;`	188	`const int k1 = (step_dir + 1) % 3;`
189	`const int k2 = (step_dir + 2) % 3;`	189	`const int k2 = (step_dir + 2) % 3;`
190	`cell[step_dir] += sgn[step_dir];`	190	`cell[step_dir] += sgn[step_dir];`
191	`disc[step_dir][k1] += dir_pre_mul[k1];`	191	`disc[step_dir][k1] += dir_pre_mul[k1];`
192	`disc[step_dir][k2] += dir_pre_mul[k2];`	192	`disc[step_dir][k2] += dir_pre_mul[k2];`
193	`}`	193	`}`
194		194
195	`/// Get the initial point containing the ray.`	195	`/// Get the initial point containing the ray.`
196	`const CGLA::Vec3f step();`	196	`const CGLA::Vec3f step();`
197		197
198		198
199	`};`	199	`};`
200	`}`	200	`}`
201		201
202	`#endif`	202	`#endif`
203		203