Context Navigation

source: proiecte/PPPP/eigenface2/eigenface_p.c @ 132

Last change on this file since 132 was 132, checked in by (none), 14 years ago
PPPP
File size: 10.6 KB

Line
1	#include <stdio.h>
2	#include <string.h>
3	#include <math.h>
4	#include "cv.h"
5	#include "cvaux.h"
6	#include "highgui.h"
7	#include <omp.h>
8
9	#define ROTATE(a,i,j,k,l) g=a[in + j];h=a[kn + l];a[in + j]=g-s(h+g*tau);\
10	a[kn + l]=h+s(g-h*tau);
11
12	/*
13	Computes all eigenvalues and eigenvectors of a real symmetric matrix a[1..n][1..n]. On
14	output, elements of a above the diagonal are destroyed. d[1..n] returns the eigenvalues of a.
15	v[1..n][1..n] is a matrix whose columns contain, on output, the normalized eigenvectors of
16	a. nrot returns the number of Jacobi rotations that were required.
17	*/
18	void jacobi(float a, int n, float d[], float v, int *nrot)
19	{
20	int j,iq,ip,i;
21	float tresh,theta,tau,t,sm,s,h,g,c,b,z;
22	int tid, nthreads;
23
24	b = (float ) malloc(n sizeof(float));
25	z = (float ) malloc(n sizeof(float));
26
27	/* Fork a team of threads giving them their own copies of variables */
28	#pragma omp parallel private(nthreads, tid, i, j, ip, iq) shared(tresh, theta, tau, sm, s, h, g, c, b, z, a, n, d, v, nrot)
29	{
30
31	/* Obtain thread number */
32	tid = omp_get_thread_num();
33	printf("Hello World from thread = %d\n", tid);
34
35	/* Only master thread does this */
36	if (tid == 0)
37	{
38	nthreads = omp_get_num_threads();
39	printf("Number of threads = %d\n", nthreads);
40	}
41
42
43	#pragma omp parallel for
44	for (ip=0;ip<n;ip++) { // Initialize to the identity matrix.
45	for (iq=0;iq<n;iq++) v[ip*n + iq]=0.0;
46	v[ip*n + ip]=1.0;
47	}
48	#pragma omp parallel for
49	for (ip=0;ip<n;ip++) { // Initialize b and d to the diagonal of a.
50	b[ip]=d[ip]=a[ip*n + ip];
51	z[ip]=0.0; // This vector will accumulate terms of the form t*a[pq] as in equation (11.1.14).
52	}
53	}
54	*nrot=0;
55	for (i=1;i<=50;i++) {
56	sm=0.0;
57	#pragma omp parallel for reduction(+:sm)
58	for (ip=0;ip<n-1;ip++) { // Sum off-diagonal elements.
59	for (iq=ip+1;iq<n;iq++)
60	sm += fabs(a[ip*n + iq]);
61	}
62	if (sm == 0.0) { // The normal return, which relies on quadratic convergence to machine underflow.
63	free(z);
64	free(b);
65	return;
66	}
67	if (i < 4)
68	tresh=0.2sm/(nn); // ...on the first three sweeps.
69	else
70	tresh=0.0; // ...thereafter.
71	for (ip=0;ip<n-1;ip++) {
72	for (iq=ip+1;iq<n;iq++) {
73	g=100.0fabs(a[ipn + iq]);
74	if (i > 4 && (float)(fabs(d[ip])+g) == (float)fabs(d[ip]) // After four sweeps, skip the rotation if the off-diagonal element is small.
75	&& (float)(fabs(d[iq])+g) == (float)fabs(d[iq]))
76	a[ip*n + iq]=0.0;
77	else if (fabs(a[ip*n + iq]) > tresh) {
78	h=d[iq]-d[ip];
79	if ((float)(fabs(h)+g) == (float)fabs(h))
80	t=(a[ipn + iq])/h; // t = 1/(2theta)
81	else {
82	theta=0.5h/(a[ipn + iq]); // Equation (11.1.10).
83	t=1.0/(fabs(theta)+sqrt(1.0+theta*theta));
84	if (theta < 0.0) t = -t;
85	}
86
87	c=1.0/sqrt(1+t*t);
88	s=t*c;
89	tau=s/(1.0+c);
90	h=ta[ipn + iq];
91	z[ip] -= h;
92	z[iq] += h;
93	d[ip] -= h;
94	d[iq] += h;
95	a[ip*n + iq]=0.0;
96
97	#pragma omp parallel for
98	for (j=0;j<=ip-1;j++) { // Case of rotations 1 <= j < p.
99	ROTATE(a,j,ip,j,iq)
100	}
101	#pragma omp parallel for
102	for (j=ip+1;j<=iq-1;j++) { // Case of rotations p < j < q.
103	ROTATE(a,ip,j,j,iq)
104	}
105	#pragma omp parallel for
106	for (j=iq+1;j<n;j++) { // Case of rotations q < j <= n.
107	ROTATE(a,ip,j,iq,j)
108	}
109	#pragma omp parallel for
110	for (j=0;j<n;j++) {
111	ROTATE(v,j,ip,j,iq)
112	}
113
114	++(*nrot);
115	}
116	}
117	}
118
119	#pragma omp parallel for
120	for (ip=0;ip<n;ip++) {
121	b[ip] += z[ip];
122	d[ip]=b[ip]; // Updte d with the sum of t*a[pq],
123	z[ip]=0.0; // and reinitialize z.
124	}
125	}
126	}
127	/*
128	Given the eigenvalues d[1..n] and eigenvectors v[1..n][1..n] as output from jacobi
129	(x11.1) or tqli (x11.3), this routine sorts the eigenvalues into descending order, and rearranges
130	the columns of v correspondingly. The method is straight insertion.
131	*/
132	void eigsrt(float d[], float *v, int n)
133	{
134	int k,j,i;
135	float p;
136	for (i=0;i<n-1;i++) {
137	p=d[k=i];
138	for (j=i+1;j<n;j++)
139	if (d[j] >= p) p=d[k=j];
140	if (k != i) {
141	d[k]=d[i];
142	d[i]=p;
143	for (j=0;j<n;j++) {
144	p=v[j*n + i];
145	v[jn + i]=v[jn + k];
146	v[j*n + k]=p;
147	}
148	}
149	}
150	}
151
152	void calcMeanImage(int nFaces, uchar** faceArr, int faceStep,
153	CvSize size, float* avg, int avgStep)
154	{
155	int i,j,k;
156	float m = 1.0f / (float) nFaces;
157	float* bf = avg;
158
159	#pragma omp parallel for
160	for( i = 0; i < size.height; i++, bf += avgStep)
161	for( j = 0; j < size.width; j++ )
162	bf[j] = 0.f;
163
164	#pragma omp parallel for
165	for( i = 0; i < nFaces; i++ )
166	{
167	uchar* bu = faceArr[i];
168	bf = avg;
169	for( k = 0; k < size.height; k++, bf += avgStep, bu += faceStep )
170	for( j = 0; j < size.width; j++ )
171	bf[j] += bu[j];
172	}
173
174	bf = avg;
175	for( i = 0; i < size.height; i++, bf += avgStep)
176	#pragma omp parallel for
177	for( j = 0; j < size.width; j++ ) {
178	bf[j] *= m;
179	}
180	}
181
182	void calcCovarMatrix(int nFaces, uchar** faceArr, int faceStep, CvSize size,
183	float* avg, int avgStep, float *covarMatrix)
184	{
185	int i, j;
186
187	for( i = 0; i < nFaces; i++ )
188	{
189	uchar *bu = faceArr[i];
190
191	for( j = i; j < nFaces; j++ )
192	{
193	int k, l;
194	float w = 0.f;
195	float *a = avg;
196	uchar *bu1 = bu;
197	uchar *bu2 = faceArr[j];
198
199	#pragma omp parallel for
200	for( k = 0; k < size.height;
201	k++, bu1 += faceStep, bu2 += faceStep, a += avgStep )
202	{
203	for( l = 0; l < size.width - 3; l += 4 )
204	{
205	float f = a[l];
206	uchar u1 = bu1[l];
207	uchar u2 = bu2[l];
208
209	w += (u1 - f) * (u2 - f);
210	f = a[l + 1];
211	u1 = bu1[l + 1];
212	u2 = bu2[l + 1];
213	w += (u1 - f) * (u2 - f);
214	f = a[l + 2];
215	u1 = bu1[l + 2];
216	u2 = bu2[l + 2];
217	w += (u1 - f) * (u2 - f);
218	f = a[l + 3];
219	u1 = bu1[l + 3];
220	u2 = bu2[l + 3];
221	w += (u1 - f) * (u2 - f);
222	}
223	for( ; l < size.width; l++ )
224	{
225	float f = a[l];
226	uchar u1 = bu1[l];
227	uchar u2 = bu2[l];
228
229	w += (u1 - f) * (u2 - f);
230	}
231	}
232
233	covarMatrix[i * nFaces + j] = covarMatrix[j * nFaces + i] = w;
234	}
235	}
236	}
237
238	void calcEigenFaces(int nFaces, IplImage facesArr, IplImage eigArr, int iter,
239	IplImage avg, float eigVals)
240	{
241	int i,j,k,l, p;
242	float covarMat, ev;
243
244	float *avg_data;
245	int avg_step = 0, eig_step = 0;
246	CvSize size;
247
248	cvGetImageRawData( avg, (uchar **) & avg_data, &avg_step, &size );
249
250	avg_step = avg_step/4;
251	uchar faces = (uchar ) cvAlloc( sizeof( uchar * ) * nFaces );
252	float eigs = (float ) cvAlloc( sizeof( float * ) * iter );
253	int face_step = 0;
254
255	for( i = 0; i < nFaces; i++ )
256	{
257	IplImage *face = facesArr[i];
258	uchar *face_data;
259
260	cvGetImageRawData( face, (uchar **) &face_data, &face_step, NULL);
261	faces[i] = face_data;
262	}
263
264
265	for( i = 0; i < iter; i++ )
266	{
267	IplImage *eig = eigArr[i];
268	float *eig_data;
269
270	cvGetImageRawData( eig, (uchar **) & eig_data, NULL, NULL);
271	eigs[i] = eig_data;
272	}
273
274	calcMeanImage( nFaces,
275	faces,
276	face_step,
277	size,
278	avg_data,
279	avg_step );
280
281	covarMat = (float ) cvAlloc( sizeof( float ) nFaces * nFaces );
282
283	//~ calcCovarMatrix( nFaces,
284	//~ faces,
285	//~ avg_step,
286	//~ size,
287	//~ avg_data,
288	//~ avg_step,
289	//~ covarMat );
290
291	//~ for ( i = 0; i < nFaces; i++ )
292	//~ {
293	//~ for ( j = 0; j < nFaces; j++ )
294	//~ {
295	//~ printf("%f ", covarMat[i*nFaces+j]);
296	//~ }
297	//~ printf("\n");
298	//~ }
299
300	cvCalcCovarMatrixEx( nFaces, facesArr, 0, 0, NULL, NULL, avg, covarMat );
301
302	printf("\n");
303	for ( i = 0; i < nFaces; i++ )
304	{
305	for ( j = 0; j < nFaces; j++ )
306	{
307	printf("%f ", covarMat[i*nFaces+j]);
308	}
309	printf("\n");
310	}
311	printf("\n");
312
313	ev = (float ) cvAlloc( sizeof( float ) nFaces * nFaces );
314
315	int nrot=0;
316	jacobi(covarMat, nFaces, eigVals, ev, &nrot);
317	eigsrt(eigVals, ev, nFaces);
318	//JacobiEigens_32f(covarMat, ev, eigVals, nFaces, 0);
319
320	for ( j = 0; j < nFaces; j++ )
321	{
322	printf("%f ", eigVals[j]);
323	}
324	printf("\n\n");
325
326	for ( i = 0; i < nFaces; i++ )
327	{
328	for ( j = 0; j < nFaces; j++ )
329	{
330	printf("%f ", ev[i*nFaces+j]);
331	}
332	printf("\n");
333	}
334
335	#pragma omp parallel for
336	for( i = 0; i < iter; i++ )
337	eigVals[i] = (float) (1.0 / sqrt( (double)eigVals[i] ));
338
339	for( i = 0; i < iter; i++ )
340	{
341	float *be = eigs[i];
342
343	#pragma omp parallel for
344	for( k = 0; k < size.height; k++, be += avg_step )
345	for( l = 0; l < size.width; l++ )
346	be[l] = 0.0f;
347	}
348
349	for( k = 0; k < nFaces; k++ )
350	{
351	uchar *bv = faces[k];
352
353	for( i = 0; i < iter; i++ )
354	{
355	float v = eigVals[i] * ev[k * nFaces + i];
356	// float v = ev[i * nFaces + k];
357	float *be = eigs[i];
358	uchar *bu = bv;
359
360	float *bf = avg_data;
361
362	#pragma omp parallel for
363	for( p = 0; p < size.height; p++, bu += face_step,
364	bf += avg_step, be += avg_step )
365	{
366	for( l = 0; l < size.width - 3; l += 4 )
367	{
368	float f = bf[l];
369	uchar u = bu[l];
370
371	be[l] += v * (u - f);
372	f = bf[l + 1];
373	u = bu[l + 1];
374	be[l + 1] += v * (u - f);
375	f = bf[l + 2];
376	u = bu[l + 2];
377	be[l + 2] += v * (u - f);
378	f = bf[l + 3];
379	u = bu[l + 3];
380	be[l + 3] += v * (u - f);
381	}
382	for( ; l < size.width; l++ )
383	be[l] += v * (bu[l] - bf[l]);
384	}
385	}
386	}
387
388	#pragma omp parallel for
389	for( i = 0; i < iter; i++ )
390	eigVals[i] = 1.f / (eigVals[i] * eigVals[i]);
391	}
392
393	void calcDecomp( IplImage* face, int nEigens, IplImage** eigArr, IplImage avg, float coeffs )
394	{
395	int i, j, k;
396	float w = 0.0f;
397
398	float *avg_data;
399	uchar *face_data;
400	int avg_step = 0, face_step = 0;
401	CvSize size;
402
403	cvGetImageRawData( avg, (uchar **) & avg_data, &avg_step, &size );
404	cvGetImageRawData( face, &face_data, &face_step, NULL);
405	avg_step = avg_step/4;
406
407	float eigs = (float ) cvAlloc( sizeof( float * ) * nEigens );
408	for( i = 0; i < nEigens; i++ )
409	{
410	IplImage *eig = eigArr[i];
411	float *eig_data;
412
413	cvGetImageRawData( eig, (uchar **) & eig_data, NULL, NULL );
414	eigs[i] = eig_data;
415	}
416
417	for( k = 0; k < nEigens; k++ )
418	{
419	float *be = eigs[k];
420	uchar *bu = face_data;
421	float *bf = avg_data;
422
423	for( i = 0; i < size.height; i++, bu+= face_step,
424	be += avg_step, bf += avg_step )
425	{
426	for( j = 0; j < size.width - 4; j += 4 )
427	{
428	float o = (float) bu[j];
429	float e = be[j];
430	float a = bf[j];
431
432	w += e * (o - a);
433	o = (float) bu[j + 1];
434	e = be[j + 1];
435	a = bf[j + 1];
436	w += e * (o - a);
437	o = (float) bu[j + 2];
438	e = be[j + 2];
439	a = bf[j + 2];
440	w += e * (o - a);
441	o = (float) bu[j + 3];
442	e = be[j + 3];
443	a = bf[j + 3];
444	w += e * (o - a);
445	}
446	for( ; j < size.width; j++ )
447	w += be[j] * ((float) bu[j] - bf[j]);
448	}
449
450	//~ if( w < -1.0e29f )
451	//~ return CV_NOTDEFINED_ERR;
452	coeffs[i] = w;
453	}
454	}

Note: See TracBrowser for help on using the repository browser.

Download in other formats: