Home My Page Projects Code Snippets Project Openings diderot
Summary Activity Tracker Tasks SCM

SCM Repository

[diderot] Annotation of /trunk/src/lib/common/eigen3x3.c
ViewVC logotype

Annotation of /trunk/src/lib/common/eigen3x3.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1679 - (view) (download) (as text)

1 : jhr 1640 /*! \file eigen.c
2 :     *
3 :     * \author Gordon Kindlmann
4 :     */
5 :    
6 :     /*
7 :     * COPYRIGHT (c) 2011 The Diderot Project (http://diderot-language.cs.uchicago.edu)
8 :     * All rights reserved.
9 :     */
10 :    
11 :     /*
12 :     Teem: Tools to process and visualize scientific data and images
13 :     Copyright (C) 2011, 2010, 2009, University of Chicago
14 :     Copyright (C) 2008, 2007, 2006, 2005 Gordon Kindlmann
15 :     Copyright (C) 2004, 2003, 2002, 2001, 2000, 1999, 1998 University of Utah
16 :    
17 :     This library is free software; you can redistribute it and/or
18 :     modify it under the terms of the GNU Lesser General Public License
19 :     (LGPL) as published by the Free Software Foundation; either
20 :     version 2.1 of the License, or (at your option) any later version.
21 :     The terms of redistributing and/or modifying this software also
22 :     include exceptions to the LGPL that facilitate static linking.
23 :    
24 :     This library is distributed in the hope that it will be useful,
25 :     but WITHOUT ANY WARRANTY; without even the implied warranty of
26 :     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
27 :     Lesser General Public License for more details.
28 :    
29 :     You should have received a copy of the GNU Lesser General Public License
30 :     along with this library; if not, write to Free Software Foundation, Inc.,
31 :     51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
32 :     */
33 :    
34 :     #include <Diderot/diderot.h>
35 :    
36 :     #define ROOT_TRIPLE 2 /* ell_cubic_root_triple */
37 :     #define ROOT_SINGLE_DOUBLE 3 /* ell_cubic_root_single_double */
38 :     #define ROOT_THREE 4 /* ell_cubic_root_three */
39 :    
40 :     // from http://en.wikipedia.org/wiki/Square_root_of_3
41 :     #define M_SQRT3 1.732050807568877293527446341506
42 :    
43 :     #define SUB(v,i) (((Diderot_union3_t)(v)).r[i])
44 :    
45 :     // OpenCL 1.0 does not specify the C representation of the host types
46 :     // for OpenCL vector types (e.g., cl_float4), so we have to handle this
47 :     // mechanism with a macro.
48 :     #if defined(CL_HOST_VECTORS_ARE_STRUCTS)
49 :     # define VSUB(v,i) (v).s[i]
50 :     #elif defined(CL_HOST_VECTORS_ARE_ARRAYS)
51 :     # define VSUB(v,i) (v)[i]
52 :     #else
53 :     # error unable to access elements of host vectors
54 :     #endif
55 :    
56 :     /*
57 :     ** All the three given vectors span only a 2D space, and this finds
58 :     ** the normal to that plane. Simply sums up all the pair-wise
59 :     ** cross-products to get a good estimate. Trick is getting the cross
60 :     ** products to line up before summing.
61 :     */
62 :     STATIC_INLINE Diderot_vec3_t nullspace1 (
63 :     const Diderot_vec3_t r0,
64 :     const Diderot_vec3_t r1,
65 :     const Diderot_vec3_t r2)
66 :     {
67 :     Diderot_vec3_t ret, crs;
68 :    
69 :     ret = cross3(r0, r1);
70 :     crs = cross3(r1, r2);
71 :     /* ret += crs or ret -= crs; whichever makes ret longer */
72 :     if (dot3(ret, crs) > 0.0)
73 :     ret += crs;
74 :     else
75 :     ret -= crs;
76 :    
77 :     crs = cross3(r0, r2);
78 :     /* ret += crs or ret -= crs; whichever makes ret longer */
79 :     if (dot3(ret, crs) > 0.0)
80 :     ret += crs;
81 :     else
82 :     ret -= crs;
83 :    
84 :     return ret;
85 :     }
86 :    
87 :     /*
88 :     ** All vectors are in the same 1D space, we have to find two
89 :     ** mutually vectors perpendicular to that span
90 :     */
91 :     static void nullspace2 (Diderot_vec3_t rets[2],
92 :     const Diderot_vec3_t r0,
93 :     const Diderot_vec3_t r1,
94 :     const Diderot_vec3_t r2)
95 :     {
96 :     Diderot_vec3_t sqr, sum;
97 :     int idx;
98 :    
99 :     sum = r0;
100 :     if (dot3(sum, r1) > 0)
101 :     sum += r1;
102 :     else
103 :     sum -= r1;
104 :     if (dot3(sum, r2) > 0)
105 :     sum += r2;
106 :     else
107 :     sum -= r2;
108 :     // find largest component, to get most stable expression for a perpendicular vector
109 :     sqr = sum*sum;
110 :     idx = 0;
111 :     if (SUB(sqr,0) < SUB(sqr,1))
112 :     idx = 1;
113 :     if (SUB(sqr, idx) < SUB(sqr,2))
114 :     idx = 2;
115 :    
116 :     if (0 == idx) {
117 :     rets[0] = vec3(SUB(sum,1) - SUB(sum,2), -SUB(sum,0), SUB(sum,0));
118 :     } else if (1 == idx) {
119 :     rets[0] = vec3(-SUB(sum,1), SUB(sum,0) - SUB(sum,2), SUB(sum,1));
120 :     } else {
121 :     rets[0] = vec3(-SUB(sum,2), SUB(sum,2), SUB(sum,0) - SUB(sum,1));
122 :     }
123 :    
124 :     rets[1] = cross3(rets[0], sum);
125 :     return;
126 :     }
127 :    
128 :     /*
129 :     ** Eigensolver for symmetric 3x3 matrix:
130 :     **
131 :     ** M00 M01 M02
132 :     ** M01 M11 M12
133 :     ** M02 M12 M22
134 :     **
135 :     ** Must be passed pointer to eval vector, and will compute eigenvalues
136 :     ** only if evec[9] is non-NULL. Computed eigenvectors are at evec+0,
137 :     ** evec+3, and evec+6.
138 :     **
139 :     ** Return value indicates something about the eigenvalue solution to
140 :     ** the cubic characteristic equation; see ROOT_ #defines above
141 :     **
142 :     ** HEY: the numerical precision issues here are very subtle, and
143 :     ** merit some more scrutiny. With evals (1.000001, 1, 1), for example,
144 :     ** whether it comes back as a single/double root, vs three distinct roots,
145 :     ** is determines by the comparison between "D" and "epsilon", and the
146 :     ** setting of epsilon seems pretty arbitrary at this point...
147 :     **
148 :     */
149 :     int Diderot_evals3x3 (
150 :     Diderot_real_t eval[3],
151 :     const Diderot_real_t _M00, const Diderot_real_t _M01, const Diderot_real_t _M02,
152 :     const Diderot_real_t _M11, const Diderot_real_t _M12,
153 :     const Diderot_real_t _M22)
154 :     {
155 :     Diderot_real_t mean, norm, rnorm, Q, R, QQQ, D, theta, M00, M01, M02, M11, M12, M22;
156 :     // FIXME: probably want a different value of epsilon for single-precision
157 :     Diderot_real_t epsilon = 1.0E-12;
158 :     int roots;
159 :    
160 :     /* copy the given matrix elements */
161 :     M00 = _M00;
162 :     M01 = _M01;
163 :     M02 = _M02;
164 :     M11 = _M11;
165 :     M12 = _M12;
166 :     M22 = _M22;
167 :    
168 :     /*
169 :     ** subtract out the eigenvalue mean (will add back to evals later);
170 :     ** helps with numerical stability
171 :     */
172 :     mean = (M00 + M11 + M22)/3.0;
173 :     M00 -= mean;
174 :     M11 -= mean;
175 :     M22 -= mean;
176 :    
177 :     /*
178 :     ** divide out L2 norm of eigenvalues (will multiply back later);
179 :     ** this too seems to help with stability
180 :     */
181 :     norm = SQRT(M00*M00 + 2*M01*M01 + 2*M02*M02 +
182 :     M11*M11 + 2*M12*M12 +
183 :     M22*M22);
184 :     // QUESTION: what if norm is very close to zero?
185 :     rnorm = norm ? 1.0/norm : 1.0;
186 :     M00 *= rnorm;
187 :     M01 *= rnorm;
188 :     M02 *= rnorm;
189 :     M11 *= rnorm;
190 :     M12 *= rnorm;
191 :     M22 *= rnorm;
192 :    
193 :     /* this code is a mix of prior Teem code and ideas from Eberly's
194 :     * "Eigensystems for 3 x 3 Symmetric Matrices (Revisited)"
195 :     */
196 :     Q = (M01*M01 + M02*M02 + M12*M12 - M00*M11 - M00*M22 - M11*M22)/3.0;
197 :     QQQ = Q*Q*Q;
198 :     R = (M00*M11*M22 + M02*(2*M01*M12 - M02*M11)
199 :     - M00*M12*M12 - M01*M01*M22)/2.0;
200 :     D = QQQ - R*R;
201 :     if (D > epsilon) {
202 :     /* three distinct roots- this is the most common case */
203 : jhr 1679 Diderot_real_t mm, ss, cc;
204 : jhr 1640 theta = ATAN2(SQRT(D), R)/3.0;
205 :     mm = SQRT(Q);
206 :     ss = SIN(theta);
207 :     cc = COS(theta);
208 :     eval[0] = 2*mm*cc;
209 :     eval[1] = mm*(-cc + M_SQRT3*ss);
210 :     eval[2] = mm*(-cc - M_SQRT3*ss);
211 :     roots = ROOT_THREE;
212 :     }
213 :     /* else D is near enough to zero */
214 :     else if (R < -epsilon || epsilon < R) {
215 : jhr 1679 Diderot_real_t U;
216 : jhr 1640 /* one double root and one single root */
217 :     U = CBRT(R); /* cube root function */
218 :     if (U > 0) {
219 :     eval[0] = 2*U;
220 :     eval[1] = -U;
221 :     eval[2] = -U;
222 :     }
223 :     else {
224 :     eval[0] = -U;
225 :     eval[1] = -U;
226 :     eval[2] = 2*U;
227 :     }
228 :     roots = ROOT_SINGLE_DOUBLE;
229 :     }
230 :     else {
231 :     /* a triple root! */
232 :     eval[0] = eval[1] = eval[2] = 0.0;
233 :     roots = ROOT_TRIPLE;
234 :     }
235 :    
236 :     /* multiply back by eigenvalue L2 norm */
237 :     eval[0] /= rnorm;
238 :     eval[1] /= rnorm;
239 :     eval[2] /= rnorm;
240 :    
241 :     /* add back in the eigenvalue mean */
242 :     eval[0] += mean;
243 :     eval[1] += mean;
244 :     eval[2] += mean;
245 :    
246 :     return roots;
247 :     }
248 :    
249 :     int Diderot_evecs3x3 (
250 :     Diderot_real_t eval[3], Diderot_vec3_t evecs[3],
251 :     const Diderot_real_t _M00, const Diderot_real_t _M01, const Diderot_real_t _M02,
252 :     const Diderot_real_t _M11, const Diderot_real_t _M12,
253 :     const Diderot_real_t _M22)
254 :     {
255 :     Diderot_real_t len, dot;
256 :    
257 :     /* BEGIN #include "teigen-evals-A.c" */
258 :     Diderot_real_t mean, norm, rnorm, Q, R, QQQ, D, theta, M00, M01, M02, M11, M12, M22;
259 :     Diderot_real_t epsilon = 1.0E-12;
260 :     int roots;
261 :    
262 :     /* copy the given matrix elements */
263 :     M00 = _M00;
264 :     M01 = _M01;
265 :     M02 = _M02;
266 :     M11 = _M11;
267 :     M12 = _M12;
268 :     M22 = _M22;
269 :    
270 :     /*
271 :     ** subtract out the eigenvalue mean (will add back to evals later);
272 :     ** helps with numerical stability
273 :     */
274 :     mean = (M00 + M11 + M22)/3.0;
275 :     M00 -= mean;
276 :     M11 -= mean;
277 :     M22 -= mean;
278 :    
279 :     /*
280 :     ** divide out L2 norm of eigenvalues (will multiply back later);
281 :     ** this too seems to help with stability
282 :     */
283 :     norm = SQRT(M00*M00 + 2*M01*M01 + 2*M02*M02 +
284 :     M11*M11 + 2*M12*M12 +
285 :     M22*M22);
286 :     rnorm = norm ? 1.0/norm : 1.0;
287 :     M00 *= rnorm;
288 :     M01 *= rnorm;
289 :     M02 *= rnorm;
290 :     M11 *= rnorm;
291 :     M12 *= rnorm;
292 :     M22 *= rnorm;
293 :    
294 :     /* this code is a mix of prior Teem code and ideas from Eberly's
295 :     * "Eigensystems for 3 x 3 Symmetric Matrices (Revisited)"
296 :     */
297 :     Q = (M01*M01 + M02*M02 + M12*M12 - M00*M11 - M00*M22 - M11*M22)/3.0;
298 :     QQQ = Q*Q*Q;
299 :     R = (M00*M11*M22 + M02*(2*M01*M12 - M02*M11) - M00*M12*M12 - M01*M01*M22)/2.0;
300 :     D = QQQ - R*R;
301 :     if (D > epsilon) {
302 :     /* three distinct roots- this is the most common case */
303 : jhr 1679 Diderot_real_t mm, ss, cc;
304 : jhr 1640 theta = ATAN2(SQRT(D), R)/3.0;
305 :     mm = SQRT(Q);
306 :     ss = SIN(theta);
307 :     cc = COS(theta);
308 :     eval[0] = 2*mm*cc;
309 :     eval[1] = mm*(-cc + SQRT(3.0)*ss);
310 :     eval[2] = mm*(-cc - SQRT(3.0)*ss);
311 :     roots = ROOT_THREE;
312 :     }
313 :     /* else D is near enough to zero */
314 :     else if (R < -epsilon || epsilon < R) {
315 : jhr 1679 Diderot_real_t U;
316 : jhr 1640 /* one double root and one single root */
317 :     U = CBRT(R); /* cube root function */
318 :     if (U > 0) {
319 :     eval[0] = 2*U;
320 :     eval[1] = -U;
321 :     eval[2] = -U;
322 :     } else {
323 :     eval[0] = -U;
324 :     eval[1] = -U;
325 :     eval[2] = 2*U;
326 :     }
327 :     roots = ROOT_SINGLE_DOUBLE;
328 :     }
329 :     else {
330 :     /* a triple root! */
331 :     eval[0] = eval[1] = eval[2] = 0.0;
332 :     roots = ROOT_TRIPLE;
333 :     }
334 :     /* END #include "teigen-evals-A.c" */
335 :    
336 :     Diderot_vec3_t ev = vec3(eval[0], eval[1], eval[2]);
337 :     if (ROOT_THREE == roots) {
338 :     evecs[0] = nullspace1 (
339 :     vec3(M00 - eval[0], M01, M02),
340 :     vec3(M01, M11 - eval[0], M12),
341 :     vec3(M02, M12, M22 - eval[0]));
342 :     evecs[1] = nullspace1 (
343 :     vec3(M00 - eval[1], M01, M02),
344 :     vec3(M01, M11 - eval[1], M12),
345 :     vec3(M02, M12, M22 - eval[1]));
346 :     evecs[2] = nullspace1 (
347 :     vec3(M00 - eval[2], M01, M02),
348 :     vec3(M01, M11 - eval[2], M12),
349 :     vec3(M02, M12, M22 - eval[2]));
350 :     }
351 :     else if (ROOT_SINGLE_DOUBLE == roots) {
352 :     if (eval[1] == eval[2]) {
353 :     /* one big (eval[0]) , two small (eval[1,2]) */
354 :     evecs[0] = nullspace1 (
355 :     vec3(M00 - eval[0], M01, M02),
356 :     vec3(M01, M11 - eval[0], M12),
357 :     vec3(M02, M12, M22 - eval[0]));
358 :     nullspace2 (evecs+1,
359 :     vec3(M00 - eval[1], M01, M02),
360 :     vec3(M01, M11 - eval[1], M12),
361 :     vec3(M02, M12, M22 - eval[1]));
362 :     }
363 :     else {
364 :     /* two big (eval[0,1]), one small (eval[2]) */
365 :     nullspace2 (evecs,
366 :     vec3(M00 - eval[0], M01, M02),
367 :     vec3(M01, M11 - eval[0], M12),
368 :     vec3(M02, M12, M22 - eval[0]));
369 :     evecs[2] = nullspace1 (
370 :     vec3(M00 - eval[2], M01, M02),
371 :     vec3(M01, M11 - eval[2], M12),
372 :     vec3(M02, M12, M22 - eval[2]));
373 :     }
374 :     }
375 :     else {
376 :     /* ROOT_TRIPLE == roots; use any basis for eigenvectors */
377 :     evecs[0] = vec3(1, 0, 0);
378 :     evecs[1] = vec3(0, 1, 0);
379 :     evecs[2] = vec3(0, 0, 1);
380 :     }
381 :     /* we always make sure it's really orthonormal; keeping fixed the
382 :     * eigenvector associated with the largest-magnitude eigenvalue
383 :     */
384 :     if (ABS(eval[0]) > ABS(eval[2])) {
385 :     /* normalize evec[0] but don't move it */
386 :     evecs[0] = normalize3(evecs[0]);
387 :     evecs[1] -= scale3(dot3(evecs[1], evecs[0]), evecs[0]);
388 :     evecs[1] = normalize3(evecs[1]);
389 :     evecs[2] = cross3(evecs[0], evecs[1]);
390 :     }
391 :     else {
392 :     /* normalize evec[2] but don't move it */
393 :     evecs[2] = normalize3(evecs[2]);
394 :     evecs[1] -= scale3(dot3(evecs[1], evecs[2]), evecs[2]);
395 :     evecs[1] = normalize3(evecs[1]);
396 :     evecs[0] = cross3(evecs[1], evecs[2]);
397 :     }
398 :     /* note that the right-handedness check has been folded into
399 :     the code above to enforce orthogonality. Indeed, some work
400 :     could be removed by never really bothering to find all three
401 :     eigenvectors; just find two and then use the cross-product.
402 :     The next iteration of the code will do that */
403 :    
404 :     /* BEGIN #include "teigen-evals-B.c" */
405 :     /* multiply back by eigenvalue L2 norm */
406 :     eval[0] /= rnorm;
407 :     eval[1] /= rnorm;
408 :     eval[2] /= rnorm;
409 :    
410 :     /* add back in the eigenvalue mean */
411 :     eval[0] += mean;
412 :     eval[1] += mean;
413 :     eval[2] += mean;
414 :     /* END #include "teigen-evals-B.c" */
415 :    
416 :     return roots;
417 :     }
418 :    
419 :    
420 :     #if 0
421 :     /* command-line tester for above*/
422 :    
423 :     #include <math.h>
424 :     #include <teem/hest.h>
425 :     #include <teem/nrrd.h>
426 :     #include <teem/ell.h>
427 :     char *info = ("tests symmetric 3x3 eigensolve.");
428 :    
429 :     int
430 :     main(int argc, const char *argv[0]) {
431 :     const char *me;
432 :     hestOpt *hopt=NULL;
433 :     airArray *mop;
434 :    
435 :     Nrrd *nbhess;
436 :     float *bhess, mat[9];
437 :     Diderot_real_t eval[3];
438 :     Diderot_vec3_t evec[3];
439 :     int roots, ri;
440 :    
441 :     mop = airMopNew();
442 :     me = argv[0];
443 :     hestOptAdd(&hopt, "i", "hess", airTypeOther, 1, 1, &nbhess, NULL,
444 :     "bad hessian",
445 :     NULL, NULL, nrrdHestNrrd);
446 :     hestParseOrDie(hopt, argc-1, argv+1, NULL,
447 :     me, info, AIR_TRUE, AIR_TRUE, AIR_TRUE);
448 :     airMopAdd(mop, hopt, (airMopper)hestOptFree, airMopAlways);
449 :     airMopAdd(mop, hopt, (airMopper)hestParseFree, airMopAlways);
450 :    
451 :     bhess = AIR_CAST(float *, nbhess->data);
452 :     roots = Diderot_evals3x3(eval,
453 :     bhess[0], bhess[1], bhess[2],
454 :     bhess[4], bhess[5],
455 :     bhess[8]);
456 :     printf("%s: roots(%d) = %g %g %g\n", me, roots, eval[0], eval[1], eval[2]);
457 :     ELL_3V_SET(mat + 0, bhess[0], bhess[1], bhess[2]);
458 :     ELL_3V_SET(mat + 3, bhess[3], bhess[4], bhess[5]);
459 :     ELL_3V_SET(mat + 6, bhess[6], bhess[7], bhess[8]);
460 :     printf("%s: det(H) = %g\n", me, ELL_3M_DET(mat));
461 :     for (ri=0; ri<3; ri++) {
462 :     ELL_3V_SET(mat + 0, bhess[0] - eval[ri], bhess[1], bhess[2]);
463 :     ELL_3V_SET(mat + 3, bhess[3], bhess[4] - eval[ri], bhess[5]);
464 :     ELL_3V_SET(mat + 6, bhess[6], bhess[7], bhess[8] - eval[ri]);
465 :     printf("%s: det(H - lI) = %g\n", me, ELL_3M_DET(mat));
466 :     }
467 :    
468 :     roots = Diderot_evecs3x3(eval, evec,
469 :     bhess[0], bhess[1], bhess[2],
470 :     bhess[4], bhess[5],
471 :     bhess[8]);
472 :     printf("%s: roots(%d) = %g %g %g\n", me, roots, eval[0], eval[1], eval[2]);
473 :     printf("%s: evec0 . evec0 = %g\n", me, ELL_3V_DOT(evec[0], evec[0]));
474 :     printf("%s: evec1 . evec1 = %g\n", me, ELL_3V_DOT(evec[1], evec[1]));
475 :     printf("%s: evec2 . evec2 = %g\n", me, ELL_3V_DOT(evec[2], evec[2]));
476 :     printf("%s: evec0 . evec1 = %g\n", me, ELL_3V_DOT(evec[0], evec[1]));
477 :     printf("%s: evec0 . evec2 = %g\n", me, ELL_3V_DOT(evec[0], evec[2]));
478 :     printf("%s: evec1 . evec2 = %g\n", me, ELL_3V_DOT(evec[1], evec[2]));
479 :    
480 :     return 0;
481 :     }
482 :    
483 :     /*
484 :    
485 :     clang -m64 -I/Users/gk/diderot/diderot/branches/pure-cfg/src/include \
486 :     -I../../include -I../../../../../../teem/include/ \
487 :     -L../../../../../../teem/lib \
488 :     -DDIDEROT_TARGET_C -DNDEBUG -DDIDEROT_SINGLE_PRECISION \
489 :     -o eigen3 eigen3x3.c -lteem -lm
490 :    
491 :     */
492 :    
493 :     #endif

root@smlnj-gforge.cs.uchicago.edu
ViewVC Help
Powered by ViewVC 1.0.0