aboutsummaryrefslogtreecommitdiff
path: root/dep/include/g3dlite/G3D/Matrix.h
diff options
context:
space:
mode:
Diffstat (limited to 'dep/include/g3dlite/G3D/Matrix.h')
-rw-r--r--dep/include/g3dlite/G3D/Matrix.h634
1 files changed, 634 insertions, 0 deletions
diff --git a/dep/include/g3dlite/G3D/Matrix.h b/dep/include/g3dlite/G3D/Matrix.h
new file mode 100644
index 00000000000..3c5394d9a76
--- /dev/null
+++ b/dep/include/g3dlite/G3D/Matrix.h
@@ -0,0 +1,634 @@
+/**
+ @file Matrix.h
+ @author Morgan McGuire, http://graphics.cs.williams.edu
+
+ @created 2005-10-23
+ @edited 2007-07-18
+ */
+
+#ifndef G3D_MATRIX_H
+#define G3D_MATRIX_H
+
+#include "G3D/g3dmath.h"
+#include "G3D/Vector3.h"
+#include "G3D/Vector4.h"
+#include "G3D/Matrix3.h"
+#include "G3D/Matrix4.h"
+#include "G3D/ReferenceCount.h"
+
+namespace G3D {
+
+/**
+ N x M matrix.
+
+ The actual data is tracked internally by a reference counted pointer;
+ it is efficient to pass and assign Matrix objects because no data is actually copied.
+ This avoids the headache of pointers and allows natural math notation:
+
+ <PRE>
+ Matrix A, B, C;
+ // ...
+
+ C = A * f(B);
+ C = C.inverse();
+
+ A = Matrix::identity(4);
+ C = A;
+ C.set(0, 0, 2.0); // Triggers a copy of the data so that A remains unchanged.
+
+ // etc.
+
+ </PRE>
+
+ The Matrix::debugNumCopyOps and Matrix::debugNumAllocOps counters
+ increment every time an operation forces the copy and allocation of matrices. You
+ can use these to detect slow operations when efficiency is a major concern.
+
+ Some methods accept an output argument instead of returning a value. For example,
+ <CODE>A = B.transpose()</CODE> can also be invoked as <CODE>B.transpose(A)</CODE>.
+ The latter may be more efficient, since Matrix may be able to re-use the storage of
+ A (if it has approximatly the right size and isn't currently shared with another matrix).
+
+ @sa G3D::Matrix3, G3D::Matrix4, G3D::Vector2, G3D::Vector3, G3D::Vector4, G3D::CoordinateFrame
+
+ @beta
+ */
+class Matrix {
+public:
+ /**
+ Internal precision. Currently float, but this may become a templated class in the future
+ to allow operations like Matrix<double> and Matrix<ComplexFloat>.
+
+ Not necessarily a plain-old-data type (e.g., could ComplexFloat), but must be something
+ with no constructor, that can be safely memcpyd, and that has a bit pattern of all zeros
+ when zero.*/
+ typedef float T;
+
+ /** Incremented every time the elements of a matrix are copied. Useful for profiling your
+ own code that uses Matrix to determine when it is slow due to copying.*/
+ static int debugNumCopyOps;
+
+ /** Incremented every time a new matrix object is allocated. Useful for profiling your
+ own code that uses Matrix to determine when it is slow due to allocation.*/
+ static int debugNumAllocOps;
+
+private:
+public:
+
+ /** Used internally by Matrix.
+
+ Does not throw exceptions-- assumes the caller has taken care of
+ argument checking. */
+ class Impl : public ReferenceCountedObject {
+ public:
+
+ static void* operator new(size_t size) {
+ return System::malloc(size);
+ }
+
+ static void operator delete(void* p) {
+ System::free(p);
+ }
+
+ ~Impl();
+
+ private:
+ friend class Matrix;
+
+ /** elt[r][c] = the element. Pointers into data.*/
+ T** elt;
+
+ /** Row major data for the entire matrix. */
+ T* data;
+
+ /** The number of rows */
+ int R;
+
+ /** The number of columns */
+ int C;
+
+ int dataSize;
+
+ /** If R*C is much larger or smaller than the current, deletes all previous data
+ and resets to random data. Otherwise re-uses existing memory and just resets
+ R, C, and the row pointers. */
+ void setSize(int newRows, int newCols);
+
+ inline Impl() : elt(NULL), data(NULL), R(0), C(0), dataSize(0) {}
+
+ Impl(const Matrix3& M);
+
+ Impl(const Matrix4& M);
+
+ inline Impl(int r, int c) : elt(NULL), data(NULL), R(0), C(0), dataSize(0) {
+ setSize(r, c);
+ }
+
+ Impl& operator=(const Impl& m);
+
+ inline Impl(const Impl& B) : elt(NULL), data(NULL), R(0), C(0), dataSize(0) {
+ // Use the assignment operator
+ *this = B;
+ }
+
+ void setZero();
+
+ inline void set(int r, int c, T v) {
+ debugAssert(r < R);
+ debugAssert(c < C);
+ elt[r][c] = v;
+ }
+
+ inline const T& get(int r, int c) const {
+ debugAssert(r < R);
+ debugAssert(c < C);
+ return elt[r][c];
+ }
+
+ /** Multiplies this by B and puts the result in out. */
+ void mul(const Impl& B, Impl& out) const;
+
+ /** Ok if out == this or out == B */
+ void add(const Impl& B, Impl& out) const;
+
+ /** Ok if out == this or out == B */
+ void add(T B, Impl& out) const;
+
+ /** Ok if out == this or out == B */
+ void sub(const Impl& B, Impl& out) const;
+
+ /** Ok if out == this or out == B */
+ void sub(T B, Impl& out) const;
+
+ /** B - this */
+ void lsub(T B, Impl& out) const;
+
+ /** Ok if out == this or out == B */
+ void arrayMul(const Impl& B, Impl& out) const;
+
+ /** Ok if out == this or out == B */
+ void mul(T B, Impl& out) const;
+
+ /** Ok if out == this or out == B */
+ void arrayDiv(const Impl& B, Impl& out) const;
+
+ /** Ok if out == this or out == B */
+ void div(T B, Impl& out) const;
+
+ void negate(Impl& out) const;
+
+ /** Slow way of computing an inverse; for reference */
+ void inverseViaAdjoint(Impl& out) const;
+
+ /** Use Gaussian elimination with pivots to solve for the inverse destructively in place. */
+ void inverseInPlaceGaussJordan();
+
+ void adjoint(Impl& out) const;
+
+ /** Matrix of all cofactors */
+ void cofactor(Impl& out) const;
+
+ /**
+ Cofactor [r][c] is defined as C[r][c] = -1 ^(r+c) * det(A[r][c]),
+ where A[r][c] is the (R-1)x(C-1) matrix formed by removing row r and
+ column c from the original matrix.
+ */
+ T cofactor(int r, int c) const;
+
+ /** Ok if out == this or out == B */
+ void transpose(Impl& out) const;
+
+ T determinant() const;
+
+ /** Determinant computed without the given row and column */
+ T determinant(int r, int c) const;
+
+ void arrayLog(Impl& out) const;
+
+ void arrayExp(Impl& out) const;
+
+ void arraySqrt(Impl& out) const;
+
+ void arrayCos(Impl& out) const;
+
+ void arraySin(Impl& out) const;
+
+ void swapRows(int r0, int r1);
+
+ void swapAndNegateCols(int c0, int c1);
+
+ void mulRow(int r, const T& v);
+
+ void abs(Impl& out) const;
+
+ /** Makes a (R-1)x(C-1) copy of this matrix */
+ void withoutRowAndCol(int excludeRow, int excludeCol, Impl& out) const;
+
+ bool anyNonZero() const;
+
+ bool allNonZero() const;
+
+ void setRow(int r, const T* vals);
+
+ void setCol(int c, const T* vals);
+ };
+private:
+
+ typedef ReferenceCountedPointer<Impl> ImplRef;
+
+ ImplRef impl;
+
+ inline Matrix(ImplRef i) : impl(i) {}
+ inline Matrix(Impl* i) : impl(ImplRef(i)) {}
+
+ /** Used by SVD */
+ class SortRank {
+ public:
+ T value;
+ int col;
+
+ inline bool operator>(const SortRank& x) const {
+ return x.value > value;
+ }
+
+ inline bool operator<(const SortRank& x) const {
+ return x.value < value;
+ }
+
+ inline bool operator>=(const SortRank& x) const {
+ return x.value >= value;
+ }
+
+ inline bool operator<=(const SortRank& x) const {
+ return x.value <= value;
+ }
+
+ inline bool operator==(const SortRank& x) const {
+ return x.value == value;
+ }
+
+ inline bool operator!=(const SortRank& x) const {
+ return x.value != value;
+ }
+ };
+
+ Matrix vectorPseudoInverse() const;
+ Matrix partitionPseudoInverse() const;
+ Matrix colPartPseudoInverse() const;
+ Matrix rowPartPseudoInverse() const;
+
+ Matrix col2PseudoInverse(const Matrix& B) const;
+ Matrix col3PseudoInverse(const Matrix& B) const;
+ Matrix col4PseudoInverse(const Matrix& B) const;
+ Matrix row2PseudoInverse(const Matrix& B) const;
+ Matrix row3PseudoInverse(const Matrix& B) const;
+ Matrix row4PseudoInverse(const Matrix& B) const;
+
+public:
+
+ Matrix() : impl(new Impl(0, 0)) {}
+
+ Matrix(const Matrix3& M) : impl(new Impl(M)) {}
+
+ Matrix(const Matrix4& M) : impl(new Impl(M)) {}
+
+ template<class S>
+ static Matrix fromDiagonal(const Array<S>& d) {
+ Matrix D = zero(d.length(), d.length());
+ for (int i = 0; i < d.length(); ++i) {
+ D.set(i, i, d[i]);
+ }
+ return D;
+ }
+
+ static Matrix fromDiagonal(const Matrix& d);
+
+ /** Returns a new matrix that is all zero. */
+ Matrix(int R, int C) : impl(new Impl(R, C)) {
+ impl->setZero();
+ }
+
+ /** Returns a new matrix that is all zero. */
+ static Matrix zero(int R, int C);
+
+ /** Returns a new matrix that is all one. */
+ static Matrix one(int R, int C);
+
+ /** Returns a new identity matrix */
+ static Matrix identity(int N);
+
+ /** Uniformly distributed values between zero and one. */
+ static Matrix random(int R, int C);
+
+ /** The number of rows */
+ inline int rows() const {
+ return impl->R;
+ }
+
+ /** Number of columns */
+ inline int cols() const {
+ return impl->C;
+ }
+
+ /** Generally more efficient than A * B */
+ Matrix& operator*=(const T& B);
+
+ /** Generally more efficient than A / B */
+ Matrix& operator/=(const T& B);
+
+ /** Generally more efficient than A + B */
+ Matrix& operator+=(const T& B);
+
+ /** Generally more efficient than A - B */
+ Matrix& operator-=(const T& B);
+
+ /** No performance advantage over A * B because
+ matrix multiplication requires intermediate
+ storage. */
+ Matrix& operator*=(const Matrix& B);
+
+ /** Generally more efficient than A + B */
+ Matrix& operator+=(const Matrix& B);
+
+ /** Generally more efficient than A - B */
+ Matrix& operator-=(const Matrix& B);
+
+ /** Returns a new matrix that is a subset of this one,
+ from r1:r2 to c1:c2, inclusive.*/
+ Matrix subMatrix(int r1, int r2, int c1, int c2) const;
+
+ /** Matrix multiplication. To perform element-by-element multiplication,
+ see arrayMul. */
+ inline Matrix operator*(const Matrix& B) const {
+ Matrix C(impl->R, B.impl->C);
+ impl->mul(*B.impl, *C.impl);
+ return C;
+ }
+
+ /** See also A *= B, which is more efficient in many cases */
+ inline Matrix operator*(const T& B) const {
+ Matrix C(impl->R, impl->C);
+ impl->mul(B, *C.impl);
+ return C;
+ }
+
+ /** See also A += B, which is more efficient in many cases */
+ inline Matrix operator+(const Matrix& B) const {
+ Matrix C(impl->R, impl->C);
+ impl->add(*B.impl, *C.impl);
+ return C;
+ }
+
+ /** See also A -= B, which is more efficient in many cases */
+ inline Matrix operator-(const Matrix& B) const {
+ Matrix C(impl->R, impl->C);
+ impl->sub(*B.impl, *C.impl);
+ return C;
+ }
+
+ /** See also A += B, which is more efficient in many cases */
+ inline Matrix operator+(const T& v) const {
+ Matrix C(impl->R, impl->C);
+ impl->add(v, *C.impl);
+ return C;
+ }
+
+ /** See also A -= B, which is more efficient in many cases */
+ inline Matrix operator-(const T& v) const {
+ Matrix C(impl->R, impl->C);
+ impl->sub(v, *C.impl);
+ return C;
+ }
+
+
+ Matrix operator>(const T& scalar) const;
+
+ Matrix operator<(const T& scalar) const;
+
+ Matrix operator>=(const T& scalar) const;
+
+ Matrix operator<=(const T& scalar) const;
+
+ Matrix operator==(const T& scalar) const;
+
+ Matrix operator!=(const T& scalar) const;
+
+ /** scalar B - this */
+ inline Matrix lsub(const T& B) const {
+ Matrix C(impl->R, impl->C);
+ impl->lsub(B, *C.impl);
+ return C;
+ }
+
+ inline Matrix arrayMul(const Matrix& B) const {
+ Matrix C(impl->R, impl->C);
+ impl->arrayMul(*B.impl, *C.impl);
+ return C;
+ }
+
+ Matrix3 toMatrix3() const;
+
+ Matrix4 toMatrix4() const;
+
+ Vector2 toVector2() const;
+
+ Vector3 toVector3() const;
+
+ Vector4 toVector4() const;
+
+ /** Mutates this */
+ void arrayMulInPlace(const Matrix& B);
+
+ /** Mutates this */
+ void arrayDivInPlace(const Matrix& B);
+
+ // Declares an array unary method and its explicit-argument counterpart
+# define DECLARE_METHODS_1(method)\
+ inline Matrix method() const {\
+ Matrix C(impl->R, impl->C);\
+ impl->method(*C.impl);\
+ return C;\
+ }\
+ void method(Matrix& out) const;
+
+
+ DECLARE_METHODS_1(abs)
+ DECLARE_METHODS_1(arrayLog)
+ DECLARE_METHODS_1(arrayExp)
+ DECLARE_METHODS_1(arraySqrt)
+ DECLARE_METHODS_1(arrayCos)
+ DECLARE_METHODS_1(arraySin)
+ DECLARE_METHODS_1(negate)
+
+# undef DECLARE_METHODS_1
+
+ inline Matrix operator-() const {
+ return negate();
+ }
+
+ /**
+ A<SUP>-1</SUP> computed using the Gauss-Jordan algorithm,
+ for square matrices.
+ Run time is <I>O(R<sup>3</sup>)</I>, where <I>R</i> is the
+ number of rows.
+ */
+ inline Matrix inverse() const {
+ Impl* A = new Impl(*impl);
+ A->inverseInPlaceGaussJordan();
+ return Matrix(A);
+ }
+
+ inline T determinant() const {
+ return impl->determinant();
+ }
+
+ /**
+ A<SUP>T</SUP>
+ */
+ inline Matrix transpose() const {
+ Impl* A = new Impl(cols(), rows());
+ impl->transpose(*A);
+ return Matrix(A);
+ }
+
+ /** Transpose in place; more efficient than transpose */
+ void transpose(Matrix& out) const;
+
+ inline Matrix adjoint() const {
+ Impl* A = new Impl(cols(), rows());
+ impl->adjoint(*A);
+ return Matrix(A);
+ }
+
+ /**
+ (A<SUP>T</SUP>A)<SUP>-1</SUP>A<SUP>T</SUP>) computed
+ using SVD.
+
+ @param tolerance Use -1 for automatic tolerance.
+ */
+ Matrix pseudoInverse(float tolerance = -1) const;
+
+ /** Called from pseudoInverse when the matrix has size > 4 along some dimension.*/
+ Matrix svdPseudoInverse(float tolerance = -1) const;
+
+ /**
+ (A<SUP>T</SUP>A)<SUP>-1</SUP>A<SUP>T</SUP>) computed
+ using Gauss-Jordan elimination.
+ */
+ inline Matrix gaussJordanPseudoInverse() const {
+ Matrix trans = transpose();
+ return (trans * (*this)).inverse() * trans;
+ }
+
+ /** Singular value decomposition. Factors into three matrices
+ such that @a this = @a U * fromDiagonal(@a d) * @a V.transpose().
+
+ The matrix must have at least as many rows as columns.
+
+ Run time is <I>O(C<sup>2</sup>*R)</I>.
+
+ @param sort If true (default), the singular values
+ are arranged so that D is sorted from largest to smallest.
+ */
+ void svd(Matrix& U, Array<T>& d, Matrix& V, bool sort = true) const;
+
+ void set(int r, int c, T v);
+
+ void setCol(int c, const Matrix& vec);
+
+ void setRow(int r, const Matrix& vec);
+
+ Matrix col(int c) const;
+
+ Matrix row(int r) const;
+
+ T get(int r, int c) const;
+
+ Vector2int16 size() const {
+ return Vector2int16(rows(), cols());
+ }
+
+ int numElements() const {
+ return rows() * cols();
+ }
+
+ void swapRows(int r0, int r1);
+
+ /** Swaps columns c0 and c1 and negates both */
+ void swapAndNegateCols(int c0, int c1);
+
+ void mulRow(int r, const T& v);
+
+ /** Returns true if any element is non-zero */
+ bool anyNonZero() const;
+
+ /** Returns true if all elements are non-zero */
+ bool allNonZero() const;
+
+ inline bool allZero() const {
+ return !anyNonZero();
+ }
+
+ inline bool anyZero() const {
+ return !allNonZero();
+ }
+
+ /** Serializes in Matlab source format */
+ void serialize(TextOutput& t) const;
+
+ std::string toString(const std::string& name) const;
+
+ std::string toString() const {
+ static const std::string name = "";
+ return toString(name);
+ }
+
+ /** 2-norm squared: sum(squares). (i.e., dot product with itself) */
+ double normSquared() const;
+
+ /** 2-norm (sqrt(sum(squares)) */
+ double norm() const;
+
+ /**
+ Low-level SVD functionality. Useful for applications that do not want
+ to construct a Matrix but need to perform the SVD operation.
+
+ this = U * D * V'
+
+ Assumes that rows >= cols
+
+ @return NULL on success, a string describing the error on failure.
+ @param U rows x cols matrix to be decomposed, gets overwritten with U, a rows x cols matrix with orthogonal columns.
+ @param D vector of singular values of a (diagonal of the D matrix). Length cols.
+ @param V returns the right orthonormal transformation matrix, size cols x cols
+
+ @cite Based on Dianne Cook's implementation, which is adapted from
+ svdecomp.c in XLISP-STAT 2.1, which is code from Numerical Recipes
+ adapted by Luke Tierney and David Betz. The Numerical Recipes code
+ is adapted from Forsythe et al, who based their code on Golub and
+ Reinsch's original implementation.
+ */
+ static const char* svdCore(float** U, int rows, int cols, float* D, float** V);
+
+};
+
+}
+
+inline G3D::Matrix operator-(const G3D::Matrix::T& v, const G3D::Matrix& M) {
+ return M.lsub(v);
+}
+
+inline G3D::Matrix operator*(const G3D::Matrix::T& v, const G3D::Matrix& M) {
+ return M * v;
+}
+
+inline G3D::Matrix operator+(const G3D::Matrix::T& v, const G3D::Matrix& M) {
+ return M + v;
+}
+
+inline G3D::Matrix abs(const G3D::Matrix& M) {
+ return M.abs();
+}
+
+#endif
+