00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026 #include "stdafx.h"
00027 #include "transform.h"
00028 #include "shape.h"
00029
00030
00031 bool SolveLinearSystem2x2(const float A[2][2],
00032 const float B[2], float *x0, float *x1) {
00033 float det = A[0][0]*A[1][1] - A[0][1]*A[1][0];
00034 if (fabsf(det) < 1e-10f)
00035 return false;
00036 *x0 = (A[1][1]*B[0] - A[0][1]*B[1]) / det;
00037 *x1 = (A[0][0]*B[1] - A[1][0]*B[0]) / det;
00038 if (isnan(*x0) || isnan(*x1))
00039 return false;
00040 return true;
00041 }
00042
00043
00044 Matrix4x4::Matrix4x4(float mat[4][4]) {
00045 memcpy(m, mat, 16*sizeof(float));
00046 }
00047
00048
00049 Matrix4x4::Matrix4x4(float t00, float t01, float t02, float t03,
00050 float t10, float t11, float t12, float t13,
00051 float t20, float t21, float t22, float t23,
00052 float t30, float t31, float t32, float t33) {
00053 m[0][0] = t00; m[0][1] = t01; m[0][2] = t02; m[0][3] = t03;
00054 m[1][0] = t10; m[1][1] = t11; m[1][2] = t12; m[1][3] = t13;
00055 m[2][0] = t20; m[2][1] = t21; m[2][2] = t22; m[2][3] = t23;
00056 m[3][0] = t30; m[3][1] = t31; m[3][2] = t32; m[3][3] = t33;
00057 }
00058
00059
00060 Matrix4x4 Transpose(const Matrix4x4 &m) {
00061 return Matrix4x4(m.m[0][0], m.m[1][0], m.m[2][0], m.m[3][0],
00062 m.m[0][1], m.m[1][1], m.m[2][1], m.m[3][1],
00063 m.m[0][2], m.m[1][2], m.m[2][2], m.m[3][2],
00064 m.m[0][3], m.m[1][3], m.m[2][3], m.m[3][3]);
00065 }
00066
00067
00068 Matrix4x4 Inverse(const Matrix4x4 &m) {
00069 int indxc[4], indxr[4];
00070 int ipiv[4] = { 0, 0, 0, 0 };
00071 float minv[4][4];
00072 memcpy(minv, m.m, 4*4*sizeof(float));
00073 for (int i = 0; i < 4; i++) {
00074 int irow = -1, icol = -1;
00075 float big = 0.;
00076
00077 for (int j = 0; j < 4; j++) {
00078 if (ipiv[j] != 1) {
00079 for (int k = 0; k < 4; k++) {
00080 if (ipiv[k] == 0) {
00081 if (fabsf(minv[j][k]) >= big) {
00082 big = float(fabsf(minv[j][k]));
00083 irow = j;
00084 icol = k;
00085 }
00086 }
00087 else if (ipiv[k] > 1)
00088 Error("Singular matrix in MatrixInvert");
00089 }
00090 }
00091 }
00092 ++ipiv[icol];
00093
00094 if (irow != icol) {
00095 for (int k = 0; k < 4; ++k)
00096 swap(minv[irow][k], minv[icol][k]);
00097 }
00098 indxr[i] = irow;
00099 indxc[i] = icol;
00100 if (minv[icol][icol] == 0.)
00101 Error("Singular matrix in MatrixInvert");
00102
00103
00104 float pivinv = 1.f / minv[icol][icol];
00105 minv[icol][icol] = 1.f;
00106 for (int j = 0; j < 4; j++)
00107 minv[icol][j] *= pivinv;
00108
00109
00110 for (int j = 0; j < 4; j++) {
00111 if (j != icol) {
00112 float save = minv[j][icol];
00113 minv[j][icol] = 0;
00114 for (int k = 0; k < 4; k++)
00115 minv[j][k] -= minv[icol][k]*save;
00116 }
00117 }
00118 }
00119
00120 for (int j = 3; j >= 0; j--) {
00121 if (indxr[j] != indxc[j]) {
00122 for (int k = 0; k < 4; k++)
00123 swap(minv[k][indxr[j]], minv[k][indxc[j]]);
00124 }
00125 }
00126 return Matrix4x4(minv);
00127 }
00128
00129
00130
00131
00132 void Transform::Print(FILE *f) const {
00133 m.Print(f);
00134 }
00135
00136
00137 Transform Translate(const Vector &delta) {
00138 Matrix4x4 m(1, 0, 0, delta.x,
00139 0, 1, 0, delta.y,
00140 0, 0, 1, delta.z,
00141 0, 0, 0, 1);
00142 Matrix4x4 minv(1, 0, 0, -delta.x,
00143 0, 1, 0, -delta.y,
00144 0, 0, 1, -delta.z,
00145 0, 0, 0, 1);
00146 return Transform(m, minv);
00147 }
00148
00149
00150 Transform Scale(float x, float y, float z) {
00151 Matrix4x4 m(x, 0, 0, 0,
00152 0, y, 0, 0,
00153 0, 0, z, 0,
00154 0, 0, 0, 1);
00155 Matrix4x4 minv(1.f/x, 0, 0, 0,
00156 0, 1.f/y, 0, 0,
00157 0, 0, 1.f/z, 0,
00158 0, 0, 0, 1);
00159 return Transform(m, minv);
00160 }
00161
00162
00163 Transform RotateX(float angle) {
00164 float sin_t = sinf(Radians(angle));
00165 float cos_t = cosf(Radians(angle));
00166 Matrix4x4 m(1, 0, 0, 0,
00167 0, cos_t, -sin_t, 0,
00168 0, sin_t, cos_t, 0,
00169 0, 0, 0, 1);
00170 return Transform(m, Transpose(m));
00171 }
00172
00173
00174 Transform RotateY(float angle) {
00175 float sin_t = sinf(Radians(angle));
00176 float cos_t = cosf(Radians(angle));
00177 Matrix4x4 m( cos_t, 0, sin_t, 0,
00178 0, 1, 0, 0,
00179 -sin_t, 0, cos_t, 0,
00180 0, 0, 0, 1);
00181 return Transform(m, Transpose(m));
00182 }
00183
00184
00185
00186 Transform RotateZ(float angle) {
00187 float sin_t = sinf(Radians(angle));
00188 float cos_t = cosf(Radians(angle));
00189 Matrix4x4 m(cos_t, -sin_t, 0, 0,
00190 sin_t, cos_t, 0, 0,
00191 0, 0, 1, 0,
00192 0, 0, 0, 1);
00193 return Transform(m, Transpose(m));
00194 }
00195
00196
00197 Transform Rotate(float angle, const Vector &axis) {
00198 Vector a = Normalize(axis);
00199 float s = sinf(Radians(angle));
00200 float c = cosf(Radians(angle));
00201 float m[4][4];
00202
00203 m[0][0] = a.x * a.x + (1.f - a.x * a.x) * c;
00204 m[0][1] = a.x * a.y * (1.f - c) - a.z * s;
00205 m[0][2] = a.x * a.z * (1.f - c) + a.y * s;
00206 m[0][3] = 0;
00207
00208 m[1][0] = a.x * a.y * (1.f - c) + a.z * s;
00209 m[1][1] = a.y * a.y + (1.f - a.y * a.y) * c;
00210 m[1][2] = a.y * a.z * (1.f - c) - a.x * s;
00211 m[1][3] = 0;
00212
00213 m[2][0] = a.x * a.z * (1.f - c) - a.y * s;
00214 m[2][1] = a.y * a.z * (1.f - c) + a.x * s;
00215 m[2][2] = a.z * a.z + (1.f - a.z * a.z) * c;
00216 m[2][3] = 0;
00217
00218 m[3][0] = 0;
00219 m[3][1] = 0;
00220 m[3][2] = 0;
00221 m[3][3] = 1;
00222
00223 Matrix4x4 mat(m);
00224 return Transform(mat, Transpose(mat));
00225 }
00226
00227
00228 Transform LookAt(const Point &pos, const Point &look, const Vector &up) {
00229 float m[4][4];
00230
00231 m[0][3] = pos.x;
00232 m[1][3] = pos.y;
00233 m[2][3] = pos.z;
00234 m[3][3] = 1;
00235
00236
00237 Vector dir = Normalize(look - pos);
00238 Vector left = Normalize(Cross(Normalize(up), dir));
00239 Vector newUp = Cross(dir, left);
00240 m[0][0] = left.x;
00241 m[1][0] = left.y;
00242 m[2][0] = left.z;
00243 m[3][0] = 0.;
00244 m[0][1] = newUp.x;
00245 m[1][1] = newUp.y;
00246 m[2][1] = newUp.z;
00247 m[3][1] = 0.;
00248 m[0][2] = dir.x;
00249 m[1][2] = dir.y;
00250 m[2][2] = dir.z;
00251 m[3][2] = 0.;
00252 Matrix4x4 camToWorld(m);
00253 return Transform(Inverse(camToWorld), camToWorld);
00254 }
00255
00256
00257 BBox Transform::operator()(const BBox &b) const {
00258 const Transform &M = *this;
00259 BBox ret( M(Point(b.pMin.x, b.pMin.y, b.pMin.z)));
00260 ret = Union(ret, M(Point(b.pMax.x, b.pMin.y, b.pMin.z)));
00261 ret = Union(ret, M(Point(b.pMin.x, b.pMax.y, b.pMin.z)));
00262 ret = Union(ret, M(Point(b.pMin.x, b.pMin.y, b.pMax.z)));
00263 ret = Union(ret, M(Point(b.pMin.x, b.pMax.y, b.pMax.z)));
00264 ret = Union(ret, M(Point(b.pMax.x, b.pMax.y, b.pMin.z)));
00265 ret = Union(ret, M(Point(b.pMax.x, b.pMin.y, b.pMax.z)));
00266 ret = Union(ret, M(Point(b.pMax.x, b.pMax.y, b.pMax.z)));
00267 return ret;
00268 }
00269
00270
00271 Transform Transform::operator*(const Transform &t2) const {
00272 Matrix4x4 m1 = Matrix4x4::Mul(m, t2.m);
00273 Matrix4x4 m2 = Matrix4x4::Mul(t2.mInv, mInv);
00274 return Transform(m1, m2);
00275 }
00276
00277
00278 bool Transform::SwapsHandedness() const {
00279 float det = ((m.m[0][0] *
00280 (m.m[1][1] * m.m[2][2] -
00281 m.m[1][2] * m.m[2][1])) -
00282 (m.m[0][1] *
00283 (m.m[1][0] * m.m[2][2] -
00284 m.m[1][2] * m.m[2][0])) +
00285 (m.m[0][2] *
00286 (m.m[1][0] * m.m[2][1] -
00287 m.m[1][1] * m.m[2][0])));
00288 return det < 0.f;
00289 }
00290
00291
00292 Transform Orthographic(float znear, float zfar) {
00293 return Scale(1.f, 1.f, 1.f / (zfar-znear)) *
00294 Translate(Vector(0.f, 0.f, -znear));
00295 }
00296
00297
00298 Transform Perspective(float fov, float n, float f) {
00299
00300 Matrix4x4 persp = Matrix4x4(1, 0, 0, 0,
00301 0, 1, 0, 0,
00302 0, 0, f / (f - n), -f*n / (f - n),
00303 0, 0, 1, 0);
00304
00305
00306 float invTanAng = 1.f / tanf(Radians(fov) / 2.f);
00307 return Scale(invTanAng, invTanAng, 1) * Transform(persp);
00308 }
00309
00310
00311
00312
00313 void AnimatedTransform::Decompose(const Matrix4x4 &m, Vector *T,
00314 Quaternion *Rquat, Matrix4x4 *S) {
00315
00316 T->x = m.m[0][3];
00317 T->y = m.m[1][3];
00318 T->z = m.m[2][3];
00319
00320
00321 Matrix4x4 M = m;
00322 for (int i = 0; i < 3; ++i)
00323 M.m[i][3] = M.m[3][i] = 0.f;
00324 M.m[3][3] = 1.f;
00325
00326
00327 float norm;
00328 int count = 0;
00329 Matrix4x4 R = M;
00330 do {
00331
00332 Matrix4x4 Rnext;
00333 Matrix4x4 Rit = Inverse(Transpose(R));
00334 for (int i = 0; i < 4; ++i)
00335 for (int j = 0; j < 4; ++j)
00336 Rnext.m[i][j] = 0.5f * (R.m[i][j] + Rit.m[i][j]);
00337
00338
00339 norm = 0.f;
00340 for (int i = 0; i < 3; ++i) {
00341 float n = fabsf(R.m[i][0] - Rnext.m[i][0]) +
00342 fabsf(R.m[i][1] - Rnext.m[i][1]) +
00343 fabsf(R.m[i][2] - Rnext.m[i][2]);
00344 norm = max(norm, n);
00345 }
00346 R = Rnext;
00347 } while (++count < 100 && norm > .0001f);
00348
00349 *Rquat = Quaternion(R);
00350
00351
00352 *S = Matrix4x4::Mul(Inverse(R), M);
00353 }
00354
00355
00356 void AnimatedTransform::Interpolate(float time, Transform *t) const {
00357
00358 if (!actuallyAnimated || time <= startTime) {
00359 *t = *startTransform;
00360 return;
00361 }
00362 if (time >= endTime) {
00363 *t = *endTransform;
00364 return;
00365 }
00366 float dt = (time - startTime) / (endTime - startTime);
00367
00368 Vector trans = (1.f - dt) * T[0] + dt * T[1];
00369
00370
00371 Quaternion rotate = Slerp(dt, R[0], R[1]);
00372
00373
00374 Matrix4x4 scale;
00375 for (int i = 0; i < 3; ++i)
00376 for (int j = 0; j < 3; ++j)
00377 scale.m[i][j] = Lerp(dt, S[0].m[i][j], S[1].m[i][j]);
00378
00379
00380 *t = Translate(trans) * rotate.ToTransform() * Transform(scale);
00381 }
00382
00383
00384 BBox AnimatedTransform::MotionBounds(const BBox &b,
00385 bool useInverse) const {
00386 if (!actuallyAnimated) return Inverse(*startTransform)(b);
00387 BBox ret;
00388 const int nSteps = 128;
00389 for (int i = 0; i < nSteps; ++i) {
00390 Transform t;
00391 float time = Lerp(float(i)/float(nSteps-1), startTime, endTime);
00392 Interpolate(time, &t);
00393 if (useInverse) t = Inverse(t);
00394 ret = Union(ret, t(b));
00395 }
00396 return ret;
00397 }
00398
00399
00400 void AnimatedTransform::operator()(const Ray &r, Ray *tr) const {
00401 if (!actuallyAnimated || r.time <= startTime)
00402 (*startTransform)(r, tr);
00403 else if (r.time >= endTime)
00404 (*endTransform)(r, tr);
00405 else {
00406 Transform t;
00407 Interpolate(r.time, &t);
00408 t(r, tr);
00409 }
00410 tr->time = r.time;
00411 }
00412
00413
00414 void AnimatedTransform::operator()(const RayDifferential &r,
00415 RayDifferential *tr) const {
00416 if (!actuallyAnimated || r.time <= startTime)
00417 (*startTransform)(r, tr);
00418 else if (r.time >= endTime)
00419 (*endTransform)(r, tr);
00420 else {
00421 Transform t;
00422 Interpolate(r.time, &t);
00423 t(r, tr);
00424 }
00425 tr->time = r.time;
00426 }
00427
00428
00429 Point AnimatedTransform::operator()(float time, const Point &p) const {
00430 if (!actuallyAnimated || time <= startTime)
00431 return (*startTransform)(p);
00432 else if (time >= endTime)
00433 return (*endTransform)(p);
00434 Transform t;
00435 Interpolate(time, &t);
00436 return t(p);
00437 }
00438
00439
00440 Vector AnimatedTransform::operator()(float time, const Vector &v) const {
00441 if (!actuallyAnimated || time <= startTime)
00442 return (*startTransform)(v);
00443 else if (time >= endTime)
00444 return (*endTransform)(v);
00445 Transform t;
00446 Interpolate(time, &t);
00447 return t(v);
00448 }
00449
00450
00451 Ray AnimatedTransform::operator()(const Ray &r) const {
00452 Ray ret;
00453 (*this)(r, &ret);
00454 return ret;
00455 }
00456
00457