source: foam/trunk/vision/src/Matrix.h @ 90

Revision 90, 10.1 KB checked in by dave, 10 years ago (diff)

PCA updated - projection and resynthesis seem to work

RevLine 
[85]1// Copyright (C) 2009 foam
2//
3// This program is free software; you can redistribute it and/or modify
4// it under the terms of the GNU General Public License as published by
5// the Free Software Foundation; either version 2 of the License, or
6// (at your option) any later version.
7//
8// This program is distributed in the hope that it will be useful,
9// but WITHOUT ANY WARRANTY; without even the implied warranty of
10// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
11// GNU General Public License for more details.
12//
13// You should have received a copy of the GNU General Public License
14// along with this program; if not, write to the Free Software
15// Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
16
17#include <assert.h>
18#include <iostream>
[86]19#include "Vector.h"
[85]20
21#ifndef FOAM_MATRIX
22#define FOAM_MATRIX
23
[86]24template<class T>
[85]25class Matrix
26{
27public:
[86]28        Matrix();
[85]29        Matrix(unsigned int r, unsigned int c);
30        ~Matrix();
31        Matrix(const Matrix &other);
[86]32        Matrix(unsigned int r, unsigned int c, float *data);
33
[85]34        // Row proxy classes to allow matrix[r][c] notation
35        class Row
36        {
37        public:
38                Row(Matrix *owner, unsigned int r)
39                {
40                        m_Data=&owner->GetRawData()[r*owner->GetCols()];
41                        m_Cols=owner->GetCols();
42                }
43               
44                T &operator[](unsigned int c)
45                {
46                        assert(c<m_Cols);
47                        return m_Data[c];
48                }
49               
50        private:
51                T *m_Data;
52                unsigned int m_Cols;
53        };
54
55        class ConstRow
56        {
57        public:
58                ConstRow(const Matrix *owner, unsigned int r)
59                {
60                        m_Data=&owner->GetRawDataConst()[r*owner->GetCols()];
61                        m_Cols=owner->GetCols();
62                }
63               
64                const T &operator[](unsigned int c) const
65                {
66                        assert(c<m_Cols);
67                        return m_Data[c];
68                }
69               
70        private:
71                const T *m_Data;
72                unsigned int m_Cols;
73        };
74
75       
76        Row operator[](unsigned int r)
77        {
78                assert(r<m_Rows);
79                return Row(this,r);
80        }
81
82        ConstRow operator[](unsigned int r) const
83        {
84                assert(r<m_Rows);
85                return ConstRow(this,r);
86        }
87 
88        unsigned int GetRows() const { return m_Rows; }
89        unsigned int GetCols() const { return m_Cols; }
90        T *GetRawData() { return m_Data; }
91        const T *GetRawDataConst() const { return m_Data; }
[86]92        Vector<T> GetRowVector(unsigned int r) const;
93        Vector<T> GetColVector(unsigned int c) const;
[89]94        void SetRowVector(unsigned int r, const Vector<T> &row);
95        void SetColVector(unsigned int c, const Vector<T> &col);
[85]96         
97        void Print() const;
98        void SetAll(T s);
[86]99        void Zero() { SetAll(0); }
100        bool IsInf();
[90]101        Matrix Transposed() const;
[85]102
103        Matrix &operator=(const Matrix &other);
104        Matrix operator+(const Matrix &other) const;
105        Matrix operator-(const Matrix &other) const;
106        Matrix operator*(const Matrix &other) const;
[90]107        Vector<T> operator*(const Vector<T> &other) const;
108        Vector<T> VecMulTransposed(const Vector<T> &other) const;
[85]109        Matrix &operator+=(const Matrix &other);
110        Matrix &operator-=(const Matrix &other);
111        Matrix &operator*=(const Matrix &other);
[90]112
[89]113        void SortRows(Vector<T> &v);
114        void SortCols(Vector<T> &v);
[90]115       
116        Matrix CropRows(unsigned int s, unsigned int e);
117        Matrix CropCols(unsigned int s, unsigned int e);
[89]118
[90]119        void Save(FILE *f);
120        void Load(FILE *f);
121
[85]122        static void RunTests();
123       
124private:
125
126        unsigned int m_Rows;
127        unsigned int m_Cols;
128       
129        T *m_Data;
130       
131};
132
133template<class T>
134Matrix<T>::Matrix(unsigned int r, unsigned int c) :
135m_Rows(r),
136m_Cols(c)
137{
138        m_Data=new T[r*c];
139}
140
141template<class T>
[86]142Matrix<T>::Matrix() :
143m_Rows(0),
144m_Cols(0),
145m_Data(NULL)
146{
147}
148
149template<class T>
150Matrix<T>::Matrix(unsigned int r, unsigned int c, float *data) :
151m_Rows(r),
152m_Cols(c),
153m_Data(data)
154{
155}
156
157template<class T>
[85]158Matrix<T>::~Matrix()
159{
160        delete[] m_Data;
161}
162
163template<class T>
164Matrix<T>::Matrix(const Matrix &other)
165{
166        m_Rows = other.m_Rows;
167        m_Cols = other.m_Cols;
168        m_Data=new T[m_Rows*m_Cols];
169        memcpy(m_Data,other.m_Data,m_Rows*m_Cols*sizeof(T));
170}
171
172template<class T>
173Matrix<T> &Matrix<T>::operator=(const Matrix &other)
174{
175        if (m_Data!=NULL)
176        {
177                delete[] m_Data;
178        }
179       
180        m_Rows = other.m_Rows;
181        m_Cols = other.m_Cols;
182        m_Data=new T[m_Rows*m_Cols];
183        memcpy(m_Data,other.m_Data,m_Rows*m_Cols*sizeof(T));
184       
185        return *this;
186}
187
188template<class T>
189void Matrix<T>::Print() const
190{
191        for (unsigned int i=0; i<m_Rows; i++)
192        {
193                for (unsigned int j=0; j<m_Cols; j++)
194                {
195                        std::cerr<<(*this)[i][j]<<" ";
196                }
197                std::cerr<<std::endl;
198        }
199}
200
201template<class T>
202void Matrix<T>::SetAll(T s)
203{
204        for (unsigned int i=0; i<m_Rows; i++)
205        {
206                for (unsigned int j=0; j<m_Cols; j++)
207                {
208                        (*this)[i][j]=s;
209                }
210        }
211}
212
213template<class T>
[86]214bool Matrix<T>::IsInf()
215{
216        for (unsigned int i=0; i<m_Rows; i++)
217        {
218                for (unsigned int j=0; j<m_Cols; j++)
219                {
220                        if (isinf((*this)[i][j])) return true;
221                        if (isnan((*this)[i][j])) return true;
222                }
223        }
224        return false;
225}
226
227template<class T>
[90]228Matrix<T> Matrix<T>::Transposed() const
[85]229{
230        Matrix<T> copy(*this);
231        for (unsigned int i=0; i<m_Rows; i++)
232        {
233                for (unsigned int j=0; j<m_Cols; j++)
234                {
235                        copy[i][j]=(*this)[j][i];
236                }
237        }
238        return copy;
239}
240
241
242template<class T>
243Matrix<T> Matrix<T>::operator+(const Matrix &other) const
244{
245        assert(m_Rows=other.m_Rows);
246        assert(m_Cols=other.m_Cols);
247       
248        Matrix<T> ret(m_Rows,m_Cols);
249        for (unsigned int i=0; i<m_Rows; i++)
250        {
251                for (unsigned int j=0; j<m_Cols; j++)
252                {
253                        ret[i][j]=(*this)[i][j]+other[i][j];
254                }
255        }
256        return ret;
257}
258
259template<class T>
260Matrix<T> Matrix<T>::operator-(const Matrix &other) const
261{
262        assert(m_Rows=other.m_Rows);
263        assert(m_Cols=other.m_Cols);
264       
265        Matrix<T> ret(m_Rows,m_Cols);
266        for (unsigned int i=0; i<m_Rows; i++)
267        {
268                for (unsigned int j=0; j<m_Cols; j++)
269                {
270                        ret[i][j]=(*this)[i][j]-other[i][j];
271                }
272        }
273        return ret;
274}
275
276template<class T>
277Matrix<T> Matrix<T>::operator*(const Matrix &other) const
278{
279        assert(m_Cols==other.m_Rows);
280       
281        Matrix<T> ret(m_Rows,other.m_Cols);
282       
283        for (unsigned int i=0; i<m_Rows; i++)
284        {
285                for (unsigned int j=0; j<other.m_Cols; j++)
286                {
287                        ret[i][j]=0;
288                        for (unsigned int k=0; k<m_Cols; k++)
289                        {
290                                ret[i][j]+=(*this)[i][k]*other[k][j];
291                        }
292                }
293        }
294        return ret;
295}
296
297template<class T>
[90]298Vector<T> Matrix<T>::operator*(const Vector<T> &other) const
299{
300        assert(m_Cols==other.Size());
301       
302        Vector<T> ret(m_Rows);
303       
304        for (unsigned int i=0; i<m_Rows; i++)
305        {
306                for (unsigned int j=0; j<other.Size(); j++)
307                {
308                        ret[i]=0;
309                        for (unsigned int k=0; k<m_Cols; k++)
310                        {
311                                ret[i]+=(*this)[i][k]*other[k];
312                        }
313                }
314        }
315        return ret;
316}
317
318template<class T>
319Vector<T> Matrix<T>::VecMulTransposed(const Vector<T> &other) const
320{
321        assert(m_Rows==other.Size());
322       
323        Vector<T> ret(m_Cols);
324       
325        for (unsigned int i=0; i<m_Cols; i++)
326        {
327                for (unsigned int j=0; j<other.Size(); j++)
328                {
329                        ret[i]=0;
330                        for (unsigned int k=0; k<m_Rows; k++)
331                        {
332                                ret[i]+=(*this)[k][i]*other[k];
333                        }
334                }
335        }
336        return ret;
337}
338
339template<class T>
[85]340Matrix<T> &Matrix<T>::operator+=(const Matrix &other)
341{
342        (*this)=(*this)+other;
343        return *this;
344}
345
346template<class T>
347Matrix<T> &Matrix<T>::operator-=(const Matrix &other)
348{
349        (*this)=(*this)-other;
350        return *this;
351}
352
353template<class T>
354Matrix<T> &Matrix<T>::operator*=(const Matrix &other)
355{
356        (*this)=(*this)*other;
357        return *this;
358}
359
[89]360
361//todo: use memcpy for these 4 functions
[86]362template<class T>
363Vector<T> Matrix<T>::GetRowVector(unsigned int r) const
364{
365        assert(r<m_Rows);
366        Vector<T> ret(m_Cols);
367        for (unsigned int j=0; j<m_Cols; j++)
368        {
369                ret[j]=(*this)[r][j];
370        }
371        return ret;
372}
[85]373
374template<class T>
[86]375Vector<T> Matrix<T>::GetColVector(unsigned int c) const
376{
377        assert(c<m_Cols);
378        Vector<T> ret(m_Rows);
379        for (unsigned int i=0; i<m_Rows; i++)
380        {
381                ret[i]=(*this)[i][c];
382        }
383        return ret;
384}
385
386template<class T>
[89]387void Matrix<T>::SetRowVector(unsigned int r, const Vector<T> &row)
388{
389        assert(r<m_Rows);
390        assert(row.Size()==m_Cols);
391        for (unsigned int j=0; j<m_Cols; j++)
392        {
393                (*this)[r][j]=row[j];
394        }
395}
396
397template<class T>
398void Matrix<T>::SetColVector(unsigned int c, const Vector<T> &col)
399{
400        assert(c<m_Cols);
401        assert(col.Size()==m_Rows);
402        for (unsigned int i=0; i<m_Rows; i++)
403        {
404                (*this)[i][c]=col[i];
405        }
406}
407
408// sort rows by v
409template<class T>
410void Matrix<T>::SortRows(Vector<T> &v)
411{
412        assert(v.Size()==m_Rows);
413       
414        bool sorted=false;
415        while(!sorted)
416        {
417                sorted=true;
418               
419                for (unsigned int i=0; i<v.Size()-1; i++)
420                {
421                        if (v[i]<v[i+1])
422                        {
423                                sorted=false;
424                                float vtmp = v[i];
425                                v[i]=v[i+1];
426                                v[i+1]=vtmp;
427                               
428                                Vector<float> rtmp = GetRowVector(i);
429                                SetRowVector(i,GetRowVector(i+1));
430                                SetRowVector(i+1,rtmp);
431                        }
432                }
433        }
434}
435
436// sort cols by v
437template<class T>
438void Matrix<T>::SortCols(Vector<T> &v)
439{
440        assert(v.Size()==m_Cols);
441       
442        bool sorted=false;
443        while(!sorted)
444        {
445                sorted=true;
446               
447                for (unsigned int i=0; i<v.Size()-1; i++)
448                {
449                        if (v[i]<v[i+1])
450                        {
451                                sorted=false;
452                                float vtmp = v[i];
453                                v[i]=v[i+1];
454                                v[i+1]=vtmp;
455                               
456                                Vector<float> rtmp = GetColVector(i);
457                                SetColVector(i,GetColVector(i+1));
458                                SetColVector(i+1,rtmp);
459                        }
460                }
461        }
462}
[90]463       
464template<class T>
465Matrix<T> Matrix<T>::CropRows(unsigned int s, unsigned int e)
466{
467        assert(s<e);
468        assert(s<m_Rows);
469        assert(e<=m_Rows);
470       
471        Matrix r(e-s,m_Cols);
472        unsigned int c=0;
473        for(unsigned int i=s; i<e; i++)
474        {
475                r.SetRowVector(c,GetRowVector(i));
476                c++;
477        }
478       
479        return r;
480}
[89]481
[90]482template<class T>
483Matrix<T> Matrix<T>::CropCols(unsigned int s, unsigned int e)
484{
485        assert(s<e);
486        assert(s<m_Cols);
487        assert(e<=m_Cols);
488       
489        Matrix r(m_Rows,e-s);
490        unsigned int c=0;
491        for(unsigned int i=s; i<e; i++)
492        {
493                r.SetColVector(c,GetColVector(i));
494                c++;
495        }
496       
497        return r;
498}
[89]499
500template<class T>
[90]501void Matrix<T>::Save(FILE* f)
502{
503        int version = 1;       
504        fwrite(&version,1,sizeof(version),f);
505        fwrite(&m_Rows,1,sizeof(m_Rows),f);
506        fwrite(&m_Cols,1,sizeof(m_Cols),f);
507        fwrite(m_Data,1,sizeof(T)*m_Rows*m_Cols,f);
508}
509
510template<class T>
511void Matrix<T>::Load(FILE* f)
512{
513        int version;   
514        fread(&version,sizeof(version),1,f);
515        fread(&m_Rows,sizeof(m_Rows),1,f);
516        fread(&m_Cols,sizeof(m_Cols),1,f);
517        m_Data=new T[m_Rows*m_Cols];
518        fread(m_Data,sizeof(T)*m_Rows*m_Cols,1,f);
519}
520
521template<class T>
[85]522void Matrix<T>::RunTests()
523{
[90]524        Vector<T>::RunTests();
525
[85]526        Matrix<T> m(10,10);
527        m.SetAll(0);
528        assert(m[0][0]==0);
529        m[5][2]=0.5;
530        assert(m[5][2]==0.5);
531        Matrix<T> om(m);
532        assert(om[5][2]==0.5);
533        Matrix<T> a(2,3);
534        a[0][0]=1; a[0][1]=2; a[0][2]=3;
535        a[1][0]=4; a[1][1]=5; a[1][2]=6;
536        Matrix<T> b(3,1);
537        b[0][0]=3;
538        b[1][0]=1;
539        b[2][0]=2;
540        Matrix<T> c=a*b;
541        assert(c[0][0]==11 && c[1][0]==29);
[90]542       
543        // test matrix * vector
544        Vector<T> d(3);
545        d[0]=3;
546        d[1]=1;
547        d[2]=2;
548        Vector<T> e=a*d;
549        assert(e[0]==11 && e[1]==29);
550       
551        Matrix<T> f=a.CropCols(1,3);
552        assert(f.GetRows()==2 && f.GetCols()==2 && f[0][0]==2);
553        Matrix<T> g=a.CropRows(0,1);
554        assert(g.GetRows()==1 && g.GetCols()==3 && g[0][0]==1);
555
[85]556}
557
558#endif
Note: See TracBrowser for help on using the repository browser.