WebSVN – gelsvn – Diff – /trunk/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)/trunk/src/Geometry/ThreeDDDA.h – Rev 443 → 595

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