FFmpeg
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  if (!pca)
45  return NULL;
46 
47  pca->n= n;
48  pca->z = av_malloc_array(n, sizeof(*pca->z));
49  pca->count=0;
50  pca->covariance= av_calloc(n*n, sizeof(double));
51  pca->mean= av_calloc(n, sizeof(double));
52 
53  if (!pca->z || !pca->covariance || !pca->mean) {
54  ff_pca_free(pca);
55  return NULL;
56  }
57 
58  return pca;
59 }
60 
61 void ff_pca_free(PCA *pca){
62  av_freep(&pca->covariance);
63  av_freep(&pca->mean);
64  av_freep(&pca->z);
65  av_free(pca);
66 }
67 
68 void ff_pca_add(PCA *pca, const double *v){
69  int i, j;
70  const int n= pca->n;
71 
72  for(i=0; i<n; i++){
73  pca->mean[i] += v[i];
74  for(j=i; j<n; j++)
75  pca->covariance[j + i*n] += v[i]*v[j];
76  }
77  pca->count++;
78 }
79 
80 int ff_pca(PCA *pca, double *eigenvector, double *eigenvalue){
81  int i, j, pass;
82  int k=0;
83  const int n= pca->n;
84  double *z = pca->z;
85 
86  memset(eigenvector, 0, sizeof(double)*n*n);
87 
88  for(j=0; j<n; j++){
89  pca->mean[j] /= pca->count;
90  eigenvector[j + j*n] = 1.0;
91  for(i=0; i<=j; i++){
92  pca->covariance[j + i*n] /= pca->count;
93  pca->covariance[j + i*n] -= pca->mean[i] * pca->mean[j];
94  pca->covariance[i + j*n] = pca->covariance[j + i*n];
95  }
96  eigenvalue[j]= pca->covariance[j + j*n];
97  z[j]= 0;
98  }
99 
100  for(pass=0; pass < 50; pass++){
101  double sum=0;
102 
103  for(i=0; i<n; i++)
104  for(j=i+1; j<n; j++)
105  sum += fabs(pca->covariance[j + i*n]);
106 
107  if(sum == 0){
108  for(i=0; i<n; i++){
109  double maxvalue= -1;
110  for(j=i; j<n; j++){
111  if(eigenvalue[j] > maxvalue){
112  maxvalue= eigenvalue[j];
113  k= j;
114  }
115  }
116  eigenvalue[k]= eigenvalue[i];
117  eigenvalue[i]= maxvalue;
118  for(j=0; j<n; j++){
119  double tmp= eigenvector[k + j*n];
120  eigenvector[k + j*n]= eigenvector[i + j*n];
121  eigenvector[i + j*n]= tmp;
122  }
123  }
124  return pass;
125  }
126 
127  for(i=0; i<n; i++){
128  for(j=i+1; j<n; j++){
129  double covar= pca->covariance[j + i*n];
130  double t,c,s,tau,theta, h;
131 
132  if(pass < 3 && fabs(covar) < sum / (5*n*n)) //FIXME why pass < 3
133  continue;
134  if(fabs(covar) == 0.0) //FIXME should not be needed
135  continue;
136  if(pass >=3 && fabs((eigenvalue[j]+z[j])/covar) > (1LL<<32) && fabs((eigenvalue[i]+z[i])/covar) > (1LL<<32)){
137  pca->covariance[j + i*n]=0.0;
138  continue;
139  }
140 
141  h= (eigenvalue[j]+z[j]) - (eigenvalue[i]+z[i]);
142  theta=0.5*h/covar;
143  t=1.0/(fabs(theta)+sqrt(1.0+theta*theta));
144  if(theta < 0.0) t = -t;
145 
146  c=1.0/sqrt(1+t*t);
147  s=t*c;
148  tau=s/(1.0+c);
149  z[i] -= t*covar;
150  z[j] += t*covar;
151 
152 #define ROTATE(a,i,j,k,l) {\
153  double g=a[j + i*n];\
154  double h=a[l + k*n];\
155  a[j + i*n]=g-s*(h+g*tau);\
156  a[l + k*n]=h+s*(g-h*tau); }
157  for(k=0; k<n; k++) {
158  if(k!=i && k!=j){
159  ROTATE(pca->covariance,FFMIN(k,i),FFMAX(k,i),FFMIN(k,j),FFMAX(k,j))
160  }
161  ROTATE(eigenvector,k,i,k,j)
162  }
163  pca->covariance[j + i*n]=0.0;
164  }
165  }
166  for (i=0; i<n; i++) {
167  eigenvalue[i] += z[i];
168  z[i]=0.0;
169  }
170  }
171 
172  return -1;
173 }
Definition: pca.c:30
#define NULL
Definition: coverity.c:32
int count
Definition: pca.c:31
void * av_mallocz(size_t size)
Allocate a memory block with alignment suitable for all memory accesses (including vectors if availab...
Definition: mem.c:236
PCA * ff_pca_init(int n)
Definition: pca.c:38
double * z
Definition: pca.c:35
void * av_calloc(size_t nmemb, size_t size)
Non-inlined equivalent of av_mallocz_array().
Definition: mem.c:244
Undefined Behavior In the C some operations are like signed integer dereferencing freed accessing outside allocated Undefined Behavior must not occur in a C it is not safe even if the output of undefined operations is unused The unsafety may seem nit picking but Optimizing compilers have in fact optimized code on the assumption that no undefined Behavior occurs Optimizing code based on wrong assumptions can and has in some cases lead to effects beyond the output of computations The signed integer overflow problem in speed critical code Code which is highly optimized and works with signed integers sometimes has the problem that often the output of the computation does not c
Definition: undefined.txt:32
int n
Definition: pca.c:32
#define i(width, name, range_min, range_max)
Definition: cbs_h2645.c:259
#define FFMAX(a, b)
Definition: common.h:94
#define pass
Definition: fft_template.c:619
#define ROTATE(a, i, j, k, l)
#define FFMIN(a, b)
Definition: common.h:96
#define s(width, name)
Definition: cbs_vp9.c:257
int ff_pca(PCA *pca, double *eigenvector, double *eigenvalue)
Definition: pca.c:80
void ff_pca_free(PCA *pca)
Definition: pca.c:61
void ff_pca_add(PCA *pca, const double *v)
Definition: pca.c:68
principal component analysis (PCA)
double * mean
Definition: pca.c:34
common internal and external API header
#define av_free(p)
#define av_freep(p)
#define av_malloc_array(a, b)
double * covariance
Definition: pca.c:33
static uint8_t tmp[11]
Definition: aes_ctr.c:26