VRS - The Virtual Rendering System

version 3.3

mat.h

Go to the documentation of this file.

/******************************************************************************
 * VRS - The Virtual Rendering System
 * Copyright (C) 2000-2008 Computer Graphics Systems Group at the 
 * Hasso-Plattner-Institute (HPI), Potsdam, Germany.
 *
 * This  library  is  free  software; you can redistribute it and/or  modify it
 * under the terms of the GNU Lesser General Public  License  as  published  by
 * the  Free  Software  Foundation;  either  version 2.1 of the License, or (at
 * your option) any later version. This library  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 Lesser General Public License for more details. You should have received
 * a copy of the GNU Lesser+ General Public License along with this library; if
 * not, write to the FreeSoftware Foundation, Inc., 59 Temple Place, Suite 330,
 * Boston, MA, 02111-1307, USA.
 ******************************************************************************/

#ifndef VRS_MAT_H
#define VRS_MAT_H

#include <vrs/vec.h>

namespace VRS {

    template<unsigned int DIM, typename T>
    class Mat {
    public:
        typedef T value_type;

    public:
        explicit Mat(T diag = T(1)) { for(size_t i = 0; i < DIM; ++i) { m_columns[i][i] = diag; } }
        explicit Mat(Vec<DIM, T> const& diag) { for(size_t i = 0; i < DIM; ++i) { m_columns[i][i] = diag[i]; } }

    public:
        Mat(Vec<4, T> const& column0, Vec<4, T> const& column1, Vec<4, T> const& column2, Vec<4, T> const& column3) {
            m_columns[0] = column0; m_columns[1] = column1; m_columns[2] = column2; m_columns[3] = column3;
            VRS_STATIC_ASSERT(DIM == 4);
        }

        Mat(
            T c0r0, T c0r1, T c0r2, T c0r3, // first column
            T c1r0, T c1r1, T c1r2, T c1r3, // second column
            T c2r0, T c2r1, T c2r2, T c2r3, // third column
            T c3r0, T c3r1, T c3r2, T c3r3  // fourth column
        ) {
            m_columns[0] = Vec<4, T>(c0r0, c0r1, c0r2, c0r3);
            m_columns[1] = Vec<4, T>(c1r0, c1r1, c1r2, c1r3);
            m_columns[2] = Vec<4, T>(c2r0, c2r1, c2r2, c2r3);
            m_columns[3] = Vec<4, T>(c3r0, c3r1, c3r2, c3r3);
            VRS_STATIC_ASSERT(DIM == 4);
        }

    public:
        Mat(Vec<3, T> const& column0, Vec<3, T> const& column1, Vec<3, T> const& column2) {
            m_columns[0] = column0; m_columns[1] = column1; m_columns[2] = column2;
            VRS_STATIC_ASSERT(DIM == 3);
        }

        Mat(
            T c0r0, T c0r1, T c0r2, // first column
            T c1r0, T c1r1, T c1r2, // second column
            T c2r0, T c2r1, T c2r2  // third column
        ) {
            m_columns[0] = Vec<3, T>(c0r0, c0r1, c0r2);
            m_columns[1] = Vec<3, T>(c1r0, c1r1, c1r2);
            m_columns[2] = Vec<3, T>(c2r0, c2r1, c2r2);
            VRS_STATIC_ASSERT(DIM == 3);
        }

    public:
        Mat(Vec<2, T> const& column0, Vec<2, T> const& column1) {
            m_columns[0] = column0; m_columns[1] = column1;
            VRS_STATIC_ASSERT(DIM == 2);
        }

        Mat(
            T c0r0, T c0r1, // first column
            T c1r0, T c1r1  // second column
        ) {
            m_columns[0] = Vec<2, T>(c0r0, c0r1);
            m_columns[1] = Vec<2, T>(c1r0, c1r1);
            VRS_STATIC_ASSERT(DIM == 2);
        }

    public:
        template<unsigned int DIM2, typename S>
        explicit Mat(Mat<DIM2, S> const& other) {
            for(size_t i = 0; i < DIM; ++i) {
                for(size_t j = 0; j < DIM; ++j) {
                    if((i < DIM2) && (j < DIM2)) {
                        m_columns[i][j] = other[i][j];
                    } else {
                        m_columns[i][j] = (i == j) ? T(1) : T(0);
                    }
                }
            }
        }

    public:
        Vec<DIM, T> const& operator[](size_t column) const { assert(column < DIM); return m_columns[column]; }
        Vec<DIM, T>      & operator[](size_t column)       { assert(column < DIM); return m_columns[column]; }

        T  element(size_t column, size_t row) const { assert(column < DIM); assert(row < DIM); return m_columns[column][row]; }
        T& element(size_t column, size_t row)       { assert(column < DIM); assert(row < DIM); return m_columns[column][row]; }

        Vec<DIM, T> const& column(size_t column) const { assert(column < DIM); return m_columns[column]; }
        Vec<DIM, T> row(size_t row) const { assert(row < DIM); Vec<DIM, T> result; for(size_t i = 0; i < DIM; ++i) { result[i] = m_columns[i][row]; } return result; }

        Mat<DIM - 1, T> minor(size_t column, size_t row) const { assert(column < DIM); assert(row < DIM); Mat<DIM - 1, T> result; for(size_t i = 0; i < DIM - 1; ++i) { for(size_t j = 0; j < DIM - 1; ++j) { result[i][j] = m_columns[(i < column) ? i : i+1][(j < row) ? j : j+1]; } } return result; }

    public:
        bool operator==(Mat<DIM, T> const& other) const { for(size_t i = 0; i < DIM; ++i) { if(m_columns[i] != other[i]) { return false; } } return true; }
        bool operator<(Mat<DIM, T> const& other) const { for(size_t i = 0; i < DIM; ++i) { if(m_columns[i] < other[i]) { return true; } if(m_columns[i] > other[i]) { return false; } } return false; }

    public:
        Mat<DIM, T>& operator+=(Mat<DIM, T> const& other) { for(size_t i = 0; i < DIM; ++i) { m_columns[i] += other[i]; } return *this; }
        Mat<DIM, T>& operator-=(Mat<DIM, T> const& other) { for(size_t i = 0; i < DIM; ++i) { m_columns[i] -= other[i]; } return *this; }
        Mat<DIM, T>& operator*=(Mat<DIM, T> const& mat) { return operator=(*this * mat); }

    private:
        Vec<DIM, T> m_columns[DIM];
    }; // class Mat<DIM, T>

    template<unsigned int DIM, typename T>
    inline std::ostream& operator<<(std::ostream& out, const Mat<DIM, T>& mat) {
        for(size_t i = 0; i < DIM; ++i) {
            out << ((i > 0) ? ",   " : "") << mat[i];
        }
        return out;
    }

    template<unsigned int DIM, typename T>
    bool operator!=(Mat<DIM, T> const& m0, Mat<DIM, T> const& m1) { return !(m0 == m1); }

    template<unsigned int DIM, typename T>
    bool operator>(Mat<DIM, T> const& m0, Mat<DIM, T> const& m1) { return (m1 < m0); }
    template<unsigned int DIM, typename T>
    bool operator<=(Mat<DIM, T> const& m0, Mat<DIM, T> const& m1) { return !(m0 > m1); }
    template<unsigned int DIM, typename T>
    bool operator>=(Mat<DIM, T> const& m0, Mat<DIM, T> const& m1) { return !(m0 < m1); }

    template<unsigned int DIM, typename T>
    inline Mat<DIM, T> operator+(Mat<DIM, T> const& m0, Mat<DIM, T> const& m1) { Mat<DIM, T> result = m0; return (result += m1); }
    template<unsigned int DIM, typename T>
    inline Mat<DIM, T> operator-(Mat<DIM, T> const& m0, Mat<DIM, T> const& m1) { Mat<DIM, T> result = m0; return (result -= m1); }
    template<unsigned int DIM, typename T>
    inline Mat<DIM, T> operator*(Mat<DIM, T> const& m0, Mat<DIM, T> const& m1) {
        Mat<DIM, T> result;
        for(size_t j = 0; j < DIM; ++j) {
            const Vec<DIM, T> m0rowJ(m0.row(j));
            for(size_t i = 0; i < DIM; ++i) {
                result[i][j] = dot(m0rowJ, m1.column(i));
            }
        }
        return result;
    }

    template<unsigned int DIM, typename T>
    inline Mat<DIM, T> operator+(Mat<DIM, T> const& m0, T s1) { Mat<DIM, T> result; for(size_t i = 0; i < DIM; ++i) { for(size_t j = 0; j < DIM; ++j) { result[i][j] = m0[i][j] + s1; } } return result; }
    template<unsigned int DIM, typename T>
    inline Mat<DIM, T> operator-(Mat<DIM, T> const& m0, T s1) { Mat<DIM, T> result; for(size_t i = 0; i < DIM; ++i) { for(size_t j = 0; j < DIM; ++j) { result[i][j] = m0[i][j] - s1; } } return result; }
    template<unsigned int DIM, typename T>
    inline Mat<DIM, T> operator*(Mat<DIM, T> const& m0, T s1) { Mat<DIM, T> result; for(size_t i = 0; i < DIM; ++i) { for(size_t j = 0; j < DIM; ++j) { result[i][j] = m0[i][j] * s1; } } return result; }

    template<unsigned int DIM, typename T>
    inline Mat<DIM, T> operator+(T s0, Mat<DIM, T> const& m1) { Mat<DIM, T> result; for(size_t i = 0; i < DIM; ++i) { for(size_t j = 0; j < DIM; ++j) { result[i][j] = s0 + m1[i][j]; } } return result; }
    template<unsigned int DIM, typename T>
    inline Mat<DIM, T> operator-(T s0, Mat<DIM, T> const& m1) { Mat<DIM, T> result; for(size_t i = 0; i < DIM; ++i) { for(size_t j = 0; j < DIM; ++j) { result[i][j] = s0 - m1[i][j]; } } return result; }
    template<unsigned int DIM, typename T>
    inline Mat<DIM, T> operator*(T s0, Mat<DIM, T> const& m1) { Mat<DIM, T> result; for(size_t i = 0; i < DIM; ++i) { for(size_t j = 0; j < DIM; ++j) { result[i][j] = s0 * m1[i][j]; } } return result; }

    template<unsigned int DIM, typename T>
    inline Mat<DIM, T> operator*(Mat<DIM, T> const& mat) { Mat<DIM, T> result; for(size_t i = 0; i < DIM; ++i) { for(size_t j = 0; j < DIM; ++j) { result[i][j] = -mat[i][j]; } } return result; }

    template<unsigned int DIM, typename T>
    inline Mat<DIM, T>& operator++(Mat<DIM, T>& v) { return (v += T(1)); }
    template<unsigned int DIM, typename T>
    inline Mat<DIM, T> operator++(Mat<DIM, T>& v, int) { const Mat<DIM, T> result = v; v += T(1); return result; }
    template<unsigned int DIM, typename T>
    inline Mat<DIM, T>& operator--(Mat<DIM, T>& v) { return (v -= T(1)); }
    template<unsigned int DIM, typename T>
    inline Mat<DIM, T> operator--(Mat<DIM, T>& v, int) { const Mat<DIM, T> result = v; v -= T(1); return result; }

    template<unsigned int DIM, typename T>
    inline Mat<DIM, T> transpose(Mat<DIM, T> const& x) { Mat<DIM, T> result; for(size_t i = 0; i < DIM; ++i) { for(size_t j = 0; j < DIM; ++j) { result[i][j] = x[j][i]; } } return result; }

    template<typename T>
    inline T det(Mat<1, T> const& x) { return x[0][0]; }

    template<typename T>
    inline T det(Mat<2, T> const& x) { return (x[0][0]*x[1][1] - x[1][0]*x[0][1]); }

    template<typename T>
    inline T det(Mat<3, T> const& x) { return (x[0][0]*x[1][1]*x[2][2] - x[0][0]*x[2][1]*x[1][2] - x[1][0]*x[0][1]*x[2][2] - x[1][0]*x[2][1]*x[0][2] - x[2][0]*x[0][1]*x[1][2] - x[2][0]*x[1][1]*x[0][2]); }

    template<unsigned int DIM, typename T>
    inline T det(Mat<DIM, T> const& x) {
        T result = T(0);
        for(size_t i = 0; i < DIM; ++i) {
            const T d = x[i][0] * det(x.minor(i, 0));
            result += ((i % 2) == 0) ? d : -d;
        }
        return result;
        VRS_STATIC_ASSERT(DIM >= 4);
    }

    template<typename T>
    inline Mat<1, T> inverse(Mat<1, T> const& x) {
        if(x[0][0] == T(0)) { return Mat<1, T>(T(0)); }
        return Mat<1, T>(T(1) / x[0][0]);
    }

    template<typename T>
    inline Mat<2, T> inverse(Mat<2, T> const& x) {
        const T d = det(x);
        if(d == T(0)) { return Mat<2, T>(T(0)); }
        const T invD = T(1) / d;
        return Mat<2, T>(invD * x[1][1], -invD * x[0][1], -invD * x[1][0], invD * x[0][0]);
    }

    template<unsigned int DIM, typename T>
    inline Mat<DIM, T> inverse(Mat<DIM, T> const& x) {
        Mat<DIM, T> inv;
        Mat<DIM, T> xx(x);

        // iterate over all columns
        for(size_t i = 0; i < DIM; ++i) {
            // look for column with largest pivot element
            T temp = T(0);
            size_t pivot = DIM;
            for(size_t j = i; j < DIM; ++j) {
                const T temp2 = fabs(xx[j][i]);
                if(temp2 > temp) {
                    temp = temp2;
                    pivot = j;
                }
            }
            if(pivot == DIM) {
                // No replacing column => no inverse
                return Mat<DIM, T>(T(0));
            }

            if(pivot != i) {
                // do the column swap
                std::swap(inv[i], inv[pivot]);
                std::swap(xx[i], xx[pivot]);
            }

            // change pivot element to 1.
            const T div = T(1) / xx[i][i];
            inv[i] *= div;
            xx[i] *= div;

            // do the gaussian elimination column wise
            for(size_t j = 0; j < DIM; ++j) {
                if(i == j) { continue; }
                const Vec<DIM, T> temp2 = xx[j][i] * xx[i];
                inv[j] -= temp2;
                xx[j] -= temp2;
            }
        }

        return inv;
        VRS_STATIC_ASSERT(DIM >= 3);
    }

    template<unsigned int DIM, typename T>
    inline Mat<DIM + 1, T> scaling(Vec<DIM, T> const& x) {
        return Mat<DIM + 1, T>(Vec<DIM + 1, T>(x, T(1)));
    }

    template<unsigned int DIM, typename T>
    inline Mat<DIM + 1, T> translation(Vec<DIM, T> const& x) {
        Mat<DIM + 1, T> result; // identity matrix
        result[DIM] = Vec<DIM + 1, T>(x, T(1)); // set "translation" column
        return result;
    }

} // namespace VRS

#include <vrs/vec_mat_ops.h>

#endif // VRS_MAT_H

---

Generated on Wed May 23 06:00:13 2012 by Doxygen 1.5.6

© 2001-2010 Hasso-Plattner-Institut | Impressum | Contact