[176] | 1 | #include "Matrix4x4.h"
|
---|
[162] | 2 | #include "Vector3.h"
|
---|
[354] | 3 | #include "Halton.h"
|
---|
[162] | 4 |
|
---|
[349] | 5 | float Vector3::sDistTolerance = 0.002f; |
---|
| 6 | float Vector3::sDistToleranceSqrt = 0.000004f;
|
---|
[162] | 7 |
|
---|
| 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
|
---|
| 127 | #ifdef _DEBUG
|
---|
| 128 | const float eps = 1e-3;
|
---|
| 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 | }
|
---|
| 138 | #endif // _DEBUG
|
---|
| 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 | {
|
---|
| 155 | #ifdef _DEBUG
|
---|
| 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);
|
---|
| 185 | #ifdef _DEBUG
|
---|
| 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 |
|
---|
| 195 | Vector3
|
---|
| 196 | UniformRandomVector(const Vector3 &normal)
|
---|
| 197 | {
|
---|
[354] | 198 | // float r1 = RandomValue(0.0f, 1.0f);
|
---|
| 199 | // float r2 = RandomValue(0.0f, 1.0f);
|
---|
| 200 | float r1, r2;
|
---|
| 201 |
|
---|
| 202 | halton2.GetNext(r1, r2);
|
---|
| 203 |
|
---|
| 204 |
|
---|
| 205 | float cosTheta = 1.0f - r1;
|
---|
[176] | 206 | float sinTheta = sqrt(1 - sqr(cosTheta));
|
---|
| 207 | float fi = 2.0f*M_PI*r2;
|
---|
| 208 |
|
---|
| 209 | Vector3 dir(sinTheta*sin(fi),
|
---|
[340] | 210 | cosTheta,
|
---|
| 211 | sinTheta*cos(fi));
|
---|
| 212 |
|
---|
[176] | 213 |
|
---|
| 214 | // return Normalize(dir);
|
---|
[340] | 215 |
|
---|
[176] | 216 | Matrix4x4 m = RotationVectorsMatrix(
|
---|
[340] | 217 | normal,
|
---|
| 218 | Vector3(0,1,0));
|
---|
[176] | 219 | Matrix4x4 mi = Invert(m);
|
---|
| 220 | m = m*RotationVectorsMatrix(
|
---|
[340] | 221 | Vector3(0,1,0),
|
---|
| 222 | Normalize(dir))*mi;
|
---|
[176] | 223 |
|
---|
| 224 | return TransformNormal(m, normal);
|
---|
| 225 |
|
---|
[340] | 226 | // return TransformNormal(
|
---|
| 227 | // RotationVectorsMatrix(
|
---|
| 228 | // Vector3(0,1,0),
|
---|
| 229 | // Normalize(dir)
|
---|
| 230 | // ),
|
---|
| 231 | // normal
|
---|
| 232 | // );
|
---|
[176] | 233 | }
|
---|