diff engine/core/util/math/matrix.h @ 0:4a0efb7baf70

* Datasets becomes the new trunk and retires after that :-)
author mvbarracuda@33b003aa-7bff-0310-803a-e67f0ece8222
date Sun, 29 Jun 2008 18:44:17 +0000
parents
children 90005975cdbb
line wrap: on
line diff
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/engine/core/util/math/matrix.h	Sun Jun 29 18:44:17 2008 +0000
@@ -0,0 +1,401 @@
+/***************************************************************************
+ *   Copyright (C) 2005-2008 by the FIFE team                              *
+ *   http://www.fifengine.de                                               *
+ *   This file is part of FIFE.                                            *
+ *                                                                         *
+ *   FIFE is free software; you can redistribute it and/or modify          *
+ *   it under the terms of the GNU General Public License as published by  *
+ *   the Free Software Foundation; either version 2 of the License, or     *
+ *   (at your option) any later version.                                   *
+ *                                                                         *
+ *   This program is distributed in the hope that it will be useful,       *
+ *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
+ *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
+ *   GNU General Public License for more details.                          *
+ *                                                                         *
+ *   You should have received a copy of the GNU General Public License     *
+ *   along with this program; if not, write to the                         *
+ *   Free Software Foundation, Inc.,                                       *
+ *   51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA              *
+ ***************************************************************************/
+/***************************************************************************
+ * Includes some heavy copying from mathgl-pp project                      *
+ * (http://sourceforge.net/projects/mathgl-pp/)                            *
+ ***************************************************************************/
+
+#ifndef FIFE_UTIL_MATRIX_H
+#define FIFE_UTIL_MATRIX_H
+
+// Standard C++ library includes
+#include <cassert>
+#include <iostream>
+
+// Platform specific includes
+
+// 3rd party library includes
+
+// FIFE includes
+// These includes are split up in two parts, separated by one empty line
+// First block: files included from the FIFE root src directory
+// Second block: files included from the same folder
+#include "util/base/fife_stdint.h"
+#include "util/structures/point.h"
+
+#include "fife_math.h"
+
+namespace FIFE {
+
+
+	/** Minimal matrix class to assist in view 3d calculations
+	 */
+	template <typename T>
+	class Matrix {
+		public:
+			Matrix<T>() {}
+			template <typename U> friend class Matrix;
+			template <typename U> Matrix<T>(const Matrix<U>& mat) {
+				memmove(m, mat.m, 16*sizeof(T));
+			}
+			~Matrix() {}
+
+			/** Adjoint method inverse, constant time inversion implementation
+			 */
+			Matrix inverse() const {
+				Matrix ret(adjoint());
+
+				T determinant = m0*ret[0] + m1*ret[4] + m2*ret[8] + m3*ret[12];
+				assert(determinant!=0 && "Singular matrix has no inverse");
+
+				ret/=determinant;
+				return ret;
+			}
+
+			/** Divide this matrix by a scalar
+			 */
+			inline Matrix& operator/= (T val) {
+				for (register unsigned i = 0; i < 16; ++i)
+					m[i] /= val;
+				return *this;
+			}
+
+			/** Get the adjoint matrix
+			 */
+			Matrix adjoint() const {
+				Matrix ret;
+
+				ret[0] = cofactorm0();
+				ret[1] = -cofactorm4();
+				ret[2] = cofactorm8();
+				ret[3] = -cofactorm12();
+
+				ret[4] = -cofactorm1();
+				ret[5] = cofactorm5();
+				ret[6] = -cofactorm9();
+				ret[7] = cofactorm13();
+
+				ret[8] = cofactorm2();
+				ret[9] = -cofactorm6();
+				ret[10] = cofactorm10();
+				ret[11] = -cofactorm14();
+
+				ret[12] = -cofactorm3();
+				ret[13] = cofactorm7();
+				ret[14] = -cofactorm11();
+				ret[15] = cofactorm15();
+
+				return ret;
+			}
+
+
+			/** Make this a rotation matrix
+			 */
+
+			inline Matrix& loadRotate(T angle, T x, T y, T z) {
+				register T magSqr = x*x + y*y + z*z;
+				if (magSqr != 1.0) {
+					register T mag = sqrt(magSqr);
+					x/=mag;
+					y/=mag;
+					z/=mag;
+				}
+				T c = cos(angle*M_PI/180);
+				T s = sin(angle*M_PI/180);
+				m0 = x*x*(1-c)+c;
+				m1 = y*x*(1-c)+z*s;
+				m2 = z*x*(1-c)-y*s;
+				m3 = 0;
+
+				m4 = x*y*(1-c)-z*s;
+				m5 = y*y*(1-c)+c;
+				m6 = z*y*(1-c)+x*s;
+				m7 = 0;
+
+				m8 = x*z*(1-c)+y*s;
+				m9 = y*z*(1-c)-x*s;
+				m10 = z*z*(1-c)+c;
+				m11 = 0;
+
+				m12 = 0;
+				m13 = 0;
+				m14 = 0;
+				m15 = 1;
+
+				return *this;
+			}
+
+			/** Apply scale into this matrix
+			 */
+			inline Matrix& applyScale(T x, T y, T z) {
+				static Matrix<T> temp;
+				temp.loadScale(x,y,z);
+				*this = temp.mult4by4(*this);
+				return  *this;
+			}
+
+			/** Make this a scale matrix
+			 */
+			inline Matrix& loadScale(T x, T y, T z = 1) {
+				m0 = x;
+				m4 = 0;
+				m8  = 0;
+				m12 = 0;
+				m1 = 0;
+				m5 = y;
+				m9  = 0;
+				m13 = 0;
+				m2 = 0;
+				m6 = 0;
+				m10 = z;
+				m14 = 0;
+				m3 = 0;
+				m7 = 0;
+				m11 = 0;
+				m15 = 1;
+
+				return *this;
+			}
+
+			/** Apply translation into this matrix
+			 */
+			inline Matrix& applyTranslate(T x, T y, T z) {
+				static Matrix<T> temp;
+				temp.loadTranslate(x,y,z);
+				*this = temp.mult4by4(*this);
+				return  *this;
+			}
+
+			/** Make this a translation matrix
+			 */
+			inline Matrix& loadTranslate( const T x, const T y, const T z) {
+				m0 = 1;
+				m4 = 0;
+				m8  = 0;
+				m12 = x;
+				m1 = 0;
+				m5 = 1;
+				m9  = 0;
+				m13 = y;
+				m2 = 0;
+				m6 = 0;
+				m10 = 1;
+				m14 = z;
+				m3 = 0;
+				m7 = 0;
+				m11 = 0;
+				m15 = 1;
+
+				return *this;
+			}
+
+			/** Transform given point using this matrix
+			 */
+			inline PointType3D<T> operator* (const PointType3D<T>& vec) {
+				PointType3D<T> ret(T(0));
+				for (register unsigned j = 0; j < 3; ++j)
+					for (register unsigned i = 0; i < 3; ++i)
+						ret.val[j] += vec.val[i]*m[j+i*4]; //scale and rotate disregarding w scaling
+
+				for (register unsigned i = 0; i < 3; ++i)
+					ret.val[i] += m[i+3*4]; //translate
+
+				//do w scaling
+				T w = m[15];
+				for (register unsigned i = 0; i < 3; ++i)
+					w += vec.val[i]*m[3+i*4];
+
+				register T resip = 1/w;
+
+				for (register unsigned i = 0; i < 3; ++i)
+					ret[i] *= resip;
+
+				return ret;
+			}
+
+			/** Direct access to the matrix elements, just remember they are in column major format!!
+			 */
+			inline T& operator[] (int ind) {
+				assert(ind > -1 && ind < 16);
+				return m[ind];
+			}
+
+			/** Apply the matrix dot product to this matrix
+			 */
+			inline Matrix& mult3by3(const Matrix& mat) {
+				Matrix temp(*this);
+				m0 = temp.m0*mat.m0+temp.m4*mat.m1+temp.m8*mat.m2;
+				m4 = temp.m0*mat.m4+temp.m4*mat.m5+temp.m8*mat.m6;
+				m8 = temp.m0*mat.m8+temp.m4*mat.m9+temp.m8*mat.m10;
+
+				m1 = temp.m1*mat.m0+temp.m5*mat.m1+temp.m9*mat.m2;
+				m5 = temp.m1*mat.m4+temp.m5*mat.m5+temp.m9*mat.m6;
+				m9 = temp.m1*mat.m8+temp.m5*mat.m9+temp.m9*mat.m10;
+
+				m2 = temp.m2*mat.m0+temp.m6*mat.m1+temp.m10*mat.m2;
+				m6 = temp.m2*mat.m4+temp.m6*mat.m5+temp.m10*mat.m6;
+				m10 = temp.m2*mat.m8+temp.m6*mat.m9+temp.m10*mat.m10;
+
+				m3 = temp.m3*mat.m0+temp.m7*mat.m1+temp.m11*mat.m2;
+				m7 = temp.m3*mat.m4+temp.m7*mat.m5+temp.m11*mat.m6;
+				m11 = temp.m3*mat.m8+temp.m7*mat.m9+temp.m11*mat.m10;
+				return *this;
+			}
+
+			/** this->Rmult4by4(temp) == [temp] X [*this] **/
+			/** also equal to temp->mult4by4(*this) **/
+			inline Matrix<T>& Rmult4by4(const Matrix<T>& mat) {
+				Matrix temp(*this);
+
+				m0 = mat.m0*temp.m0+mat.m4*temp.m1+mat.m8*temp.m2+mat.m12*temp.m3;
+				m4 = mat.m0*temp.m4+mat.m4*temp.m5+mat.m8*temp.m6+mat.m12*temp.m7;
+				m8 = mat.m0*temp.m8+mat.m4*temp.m9+mat.m8*temp.m10+mat.m12*temp.m11;
+				m12 = mat.m0*temp.m12+mat.m4*temp.m13+mat.m8*temp.m14+mat.m12*temp.m15;
+
+				m1 =  mat.m1*temp.m0 + mat.m5*temp.m1  + mat.m9*temp.m2+mat.m13*temp.m3;
+				m5 =  mat.m1*temp.m4 + mat.m5*temp.m5  + mat.m9*temp.m6+mat.m13*temp.m7;
+				m9 =  mat.m1*temp.m8 + mat.m5*temp.m9  + mat.m9*temp.m10+mat.m13*temp.m11;
+				m13 = mat.m1*temp.m12+ mat.m5*temp.m13 + mat.m9*temp.m14+mat.m13*temp.m15;
+
+				m2 = mat.m2*temp.m0+mat.m6*temp.m1+mat.m10*temp.m2+mat.m14*temp.m3;
+				m6 = mat.m2*temp.m4+mat.m6*temp.m5+mat.m10*temp.m6+mat.m14*temp.m7;
+				m10 = mat.m2*temp.m8+mat.m6*temp.m9+mat.m10*temp.m10+mat.m14*temp.m11;
+				m14 = mat.m2*temp.m12+mat.m6*temp.m13+mat.m10*temp.m14+mat.m14*temp.m15;
+
+				m3 = mat.m3*temp.m0+mat.m7*temp.m1+mat.m11*temp.m2+mat.m15*temp.m3;
+				m7 = mat.m3*temp.m4+mat.m7*temp.m5+mat.m11*temp.m6+mat.m15*temp.m7;
+				m11 = mat.m3*temp.m8+mat.m7*temp.m9+mat.m11*temp.m10+mat.m15*temp.m11;
+				m15 = mat.m3*temp.m12+mat.m7*temp.m13+mat.m11*temp.m14+mat.m15*temp.m15;
+				return *this;
+			}
+
+
+			inline Matrix<T>& mult4by4(const Matrix<T>& mat) {
+				Matrix temp(*this);
+
+				m0 = temp.m0*mat.m0+temp.m4*mat.m1+temp.m8*mat.m2+temp.m12*mat.m3;
+				m4 = temp.m0*mat.m4+temp.m4*mat.m5+temp.m8*mat.m6+temp.m12*mat.m7;
+				m8 = temp.m0*mat.m8+temp.m4*mat.m9+temp.m8*mat.m10+temp.m12*mat.m11;
+				m12 = temp.m0*mat.m12+temp.m4*mat.m13+temp.m8*mat.m14+temp.m12*mat.m15;
+
+				m1 =  temp.m1*mat.m0 + temp.m5*mat.m1  + temp.m9*mat.m2+temp.m13*mat.m3;
+				m5 =  temp.m1*mat.m4 + temp.m5*mat.m5  + temp.m9*mat.m6+temp.m13*mat.m7;
+				m9 =  temp.m1*mat.m8 + temp.m5*mat.m9  + temp.m9*mat.m10+temp.m13*mat.m11;
+				m13 = temp.m1*mat.m12+ temp.m5*mat.m13 + temp.m9*mat.m14+temp.m13*mat.m15;
+
+				m2 = temp.m2*mat.m0+temp.m6*mat.m1+temp.m10*mat.m2+temp.m14*mat.m3;
+				m6 = temp.m2*mat.m4+temp.m6*mat.m5+temp.m10*mat.m6+temp.m14*mat.m7;
+				m10 = temp.m2*mat.m8+temp.m6*mat.m9+temp.m10*mat.m10+temp.m14*mat.m11;
+				m14 = temp.m2*mat.m12+temp.m6*mat.m13+temp.m10*mat.m14+temp.m14*mat.m15;
+
+				m3 = temp.m3*mat.m0+temp.m7*mat.m1+temp.m11*mat.m2+temp.m15*mat.m3;
+				m7 = temp.m3*mat.m4+temp.m7*mat.m5+temp.m11*mat.m6+temp.m15*mat.m7;
+				m11 = temp.m3*mat.m8+temp.m7*mat.m9+temp.m11*mat.m10+temp.m15*mat.m11;
+				m15 = temp.m3*mat.m12+temp.m7*mat.m13+temp.m11*mat.m14+temp.m15*mat.m15;
+				return *this;
+			}
+
+			Matrix& applyRotate(T angle, T x, T y, T z) {
+				static Matrix<T> temp;
+				temp.loadRotate(angle,x,y,z);
+				*this = temp.mult4by4(*this);
+				return  *this;
+			}
+
+
+		private:
+#define cofactor_maker(f1,mj1,mi1, f2,mj2,mi2, f3,mj3,mi3) \
+		f1*(mj1*mi1-mj2*mi3) + f2*(mj2*mi2-mj3*mi1) + f3*(mj3*mi3-mj1*mi2)
+
+			inline T cofactorm0() const {
+				return cofactor_maker(m5,m10,m15, m6,m11,m13, m7,m9,m14);
+			}
+			inline T cofactorm1() const {
+				return cofactor_maker(m6,m11,m12, m7,m8,m14, m4,m10,m15);
+			}
+			inline T cofactorm2() const {
+				return cofactor_maker(m7,m8,m13, m4,m9,m15, m5,m11,m12);
+			}
+			inline T cofactorm3() const {
+				return cofactor_maker(m4,m9,m14, m5,m10,m12, m6,m8,m13);
+			}
+			inline T cofactorm4() const {
+				return cofactor_maker(m9,m14,m3, m10,m15,m1, m11,m13,m2);
+			}
+			inline T cofactorm5() const {
+				return cofactor_maker(m10,m15,m0, m11,m12,m2, m8,m14,m3);
+			}
+			inline T cofactorm6() const {
+				return cofactor_maker(m11,m12,m1, m8,m13,m3, m9,m15,m0);
+			}
+			inline T cofactorm7() const {
+				return cofactor_maker(m8,m13,m2, m9,m14,m0, m10,m12,m1);
+			}
+			inline T cofactorm8() const {
+				return cofactor_maker(m13,m2,m7, m14,m3,m5, m15,m1,m6);
+			}
+			inline T cofactorm9() const {
+				return cofactor_maker(m14,m13,m4, m15,m0,m6, m12,m2,m7);
+			}
+			inline T cofactorm10() const {
+				return cofactor_maker(m15,m0,m5, m12,m1,m7, m13,m3,m4);
+			}
+			inline T cofactorm11() const {
+				return cofactor_maker(m12,m1,m6, m13,m2,m4, m14,m0,m5);
+			}
+			inline T cofactorm12() const {
+				return cofactor_maker(m1,m6,m11, m2,m7,m9, m3,m5,m10);
+			}
+			inline T cofactorm13() const {
+				return cofactor_maker(m2,m7,m8, m3,m4,m10, m10,m6,m11);
+			}
+			inline T cofactorm14() const {
+				return cofactor_maker(m3,m4,m9, m0,m5,m11, m1,m7,m8);
+			}
+			inline T cofactorm15() const {
+				return cofactor_maker(m0,m5,m10, m1,m6,m8, m2,m4,m9);
+			}
+
+			union {
+				T m[16];
+				struct {
+					T m0,m1,m2,m3,m4,m5,m6,m7,m8,m9,m10,m11,m12,m13,m14,m15;
+				};
+			};
+	};
+
+	typedef Matrix<double> DoubleMatrix;
+	typedef Matrix<int> IntMatrix;
+
+	/** Print coords of the Matrix to a stream
+	 */
+	template<typename T>
+	std::ostream& operator<<(std::ostream& os, Matrix<T>& m) {
+		return os << "\n|" << m[0] << "," << m[4] << "," << m[8] << ","  << m[12] << "|\n" << \
+		               "|" << m[1] << "," << m[5] << "," << m[9] << ","  << m[13] << "|\n" << \
+		               "|" << m[2] << "," << m[6] << "," << m[10] << "," << m[14] << "|\n" << \
+		               "|" << m[3] << "," << m[7] << "," << m[11] << "," << m[15] << "|\n";
+	}
+
+
+}
+
+#endif