add phong interpolation (lighting)

[chaz/rasterize] / mat.h

diff --git a/mat.h b/mat.h
index 499d09261697a8c02261c7ac94844475ef507f07..c188723346fce7faeec94210a2eca2d9a61fa042 100644 (file)
--- a/mat.h
+++ b/mat.h
@@ -237,22 +237,16 @@ mat_t MAT_ROTATE_Z(scal_t theta)
 INLINE_MAYBE
 mat_t MAT_ROTATE(scal_t theta, scal_t x, scal_t y, scal_t z)
 {
-    /*
-     * This code is an implementation of an algorithm described by Glenn
-     * Murray at http://inside.mines.edu/~gmurray/ArbitraryAxisRotation/
-     */
     vec_t v = vec_normalize(vec_new(x, y, z));
     x = v.x;
     y = v.y;
     z = v.z;
-    scal_t sin_a = scal_sin(theta);
-    scal_t cos_a = scal_cos(theta);
-    scal_t x2 = x * x;
-    scal_t y2 = y * y;
-    scal_t z2 = z * z;
-    return mat_new(x2+(S(1.0)-x2)*cos_a, x*y*(S(1.0)-cos_a)-z*sin_a, x*z*(S(1.0)-cos_a)+y*sin_a, S(0.0),
-                   x*y*(S(1.0)-cos_a)+z*sin_a, y2+(S(1.0)-y2)*cos_a, y*z*(S(1.0)-cos_a)-y*sin_a, S(0.0),
-                   x*z*(S(1.0)-cos_a)-y*sin_a, y*z*(S(1.0)-cos_a)+x*sin_a, z2+(S(1.0)-z2)*cos_a, S(0.0),
+    scal_t c = scal_cos(-theta);
+    scal_t s = scal_sin(-theta);
+    scal_t t = S(1.0) - c;
+    return mat_new(t * x * x + c, t * x * y + s * z, t * x * z - s * y, S(0.0),
+                   t * x * y - s * z, t * y * y + c, t * y * z + s * x, S(0.0),
+                   t * x * z + s * y, t * y * z - s * x, t * z * z + c, S(0.0),
                    S(0.0), S(0.0), S(0.0), S(1.0));
 }