summaryrefslogtreecommitdiffstats
path: root/libdimension
diff options
context:
space:
mode:
Diffstat (limited to 'libdimension')
-rw-r--r--libdimension/dimension/geometry.h24
-rw-r--r--libdimension/geometry.c107
-rw-r--r--libdimension/object.c5
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])
)
);
}