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

SCM Repository

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

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

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1589, Mon Oct 31 14:50:15 2011 UTC revision 1591, Mon Oct 31 19:31:50 2011 UTC
# Line 9  Line 9 
9  #include "teem/gage.h"  #include "teem/gage.h"
10  #include "teem/ell.h"  #include "teem/ell.h"
11    
12    extern int
13    _gageProbe(gageContext *ctx, double _xi, double _yi, double _zi, double _si);
14    
15    
16  #include "teem-defs.h"  #include "teem-defs.h"
17    
18  int inside (double pos[3]) {  
19    if(pos[0] < -0.5 || pos[0] > 175.5) return 0;  #define ROOT_TRIPLE 2           /* ell_cubic_root_triple */
20    if(pos[1] < -0.5 || pos[1] > 186.5) return 0;  #define ROOT_SINGLE_DOUBLE 3    /* ell_cubic_root_single_double */
21    if(pos[2] < -0.5 || pos[2] > 189.5) return 0;  #define ROOT_THREE 4            /* ell_cubic_root_three */
22    return 1;  
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    }
102    
103    /*
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  int main (int argc, const char *argv[]) {  int main (int argc, const char *argv[]) {
294    hestOpt *hopt=NULL;    hestOpt *hopt=NULL;
295    hestParm *hparm;    hestParm *hparm;
# Line 80  Line 352 
352        || gageKernelSet(ctx, gageKernel22, nrrdKernelBSpline3DD, kparm)        || gageKernelSet(ctx, gageKernel22, nrrdKernelBSpline3DD, kparm)
353        || gageQueryItemOn(ctx, pvl, gageSclGradVec)        || gageQueryItemOn(ctx, pvl, gageSclGradVec)
354        || gageQueryItemOn(ctx, pvl, gageSclHessian)        || gageQueryItemOn(ctx, pvl, gageSclHessian)
       || gageQueryItemOn(ctx, pvl, gageSclHessEval)  
       || gageQueryItemOn(ctx, pvl, gageSclHessEvec)  
355        || gageUpdate(ctx)) {        || gageUpdate(ctx)) {
356      airMopAdd(mop, err = biffGetDone(GAGE), airFree, airMopAlways);      airMopAdd(mop, err = biffGetDone(GAGE), airFree, airMopAlways);
357      fprintf(stderr, "%s: Trouble setting up gage context:\n%s", me, err);      fprintf(stderr, "%s: Trouble setting up gage context:\n%s", me, err);
# Line 101  Line 371 
371    }    }
372    float *pos = AIR_CAST(float *, npos->data);    float *pos = AIR_CAST(float *, npos->data);
373    
   double t0 = GetTime(); // start timing  
   
374    int ui, vi, wi, si;    int ui, vi, wi, si;
375    float xx, yy, zz, pp[3];    float xx, yy, zz, wpos[4], ipos[4];
376    const double *grad, *hess, *eval, *evec;    const double *grad, *hess;
377      REAL eval[3], evec[9];
378    unsigned int posIdx = 0;    unsigned int posIdx = 0;
379    if (!( (grad = gageAnswerPointer(ctx, pvl, gageSclGradVec))    if (!( (grad = gageAnswerPointer(ctx, pvl, gageSclGradVec))
380           && (hess = gageAnswerPointer(ctx, pvl, gageSclHessian))           && (hess = gageAnswerPointer(ctx, pvl, gageSclHessian))
          && (eval = gageAnswerPointer(ctx, pvl, gageSclHessEval))  
          && (evec = gageAnswerPointer(ctx, pvl, gageSclHessEvec))  
381           )) {           )) {
382      fprintf(stderr, "%s: got null answer pointers", me);      fprintf(stderr, "%s: got null answer pointers", me);
383      airMopError(mop);      airMopError(mop);
384      return 1;      return 1;
385    }    }
386      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    for (wi=0; wi<gridPoints; wi++) {    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 */      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++) {      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 */        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++) {        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 */          xx = AIR_AFFINE(-0.5, ui, gridPoints-0.5, 21.6401, 21.6401 + 140.0*0.5); /* HEY */
401          ELL_3V_SET(pp, xx, yy, zz);          ELL_3V_SET(wpos, xx, yy, zz);
402          double travel = 0.0, proj1[9], proj2[9], proj[9],  
403            REAL travel = 0.0, proj1[9], proj2[9], proj[9],
404            dir[3], len, fdd, sdd, del;            dir[3], len, fdd, sdd, del;
405    
406  #define PROBE(posvec)                                                \  #define PROBE(worldPos)                                              \
407          gageProbeSpace(ctx, posvec[0], posvec[1], posvec[2],         \          ELL_4MV_MUL(ipos, ctx->shape->WtoI, worldPos);               \
408                         AIR_FALSE /* index */, AIR_TRUE /* clamp */); \          if ( ipos[0] <= 2 || ipos[0] >= 137                          \
409          if (0.0 != ctx->edgeFrac) {                                  \               || ipos[1] <= 2 || ipos[1] >= 197                       \
410                 || ipos[2] <= 2 || ipos[2] >= 137 ) {                   \
411            /* if not inside */                                        \            /* if not inside */                                        \
412            break;                                                     \            break;                                                     \
413          }          }                                                            \
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          for (si=0; si<stepsMax; si++) {          for (si=0; si<stepsMax; si++) {
421            if (travel > travelMax) {            if (travel > travelMax) {
422              break;              break;
423            }            }
424            PROBE(pp);            PROBE(wpos);
425            if (!ELL_3V_LEN(grad)) {            if (!ELL_3V_LEN(grad)) {
426              break;              break;
427            }            }
# Line 146  Line 429 
429            ELL_3MV_OUTER(proj2, evec + 2*3, evec + 2*3);            ELL_3MV_OUTER(proj2, evec + 2*3, evec + 2*3);
430            ELL_3M_ADD2(proj, proj1, proj2);            ELL_3M_ADD2(proj, proj1, proj2);
431            ELL_3MV_MUL(dir, proj, grad);            ELL_3MV_MUL(dir, proj, grad);
432            ELL_3V_NORM(dir, dir, len);            VEC_NORM(dir, len);
433            if (!len) {            if (!len) {
434              break;              break;
435            }            }
# Line 156  Line 439 
439            if (del < epsilon) {            if (del < epsilon) {
440              if (-eval[1] >= strnMin) {              if (-eval[1] >= strnMin) {
441                /* converged! save point */                /* converged! save point */
442                ELL_3V_COPY(pos + 3*posIdx, pp);                ELL_3V_COPY(pos + 3*posIdx, wpos);
443                posIdx++;                posIdx++;
444                break;                break;
445              }              }
446              break;              break;
447            }            }
448            ELL_3V_SCALE_INCR(pp, del, dir);            ELL_3V_SCALE_INCR(wpos, del, dir);
449            travel += del;            travel += del;
450          }          }
451        }        }

Legend:
Removed from v.1589  
changed lines
  Added in v.1591

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