WebSVN – gelsvn – Diff – /branches/SPD/src/Geometry/ThreeDDDA.cpp


#include <cmath>
#include <stdlib.h>
#include <iostream>
#include "CGLA/Vec3d.h"
#include "ThreeDDDA.h"
 
using namespace std;
using namespace CGLA;
 
namespace Geometry
{
 
 
	/** 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 Vec3i ThreeDDDA::grid_pos(const Vec3f& v) const
	{
		Vec3i g;
		for(int i=0;i<3;++i)
			g[i] = static_cast<int>(floor(v[i]*PRECISION));
		return g;
	}
 
 
	/// Compute the cell containing a given fine grid position.
	Vec3i ThreeDDDA::get_cell(const Vec3i& v) const
	{
		Vec3i c;
		for(int i=0;i<3;++i)
			{
				// If we are on the right side of the axis
				if(v[i]>=0)
					// We simply divide by the precision.
					c[i] = (v[i])/PRECISION;
				else
					// If we are on the left side, we have to add one
					// before and subtract one after dividing. This reflects
					// the lack of symmetry.
					c[i] = (v[i]+1)/PRECISION-1;
			}
		return c;
	}
 
	/// Mirror a fine grid position
	Vec3i ThreeDDDA::mirror(const Vec3i& v) const
	{
		Vec3i m;
		for(int i=0;i<3;++i)
			{
				// Mirroring simply inverts the sign, however, 
				// to reflect the fact that the cells on either side of the
				// axis are numbered differently, we have to add one before flipping 
				// the sign.
				if(sgn[i]==-1)
					m[i] = -(v[i]+1);
				else
					m[i] = v[i];
			}
		return m;
	}
 
 
	/** Create a 3DDDA based on the two end-points of the ray. An optional
			last argument indicates the resolution of the grid. */
	ThreeDDDA::ThreeDDDA(const CGLA::Vec3f& p0, const CGLA::Vec3f& p1, 
											 int _PRECISION):
		PRECISION(_PRECISION),
		sgn(p1[0]>=p0[0] ? 1 : -1,
				p1[1]>=p0[1] ? 1 : -1,
				p1[2]>=p0[2] ? 1 : -1),
		p0_int(grid_pos(p0)),
		p1_int(grid_pos(p1)),
		dir(sgn*(p1_int-p0_int)),
		first_cell(get_cell(p0_int)),
		last_cell(get_cell(p1_int)),
		no_steps(dot(sgn, last_cell-first_cell)),
		dir_pre_mul(2*dir*PRECISION)
	{
		// Compute the mirrored position of the ray starting point
		const Vec3i p0_int_m = mirror(p0_int);
 
		// Compute the cell of the mirrored position.
		const Vec3i first_cell_m = get_cell(p0_int_m);
 
		/* Finally, compute delta = p - p0 where 
			 "p" is the position of the far,upper,right corner of the cell containing
			 the ray starting point and "p0" is the starting point (mirrored). 
			 Since the ray starts at a grid position + (0.5,0.5,0.5) we have to
			 add one half to each grid coordinate of p0. Since we are doing integer
			 arithmetic, we simply multiply everything by two instead.
		*/
		const Vec3i delta = 2*(first_cell_m*PRECISION+Vec3i(PRECISION))-
			(2*p0_int_m  + Vec3i(1));
 
		/* Compute the discriminators. These are (x-x0) * dir_y, (y-y0)*dir_x
			 and so on for all permutations. These will be updated in the course
			 of the algorithm. */ 
		for(int i=0;i<3;++i)
			{
				const int k1 = (i + 1) % 3;
				const int k2 = (i + 2) % 3;
				disc[i][k1] = static_cast<long long int>(delta[i])*dir[k1];
				disc[i][k2] = static_cast<long long int>(delta[i])*dir[k2];
			}
 
		// Finally, we set the cell position equal to the first_cell and we
		// are ready to iterate.
		cell = first_cell;
	}
 
	/// Get the initial point containing the ray.
	const CGLA::Vec3f ThreeDDDA::step() 
	{
		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;
		
		double delta;
		delta=cell[step_dir]+sgn[step_dir]-double(p0_int[step_dir]+0.5f)/PRECISION;
		Vec3f real_dir = get_last_point() - get_first_point();
		double ck1 = real_dir[k1]/real_dir[step_dir];
		double ck2 = real_dir[k2]/real_dir[step_dir];
 
		Vec3f p;
		p[step_dir] = delta;
		p[k1] = ck1 * delta;
		p[k2] = ck2 * delta;
 		p += get_first_point();
 
		cell[step_dir] += sgn[step_dir];
		disc[step_dir][k1] += dir_pre_mul[k1];
		disc[step_dir][k2] += dir_pre_mul[k2];
		return Vec3f(p);
	}
 
	
}
 

Rev 61	Rev 207
1	`#include <cmath>`	1	`#include <cmath>`
2	`#include <stdlib.h>`	2	`#include <stdlib.h>`
3	`#include <iostream>`	3	`#include <iostream>`
4	`#include "CGLA/Vec3d.h"`	4	`#include "CGLA/Vec3d.h"`
5	`#include "ThreeDDDA.h"`	5	`#include "ThreeDDDA.h"`
6		6
7	`using namespace std;`	7	`using namespace std;`
8	`using namespace CGLA;`	8	`using namespace CGLA;`
9		9
10	`namespace Geometry`	10	`namespace Geometry`
11	`{`	11	`{`
12		12
13		13
14	`/** Return the position on the fine grid of a vector v.`	14	`/** Return the position on the fine grid of a vector v.`
15	`We simply multiply by the precision and round down. Note that`	15	`We simply multiply by the precision and round down. Note that`
16	`it is necessary to use floor since simply clamping the value`	16	`it is necessary to use floor since simply clamping the value`
17	`will always round toward zero.`	17	`will always round toward zero.`
18	`*/`	18	`*/`
19	`const Vec3i ThreeDDDA::grid_pos(const Vec3f& v) const`	19	`const Vec3i ThreeDDDA::grid_pos(const Vec3f& v) const`
20	`{`	20	`{`
21	`Vec3i g;`	21	`Vec3i g;`
22	`for(int i=0;i<3;++i)`	22	`for(int i=0;i<3;++i)`
23	`g[i] = static_cast<int>(floor(v[i]*PRECISION));`	23	`g[i] = static_cast<int>(floor(v[i]*PRECISION));`
24	`return g;`	24	`return g;`
25	`}`	25	`}`
26		26
27		27
28	`/// Compute the cell containing a given fine grid position.`	28	`/// Compute the cell containing a given fine grid position.`
29	`Vec3i ThreeDDDA::get_cell(const Vec3i& v) const`	29	`Vec3i ThreeDDDA::get_cell(const Vec3i& v) const`
30	`{`	30	`{`
31	`Vec3i c;`	31	`Vec3i c;`
32	`for(int i=0;i<3;++i)`	32	`for(int i=0;i<3;++i)`
33	`{`	33	`{`
34	`// If we are on the right side of the axis`	34	`// If we are on the right side of the axis`
35	`if(v[i]>=0)`	35	`if(v[i]>=0)`
36	`// We simply divide by the precision.`	36	`// We simply divide by the precision.`
37	`c[i] = (v[i])/PRECISION;`	37	`c[i] = (v[i])/PRECISION;`
38	`else`	38	`else`
39	`// If we are on the left side, we have to add one`	39	`// If we are on the left side, we have to add one`
40	`// before and subtract one after dividing. This reflects`	40	`// before and subtract one after dividing. This reflects`
41	`// the lack of symmetry.`	41	`// the lack of symmetry.`
42	`c[i] = (v[i]+1)/PRECISION-1;`	42	`c[i] = (v[i]+1)/PRECISION-1;`
43	`}`	43	`}`
44	`return c;`	44	`return c;`
45	`}`	45	`}`
46		46
47	`/// Mirror a fine grid position`	47	`/// Mirror a fine grid position`
48	`Vec3i ThreeDDDA::mirror(const Vec3i& v) const`	48	`Vec3i ThreeDDDA::mirror(const Vec3i& v) const`
49	`{`	49	`{`
50	`Vec3i m;`	50	`Vec3i m;`
51	`for(int i=0;i<3;++i)`	51	`for(int i=0;i<3;++i)`
52	`{`	52	`{`
53	`// Mirroring simply inverts the sign, however,`	53	`// Mirroring simply inverts the sign, however,`
54	`// to reflect the fact that the cells on either side of the`	54	`// to reflect the fact that the cells on either side of the`
55	`// axis are numbered differently, we have to add one before flipping`	55	`// axis are numbered differently, we have to add one before flipping`
56	`// the sign.`	56	`// the sign.`
57	`if(sgn[i]==-1)`	57	`if(sgn[i]==-1)`
58	`m[i] = -(v[i]+1);`	58	`m[i] = -(v[i]+1);`
59	`else`	59	`else`
60	`m[i] = v[i];`	60	`m[i] = v[i];`
61	`}`	61	`}`
62	`return m;`	62	`return m;`
63	`}`	63	`}`
64		64
65		65
66	`/** Create a 3DDDA based on the two end-points of the ray. An optional`	66	`/** Create a 3DDDA based on the two end-points of the ray. An optional`
67	`last argument indicates the resolution of the grid. */`	67	`last argument indicates the resolution of the grid. */`
68	`ThreeDDDA::ThreeDDDA(const CGLA::Vec3f& p0, const CGLA::Vec3f& p1,`	68	`ThreeDDDA::ThreeDDDA(const CGLA::Vec3f& p0, const CGLA::Vec3f& p1,`
69	`int _PRECISION):`	69	`int _PRECISION):`
70	`PRECISION(_PRECISION),`	70	`PRECISION(_PRECISION),`
71	`sgn(p1[0]>=p0[0] ? 1 : -1,`	71	`sgn(p1[0]>=p0[0] ? 1 : -1,`
72	`p1[1]>=p0[1] ? 1 : -1,`	72	`p1[1]>=p0[1] ? 1 : -1,`
73	`p1[2]>=p0[2] ? 1 : -1),`	73	`p1[2]>=p0[2] ? 1 : -1),`
74	`p0_int(grid_pos(p0)),`	74	`p0_int(grid_pos(p0)),`
75	`p1_int(grid_pos(p1)),`	75	`p1_int(grid_pos(p1)),`
76	`dir(sgn*(p1_int-p0_int)),`	76	`dir(sgn*(p1_int-p0_int)),`
77	`first_cell(get_cell(p0_int)),`	77	`first_cell(get_cell(p0_int)),`
78	`last_cell(get_cell(p1_int)),`	78	`last_cell(get_cell(p1_int)),`
79	`no_steps(dot(sgn, last_cell-first_cell)),`	79	`no_steps(dot(sgn, last_cell-first_cell)),`
80	`dir_pre_mul(2dirPRECISION)`	80	`dir_pre_mul(2dirPRECISION)`
81	`{`	81	`{`
82	`// Compute the mirrored position of the ray starting point`	82	`// Compute the mirrored position of the ray starting point`
83	`const Vec3i p0_int_m = mirror(p0_int);`	83	`const Vec3i p0_int_m = mirror(p0_int);`
84		84
85	`// Compute the cell of the mirrored position.`	85	`// Compute the cell of the mirrored position.`
86	`const Vec3i first_cell_m = get_cell(p0_int_m);`	86	`const Vec3i first_cell_m = get_cell(p0_int_m);`
87		87
88	`/* Finally, compute delta = p - p0 where`	88	`/* Finally, compute delta = p - p0 where`
89	`"p" is the position of the far,upper,right corner of the cell containing`	89	`"p" is the position of the far,upper,right corner of the cell containing`
90	`the ray starting point and "p0" is the starting point (mirrored).`	90	`the ray starting point and "p0" is the starting point (mirrored).`
91	`Since the ray starts at a grid position + (0.5,0.5,0.5) we have to`	91	`Since the ray starts at a grid position + (0.5,0.5,0.5) we have to`
92	`add one half to each grid coordinate of p0. Since we are doing integer`	92	`add one half to each grid coordinate of p0. Since we are doing integer`
93	`arithmetic, we simply multiply everything by two instead.`	93	`arithmetic, we simply multiply everything by two instead.`
94	`*/`	94	`*/`
95	`const Vec3i delta = 2(first_cell_mPRECISION+Vec3i(PRECISION))-`	95	`const Vec3i delta = 2(first_cell_mPRECISION+Vec3i(PRECISION))-`
96	`(2*p0_int_m + Vec3i(1));`	96	`(2*p0_int_m + Vec3i(1));`
97		97
98	`/* Compute the discriminators. These are (x-x0) * dir_y, (y-y0)*dir_x`	98	`/* Compute the discriminators. These are (x-x0) * dir_y, (y-y0)*dir_x`
99	`and so on for all permutations. These will be updated in the course`	99	`and so on for all permutations. These will be updated in the course`
100	`of the algorithm. */`	100	`of the algorithm. */`
101	`for(int i=0;i<3;++i)`	101	`for(int i=0;i<3;++i)`
102	`{`	102	`{`
103	`const int k1 = (i + 1) % 3;`	103	`const int k1 = (i + 1) % 3;`
104	`const int k2 = (i + 2) % 3;`	104	`const int k2 = (i + 2) % 3;`
105	`disc[i][k1] = static_cast<long long int>(delta[i])*dir[k1];`	105	`disc[i][k1] = static_cast<long long int>(delta[i])*dir[k1];`
106	`disc[i][k2] = static_cast<long long int>(delta[i])*dir[k2];`	106	`disc[i][k2] = static_cast<long long int>(delta[i])*dir[k2];`
107	`}`	107	`}`
108		108
109	`// Finally, we set the cell position equal to the first_cell and we`	109	`// Finally, we set the cell position equal to the first_cell and we`
110	`// are ready to iterate.`	110	`// are ready to iterate.`
111	`cell = first_cell;`	111	`cell = first_cell;`
112	`}`	112	`}`
113		113
114	`/// Get the initial point containing the ray.`	114	`/// Get the initial point containing the ray.`
115	`const CGLA::Vec3f ThreeDDDA::step()`	115	`const CGLA::Vec3f ThreeDDDA::step()`
116	`{`	116	`{`
117	`int step_dir;`	117	`int step_dir;`
118	`if(disc[1][0] > disc[0][1])`	118	`if(disc[1][0] > disc[0][1])`
119	`if(disc[2][0] > disc[0][2])`	119	`if(disc[2][0] > disc[0][2])`
120	`step_dir = 0;`	120	`step_dir = 0;`
121	`else`	121	`else`
122	`step_dir = 2;`	122	`step_dir = 2;`
123	`else`	123	`else`
124	`if(disc[2][1] > disc[1][2])`	124	`if(disc[2][1] > disc[1][2])`
125	`step_dir = 1;`	125	`step_dir = 1;`
126	`else`	126	`else`
127	`step_dir = 2;`	127	`step_dir = 2;`
128		128
129	`const int k1 = (step_dir + 1) % 3;`	129	`const int k1 = (step_dir + 1) % 3;`
130	`const int k2 = (step_dir + 2) % 3;`	130	`const int k2 = (step_dir + 2) % 3;`
131		131
132	`double delta;`	132	`double delta;`
133	`delta=cell[step_dir]+sgn[step_dir]-double(p0_int[step_dir]+0.5f)/PRECISION;`	133	`delta=cell[step_dir]+sgn[step_dir]-double(p0_int[step_dir]+0.5f)/PRECISION;`
134	`Vec3f real_dir = get_last_point() - get_first_point();`	134	`Vec3f real_dir = get_last_point() - get_first_point();`
135	`double ck1 = real_dir[k1]/real_dir[step_dir];`	135	`double ck1 = real_dir[k1]/real_dir[step_dir];`
136	`double ck2 = real_dir[k2]/real_dir[step_dir];`	136	`double ck2 = real_dir[k2]/real_dir[step_dir];`
137		137
138	`Vec3f p;`	138	`Vec3f p;`
139	`p[step_dir] = delta;`	139	`p[step_dir] = delta;`
140	`p[k1] = ck1 * delta;`	140	`p[k1] = ck1 * delta;`
141	`p[k2] = ck2 * delta;`	141	`p[k2] = ck2 * delta;`
142	`p += get_first_point();`	142	`p += get_first_point();`
143		143
144	`cell[step_dir] += sgn[step_dir];`	144	`cell[step_dir] += sgn[step_dir];`
145	`disc[step_dir][k1] += dir_pre_mul[k1];`	145	`disc[step_dir][k1] += dir_pre_mul[k1];`
146	`disc[step_dir][k2] += dir_pre_mul[k2];`	146	`disc[step_dir][k2] += dir_pre_mul[k2];`
147	`return Vec3f(p);`	147	`return Vec3f(p);`
148	`}`	148	`}`
149		149
150		150
151	`}`	151	`}`
152		152

Subversion Repositories gelsvn

(root)/branches/SPD/src/Geometry/ThreeDDDA.cpp – Rev 61 → 207