old_deprecated_matrix.h

Go to the documentation of this file.

/****************************************************************************
* VCGLib                                                            o o     *
* Visual and Computer Graphics Library                            o     o   *
*                                                                _   O  _   *
* Copyright(C) 2004                                                \/)\/    *
* Visual Computing Lab                                            /\/|      *
* ISTI - Italian National Research Council                           |      *
*                                                                    \      *
* All rights reserved.                                                      *
*                                                                           *
* This program 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 (http://www.gnu.org/licenses/gpl.txt)          *
* for more details.                                                         *
*                                                                           *
****************************************************************************/
/***************************************************************************
$Log: not supported by cvs2svn $
Revision 1.9  2006/09/11 16:11:39  marfr960
Added const to declarations of the overloaded (operators *).
Otherwise the  * operator would always attempt to convert any type of data passed as an argument to Point3<TYPE>

Revision 1.8  2006/08/23 15:24:45  marfr960
Copy constructor : faster memcpy instead of slow 'for' cycle
empty constructor

Revision 1.7  2006/04/29 10:26:04  fiorin
Added some utility methods (swapping of columns and rows, matrix-vector multiplication)

Revision 1.6  2006/04/11 08:09:35  zifnab1974
changes necessary for gcc 3.4.5 on linux 64bit. Please take note of case-sensitivity of filenames

Revision 1.5  2005/12/12 11:25:00  ganovelli
added diagonal matrix, outer produce and namespace

***************************************************************************/

#ifndef MATRIX_VCGLIB
#define MATRIX_VCGLIB

#include <stdio.h>
#include <math.h>
#include <memory.h>
#include <assert.h>
#include <algorithm>
#include <vcg/space/point.h>

namespace vcg{
        namespace ndim{

   /* @{ */

        class MatrixDiagBase{public: 
        virtual const int & Dimension()const =0;
        virtual float operator[](const int & i)const = 0;
        };
        template<int N, class S> 
        class MatrixDiag: public Point<N,S>, public MatrixDiagBase{
        public:
                const int & Dimension() const {return N;}
                MatrixDiag(const Point<N,S>&p):Point<N,S>(p){}
        };

        template<class TYPE> 
        class Matrix
                {
                        
                public:
                        typedef TYPE ScalarType;

                        Matrix(unsigned int m, unsigned int n)
                        {
                                _rows = m;
                                _columns = n;
                                _data = new ScalarType[m*n];
                                memset(_data, 0, m*n*sizeof(ScalarType));
                        };

                        Matrix(unsigned int m, unsigned int n, TYPE *values)
                        {
                                _rows = m;
                                _columns = n;
                                unsigned int dim = m*n;
                                _data = new ScalarType[dim];
                                memcpy(_data, values, dim*sizeof(ScalarType));
                                //unsigned int i;
                                //for (i=0; i<_rows*_columns; i++)
                                //      _data[i] = values[i];
                        };

                        Matrix()
                        {
                                _rows = 0;
                                _columns = 0;
                                _data = NULL;
                        };

                        Matrix(const Matrix<TYPE> &m)
                        {
                                _rows = m._rows;
                                _columns = m._columns;
                                _data = new ScalarType[_rows*_columns];

                                unsigned int dim = _rows * _columns;
                                memcpy(_data, m._data, dim * sizeof(ScalarType));

//                              for (unsigned int i=0; i<_rows*_columns; i++)
//                                      _data[i] = m._data[i];
                        };

                        ~Matrix()
                        {
                                delete []_data;
                        };

                        inline unsigned int ColumnsNumber() const
                        {
                                return _columns;
                        };


                        inline unsigned int RowsNumber() const 
                        {
                                return _rows;
                        };

                        bool operator==(const Matrix<TYPE> &m) const
                        {
                                if (_rows==m._rows && _columns==m._columns)
                                {
                                        bool result = true;
                                        for (unsigned int i=0; i<_rows*_columns && result; i++)
                                                result = (_data[i]==m._data[i]);
                                        return result;
                                }
                                return false;
                        };

                        bool operator!=(const Matrix<TYPE> &m) const
                        {
                                if (_rows==m._rows && _columns==m._columns)
                                {
                                        bool result = false;
                                        for (unsigned int i=0; i<_rows*_columns && !result; i++)
                                                result = (_data[i]!=m._data[i]);
                                        return result;
                                }
                                return true;
                        };

                        inline TYPE ElementAt(unsigned int i, unsigned int j)
                        {
                                assert(i>=0 && i<_rows);
                                assert(j>=0 && j<_columns);
                                return _data[i*_columns+j];
                        };

                        TYPE Determinant() const
                        {
                                assert(_rows == _columns);
                                switch (_rows)
                                {
                                case 2:
                                        {
                                                return _data[0]*_data[3]-_data[1]*_data[2];
                                                break;
                                        };
                                case 3:
                                        {
                                                return  _data[0]*(_data[4]*_data[8]-_data[5]*_data[7]) - 
                                                                                _data[1]*(_data[3]*_data[8]-_data[5]*_data[6]) + 
                                                                                _data[2]*(_data[3]*_data[7]-_data[4]*_data[6]) ;
                                                break;
                                        };
                                default:
                                        {
                                                // da migliorare: si puo' cercare la riga/colonna con maggior numero di zeri
                                                ScalarType det = 0;
                                                for (unsigned int j=0; j<_columns; j++)
                                                        if (_data[j]!=0)
                                                                det += _data[j]*this->Cofactor(0, j);

                                                return det;
                                        }
                                };
                        };

                        TYPE Cofactor(unsigned int i, unsigned int j) const
                        {
                                assert(_rows == _columns);
                                assert(_rows>2);
                                TYPE *values = new TYPE[(_rows-1)*(_columns-1)];
                                unsigned int u, v, p, q, s, t;
                                for (u=0, p=0, s=0, t=0; u<_rows; u++, t+=_rows)
                                {
                                        if (i==u)
                                                continue;

                                        for (v=0, q=0; v<_columns; v++)
                                        {
                                                if (j==v)
                                                        continue;
                                                values[s+q] = _data[t+v];
                                                q++;
                                        }
                                        p++;
                                        s+=(_rows-1);
                                }
                                Matrix<TYPE> temp(_rows-1, _columns-1, values);
        return (pow(TYPE(-1.0), TYPE(i+j))*temp.Determinant());
                        };

                        inline TYPE* operator[](const unsigned int i)
                        {
        assert(i<_rows);
                                return _data + i*_columns;
                        };

                        inline const TYPE* operator[](const unsigned int i) const 
                        {
        assert(i<_rows);
                                return _data + i*_columns;
                        };

                        TYPE* GetColumn(const unsigned int j)
                        {
                                assert(j>=0 && j<_columns);
                                ScalarType *v = new ScalarType[_columns];
                                unsigned int i, p;
                                for (i=0, p=j; i<_rows; i++, p+=_columns)
                                        v[i] = _data[p];
                                return v;
                        };

                        TYPE* GetRow(const unsigned int i)
                        {
                                assert(i>=0 && i<_rows);
                                ScalarType *v = new ScalarType[_rows];
                                unsigned int j, p;
                                for (j=0, p=i*_columns; j<_columns; j++, p++)
                                        v[j] = _data[p];
                                return v;
                        };

                        void SwapColumns(const unsigned int i, const unsigned int j)
                        {
                                assert(0<=i && i<_columns);
                                assert(0<=j && j<_columns);
                                if (i==j)
                                        return;

                                unsigned int r, e0, e1;
                                for (r=0, e0=i, e1=j; r<_rows; r++, e0+=_columns, e1+=_columns)
                                        std::swap(_data[e0], _data[e1]);
                        };

                        void SwapRows(const unsigned int i, const unsigned int j)
                        {
                                assert(0<=i && i<_rows);
                                assert(0<=j && j<_rows);
                                if (i==j)
                                        return;

                                unsigned int r, e0, e1;
                                for (r=0, e0=i*_columns, e1=j*_columns; r<_columns; r++, e0++, e1++)
                                        std::swap(_data[e0], _data[e1]);
                        };

                        Matrix<TYPE>& operator=(const Matrix<TYPE> &m)
                        {
                                if (this != &m)
                                {
                                        assert(_rows == m._rows);
                                        assert(_columns == m._columns);
                                        for (unsigned int i=0; i<_rows*_columns; i++)
                                                _data[i] = m._data[i];
                                }
                                return *this;
                        };

                        Matrix<TYPE>& operator+=(const Matrix<TYPE> &m)
                        {
                                assert(_rows == m._rows);
                                assert(_columns == m._columns);
                                for (unsigned int i=0; i<_rows*_columns; i++)
                                        _data[i] += m._data[i];
                                return *this;
                        };

                        Matrix<TYPE>& operator-=(const Matrix<TYPE> &m)
                        {
                                assert(_rows == m._rows);
                                assert(_columns == m._columns);
                                for (unsigned int i=0; i<_rows*_columns; i++)
                                        _data[i] -= m._data[i];
                                return *this;
                        };

                        Matrix<TYPE>& operator+=(const TYPE k)
                        {
                                for (unsigned int i=0; i<_rows*_columns; i++)
                                        _data[i] += k;
                                return *this;
                        };

                        Matrix<TYPE>& operator-=(const TYPE k)
                        {
                                for (unsigned int i=0; i<_rows*_columns; i++)
                                        _data[i] -= k;
                                return *this;
                        };

                        Matrix<TYPE>& operator*=(const TYPE k)
                        {
                                for (unsigned int i=0; i<_rows*_columns; i++)
                                        _data[i] *= k;
                                return *this;
                        };

                        Matrix<TYPE>& operator/=(const TYPE k)
                        {
                                assert(k!=0);
                                for (unsigned int i=0; i<_rows*_columns; i++)
                                        _data[i] /= k;
                                return *this;
                        };

                        Matrix<TYPE> operator*(const Matrix<TYPE> &m) const
                        {
                                assert(_columns == m._rows);
                                Matrix<TYPE> result(_rows, m._columns);
                                unsigned int i, j, k, p, q, r;
                                for (i=0, p=0, r=0; i<result._rows; i++, p+=_columns, r+=result._columns)
                                        for (j=0; j<result._columns; j++)
                                        {
                                                ScalarType temp = 0;
                                                for (k=0, q=0; k<_columns; k++, q+=m._columns)
                                                        temp+=(_data[p+k]*m._data[q+j]);
                                                result._data[r+j] = temp;
                                        }
                                
                                return result;
                        };

                        ScalarType* operator*(const ScalarType v[]) const
                        {
                                ScalarType *result = new ScalarType[_rows];
                                memset(result, 0, _rows*sizeof(ScalarType));
                                unsigned int r, c, i;
                                for (r=0, i=0; r<_rows; r++)
                                        for (c=0; c<_columns; c++, i++)
                                                result[r] += _data[i]*v[c];

                                return result;
                        };

                        template <int N,int M>
                        void DotProduct(Point<N,TYPE> &m,Point<M,TYPE> &result)
                        {
                                unsigned int i, j,  p,  r;
                                for (i=0, p=0, r=0; i<M; i++)
                                { result[i]=0;
                                        for (j=0; j<N; j++)
                                                result[i]+=(*this)[i][j]*m[j];
                                }
                        };

                        Matrix<TYPE> operator*(const MatrixDiagBase &m) const
                        {
                                assert(_columns == _rows);
                                assert(_columns == m.Dimension());
                                int i,j;
                                Matrix<TYPE> result(_rows, _columns);

                                for (i=0; i<result._rows; i++)
                                        for (j=0; j<result._columns; j++)
                                                result[i][j]*= m[j];
                                
                                return result;
                        };
                        template <int N, int M>
                        void OuterProduct(const Point<N,TYPE> a, const Point< M,TYPE> b)
                        {
                                assert(N == _rows);
                                assert(M == _columns);
                                Matrix<TYPE> result(_rows,_columns);
                                unsigned int i, j;

                                for (i=0; i<result._rows; i++)
                                        for (j=0; j<result._columns; j++)
                                                (*this)[i][j] = a[i] * b[j];
                        };


                        Point3<TYPE> operator*(Point3<TYPE> &p) const
                        {
                                assert(_columns==3 && _rows==3);
                                vcg::Point3<TYPE> result;
                                result[0] = _data[0]*p[0]+_data[1]*p[1]+_data[2]*p[2];
                                result[1] = _data[3]*p[0]+_data[4]*p[1]+_data[5]*p[2];
                                result[2] = _data[6]*p[0]+_data[7]*p[1]+_data[8]*p[2];
                                return result;
                        };


                        Matrix<TYPE> operator+(const TYPE k)
                        {
                                Matrix<TYPE> result(_rows, _columns);
                                for (unsigned int i=0; i<_rows*_columns; i++)
                                        result._data[i] =  _data[i]+k;
                                return result;
                        };

                        Matrix<TYPE> operator-(const TYPE k)
                        {
                                Matrix<TYPE> result(_rows, _columns);
                                for (unsigned int i=0; i<_rows*_columns; i++)
                                        result._data[i] =  _data[i]-k;
                                return result;
                        };

                        Matrix<TYPE> operator-() const
                        {
                                Matrix<TYPE> result(_rows, _columns, _data);
                                for (unsigned int i=0; i<_columns*_rows; i++)
                                        result._data[i] = -1*_data[i];
                                return result;
                        };

                        Matrix<TYPE> operator*(const TYPE k) const
                        {
                                Matrix<TYPE> result(_rows, _columns);
                                for (unsigned int i=0; i<_rows*_columns; i++)
                                        result._data[i] =  _data[i]*k;
                                return result;
                        };

                        Matrix<TYPE> operator/(const TYPE k)
                        {
                                Matrix<TYPE> result(_rows, _columns);
                                for (unsigned int i=0; i<_rows*_columns; i++)
                                        result._data[i] =  _data[i]/k;
                                return result;
                        };


                        void SetZero()
                        {
                                for (unsigned int i=0; i<_rows*_columns; i++)
                                        _data[i] = ScalarType(0.0);
                        };

                        void SetIdentity()
                        {
                                assert(_rows==_columns);
                                for (unsigned int i=0; i<_rows; i++)
                                        for (unsigned int j=0; j<_columns; j++)
                                                _data[i] = (i==j) ? ScalarType(1.0) : ScalarType(0.0f);
                        };

                        void SetColumn(const unsigned int j, TYPE* v)
                        {
                                assert(j>=0 && j<_columns);
                                unsigned int i, p;
                                for (i=0, p=j; i<_rows; i++, p+=_columns)
                                        _data[p] = v[i];
                        };

                        void SetRow(const unsigned int i, TYPE* v)
                        {
                                assert(i>=0 && i<_rows);
                                unsigned int j, p;
                                for (j=0, p=i*_rows; j<_columns; j++, p++)
                                        _data[p] = v[j];
                        };

                        void SetDiagonal(TYPE *v)
                        {
                                assert(_rows == _columns);
                                for (unsigned int i=0, p=0; i<_rows; i++, p+=_rows)
                                        _data[p+i] = v[i];
                        };

                        void Resize(const unsigned int m, const unsigned int n)
                        {
                                assert(m>=2);
                                assert(n>=2);
                                _rows = m;
                                _columns = n;
                                delete []_data;
                                _data = new ScalarType[m*n];
                                for (unsigned int i=0; i<m*n; i++)
                                        _data[i] = 0;
                        };


                        void Transpose()
                        {
                                ScalarType *temp = new ScalarType[_rows*_columns];
                                unsigned int i, j, p, q;
                                for (i=0, p=0; i<_rows; i++, p+=_columns)
                                        for (j=0, q=0; j<_columns; j++, q+=_rows)
                                                temp[q+i] = _data[p+j];
                                
                                std::swap(_columns, _rows);
                                std::swap(_data, temp);
                                delete []temp;
                        };

                        // for the transistion to eigen
                        Matrix transpose()
                        {
                                Matrix res = *this;
                                res.Transpose();
                                return res;
                        }
                        void transposeInPlace() { Transpose(); }
                        // for the transistion to eigen

                        void Dump()
                        {
                                unsigned int i, j, p;
                                for (i=0, p=0; i<_rows; i++, p+=_columns)
                                {
                                        printf("[\t");
                                        for (j=0; j<_columns; j++)
                                                printf("%f\t", _data[p+j]);
                                        
                                        printf("]\n");
                                }
                                printf("\n");
                        };

                protected:
                        unsigned int _rows;

                        unsigned int _columns;

                        ScalarType *_data;
                };

                typedef vcg::ndim::Matrix<double> MatrixMNd;
                typedef vcg::ndim::Matrix<float>  MatrixMNf;

//      template <class MatrixType>
//      void Invert(MatrixType & m){
//              typedef typename MatrixType::ScalarType X;
//              X  *diag;
//              diag = new  X [m.ColumnsNumber()];

//              MatrixType res(m.RowsNumber(),m.ColumnsNumber());
//              vcg::SingularValueDecomposition<MatrixType > (m,&diag[0],res,LeaveUnsorted,50 );
//              m.Transpose();
//              // prodotto per la diagonale
//              unsigned  int i,j;
//              for (i=0; i<m.RowsNumber(); i++)
//                                      for (j=0; j<m.ColumnsNumber(); j++)
//                                              res[i][j]/= diag[j];
//              m = res *m;
//              }

        }
}; // end of namespace

#endif