[176] | 1 | #include "Matrix4x4.h"
|
---|
[162] | 2 | #include "Vector3.h"
|
---|
[354] | 3 | #include "Halton.h"
|
---|
[1420] | 4 | #include "float.h"
|
---|
[162] | 5 |
|
---|
[1420] | 6 |
|
---|
[863] | 7 | namespace GtpVisibilityPreprocessor {
|
---|
[860] | 8 |
|
---|
[162] | 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.
|
---|
| 13 | void
|
---|
| 14 | Minimize(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.
|
---|
| 28 | void
|
---|
| 29 | Maximize(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 |
|
---|
| 46 | void
|
---|
| 47 | Vector3::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
|
---|
| 66 | int
|
---|
| 67 | Vector3::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
|
---|
| 84 | int
|
---|
| 85 | Vector3::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
|
---|
| 110 | void
|
---|
| 111 | ViewVectors(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
|
---|
[1715] | 128 | #ifdef GTP_DEBUG
|
---|
[728] | 129 | const float eps = 1e-3f;
|
---|
[162] | 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 | }
|
---|
[1715] | 139 | #endif // GTP_DEBUG
|
---|
[162] | 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.
|
---|
| 150 | void
|
---|
| 151 | TangentVectors(Vector3 &U,
|
---|
| 152 | Vector3 &V, // output
|
---|
| 153 | const Vector3 &normal, // input
|
---|
| 154 | const Vector3 &dirIncoming)
|
---|
| 155 | {
|
---|
[1715] | 156 | #ifdef GTP_DEBUG
|
---|
[162] | 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);
|
---|
[1715] | 186 | #ifdef GTP_DEBUG
|
---|
[162] | 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 | }
|
---|
[176] | 195 |
|
---|
[1715] | 196 | #define USE_HALTON 0
|
---|
| 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
|
---|
| 211 |
|
---|
| 212 | float cosTheta = 1.0f - 2*r1;
|
---|
[176] | 213 | float sinTheta = sqrt(1 - sqr(cosTheta));
|
---|
| 214 | float fi = 2.0f*M_PI*r2;
|
---|
| 215 |
|
---|
| 216 | Vector3 dir(sinTheta*sin(fi),
|
---|
[492] | 217 | cosTheta,
|
---|
| 218 | sinTheta*cos(fi));
|
---|
[176] | 219 |
|
---|
[492] | 220 | return dir;
|
---|
| 221 | }
|
---|
| 222 |
|
---|
| 223 | Vector3
|
---|
| 224 | UniformRandomVector(const Vector3 &normal)
|
---|
| 225 | {
|
---|
| 226 | // float r1 = RandomValue(0.0f, 1.0f);
|
---|
| 227 | // float r2 = RandomValue(0.0f, 1.0f);
|
---|
| 228 | float r1, r2;
|
---|
| 229 |
|
---|
[1715] | 230 | #if USE_HALTON
|
---|
[492] | 231 | halton2.GetNext(r1, r2);
|
---|
[1715] | 232 | #else
|
---|
| 233 | r1 = RandomValue(0.0f, 1.0f);
|
---|
| 234 | r2 = RandomValue(0.0f, 1.0f);
|
---|
| 235 | #endif
|
---|
[492] | 236 |
|
---|
| 237 |
|
---|
| 238 | float cosTheta = 1.0f - r1;
|
---|
| 239 | float sinTheta = sqrt(1 - sqr(cosTheta));
|
---|
| 240 | float fi = 2.0f*M_PI*r2;
|
---|
| 241 |
|
---|
| 242 | Vector3 dir(sinTheta*sin(fi),
|
---|
| 243 | cosTheta,
|
---|
| 244 | sinTheta*cos(fi));
|
---|
| 245 |
|
---|
| 246 |
|
---|
[176] | 247 | // return Normalize(dir);
|
---|
[492] | 248 |
|
---|
[176] | 249 | Matrix4x4 m = RotationVectorsMatrix(
|
---|
[492] | 250 | normal,
|
---|
| 251 | Vector3(0,1,0));
|
---|
[176] | 252 | Matrix4x4 mi = Invert(m);
|
---|
| 253 | m = m*RotationVectorsMatrix(
|
---|
[492] | 254 | Vector3(0,1,0),
|
---|
| 255 | Normalize(dir))*mi;
|
---|
[176] | 256 |
|
---|
| 257 | return TransformNormal(m, normal);
|
---|
| 258 |
|
---|
[492] | 259 | // return TransformNormal(
|
---|
| 260 | // RotationVectorsMatrix(
|
---|
| 261 | // Vector3(0,1,0),
|
---|
| 262 | // Normalize(dir)
|
---|
| 263 | // ),
|
---|
| 264 | // normal
|
---|
| 265 | // );
|
---|
[176] | 266 | }
|
---|
[492] | 267 |
|
---|
[1420] | 268 |
|
---|
| 269 | bool Vector3::CheckValidity() const
|
---|
| 270 | {
|
---|
| 271 | return !(_isnan(x) || _isnan(y) || _isnan(z));
|
---|
| 272 | //return ((x != x) || (y != y) || (z != z));
|
---|
| 273 | }
|
---|
| 274 |
|
---|
[1715] | 275 | }
|
---|