1 | #include "Matrix4x4.h"
|
---|
2 | #include "Vector3.h"
|
---|
3 | #include "Halton.h"
|
---|
4 | #include <float.h>
|
---|
5 |
|
---|
6 | namespace GtpVisibilityPreprocessor {
|
---|
7 |
|
---|
8 | // Given min a vector to minimize and a candidate vector, replace
|
---|
9 | // elements of min whose corresponding elements in Candidate are
|
---|
10 | // smaller. This function is used for finding objects' bounds,
|
---|
11 | // among other things.
|
---|
12 | void
|
---|
13 | Minimize(Vector3 &min, const Vector3 &Candidate)
|
---|
14 | {
|
---|
15 | if (Candidate.x < min.x)
|
---|
16 | min.x = Candidate.x;
|
---|
17 | if (Candidate.y < min.y)
|
---|
18 | min.y = Candidate.y;
|
---|
19 | if (Candidate.z < min.z)
|
---|
20 | min.z = Candidate.z;
|
---|
21 | }
|
---|
22 |
|
---|
23 | // Given min a vector to minimize and a candidate vector, replace
|
---|
24 | // elements of min whose corresponding elements in Candidate are
|
---|
25 | // larger. This function is used for finding objects' bounds,
|
---|
26 | // among other things.
|
---|
27 | void
|
---|
28 | Maximize(Vector3 &max, const Vector3 &Candidate)
|
---|
29 | {
|
---|
30 | if (Candidate.x > max.x)
|
---|
31 | max.x = Candidate.x;
|
---|
32 | if (Candidate.y > max.y)
|
---|
33 | max.y = Candidate.y;
|
---|
34 | if (Candidate.z > max.z)
|
---|
35 | max.z = Candidate.z;
|
---|
36 | }
|
---|
37 |
|
---|
38 | // Project the vector onto the YZ, XZ, or XY plane depending on which.
|
---|
39 | // which Coordinate plane to project onto
|
---|
40 | // 0 YZ
|
---|
41 | // 1 XZ
|
---|
42 | // 2 XY
|
---|
43 | // This function is used by the polygon intersection code.
|
---|
44 |
|
---|
45 | void
|
---|
46 | Vector3::ExtractVerts(float *px, float *py, int which) const
|
---|
47 | {
|
---|
48 | switch (which) {
|
---|
49 | case 0:
|
---|
50 | *px = y;
|
---|
51 | *py = z;
|
---|
52 | break;
|
---|
53 | case 1:
|
---|
54 | *px = x;
|
---|
55 | *py = z;
|
---|
56 | break;
|
---|
57 | case 2:
|
---|
58 | *px = x;
|
---|
59 | *py = y;
|
---|
60 | break;
|
---|
61 | }
|
---|
62 | }
|
---|
63 |
|
---|
64 | // returns the axis, where the vector has the largest value
|
---|
65 | int
|
---|
66 | Vector3::DrivingAxis(void) const
|
---|
67 | {
|
---|
68 | int axis = 0;
|
---|
69 | float val = fabs(x);
|
---|
70 |
|
---|
71 | if (fabs(y) > val) {
|
---|
72 | val = fabs(y);
|
---|
73 | axis = 1;
|
---|
74 | }
|
---|
75 |
|
---|
76 | if (fabs(z) > val)
|
---|
77 | axis = 2;
|
---|
78 |
|
---|
79 | return axis;
|
---|
80 | }
|
---|
81 |
|
---|
82 | // returns the axis, where the vector has the smallest value
|
---|
83 | int
|
---|
84 | Vector3::TinyAxis(void) const
|
---|
85 | {
|
---|
86 | int axis = 0;
|
---|
87 | float val = fabs(x);
|
---|
88 |
|
---|
89 | if (fabs(y) < val) {
|
---|
90 | val = fabs(y);
|
---|
91 | axis = 1;
|
---|
92 | }
|
---|
93 |
|
---|
94 | if (fabs(z) < val)
|
---|
95 | axis = 2;
|
---|
96 |
|
---|
97 | return axis;
|
---|
98 | }
|
---|
99 |
|
---|
100 |
|
---|
101 | // Construct a view vector ViewN, and the vector ViewU perpendicular
|
---|
102 | // to ViewN and lying in the plane given by ViewNxUpl
|
---|
103 | // the last vector of ortogonal system is ViewV, that is
|
---|
104 | // perpendicular to both ViewN and ViewU.
|
---|
105 | // |ViewN| = |ViewU| = |ViewV| = 1
|
---|
106 | // The ViewN vector pierces the center of the synthesized image
|
---|
107 | // ViewU vector goes from the center image rightwards
|
---|
108 | // ViewV vector goes from the center image upwards
|
---|
109 | void
|
---|
110 | ViewVectors(const Vector3 &DirAt, const Vector3 &Viewer,
|
---|
111 | const Vector3 &UpL, Vector3 &ViewV, Vector3 &ViewU,
|
---|
112 | Vector3 &ViewN)
|
---|
113 | {
|
---|
114 | Vector3 U, V, N;
|
---|
115 | Vector3 Up = Normalize(UpL);
|
---|
116 |
|
---|
117 | N = -Normalize(DirAt);
|
---|
118 |
|
---|
119 | V = Normalize(Up - DirAt);
|
---|
120 | V -= N * DotProd(V, N);
|
---|
121 | V = Normalize(V);
|
---|
122 | U = CrossProd(V, N);
|
---|
123 |
|
---|
124 | ViewU = U; // rightwards
|
---|
125 | ViewV = V; // upwards
|
---|
126 | ViewN = -N; // forwards
|
---|
127 | #ifdef GTP_DEBUG
|
---|
128 | const float eps = 1e-3f;
|
---|
129 | if (fabs(Magnitude(ViewU) - 1.0) > eps) {
|
---|
130 | Debug << "ViewU magnitude error= " << Magnitude(ViewU) << "\n";
|
---|
131 | }
|
---|
132 | if (fabs(Magnitude(ViewV) - 1.0) > eps) {
|
---|
133 | Debug << "ViewU magnitude error= " << Magnitude(ViewV) << "\n";
|
---|
134 | }
|
---|
135 | if (fabs(Magnitude(ViewN) - 1.0) > eps) {
|
---|
136 | Debug << "ViewU magnitude error= " << Magnitude(ViewN) << "\n";
|
---|
137 | }
|
---|
138 | #endif // GTP_DEBUG
|
---|
139 |
|
---|
140 | return;
|
---|
141 | }
|
---|
142 |
|
---|
143 | // Given the intersection point `P', you have available normal `N'
|
---|
144 | // of unit length. Let us suppose the incoming ray has direction `D'.
|
---|
145 | // Then we can construct such two vectors `U' and `V' that
|
---|
146 | // `U',`N', and `D' are coplanar, and `V' is perpendicular
|
---|
147 | // to the vectors `N','D', and `V'. Then 'N', 'U', and 'V' create
|
---|
148 | // the orthonormal base in space R3.
|
---|
149 | void
|
---|
150 | TangentVectors(Vector3 &U,
|
---|
151 | Vector3 &V, // output
|
---|
152 | const Vector3 &normal, // input
|
---|
153 | const Vector3 &dirIncoming)
|
---|
154 | {
|
---|
155 | #ifdef GTP_DEBUG
|
---|
156 | float d = Magnitude(normal);
|
---|
157 | if ( (d < 0.99) ||
|
---|
158 | (d > 1.01) ) {
|
---|
159 | Debug << " The normal has not unit length = " << d << endl;
|
---|
160 | }
|
---|
161 | d = Magnitude(dirIncoming);
|
---|
162 | if ( (d < 0.99) ||
|
---|
163 | (d > 1.01) ) {
|
---|
164 | Debug << " The incoming dir has not unit length = " << d << endl;
|
---|
165 | }
|
---|
166 | #endif
|
---|
167 |
|
---|
168 | V = CrossProd(normal, dirIncoming);
|
---|
169 |
|
---|
170 | if (SqrMagnitude(V) < 1e-3) {
|
---|
171 | // the normal and dirIncoming are colinear
|
---|
172 | // we can/have to generate arbitrary perpendicular vector to normal.
|
---|
173 | if (fabs(normal.x) < 0.6)
|
---|
174 | V.SetValue(0.0, -normal.z, normal.y);
|
---|
175 | else {
|
---|
176 | if (fabs(normal.y) < 0.6)
|
---|
177 | V.SetValue(-normal.z, 0.0, normal.x);
|
---|
178 | else
|
---|
179 | V.SetValue(-normal.y, normal.x, 0.0);
|
---|
180 | }
|
---|
181 | }
|
---|
182 | V = Normalize(V);
|
---|
183 |
|
---|
184 | U = CrossProd(normal, V);
|
---|
185 | #ifdef GTP_DEBUG
|
---|
186 | d = SqrMagnitude(U);
|
---|
187 | if ( (d < 0.99) ||
|
---|
188 | (d > 1.01) ) {
|
---|
189 | Debug << "The size of U vector incorrect\n";
|
---|
190 | }
|
---|
191 | #endif
|
---|
192 | return;
|
---|
193 | }
|
---|
194 |
|
---|
195 | #define USE_HALTON 0
|
---|
196 |
|
---|
197 |
|
---|
198 | Vector3
|
---|
199 | UniformRandomVector()
|
---|
200 | {
|
---|
201 | // float r1 = RandomValue(0.0f, 1.0f);
|
---|
202 | // float r2 = RandomValue(0.0f, 1.0f);
|
---|
203 | float r1, r2;
|
---|
204 |
|
---|
205 | #if USE_HALTON
|
---|
206 | halton2.GetNext(r1, r2);
|
---|
207 | #else
|
---|
208 | r1 = RandomValue(0.0f, 1.0f);
|
---|
209 | r2 = RandomValue(0.0f, 1.0f);
|
---|
210 | #endif
|
---|
211 |
|
---|
212 | return UniformRandomVector(r1, r2);
|
---|
213 | }
|
---|
214 |
|
---|
215 | Vector3
|
---|
216 | UniformRandomVector(const float r1, const float r2)
|
---|
217 | {
|
---|
218 | float cosTheta = 1.0f - 2*r1;
|
---|
219 | float sinTheta = sqrt(1 - sqr(cosTheta));
|
---|
220 | float fi = 2.0f*M_PI*r2;
|
---|
221 |
|
---|
222 | Vector3 dir(sinTheta*sin(fi),
|
---|
223 | cosTheta,
|
---|
224 | sinTheta*cos(fi));
|
---|
225 |
|
---|
226 | return dir;
|
---|
227 | }
|
---|
228 |
|
---|
229 | Vector3
|
---|
230 | CosineRandomVector(const Vector3 &normal)
|
---|
231 | {
|
---|
232 | float r1, r2;
|
---|
233 |
|
---|
234 | r1 = RandomValue(0.0f, 1.0f);
|
---|
235 | r2 = RandomValue(0.0f, 1.0f);
|
---|
236 |
|
---|
237 | return CosineRandomVector(r1, r2, normal);
|
---|
238 | }
|
---|
239 |
|
---|
240 | Vector3
|
---|
241 | CosineRandomVector(const float r1,
|
---|
242 | const float r2,
|
---|
243 | const Vector3 &normal)
|
---|
244 | {
|
---|
245 | float theta = 2.0f*M_PI*r2;
|
---|
246 | float radius = sqrt(r1);
|
---|
247 | float x = radius*sin(theta);
|
---|
248 | float z = radius*cos(theta);
|
---|
249 | float y = sqrt(1.0f - x*x - z*z);
|
---|
250 |
|
---|
251 | Vector3 dir(x,
|
---|
252 | y,
|
---|
253 | z);
|
---|
254 |
|
---|
255 | // return Normalize(dir);
|
---|
256 |
|
---|
257 | Matrix4x4 m = RotationVectorsMatrix(
|
---|
258 | normal,
|
---|
259 | Vector3(0,1,0)
|
---|
260 | );
|
---|
261 | Matrix4x4 mi = Invert(m);
|
---|
262 | m = m*RotationVectorsMatrix(
|
---|
263 | Vector3(0,1,0),
|
---|
264 | Normalize(dir))*mi;
|
---|
265 |
|
---|
266 | return TransformNormal(m, normal);
|
---|
267 | }
|
---|
268 |
|
---|
269 | Vector3
|
---|
270 | UniformRandomVector(const Vector3 &normal)
|
---|
271 | {
|
---|
272 | // float r1 = RandomValue(0.0f, 1.0f);
|
---|
273 | // float r2 = RandomValue(0.0f, 1.0f);
|
---|
274 | float r1, r2;
|
---|
275 |
|
---|
276 | #if USE_HALTON
|
---|
277 | halton2.GetNext(r1, r2);
|
---|
278 | #else
|
---|
279 | r1 = RandomValue(0.0f, 1.0f);
|
---|
280 | r2 = RandomValue(0.0f, 1.0f);
|
---|
281 | #endif
|
---|
282 |
|
---|
283 |
|
---|
284 | float cosTheta = 1.0f - r1;
|
---|
285 | float sinTheta = sqrt(1 - sqr(cosTheta));
|
---|
286 | float fi = 2.0f*M_PI*r2;
|
---|
287 |
|
---|
288 | Vector3 dir(sinTheta*sin(fi),
|
---|
289 | cosTheta,
|
---|
290 | sinTheta*cos(fi));
|
---|
291 |
|
---|
292 |
|
---|
293 | // return Normalize(dir);
|
---|
294 |
|
---|
295 | Matrix4x4 m = RotationVectorsMatrix(
|
---|
296 | normal,
|
---|
297 | Vector3(0,1,0));
|
---|
298 | Matrix4x4 mi = Invert(m);
|
---|
299 | m = m*RotationVectorsMatrix(
|
---|
300 | Vector3(0,1,0),
|
---|
301 | Normalize(dir))*mi;
|
---|
302 |
|
---|
303 | return TransformNormal(m, normal);
|
---|
304 |
|
---|
305 | // return TransformNormal(
|
---|
306 | // RotationVectorsMatrix(
|
---|
307 | // Vector3(0,1,0),
|
---|
308 | // Normalize(dir)
|
---|
309 | // ),
|
---|
310 | // normal
|
---|
311 | // );
|
---|
312 | }
|
---|
313 |
|
---|
314 |
|
---|
315 | bool Vector3::CheckValidity() const
|
---|
316 | {
|
---|
317 | #ifdef _WIN32
|
---|
318 | return !(_isnan(x) || _isnan(y) || _isnan(z));
|
---|
319 | #else
|
---|
320 | return !(isnan(x) || isnan(y) || isnan(z));
|
---|
321 | #endif
|
---|
322 | //return ((x != x) || (y != y) || (z != z));
|
---|
323 | }
|
---|
324 |
|
---|
325 | }
|
---|