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

Revision 1420, 6.1 KB checked in by mattausch, 18 years ago (diff)

corrected raycasting bug for triangles because of ill defined triangles

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