diff options
Diffstat (limited to 'libdimension')
| -rw-r--r-- | libdimension/dimension/geometry.h | 24 | ||||
| -rw-r--r-- | libdimension/geometry.c | 107 | ||||
| -rw-r--r-- | libdimension/object.c | 5 |
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]) ) ); } |