WebSVN – gelsvn – Diff – /trunk/src/Geometry/ThreeDDDA.cpp


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

Subversion Repositories gelsvn

(root)/trunk/src/Geometry/ThreeDDDA.cpp – Rev 304 → 595