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 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