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

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