FFmpeg
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Groups Pages
pca.c
Go to the documentation of this file.
1 /*
2  * principal component analysis (PCA)
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 02110-1301 USA
20  */
21 
22 /**
23  * @file
24  * principal component analysis (PCA)
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  double *z;
36 }PCA;
37 
39  PCA *pca;
40  if(n<=0)
41  return NULL;
42 
43  pca= av_mallocz(sizeof(*pca));
44  pca->n= n;
45  pca->z = av_malloc_array(n, sizeof(*pca->z));
46  pca->count=0;
47  pca->covariance= av_calloc(n*n, sizeof(double));
48  pca->mean= av_calloc(n, sizeof(double));
49 
50  return pca;
51 }
52 
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 }
59 
60 void ff_pca_add(PCA *pca, const double *v){
61  int i, j;
62  const int n= pca->n;
63 
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 }
71 
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;
77 
78  memset(eigenvector, 0, sizeof(double)*n*n);
79 
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  }
91 
92  for(pass=0; pass < 50; pass++){
93  double sum=0;
94 
95  for(i=0; i<n; i++)
96  for(j=i+1; j<n; j++)
97  sum += fabs(pca->covariance[j + i*n]);
98 
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