自学教程：C++ AIR_AFFINE函数代码示例

voidplotDots(plotPS *pps, plotParm *pparm, Nrrd **ndata, int nidx) {  int ii, npts;  double xx, yy, orad, irad, *data, val;  if (!( pparm->dotDiameter[nidx] > 0 )) {    return;  }  fprintf(pps->file, "gsave/n");  fprintf(pps->file, "newpath/n");  plotWidth(pps, pparm, 0);  data = (double *)(ndata[nidx]->data);  npts = ndata[nidx]->axis[1].size;  orad = pparm->dotDiameter[nidx]/2;  irad = pparm->dotInnerDiameterFraction*orad;  for (ii=0; ii<npts; ii++) {    val = data[ii];    xx = AIR_AFFINE(0, ii, npts-1,                    ndata[nidx]->axis[1].min, ndata[nidx]->axis[1].max);    xx = AIR_AFFINE(pparm->dbox[0], xx, pparm->dbox[2],                    pps->bbox[0], pps->bbox[2]);    yy = AIR_AFFINE(pparm->dbox[1], val, pparm->dbox[3],                    pps->bbox[1], pps->bbox[3]);    plotGray(pps, pparm, pparm->dotGray[nidx]);    fprintf(pps->file, "%g %g %g 0 360 arc closepath fill/n",             PPS_X(xx), PPS_Y(yy), PPS_S(orad));    if (irad) {      plotGray(pps, pparm, 1.0);      fprintf(pps->file, "%g %g %g 0 360 arc closepath fill/n",               PPS_X(xx), PPS_Y(yy), PPS_S(irad));    }  }  fprintf(pps->file, "grestore/n");}

/********* tenGradientRandom**** generates num random unit vectors of type double*/inttenGradientRandom(Nrrd *ngrad, unsigned int num, unsigned int seed) {  static const char me[]="tenGradientRandom";  double *grad, len;  unsigned int gi;  if (nrrdMaybeAlloc_va(ngrad, nrrdTypeDouble, 2,                        AIR_CAST(size_t, 3), AIR_CAST(size_t, num))) {    biffMovef(TEN, NRRD, "%s: couldn't allocate output", me);    return 1;  }  airSrandMT(seed);  grad = AIR_CAST(double*, ngrad->data);  for (gi=0; gi<num; gi++) {    do {      grad[0] = AIR_AFFINE(0, airDrandMT(), 1, -1, 1);      grad[1] = AIR_AFFINE(0, airDrandMT(), 1, -1, 1);      grad[2] = AIR_AFFINE(0, airDrandMT(), 1, -1, 1);      len = ELL_3V_LEN(grad);    } while (len > 1 || !len);    ELL_3V_SCALE(grad, 1.0/len, grad);    grad += 3;  }  return 0;}

voidplotGraph(plotPS *pps, plotParm *pparm, Nrrd **ndata, int nidx) {  int ii, npts;  double xx, yy, *data, val;  if (!( pparm->graphThick[nidx] > 0 )) {    return;  }  data = (double *)(ndata[nidx]->data);  npts = ndata[nidx]->axis[1].size;  plotGray(pps, pparm, pparm->graphGray[nidx]);  fprintf(pps->file, "%g W/n", pps->psc*pparm->graphThick[nidx]);  for (ii=0; ii<npts; ii++) {    val = data[ii];    xx = AIR_AFFINE(0, ii, npts-1,                    ndata[nidx]->axis[1].min, ndata[nidx]->axis[1].max);    xx = AIR_AFFINE(pparm->dbox[0], xx, pparm->dbox[2],                    pps->bbox[0], pps->bbox[2]);    yy = AIR_AFFINE(pparm->dbox[1], val, pparm->dbox[3],                    pps->bbox[1], pps->bbox[3]);    fprintf(pps->file, "%g %g %s/n", PPS_X(xx), PPS_Y(yy),            ii ? "L" : "M");  }  fprintf(pps->file, "S/n");}

void_cap2xyz(double xyz[3], double ca, double cp, int version, int whole) {  double cl, cs, mean;  cs = 1 - ca;  cl = 1 - cs - cp;  mean = (cs + cp + cl)/3;  /*    xyz[0] = cs*0.333 + cl*1.0 + cp*0.5;    xyz[1] = cs*0.333 + cl*0.0 + cp*0.5;    xyz[2] = cs*0.333 + cl*0.0 + cp*0.0;    xyz[0] = AIR_AFFINE(0, ca, 1, 1.1*xyz[0], 0.86*xyz[0]);    xyz[1] = AIR_AFFINE(0, ca, 1, 1.1*xyz[1], 0.86*xyz[1]);    xyz[2] = AIR_AFFINE(0, ca, 1, 1.1*xyz[2], 0.86*xyz[2]);  */  if (whole) {    ELL_3V_SET(xyz,               AIR_AFFINE(0.0, 0.9, 1.0, mean, cl),               AIR_AFFINE(0.0, 0.9, 1.0, mean, cp),               AIR_AFFINE(0.0, 0.9, 1.0, mean, cs));    ELL_3V_SET(xyz, cl, cp, cs);  } else {    if (1 == version) {      ELL_3V_SET(xyz,                 (3 + 3*cl - cs)/6,                 (2 - 2*cl + cp)/6,                 2*cs/6);    } else {      ELL_3V_SET(xyz, 1, 1 - cl, cs);    }  }}

/********* nrrdAxisInfoPosRange()**** given a nrrd, an axis, and two (floating point) index space positions,** return the range of positions implied the axis's min, max, and center** The opposite of nrrdAxisIdxRange()*/voidnrrdAxisInfoPosRange(double *loP, double *hiP,                     const Nrrd *nrrd, unsigned int ax,                      double loIdx, double hiIdx) {  int center, flip = 0;  size_t size;  double min, max, tmp;  if (!( loP && hiP && nrrd && ax <= nrrd->dim-1 )) {    *loP = *hiP = AIR_NAN;    return;  }  center = _nrrdCenter(nrrd->axis[ax].center);  min = nrrd->axis[ax].min;  max = nrrd->axis[ax].max;  size = nrrd->axis[ax].size;  if (loIdx > hiIdx) {    flip = 1;    tmp = loIdx; loIdx = hiIdx; hiIdx = tmp;  }  if (nrrdCenterCell == center) {    *loP = AIR_AFFINE(0, loIdx, size, min, max);    *hiP = AIR_AFFINE(0, hiIdx+1, size, min, max);  } else {    *loP = AIR_AFFINE(0, loIdx, size-1, min, max);    *hiP = AIR_AFFINE(0, hiIdx, size-1, min, max);  }  if (flip) {    tmp = *loP; *loP = *hiP; *hiP = tmp;  }  return;}

voidmakeSceneRainLights(limnCamera *cam, echoRTParm *parm, echoScene *scene) {  echoObject *sphere, *rect;  int i, N;  echoPos_t w;  float r, g, b;  ELL_3V_SET(cam->from, 2.5, 0, 5);  ELL_3V_SET(cam->at,   0, 0, 0);  ELL_3V_SET(cam->up,   0, 0, 1);  cam->uRange[0] = -1.7;  cam->uRange[1] = 1.7;  cam->vRange[0] = -1.7;  cam->vRange[1] = 1.7;  parm->jitterType = echoJitterJitter;  parm->numSamples = 36;  parm->imgResU = 1000;  parm->imgResV = 1000;  parm->numSamples = 16;  parm->imgResU = 200;  parm->imgResV = 200;  parm->aperture = 0.0;  parm->renderLights = AIR_TRUE;  parm->shadow = 0.0;  ELL_3V_SET(scene->bkgr, 0.1, 0.1, 0.1);  /* create scene */  sphere = echoObjectNew(scene, echoTypeSphere);  echoSphereSet(sphere, 0, 0, 0, 1.0);  echoColorSet(sphere, 1.0, 1.0, 1.0, 1.0);  echoMatterPhongSet(scene, sphere, 0.02, 0.2, 1.0, 400);  echoObjectAdd(scene, sphere);  N = 8;  w = 1.7/N;  for (i=0; i<N; i++) {    rect = echoObjectNew(scene, echoTypeRectangle);    echoRectangleSet(rect,                     w/2, AIR_AFFINE(0, i, N-1, -1-w/2, 1-w/2), 1.5,                     0, w, 0,                     w, 0, 0);    _dyeHSVtoRGB(&r, &g, &b, AIR_AFFINE(0, i, N, 0.0, 1.0), 1.0, 1.0);    echoColorSet(rect, r, g, b, 1);    echoMatterLightSet(scene, rect, 1, 0);    echoObjectAdd(scene, rect);  }}

echoObject *echoRoughSphereNew(echoScene *scene, int theRes, int phiRes, echoPos_t *matx) {  echoObject *trim;  echoPos_t *_pos, *pos, tmp[3];  int *_vert, *vert, thidx, phidx, n;  echoPos_t th, ph;  trim = echoObjectNew(scene, echoTypeTriMesh);  TRIMESH(trim)->numV = 2 + (phiRes-1)*theRes;  TRIMESH(trim)->numF = (2 + 2*(phiRes-2))*theRes;  _pos = pos = (echoPos_t *)calloc(3*TRIMESH(trim)->numV, sizeof(echoPos_t));  _vert = vert = (int *)calloc(3*TRIMESH(trim)->numF, sizeof(int));  ELL_3V_SET(tmp, 0, 0, 1); _echoPosSet(pos, matx, tmp); pos += 3;  for (phidx=1; phidx<phiRes; phidx++) {    ph = AIR_AFFINE(0, phidx, phiRes, 0.0, AIR_PI);    for (thidx=0; thidx<theRes; thidx++) {      th = AIR_AFFINE(0, thidx, theRes, 0.0, 2*AIR_PI);      ELL_3V_SET(tmp, cos(th)*sin(ph), sin(th)*sin(ph), cos(ph));      _echoPosSet(pos, matx, tmp); pos += 3;    }  }  ELL_3V_SET(tmp, 0, 0, -1); _echoPosSet(pos, matx, tmp);  for (thidx=0; thidx<theRes; thidx++) {    n = AIR_MOD(thidx+1, theRes);    ELL_3V_SET(vert, 0, 1+thidx, 1+n); vert += 3;  }  for (phidx=0; phidx<phiRes-2; phidx++) {    for (thidx=0; thidx<theRes; thidx++) {      n = AIR_MOD(thidx+1, theRes);      ELL_3V_SET(vert, 1+phidx*theRes+thidx, 1+(1+phidx)*theRes+thidx,                 1+phidx*theRes+n); vert += 3;      ELL_3V_SET(vert, 1+(1+phidx)*theRes+thidx, 1+(1+phidx)*theRes+n,                  1+phidx*theRes+n); vert += 3;    }  }  for (thidx=0; thidx<theRes; thidx++) {    n = AIR_MOD(thidx+1, theRes);    ELL_3V_SET(vert, 1+(phiRes-2)*theRes+thidx, TRIMESH(trim)->numV-1,               1+(phiRes-2)*theRes+n);     vert += 3;  }  echoTriMeshSet(trim, TRIMESH(trim)->numV, _pos, TRIMESH(trim)->numF, _vert);  return(trim);}

/* Calculate the Q-ball profile from DWIs */void_tenQball(const double b, const int gradcount, const double svals[],          const double grads[], double qvals[] ) {  /* Not an optimal Q-ball implementation! (Taken from submission to     MICCAI 2006) Should be solved analytically in the future,     implemented from recent papers. */  int i,j;  double d, dist, weight, min, max;  AIR_UNUSED(b);  min = max = svals[1] / svals[0];  for( i = 0; i < gradcount; i++ ) {    d = svals[i+1] / svals[0];    if( d > max )      max = d;    else if( d < min )      min = d;  }  for( i = 0; i < gradcount; i++ ) {    qvals[i] = 0;    for( j = 0; j < gradcount; j++ ) {      d = AIR_AFFINE( min, svals[j+1] / svals[0], max, 0,1 );      dist = ELL_3V_DOT(grads + 3*i, grads + 3*j);      dist = AIR_ABS(dist);      weight = cos( 0.5 * AIR_PI * dist );      qvals[i] += d * weight*weight*weight*weight;    }  }  return;}

/* XX != AB */static intlattABtoXX(int dstLatt, double *dstParm, const double *srcParm) {  double AA[2], BB[2], theta, phase, radi, area, len;  int ret = 0;  /* we have to reduce the DOF, which always starts with the same     loss of orientation information */  ELL_2V_COPY(AA, srcParm + 0);  ELL_2V_COPY(BB, srcParm + 2);  getToPosY(AA, BB);  switch(dstLatt) {  case rvaLattPRA:   /* AB -> PRA (loss off orientation) */    theta = atan2(BB[1], BB[0]);    phase = AIR_AFFINE(AIR_PI/2, theta, AIR_PI/3, 0.0, 1.0);    radi = _rvaLen2(BB)/_rvaLen2(AA);    area = _rvaLen2(AA)*BB[1];    ELL_3V_SET(dstParm, phase, radi, area);    break;  case rvaLattUVW:  /* AB -> UVW (loss of orientation) */    ELL_3V_SET(dstParm, BB[0], BB[1], AA[0]);    break;  case rvaLattXY:   /* AB -> XY (loss of orientation and scale) */    len = _rvaLen2(AA);    ELL_2V_SET(dstParm, BB[0]/len, BB[1]/len);    break;  default: ret = 1; break; /* unimplemented */  }  return ret;}

/* XX != AB */static intlattXXtoAB(double *dstParm, int srcLatt, const double *srcParm) {  double AA[2], BB[2], area, theta, radi, scl;  int ret = 0;  switch(srcLatt) {  case rvaLattPRA:  /* PRA -> AB */    area = AIR_ABS(srcParm[2]);    theta = AIR_AFFINE(0, srcParm[0], 1, AIR_PI/2, AIR_PI/3);    radi = srcParm[1];    ELL_2V_SET(AA, 1.0, 0.0);    ELL_2V_SET(BB, radi*cos(theta), radi*sin(theta));    /* area from AA and BB is BB[1], but need to scale       these to get to requested area */    scl = sqrt(area/BB[1]);    ELL_4V_SET(dstParm, scl*AA[0], scl*AA[1], scl*BB[0], scl*BB[1]);    break;  case rvaLattAB:  /* UVW -> AB */    ELL_4V_SET(dstParm, srcParm[2], 0.0, srcParm[0], srcParm[1]);    break;  case rvaLattXY:  /* XY -> AB */    ELL_4V_SET(dstParm, 1.0, 0.0, srcParm[0], srcParm[1]);    break;  default: ret = 1; break; /* unimplemented */  }  return ret;}

/* ----------------------------------------------------------------** ------------------------------ SPRING --------------------------** ----------------------------------------------------------------** 1 parms:** parm[0]: width of pull region (beyond 1.0)**** learned: "1/2" is not 0.5 !!!!!*/double_pullEnergySpringEval(double *frc, double dist, const double *parm) {  /* char me[]="_pullEnergySpringEval"; */  double enr, xx, pull;  pull = parm[0];  /* support used to be [0,1 + pull], but now is scrunched to [0,1],     so hack "dist" to match old parameterization */  dist = AIR_AFFINE(0, dist, 1, 0, 1+pull);  xx = dist - 1.0;  if (xx > pull) {    enr = 0;    *frc = 0;  } else if (xx > 0) {    enr = xx*xx*(xx*xx/(4*pull*pull) - 2*xx/(3*pull) + 1.0/2.0);    *frc = xx*(xx*xx/(pull*pull) - 2*xx/pull + 1);  } else {    enr = xx*xx/2;    *frc = xx;  }  /*  if (!AIR_EXISTS(ret)) {    fprintf(stderr, "!%s: dist=%g, pull=%g, blah=%d --> ret=%g/n",            me, dist, pull, blah, ret);  }  */  return enr;}

intlimnSplineSample(Nrrd *nout, limnSpline *spline,                 double minT, size_t M, double maxT) {  char me[]="limnSplineSample", err[BIFF_STRLEN];  airArray *mop;  Nrrd *ntt;  double *tt;  size_t I;  if (!(nout && spline)) {    sprintf(err, "%s: got NULL pointer", me);    biffAdd(LIMN, err); return 1;  }  mop = airMopNew();  airMopAdd(mop, ntt=nrrdNew(), (airMopper)nrrdNuke, airMopAlways);  if (nrrdMaybeAlloc_va(ntt, nrrdTypeDouble, 1,                        M)) {    sprintf(err, "%s: trouble allocating tmp nrrd", me);    biffMove(LIMN, err, NRRD); airMopError(mop); return 1;  }  tt = (double*)(ntt->data);  for (I=0; I<M; I++) {    tt[I] = AIR_AFFINE(0, I, M-1, minT, maxT);  }  if (limnSplineNrrdEvaluate(nout, spline, ntt)) {    sprintf(err, "%s: trouble", me);    biffAdd(LIMN, err); airMopError(mop); return 1;  }  airMopOkay(mop);  return 0;}

voidincTest(char *me, int num, baneRange *range) {  double *val, tmp, incParm[BANE_PARM_NUM], omin, omax, rmin, rmax;  baneInc *inc;  Nrrd *hist;  int i, j;  airSrand48();  val = (double*)malloc(num*sizeof(double));  /* from <http://www.itl.nist.gov/div898/handbook/index.htm>:     the standard dev of a uniform distribution between A and B is     sqrt((B-A)^2/12) */  for (j=0; j<num; j++) {    tmp = AIR_AFFINE(0.0, airDrand48(), 1.0, -1.0, 1.0);    /* val[j] = tmp*tmp*tmp; */    val[j] = tmp;  }  rmin = rmax = val[0];  for (j=0; j<num; j++) {    rmin = AIR_MIN(rmin, val[j]);    rmax = AIR_MAX(rmax, val[j]);  }  fprintf(stderr, "incTest: real min,max = %g,%g/n", rmin, rmax);    for (i=1; i<baneIncLast; i++) {    /* NOTE: THIS IS BROKEN !!! */    inc = baneIncNew(i, NULL, incParm);    printf("%s: inclusion %s ------/n", me, inc->name);    switch(i) {    case baneIncAbsolute:      ELL_3V_SET(incParm, -0.8, 1.5, AIR_NAN);      break;    case baneIncRangeRatio:      ELL_3V_SET(incParm, 0.99, AIR_NAN, AIR_NAN);      break;    case baneIncPercentile:      ELL_3V_SET(incParm, 1024, 10, AIR_NAN);      break;    case baneIncStdv:      ELL_3V_SET(incParm, 1.0, AIR_NAN, AIR_NAN);      break;    }    fprintf(stderr, "!%s: THIS IS BROKEN!!!/n", "incTest");    /*    if (inc->passA) {      for (j=0; j<num; j++)        inc->process[0](hist, val[j], incParm);    }    if (inc->passB) {      for (j=0; j<num; j++)        inc->process[1](hist, val[j], incParm);    }    inc->ans(&omin, &omax, hist, incParm, range);    */    printf(" --> (%g,%g)/n", omin, omax);  }  free(val);}

/********* tenGradientJitter**** moves all gradients by amount dist on tangent plane, in a random** direction, and then renormalizes. The distance is a fraction** of the ideal edge length (via tenGradientIdealEdge)*/inttenGradientJitter(Nrrd *nout, const Nrrd *nin, double dist) {  static const char me[]="tenGradientJitter";  double *grad, perp0[3], perp1[3], len, theta, cc, ss, edge;  unsigned int gi, num;  if (nrrdConvert(nout, nin, nrrdTypeDouble)) {    biffMovef(TEN, NRRD, "%s: trouble converting input to double", me);    return 1;  }  if (tenGradientCheck(nout, nrrdTypeDouble, 3)) {    biffAddf(TEN, "%s: didn't get valid gradients", me);    return 1;  }  grad = AIR_CAST(double*, nout->data);  num = AIR_UINT(nout->axis[1].size);  /* HEY: possible confusion between single and not */  edge = tenGradientIdealEdge(num, AIR_FALSE);  for (gi=0; gi<num; gi++) {    ELL_3V_NORM(grad, grad, len);    ell_3v_perp_d(perp0, grad);    ELL_3V_CROSS(perp1, perp0, grad);    theta = AIR_AFFINE(0, airDrandMT(), 1, 0, 2*AIR_PI);    cc = dist*edge*cos(theta);    ss = dist*edge*sin(theta);    ELL_3V_SCALE_ADD3(grad, 1.0, grad, cc, perp0, ss, perp1);    ELL_3V_NORM(grad, grad, len);    grad += 3;  }  return 0;}

/*** Eval from UV*/voidtenGlyphBqdEvalUv(double eval[3], const double uv[2]) {  double xx, yy, zz, ll;  yy = AIR_AFFINE(0, uv[0], 1, -1, 1);  if (uv[0] + uv[1] > 1) {    zz = AIR_AFFINE(0, uv[1], 1, -1, 1) - 1 + yy;    xx = 1;  } else {    xx = AIR_AFFINE(0, uv[1], 1, -1, 1) + yy + 1;    zz = -1;  }  ELL_3V_SET(eval, xx, yy, zz);  ELL_3V_NORM(eval, eval, ll);  return;}

void_echoRayIntxUV_Sphere(echoIntx *intx) {  echoPos_t u, v;  if (intx->norm[0] || intx->norm[1]) {    u = atan2(intx->norm[1], intx->norm[0]);    intx->u = AIR_AFFINE(-AIR_PI, u, AIR_PI, 0.0, 1.0);    v = -asin(intx->norm[2]);    intx->v = AIR_AFFINE(-AIR_PI/2, v, AIR_PI/2, 0.0, 1.0);  }  else {    intx->u = 0;    /* this is valid because if we're here, then intx->norm[2]       is either 1.0 or -1.0 */    intx->v = AIR_AFFINE(1.0, intx->norm[2], -1.0, 0.0, 1.0);  }}

/*** hist[0]: hue vs sat** hist[1]: hue vs val*/voidimageProc(Nrrd *nhproj[3], Nrrd *nhist[2], unsigned int sH,          float *rgb, unsigned int size0, unsigned int sXY,          unsigned int overSampleNum, float overSampleScale) {  unsigned int xyi, hi, si, vi, oi;  float rr, gg, bb, hh, ss, vv, *hist[2];  double rndA, rndB;  nrrdZeroSet(nhist[0]);  nrrdZeroSet(nhist[1]);  hist[0] = AIR_CAST(float *, nhist[0]->data);  hist[1] = AIR_CAST(float *, nhist[1]->data);  for (xyi=0; xyi<sXY; xyi++) {    rr = AIR_CLAMP(0, rgb[0], 255);    gg = AIR_CLAMP(0, rgb[1], 255);    bb = AIR_CLAMP(0, rgb[2], 255);    rr = AIR_AFFINE(-1, rr, 256, 0, 1);    gg = AIR_AFFINE(-1, gg, 256, 0, 1);    bb = AIR_AFFINE(-1, bb, 256, 0, 1);    dyeRGBtoHSV(&hh, &ss, &vv, rr, gg, bb);    si = airIndexClamp(0, ss, 1, sH);    vi = airIndexClamp(0, vv, 1, sH);#define UPDATE_HIST(rnd)                                                /    hi = airIndexClamp(0, hh + overSampleScale*(1-ss)*(rnd), 1, sH);    /    hist[0][hi + sH*si] += 1.0/overSampleNum;                           /    hist[1][hi + sH*vi] += 1.0/overSampleNum    if (overSampleNum % 2 == 1) {      airNormalRand(&rndA, NULL);      UPDATE_HIST(rndA);      overSampleNum -= 1;    }    for (oi=0; oi<overSampleNum; oi+=2) {      airNormalRand(&rndA, &rndB);      UPDATE_HIST(rndA);      UPDATE_HIST(rndB);    }    rgb += size0;  }  nrrdProject(nhproj[0], nhist[0], 1, nrrdMeasureHistoMean, nrrdTypeFloat);  nrrdProject(nhproj[1], nhist[1], 1, nrrdMeasureHistoMean, nrrdTypeFloat);  nrrdProject(nhproj[2], nhist[1], 1, nrrdMeasureSum, nrrdTypeFloat);}

intmain(int argc, char *argv[]) {  char *me;  hestOpt *hopt=NULL;  airArray *mop;    pushEnergySpec *ensp;  unsigned int pi, xi, nn;  double xx, supp, del;  mop = airMopNew();  me = argv[0];  hestOptAdd(&hopt, "energy", "spec", airTypeOther, 1, 1, &ensp, NULL,             "specification of force function to use",             NULL, NULL, pushHestEnergySpec);  hestParseOrDie(hopt, argc-1, argv+1, NULL,                 me, info, AIR_TRUE, AIR_TRUE, AIR_TRUE);  airMopAdd(mop, hopt, (airMopper)hestOptFree, airMopAlways);  airMopAdd(mop, hopt, (airMopper)hestParseFree, airMopAlways);  fprintf(stderr, "%s: parsed energy /"%s/", with %u parms./n", me,          ensp->energy->name, ensp->energy->parmNum);  for (pi=0; pi<ensp->energy->parmNum; pi++) {    fprintf(stderr, "%u: %g/n", pi, ensp->parm[pi]);  }  fprintf(stderr, "/n");  nn = 600;  supp = ensp->energy->support(ensp->parm);  del = AIR_DELTA(0, 2, nn, 0, supp);  for (xi=1; xi<nn; xi++) {    double x0, x1, ee, ff, e0, e1, dummy;    xx = AIR_AFFINE(0, xi,   nn, 0, supp);    x1 = AIR_AFFINE(0, xi+1, nn, 0, supp);    x0 = AIR_AFFINE(0, xi-1, nn, 0, supp);    ensp->energy->eval(&e1, &dummy, x1, ensp->parm);    ensp->energy->eval(&e0, &dummy, x0, ensp->parm);    ensp->energy->eval(&ee, &ff, xx, ensp->parm);    printf("%g %g %g %g/n", xx, ee, ff, (e1 - e0)/del);  }  airMopOkay(mop);  return 0;}

void_echoRayIntxUV_TriMesh(echoIntx *intx) {  echoPos_t u, v, norm[3], len;  echoTriMesh *trim;  trim = TRIMESH(intx->obj);  ELL_3V_SUB(norm, intx->pos, trim->meanvert);  ELL_3V_NORM(norm, norm, len);  if (norm[0] || norm[1]) {    u = atan2(norm[1], norm[0]);    intx->u = AIR_AFFINE(-AIR_PI, u, AIR_PI, 0.0, 1.0);    v = -asin(norm[2]);    intx->v = AIR_AFFINE(-AIR_PI/2, v, AIR_PI/2, 0.0, 1.0);  }  else {    intx->u = 0;    intx->v = AIR_AFFINE(1.0, norm[2], -1.0, 0.0, 1.0);  }}

intmain(int argc, char *argv[]) {  char *me;  unsigned int ii, NN;  double min, max, *out;  Nrrd *nout;  me = argv[0];  if (5 != argc) {    usage(me);  }  if (3 != (sscanf(argv[2], "%u", &NN)            + sscanf(argv[1], "%lf", &min)            + sscanf(argv[3], "%lf", &max))) {    fprintf(stderr, "%s: couldn't parse %s %s %s double uint double",            me, argv[1], argv[2], argv[3]);    usage(me);  }  nout = nrrdNew();  if (nrrdAlloc_va(nout, nrrdTypeDouble, 2,                   AIR_CAST(size_t, 5),                   AIR_CAST(size_t, NN))) {    fprintf(stderr, "%s: trouble allocating:/n%s", me,            biffGetDone(NRRD));    exit(1);  }  out = AIR_CAST(double *, nout->data);  for (ii=0; ii<NN; ii++) {    double xx, rr, ff, gg;    xx = AIR_AFFINE(0, ii, NN-1, min, max);    rr = exp(xx);    ff = airFastExp(xx);    gg = airExp(xx);    if (rr < 0 || ff < 0 || gg < 0        || !AIR_EXISTS(rr) || !AIR_EXISTS(ff) || !AIR_EXISTS(gg)) {      fprintf(stderr, "%s: problem: %f -> real %f, approx %f %f/n",              me, xx, rr, ff, gg);      exit(1);    }    out[0] = rr;    out[1] = ff;    out[2] = (ff-rr)/rr;    out[3] = gg;    out[4] = (gg-rr)/rr;    out += 5;  }  if (nrrdSave(argv[4], nout, NULL)) {    fprintf(stderr, "%s: trouble saving:/n%s", me,            biffGetDone(NRRD));    exit(1);  }  exit(0);}

/********* nrrd1DIrregAclGenerate()**** Generates the "acl" that is used to speed up the action of** nrrdApply1DIrregMap().  Basically, the domain of the map** is quantized into "acllen" bins, and for each bin, the** lowest and highest possible map interval is stored. This** either obviates or speeds up the task of finding which** interval contains a given value.**** Assumes that nrrd1DIrregMapCheck has been called on "nmap".*/intnrrd1DIrregAclGenerate(Nrrd *nacl, const Nrrd *nmap, size_t aclLen) {  char me[]="nrrd1DIrregAclGenerate", err[BIFF_STRLEN];  int posLen;  unsigned int ii;  unsigned short *acl;  double lo, hi, min, max, *pos;  if (!(nacl && nmap)) {    sprintf(err, "%s: got NULL pointer", me);    biffAdd(NRRD, err); return 1;  }  if (!(aclLen >= 2)) {    sprintf(err, "%s: given acl length (" _AIR_SIZE_T_CNV             ") is too small", me, aclLen);    biffAdd(NRRD, err); return 1;  }  if (nrrdMaybeAlloc_va(nacl, nrrdTypeUShort, 2,                         AIR_CAST(size_t, 2), AIR_CAST(size_t, aclLen))) {    sprintf(err, "%s: ", me);    biffAdd(NRRD, err); return 1;  }  acl = (unsigned short *)nacl->data;  pos = _nrrd1DIrregMapDomain(&posLen, NULL, nmap);  if (!pos) {    sprintf(err, "%s: couldn't determine domain", me);     biffAdd(NRRD, err); return 1;  }  nacl->axis[1].min = min = pos[0];  nacl->axis[1].max = max = pos[posLen-1];  for (ii=0; ii<=aclLen-1; ii++) {    lo = AIR_AFFINE(0, ii, aclLen, min, max);    hi = AIR_AFFINE(0, ii+1, aclLen, min, max);    acl[0 + 2*ii] = _nrrd1DIrregFindInterval(pos, lo, 0, posLen-2);    acl[1 + 2*ii] = _nrrd1DIrregFindInterval(pos, hi, 0, posLen-2);  }  free(pos);  return 0;}

/*** UV from Eval*/voidtenGlyphBqdUvEval(double uv[2], const double eval[3]) {  double xx, yy, zz, ax, ay, az, mm;  ax = AIR_ABS(eval[0]);  ay = AIR_ABS(eval[1]);  az = AIR_ABS(eval[2]);  mm = AIR_MAX(ax, AIR_MAX(ay, az));  if (mm==0) { /* do not divide */    uv[0]=uv[1]=0;    return;  }  xx = eval[0]/mm;  yy = eval[1]/mm;  zz = eval[2]/mm;  uv[0] = AIR_AFFINE(-1, yy, 1, 0, 1);  if (xx > -zz) {    uv[1] = AIR_AFFINE(-1, zz, 1, 0, 1) - uv[0] + 1;  } else {    uv[1] = AIR_AFFINE(-1, xx, 1, -1, 0) - uv[0] + 1;  }  return;}

voidmakeSceneBVH(limnCamera *cam, echoRTParm *parm, echoObject **sceneP) {  echoObject *sphere;  int i, N;  float r, g, b;  echoObject *scene;  double time0, time1;    *sceneP = scene = echoObjectNew(echoList);  ELL_3V_SET(cam->from, 9, 6, 0);  ELL_3V_SET(cam->at,   0, 0, 0);  ELL_3V_SET(cam->up,   0, 0, 1);  cam->uRange[0] = -3;  cam->uRange[1] = 3;  cam->vRange[0] = -3;  cam->vRange[1] = 3;  parm->jitterType = echoJitterNone;  parm->numSamples = 1;  parm->imgResU = 500;  parm->imgResV = 500;  parm->aperture = 0.0;  parm->renderLights = AIR_TRUE;  parm->renderBoxes = AIR_FALSE;  parm->seedRand = AIR_FALSE;  parm->maxRecDepth = 10;  parm->shadow = 0.0;  N = 1000000;  airArrayLenSet(LIST(scene)->objArr, N);  for (i=0; i<N; i++) {    sphere = echoObjectNew(echoSphere);    echoSphereSet(sphere,                  4*airDrandMT()-2, 4*airDrandMT()-2, 4*airDrandMT()-2, 0.005);    _dyeHSVtoRGB(&r, &g, &b, AIR_AFFINE(0, i, N, 0.0, 1.0), 1.0, 1.0);    echoMatterPhongSet(sphere, r, g, b, 1.0,                       1.0, 0.0, 0.0, 50);    LIST(scene)->obj[i] = sphere;  }  time0 = airTime();  *sceneP = scene = echoListSplit3(scene, 8);  time1 = airTime();  printf("BVH build time = %g seconds/n", time1 - time0);}

/***  0    1       (2)** texp  N*/intmain(int argc, char *argv[]) {  char *me;  unsigned int ii, NN;  me = argv[0];  if (2 != argc || 1 != sscanf(argv[1], "%u", &NN)) {    fprintf(stderr, "%s: need one uint as argument/n", me);    exit(1);  }  for (ii=0; ii<NN; ii++) {    double xx;    xx = AIR_AFFINE(0.0, airDrandMT(), 1.0, -10, 10);    printf("%f %f/n", exp(xx), airFastExp(xx));  }  exit(0);}

voidwheelGraph(wheelPS *wps, double a, double d, double f) {  double A, B, C;  int xi;  double x, y;    A = -a - d - f;  B = a*d + a*f + d*f;  C = -a*d*f;  for (xi=0; xi<=99; xi++) {    x = AIR_AFFINE(0, xi, 99, wps->bbox[0], wps->bbox[2]);    y = (((x + A)*x + B)*x + C)/2;    fprintf(wps->file, "%g %g %s/n", WPS_X(x), WPS_Y(wps->yscale*y),            xi ? "L" : "M");  }  fprintf(wps->file, "S/n");  return;}

intnrrdArithAffine(Nrrd *nout, double minIn,                const Nrrd *nin, double maxIn,                double minOut, double maxOut, int clamp) {  static const char me[]="nrrdArithAffine";  size_t I, N;  double (*ins)(void *v, size_t I, double d),    (*lup)(const void *v, size_t I), mmin, mmax;  if ( !nout || nrrdCheck(nin) ) {    biffAddf(NRRD, "%s: got NULL pointer or invalid input", me);    return 1;  }  if (nout != nin) {    if (nrrdCopy(nout, nin)) {      biffAddf(NRRD, "%s: couldn't initialize output", me);      return 1;    }  }  N = nrrdElementNumber(nin);  ins = nrrdDInsert[nout->type];  lup = nrrdDLookup[nin->type];  mmin = AIR_MIN(minOut, maxOut);  mmax = AIR_MAX(minOut, maxOut);  for (I=0; I<N; I++) {    double val;    val = lup(nin->data, I);    val = AIR_AFFINE(minIn, val, maxIn, minOut, maxOut);    if (clamp) {      val = AIR_CLAMP(mmin, val, mmax);    }    ins(nout->data, I, val);  }  /* HEY: it would be much better if the ordering here was the same as in     AIR_AFFINE, but that's not easy with the way the content functions are     now set up */  if (nrrdContentSet_va(nout, "affine", nin,                        "%g,%g,%g,%g", minIn, maxIn,                        minOut, maxOut)) {    biffAddf(NRRD, "%s:", me);  }  return 0;}

intmain(int argc, char *argv[]) {  char *resS, *scS, *outS;  int i, res;  FILE *out;  float hue, R, G, B, sc;    me = argv[0];  if (4 != argc)    usage();  resS = argv[1];  scS = argv[2];  outS = argv[3];    if (1 != sscanf(resS, "%d", &res)) {    fprintf(stderr, "%s: couldn't parse /"%s/" as int/n", me, resS);    exit(1);  }  if (1 != sscanf(scS, "%f", &sc)) {    fprintf(stderr, "%s: couldn't parse /"%s/" as float/n", me, scS);    exit(1);  }  if (!(out = fopen(outS, "wa"))) {    fprintf(stderr, "%s: couldn't open /"%s/" for writing/n", me, outS);    exit(1);  }  fprintf(out, "# space: RGB/n");  for (i=0; i<=res-1; i++) {    hue = AIR_AFFINE(0, i, res, 0.0, 1.0);    dyeHSVtoRGB(&R, &G, &B, hue, 1, 1);    fprintf(out, "%g %g %g %g/n", hue, sc*R, sc*G, sc*B);  }  fclose(out);  exit(0);}

doublegageStackWtoI(gageContext *ctx, double swrl, int *outside) {  double si;  if (ctx && ctx->parm.stackUse && outside) {    unsigned int sidx;    if (swrl < ctx->stackPos[0]) {      /* we'll extrapolate from stackPos[0] and [1] */      sidx = 0;      *outside = AIR_TRUE;    } else if (swrl > ctx->stackPos[ctx->pvlNum-2]) {      /* extrapolate from stackPos[ctx->pvlNum-3] and [ctx->pvlNum-2];         gageStackPerVolumeAttach ensures that we there are at least two         blurrings pvls & one base pvl ==> pvlNum >= 3 ==> pvlNum-3 >= 0 */      sidx = ctx->pvlNum-3;      *outside = AIR_TRUE;    } else {      /* HEY: stupid linear search */      for (sidx=0; sidx<ctx->pvlNum-2; sidx++) {        if (AIR_IN_CL(ctx->stackPos[sidx], swrl, ctx->stackPos[sidx+1])) {          break;        }      }      if (sidx == ctx->pvlNum-2) {        /* search failure */        *outside = AIR_FALSE;        return AIR_NAN;      }      *outside = AIR_FALSE;    }    si = AIR_AFFINE(ctx->stackPos[sidx], swrl, ctx->stackPos[sidx+1],                    sidx, sidx+1);  } else {    si = AIR_NAN;  }  return si;}

doublegageStackItoW(gageContext *ctx, double si, int *outside) {  unsigned int sidx;  double swrl, sfrac;  if (ctx && ctx->parm.stackUse && outside) {    if (si < 0) {      sidx = 0;      *outside = AIR_TRUE;    } else if (si > ctx->pvlNum-2) {      sidx = ctx->pvlNum-3;      *outside = AIR_TRUE;    } else {      sidx = AIR_CAST(unsigned int, si);      *outside = AIR_FALSE;    }    sfrac = si - sidx;    swrl = AIR_AFFINE(0, sfrac, 1, ctx->stackPos[sidx], ctx->stackPos[sidx+1]);    /*    fprintf(stderr, "!%s: si %g (%u) --> %u + %g --> [%g,%g] -> %g/n", me,            si, ctx->pvlNum, sidx, sfrac,            ctx->stackPos[sidx], ctx->stackPos[sidx+1], swrl);    */  } else {

intmain(int argc, char *argv[]) {    char *me;    double minIn, valIn, maxIn, minOut, maxOut, valOut;    me = argv[0];    if (6 != argc) {        usage(me);    }    if (5 != (sscanf(argv[1], "%lg", &minIn) +              sscanf(argv[2], "%lg", &valIn) +              sscanf(argv[3], "%lg", &maxIn) +              sscanf(argv[4], "%lg", &minOut) +              sscanf(argv[5], "%lg", &maxOut))) {        fprintf(stderr, "%s: couldn't parse all args as doubles/n", me);        usage(me);    }    valOut = AIR_AFFINE(minIn, valIn, maxIn, minOut, maxOut);    printf("%g/n", valOut);    exit(0);}