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

Subversion Repositories gelsvn

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