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 3349 - (view) (download) (as text)

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

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