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