Annotation of /benchmarks/programs/ridge3d/bmark-teem.c

Parent Directory | Revision Log

---

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

1 :	jhr	1584	/*! \file bmark-teem.c
2 :			*
3 :			* \author Gordon Kindlmann
4 :			*/
5 :
6 :	glk	1582	#include <stdio.h>
7 :			#include <math.h>
8 :			#include "teem/nrrd.h"
9 :			#include "teem/gage.h"
10 :			#include "teem/ell.h"
11 :
12 :	jhr	1597	extern int _gageProbe(gageContext *ctx, double _xi, double _yi, double _zi, double _si);
13 :	glk	1591
14 :	jhr	1586	#include "teem-defs.h"
15 :	jhr	1584
16 :	glk	1591	#define ROOT_TRIPLE 2 /* ell_cubic_root_triple */
17 :			#define ROOT_SINGLE_DOUBLE 3 /* ell_cubic_root_single_double */
18 :			#define ROOT_THREE 4 /* ell_cubic_root_three */
19 :
20 :			#define ABS(a) (((a) > 0.0f ? (a) : -(a)))
21 :			#define VEC_SET(v, a, b, c) \
22 :			((v)[0] = (a), (v)[1] = (b), (v)[2] = (c))
23 :			#define VEC_DOT(v1, v2) \
24 :			((v1)[0]*(v2)[0] + (v1)[1]*(v2)[1] + (v1)[2]*(v2)[2])
25 :			#define VEC_CROSS(v3, v1, v2) \
26 :			((v3)[0] = (v1)[1]*(v2)[2] - (v1)[2]*(v2)[1], \
27 :			(v3)[1] = (v1)[2]*(v2)[0] - (v1)[0]*(v2)[2], \
28 :			(v3)[2] = (v1)[0]*(v2)[1] - (v1)[1]*(v2)[0])
29 :			#define VEC_ADD(v1, v2) \
30 :			((v1)[0] += (v2)[0], \
31 :			(v1)[1] += (v2)[1], \
32 :			(v1)[2] += (v2)[2])
33 :			#define VEC_SUB(v1, v2) \
34 :			((v1)[0] -= (v2)[0], \
35 :			(v1)[1] -= (v2)[1], \
36 :			(v1)[2] -= (v2)[2])
37 :			#define VEC_SCL(v1, s) \
38 :			((v1)[0] *= (s), \
39 :			(v1)[1] *= (s), \
40 :			(v1)[2] *= (s))
41 :			#define VEC_LEN(v) (SQRT(VEC_DOT(v,v)))
42 :			#define VEC_NORM(v, len) ((len) = VEC_LEN(v), VEC_SCL(v, 1.0/len))
43 :			#define VEC_SCL_SUB(v1, s, v2) \
44 :			((v1)[0] -= (s)*(v2)[0], \
45 :			(v1)[1] -= (s)*(v2)[1], \
46 :			(v1)[2] -= (s)*(v2)[2])
47 :			#define VEC_COPY(v1, v2) \
48 :			((v1)[0] = (v2)[0], \
49 :			(v1)[1] = (v2)[1], \
50 :			(v1)[2] = (v2)[2])
51 :
52 :			#if 1
53 :			# define REAL float
54 :			# define ATAN2 atan2f
55 :			# define SQRT sqrtf
56 :			# define COS cosf
57 :			# define SIN sinf
58 :			#else
59 :			# define REAL double
60 :			# define ATAN2 atan2
61 :			# define SQRT sqrt
62 :			# define COS cos
63 :			# define SIN sin
64 :			#endif
65 :
66 :			/*
67 :			** All the three given vectors span only a 2D space, and this finds
68 :			** the normal to that plane. Simply sums up all the pair-wise
69 :			** cross-products to get a good estimate. Trick is getting the cross
70 :			** products to line up before summing.
71 :			*/
72 :	jhr	1597	static inline void nullspace1(REAL ret[3], const REAL r0[3], const REAL r1[3], const REAL r2[3])
73 :			{
74 :	glk	1591	REAL crs[3];
75 :
76 :			/* ret = r0 x r1 */
77 :			VEC_CROSS(ret, r0, r1);
78 :			/* crs = r1 x r2 */
79 :			VEC_CROSS(crs, r1, r2);
80 :			/* ret += crs or ret -= crs; whichever makes ret longer */
81 :			if (VEC_DOT(ret, crs) > 0) {
82 :			VEC_ADD(ret, crs);
83 :			} else {
84 :			VEC_SUB(ret, crs);
85 :			}
86 :			/* crs = r0 x r2 */
87 :			VEC_CROSS(crs, r0, r2);
88 :			/* ret += crs or ret -= crs; whichever makes ret longer */
89 :			if (VEC_DOT(ret, crs) > 0) {
90 :			VEC_ADD(ret, crs);
91 :			} else {
92 :			VEC_SUB(ret, crs);
93 :			}
94 :
95 :			return;
96 :	glk	1582	}
97 :
98 :	glk	1591	/*
99 :			** All vectors are in the same 1D space, we have to find two
100 :			** mutually vectors perpendicular to that span
101 :			*/
102 :	jhr	1597	static inline void nullspace2(REAL reta[3], REAL retb[3], const REAL r0[3], const REAL r1[3], const REAL r2[3])
103 :			{
104 :	glk	1591	REAL sqr[3], sum[3];
105 :			int idx;
106 :
107 :			VEC_COPY(sum, r0);
108 :			if (VEC_DOT(sum, r1) > 0) {
109 :			VEC_ADD(sum, r1);
110 :			} else {
111 :			VEC_SUB(sum, r1);
112 :			}
113 :			if (VEC_DOT(sum, r2) > 0) {
114 :			VEC_ADD(sum, r2);
115 :			} else {
116 :			VEC_SUB(sum, r2);
117 :			}
118 :			/* find largest component, to get most stable expression for a
119 :			perpendicular vector */
120 :			sqr[0] = sum[0]*sum[0];
121 :			sqr[1] = sum[1]*sum[1];
122 :			sqr[2] = sum[2]*sum[2];
123 :			idx = 0;
124 :			if (sqr[0] < sqr[1])
125 :			idx = 1;
126 :			if (sqr[idx] < sqr[2])
127 :			idx = 2;
128 :			/* reta will be perpendicular to sum */
129 :			if (0 == idx) {
130 :			VEC_SET(reta, sum[1] - sum[2], -sum[0], sum[0]);
131 :			} else if (1 == idx) {
132 :			VEC_SET(reta, -sum[1], sum[0] - sum[2], sum[1]);
133 :			} else {
134 :			VEC_SET(reta, -sum[2], sum[2], sum[0] - sum[1]);
135 :			}
136 :			/* and now retb will be perpendicular to both reta and sum */
137 :			VEC_CROSS(retb, reta, sum);
138 :			return;
139 :			}
140 :
141 :			int
142 :			evals_evecs(REAL eval[3], REAL evec[9],
143 :			const REAL _M00, const REAL _M01, const REAL _M02,
144 :			const REAL _M11, const REAL _M12,
145 :			const REAL _M22) {
146 :			REAL r0[3], r1[3], r2[3], crs[3], len, dot;
147 :
148 :			REAL mean, norm, rnorm, Q, R, QQQ, D, theta,
149 :			M00, M01, M02, M11, M12, M22;
150 :			REAL epsilon = 1.0E-12;
151 :			int roots;
152 :
153 :			/* copy the given matrix elements */
154 :			M00 = _M00;
155 :			M01 = _M01;
156 :			M02 = _M02;
157 :			M11 = _M11;
158 :			M12 = _M12;
159 :			M22 = _M22;
160 :
161 :			/*
162 :			** subtract out the eigenvalue mean (will add back to evals later);
163 :			** helps with numerical stability
164 :			*/
165 :			mean = (M00 + M11 + M22)/3.0;
166 :			M00 -= mean;
167 :			M11 -= mean;
168 :			M22 -= mean;
169 :
170 :			/*
171 :			** divide out L2 norm of eigenvalues (will multiply back later);
172 :			** this too seems to help with stability
173 :			*/
174 :			norm = SQRT(M00*M00 + 2*M01*M01 + 2*M02*M02 +
175 :			M11*M11 + 2*M12*M12 +
176 :			M22*M22);
177 :			rnorm = norm ? 1.0/norm : 1.0;
178 :			M00 *= rnorm;
179 :			M01 *= rnorm;
180 :			M02 *= rnorm;
181 :			M11 *= rnorm;
182 :			M12 *= rnorm;
183 :			M22 *= rnorm;
184 :
185 :			/* this code is a mix of prior Teem code and ideas from Eberly's
186 :			"Eigensystems for 3 x 3 Symmetric Matrices (Revisited)" */
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 :			REAL mm, ss, cc;
195 :			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 + SQRT(3.0)*ss);
201 :			eval[2] = mm*(-cc - SQRT(3.0)*ss);
202 :			roots = ROOT_THREE;
203 :			/* else D is near enough to zero */
204 :			} else if (R < -epsilon || epsilon < R) {
205 :			REAL U;
206 :			/* one double root and one single root */
207 :			U = airCbrt(R); /* cube root function */
208 :			if (U > 0) {
209 :			eval[0] = 2*U;
210 :			eval[1] = -U;
211 :			eval[2] = -U;
212 :			} else {
213 :			eval[0] = -U;
214 :			eval[1] = -U;
215 :			eval[2] = 2*U;
216 :			}
217 :			roots = ROOT_SINGLE_DOUBLE;
218 :			} else {
219 :			/* a triple root! */
220 :			eval[0] = eval[1] = eval[2] = 0.0;
221 :			roots = ROOT_TRIPLE;
222 :			}
223 :
224 :			/* r0, r1, r2 are the vectors we manipulate to
225 :			find the nullspaces of M - lambda*I */
226 :			VEC_SET(r0, 0.0, M01, M02);
227 :			VEC_SET(r1, M01, 0.0, M12);
228 :			VEC_SET(r2, M02, M12, 0.0);
229 :			if (ROOT_THREE == roots) {
230 :			r0[0] = M00 - eval[0]; r1[1] = M11 - eval[0]; r2[2] = M22 - eval[0];
231 :			nullspace1(evec+0, r0, r1, r2);
232 :			r0[0] = M00 - eval[1]; r1[1] = M11 - eval[1]; r2[2] = M22 - eval[1];
233 :			nullspace1(evec+3, r0, r1, r2);
234 :			r0[0] = M00 - eval[2]; r1[1] = M11 - eval[2]; r2[2] = M22 - eval[2];
235 :			nullspace1(evec+6, r0, r1, r2);
236 :			} else if (ROOT_SINGLE_DOUBLE == roots) {
237 :			if (eval[1] == eval[2]) {
238 :			/* one big (eval[0]) , two small (eval[1,2]) */
239 :			r0[0] = M00 - eval[0]; r1[1] = M11 - eval[0]; r2[2] = M22 - eval[0];
240 :			nullspace1(evec+0, r0, r1, r2);
241 :			r0[0] = M00 - eval[1]; r1[1] = M11 - eval[1]; r2[2] = M22 - eval[1];
242 :			nullspace2(evec+3, evec+6, r0, r1, r2);
243 :			}
244 :			else {
245 :			/* two big (eval[0,1]), one small (eval[2]) */
246 :			r0[0] = M00 - eval[0]; r1[1] = M11 - eval[0]; r2[2] = M22 - eval[0];
247 :			nullspace2(evec+0, evec+3, r0, r1, r2);
248 :			r0[0] = M00 - eval[2]; r1[1] = M11 - eval[2]; r2[2] = M22 - eval[2];
249 :			nullspace1(evec+6, r0, r1, r2);
250 :			}
251 :			} else {
252 :			/* ROOT_TRIPLE == roots; use any basis for eigenvectors */
253 :			VEC_SET(evec+0, 1, 0, 0);
254 :			VEC_SET(evec+3, 0, 1, 0);
255 :			VEC_SET(evec+6, 0, 0, 1);
256 :			}
257 :			/* we always make sure its really orthonormal; keeping fixed the
258 :			eigenvector associated with the largest-magnitude eigenvalue */
259 :			if (ABS(eval[0]) > ABS(eval[2])) {
260 :			/* normalize evec+0 but don't move it */
261 :			VEC_NORM(evec+0, len);
262 :			dot = VEC_DOT(evec+0, evec+3); VEC_SCL_SUB(evec+3, dot, evec+0);
263 :			VEC_NORM(evec+3, len);
264 :			ELL_3V_CROSS(evec+6, evec+0, evec+3);
265 :			} else {
266 :			/* normalize evec+6 but don't move it */
267 :			VEC_NORM(evec+6, len);
268 :			dot = VEC_DOT(evec+6, evec+3); VEC_SCL_SUB(evec+3, dot, evec+6);
269 :			VEC_NORM(evec+3, len);
270 :			ELL_3V_CROSS(evec+0, evec+3, evec+6);
271 :			}
272 :
273 :			/* multiply back by eigenvalue L2 norm */
274 :			eval[0] /= rnorm;
275 :			eval[1] /= rnorm;
276 :			eval[2] /= rnorm;
277 :
278 :			/* add back in the eigenvalue mean */
279 :			eval[0] += mean;
280 :			eval[1] += mean;
281 :			eval[2] += mean;
282 :
283 :			return roots;
284 :			}
285 :
286 :	jhr	1597	int main (int argc, const char *argv[])
287 :			{
288 :	glk	1582	hestOpt *hopt=NULL;
289 :			hestParm *hparm;
290 :			airArray *mop;
291 :			const char *me;
292 :
293 :			#define DATA_IN_STR "../../data/lungcrop.nrrd"
294 :			unsigned int gridPoints, stepsMax;
295 :			Nrrd *nin, *nout, *npos;
296 :			char *err, *outS;
297 :			float epsilon;
298 :	glk	1589	double travelMax, strnMin;
299 :	glk	1582
300 :			mop = airMopNew();
301 :			hparm = hestParmNew();
302 :			airMopAdd(mop, hparm, (airMopper)hestParmFree, airMopAlways);
303 :			nout = nrrdNew();
304 :			airMopAdd(mop, nout, (airMopper)nrrdNuke, airMopAlways);
305 :
306 :			hparm->noArgsIsNoProblem = AIR_TRUE;
307 :			me = argv[0];
308 :			hestOptAdd(&hopt, "gridSize", "grid points", airTypeUInt, 1, 1,
309 :			&gridPoints, "120",
310 :			"number of points along volume edge");
311 :			hestOptAdd(&hopt, "stepsMax", "# iters", airTypeUInt, 1, 1,
312 :	glk	1589	&stepsMax, "30",
313 :	glk	1582	"# iterations allowed for feature detection");
314 :			hestOptAdd(&hopt, "epsilon", "convg", airTypeFloat, 1, 1,
315 :	glk	1589	&epsilon, "0.001", "convergence threshold");
316 :			hestOptAdd(&hopt, "travelMax", "tmax", airTypeDouble, 1, 1,
317 :			&travelMax, "5.0", "max on travel");
318 :			hestOptAdd(&hopt, "strnMin", "smin", airTypeDouble, 1, 1,
319 :			&strnMin, "150", "minimum significant features strength");
320 :			hestOptAdd(&hopt, "o", "out", airTypeString, 1, 1,
321 :			&outS, "ridge-pos.nrrd", "output filename");
322 :	glk	1582
323 :			hestParseOrDie(hopt, argc-1, argv+1, hparm,
324 :			me, "ridge3d", AIR_TRUE, AIR_TRUE, AIR_TRUE);
325 :			airMopAdd(mop, hopt, (airMopper)hestOptFree, airMopAlways);
326 :			airMopAdd(mop, hopt, (airMopper)hestParseFree, airMopAlways);
327 :
328 :			nin = nrrdNew();
329 :			airMopAdd(mop, nin, (airMopper)nrrdNuke, airMopAlways);
330 :			if (nrrdLoad(nin, DATA_IN_STR, NULL)) {
331 :			airMopAdd(mop, err = biffGetDone(NRRD), airFree, airMopAlways);
332 :			fprintf(stderr, "Trouble reading:\n%s", err);
333 :			airMopError(mop);
334 :			return 1;
335 :			}
336 :
337 :			gageContext *ctx;
338 :			gagePerVolume *pvl;
339 :			double kparm[NRRD_KERNEL_PARMS_NUM] = {0.0};
340 :			ctx = gageContextNew();
341 :			airMopAdd(mop, ctx, (airMopper)gageContextNix, airMopAlways);
342 :			pvl = gagePerVolumeNew(ctx, nin, gageKindScl);
343 :			if (gagePerVolumeAttach(ctx, pvl)
344 :			|| gageKernelSet(ctx, gageKernel00, nrrdKernelBSpline3, kparm)
345 :			|| gageKernelSet(ctx, gageKernel11, nrrdKernelBSpline3D, kparm)
346 :			|| gageKernelSet(ctx, gageKernel22, nrrdKernelBSpline3DD, kparm)
347 :			|| gageQueryItemOn(ctx, pvl, gageSclGradVec)
348 :			|| gageQueryItemOn(ctx, pvl, gageSclHessian)
349 :			|| gageUpdate(ctx)) {
350 :			airMopAdd(mop, err = biffGetDone(GAGE), airFree, airMopAlways);
351 :	glk	1589	fprintf(stderr, "%s: Trouble setting up gage context:\n%s", me, err);
352 :	glk	1582	airMopError(mop);
353 :			return 1;
354 :			}
355 :	glk	1589
356 :	glk	1582	npos = nrrdNew();
357 :			airMopAdd(mop, npos, (airMopper)nrrdNuke, airMopAlways);
358 :	glk	1589	if (nrrdAlloc_va(npos, nrrdTypeFloat, 2,
359 :			AIR_CAST(size_t, 3),
360 :			AIR_CAST(size_t, gridPoints*gridPoints*gridPoints))) {
361 :			airMopAdd(mop, err = biffGetDone(NRRD), airFree, airMopAlways);
362 :			fprintf(stderr, "%s: Trouble setting up nrrd output:\n%s", me, err);
363 :			airMopError(mop);
364 :			return 1;
365 :			}
366 :			float *pos = AIR_CAST(float *, npos->data);
367 :	glk	1582
368 :	glk	1589	int ui, vi, wi, si;
369 :	glk	1591	float xx, yy, zz, wpos[4], ipos[4];
370 :			const double *grad, *hess;
371 :			REAL eval[3], evec[9];
372 :	glk	1589	unsigned int posIdx = 0;
373 :			if (!( (grad = gageAnswerPointer(ctx, pvl, gageSclGradVec))
374 :			&& (hess = gageAnswerPointer(ctx, pvl, gageSclHessian))
375 :			)) {
376 :			fprintf(stderr, "%s: got null answer pointers", me);
377 :			airMopError(mop);
378 :			return 1;
379 :			}
380 :	glk	1591	wpos[3] = 1.0;
381 :			/*
382 :			printf("%s: hit return\n", me);
383 :			fgetc(stdin);
384 :			printf("%s: here we go\n", me);
385 :			*/
386 :
387 :	glk	1601	double t0 = airTime(); // start timing
388 :	glk	1591
389 :	glk	1589	for (wi=0; wi<gridPoints; wi++) {
390 :			zz = AIR_AFFINE(-0.5, wi, gridPoints-0.5, -188.0000, -188.0000 + 140.0*0.5); /* HEY */
391 :			for (vi=0; vi<gridPoints; vi++) {
392 :			yy = AIR_AFFINE(-0.5, vi, gridPoints-0.5, -175.8792, -175.8792 + 200.0*0.5); /* HEY */
393 :			for (ui=0; ui<gridPoints; ui++) {
394 :			xx = AIR_AFFINE(-0.5, ui, gridPoints-0.5, 21.6401, 21.6401 + 140.0*0.5); /* HEY */
395 :	glk	1591	ELL_3V_SET(wpos, xx, yy, zz);
396 :
397 :			REAL travel = 0.0, proj1[9], proj2[9], proj[9],
398 :	glk	1589	dir[3], len, fdd, sdd, del;
399 :	jhr	1584
400 :	glk	1591	#define PROBE(worldPos) \
401 :			ELL_4MV_MUL(ipos, ctx->shape->WtoI, worldPos); \
402 :			if ( ipos[0] <= 2 || ipos[0] >= 137 \
403 :			|| ipos[1] <= 2 || ipos[1] >= 197 \
404 :			|| ipos[2] <= 2 || ipos[2] >= 137 ) { \
405 :	glk	1589	/* if not inside */ \
406 :			break; \
407 :	glk	1591	} \
408 :			_gageProbe(ctx, ipos[0], ipos[1], ipos[2], 0.0); \
409 :			evals_evecs(eval, evec, \
410 :			hess[0], hess[1], hess[2], \
411 :			hess[4], hess[5], \
412 :			hess[8])
413 :
414 :	glk	1589	for (si=0; si<stepsMax; si++) {
415 :			if (travel > travelMax) {
416 :			break;
417 :			}
418 :	glk	1591	PROBE(wpos);
419 :	glk	1589	if (!ELL_3V_LEN(grad)) {
420 :			break;
421 :			}
422 :			ELL_3MV_OUTER(proj1, evec + 1*3, evec + 1*3);
423 :			ELL_3MV_OUTER(proj2, evec + 2*3, evec + 2*3);
424 :			ELL_3M_ADD2(proj, proj1, proj2);
425 :			ELL_3MV_MUL(dir, proj, grad);
426 :	glk	1591	VEC_NORM(dir, len);
427 :	glk	1589	if (!len) {
428 :			break;
429 :			}
430 :			fdd = ELL_3V_DOT(grad, dir);
431 :			sdd = ELL_3MV_CONTR(hess, dir);
432 :			del = sdd >= 0 ? 0.1 : -fdd/sdd;
433 :			if (del < epsilon) {
434 :			if (-eval[1] >= strnMin) {
435 :			/* converged! save point */
436 :	glk	1591	ELL_3V_COPY(pos + 3*posIdx, wpos);
437 :	glk	1589	posIdx++;
438 :			break;
439 :			}
440 :			break;
441 :			}
442 :	glk	1591	ELL_3V_SCALE_INCR(wpos, del, dir);
443 :	glk	1589	travel += del;
444 :			}
445 :			}
446 :			}
447 :			}
448 :			/* update number of actual output points */
449 :			npos->axis[1].size = posIdx;
450 :
451 :	glk	1601	double totalTime = airTime() - t0; // report timing
452 :	jhr	1584	printf("usr=%f\n", totalTime);
453 :
454 :	glk	1589	if (nrrdSave(outS, npos, NULL)) {
455 :			airMopAdd(mop, err = biffGetDone(NRRD), airFree, airMopAlways);
456 :			fprintf(stderr, "%s: Trouble saving output:\n%s", me, err);
457 :			airMopError(mop);
458 :			return 1;
459 :			}
460 :
461 :	glk	1582	airMopOkay(mop);
462 :			return 0;
463 :			}

---

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