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

Revision 96, 3.7 KB checked in by dave, 10 years ago (diff)

Benchmark testing added for the eigen face recognition

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) const
80{
81        return m_EigenTransform*v;
82}
83       
84Vector<float> PCA::Synth(Vector<float> v) const
85{
86        return m_Mean+m_EigenTransform.VecMulTransposed(v);
87}
88       
89void PCA::Mult(const PCA &other)
90{
91        m_EigenTransform *= other.GetEigenTransform().Inverted();
92}
93
94void PCA::Save(FILE *f)
95{
96        int version = 2;
97        fwrite(&version,sizeof(version),1,f);
98        m_EigenTransform.Save(f);
99        m_EigenValues.Save(f);
100        m_Mean.Save(f);
101        fwrite(&m_FeatureSize,sizeof(m_FeatureSize),1,f);
102        for (unsigned int i=0; i<m_Features.size(); i++)
103        {
104                m_Features[i].Save(f);
105        }
106       
107}
108
109void PCA::Load(FILE *f)
110{
111        int version;   
112        fread(&version,sizeof(version),1,f);
113        m_EigenTransform.Load(f);
114        m_EigenValues.Load(f);
115       
116        if (version == 1)
117        {
118                m_EigenTransform=m_EigenTransform.Transposed();
119        }
120       
121        if (version>1)
122        {
123                m_Mean.Load(f);
124                fread(&m_FeatureSize,sizeof(m_FeatureSize),1,f);
125                for (unsigned int i=0; i<m_Features.size(); i++)
126                {
127                        m_Features[i].Load(f);
128                }
129        }
130}
131
132void PCA::RunTests()
133{
134        Matrix<float>::RunTests();
135        PCA pca(2);
136       
137        Vector<float> in(2);
138        in[0]=2.5; in[1]=2.4;
139        pca.AddFeature(in);
140        in[0]=0.5; in[1]=0.7;
141        pca.AddFeature(in);
142        in[0]=2.2; in[1]=2.9;
143        pca.AddFeature(in);
144        in[0]=1.9; in[1]=2.2;
145        pca.AddFeature(in);
146        in[0]=3.1; in[1]=3.0;
147        pca.AddFeature(in);
148        in[0]=2.3; in[1]=2.7;
149        pca.AddFeature(in);
150        in[0]=2; in[1]=1.6;
151        pca.AddFeature(in);
152        in[0]=1; in[1]=1.1;
153        pca.AddFeature(in);
154        in[0]=1.5; in[1]=1.6;
155        pca.AddFeature(in);
156        in[0]=1.1; in[1]=0.9;
157        pca.AddFeature(in);
158       
159        pca.Calculate();
160               
161        in[0]=.69; in[1]=.49;
162        Vector<float> out = pca.Project(in);
163        assert(feq(out[0],-0.82797f) && feq(out[1],-0.175115f));
164       
165}
166
167
168
Note: See TracBrowser for help on using the repository browser.