source: foam/trunk/vision/src/PCA.cpp @ 90

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

PCA updated - projection and resynthesis seem to work

Line 
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 "PCA.h"
18#include "SVD.h"
19#include <iostream>
20
21using namespace std;
22
23PCA::PCA(unsigned int FeatureSize) :
24m_FeatureSize(FeatureSize),
25m_Mean(FeatureSize)
26{
27}
28
29PCA::~PCA()
30{
31}
32
33
34void PCA::Calculate()
35{
36        // calculate the mean
37        m_Mean.Zero();
38        for (FeatureVec::iterator vi = m_Features.begin(); vi!=m_Features.end(); ++vi)
39        {
40                m_Mean+=*vi;
41        }
42       
43        m_Mean/=m_Features.size();
44               
45        // subtract the mean
46        FeatureVec SubMean;
47        for (FeatureVec::iterator vi = m_Features.begin(); vi!=m_Features.end(); ++vi)
48        {
49                SubMean.push_back(*vi-m_Mean);
50        }
51       
52        // allocate the transform matrix (this is where it'll run out of memory)
53        cerr<<"Allocating "<<m_FeatureSize*m_FeatureSize*sizeof(float)/1024/1024.0<<" megs for covariance matrix"<<endl;
54        m_EigenTransform = Matrix<float>(m_FeatureSize,m_FeatureSize);
55        m_EigenTransform.Zero();
56       
57        // start by calculating the covariance matrix
58        for (unsigned int i=0; i<m_FeatureSize; i++)
59        {
60                for (unsigned int j=0; j<m_FeatureSize; j++)
61                {
62                        for (FeatureVec::iterator f = SubMean.begin(); f!=SubMean.end(); ++f)
63                        {
64                                m_EigenTransform[i][j]+=(*f)[i]*(*f)[j];
65                        }
66                       
67                        m_EigenTransform[i][j]/=(float)(SubMean.size()-1);
68                }
69        }
70        m_EigenValues = SVD(m_EigenTransform);
71        m_EigenTransform=m_EigenTransform.Transposed();
72}
73
74void PCA::Compress(unsigned int s, unsigned int e)
75{
76        m_EigenTransform=m_EigenTransform.CropRows(s,e);
77}
78
79Vector<float> PCA::Project(Vector<float> v)
80{
81        return m_EigenTransform*v;
82}
83       
84Vector<float> PCA::Synth(Vector<float> v)
85{
86        return m_Mean+m_EigenTransform.VecMulTransposed(v);
87}
88
89void PCA::Save(FILE *f)
90{
91        int version = 2;
92        fwrite(&version,sizeof(version),1,f);
93        m_EigenTransform.Save(f);
94        m_EigenValues.Save(f);
95        m_Mean.Save(f);
96        fwrite(&m_FeatureSize,sizeof(m_FeatureSize),1,f);
97        for (unsigned int i=0; i<m_Features.size(); i++)
98        {
99                m_Features[i].Save(f);
100        }
101       
102}
103
104void PCA::Load(FILE *f)
105{
106        int version;   
107        fread(&version,sizeof(version),1,f);
108        m_EigenTransform.Load(f);
109        m_EigenValues.Load(f);
110       
111        if (version == 1)
112        {
113                m_EigenTransform=m_EigenTransform.Transposed();
114        }
115       
116        if (version>1)
117        {
118                m_Mean.Load(f);
119                fread(&m_FeatureSize,sizeof(m_FeatureSize),1,f);
120                for (unsigned int i=0; i<m_Features.size(); i++)
121                {
122                        m_Features[i].Load(f);
123                }
124        }
125}
126
127void PCA::RunTests()
128{
129        Matrix<float>::RunTests();
130        PCA pca(2);
131       
132        Vector<float> in(2);
133        in[0]=2.5; in[1]=2.4;
134        pca.AddFeature(in);
135        in[0]=0.5; in[1]=0.7;
136        pca.AddFeature(in);
137        in[0]=2.2; in[1]=2.9;
138        pca.AddFeature(in);
139        in[0]=1.9; in[1]=2.2;
140        pca.AddFeature(in);
141        in[0]=3.1; in[1]=3.0;
142        pca.AddFeature(in);
143        in[0]=2.3; in[1]=2.7;
144        pca.AddFeature(in);
145        in[0]=2; in[1]=1.6;
146        pca.AddFeature(in);
147        in[0]=1; in[1]=1.1;
148        pca.AddFeature(in);
149        in[0]=1.5; in[1]=1.6;
150        pca.AddFeature(in);
151        in[0]=1.1; in[1]=0.9;
152        pca.AddFeature(in);
153       
154        pca.Calculate();
155               
156        in[0]=.69; in[1]=.49;
157        Vector<float> out = pca.Project(in);
158        assert(feq(out[0],-0.82797f) && feq(out[1],-0.175115f));
159       
160}
161
162
163
Note: See TracBrowser for help on using the repository browser.