/*! \file bmark-teem.c
 *
 * \author Gordon Kindlmann
 */

/*
 * COPYRIGHT (c) 2011 The Diderot Project (http://diderot-language.cs.uchicago.edu)
 * All rights reserved.
 */

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

static double GetTime ()
{
    struct timeval tv;

    gettimeofday (&tv, 0);

    return (double)tv.tv_sec + 0.000001 * (double)tv.tv_usec;
}

#define STATIC_INLINE static inline

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;

  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, "20",
             "# iterations allowed for feature detection");
  hestOptAdd(&hopt, "epsilon", "convg", airTypeFloat, 1, 1,
             &epsilon, "0.005", "convergence threshold");

  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)
      || gageUpdate(ctx)) {
    airMopAdd(mop, err = biffGetDone(GAGE), airFree, airMopAlways);
    fprintf(stderr, "Trouble setting up gage context:\n%s", err);
    airMopError(mop); 
    return 1;
  }
  
  npos = nrrdNew();
  airMopAdd(mop, npos, (airMopper)nrrdNuke, airMopAlways);

  int ui, vi, wi;
  float xx, yy, zz;

  double t0 = GetTime(); // start timing

/* COMPUTATION GOES HERE */

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

  airMopOkay(mop);
  return 0;
}

<script type="text/javascript">//<![CDATA[
		function toggle_ffErrors() {
			var errorsblock = document.getElementById("ffErrorsBlock");
			var errorsgroup = document.getElementById("ffErrors");
			if (errorsblock.style.display == "none") {
				errorsblock.style.display = "block";
				errorsgroup.style.right = "10px";
			} else {
				errorsblock.style.display = "none";
				errorsgroup.style.right = "300px";
			}
		}
//]]></script>
<div id="ffErrors">
<a href="javascript:toggle_ffErrors();">Click to toggle</a>
<div id="ffErrorsBlock">
	<div class="error">does not end with &lt;/html&gt; tag</div>
	<div class="error">does not end with &lt;/body&gt; tag</div>
	<div class="info">The output has ended thus: GetTime() - t0;  // report timing
  printf(&quot;usr=%f\n&quot;, totalTime);

  airMopOkay(mop);
  return 0;
}</div>
</div></div>
</body></html>