ffmpeg / libavutil / pca.c @ b76e3424
History  View  Annotate  Download (6.24 KB)
1 
/*


2 
* Principal component analysis

3 
* Copyright (c) 2004 Michael Niedermayer <michaelni@gmx.at>

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

20 
*/

21  
22 
/**

23 
* @file pca.c

24 
* Principal component analysis

25 
*/

26  
27 
#include "common.h" 
28 
#include "pca.h" 
29  
30 
typedef struct PCA{ 
31 
int count;

32 
int n;

33 
double *covariance;

34 
double *mean;

35 
}PCA; 
36  
37 
PCA *ff_pca_init(int n){

38 
PCA *pca; 
39 
if(n<=0) 
40 
return NULL; 
41  
42 
pca= av_mallocz(sizeof(PCA));

43 
pca>n= n; 
44 
pca>count=0;

45 
pca>covariance= av_mallocz(sizeof(double)*n*n); 
46 
pca>mean= av_mallocz(sizeof(double)*n); 
47  
48 
return pca;

49 
} 
50  
51 
void ff_pca_free(PCA *pca){

52 
av_freep(&pca>covariance); 
53 
av_freep(&pca>mean); 
54 
av_free(pca); 
55 
} 
56  
57 
void ff_pca_add(PCA *pca, double *v){ 
58 
int i, j;

59 
const int n= pca>n; 
60  
61 
for(i=0; i<n; i++){ 
62 
pca>mean[i] += v[i]; 
63 
for(j=i; j<n; j++)

64 
pca>covariance[j + i*n] += v[i]*v[j]; 
65 
} 
66 
pca>count++; 
67 
} 
68  
69 
int ff_pca(PCA *pca, double *eigenvector, double *eigenvalue){ 
70 
int i, j, k, pass;

71 
const int n= pca>n; 
72 
double z[n];

73  
74 
memset(eigenvector, 0, sizeof(double)*n*n); 
75  
76 
for(j=0; j<n; j++){ 
77 
pca>mean[j] /= pca>count; 
78 
eigenvector[j + j*n] = 1.0; 
79 
for(i=0; i<=j; i++){ 
80 
pca>covariance[j + i*n] /= pca>count; 
81 
pca>covariance[j + i*n] = pca>mean[i] * pca>mean[j]; 
82 
pca>covariance[i + j*n] = pca>covariance[j + i*n]; 
83 
} 
84 
eigenvalue[j]= pca>covariance[j + j*n]; 
85 
z[j]= 0;

86 
} 
87  
88 
for(pass=0; pass < 50; pass++){ 
89 
double sum=0; 
90  
91 
for(i=0; i<n; i++) 
92 
for(j=i+1; j<n; j++) 
93 
sum += fabs(pca>covariance[j + i*n]); 
94  
95 
if(sum == 0){ 
96 
for(i=0; i<n; i++){ 
97 
double maxvalue= 1; 
98 
for(j=i; j<n; j++){

99 
if(eigenvalue[j] > maxvalue){

100 
maxvalue= eigenvalue[j]; 
101 
k= j; 
102 
} 
103 
} 
104 
eigenvalue[k]= eigenvalue[i]; 
105 
eigenvalue[i]= maxvalue; 
106 
for(j=0; j<n; j++){ 
107 
double tmp= eigenvector[k + j*n];

108 
eigenvector[k + j*n]= eigenvector[i + j*n]; 
109 
eigenvector[i + j*n]= tmp; 
110 
} 
111 
} 
112 
return pass;

113 
} 
114  
115 
for(i=0; i<n; i++){ 
116 
for(j=i+1; j<n; j++){ 
117 
double covar= pca>covariance[j + i*n];

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

119  
120 
if(pass < 3 && fabs(covar) < sum / (5*n*n)) //FIXME why pass < 3 
121 
continue;

122 
if(fabs(covar) == 0.0) //FIXME shouldnt be needed 
123 
continue;

124 
if(pass >=3 && fabs((eigenvalue[j]+z[j])/covar) > (1LL<<32) && fabs((eigenvalue[i]+z[i])/covar) > (1LL<<32)){ 
125 
pca>covariance[j + i*n]=0.0; 
126 
continue;

127 
} 
128  
129 
h= (eigenvalue[j]+z[j])  (eigenvalue[i]+z[i]); 
130 
theta=0.5*h/covar; 
131 
t=1.0/(fabs(theta)+sqrt(1.0+theta*theta)); 
132 
if(theta < 0.0) t = t; 
133  
134 
c=1.0/sqrt(1+t*t); 
135 
s=t*c; 
136 
tau=s/(1.0+c); 
137 
z[i] = t*covar; 
138 
z[j] += t*covar; 
139  
140 
#define ROTATE(a,i,j,k,l) {\

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

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

143 
a[j + i*n]=gs*(h+g*tau);\ 
144 
a[l + k*n]=h+s*(gh*tau); } 
145 
for(k=0; k<n; k++) { 
146 
if(k!=i && k!=j){

147 
ROTATE(pca>covariance,FFMIN(k,i),FFMAX(k,i),FFMIN(k,j),FFMAX(k,j)) 
148 
} 
149 
ROTATE(eigenvector,k,i,k,j) 
150 
} 
151 
pca>covariance[j + i*n]=0.0; 
152 
} 
153 
} 
154 
for (i=0; i<n; i++) { 
155 
eigenvalue[i] += z[i]; 
156 
z[i]=0.0; 
157 
} 
158 
} 
159  
160 
return 1; 
161 
} 
162  
163 
#ifdef TEST

164  
165 
#undef printf

166 
#undef random

167 
#include <stdio.h> 
168 
#include <stdlib.h> 
169  
170 
int main(){

171 
PCA *pca; 
172 
int i, j, k;

173 
#define LEN 8 
174 
double eigenvector[LEN*LEN];

175 
double eigenvalue[LEN];

176  
177 
pca= ff_pca_init(LEN); 
178  
179 
for(i=0; i<9000000; i++){ 
180 
double v[2*LEN+100]; 
181 
double sum=0; 
182 
int pos= random()%LEN;

183 
int v2= (random()%101)  50; 
184 
v[0]= (random()%101)  50; 
185 
for(j=1; j<8; j++){ 
186 
if(j<=pos) v[j]= v[0]; 
187 
else v[j]= v2;

188 
sum += v[j]; 
189 
} 
190 
/* for(j=0; j<LEN; j++){

191 
v[j] = v[pos];

192 
}*/

193 
// sum += random()%10;

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

195 
v[j] = sum/LEN;

196 
}*/

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

198 
ff_pca_add(pca, v); 
199 
} 
200  
201  
202 
ff_pca(pca, eigenvector, eigenvalue); 
203 
for(i=0; i<LEN; i++){ 
204 
pca>count= 1;

205 
pca>mean[i]= 0;

206  
207 
// (0.5^x)^2 = 0.5^2x = 0.25^x

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

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

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

213 
} 
214 
printf("\n");

215 
} 
216  
217 
#if 1 
218 
for(i=0; i<LEN; i++){ 
219 
double v[LEN];

220 
double error=0; 
221 
memset(v, 0, sizeof(v)); 
222 
for(j=0; j<LEN; j++){ 
223 
for(k=0; k<LEN; k++){ 
224 
v[j] += pca>covariance[FFMIN(k,j) + FFMAX(k,j)*LEN] * eigenvector[i + k*LEN]; 
225 
} 
226 
v[j] /= eigenvalue[i]; 
227 
error += fabs(v[j]  eigenvector[i + j*LEN]); 
228 
} 
229 
printf("%f ", error);

230 
} 
231 
printf("\n");

232 
#endif

233 
for(i=0; i<LEN; i++){ 
234 
for(j=0; j<LEN; j++){ 
235 
printf("%9.6f ", eigenvector[i + j*LEN]);

236 
} 
237 
printf(" %9.1f %f\n", eigenvalue[i], eigenvalue[i]/eigenvalue[0]); 
238 
} 
239  
240 
return 0; 
241 
} 
242 
#endif
