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