FFmpeg: libavutil/pca.c Source File

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 
38  PCA *ff_pca_init(int n){
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 }
tmp
static uint8_t tmp[11]
Definition: aes_ctr.c:28
PCA::z
double * z
Definition: pca.c:35
FFMAX
#define FFMAX(a, b)
Definition: macros.h:47
PCA::n
int n
Definition: pca.c:32
PCA::count
int count
Definition: pca.c:31
ff_pca_init
PCA * ff_pca_init(int n)
Definition: pca.c:38
s
#define s(width, name)
Definition: cbs_vp9.c:256
PCA
Definition: pca.c:30
pass
#define pass
Definition: fft_template.c:608
ROTATE
#define ROTATE(a, i, j, k, l)
ff_pca_free
void ff_pca_free(PCA *pca)
Definition: pca.c:61
fabs
static __device__ float fabs(float a)
Definition: cuda_runtime.h:182
NULL
#define NULL
Definition: coverity.c:32
c
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
pca.h
i
#define i(width, name, range_min, range_max)
Definition: cbs_h2645.c:269
av_malloc_array
#define av_malloc_array(a, b)
Definition: tableprint_vlc.h:31
common.h
PCA::mean
double * mean
Definition: pca.c:34
FFMIN
#define FFMIN(a, b)
Definition: macros.h:49
av_mallocz
void * av_mallocz(size_t size)
Allocate a memory block with alignment suitable for all memory accesses (including vectors if availab...
Definition: mem.c:254
av_calloc
void * av_calloc(size_t nmemb, size_t size)
Definition: mem.c:262
PCA::covariance
double * covariance
Definition: pca.c:33
ff_pca
int ff_pca(PCA *pca, double *eigenvector, double *eigenvalue)
Definition: pca.c:80
av_free
#define av_free(p)
Definition: tableprint_vlc.h:33
av_freep
#define av_freep(p)
Definition: tableprint_vlc.h:34
h
h
Definition: vp9dsp_template.c:2038
ff_pca_add
void ff_pca_add(PCA *pca, const double *v)
Definition: pca.c:68
Generated on Tue Feb 28 2023 21:34:14 for FFmpeg by   doxygen 1.8.17

AltStyle によって変換されたページ (->オリジナル) /