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

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

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 | |

21 | using namespace std; |

22 | |

23 | PCA::PCA(unsigned int FeatureSize) : |

24 | m_FeatureSize(FeatureSize), |

25 | m_Mean(FeatureSize) |

26 | { |

27 | } |

28 | |

29 | PCA::~PCA() |

30 | { |

31 | } |

32 | |

33 | |

34 | void 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 | |

74 | void PCA::Compress(unsigned int s, unsigned int e) |

75 | { |

76 | m_EigenTransform=m_EigenTransform.CropRows(s,e); |

77 | } |

78 | |

79 | Vector<float> PCA::Project(Vector<float> v) const |

80 | { |

81 | return m_EigenTransform*v; |

82 | } |

83 | |

84 | Vector<float> PCA::Synth(Vector<float> v) const |

85 | { |

86 | return m_Mean+m_EigenTransform.VecMulTransposed(v); |

87 | } |

88 | |

89 | void PCA::Mult(const PCA &other) |

90 | { |

91 | m_EigenTransform *= other.GetEigenTransform().Inverted(); |

92 | } |

93 | |

94 | void 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 | |

109 | void 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 | |

132 | void 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.