Make the bottom [0 0 0 1] of affine transformation matricies implicit.

3 files changed, 60 insertions, 76 deletions

diff --git a/libdimension/dimension/geometry.h b/libdimension/dimension/geometry.h
index 831f48c..f5a49ec 100644
--- a/libdimension/dimension/geometry.h
+++ b/libdimension/dimension/geometry.h
@@ -41,9 +41,9 @@ typedef struct dmnsn_vector {
 /** The appropriate arguements to printf() a vector. */
 #define DMNSN_VECTOR_PRINTF(v) (v).x, (v).y, (v).z
 
-/** A 4x4 affine transformation matrix. */
+/** A 4x4 affine transformation matrix, with implied [0 0 0 1] bottom row. */
 typedef struct dmnsn_matrix {
-  double n[4][4]; /**< The matrix elements in row-major order. */
+  double n[3][4]; /**< The matrix elements in row-major order. */
 } dmnsn_matrix;
 
 /** A standard format string for matricies. */
@@ -57,7 +57,7 @@ typedef struct dmnsn_matrix {
   (m).n[0][0], (m).n[0][1], (m).n[0][2], (m).n[0][3],   \
   (m).n[1][0], (m).n[1][1], (m).n[1][2], (m).n[1][3],   \
   (m).n[2][0], (m).n[2][1], (m).n[2][2], (m).n[2][3],   \
-  (m).n[3][0], (m).n[3][1], (m).n[3][2], (m).n[3][3]
+  0.0, 0.0, 0.0, 1.0
 
 /** A line, or ray. */
 typedef struct dmnsn_line {
@@ -145,17 +145,15 @@ dmnsn_new_vector(double x, double y, double z)
   return v;
 }
 
-/** Construct a new matrix. */
+/** Construct a new transformation matrix. */
 DMNSN_INLINE dmnsn_matrix
 dmnsn_new_matrix(double a0, double a1, double a2, double a3,
                  double b0, double b1, double b2, double b3,
-                 double c0, double c1, double c2, double c3,
-                 double d0, double d1, double d2, double d3)
+                 double c0, double c1, double c2, double c3)
 {
   dmnsn_matrix m = { { { a0, a1, a2, a3 },
                        { b0, b1, b2, b3 },
-                       { c0, c1, c2, c3 },
-                       { d0, d1, d2, d3 } } };
+                       { c0, c1, c2, c3 } } };
   return m;
 }
 
@@ -378,16 +376,12 @@ dmnsn_matrix dmnsn_matrix_mul(dmnsn_matrix lhs, dmnsn_matrix rhs);
 DMNSN_INLINE dmnsn_vector
 dmnsn_transform_vector(dmnsn_matrix T, dmnsn_vector v)
 {
-  /* 12 multiplications, 3 divisions, 12 additions */
+  /* 9 multiplications, 9 additions */
   dmnsn_vector r;
-  double w;
-
   r.x = T.n[0][0]*v.x + T.n[0][1]*v.y + T.n[0][2]*v.z + T.n[0][3];
   r.y = T.n[1][0]*v.x + T.n[1][1]*v.y + T.n[1][2]*v.z + T.n[1][3];
   r.z = T.n[2][0]*v.x + T.n[2][1]*v.y + T.n[2][2]*v.z + T.n[2][3];
-  w   = T.n[3][0]*v.x + T.n[3][1]*v.y + T.n[3][2]*v.z + T.n[3][3];
-
-  return dmnsn_vector_div(r, w);
+  return r;
 }
 
 /** Transform a bounding box by a matrix. */
@@ -402,7 +396,7 @@ dmnsn_bounding_box dmnsn_transform_bounding_box(dmnsn_matrix T,
 DMNSN_INLINE dmnsn_line
 dmnsn_transform_line(dmnsn_matrix T, dmnsn_line l)
 {
-  /* 24 multiplications, 6 divisions, 30 additions */
+  /* 18 multiplications, 24 additions */
   dmnsn_line ret;
   ret.x0 = dmnsn_transform_vector(T, l.x0);
   ret.n  = dmnsn_vector_sub(
diff --git a/libdimension/geometry.c b/libdimension/geometry.c
index 1f0b054..11a59d3 100644
--- a/libdimension/geometry.c
+++ b/libdimension/geometry.c
@@ -32,8 +32,7 @@ dmnsn_identity_matrix()
 {
   return dmnsn_new_matrix(1.0, 0.0, 0.0, 0.0,
                           0.0, 1.0, 0.0, 0.0,
-                          0.0, 0.0, 1.0, 0.0,
-                          0.0, 0.0, 0.0, 1.0);
+                          0.0, 0.0, 1.0, 0.0);
 }
 
 /* Scaling matrix */
@@ -42,8 +41,7 @@ dmnsn_scale_matrix(dmnsn_vector s)
 {
   return dmnsn_new_matrix(s.x, 0.0, 0.0, 0.0,
                           0.0, s.y, 0.0, 0.0,
-                          0.0, 0.0, s.z, 0.0,
-                          0.0, 0.0, 0.0, 1.0);
+                          0.0, 0.0, s.z, 0.0);
 }
 
 /* Translation matrix */
@@ -52,8 +50,7 @@ dmnsn_translation_matrix(dmnsn_vector d)
 {
   return dmnsn_new_matrix(1.0, 0.0, 0.0, d.x,
                           0.0, 1.0, 0.0, d.y,
-                          0.0, 0.0, 1.0, d.z,
-                          0.0, 0.0, 0.0, 1.0);
+                          0.0, 0.0, 1.0, d.z);
 }
 
 /* Left-handed rotation matrix; theta/|theta| = axis, |theta| = angle */
@@ -80,8 +77,7 @@ dmnsn_rotation_matrix(dmnsn_vector theta)
   return dmnsn_new_matrix(
     1.0 + t*(x*x - 1.0), -z*s + t*x*y,        y*s + t*x*z,         0.0,
     z*s + t*x*y,         1.0 + t*(y*y - 1.0), -x*s + t*y*z,        0.0,
-    -y*s + t*x*z,        x*s + t*y*z,         1.0 + t*(z*z - 1.0), 0.0,
-    0.0,                 0.0,                 0.0,                 1.0
+    -y*s + t*x*z,        x*s + t*y*z,         1.0 + t*(z*z - 1.0), 0.0
   );
 }
 
@@ -138,11 +134,11 @@ dmnsn_matrix_inverse(dmnsn_matrix A)
   /*
    * Use partitioning to invert a matrix:
    *
-   *     ( P Q ) -1
-   *     ( R S )
+   *     [ P Q ] -1
+   *     [ R S ]
    *
-   *   = ( PP QQ )
-   *     ( RR SS ),
+   *   = [ PP QQ ]
+   *     [ RR SS ],
    *
    * with PP = inv(P) - inv(P)*Q*RR,
    *      QQ = -inv(P)*Q*SS,
@@ -173,9 +169,9 @@ dmnsn_matrix_inverse(dmnsn_matrix A)
   Q = dmnsn_new_matrix2(A.n[0][2], A.n[0][3],
                         A.n[1][2], A.n[1][3]);
   R = dmnsn_new_matrix2(A.n[2][0], A.n[2][1],
-                        A.n[3][0], A.n[3][1]);
+                        0.0,       0.0);
   S = dmnsn_new_matrix2(A.n[2][2], A.n[2][3],
-                        A.n[3][2], A.n[3][3]);
+                        0.0,       1.0);
 
   /* Do this inversion ourselves, since we already have the determinant */
   Pi = dmnsn_new_matrix2( P.n[1][1]/Pdet, -P.n[0][1]/Pdet,
@@ -195,8 +191,7 @@ dmnsn_matrix_inverse(dmnsn_matrix A)
   /* Reconstruct the matrix */
   return dmnsn_new_matrix(PP.n[0][0], PP.n[0][1], QQ.n[0][0], QQ.n[0][1],
                           PP.n[1][0], PP.n[1][1], QQ.n[1][0], QQ.n[1][1],
-                          RR.n[0][0], RR.n[0][1], SS.n[0][0], SS.n[0][1],
-                          RR.n[1][0], RR.n[1][1], SS.n[1][0], SS.n[1][1]);
+                          RR.n[0][0], RR.n[0][1], SS.n[0][0], SS.n[0][1]);
 }
 
 /* For nice shorthand */
@@ -255,48 +250,56 @@ static dmnsn_matrix
 dmnsn_matrix_inverse_generic(dmnsn_matrix A)
 {
   /*
-   * Simply form the matrix's adjugate and divide each element by the
-   * determinant as we go.  The routine itself has 4 additions and 16 divisions
-   * plus 16 cofactor calculations, giving 144 multiplications, 84 additions,
-   * and 16 divisions.
+   * For A = [ A'      b ], A^-1 = [ A'^-1   -(A'^-1)*b ].
+   *         [ 0 ... 0 1 ]         [ 0 ... 0      1     ]
+   *
+   * Invert A' by calculating its adjucate.
    */
   dmnsn_matrix inv;
   double det = 0.0, C;
 
   /* Perform a Laplace expansion along the first row to give us the adjugate's
      first column and the determinant */
-  for (size_t j = 0; j < 4; ++j) {
+  for (size_t j = 0; j < 3; ++j) {
     C = dmnsn_matrix_cofactor(A, 0, j);
     det += A.n[0][j]*C;
     inv.n[j][0] = C;
   }
 
   /* Divide the first column by the determinant */
-  for (size_t j = 0; j < 4; ++j) {
+  for (size_t j = 0; j < 3; ++j) {
     inv.n[j][0] /= det;
   }
 
-  /* Find columns 2 through 4 */
-  for (size_t i = 1; i < 4; ++i) {
-    for (size_t j = 0; j < 4; ++j) {
+  /* Find the rest of A' */
+  for (size_t j = 0; j < 3; ++j) {
+    for (size_t i = 1; i < 3; ++i) {
       inv.n[j][i] = dmnsn_matrix_cofactor(A, i, j)/det;
     }
+    inv.n[j][3] = 0.0;
+  }
+
+  /* Find the translational component of the inverse */
+  for (size_t i = 0; i < 3; ++i) {
+    for (size_t j = 0; j < 3; ++j) {
+      inv.n[i][3] -= inv.n[i][j]*A.n[i][3];
+    }
   }
 
   return inv;
 }
 
-/* Gives the cofactor at row, col; the determinant of the matrix formed from A
-   by ignoring row `row' and column `col', times (-1)**(row + col) */
+/* Gives the cofactor at row, col; the determinant of the matrix formed from the
+   upper-left 3x3 corner of A by ignoring row `row' and column `col',
+   times (-1)^(row + col) */
 static double
 dmnsn_matrix_cofactor(dmnsn_matrix A, unsigned int row, unsigned int col)
 {
-  /* 9 multiplications, 5 additions */
-  double n[9], C;
-  unsigned int k = 0;
-
-  for (size_t i = 0; i < 4; ++i) {
-    for (size_t j = 0; j < 4; ++j) {
+  /* 2 multiplications, 1 addition */
+  double n[4];
+  size_t k = 0;
+  for (size_t i = 0; i < 3; ++i) {
+    for (size_t j = 0; j < 3; ++j) {
       if (i != row && j != col) {
         n[k] = A.n[i][j];
         ++k;
@@ -304,8 +307,7 @@ dmnsn_matrix_cofactor(dmnsn_matrix A, unsigned int row, unsigned int col)
     }
   }
 
-  C = n[0]*(n[4]*n[8] - n[5]*n[7]) + n[1]*(n[5]*n[6] - n[3]*n[8])
-    + n[2]*(n[3]*n[7] - n[4]*n[6]);
+  double C = n[0]*n[3] - n[1]*n[2];
   if ((row + col)%2 == 0) {
     return C;
   } else {
@@ -317,44 +319,35 @@ dmnsn_matrix_cofactor(dmnsn_matrix A, unsigned int row, unsigned int col)
 dmnsn_matrix
 dmnsn_matrix_mul(dmnsn_matrix lhs, dmnsn_matrix rhs)
 {
-  /* 64 multiplications, 48 additions */
+  /* 36 multiplications, 27 additions */
   dmnsn_matrix r;
 
   r.n[0][0] = lhs.n[0][0]*rhs.n[0][0] + lhs.n[0][1]*rhs.n[1][0]
-    + lhs.n[0][2]*rhs.n[2][0] + lhs.n[0][3]*rhs.n[3][0];
+    + lhs.n[0][2]*rhs.n[2][0];
   r.n[0][1] = lhs.n[0][0]*rhs.n[0][1] + lhs.n[0][1]*rhs.n[1][1]
-    + lhs.n[0][2]*rhs.n[2][1] + lhs.n[0][3]*rhs.n[3][1];
+    + lhs.n[0][2]*rhs.n[2][1];
   r.n[0][2] = lhs.n[0][0]*rhs.n[0][2] + lhs.n[0][1]*rhs.n[1][2]
-    + lhs.n[0][2]*rhs.n[2][2] + lhs.n[0][3]*rhs.n[3][2];
+    + lhs.n[0][2]*rhs.n[2][2];
   r.n[0][3] = lhs.n[0][0]*rhs.n[0][3] + lhs.n[0][1]*rhs.n[1][3]
-    + lhs.n[0][2]*rhs.n[2][3] + lhs.n[0][3]*rhs.n[3][3];
+    + lhs.n[0][2]*rhs.n[2][3] + lhs.n[0][3];
 
   r.n[1][0] = lhs.n[1][0]*rhs.n[0][0] + lhs.n[1][1]*rhs.n[1][0]
-    + lhs.n[1][2]*rhs.n[2][0] + lhs.n[1][3]*rhs.n[3][0];
+    + lhs.n[1][2]*rhs.n[2][0];
   r.n[1][1] = lhs.n[1][0]*rhs.n[0][1] + lhs.n[1][1]*rhs.n[1][1]
-    + lhs.n[1][2]*rhs.n[2][1] + lhs.n[1][3]*rhs.n[3][1];
+    + lhs.n[1][2]*rhs.n[2][1];
   r.n[1][2] = lhs.n[1][0]*rhs.n[0][2] + lhs.n[1][1]*rhs.n[1][2]
-    + lhs.n[1][2]*rhs.n[2][2] + lhs.n[1][3]*rhs.n[3][2];
+    + lhs.n[1][2]*rhs.n[2][2];
   r.n[1][3] = lhs.n[1][0]*rhs.n[0][3] + lhs.n[1][1]*rhs.n[1][3]
-    + lhs.n[1][2]*rhs.n[2][3] + lhs.n[1][3]*rhs.n[3][3];
+    + lhs.n[1][2]*rhs.n[2][3] + lhs.n[1][3];
 
   r.n[2][0] = lhs.n[2][0]*rhs.n[0][0] + lhs.n[2][1]*rhs.n[1][0]
-    + lhs.n[2][2]*rhs.n[2][0] + lhs.n[2][3]*rhs.n[3][0];
+    + lhs.n[2][2]*rhs.n[2][0];
   r.n[2][1] = lhs.n[2][0]*rhs.n[0][1] + lhs.n[2][1]*rhs.n[1][1]
-    + lhs.n[2][2]*rhs.n[2][1] + lhs.n[2][3]*rhs.n[3][1];
+    + lhs.n[2][2]*rhs.n[2][1];
   r.n[2][2] = lhs.n[2][0]*rhs.n[0][2] + lhs.n[2][1]*rhs.n[1][2]
-    + lhs.n[2][2]*rhs.n[2][2] + lhs.n[2][3]*rhs.n[3][2];
+    + lhs.n[2][2]*rhs.n[2][2];
   r.n[2][3] = lhs.n[2][0]*rhs.n[0][3] + lhs.n[2][1]*rhs.n[1][3]
-    + lhs.n[2][2]*rhs.n[2][3] + lhs.n[2][3]*rhs.n[3][3];
-
-  r.n[3][0] = lhs.n[3][0]*rhs.n[0][0] + lhs.n[3][1]*rhs.n[1][0]
-    + lhs.n[3][2]*rhs.n[2][0] + lhs.n[3][3]*rhs.n[3][0];
-  r.n[3][1] = lhs.n[3][0]*rhs.n[0][1] + lhs.n[3][1]*rhs.n[1][1]
-    + lhs.n[3][2]*rhs.n[2][1] + lhs.n[3][3]*rhs.n[3][1];
-  r.n[3][2] = lhs.n[3][0]*rhs.n[0][2] + lhs.n[3][1]*rhs.n[1][2]
-    + lhs.n[3][2]*rhs.n[2][2] + lhs.n[3][3]*rhs.n[3][2];
-  r.n[3][3] = lhs.n[3][0]*rhs.n[0][3] + lhs.n[3][1]*rhs.n[1][3]
-    + lhs.n[3][2]*rhs.n[2][3] + lhs.n[3][3]*rhs.n[3][3];
+    + lhs.n[2][2]*rhs.n[2][3] + lhs.n[2][3];
 
   return r;
 }
diff --git a/libdimension/object.c b/libdimension/object.c
index bb4eeb1..5914b46 100644
--- a/libdimension/object.c
+++ b/libdimension/object.c
@@ -94,10 +94,7 @@ dmnsn_transform_normal(dmnsn_matrix trans, dmnsn_vector normal)
     dmnsn_vector_sub(
       dmnsn_transform_vector(trans, normal),
       /* Optimized form of dmnsn_transform_vector(trans, dmnsn_zero) */
-      dmnsn_vector_div(
-        dmnsn_new_vector(trans.n[0][3], trans.n[1][3], trans.n[2][3]),
-        trans.n[3][3]
-      )
+      dmnsn_new_vector(trans.n[0][3], trans.n[1][3], trans.n[2][3])
     )
   );
 }