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

SCM Repository

[diderot] Annotation of /benchmarks/programs/ridge3d/bmark-teem.c
ViewVC logotype

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

Parent Directory Parent Directory | Revision Log Revision Log


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

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