source: GTP/trunk/Lib/Vis/Preprocessing/src/Vector3.cpp @ 2575

Revision 2575, 7.3 KB checked in by bittner, 17 years ago (diff)

big merge: preparation for havran ray caster, check if everything works

Line 
1#include "Matrix4x4.h"
2#include "Vector3.h"
3#include "Halton.h"
4#include <float.h>
5
6namespace 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.
12void
13Minimize(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.
27void
28Maximize(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
45void
46Vector3::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
65int
66Vector3::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
83int
84Vector3::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
109void
110ViewVectors(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.
149void
150TangentVectors(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
198Vector3
199UniformRandomVector()
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
215Vector3
216UniformRandomVector(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
229Vector3
230CosineRandomVector(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
240Vector3
241CosineRandomVector(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
269Vector3
270UniformRandomVector(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
315bool 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}
Note: See TracBrowser for help on using the repository browser.