FFmpeg: libavutil/pca.c Source File

FFmpeg

[フレーム]

libavutil

pca.c

Go to the documentation of this file.

1 /*

2 * principal component analysis (PCA)

4 *

5 * This file is part of FFmpeg.

6 *

7 * FFmpeg is free software; you can redistribute it and/or

8 * modify it under the terms of the GNU Lesser General Public

9 * License as published by the Free Software Foundation; either

10 * version 2.1 of the License, or (at your option) any later version.

11 *

12 * FFmpeg is distributed in the hope that it will be useful,

13 * but WITHOUT ANY WARRANTY; without even the implied warranty of

14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU

15 * Lesser General Public License for more details.

16 *

17 * You should have received a copy of the GNU Lesser General Public

18 * License along with FFmpeg; if not, write to the Free Software

19 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA

20 */

22 /**

23 * @file

24 * principal component analysis (PCA)

25 */

27 #include "common.h"

28 #include "pca.h"

30 typedef struct PCA{

31 int count;

32 int n;

33 double *covariance;

34 double *mean;

35 double *z;

36 }PCA;

38 PCA *ff_pca_init(int n){

39 PCA *pca;

40 if(n<=0)

41 return NULL;

43 pca= av_mallocz(sizeof(*pca));

44 pca->n= n;

45 pca->z = av_malloc(sizeof(*pca->z) * n);

46 pca->count=0;

47 pca->covariance= av_calloc(n*n, sizeof(double));

48 pca->mean= av_calloc(n, sizeof(double));

50 return pca;

51 }

53 void ff_pca_free(PCA *pca){

54 av_freep(&pca->covariance);

55 av_freep(&pca->mean);

56 av_freep(&pca->z);

57 av_free(pca);

58 }

60 void ff_pca_add(PCA *pca, double *v){

61 int i, j;

62 const int n= pca->n;

64 for(i=0; i<n; i++){

65 pca->mean[i] += v[i];

66 for(j=i; j<n; j++)

67 pca->covariance[j + i*n] += v[i]*v[j];

68 }

69 pca->count++;

70 }

72 int ff_pca(PCA *pca, double *eigenvector, double *eigenvalue){

73 int i, j, pass;

74 int k=0;

75 const int n= pca->n;

76 double *z = pca->z;

78 memset(eigenvector, 0, sizeof(double)*n*n);

80 for(j=0; j<n; j++){

81 pca->mean[j] /= pca->count;

82 eigenvector[j + j*n] = 1.0;

83 for(i=0; i<=j; i++){

84 pca->covariance[j + i*n] /= pca->count;

85 pca->covariance[j + i*n] -= pca->mean[i] * pca->mean[j];

86 pca->covariance[i + j*n] = pca->covariance[j + i*n];

87 }

88 eigenvalue[j]= pca->covariance[j + j*n];

89 z[j]= 0;

90 }

92 for(pass=0; pass < 50; pass++){

93 double sum=0;

95 for(i=0; i<n; i++)

96 for(j=i+1; j<n; j++)

97 sum += fabs(pca->covariance[j + i*n]);

99 if(sum == 0){

100 for(i=0; i<n; i++){

101 double maxvalue= -1;

102 for(j=i; j<n; j++){

103 if(eigenvalue[j] > maxvalue){

104 maxvalue= eigenvalue[j];

105 k= j;

106 }

107 }

108 eigenvalue[k]= eigenvalue[i];

109 eigenvalue[i]= maxvalue;

110 for(j=0; j<n; j++){

111 double tmp= eigenvector[k + j*n];

112 eigenvector[k + j*n]= eigenvector[i + j*n];

113 eigenvector[i + j*n]= tmp;

114 }

115 }

116 return pass;

117 }

118

119 for(i=0; i<n; i++){

120 for(j=i+1; j<n; j++){

121 double covar= pca->covariance[j + i*n];

122 double t,c,s,tau,theta, h;

123

124 if(pass < 3 && fabs(covar) < sum / (5*n*n)) //FIXME why pass < 3

125 continue;

126 if(fabs(covar) == 0.0) //FIXME should not be needed

127 continue;

128 if(pass >=3 && fabs((eigenvalue[j]+z[j])/covar) > (1LL<<32) && fabs((eigenvalue[i]+z[i])/covar) > (1LL<<32)){

129 pca->covariance[j + i*n]=0.0;

130 continue;

131 }

132

133 h= (eigenvalue[j]+z[j]) - (eigenvalue[i]+z[i]);

134 theta=0.5*h/covar;

135 t=1.0/(fabs(theta)+sqrt(1.0+theta*theta));

136 if(theta < 0.0) t = -t;

137

138 c=1.0/sqrt(1+t*t);

139 s=t*c;

140 tau=s/(1.0+c);

141 z[i] -= t*covar;

142 z[j] += t*covar;

143

144 #define ROTATE(a,i,j,k,l) {\

145 double g=a[j + i*n];\

146 double h=a[l + k*n];\

147 a[j + i*n]=g-s*(h+g*tau);\

148 a[l + k*n]=h+s*(g-h*tau); }

149 for(k=0; k<n; k++) {

150 if(k!=i && k!=j){

151 ROTATE(pca->covariance,FFMIN(k,i),FFMAX(k,i),FFMIN(k,j),FFMAX(k,j))

152 }

153 ROTATE(eigenvector,k,i,k,j)

154 }

155 pca->covariance[j + i*n]=0.0;

156 }

157 }

158 for (i=0; i<n; i++) {

159 eigenvalue[i] += z[i];

160 z[i]=0.0;

161 }

162 }

163

164 return -1;

165 }

166

167 #ifdef TEST

168

169 #undef printf

170 #include <stdio.h>

171 #include <stdlib.h>

172 #include "lfg.h"

173

174 int main(void){

175 PCA *pca;

176 int i, j, k;

177 #define LEN 8

178 double eigenvector[LEN*LEN];

179 double eigenvalue[LEN];

180 AVLFG prng;

181

182 av_lfg_init(&prng, 1);

183

184 pca= ff_pca_init(LEN);

185

186 for(i=0; i<9000000; i++){

187 double v[2*LEN+100];

188 double sum=0;

189 int pos = av_lfg_get(&prng) % LEN;

190 int v2 = av_lfg_get(&prng) % 101 - 50;

191 v[0] = av_lfg_get(&prng) % 101 - 50;

192 for(j=1; j<8; j++){

193 if(j<=pos) v[j]= v[0];

194 else v[j]= v2;

195 sum += v[j];

196 }

197 /* for(j=0; j<LEN; j++){

198 v[j] -= v[pos];

199 }*/

200 // sum += av_lfg_get(&prng) % 10;

201 /* for(j=0; j<LEN; j++){

202 v[j] -= sum/LEN;

203 }*/

204 // lbt1(v+100,v+100,LEN);

205 ff_pca_add(pca, v);

206 }

207

208

209 ff_pca(pca, eigenvector, eigenvalue);

210 for(i=0; i<LEN; i++){

211 pca->count= 1;

212 pca->mean[i]= 0;

213

214 // (0.5^|x|)^2 = 0.5^2|x| = 0.25^|x|

215

216

217 // pca.covariance[i + i*LEN]= pow(0.5, fabs

218 for(j=i; j<LEN; j++){

219 printf("%f ", pca->covariance[i + j*LEN]);

220 }

221 printf("\n");

222 }

223

224 for(i=0; i<LEN; i++){

225 double v[LEN];

226 double error=0;

227 memset(v, 0, sizeof(v));

228 for(j=0; j<LEN; j++){

229 for(k=0; k<LEN; k++){

230 v[j] += pca->covariance[FFMIN(k,j) + FFMAX(k,j)*LEN] * eigenvector[i + k*LEN];

231 }

232 v[j] /= eigenvalue[i];

233 error += fabs(v[j] - eigenvector[i + j*LEN]);

234 }

235 printf("%f ", error);

236 }

237 printf("\n");

238

239 for(i=0; i<LEN; i++){

240 for(j=0; j<LEN; j++){

241 printf("%9.6f ", eigenvector[i + j*LEN]);

242 }

243 printf(" %9.1f %f\n", eigenvalue[i], eigenvalue[i]/eigenvalue[0]);

244 }

245

246 return 0;

247 }

248 #endif

Generated on Wed Jul 10 2013 23:48:16 for FFmpeg by doxygen 1.8.2