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

Revision 492, 5.9 KB checked in by bittner, 19 years ago (diff)

Large merge - viewcells seem not functional now

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