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

SCM Repository

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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1589 - (download) (as text) (annotate)
Mon Oct 31 14:50:15 2011 UTC (8 years, 9 months ago) by glk
File size: 6430 byte(s)
finished first pass at C/Teem version of ridge line detection
/*! \file bmark-teem.c
 *
 * \author Gordon Kindlmann
 */

#include <stdio.h>
#include <math.h>
#include "teem/nrrd.h"
#include "teem/gage.h"
#include "teem/ell.h"

#include "teem-defs.h"

int inside (double pos[3]) {
  if(pos[0] < -0.5 || pos[0] > 175.5) return 0;
  if(pos[1] < -0.5 || pos[1] > 186.5) return 0;
  if(pos[2] < -0.5 || pos[2] > 189.5) return 0;
  return 1;
}

int main (int argc, const char *argv[]) {
  hestOpt *hopt=NULL;
  hestParm *hparm;
  airArray *mop;
  const char *me;

#define DATA_IN_STR "../../data/lungcrop.nrrd"
  unsigned int gridPoints, stepsMax;
  Nrrd *nin, *nout, *npos;
  char *err, *outS;
  float epsilon;
  double travelMax, strnMin;

  mop = airMopNew();
  hparm = hestParmNew();
  airMopAdd(mop, hparm, (airMopper)hestParmFree, airMopAlways);
  nout = nrrdNew();
  airMopAdd(mop, nout, (airMopper)nrrdNuke, airMopAlways);

  hparm->noArgsIsNoProblem = AIR_TRUE;
  me = argv[0];
  hestOptAdd(&hopt, "gridSize", "grid points", airTypeUInt, 1, 1,
             &gridPoints, "120",
             "number of points along volume edge");
  hestOptAdd(&hopt, "stepsMax", "# iters", airTypeUInt, 1, 1,
             &stepsMax, "30",
             "# iterations allowed for feature detection");
  hestOptAdd(&hopt, "epsilon", "convg", airTypeFloat, 1, 1,
             &epsilon, "0.001", "convergence threshold");
  hestOptAdd(&hopt, "travelMax", "tmax", airTypeDouble, 1, 1,
             &travelMax, "5.0", "max on travel");
  hestOptAdd(&hopt, "strnMin", "smin", airTypeDouble, 1, 1,
             &strnMin, "150", "minimum significant features strength");
  hestOptAdd(&hopt, "o", "out", airTypeString, 1, 1,
             &outS, "ridge-pos.nrrd", "output filename");

  hestParseOrDie(hopt, argc-1, argv+1, hparm,
                 me, "ridge3d", AIR_TRUE, AIR_TRUE, AIR_TRUE);
  airMopAdd(mop, hopt, (airMopper)hestOptFree, airMopAlways);
  airMopAdd(mop, hopt, (airMopper)hestParseFree, airMopAlways);

  nin = nrrdNew();
  airMopAdd(mop, nin, (airMopper)nrrdNuke, airMopAlways);
  if (nrrdLoad(nin, DATA_IN_STR, NULL)) {
    airMopAdd(mop, err = biffGetDone(NRRD), airFree, airMopAlways);
    fprintf(stderr, "Trouble reading:\n%s", err);
    airMopError(mop); 
    return 1;
  }

  gageContext *ctx;
  gagePerVolume *pvl;
  double kparm[NRRD_KERNEL_PARMS_NUM] = {0.0};
  ctx = gageContextNew();
  airMopAdd(mop, ctx, (airMopper)gageContextNix, airMopAlways);
  pvl = gagePerVolumeNew(ctx, nin, gageKindScl);
  if (gagePerVolumeAttach(ctx, pvl)
      || gageKernelSet(ctx, gageKernel00, nrrdKernelBSpline3, kparm)
      || gageKernelSet(ctx, gageKernel11, nrrdKernelBSpline3D, kparm)
      || gageKernelSet(ctx, gageKernel22, nrrdKernelBSpline3DD, kparm)
      || gageQueryItemOn(ctx, pvl, gageSclGradVec)
      || gageQueryItemOn(ctx, pvl, gageSclHessian)
      || gageQueryItemOn(ctx, pvl, gageSclHessEval)
      || gageQueryItemOn(ctx, pvl, gageSclHessEvec)
      || gageUpdate(ctx)) {
    airMopAdd(mop, err = biffGetDone(GAGE), airFree, airMopAlways);
    fprintf(stderr, "%s: Trouble setting up gage context:\n%s", me, err);
    airMopError(mop); 
    return 1;
  }

  npos = nrrdNew();
  airMopAdd(mop, npos, (airMopper)nrrdNuke, airMopAlways);
  if (nrrdAlloc_va(npos, nrrdTypeFloat, 2, 
                   AIR_CAST(size_t, 3),
                   AIR_CAST(size_t, gridPoints*gridPoints*gridPoints))) {
    airMopAdd(mop, err = biffGetDone(NRRD), airFree, airMopAlways);
    fprintf(stderr, "%s: Trouble setting up nrrd output:\n%s", me, err);
    airMopError(mop); 
    return 1;
  }
  float *pos = AIR_CAST(float *, npos->data);

  double t0 = GetTime(); // start timing

  int ui, vi, wi, si;
  float xx, yy, zz, pp[3];
  const double *grad, *hess, *eval, *evec;
  unsigned int posIdx = 0;
  if (!( (grad = gageAnswerPointer(ctx, pvl, gageSclGradVec)) 
         && (hess = gageAnswerPointer(ctx, pvl, gageSclHessian))
         && (eval = gageAnswerPointer(ctx, pvl, gageSclHessEval))
         && (evec = gageAnswerPointer(ctx, pvl, gageSclHessEvec))
         )) {
    fprintf(stderr, "%s: got null answer pointers", me);
    airMopError(mop); 
    return 1;
  }
  for (wi=0; wi<gridPoints; wi++) {
    zz = AIR_AFFINE(-0.5, wi, gridPoints-0.5, -188.0000, -188.0000 + 140.0*0.5); /* HEY */
    for (vi=0; vi<gridPoints; vi++) {
      yy = AIR_AFFINE(-0.5, vi, gridPoints-0.5, -175.8792, -175.8792 + 200.0*0.5); /* HEY */
      for (ui=0; ui<gridPoints; ui++) {
        xx = AIR_AFFINE(-0.5, ui, gridPoints-0.5, 21.6401, 21.6401 + 140.0*0.5); /* HEY */
        ELL_3V_SET(pp, xx, yy, zz);
        double travel = 0.0, proj1[9], proj2[9], proj[9], 
          dir[3], len, fdd, sdd, del;

#define PROBE(posvec)                                                \
        gageProbeSpace(ctx, posvec[0], posvec[1], posvec[2],         \
                       AIR_FALSE /* index */, AIR_TRUE /* clamp */); \
        if (0.0 != ctx->edgeFrac) {                                  \
          /* if not inside */                                        \
          break;                                                     \
        }
        
        for (si=0; si<stepsMax; si++) {
          if (travel > travelMax) {
            break;
          }
          PROBE(pp);
          if (!ELL_3V_LEN(grad)) {
            break;
          }
          ELL_3MV_OUTER(proj1, evec + 1*3, evec + 1*3);
          ELL_3MV_OUTER(proj2, evec + 2*3, evec + 2*3);
          ELL_3M_ADD2(proj, proj1, proj2);
          ELL_3MV_MUL(dir, proj, grad);
          ELL_3V_NORM(dir, dir, len);
          if (!len) {
            break;
          }
          fdd = ELL_3V_DOT(grad, dir);
          sdd = ELL_3MV_CONTR(hess, dir);
          del = sdd >= 0 ? 0.1 : -fdd/sdd;
          if (del < epsilon) {
            if (-eval[1] >= strnMin) {
              /* converged! save point */
              ELL_3V_COPY(pos + 3*posIdx, pp);
              posIdx++;
              break;
            }
            break;
          }
          ELL_3V_SCALE_INCR(pp, del, dir);
          travel += del;
        }
      }
    }
  }
  /* update number of actual output points */
  npos->axis[1].size = posIdx;


  double totalTime = GetTime() - t0;  // report timing
  printf("usr=%f\n", totalTime);

  if (nrrdSave(outS, npos, NULL)) {
    airMopAdd(mop, err = biffGetDone(NRRD), airFree, airMopAlways);
    fprintf(stderr, "%s: Trouble saving output:\n%s", me, err);
    airMopError(mop); 
    return 1;
  }

  airMopOkay(mop);
  return 0;
}

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