path: root/libavutil/pca.c (plain)
blob: 4e52c7b362cc16f4694952ce8e4614b87cd758d0

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	}
174