Context Navigation

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

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

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

Download in other formats: