|
/Users/lyon/j4p/src/face/EigenFaceComputation.java
|
1 package face;
2
3
4 import math.linearAlgebra.EigenvalueDecomposition;
5 import math.linearAlgebra.Matrix;
6
7 /**
8 * Computes an "face space" used for face recognition.
9 *
10 * This idea/algorhitm was derieved from Matthew A. Turn and Alex P. Pentland
11 * paper, titled: "Face Recognition using Eigenfaces"
12 * (<a href="<a href="http://www.cs.ucsb.edu/~mturk/Papers/mturk-CVPR91.pdf">http://www.cs.ucsb.edu/~mturk/Papers/mturk-CVPR91.pdf</a>"><a href="http://www.cs.ucsb.edu/~mturk/Papers/mturk-CVPR91.pdf">http://www.cs.ucsb.edu/~mturk/Papers/mturk-CVPR91.pdf</a></a>)
13 *<br><br>;
14 * The way this algorhitm works is by treating face recognition as a
15 * "two-dimensional recognition problem, taking advantage of the fact that
16 * faces are normally upright and thus may be described by a small set
17 * of 2-D characterisits views. Face images are projected onto a
18 * feature space ('face space') that best encodes the variation
19 * among known face images. The face space is defined by the
20 * "eigenfaces", which are the eigenvectors of the set of faces;
21 * they do not necessarily correspond to isolated features such as eyes,
22 * ears, and noses.
23 *<br><br>
24 * This work is released to the public with no license and no warranties.
25 * Do with the code as you wish.
26 * <br><b>NOTE</b>: This package uses linearAlgebra for computing eigenvalues and eigenvectors.
27 *
28 * Which is available at:
29 * <a href="<a href="http://math.nist.gov/javanumerics/jama/">http://math.nist.gov/javanumerics/jama/</a>"><a href="http://math.nist.gov/javanumerics/jama/">http://math.nist.gov/javanumerics/jama/</a></a><br><br>
30 * And in case you are worreid about its copyright states:<br><br>
31 *
32 * This software is a cooperative product of The MathWorks and the National Institute of Standards
33 * and Technology (NIST) which has been released to the public domain. Neither The MathWorks nor NIST assumes any
34 * responsibility whatsoever for its use by other parties, and makes no guarantees, expressed or implied, about its quality,
35 * reliability, or any other characteristic.
36 *
37 *
38 *
39 *<br><br>
40 * @author Konrad Rzeszutek
41 * @version 1.0
42 */
43 public class EigenFaceComputation {
44
45 /**
46 * Our "heuresticly" determined value of how many faces we want to compute.
47 * Just picked the value at random. You can fool around with it.
48 */
49 private final static int MAGIC_NR = 11;
50
51 /**
52 * Compute the "face space" used for face recognition. The recognition is
53 * actually being carried out in the FaceBundle object, but the preparation
54 * of such object requires to do lots of computations. The steps are:
55 * <ol>
56 * <li> Compute an average face.
57 * <li> Build an covariance matrix.
58 * <li> Compute eigenvalues and eigenvector
59 * <li> Select only { MAGIC_NR} largest eigenvalues (and its corresponding eigenvectors)
60 * <li> Compute the faces using our eigenvectors
61 * <li> Compute our eigenspace for our given images.
62 * </ol>
63 * From then the rest of the algorithm (trying to match a face) has to be
64 * called in {@link face.FaceBundle}.
65 *
66 * @param face_v 2-D array. Has to have 16 rows. Each column has
67 * to have the same length. Each
68 * row contains the image in a vector representation.
69 * @param width The width of the image in the row-vector in face_v.
70 * @param height The height of the image in the row-vector in face_v.
71 * @param id The string representing each of the sixteen images.
72 *
73 * @return An FaceBundle usable for recognition.
74 *
75 */
76 public static FaceBundle submit(double[][] face_v, int width, int height, String[] id) {
77
78 int length = width * height;
79 int nrfaces = face_v.length;
80 int i, j, col,rows, pix, image;
81 double temp = 0.0;
82 double[][] faces = new double[nrfaces][length];
83
84 //TestPPM simple = new TestPPM();
85 //simple.setImage(face_v[0],width,height);
86
87 double[] avgF = new double[length];
88
89 /*
90 Compute average face of all of the faces. 1xN^2
91 */
92 for (pix = 0; pix < length; pix++) {
93 temp = 0;
94 for (image = 0; image < nrfaces; image++) {
95 temp += face_v[image][pix];
96 }
97 avgF[pix] = temp / nrfaces;
98 }
99
100 //simple.setImage(avgF, width,height);
101
102 /*
103 Compute difference.
104 */
105
106 for (image = 0; image < nrfaces; image++) {
107
108 for (pix = 0; pix < length; pix++) {
109 face_v[image][pix] = face_v[image][pix] - avgF[pix];
110 }
111 }
112 /* Copy our face vector (MxN^2). We will use it later */
113
114 //for (image = 0; image < nrfaces; image++)
115 // System.arraycopy(face_v[image],0,faces[image],0,length);
116 System.arraycopy(face_v, 0, faces, 0, face_v.length);
117
118 //simple.setImage(face_v[0],width,height);
119
120 /*
121 Build covariance matrix. MxM
122 */
123
124 Matrix faceM = new Matrix(face_v, nrfaces, length);
125 Matrix faceM_transpose = faceM.transpose();
126
127 /*
128 Covariance matrix - its MxM (nrfaces x nrfaces)
129 */
130 Matrix covarM = faceM.times(faceM_transpose);
131
132 //double[][] z = covarM.getArray();
133 //System.out.println("Covariance matrix is "
134 // +z.length+" x "+z[0].length);
135
136 /*
137 Compute eigenvalues and eigenvector. Both are MxM
138 */
139 EigenvalueDecomposition E = covarM.eig();
140
141 double[] eigValue = diag(E.getD().getArray());
142 double[][] eigVector = E.getV().getArray();
143
144 /*
145 * We only need the largest associated values of the eigenvalues.
146 * Thus we sort them (and keep an index of them)
147 */
148 int[] index = new int[nrfaces];
149 double[][] tempVector = new double[nrfaces][nrfaces]; /* Temporary new eigVector */
150
151 for (i = 0; i < nrfaces; i++) /* Enumerate all the entries */
152 index[i] = i;
153
154 doubleQuickSort(eigValue, index, 0, nrfaces - 1);
155
156 // Put the index in inverse
157 int[] tempV = new int[nrfaces];
158 for (j = 0; j < nrfaces; j++)
159 tempV[nrfaces - 1 - j] = index[j];
160 /*
161 for (int j = 0; j< nrfaces; j++) {
162 System.out.println(temp[j]+" (was: "+index[j]+") "+eigValue[temp[j]]);
163 }
164 */
165 index = tempV;
166
167 /*
168 * Put the sorted eigenvalues in the appropiate columns.
169 */
170 for (col = nrfaces - 1; col >= 0; col--) {
171 for (rows = 0; rows < nrfaces; rows++) {
172 tempVector[rows][col] = eigVector[rows][index[col]];
173 }
174 }
175 eigVector = tempVector;
176 tempVector = null;
177 eigValue = null;
178 /*
179 * From this point on we don't care about our eigenvalues anymore.
180 *
181 *
182 * We multiply our faceM (MxN^2) with our new eigenvector (MxM),
183 * getting eigenfaces (MxN^2)
184 */
185 Matrix eigVectorM = new Matrix(eigVector, nrfaces, nrfaces);
186 eigVector = eigVectorM.times(faceM).getArray();
187
188
189 /* Normalize our eigen vector matrix. */
190
191 for (image = 0; image < nrfaces; image++) {
192 temp = max(eigVector[image]); // Our max
193 for (pix = 0; pix < eigVector[0].length; pix++)
194 // Normalize
195 eigVector[image][pix] = Math.abs(eigVector[image][pix] / temp);
196 }
197
198 /*
199 * And now the computation of wk (face space), which is a vector and
200 * is of some heuristically determind length. How does 11 sound?
201 *
202 * This is where we are using our copied faces vector (look at the
203 * the beginning of the method) - the reason is b/c the faceM made the
204 * array internally an covariance matrix.
205 */
206
207 double[][] wk = new double[nrfaces][MAGIC_NR]; // M rows, 11 columns
208
209 /*
210 Compute our wk.
211 */
212
213 for (image = 0; image < nrfaces; image++) {
214 for (j = 0; j < MAGIC_NR; j++) {
215 temp = 0.0;
216 for (pix = 0; pix < length; pix++)
217 temp += eigVector[j][pix] * faces[image][pix];
218 wk[image][j] = Math.abs(temp);
219 }
220 }
221
222 /*
223 That's it!
224 */
225
226 FaceBundle b = new FaceBundle(avgF, wk, eigVector, id);
227
228 /*
229 //This is what you would use to recognize a face ...
230
231 double[] inputFace = new double[length];
232 // So we are trying to recognize the 14th image..
233 System.arraycopy(faces[14],0,inputFace,0,length);
234
235 // This is done for virgin images, not ones that we already subtracted.
236 // for ( pix = 0; pix < inputFace.length; pix++) {
237 // inputFace[pix] = inputFace[pix] - avgF[pix];
238 //}
239 */
240 /*
241 double[] input_wk = new double[MAGIC_NR];
242 for (j = 0; j < MAGIC_NR; j++) {
243 temp = 0.0;
244 for ( pix=0; pix <length; pix++)
245 temp += eigVector[j][pix] * inputFace[pix];
246
247 input_wk[j] = Math.abs( temp );
248 }
249 */
250 /*
251 * Find the minimun distance from the input_wk as compared to wk
252 */
253 /*
254 int idx = 0;
255 double[] distance = new double[MAGIC_NR];
256 double[] minDistance = new double[MAGIC_NR];
257
258 for (image = 0; image < nrfaces; image++) {
259 for (j = 0; j < MAGIC_NR; j++)
260 distance[j] = Math.abs(input_wk[j] - wk[image][j]);
261 if (image == 0)
262 System.arraycopy(distance,0,minDistance,0,MAGIC_NR);
263 if (sum(minDistance) > sum(distance)) {
264 idx = image;
265 System.arraycopy(distance,0,minDistance,0,MAGIC_NR);
266
267 }
268 }
269 */
270
271 /*
272 * Normalize our minimum distance.
273 */
274 /*
275 divide(minDistance, max(minDistance)+1);
276 double minD = sum(minDistance);
277 System.out.println("image is idx "+idx+" distance from face: "+minD);
278
279
280 //TestPPM simple = new TestPPM();
281 int[] bb = new int[length];
282
283 temp = max(eigVector[0]);
284
285 for ( i = 0; i < width*height; i++) {
286 bb[i] = (int) (255*(1 - eigVector[0][i] / temp ));
287 }
288
289 simple.setImage(bb,width,height);
290 */
291 return b;
292 }
293
294 /**
295 * Find the diagonal of an matrix.
296 *
297 * @param m the number of rows and columns must be the same
298 * @return the diagonal of the matrix
299 */
300 static double[] diag(double[][] m) {
301
302 double[] d = new double[m.length];
303 for (int i = 0; i < m.length; i++)
304 d[i] = m[i][i];
305 return d;
306 }
307
308 /**
309 * Divide each element in <code>v</code> by <code>b</code>
310 * No checking for division by zero.
311 *
312 * @param v vector containing numbers.
313 * @param b scalar used to divied each element in the v vector
314 *
315 * a vector having each element divided by <code>b</code> scalar.
316 *
317 */
318 static void divide(double[] v, double b) {
319
320 for (int i = 0; i < v.length; i++)
321 v[i] = v[i] / b;
322
323
324 }
325
326 /**
327 * The sum of the vector.
328 *
329 * @param a vector with numbers
330 * @return a scalar with the sum of each element in the <code>a</code> vector
331 */
332 static double sum(double[] a) {
333
334 double b = a[0];
335 for (int i = 0; i < a.length; i++)
336 b += a[i];
337
338 return b;
339
340 }
341
342 /**
343 * The max of the vector a.
344 *
345 * @param a the vector
346 *
347 * @return the sum of all the elements on <code>a</code>
348 */
349 static double max(double[] a) {
350 double b = a[0];
351 for (int i = 0; i < a.length; i++)
352 if (a[i] > b) b = a[i];
353
354 return b;
355 }
356
357 /**
358 * Quick sort on a vector with an index.
359 *
360 * @param a the array of numbers. This will be modified and sorted
361 * ascendingly (smalles to highest)
362 * @param index the index of the numbers as related to original
363 * location.
364 * @param lo0 the index where to start from. Usually 0.
365 * @param hi0 the index where to stop. Usually a.length()
366 */
367 static void doubleQuickSort(double a[], int index[], int lo0, int hi0) {
368 int lo = lo0;
369 int hi = hi0;
370 double mid;
371
372 if (hi0 > lo0) {
373
374 /* Arbitrarily establishing partition element as the midpoint of
375 * the array.
376 */
377 mid = a[(lo0 + hi0) / 2];
378 // loop through the array until indices cross
379 while (lo <= hi) {
380 /* find the first element that is greater than or equal to
381 * the partition element starting from the left Index.
382 */
383 while ((lo < hi0) && (a[lo] < mid)) {
384 ++lo;
385 }
386
387 /* find an element that is smaller than or equal to
388 * the partition element starting from the right Index.
389 */
390 while ((hi > lo0) && (a[hi] > mid)) {
391 --hi;
392 }
393
394 // if the indexes have not crossed, swap
395 if (lo <= hi) {
396 swap(a, index, lo, hi);
397 ++lo;
398 --hi;
399 }
400 }
401 /* If the right index has not reached the left side of array
402 * must now sort the left partition.
403 */
404 if (lo0 < hi) {
405 doubleQuickSort(a, index, lo0, hi);
406 }
407 /* If the left index has not reached the right side of array
408 * must now sort the right partition.
409 */
410 if (lo < hi0) {
411 doubleQuickSort(a, index, lo, hi0);
412 }
413 }
414 }
415
416 static private void swap(double a[], int[] index, int i, int j) {
417 double T;
418 T = a[i];
419 a[i] = a[j];
420 a[j] = T;
421 // Index
422 index[i] = i;
423 index[j] = j;
424 }
425 }