text: 
  /*
      pbrt source code Copyright(c) 1998-2012 Matt Pharr and Greg Humphreys.
  
      This file is part of pbrt.
  
      Redistribution and use in source and binary forms, with or without
      modification, are permitted provided that the following conditions are
      met:
  
      - Redistributions of source code must retain the above copyright
        notice, this list of conditions and the following disclaimer.
  
      - Redistributions in binary form must reproduce the above copyright
        notice, this list of conditions and the following disclaimer in the
        documentation and/or other materials provided with the distribution.
  
      THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
      IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
      TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
      PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
      HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
      SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
      LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
      DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
      THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
      (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
      OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
  
   */
  
  
  // core/transform.cpp*
  #include "stdafx.h"
  #include "transform.h"
  #include "shape.h"
  
  // Matrix4x4 Method Definitions
  bool SolveLinearSystem2x2(const float A[2][2],
          const float B[2], float *x0, float *x1) {
      float det = A[0][0]*A[1][1] - A[0][1]*A[1][0];
      if (fabsf(det) < 1e-10f)
          return false;
      *x0 = (A[1][1]*B[0] - A[0][1]*B[1]) / det;
      *x1 = (A[0][0]*B[1] - A[1][0]*B[0]) / det;
      if (isnan(*x0) || isnan(*x1))
          return false;
      return true;
  }
  
  
  Matrix4x4::Matrix4x4(float mat[4][4]) {
      memcpy(m, mat, 16*sizeof(float));
  }
  
  
  Matrix4x4::Matrix4x4(float t00, float t01, float t02, float t03,
                       float t10, float t11, float t12, float t13,
                       float t20, float t21, float t22, float t23,
                       float t30, float t31, float t32, float t33) {
      m[0][0] = t00; m[0][1] = t01; m[0][2] = t02; m[0][3] = t03;
      m[1][0] = t10; m[1][1] = t11; m[1][2] = t12; m[1][3] = t13;
      m[2][0] = t20; m[2][1] = t21; m[2][2] = t22; m[2][3] = t23;
      m[3][0] = t30; m[3][1] = t31; m[3][2] = t32; m[3][3] = t33;
  }
  
  
  Matrix4x4 Transpose(const Matrix4x4 &m) {
     return Matrix4x4(m.m[0][0], m.m[1][0], m.m[2][0], m.m[3][0],
                      m.m[0][1], m.m[1][1], m.m[2][1], m.m[3][1],
                      m.m[0][2], m.m[1][2], m.m[2][2], m.m[3][2],
                      m.m[0][3], m.m[1][3], m.m[2][3], m.m[3][3]);
  }
  
  
  Matrix4x4 Inverse(const Matrix4x4 &m) {
      int indxc[4], indxr[4];
      int ipiv[4] = { 0, 0, 0, 0 };
      float minv[4][4];
      memcpy(minv, m.m, 4*4*sizeof(float));
      for (int i = 0; i < 4; i++) {
          int irow = -1, icol = -1;
          float big = 0.;
          // Choose pivot
          for (int j = 0; j < 4; j++) {
              if (ipiv[j] != 1) {
                  for (int k = 0; k < 4; k++) {
                      if (ipiv[k] == 0) {
                          if (fabsf(minv[j][k]) >= big) {
                              big = float(fabsf(minv[j][k]));
                              irow = j;
                              icol = k;
                          }
                      }
                      else if (ipiv[k] > 1)
                          Error("Singular matrix in MatrixInvert");
                  }
              }
          }
          ++ipiv[icol];
          // Swap rows _irow_ and _icol_ for pivot
          if (irow != icol) {
              for (int k = 0; k < 4; ++k)
                  swap(minv[irow][k], minv[icol][k]);
          }
          indxr[i] = irow;
          indxc[i] = icol;
          if (minv[icol][icol] == 0.)
              Error("Singular matrix in MatrixInvert");
  
          // Set $m[icol][icol]$ to one by scaling row _icol_ appropriately
          float pivinv = 1.f / minv[icol][icol];
          minv[icol][icol] = 1.f;
          for (int j = 0; j < 4; j++)
              minv[icol][j] *= pivinv;
  
          // Subtract this row from others to zero out their columns
          for (int j = 0; j < 4; j++) {
              if (j != icol) {
                  float save = minv[j][icol];
                  minv[j][icol] = 0;
                  for (int k = 0; k < 4; k++)
                      minv[j][k] -= minv[icol][k]*save;
              }
          }
      }
      // Swap columns to reflect permutation
      for (int j = 3; j >= 0; j--) {
          if (indxr[j] != indxc[j]) {
              for (int k = 0; k < 4; k++)
                  swap(minv[k][indxr[j]], minv[k][indxc[j]]);
          }
      }
      return Matrix4x4(minv);
  }
  
  
  
  // Transform Method Definitions
  void Transform::Print(FILE *f) const {
      m.Print(f);
  }
  
  
  Transform Translate(const Vector &delta) {
      Matrix4x4 m(1, 0, 0, delta.x,
                  0, 1, 0, delta.y,
                  0, 0, 1, delta.z,
                  0, 0, 0,       1);
      Matrix4x4 minv(1, 0, 0, -delta.x,
                     0, 1, 0, -delta.y,
                     0, 0, 1, -delta.z,
                     0, 0, 0,        1);
      return Transform(m, minv);
  }
  
  
  Transform Scale(float x, float y, float z) {
      Matrix4x4 m(x, 0, 0, 0,
                  0, y, 0, 0,
                  0, 0, z, 0,
                  0, 0, 0, 1);
      Matrix4x4 minv(1.f/x,     0,     0,     0,
                     0,     1.f/y,     0,     0,
                     0,         0,     1.f/z, 0,
                     0,         0,     0,     1);
      return Transform(m, minv);
  }
  
  
  Transform RotateX(float angle) {
      float sin_t = sinf(Radians(angle));
      float cos_t = cosf(Radians(angle));
      Matrix4x4 m(1,     0,      0, 0,
                  0, cos_t, -sin_t, 0,
                  0, sin_t,  cos_t, 0,
                  0,     0,      0, 1);
      return Transform(m, Transpose(m));
  }
  
  
  Transform RotateY(float angle) {
      float sin_t = sinf(Radians(angle));
      float cos_t = cosf(Radians(angle));
      Matrix4x4 m( cos_t,   0,  sin_t, 0,
                   0,   1,      0, 0,
                  -sin_t,   0,  cos_t, 0,
                   0,   0,   0,    1);
      return Transform(m, Transpose(m));
  }
  
  
  
  Transform RotateZ(float angle) {
      float sin_t = sinf(Radians(angle));
      float cos_t = cosf(Radians(angle));
      Matrix4x4 m(cos_t, -sin_t, 0, 0,
                  sin_t,  cos_t, 0, 0,
                  0,      0, 1, 0,
                  0,      0, 0, 1);
      return Transform(m, Transpose(m));
  }
  
  
  Transform Rotate(float angle, const Vector &axis) {
      Vector a = Normalize(axis);
      float s = sinf(Radians(angle));
      float c = cosf(Radians(angle));
      float m[4][4];
  
      m[0][0] = a.x * a.x + (1.f - a.x * a.x) * c;
      m[0][1] = a.x * a.y * (1.f - c) - a.z * s;
      m[0][2] = a.x * a.z * (1.f - c) + a.y * s;
      m[0][3] = 0;
  
      m[1][0] = a.x * a.y * (1.f - c) + a.z * s;
      m[1][1] = a.y * a.y + (1.f - a.y * a.y) * c;
      m[1][2] = a.y * a.z * (1.f - c) - a.x * s;
      m[1][3] = 0;
  
      m[2][0] = a.x * a.z * (1.f - c) - a.y * s;
      m[2][1] = a.y * a.z * (1.f - c) + a.x * s;
      m[2][2] = a.z * a.z + (1.f - a.z * a.z) * c;
      m[2][3] = 0;
  
      m[3][0] = 0;
      m[3][1] = 0;
      m[3][2] = 0;
      m[3][3] = 1;
  
      Matrix4x4 mat(m);
      return Transform(mat, Transpose(mat));
  }
  
  
  Transform LookAt(const Point &pos, const Point &look, const Vector &up) {
      float m[4][4];
      // Initialize fourth column of viewing matrix
      m[0][3] = pos.x;
      m[1][3] = pos.y;
      m[2][3] = pos.z;
      m[3][3] = 1;
  
      // Initialize first three columns of viewing matrix
      Vector dir = Normalize(look - pos);
      if (Cross(Normalize(up), dir).Length() == 0) {
          Error("\"up\" vector (%f, %f, %f) and viewing direction (%f, %f, %f) "
                "passed to LookAt are pointing in the same direction.  Using "
                "the identity transformation.", up.x, up.y, up.z, dir.x, dir.y,
                dir.z);
          return Transform();
      }
      Vector left = Normalize(Cross(Normalize(up), dir));
      Vector newUp = Cross(dir, left);
      m[0][0] = left.x;
      m[1][0] = left.y;
      m[2][0] = left.z;
      m[3][0] = 0.;
      m[0][1] = newUp.x;
      m[1][1] = newUp.y;
      m[2][1] = newUp.z;
      m[3][1] = 0.;
      m[0][2] = dir.x;
      m[1][2] = dir.y;
      m[2][2] = dir.z;
      m[3][2] = 0.;
      Matrix4x4 camToWorld(m);
      return Transform(Inverse(camToWorld), camToWorld);
  }
  
  
  BBox Transform::operator()(const BBox &b) const {
      const Transform &M = *this;
      BBox ret(        M(Point(b.pMin.x, b.pMin.y, b.pMin.z)));
      ret = Union(ret, M(Point(b.pMax.x, b.pMin.y, b.pMin.z)));
      ret = Union(ret, M(Point(b.pMin.x, b.pMax.y, b.pMin.z)));
      ret = Union(ret, M(Point(b.pMin.x, b.pMin.y, b.pMax.z)));
      ret = Union(ret, M(Point(b.pMin.x, b.pMax.y, b.pMax.z)));
      ret = Union(ret, M(Point(b.pMax.x, b.pMax.y, b.pMin.z)));
      ret = Union(ret, M(Point(b.pMax.x, b.pMin.y, b.pMax.z)));
      ret = Union(ret, M(Point(b.pMax.x, b.pMax.y, b.pMax.z)));
      return ret;
  }
  
  
  Transform Transform::operator*(const Transform &t2) const {
      Matrix4x4 m1 = Matrix4x4::Mul(m, t2.m);
      Matrix4x4 m2 = Matrix4x4::Mul(t2.mInv, mInv);
      return Transform(m1, m2);
  }
  
  
  bool Transform::SwapsHandedness() const {
      float det = ((m.m[0][0] *
                    (m.m[1][1] * m.m[2][2] -
                     m.m[1][2] * m.m[2][1])) -
                   (m.m[0][1] *
                    (m.m[1][0] * m.m[2][2] -
                     m.m[1][2] * m.m[2][0])) +
                   (m.m[0][2] *
                    (m.m[1][0] * m.m[2][1] -
                     m.m[1][1] * m.m[2][0])));
      return det < 0.f;
  }
  
  
  Transform Orthographic(float znear, float zfar) {
      return Scale(1.f, 1.f, 1.f / (zfar-znear)) *
             Translate(Vector(0.f, 0.f, -znear));
  }
  
  
  Transform Perspective(float fov, float n, float f) {
      // Perform projective divide
      Matrix4x4 persp = Matrix4x4(1, 0,           0,              0,
                                  0, 1,           0,              0,
                                  0, 0, f / (f - n), -f*n / (f - n),
                                  0, 0,           1,              0);
  
      // Scale to canonical viewing volume
      float invTanAng = 1.f / tanf(Radians(fov) / 2.f);
      return Scale(invTanAng, invTanAng, 1) * Transform(persp);
  }
  
  
  
  // AnimatedTransform Method Definitions
  void AnimatedTransform::Decompose(const Matrix4x4 &m, Vector *T,
                                    Quaternion *Rquat, Matrix4x4 *S) {
      // Extract translation _T_ from transformation matrix
      T->x = m.m[0][3];
      T->y = m.m[1][3];
      T->z = m.m[2][3];
  
      // Compute new transformation matrix _M_ without translation
      Matrix4x4 M = m;
      for (int i = 0; i < 3; ++i)
          M.m[i][3] = M.m[3][i] = 0.f;
      M.m[3][3] = 1.f;
  
      // Extract rotation _R_ from transformation matrix
      float norm;
      int count = 0;
      Matrix4x4 R = M;
      do {
          // Compute next matrix _Rnext_ in series
          Matrix4x4 Rnext;
          Matrix4x4 Rit = Inverse(Transpose(R));
          for (int i = 0; i < 4; ++i)
              for (int j = 0; j < 4; ++j)
                  Rnext.m[i][j] = 0.5f * (R.m[i][j] + Rit.m[i][j]);
  
          // Compute norm of difference between _R_ and _Rnext_
          norm = 0.f;
          for (int i = 0; i < 3; ++i) {
              float n = fabsf(R.m[i][0] - Rnext.m[i][0]) +
                        fabsf(R.m[i][1] - Rnext.m[i][1]) +
                        fabsf(R.m[i][2] - Rnext.m[i][2]);
              norm = max(norm, n);
          }
          R = Rnext;
      } while (++count < 100 && norm > .0001f);
      // XXX TODO FIXME deal with flip...
      *Rquat = Quaternion(R);
  
      // Compute scale _S_ using rotation and original matrix
      *S = Matrix4x4::Mul(Inverse(R), M);
  }
  
  
  void AnimatedTransform::Interpolate(float time, Transform *t) const {
      // Handle boundary conditions for matrix interpolation
      if (!actuallyAnimated || time <= startTime) {
          *t = *startTransform;
          return;
      }
      if (time >= endTime) {
          *t = *endTransform;
          return;
      }
      float dt = (time - startTime) / (endTime - startTime);
      // Interpolate translation at _dt_
      Vector trans = (1.f - dt) * T[0] + dt * T[1];
  
      // Interpolate rotation at _dt_
      Quaternion rotate = Slerp(dt, R[0], R[1]);
  
      // Interpolate scale at _dt_
      Matrix4x4 scale;
      for (int i = 0; i < 3; ++i)
          for (int j = 0; j < 3; ++j)
              scale.m[i][j] = Lerp(dt, S[0].m[i][j], S[1].m[i][j]);
  
      // Compute interpolated matrix as product of interpolated components
      *t = Translate(trans) * rotate.ToTransform() * Transform(scale);
  }
  
  
  BBox AnimatedTransform::MotionBounds(const BBox &b,
                                       bool useInverse) const {
      if (!actuallyAnimated) 
        return useInverse ? Inverse(*startTransform)(b) : (*startTransform)(b);
      BBox ret;
      const int nSteps = 128;
      for (int i = 0; i < nSteps; ++i) {
          Transform t;
          float time = Lerp(float(i)/float(nSteps-1), startTime, endTime);
          Interpolate(time, &t);
          if (useInverse) t = Inverse(t);
          ret = Union(ret, t(b));
      }
      return ret;
  }
  
  
  void AnimatedTransform::operator()(const Ray &r, Ray *tr) const {
      if (!actuallyAnimated || r.time <= startTime)
          (*startTransform)(r, tr);
      else if (r.time >= endTime)
          (*endTransform)(r, tr);
      else {
          Transform t;
          Interpolate(r.time, &t);
          t(r, tr);
      }
      tr->time = r.time;
  }
  
  
  void AnimatedTransform::operator()(const RayDifferential &r,
      RayDifferential *tr) const {
      if (!actuallyAnimated || r.time <= startTime)
          (*startTransform)(r, tr);
      else if (r.time >= endTime)
          (*endTransform)(r, tr);
      else {
          Transform t;
          Interpolate(r.time, &t);
          t(r, tr);
      }
      tr->time = r.time;
  }
  
  
  Point AnimatedTransform::operator()(float time, const Point &p) const {
      if (!actuallyAnimated || time <= startTime)
          return (*startTransform)(p);
      else if (time >= endTime)
          return (*endTransform)(p);
      Transform t;
      Interpolate(time, &t);
      return t(p);
  }
  
  
  Vector AnimatedTransform::operator()(float time, const Vector &v) const {
      if (!actuallyAnimated || time <= startTime)
          return (*startTransform)(v);
      else if (time >= endTime)
          return (*endTransform)(v);
      Transform t;
      Interpolate(time, &t);
      return t(v);
  }
  
  
  Ray AnimatedTransform::operator()(const Ray &r) const {
      Ray ret;
      (*this)(r, &ret);
      return ret;
  }
  
  
  