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 1586, Mon Oct 31 10:55:08 2011 UTC revision 1589, Mon Oct 31 14:50:15 2011 UTC
# Line 3  Line 3 
3   * \author Gordon Kindlmann   * \author Gordon Kindlmann
4   */   */
5    
 /*  
  * COPYRIGHT (c) 2011 The Diderot Project (http://diderot-language.cs.uchicago.edu)  
  * All rights reserved.  
  */  
   
6  #include <stdio.h>  #include <stdio.h>
7  #include <math.h>  #include <math.h>
8  #include "teem/nrrd.h"  #include "teem/nrrd.h"
# Line 34  Line 29 
29    Nrrd *nin, *nout, *npos;    Nrrd *nin, *nout, *npos;
30    char *err, *outS;    char *err, *outS;
31    float epsilon;    float epsilon;
32      double travelMax, strnMin;
33    
34    mop = airMopNew();    mop = airMopNew();
35    hparm = hestParmNew();    hparm = hestParmNew();
# Line 47  Line 43 
43               &gridPoints, "120",               &gridPoints, "120",
44               "number of points along volume edge");               "number of points along volume edge");
45    hestOptAdd(&hopt, "stepsMax", "# iters", airTypeUInt, 1, 1,    hestOptAdd(&hopt, "stepsMax", "# iters", airTypeUInt, 1, 1,
46               &stepsMax, "20",               &stepsMax, "30",
47               "# iterations allowed for feature detection");               "# iterations allowed for feature detection");
48    hestOptAdd(&hopt, "epsilon", "convg", airTypeFloat, 1, 1,    hestOptAdd(&hopt, "epsilon", "convg", airTypeFloat, 1, 1,
49               &epsilon, "0.005", "convergence threshold");               &epsilon, "0.001", "convergence threshold");
50      hestOptAdd(&hopt, "travelMax", "tmax", airTypeDouble, 1, 1,
51                 &travelMax, "5.0", "max on travel");
52      hestOptAdd(&hopt, "strnMin", "smin", airTypeDouble, 1, 1,
53                 &strnMin, "150", "minimum significant features strength");
54      hestOptAdd(&hopt, "o", "out", airTypeString, 1, 1,
55                 &outS, "ridge-pos.nrrd", "output filename");
56    
57    hestParseOrDie(hopt, argc-1, argv+1, hparm,    hestParseOrDie(hopt, argc-1, argv+1, hparm,
58                   me, "ridge3d", AIR_TRUE, AIR_TRUE, AIR_TRUE);                   me, "ridge3d", AIR_TRUE, AIR_TRUE, AIR_TRUE);
# Line 78  Line 80 
80        || gageKernelSet(ctx, gageKernel22, nrrdKernelBSpline3DD, kparm)        || gageKernelSet(ctx, gageKernel22, nrrdKernelBSpline3DD, kparm)
81        || gageQueryItemOn(ctx, pvl, gageSclGradVec)        || gageQueryItemOn(ctx, pvl, gageSclGradVec)
82        || gageQueryItemOn(ctx, pvl, gageSclHessian)        || gageQueryItemOn(ctx, pvl, gageSclHessian)
83          || gageQueryItemOn(ctx, pvl, gageSclHessEval)
84          || gageQueryItemOn(ctx, pvl, gageSclHessEvec)
85        || gageUpdate(ctx)) {        || gageUpdate(ctx)) {
86      airMopAdd(mop, err = biffGetDone(GAGE), airFree, airMopAlways);      airMopAdd(mop, err = biffGetDone(GAGE), airFree, airMopAlways);
87      fprintf(stderr, "Trouble setting up gage context:\n%s", err);      fprintf(stderr, "%s: Trouble setting up gage context:\n%s", me, err);
88      airMopError(mop);      airMopError(mop);
89      return 1;      return 1;
90    }    }
91    
92    npos = nrrdNew();    npos = nrrdNew();
93    airMopAdd(mop, npos, (airMopper)nrrdNuke, airMopAlways);    airMopAdd(mop, npos, (airMopper)nrrdNuke, airMopAlways);
94      if (nrrdAlloc_va(npos, nrrdTypeFloat, 2,
95    int ui, vi, wi;                     AIR_CAST(size_t, 3),
96    float xx, yy, zz;                     AIR_CAST(size_t, gridPoints*gridPoints*gridPoints))) {
97        airMopAdd(mop, err = biffGetDone(NRRD), airFree, airMopAlways);
98        fprintf(stderr, "%s: Trouble setting up nrrd output:\n%s", me, err);
99        airMopError(mop);
100        return 1;
101      }
102      float *pos = AIR_CAST(float *, npos->data);
103    
104    double t0 = GetTime(); // start timing    double t0 = GetTime(); // start timing
105    
106  /* COMPUTATION GOES HERE */    int ui, vi, wi, si;
107      float xx, yy, zz, pp[3];
108      const double *grad, *hess, *eval, *evec;
109      unsigned int posIdx = 0;
110      if (!( (grad = gageAnswerPointer(ctx, pvl, gageSclGradVec))
111             && (hess = gageAnswerPointer(ctx, pvl, gageSclHessian))
112             && (eval = gageAnswerPointer(ctx, pvl, gageSclHessEval))
113             && (evec = gageAnswerPointer(ctx, pvl, gageSclHessEvec))
114             )) {
115        fprintf(stderr, "%s: got null answer pointers", me);
116        airMopError(mop);
117        return 1;
118      }
119      for (wi=0; wi<gridPoints; wi++) {
120        zz = AIR_AFFINE(-0.5, wi, gridPoints-0.5, -188.0000, -188.0000 + 140.0*0.5); /* HEY */
121        for (vi=0; vi<gridPoints; vi++) {
122          yy = AIR_AFFINE(-0.5, vi, gridPoints-0.5, -175.8792, -175.8792 + 200.0*0.5); /* HEY */
123          for (ui=0; ui<gridPoints; ui++) {
124            xx = AIR_AFFINE(-0.5, ui, gridPoints-0.5, 21.6401, 21.6401 + 140.0*0.5); /* HEY */
125            ELL_3V_SET(pp, xx, yy, zz);
126            double travel = 0.0, proj1[9], proj2[9], proj[9],
127              dir[3], len, fdd, sdd, del;
128    
129    #define PROBE(posvec)                                                \
130            gageProbeSpace(ctx, posvec[0], posvec[1], posvec[2],         \
131                           AIR_FALSE /* index */, AIR_TRUE /* clamp */); \
132            if (0.0 != ctx->edgeFrac) {                                  \
133              /* if not inside */                                        \
134              break;                                                     \
135            }
136    
137            for (si=0; si<stepsMax; si++) {
138              if (travel > travelMax) {
139                break;
140              }
141              PROBE(pp);
142              if (!ELL_3V_LEN(grad)) {
143                break;
144              }
145              ELL_3MV_OUTER(proj1, evec + 1*3, evec + 1*3);
146              ELL_3MV_OUTER(proj2, evec + 2*3, evec + 2*3);
147              ELL_3M_ADD2(proj, proj1, proj2);
148              ELL_3MV_MUL(dir, proj, grad);
149              ELL_3V_NORM(dir, dir, len);
150              if (!len) {
151                break;
152              }
153              fdd = ELL_3V_DOT(grad, dir);
154              sdd = ELL_3MV_CONTR(hess, dir);
155              del = sdd >= 0 ? 0.1 : -fdd/sdd;
156              if (del < epsilon) {
157                if (-eval[1] >= strnMin) {
158                  /* converged! save point */
159                  ELL_3V_COPY(pos + 3*posIdx, pp);
160                  posIdx++;
161                  break;
162                }
163                break;
164              }
165              ELL_3V_SCALE_INCR(pp, del, dir);
166              travel += del;
167            }
168          }
169        }
170      }
171      /* update number of actual output points */
172      npos->axis[1].size = posIdx;
173    
174    
175    double totalTime = GetTime() - t0;  // report timing    double totalTime = GetTime() - t0;  // report timing
176    printf("usr=%f\n", totalTime);    printf("usr=%f\n", totalTime);
177    
178      if (nrrdSave(outS, npos, NULL)) {
179        airMopAdd(mop, err = biffGetDone(NRRD), airFree, airMopAlways);
180        fprintf(stderr, "%s: Trouble saving output:\n%s", me, err);
181        airMopError(mop);
182        return 1;
183      }
184    
185    airMopOkay(mop);    airMopOkay(mop);
186    return 0;    return 0;
187  }  }

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

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