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

Revision 90, 3.6 KB checked in by dave, 11 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) |

80 | { |

81 | return m_EigenTransform*v; |

82 | } |

83 | |

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

85 | { |

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

87 | } |

88 | |

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

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

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