void
nrrdRangeSafeSet(NrrdRange *range, const Nrrd *nrrd, int blind8BitRange) {
  double minIn, maxIn;

  if (!range) {
    return;
  }
  minIn = range->min;
  maxIn = range->max;
  nrrdRangeSet(range, nrrd, blind8BitRange);
  if (AIR_EXISTS(minIn)) {
    range->min = minIn;
  }
  if (AIR_EXISTS(maxIn)) {
    range->max = maxIn;
  }
  return;
}

/*
** its in here that we scale from "energy gradient" to "force"
*/
double
_pullPointEnergyTotal(pullTask *task, pullBin *bin, pullPoint *point,
                      /* output */
                      double force[4]) {
  char me[]="_pullPointEnergyTotal";
  double enrIm, enrPt, egradIm[4], egradPt[4], energy;
    
  ELL_4V_SET(egradIm, 0, 0, 0, 0); /* sssh */
  ELL_4V_SET(egradPt, 0, 0, 0, 0); /* sssh */
  enrIm = _energyFromImage(task, point, 
                           force ? egradIm : NULL);
  enrPt = _energyFromPoints(task, bin, point,
                            force ? egradPt : NULL);
  energy = AIR_LERP(task->pctx->alpha, enrIm, enrPt);
  /*
  fprintf(stderr, "!%s(%u): energy = lerp(%g, im %g, pt %g) = %g\n", me,
          point->idtag, task->pctx->alpha, enrIm, enrPt, energy);
  */
  if (force) {
    ELL_4V_LERP(force, task->pctx->alpha, egradIm, egradPt);
    ELL_4V_SCALE(force, -1, force);
    /*
    fprintf(stderr, "!%s(%u): egradIm = %g %g %g %g\n", me, point->idtag,
            egradIm[0], egradIm[1], egradIm[2], egradIm[3]);
    fprintf(stderr, "!%s(%u): egradPt = %g %g %g %g\n", me, point->idtag,
            egradPt[0], egradPt[1], egradPt[2], egradPt[3]);
    fprintf(stderr, "!%s(%u): ---> force = %g %g %g %g\n", me,
            point->idtag, force[0], force[1], force[2], force[3]);
    */
  }
  if (task->pctx->wall) {
    unsigned int axi;
    double frc;
    for (axi=0; axi<4; axi++) {
      frc = task->pctx->bboxMin[axi] - point->pos[axi];
      if (frc < 0) {
        /* not below min */
        frc = task->pctx->bboxMax[axi] - point->pos[axi];
        if (frc > 0) {
          /* not above max */
          frc = 0;
        }
      } 
      energy += task->pctx->wall*frc*frc/2;
      if (force) {
        force[axi] += task->pctx->wall*frc;
      }
    }
  }
  if (!AIR_EXISTS(energy)) {
    fprintf(stderr, "!%s(%u): HEY! non-exist energy %g\n", me, point->idtag,
            energy);
  }
  return energy;
}

double
_nrrdApplyDomainMax(const Nrrd *nmap, int ramps, int mapAxis) {
  double ret;
  
  ret = nmap->axis[mapAxis].max;
  if (!AIR_EXISTS(ret)) {
    ret = nmap->axis[mapAxis].size;
    ret = ramps ? ret-1 : ret;
  }
  return ret;
}

void
_nrrdMeasureHistoSum(void *ans, int ansType,
                     const void *line, int lineType, size_t len, 
                     double axmin, double axmax) {
  double sum, hits, val, (*lup)(const void*, size_t);
  size_t ii;
  
  if (!(AIR_EXISTS(axmin) && AIR_EXISTS(axmax))) {
    axmin = -0.5;
    axmax = len-0.5;
  }
  lup = nrrdDLookup[lineType];
  sum = 0;
  for (ii=0; ii<len; ii++) {
    val = NRRD_CELL_POS(axmin, axmax, len, ii);
    hits = lup(line, ii);
    hits = AIR_MAX(hits, 0);
    sum += hits*val;
  }
  nrrdDStore[ansType](ans, sum);
}

int
pullScaleTracePlotAdd(pullContext *pctx, Nrrd *nwild, Nrrd *nccd,
                      Nrrd *nmask, double mth, airArray *insideArr,
                      double velHalf, pullTrace *pts) {
  static const char me[]="pullScaleTracePlotAdd";
  double ssr[2], *pos, *velo, *wild, *ccd, *mask;
  unsigned int pnum, pidx, sizeS, sizeV;

  if (!(pctx && nwild && nccd && pts)) {
    biffAddf(PULL, "%s: got NULL pointer", me);
    return 1;
  }
  if (!nrrdSameSize(nwild, nccd, AIR_TRUE)) {
    biffMovef(PULL, NRRD, "%s: nwild not same size as nccd", me);
    return 1;
  }
  if (nmask || insideArr) {
    if (!insideArr) {
      biffAddf(PULL, "%s: got nmask but not insideArr", me);
      return 1;
    }
    if (!nmask) {
      biffAddf(PULL, "%s: got insideArr but not nmask", me);
      return 1;
    }
    if (nrrdTypeDouble != nmask->type) {
      biffAddf(PULL, "%s: nmask has type %s but want %s", me,
               airEnumStr(nrrdType, nmask->type),
               airEnumStr(nrrdType, nrrdTypeDouble));
      return 1;
    }
    if (!nrrdSameSize(nwild, nmask, AIR_TRUE)) {
      biffMovef(PULL, NRRD, "%s: nwild not same size as nmask", me);
      return 1;
    }
    if (!nrrdSameSize(nccd, nmask, AIR_TRUE)) {
      biffMovef(PULL, NRRD, "%s: nccd not same size as nmask", me);
      return 1;
    }
    if (!AIR_EXISTS(mth)) {
      biffAddf(PULL, "%s: got non-existent mask thresh %g", me, mth);
      return 1;
    }
  }
  ssr[0] = nwild->axis[0].min;
  ssr[1] = nwild->axis[0].max;
  sizeS = AIR_CAST(unsigned int, nwild->axis[0].size);
  sizeV = AIR_CAST(unsigned int, nwild->axis[1].size);
  wild = AIR_CAST(double *, nwild->data);
  ccd = AIR_CAST(double *, nccd->data);
  if (nmask) {
    mask = AIR_CAST(double *, nmask->data);
  } else {

/*
** _nrrdResampleCheckInfo()
**
** checks validity of given NrrdResampleInfo *info: 
** - all required parameters exist
** - both min[d] and max[d] for all axes d
*/
int
_nrrdResampleCheckInfo(const Nrrd *nin, const NrrdResampleInfo *info) {
  char me[] = "_nrrdResampleCheckInfo", err[BIFF_STRLEN];
  const NrrdKernel *k;
  int center, p, np;
  unsigned int ai, minsmp;

  if (nrrdTypeBlock == nin->type || nrrdTypeBlock == info->type) {
    sprintf(err, "%s: can't resample to or from type %s", me,
            airEnumStr(nrrdType, nrrdTypeBlock));
    biffAdd(NRRD, err); return 1;
  }
  if (nrrdBoundaryUnknown == info->boundary) {
    sprintf(err, "%s: didn't set boundary behavior\n", me);
    biffAdd(NRRD, err); return 1;
  }
  if (nrrdBoundaryPad == info->boundary && !AIR_EXISTS(info->padValue)) {
    sprintf(err, "%s: asked for boundary padding, but no pad value set\n", me);
    biffAdd(NRRD, err); return 1;
  }
  for (ai=0; ai<nin->dim; ai++) {
    k = info->kernel[ai];
    /* we only care about the axes being resampled */
    if (!k)
      continue;
    if (!(info->samples[ai] > 0)) {
      sprintf(err, "%s: axis %d # samples (" _AIR_SIZE_T_CNV ") invalid", 
              me, ai, info->samples[ai]);
      biffAdd(NRRD, err); return 1;
    }
    if (!( AIR_EXISTS(nin->axis[ai].min) && AIR_EXISTS(nin->axis[ai].max) )) {
      sprintf(err, "%s: input nrrd's axis %d min,max have not both been set",
              me, ai);
      biffAdd(NRRD, err); return 1;
    }
    if (!( AIR_EXISTS(info->min[ai]) && AIR_EXISTS(info->max[ai]) )) {
      sprintf(err, "%s: info's axis %d min,max not both set", me, ai);
      biffAdd(NRRD, err); return 1;
    }
    np = k->numParm;
    for (p=0; p<np; p++) {
      if (!AIR_EXISTS(info->parm[ai][p])) {
        sprintf(err, "%s: didn't set parameter %d (of %d) for axis %d\n",
                me, p, np, ai);
        biffAdd(NRRD, err); return 1;
      }
    }
    center = _nrrdCenter(nin->axis[ai].center);
    minsmp = nrrdCenterCell == center ? 1 : 2;
    if (!( nin->axis[ai].size >= minsmp && info->samples[ai] >= minsmp )) {
      sprintf(err, "%s: axis %d # input samples (" _AIR_SIZE_T_CNV 
              ") or output samples (" _AIR_SIZE_T_CNV ") "
              " invalid for %s centering",
              me, ai, nin->axis[ai].size, info->samples[ai],
              airEnumStr(nrrdCenter, center));
      biffAdd(NRRD, err); return 1;
    }
  }
  return 0;
}

void
_nrrdMeasureLineError(void *ans, int ansType,
                      const void *line, int lineType, size_t len, 
                      double axmin, double axmax) {
  double x, y, slope, intc, tmp, err=0, (*lup)(const void*, size_t);
  size_t ii;
  
  _nrrdMeasureLineFit(&intc, &slope, line, lineType, len, axmin, axmax);

  if (!( AIR_EXISTS(axmin) && AIR_EXISTS(axmax) )) {
    axmin = 0;
    axmax = len-1;
  }
  lup = nrrdDLookup[lineType];
  for (ii=0; ii<len; ii++) {
    x = NRRD_NODE_POS(axmin, axmax, len, ii);
    y = lup(line, ii);
    tmp = slope*x + intc - y;
    err += tmp*tmp;
  }
  nrrdDStore[ansType](ans, err);
}

void
_nrrdMeasureMean(void *ans, int ansType,
                 const void *line, int lineType, size_t len, 
                 double axmin, double axmax) {
  double val, S, M, (*lup)(const void*, size_t);
  size_t ii, count;

  AIR_UNUSED(axmin);
  AIR_UNUSED(axmax);
  lup = nrrdDLookup[lineType];
  if (nrrdTypeIsIntegral[lineType]) {
    S = 0.0;
    for (ii=0; ii<len; ii++) {
      S += lup(line, ii);
    }
    M = S/len;
  } else {
    S = AIR_NAN;
    for (ii=0; !AIR_EXISTS(S) && ii<len; ii++) {
      S = lup(line, ii);
    }
    if (AIR_EXISTS(S)) {
      /* there was an existant value */
      count = 1;
      for (; ii<len; ii++) {
        val = lup(line, ii);
        if (AIR_EXISTS(val)) {
          count++;
          S += val;
        }
      }
      M = S/count;
    } else {
      /* there were NO existant values */
      M = AIR_NAN;
    }
  }
  nrrdDStore[ansType](ans, M);
}

void
_nrrdMeasureMedian(void *ans, int ansType,
                   const void *_line, int lineType, size_t len, 
                   double axmin, double axmax) {
  double M=0, (*lup)(const void*, size_t);
  size_t ii, mid;
  void *line;

  AIR_UNUSED(axmin);
  AIR_UNUSED(axmax);
  lup = nrrdDLookup[lineType];
  line = calloc(len, nrrdTypeSize[lineType]);
  if (line) {
    memcpy(line, _line, len*nrrdTypeSize[lineType]);
  
    /* yes, I know, this is not the fastest median.  I'll get to it ... */
    qsort(line, len, nrrdTypeSize[lineType], nrrdValCompare[lineType]);
    M = AIR_NAN;
    for (ii=0; !AIR_EXISTS(M) && ii<len; ii++) {
      M = lup(line, ii);
    }
    
    if (AIR_EXISTS(M)) {
      /* i is index AFTER first existant value */
      ii--;
      len -= ii;
      mid = len/2;
      if (len % 2) {
        /* len is odd, there is a middle value, its at mid */
        M = lup(line, ii+mid);
      } else {
        /* len is even, two middle values are at mid-1 and mid */
        M = (lup(line, ii+mid-1) + lup(line, ii+mid))/2;
      }
    }
  }
  nrrdDStore[ansType](ans, M);
}

int
baneInputCheck (Nrrd *nin, baneHVolParm *hvp) {
  char me[]="baneInputCheck", err[BIFF_STRLEN];
  int i;

  if (nrrdCheck(nin)) {
    sprintf(err, "%s: basic nrrd validity check failed", me);
    biffMove(BANE, err, NRRD); return 1;
  }
  if (3 != nin->dim) {
    sprintf(err, "%s: need a 3-dimensional nrrd (not %d)", me, nin->dim);
    biffAdd(BANE, err); return 1;
  }
  if (nrrdTypeBlock == nin->type) {
    sprintf(err, "%s: can't operate on block type", me);
    biffAdd(BANE, err); return 1;
  }
  if (!( AIR_EXISTS(nin->axis[0].spacing) && nin->axis[0].spacing != 0 &&
         AIR_EXISTS(nin->axis[1].spacing) && nin->axis[1].spacing != 0 &&
         AIR_EXISTS(nin->axis[2].spacing) && nin->axis[2].spacing != 0 )) {
    sprintf(err, "%s: must have non-zero existant spacing for all 3 axes", me);
    biffAdd(BANE, err); return 1;
  }
  for (i=0; i<=2; i++) {
    if (_baneAxisCheck(hvp->axis + i)) {
      sprintf(err, "%s: trouble with axis %d", me, i);
      biffAdd(BANE, err); return 1;
    }
  }
  if (!hvp->clip) {
    sprintf(err, "%s: got NULL baneClip", me);
    biffAdd(BANE, err); return 1;
  }

  /* all okay */
  return 0;
}

int
_nrrdSpaceVecExists(const Nrrd *nrrd, unsigned int axi) {
  unsigned int sai;
  int ret;

  if (!( nrrd && axi < nrrd->dim && nrrd->spaceDim )) {
    ret = AIR_FALSE;
  } else {
    ret = AIR_TRUE;
    for (sai=0; sai<nrrd->spaceDim; sai++) {
      ret &= AIR_EXISTS(nrrd->axis[axi].spaceDirection[sai]);
    }
  }
  return ret;
}

void
_nrrdMeasureHistoMin(void *ans, int ansType,
                     const void *line, int lineType, size_t len, 
                     double axmin, double axmax) {
  double val, (*lup)(const void*, size_t);
  size_t ii;

  if (!(AIR_EXISTS(axmin) && AIR_EXISTS(axmax))) {
    axmin = -0.5;
    axmax = len-0.5;
  }
  lup = nrrdDLookup[lineType];
  for (ii=0; ii<len; ii++) {
    if (lup(line, ii) > 0) {
      break;
    }
  }
  if (ii==len) {
    nrrdDStore[ansType](ans, AIR_NAN);
    return;
  }
  val = NRRD_CELL_POS(axmin, axmax, len, ii);
  nrrdDStore[ansType](ans, val);
}

baneClip *
baneClipNew(int type, double *parm) {
  static const char me[]="baneClipNew";
  baneClip *clip;

  if (!( AIR_IN_OP(baneClipUnknown, type, baneClipLast) )) {
    biffAddf(BANE, "%s: baneClip %d invalid", me, type);
    return NULL;
  }
  if (!parm) {
    biffAddf(BANE, "%s: got NULL pointer", me);
    return NULL;
  }
  if (!(AIR_EXISTS(parm[0]))) {
    biffAddf(BANE, "%s: parm[0] doesn't exist", me);
    return NULL;
  }
  clip = (baneClip*)calloc(1, sizeof(baneClip));
  if (!clip) {
    biffAddf(BANE, "%s: couldn't allocate baneClip!", me);
    return NULL;
  }
  clip->parm[0] = parm[0];
  clip->type = type;
  switch(type) {
  case baneClipAbsolute:
    sprintf(clip->name, "absolute");
    clip->answer = _baneClipAnswer_Absolute;
    break;
  case baneClipPeakRatio:
    sprintf(clip->name, "peak ratio");
    clip->answer = _baneClipAnswer_PeakRatio;
    break;
  case baneClipPercentile:
    sprintf(clip->name, "percentile");
    clip->answer = _baneClipAnswer_Percentile;
    break;
  case baneClipTopN:
    sprintf(clip->name, "top N");
    clip->answer = _baneClipAnswer_TopN;
    break;
  default:
    biffAddf(BANE, "%s: sorry, baneClip %d not implemented", me, type);
    baneClipNix(clip); return NULL;
    break;
  }
  return clip;
}

/*
******** nrrdSpacingCalculate
**
** Determine nrrdSpacingStatus, and whatever can be calculated about
** spacing for a given axis.  Takes a nrrd, an axis, a double pointer
** (for returning a scalar), a space vector, and an int pointer for
** returning the known length of the space vector.
**
** The behavior of what has been set by the function is determined by
** the return value, which takes values from the nrrdSpacingStatus*
** enum, as follows:
**
** returned status value:            what it means, and what it set
** ---------------------------------------------------------------------------
** nrrdSpacingStatusUnknown          Something about the given arguments is
**                                   invalid.
**                                   *spacing = NaN,
**                                   vector = all NaNs
**
** nrrdSpacingStatusNone             There is no spacing info at all:
**                                   *spacing = NaN,
**                                   vector = all NaNs
**
** nrrdSpacingStatusScalarNoSpace    There is no surrounding space, but the
**                                   axis's spacing was known.
**                                   *spacing = axis->spacing,
**                                   vector = all NaNs
**
** nrrdSpacingStatusScalarWithSpace  There *is* a surrounding space, but the
**                                   given axis does not live in that space,
**                                   because it has no space direction.  Caller
**                                   may want to think about what's going on.
**                                   *spacing = axis->spacing,
**                                   vector = all NaNs
**
** nrrdSpacingStatusDirection        There is a surrounding space, in which
**                                   this axis has a direction V:
**                                   *spacing = |V| (length of direction),
**                                   vector = V/|V| (normalized direction)
**                                   NOTE: it is still possible for both
**                                   *spacing and vector to be all NaNs!!
*/
int
nrrdSpacingCalculate(const Nrrd *nrrd, unsigned int ax,
                     double *spacing, double vector[NRRD_SPACE_DIM_MAX]) {
  int ret;
  
  if (!( nrrd && spacing && vector
         && ax <= nrrd->dim-1
         && !_nrrdCheck(nrrd, AIR_FALSE, AIR_FALSE) )) {
    /* there's a problem with the arguments.  Note: the _nrrdCheck()
       call does not check on non-NULL-ity of nrrd->data */
    ret = nrrdSpacingStatusUnknown;
    if (spacing) { 
      *spacing = AIR_NAN;
    }
    if (vector) {
      nrrdSpaceVecSetNaN(vector);
    }
  } else {
    if (AIR_EXISTS(nrrd->axis[ax].spacing)) {
      if (nrrd->spaceDim > 0) {
        ret = nrrdSpacingStatusScalarWithSpace;
      } else {
        ret = nrrdSpacingStatusScalarNoSpace;
      }
      *spacing = nrrd->axis[ax].spacing;
      nrrdSpaceVecSetNaN(vector);      
    } else {
      if (nrrd->spaceDim > 0 && _nrrdSpaceVecExists(nrrd, ax)) {
        ret = nrrdSpacingStatusDirection;
        *spacing = nrrdSpaceVecNorm(nrrd->spaceDim, 
                                    nrrd->axis[ax].spaceDirection);
        nrrdSpaceVecScale(vector, 1.0/(*spacing),
                          nrrd->axis[ax].spaceDirection);
      } else {
        ret = nrrdSpacingStatusNone;
        *spacing = AIR_NAN;
        nrrdSpaceVecSetNaN(vector);
      }
    }
  }
  return ret;
}

dyeColor *
dyeColorSet(dyeColor *col, int space, float v0, float v1, float v2) {

  if (col && DYE_VALID_SPACE(space)) {
    col->ii = AIR_CLAMP(0, col->ii, 1);

    /* We switch to the other one if the current one seems to be used,
       but we don't switch if new and current colorspaces are the same.
       If the other one is being used too, oh well.  */
    if (dyeSpaceUnknown != col->spc[col->ii] &&
        AIR_EXISTS(col->val[col->ii][0]) &&
        col->spc[col->ii] != space) {
      col->ii = 1 - col->ii;
    }

    ELL_3V_SET(col->val[col->ii], v0, v1, v2);
    col->spc[col->ii] = space;
  }
  return col;
}

void
_nrrdStrcatSpaceVector(char *str, unsigned int spaceDim,
                       const double val[NRRD_SPACE_DIM_MAX]) {
  char buff[AIR_STRLEN_MED];  /* bad Gordon */
  unsigned int dd;

  if (AIR_EXISTS(val[0])) {
    strcat(str, "(");
    for (dd=0; dd<spaceDim; dd++) {
      strcpy(buff, "");
      airSinglePrintf(NULL, buff, "%.17g", val[dd]);
      strcat(str, buff);
      sprintf(buff, "%s", dd+1 < spaceDim ? "," : ")");
      strcat(str, buff);
    }
  } else {
    strcat(str, _nrrdNoSpaceVector);
  }
  return;
}

void
_nrrdMeasureVariance(void *ans, int ansType,
                     const void *line, int lineType, size_t len, 
                     double axmin, double axmax) {
  double val, S, SS, (*lup)(const void*, size_t);
  size_t ii, count;

  AIR_UNUSED(axmin);
  AIR_UNUSED(axmax);
  SS = S = 0.0;
  lup = nrrdDLookup[lineType];
  if (nrrdTypeIsIntegral[lineType]) {
    for (ii=0; ii<len; ii++) {
      val = lup(line, ii);
      S += val;
      SS += val*val;
    }
    S /= len;
    SS /= len;
  } else {
    count = 0;
    for (ii=0; ii<len; ii++) {
      val = lup(line, ii);
      if (AIR_EXISTS(val)) {
        count++;
        S += val;
        SS += val*val;
      }
    }
    if (count) {
      S /= count;
      SS /= count;
    } else {
      S = SS = AIR_NAN;
    }
  }
  nrrdDStore[ansType](ans, SS - S*S);
}

void
nrrdAxisInfoMinMaxSet(Nrrd *nrrd, unsigned int ax, int defCenter) {
  int center;
  double spacing;

  if (!( nrrd && ax <= nrrd->dim-1 )) {
    return;
  }
  
  center = _nrrdCenter2(nrrd->axis[ax].center, defCenter);
  spacing = nrrd->axis[ax].spacing;
  if (!AIR_EXISTS(spacing))
    spacing = nrrdDefaultSpacing;
  if (nrrdCenterCell == center) {
    nrrd->axis[ax].min = 0;
    nrrd->axis[ax].max = spacing*nrrd->axis[ax].size;
  } else {
    nrrd->axis[ax].min = 0;
    nrrd->axis[ax].max = spacing*(nrrd->axis[ax].size - 1);
  }
  
  return;
}

int
baneHVolCheck (Nrrd *hvol) {
  char me[]="baneHVolCheck", err[BIFF_STRLEN];

  if (3 != hvol->dim) {
    sprintf(err, "%s: need dimension to be 3 (not %d)", me, hvol->dim);
    biffAdd(BANE, err); return 1;
  }
  if (nrrdTypeUChar != hvol->type) {
    sprintf(err, "%s: need type to be %s (not %s)", 
            me, airEnumStr(nrrdType, nrrdTypeUChar),
            airEnumStr(nrrdType, hvol->type));
    biffAdd(BANE, err); return 1;
  }
  if (!( AIR_EXISTS(hvol->axis[0].min) && AIR_EXISTS(hvol->axis[0].max) && 
         AIR_EXISTS(hvol->axis[1].min) && AIR_EXISTS(hvol->axis[1].max) && 
         AIR_EXISTS(hvol->axis[2].min) && AIR_EXISTS(hvol->axis[2].max) )) {
    sprintf(err, "%s: axisMin and axisMax must be set for all axes", me);
    biffAdd(BANE, err); return 1;
  }
  /* 
  ** NOTE: For the time being, I'm giving up on enforcing a 
  ** particular kind of histogram volume ...
  if (strcmp(hvol->axis[0].label, baneMeasrGradMag->name)) {
    sprintf(err, "%s: expected \"%s\" on axis 0 label",
            me, baneMeasrGradMag->name);
    biffAdd(BANE, err); return 1;
  }
  if (strcmp(hvol->axis[1].label, baneMeasrLapl->name) &&
      strcmp(hvol->axis[1].label, baneMeasrHess->name)) {
    sprintf(err, "%s: expected a 2nd deriv. measr on axis 1 (%s or %s)",
            me, baneMeasrHess->name, baneMeasrLapl->name);
    biffAdd(BANE, err); return 1;    
  }
  if (strcmp(hvol->axis[2].label, baneMeasrVal->name)) {
    sprintf(err, "%s: expected \"%s\" on axis 2",
            me, baneMeasrVal->name);
    biffAdd(BANE, err); return 1;
  }
  */
  return 0;
}

/*
******** nrrdAxisInfoSet_nva()
**
** Simple means of setting fields of the axis array in the nrrd.
**
** type to pass for third argument:
**           nrrdAxisInfoSize: size_t*
**        nrrdAxisInfoSpacing: double*
**      nrrdAxisInfoThickness: double*
**            nrrdAxisInfoMin: double*
**            nrrdAxisInfoMax: double*
** nrrdAxisInfoSpaceDirection: double (*var)[NRRD_SPACE_DIM_MAX]
**         nrrdAxisInfoCenter: int*
**           nrrdAxisInfoKind: int*
**          nrrdAxisInfoLabel: char**
**          nrrdAxisInfoUnits: char**
*/
void
nrrdAxisInfoSet_nva(Nrrd *nrrd, int axInfo, const void *_info) {
  _nrrdAxisInfoSetPtrs info;
  int exists;
  unsigned int ai, si, minsi;
  
  if (!( nrrd 
         && AIR_IN_CL(1, nrrd->dim, NRRD_DIM_MAX) 
         && AIR_IN_OP(nrrdAxisInfoUnknown, axInfo, nrrdAxisInfoLast) 
         && _info )) {
    return;
  }
  info.P = _info;

  for (ai=0; ai<nrrd->dim; ai++) {
    switch (axInfo) {
    case nrrdAxisInfoSize:
      nrrd->axis[ai].size = info.ST[ai];
      break;
    case nrrdAxisInfoSpacing:
      nrrd->axis[ai].spacing = info.D[ai];
      break;
    case nrrdAxisInfoThickness:
      nrrd->axis[ai].thickness = info.D[ai];
      break;
    case nrrdAxisInfoMin:
      nrrd->axis[ai].min = info.D[ai];
      break;
    case nrrdAxisInfoMax:
      nrrd->axis[ai].max = info.D[ai];
      break;
    case nrrdAxisInfoSpaceDirection:
      /* we won't allow setting an invalid direction */
      exists = AIR_EXISTS(info.V[ai][0]);
      minsi = nrrd->spaceDim;
      for (si=0; si<nrrd->spaceDim; si++) {
        nrrd->axis[ai].spaceDirection[si] = info.V[ai][si];
        if (exists ^ AIR_EXISTS(info.V[ai][si])) {
          minsi = 0;
          break;
        }
      }
      for (si=minsi; si<NRRD_SPACE_DIM_MAX; si++) {
        nrrd->axis[ai].spaceDirection[si] = AIR_NAN;
      }
      break;
    case nrrdAxisInfoCenter:
      nrrd->axis[ai].center = info.I[ai];
      break;
    case nrrdAxisInfoKind:
      nrrd->axis[ai].kind = info.I[ai];
      break;
    case nrrdAxisInfoLabel:
      nrrd->axis[ai].label = (char *)airFree(nrrd->axis[ai].label);
      nrrd->axis[ai].label = (char *)airStrdup(info.CP[ai]);
      break;
    case nrrdAxisInfoUnits:
      nrrd->axis[ai].units = (char *)airFree(nrrd->axis[ai].units);
      nrrd->axis[ai].units = (char *)airStrdup(info.CP[ai]);
      break;
    }
  }
  if (nrrdAxisInfoSpaceDirection == axInfo) {
    for (ai=nrrd->dim; ai<NRRD_DIM_MAX; ai++) {
      for (si=0; si<NRRD_SPACE_DIM_MAX; si++) {
        nrrd->axis[ai].spaceDirection[si] = AIR_NAN;
      }
    }    
  }
  return;
}