X-Git-Url: https://git.xonotic.org/?a=blobdiff_plain;f=matrixlib.c;h=cfc4627a5b1e544cf01e32417af00179866130fd;hb=c9a3afc3de13d6956165e41cb93ceb6a0e8c2763;hp=bca33332e32940e29a1c04e290aafe4be9a5c0d2;hpb=fb7e8ff7102e5e873a8b9eac4cfbf4bab557a1d4;p=xonotic%2Fdarkplaces.git diff --git a/matrixlib.c b/matrixlib.c index bca33332..cfc4627a 100644 --- a/matrixlib.c +++ b/matrixlib.c @@ -143,6 +143,60 @@ void Matrix4x4_Transpose (matrix4x4_t *out, const matrix4x4_t *in1) #if 1 // Adapted from code contributed to Mesa by David Moore (Mesa 7.6 under SGI Free License B - which is MIT/X11-type) +// added helper for common subexpression elimination by eihrul, and other optimizations by div0 +int Matrix4x4_Invert_Full (matrix4x4_t *out, const matrix4x4_t *in1) +{ + float det; + + // note: orientation does not matter, as transpose(invert(transpose(m))) == invert(m), proof: + // transpose(invert(transpose(m))) * m + // = transpose(invert(transpose(m))) * transpose(transpose(m)) + // = transpose(transpose(m) * invert(transpose(m))) + // = transpose(identity) + // = identity + + // this seems to help gcc's common subexpression elimination, and also makes the code look nicer + float m00 = in1->m[0][0], m01 = in1->m[0][1], m02 = in1->m[0][2], m03 = in1->m[0][3], + m10 = in1->m[1][0], m11 = in1->m[1][1], m12 = in1->m[1][2], m13 = in1->m[1][3], + m20 = in1->m[2][0], m21 = in1->m[2][1], m22 = in1->m[2][2], m23 = in1->m[2][3], + m30 = in1->m[3][0], m31 = in1->m[3][1], m32 = in1->m[3][2], m33 = in1->m[3][3]; + + // calculate the adjoint + out->m[0][0] = (m11*(m22*m33 - m23*m32) - m21*(m12*m33 - m13*m32) + m31*(m12*m23 - m13*m22)); + out->m[0][1] = -(m01*(m22*m33 - m23*m32) - m21*(m02*m33 - m03*m32) + m31*(m02*m23 - m03*m22)); + out->m[0][2] = (m01*(m12*m33 - m13*m32) - m11*(m02*m33 - m03*m32) + m31*(m02*m13 - m03*m12)); + out->m[0][3] = -(m01*(m12*m23 - m13*m22) - m11*(m02*m23 - m03*m22) + m21*(m02*m13 - m03*m12)); + out->m[1][0] = -(m10*(m22*m33 - m23*m32) - m20*(m12*m33 - m13*m32) + m30*(m12*m23 - m13*m22)); + out->m[1][1] = (m00*(m22*m33 - m23*m32) - m20*(m02*m33 - m03*m32) + m30*(m02*m23 - m03*m22)); + out->m[1][2] = -(m00*(m12*m33 - m13*m32) - m10*(m02*m33 - m03*m32) + m30*(m02*m13 - m03*m12)); + out->m[1][3] = (m00*(m12*m23 - m13*m22) - m10*(m02*m23 - m03*m22) + m20*(m02*m13 - m03*m12)); + out->m[2][0] = (m10*(m21*m33 - m23*m31) - m20*(m11*m33 - m13*m31) + m30*(m11*m23 - m13*m21)); + out->m[2][1] = -(m00*(m21*m33 - m23*m31) - m20*(m01*m33 - m03*m31) + m30*(m01*m23 - m03*m21)); + out->m[2][2] = (m00*(m11*m33 - m13*m31) - m10*(m01*m33 - m03*m31) + m30*(m01*m13 - m03*m11)); + out->m[2][3] = -(m00*(m11*m23 - m13*m21) - m10*(m01*m23 - m03*m21) + m20*(m01*m13 - m03*m11)); + out->m[3][0] = -(m10*(m21*m32 - m22*m31) - m20*(m11*m32 - m12*m31) + m30*(m11*m22 - m12*m21)); + out->m[3][1] = (m00*(m21*m32 - m22*m31) - m20*(m01*m32 - m02*m31) + m30*(m01*m22 - m02*m21)); + out->m[3][2] = -(m00*(m11*m32 - m12*m31) - m10*(m01*m32 - m02*m31) + m30*(m01*m12 - m02*m11)); + out->m[3][3] = (m00*(m11*m22 - m12*m21) - m10*(m01*m22 - m02*m21) + m20*(m01*m12 - m02*m11)); + + // calculate the determinant (as inverse == 1/det * adjoint, adjoint * m == identity * det, so this calculates the det) + det = m00*out->m[0][0] + m10*out->m[0][1] + m20*out->m[0][2] + m30*out->m[0][3]; + if (det == 0.0f) + return 0; + + // multiplications are faster than divisions, usually + det = 1.0f / det; + + // manually unrolled loop to multiply all matrix elements by 1/det + out->m[0][0] *= det; out->m[0][1] *= det; out->m[0][2] *= det; out->m[0][3] *= det; + out->m[1][0] *= det; out->m[1][1] *= det; out->m[1][2] *= det; out->m[1][3] *= det; + out->m[2][0] *= det; out->m[2][1] *= det; out->m[2][2] *= det; out->m[2][3] *= det; + out->m[3][0] *= det; out->m[3][1] *= det; out->m[3][2] *= det; out->m[3][3] *= det; + + return 1; +} +#elif 1 +// Adapted from code contributed to Mesa by David Moore (Mesa 7.6 under SGI Free License B - which is MIT/X11-type) int Matrix4x4_Invert_Full (matrix4x4_t *out, const matrix4x4_t *in1) { matrix4x4_t temp; @@ -1457,31 +1511,35 @@ void Matrix4x4_Transform3x3 (const matrix4x4_t *in, const float v[3], float out[ void Matrix4x4_TransformPositivePlane(const matrix4x4_t *in, float x, float y, float z, float d, float *o) { + float scale = sqrt(in->m[0][0] * in->m[0][0] + in->m[0][1] * in->m[0][1] + in->m[0][2] * in->m[0][2]); + float iscale = 1.0f / scale; #ifdef MATRIX4x4_OPENGLORIENTATION - o[0] = x * in->m[0][0] + y * in->m[1][0] + z * in->m[2][0]; - o[1] = x * in->m[0][1] + y * in->m[1][1] + z * in->m[2][1]; - o[2] = x * in->m[0][2] + y * in->m[1][2] + z * in->m[2][2]; - o[3] = d + (x * in->m[3][0] + y * in->m[3][1] + z * in->m[3][2]); + o[0] = (x * in->m[0][0] + y * in->m[1][0] + z * in->m[2][0]) * iscale; + o[1] = (x * in->m[0][1] + y * in->m[1][1] + z * in->m[2][1]) * iscale; + o[2] = (x * in->m[0][2] + y * in->m[1][2] + z * in->m[2][2]) * iscale; + o[3] = d * scale + (o[0] * in->m[3][0] + o[1] * in->m[3][1] + o[2] * in->m[3][2]); #else - o[0] = x * in->m[0][0] + y * in->m[0][1] + z * in->m[0][2]; - o[1] = x * in->m[1][0] + y * in->m[1][1] + z * in->m[1][2]; - o[2] = x * in->m[2][0] + y * in->m[2][1] + z * in->m[2][2]; - o[3] = d + (x * in->m[0][3] + y * in->m[1][3] + z * in->m[2][3]); + o[0] = (x * in->m[0][0] + y * in->m[0][1] + z * in->m[0][2]) * iscale; + o[1] = (x * in->m[1][0] + y * in->m[1][1] + z * in->m[1][2]) * iscale; + o[2] = (x * in->m[2][0] + y * in->m[2][1] + z * in->m[2][2]) * iscale; + o[3] = d * scale + (o[0] * in->m[0][3] + o[1] * in->m[1][3] + o[2] * in->m[2][3]); #endif } void Matrix4x4_TransformStandardPlane(const matrix4x4_t *in, float x, float y, float z, float d, float *o) { + float scale = sqrt(in->m[0][0] * in->m[0][0] + in->m[0][1] * in->m[0][1] + in->m[0][2] * in->m[0][2]); + float iscale = 1.0f / scale; #ifdef MATRIX4x4_OPENGLORIENTATION - o[0] = x * in->m[0][0] + y * in->m[1][0] + z * in->m[2][0]; - o[1] = x * in->m[0][1] + y * in->m[1][1] + z * in->m[2][1]; - o[2] = x * in->m[0][2] + y * in->m[1][2] + z * in->m[2][2]; - o[3] = d - (x * in->m[3][0] + y * in->m[3][1] + z * in->m[3][2]); + o[0] = (x * in->m[0][0] + y * in->m[1][0] + z * in->m[2][0]) * iscale; + o[1] = (x * in->m[0][1] + y * in->m[1][1] + z * in->m[2][1]) * iscale; + o[2] = (x * in->m[0][2] + y * in->m[1][2] + z * in->m[2][2]) * iscale; + o[3] = d * scale - (o[0] * in->m[3][0] + o[1] * in->m[3][1] + o[2] * in->m[3][2]); #else - o[0] = x * in->m[0][0] + y * in->m[0][1] + z * in->m[0][2]; - o[1] = x * in->m[1][0] + y * in->m[1][1] + z * in->m[1][2]; - o[2] = x * in->m[2][0] + y * in->m[2][1] + z * in->m[2][2]; - o[3] = d - (x * in->m[0][3] + y * in->m[1][3] + z * in->m[2][3]); + o[0] = (x * in->m[0][0] + y * in->m[0][1] + z * in->m[0][2]) * iscale; + o[1] = (x * in->m[1][0] + y * in->m[1][1] + z * in->m[1][2]) * iscale; + o[2] = (x * in->m[2][0] + y * in->m[2][1] + z * in->m[2][2]) * iscale; + o[3] = d * scale - (o[0] * in->m[0][3] + o[1] * in->m[1][3] + o[2] * in->m[2][3]); #endif } @@ -1609,3 +1667,17 @@ void Matrix4x4_Scale (matrix4x4_t *out, double rotatescale, double originscale) out->m[2][3] *= originscale; #endif } + +void Matrix4x4_Abs (matrix4x4_t *out) +{ + out->m[0][0] = fabs(out->m[0][0]); + out->m[0][1] = fabs(out->m[0][1]); + out->m[0][2] = fabs(out->m[0][2]); + out->m[1][0] = fabs(out->m[1][0]); + out->m[1][1] = fabs(out->m[1][1]); + out->m[1][2] = fabs(out->m[1][2]); + out->m[2][0] = fabs(out->m[2][0]); + out->m[2][1] = fabs(out->m[2][1]); + out->m[2][2] = fabs(out->m[2][2]); +} +