00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044 #ifndef MATRIX_VCGLIB
00045 #define MATRIX_VCGLIB
00046
00047 #include <stdio.h>
00048 #include <math.h>
00049 #include <memory.h>
00050 #include <assert.h>
00051 #include <algorithm>
00052 #include <vcg/space/point.h>
00053 #include <vcg/math/lin_algebra.h>
00054
00055 namespace vcg{
00056 namespace ndim{
00057
00059
00060
00065 class MatrixDiagBase{public:
00066 virtual const int & Dimension()const =0;
00067 virtual const float operator[](const int & i)const = 0;
00068 };
00069 template<int N, class S>
00070 class MatrixDiag: public Point<N,S>, public MatrixDiagBase{
00071 public:
00072 const int & Dimension() const {return N;}
00073 MatrixDiag(const Point<N,S>&p):Point<N,S>(p){}
00074 };
00075
00080 template<class TYPE>
00081 class Matrix
00082 {
00083
00084 public:
00085 typedef TYPE ScalarType;
00086
00093 Matrix(unsigned int m, unsigned int n)
00094 {
00095 _rows = m;
00096 _columns = n;
00097 _data = new ScalarType[m*n];
00098 memset(_data, 0, m*n*sizeof(ScalarType));
00099 };
00100
00108 Matrix(unsigned int m, unsigned int n, TYPE *values)
00109 {
00110 _rows = m;
00111 _columns = n;
00112 unsigned int dim = m*n;
00113 _data = new ScalarType[dim];
00114 memcpy(_data, values, dim*sizeof(ScalarType));
00115
00116
00117
00118 };
00119
00124 Matrix()
00125 {
00126 _rows = 0;
00127 _columns = 0;
00128 _data = NULL;
00129 };
00130
00136 Matrix(const Matrix<TYPE> &m)
00137 {
00138 _rows = m._rows;
00139 _columns = m._columns;
00140 _data = new ScalarType[_rows*_columns];
00141
00142 unsigned int dim = _rows * _columns;
00143 memcpy(_data, m._data, dim * sizeof(ScalarType));
00144
00145
00146
00147 };
00148
00152 ~Matrix()
00153 {
00154 delete []_data;
00155 };
00156
00160 inline unsigned int ColumnsNumber() const
00161 {
00162 return _columns;
00163 };
00164
00165
00169 inline unsigned int RowsNumber() const
00170 {
00171 return _rows;
00172 };
00173
00179 bool operator==(const Matrix<TYPE> &m) const
00180 {
00181 if (_rows==m._rows && _columns==m._columns)
00182 {
00183 bool result = true;
00184 for (unsigned int i=0; i<_rows*_columns && result; i++)
00185 result = (_data[i]==m._data[i]);
00186 return result;
00187 }
00188 return false;
00189 };
00190
00196 bool operator!=(const Matrix<TYPE> &m) const
00197 {
00198 if (_rows==m._rows && _columns==m._columns)
00199 {
00200 bool result = false;
00201 for (unsigned int i=0; i<_rows*_columns && !result; i++)
00202 result = (_data[i]!=m._data[i]);
00203 return result;
00204 }
00205 return true;
00206 };
00207
00214 inline TYPE ElementAt(unsigned int i, unsigned int j)
00215 {
00216 assert(i>=0 && i<_rows);
00217 assert(j>=0 && j<_columns);
00218 return _data[i*_columns+j];
00219 };
00220
00225 TYPE Determinant() const
00226 {
00227 assert(_rows == _columns);
00228 switch (_rows)
00229 {
00230 case 2:
00231 {
00232 return _data[0]*_data[3]-_data[1]*_data[2];
00233 break;
00234 };
00235 case 3:
00236 {
00237 return _data[0]*(_data[4]*_data[8]-_data[5]*_data[7]) -
00238 _data[1]*(_data[3]*_data[8]-_data[5]*_data[6]) +
00239 _data[2]*(_data[3]*_data[7]-_data[4]*_data[6]) ;
00240 break;
00241 };
00242 default:
00243 {
00244
00245 ScalarType det = 0;
00246 for (unsigned int j=0; j<_columns; j++)
00247 if (_data[j]!=0)
00248 det += _data[j]*this->Cofactor(0, j);
00249
00250 return det;
00251 }
00252 };
00253 };
00254
00259 TYPE Cofactor(unsigned int i, unsigned int j) const
00260 {
00261 assert(_rows == _columns);
00262 assert(_rows>2);
00263 TYPE *values = new TYPE[(_rows-1)*(_columns-1)];
00264 unsigned int u, v, p, q, s, t;
00265 for (u=0, p=0, s=0, t=0; u<_rows; u++, t+=_rows)
00266 {
00267 if (i==u)
00268 continue;
00269
00270 for (v=0, q=0; v<_columns; v++)
00271 {
00272 if (j==v)
00273 continue;
00274 values[s+q] = _data[t+v];
00275 q++;
00276 }
00277 p++;
00278 s+=(_rows-1);
00279 }
00280 Matrix<TYPE> temp(_rows-1, _columns-1, values);
00281 return (pow(-1, i+j)*temp.Determinant());
00282 };
00283
00289 inline TYPE* operator[](const unsigned int i)
00290 {
00291 assert(i<_rows);
00292 return _data + i*_columns;
00293 };
00294
00300 inline const TYPE* operator[](const unsigned int i) const
00301 {
00302 assert(i>=0 && i<_rows);
00303 return _data + i*_columns;
00304 };
00305
00311 TYPE* GetColumn(const unsigned int j)
00312 {
00313 assert(j>=0 && j<_columns);
00314 ScalarType *v = new ScalarType[_columns];
00315 unsigned int i, p;
00316 for (i=0, p=j; i<_rows; i++, p+=_columns)
00317 v[i] = _data[p];
00318 return v;
00319 };
00320
00326 TYPE* GetRow(const unsigned int i)
00327 {
00328 assert(i>=0 && i<_rows);
00329 ScalarType *v = new ScalarType[_rows];
00330 unsigned int j, p;
00331 for (j=0, p=i*_columns; j<_columns; j++, p++)
00332 v[j] = _data[p];
00333 return v;
00334 };
00335
00341 void SwapColumns(const unsigned int i, const unsigned int j)
00342 {
00343 assert(0<=i && i<_columns);
00344 assert(0<=j && j<_columns);
00345 if (i==j)
00346 return;
00347
00348 unsigned int r, e0, e1;
00349 for (r=0, e0=i, e1=j; r<_rows; r++, e0+=_columns, e1+=_columns)
00350 std::swap(_data[e0], _data[e1]);
00351 };
00352
00358 void SwapRows(const unsigned int i, const unsigned int j)
00359 {
00360 assert(0<=i && i<_rows);
00361 assert(0<=j && j<_rows);
00362 if (i==j)
00363 return;
00364
00365 unsigned int r, e0, e1;
00366 for (r=0, e0=i*_columns, e1=j*_columns; r<_columns; r++, e0++, e1++)
00367 std::swap(_data[e0], _data[e1]);
00368 };
00369
00374 Matrix<TYPE>& operator=(const Matrix<TYPE> &m)
00375 {
00376 if (this != &m)
00377 {
00378 assert(_rows == m._rows);
00379 assert(_columns == m._columns);
00380 for (unsigned int i=0; i<_rows*_columns; i++)
00381 _data[i] = m._data[i];
00382 }
00383 return *this;
00384 };
00385
00391 Matrix<TYPE>& operator+=(const Matrix<TYPE> &m)
00392 {
00393 assert(_rows == m._rows);
00394 assert(_columns == m._columns);
00395 for (unsigned int i=0; i<_rows*_columns; i++)
00396 _data[i] += m._data[i];
00397 return *this;
00398 };
00399
00405 Matrix<TYPE>& operator-=(const Matrix<TYPE> &m)
00406 {
00407 assert(_rows == m._rows);
00408 assert(_columns == m._columns);
00409 for (unsigned int i=0; i<_rows*_columns; i++)
00410 _data[i] -= m._data[i];
00411 return *this;
00412 };
00413
00419 Matrix<TYPE>& operator+=(const TYPE k)
00420 {
00421 for (unsigned int i=0; i<_rows*_columns; i++)
00422 _data[i] += k;
00423 return *this;
00424 };
00425
00431 Matrix<TYPE>& operator-=(const TYPE k)
00432 {
00433 for (unsigned int i=0; i<_rows*_columns; i++)
00434 _data[i] -= k;
00435 return *this;
00436 };
00437
00443 Matrix<TYPE>& operator*=(const TYPE k)
00444 {
00445 for (unsigned int i=0; i<_rows*_columns; i++)
00446 _data[i] *= k;
00447 return *this;
00448 };
00449
00455 Matrix<TYPE>& operator/=(const TYPE k)
00456 {
00457 assert(k!=0);
00458 for (unsigned int i=0; i<_rows*_columns; i++)
00459 _data[i] /= k;
00460 return *this;
00461 };
00462
00468 Matrix<TYPE> operator*(const Matrix<TYPE> &m) const
00469 {
00470 assert(_columns == m._rows);
00471 Matrix<TYPE> result(_rows, m._columns);
00472 unsigned int i, j, k, p, q, r;
00473 for (i=0, p=0, r=0; i<result._rows; i++, p+=_columns, r+=result._columns)
00474 for (j=0; j<result._columns; j++)
00475 {
00476 ScalarType temp = 0;
00477 for (k=0, q=0; k<_columns; k++, q+=m._columns)
00478 temp+=(_data[p+k]*m._data[q+j]);
00479 result._data[r+j] = temp;
00480 }
00481
00482 return result;
00483 };
00484
00490 ScalarType* operator*(const ScalarType v[]) const
00491 {
00492 ScalarType *result = new ScalarType[_rows];
00493 memset(result, 0, _rows*sizeof(ScalarType));
00494 unsigned int r, c, i;
00495 for (r=0, i=0; r<_rows; r++)
00496 for (c=0; c<_columns; c++, i++)
00497 result[r] += _data[i]*v[c];
00498
00499 return result;
00500 };
00501
00507 template <int N,int M>
00508 void DotProduct(Point<N,TYPE> &m,Point<M,TYPE> &result)
00509 {
00510 unsigned int i, j, p, r;
00511 for (i=0, p=0, r=0; i<M; i++)
00512 { result[i]=0;
00513 for (j=0; j<N; j++)
00514 result[i]+=(*this)[i][j]*m[j];
00515 }
00516 };
00517
00521 Matrix<TYPE> operator*(const MatrixDiagBase &m) const
00522 {
00523 assert(_columns == _rows);
00524 assert(_columns == m.Dimension());
00525 int i,j;
00526 Matrix<TYPE> result(_rows, _columns);
00527
00528 for (i=0; i<result._rows; i++)
00529 for (j=0; j<result._columns; j++)
00530 result[i][j]*= m[j];
00531
00532 return result;
00533 };
00537 template <int N, int M>
00538 void OuterProduct(const Point<N,TYPE> a, const Point< M,TYPE> b)
00539 {
00540 assert(N == _rows);
00541 assert(M == _columns);
00542 Matrix<TYPE> result(_rows,_columns);
00543 unsigned int i, j;
00544
00545 for (i=0; i<result._rows; i++)
00546 for (j=0; j<result._columns; j++)
00547 (*this)[i][j] = a[i] * b[j];
00548 };
00549
00550
00557 Point3<TYPE> operator*(Point3<TYPE> &p) const
00558 {
00559 assert(_columns==3 && _rows==3);
00560 vcg::Point3<TYPE> result;
00561 result[0] = _data[0]*p[0]+_data[1]*p[1]+_data[2]*p[2];
00562 result[1] = _data[3]*p[0]+_data[4]*p[1]+_data[5]*p[2];
00563 result[2] = _data[6]*p[0]+_data[7]*p[1]+_data[8]*p[2];
00564 return result;
00565 };
00566
00567
00573 Matrix<TYPE> operator+(const TYPE k)
00574 {
00575 Matrix<TYPE> result(_rows, _columns);
00576 for (unsigned int i=0; i<_rows*_columns; i++)
00577 result._data[i] = _data[i]+k;
00578 return result;
00579 };
00580
00586 Matrix<TYPE> operator-(const TYPE k)
00587 {
00588 Matrix<TYPE> result(_rows, _columns);
00589 for (unsigned int i=0; i<_rows*_columns; i++)
00590 result._data[i] = _data[i]-k;
00591 return result;
00592 };
00593
00598 Matrix<TYPE> operator-() const
00599 {
00600 Matrix<TYPE> result(_rows, _columns, _data);
00601 for (unsigned int i=0; i<_columns*_rows; i++)
00602 result._data[i] = -1*_data[i];
00603 return result;
00604 };
00605
00611 Matrix<TYPE> operator*(const TYPE k) const
00612 {
00613 Matrix<TYPE> result(_rows, _columns);
00614 for (unsigned int i=0; i<_rows*_columns; i++)
00615 result._data[i] = _data[i]*k;
00616 return result;
00617 };
00618
00624 Matrix<TYPE> operator/(const TYPE k)
00625 {
00626 Matrix<TYPE> result(_rows, _columns);
00627 for (unsigned int i=0; i<_rows*_columns; i++)
00628 result._data[i] = _data[i]/k;
00629 return result;
00630 };
00631
00632
00636 void SetZero()
00637 {
00638 for (unsigned int i=0; i<_rows*_columns; i++)
00639 _data[i] = ScalarType(0.0);
00640 };
00641
00645 void SetIdentity()
00646 {
00647 assert(_rows==_columns);
00648 for (unsigned int i=0; i<_rows; i++)
00649 for (unsigned int j=0; j<_columns; j++)
00650 _data[i] = (i==j) ? ScalarType(1.0) : ScalarType(0.0f);
00651 };
00652
00658 void SetColumn(const unsigned int j, TYPE* v)
00659 {
00660 assert(j>=0 && j<_columns);
00661 unsigned int i, p;
00662 for (i=0, p=j; i<_rows; i++, p+=_columns)
00663 _data[p] = v[i];
00664 };
00665
00671 void SetRow(const unsigned int i, TYPE* v)
00672 {
00673 assert(i>=0 && i<_rows);
00674 unsigned int j, p;
00675 for (j=0, p=i*_rows; j<_columns; j++, p++)
00676 _data[p] = v[j];
00677 };
00678
00683 void SetDiagonal(TYPE *v)
00684 {
00685 assert(_rows == _columns);
00686 for (unsigned int i=0, p=0; i<_rows; i++, p+=_rows)
00687 _data[p+i] = v[i];
00688 };
00689
00695 void Resize(const unsigned int m, const unsigned int n)
00696 {
00697 assert(m>=2);
00698 assert(n>=2);
00699 _rows = m;
00700 _columns = n;
00701 delete []_data;
00702 _data = new ScalarType[m*n];
00703 for (unsigned int i=0; i<m*n; i++)
00704 _data[i] = 0;
00705 };
00706
00707
00711 void Transpose()
00712 {
00713 ScalarType *temp = new ScalarType[_rows*_columns];
00714 unsigned int i, j, p, q;
00715 for (i=0, p=0; i<_rows; i++, p+=_columns)
00716 for (j=0, q=0; j<_columns; j++, q+=_rows)
00717 temp[q+i] = _data[p+j];
00718
00719 std::swap(_columns, _rows);
00720 std::swap(_data, temp);
00721 delete []temp;
00722 };
00723
00724
00725 Matrix transpose()
00726 {
00727 Matrix res = *this;
00728 res.Transpose();
00729 return res;
00730 }
00731 void transposeInPlace() { Transpose(); }
00732
00733
00737 void Dump()
00738 {
00739 unsigned int i, j, p;
00740 for (i=0, p=0; i<_rows; i++, p+=_columns)
00741 {
00742 printf("[\t");
00743 for (j=0; j<_columns; j++)
00744 printf("%f\t", _data[p+j]);
00745
00746 printf("]\n");
00747 }
00748 printf("\n");
00749 };
00750
00751 protected:
00753 unsigned int _rows;
00754
00756 unsigned int _columns;
00757
00759 ScalarType *_data;
00760 };
00761
00762 typedef vcg::ndim::Matrix<double> MatrixMNd;
00763 typedef vcg::ndim::Matrix<float> MatrixMNf;
00764
00767 template <class MatrixType>
00768 void Invert(MatrixType & m){
00769 typedef typename MatrixType::ScalarType X;
00770 X *diag;
00771 diag = new X [m.ColumnsNumber()];
00772
00773 MatrixType res(m.RowsNumber(),m.ColumnsNumber());
00774 vcg::SingularValueDecomposition<MatrixType > (m,&diag[0],res,LeaveUnsorted,50 );
00775 m.Transpose();
00776
00777 unsigned int i,j;
00778 for (i=0; i<m.RowsNumber(); i++)
00779 for (j=0; j<m.ColumnsNumber(); j++)
00780 res[i][j]/= diag[j];
00781 m = res *m;
00782 }
00783
00784 }
00785 };
00786
00787 #endif
