Annotation of /branches/pure-cfg/src/lib/common/eigen3x3.c

Parent Directory | Revision Log

---

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

1 :	jhr	1172	/*! \file eigen.c
2 :			*
3 :			* \author Gordon Kindlmann
4 :			*/
5 :
6 :	jhr	1106	/*
7 :	jhr	1172	* COPYRIGHT (c) 2011 The Diderot Project (http://diderot-language.cs.uchicago.edu)
8 :			* All rights reserved.
9 :			*/
10 :
11 :			/*
12 :	jhr	1106	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 :	jhr	1172	#include <Diderot/diderot.h>
35 :	jhr	1106
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 :	jhr	1174	// from http://en.wikipedia.org/wiki/Square_root_of_3
41 :			#define M_SQRT3 1.732050807568877293527446341506
42 :	jhr	1106
43 :	jhr	1511	#define SUB(v,i) (((Diderot_union3_t)(v)).r[i])
44 :	jhr	1174
45 :	glk	1534	// 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 :	lamonts	1513	#endif
55 :	jhr	1544
56 :	jhr	1106	/*
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 :	jhr	1174	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 :	jhr	1106
69 :	jhr	1593	ret = cross3(r0, r1);
70 :			crs = cross3(r1, r2);
71 :	jhr	1106	/* ret += crs or ret -= crs; whichever makes ret longer */
72 :	jhr	1593	if (dot3(ret, crs) > 0.0)
73 :	jhr	1174	ret += crs;
74 :			else
75 :			ret -= crs;
76 :
77 :	jhr	1593	crs = cross3(r0, r2);
78 :	jhr	1106	/* ret += crs or ret -= crs; whichever makes ret longer */
79 :	jhr	1593	if (dot3(ret, crs) > 0.0)
80 :	jhr	1174	ret += crs;
81 :			else
82 :			ret -= crs;
83 :	jhr	1106
84 :	jhr	1174	return ret;
85 :	jhr	1106	}
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 :	jhr	1544	static void nullspace2 (Diderot_vec3_t rets[2],
92 :	jhr	1174	const Diderot_vec3_t r0,
93 :			const Diderot_vec3_t r1,
94 :			const Diderot_vec3_t r2)
95 :			{
96 :	jhr	1511	Diderot_vec3_t sqr, sum;
97 :			int idx;
98 :	jhr	1106
99 :	jhr	1544	sum = r0;
100 :	jhr	1593	if (dot3(sum, r1) > 0)
101 :	jhr	1184	sum += r1;
102 :			else
103 :			sum -= r1;
104 :	jhr	1593	if (dot3(sum, r2) > 0)
105 :	jhr	1184	sum += r2;
106 :			else
107 :			sum -= r2;
108 :	lamonts	1513	// find largest component, to get most stable expression for a perpendicular vector
109 :	jhr	1184	sqr = sum*sum;
110 :			idx = 0;
111 :	jhr	1511	if (SUB(sqr,0) < SUB(sqr,1))
112 :	jhr	1184	idx = 1;
113 :	jhr	1512	if (SUB(sqr, idx) < SUB(sqr,2))
114 :	jhr	1184	idx = 2;
115 :	lamonts	1513
116 :	jhr	1184	if (0 == idx) {
117 :	jhr	1593	rets[0] = vec3(SUB(sum,1) - SUB(sum,2), -SUB(sum,0), SUB(sum,0));
118 :	jhr	1184	} else if (1 == idx) {
119 :	jhr	1593	rets[0] = vec3(-SUB(sum,1), SUB(sum,0) - SUB(sum,2), SUB(sum,1));
120 :	jhr	1184	} else {
121 :	jhr	1593	rets[0] = vec3(-SUB(sum,2), SUB(sum,2), SUB(sum,0) - SUB(sum,1));
122 :	jhr	1184	}
123 :	lamonts	1513
124 :	jhr	1593	rets[1] = cross3(rets[0], sum);
125 :	jhr	1184	return;
126 :	jhr	1106	}
127 :
128 :			/*
129 :			** Eigensolver for symmetric 3x3 matrix:
130 :			**
131 :			** M00 M01 M02
132 :			** M01 M11 M12
133 :			** M02 M12 M22
134 :			**
135 :	jhr	1174	** Must be passed pointer to eval vector, and will compute eigenvalues
136 :	jhr	1106	** 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 :	jhr	1544	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 :	jhr	1174	{
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 :	jhr	1106
160 :	jhr	1174	/* copy the given matrix elements */
161 :			M00 = _M00;
162 :			M01 = _M01;
163 :			M02 = _M02;
164 :			M11 = _M11;
165 :			M12 = _M12;
166 :			M22 = _M22;
167 :	jhr	1106
168 :	jhr	1174	/*
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 :	jhr	1106
193 :	jhr	1174	/* 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 :			double mm, ss, cc;
204 :			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 :			double U;
216 :			/* 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 :	jhr	1106	}
248 :
249 :	jhr	1544	int Diderot_evecs3x3 (
250 :	jhr	1174	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 :	jhr	1106
257 :	jhr	1509	/* 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 :	jhr	1106
262 :	jhr	1509	/* 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 :			double mm, ss, cc;
304 :			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 :			double U;
316 :			/* 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 :	jhr	1593	Diderot_vec3_t ev = vec3(eval[0], eval[1], eval[2]);
337 :	jhr	1174	if (ROOT_THREE == roots) {
338 :	jhr	1509	evecs[0] = nullspace1 (
339 :	jhr	1593	vec3(M00 - eval[0], M01, M02),
340 :			vec3(M01, M11 - eval[0], M12),
341 :			vec3(M02, M12, M22 - eval[0]));
342 :	jhr	1509	evecs[1] = nullspace1 (
343 :	jhr	1593	vec3(M00 - eval[1], M01, M02),
344 :			vec3(M01, M11 - eval[1], M12),
345 :			vec3(M02, M12, M22 - eval[1]));
346 :	jhr	1509	evecs[2] = nullspace1 (
347 :	jhr	1593	vec3(M00 - eval[2], M01, M02),
348 :			vec3(M01, M11 - eval[2], M12),
349 :			vec3(M02, M12, M22 - eval[2]));
350 :	jhr	1106	}
351 :	jhr	1174	else if (ROOT_SINGLE_DOUBLE == roots) {
352 :			if (eval[1] == eval[2]) {
353 :			/* one big (eval[0]) , two small (eval[1,2]) */
354 :	jhr	1509	evecs[0] = nullspace1 (
355 :	jhr	1593	vec3(M00 - eval[0], M01, M02),
356 :			vec3(M01, M11 - eval[0], M12),
357 :			vec3(M02, M12, M22 - eval[0]));
358 :	jhr	1509	nullspace2 (evecs+1,
359 :	jhr	1593	vec3(M00 - eval[1], M01, M02),
360 :			vec3(M01, M11 - eval[1], M12),
361 :			vec3(M02, M12, M22 - eval[1]));
362 :	jhr	1174	}
363 :			else {
364 :			/* two big (eval[0,1]), one small (eval[2]) */
365 :	jhr	1509	nullspace2 (evecs,
366 :	jhr	1593	vec3(M00 - eval[0], M01, M02),
367 :			vec3(M01, M11 - eval[0], M12),
368 :			vec3(M02, M12, M22 - eval[0]));
369 :	jhr	1509	evecs[2] = nullspace1 (
370 :	jhr	1593	vec3(M00 - eval[2], M01, M02),
371 :			vec3(M01, M11 - eval[2], M12),
372 :			vec3(M02, M12, M22 - eval[2]));
373 :	jhr	1174	}
374 :			}
375 :	jhr	1106	else {
376 :	jhr	1174	/* ROOT_TRIPLE == roots; use any basis for eigenvectors */
377 :	jhr	1593	evecs[0] = vec3(1, 0, 0);
378 :			evecs[1] = vec3(0, 1, 0);
379 :			evecs[2] = vec3(0, 0, 1);
380 :	jhr	1106	}
381 :	jhr	1184	/* we always make sure it's really orthonormal; keeping fixed the
382 :	jhr	1174	* eigenvector associated with the largest-magnitude eigenvalue
383 :			*/
384 :			if (ABS(eval[0]) > ABS(eval[2])) {
385 :	jhr	1184	/* normalize evec[0] but don't move it */
386 :	jhr	1593	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 :	jhr	1174	}
391 :	jhr	1184	else {
392 :			/* normalize evec[2] but don't move it */
393 :	jhr	1593	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 :	jhr	1184	}
398 :	glk	1559	/* 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 :	jhr	1106
404 :	jhr	1509	/* BEGIN #include "teigen-evals-B.c" */
405 :	jhr	1174	/* multiply back by eigenvalue L2 norm */
406 :			eval[0] /= rnorm;
407 :			eval[1] /= rnorm;
408 :			eval[2] /= rnorm;
409 :	jhr	1106
410 :	jhr	1174	/* add back in the eigenvalue mean */
411 :			eval[0] += mean;
412 :			eval[1] += mean;
413 :			eval[2] += mean;
414 :	jhr	1509	/* END #include "teigen-evals-B.c" */
415 :	jhr	1174
416 :			return roots;
417 :	jhr	1106	}
418 :	glk	1559
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