[176] | 1 | #include "Matrix4x4.h"
|
---|
[162] | 2 | #include "Vector3.h"
|
---|
[354] | 3 | #include "Halton.h"
|
---|
[2575] | 4 | #include <float.h>
|
---|
[162] | 5 |
|
---|
[863] | 6 | namespace GtpVisibilityPreprocessor {
|
---|
[860] | 7 |
|
---|
[162] | 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
|
---|
[1715] | 127 | #ifdef GTP_DEBUG
|
---|
[728] | 128 | const float eps = 1e-3f;
|
---|
[162] | 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 | }
|
---|
[1715] | 138 | #endif // GTP_DEBUG
|
---|
[162] | 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 | {
|
---|
[1715] | 155 | #ifdef GTP_DEBUG
|
---|
[162] | 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);
|
---|
[1715] | 185 | #ifdef GTP_DEBUG
|
---|
[162] | 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 | }
|
---|
[176] | 194 |
|
---|
[1715] | 195 | #define USE_HALTON 0
|
---|
| 196 |
|
---|
[1867] | 197 |
|
---|
[176] | 198 | Vector3
|
---|
[492] | 199 | UniformRandomVector()
|
---|
[176] | 200 | {
|
---|
[492] | 201 | // float r1 = RandomValue(0.0f, 1.0f);
|
---|
| 202 | // float r2 = RandomValue(0.0f, 1.0f);
|
---|
| 203 | float r1, r2;
|
---|
[1715] | 204 |
|
---|
| 205 | #if USE_HALTON
|
---|
[492] | 206 | halton2.GetNext(r1, r2);
|
---|
[1715] | 207 | #else
|
---|
| 208 | r1 = RandomValue(0.0f, 1.0f);
|
---|
| 209 | r2 = RandomValue(0.0f, 1.0f);
|
---|
| 210 | #endif
|
---|
[1867] | 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;
|
---|
[1715] | 221 |
|
---|
[1867] | 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;
|
---|
[176] | 233 |
|
---|
[1883] | 234 | r1 = RandomValue(0.0f, 1.0f);
|
---|
| 235 | r2 = RandomValue(0.0f, 1.0f);
|
---|
| 236 |
|
---|
| 237 | return CosineRandomVector(r1, r2, normal);
|
---|
| 238 | }
|
---|
[1867] | 239 |
|
---|
[1883] | 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);
|
---|
[1867] | 256 |
|
---|
[1883] | 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;
|
---|
[176] | 265 |
|
---|
[1883] | 266 | return TransformNormal(m, normal);
|
---|
[492] | 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 |
|
---|
[1715] | 276 | #if USE_HALTON
|
---|
[492] | 277 | halton2.GetNext(r1, r2);
|
---|
[1715] | 278 | #else
|
---|
| 279 | r1 = RandomValue(0.0f, 1.0f);
|
---|
| 280 | r2 = RandomValue(0.0f, 1.0f);
|
---|
| 281 | #endif
|
---|
[492] | 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 |
|
---|
[176] | 293 | // return Normalize(dir);
|
---|
[492] | 294 |
|
---|
[176] | 295 | Matrix4x4 m = RotationVectorsMatrix(
|
---|
[492] | 296 | normal,
|
---|
| 297 | Vector3(0,1,0));
|
---|
[176] | 298 | Matrix4x4 mi = Invert(m);
|
---|
| 299 | m = m*RotationVectorsMatrix(
|
---|
[492] | 300 | Vector3(0,1,0),
|
---|
| 301 | Normalize(dir))*mi;
|
---|
[176] | 302 |
|
---|
| 303 | return TransformNormal(m, normal);
|
---|
| 304 |
|
---|
[492] | 305 | // return TransformNormal(
|
---|
| 306 | // RotationVectorsMatrix(
|
---|
| 307 | // Vector3(0,1,0),
|
---|
| 308 | // Normalize(dir)
|
---|
| 309 | // ),
|
---|
| 310 | // normal
|
---|
| 311 | // );
|
---|
[176] | 312 | }
|
---|
[492] | 313 |
|
---|
[1420] | 314 |
|
---|
| 315 | bool Vector3::CheckValidity() const
|
---|
| 316 | {
|
---|
[2575] | 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));
|
---|
[1420] | 323 | }
|
---|
| 324 |
|
---|
[1715] | 325 | }
|
---|