static void test_headers(void)
{
    /* only test the interface; the core parser is tested by the tests above */

    struct phr_header headers[4];
    size_t num_headers;

#define PARSE(s, last_len, exp, comment)                                                                                           \
    do {                                                                                                                           \
        note(comment);                                                                                                             \
        num_headers = sizeof(headers) / sizeof(headers[0]);                                                                        \
        ok(phr_parse_headers(s, strlen(s), headers, &num_headers, last_len) == (exp == 0 ? strlen(s) : exp));                      \
    } while (0)

    PARSE("Host: example.com\r\nCookie: \r\n\r\n", 0, 0, "simple");
    ok(num_headers == 2);
    ok(bufis(headers[0].name, headers[0].name_len, "Host"));
    ok(bufis(headers[0].value, headers[0].value_len, "example.com"));
    ok(bufis(headers[1].name, headers[1].name_len, "Cookie"));
    ok(bufis(headers[1].value, headers[1].value_len, ""));

    PARSE("Host: example.com\r\nCookie: \r\n\r\n", 1, 0, "slowloris");
    ok(num_headers == 2);
    ok(bufis(headers[0].name, headers[0].name_len, "Host"));
    ok(bufis(headers[0].value, headers[0].value_len, "example.com"));
    ok(bufis(headers[1].name, headers[1].name_len, "Cookie"));
    ok(bufis(headers[1].value, headers[1].value_len, ""));

    PARSE("Host: example.com\r\nCookie: \r\n\r", 0, -2, "partial");

    PARSE("Host: e\7fample.com\r\nCookie: \r\n\r", 0, -1, "error");

#undef PARSE
}

int
unrrdu_flipMain(int argc, char **argv, char *me, hestParm *hparm) {
  hestOpt *opt = NULL;
  char *out, *err;
  Nrrd *nin, *nout;
  int pret;
  unsigned int axis;
  airArray *mop;

  OPT_ADD_AXIS(axis, "axis to flip along");
  OPT_ADD_NIN(nin, "input nrrd");
  OPT_ADD_NOUT(out, "output nrrd");

  mop = airMopNew();
  airMopAdd(mop, opt, (airMopper)hestOptFree, airMopAlways);

  USAGE(_unrrdu_flipInfoL);
  PARSE();
  airMopAdd(mop, opt, (airMopper)hestParseFree, airMopAlways);

  nout = nrrdNew();
  airMopAdd(mop, nout, (airMopper)nrrdNuke, airMopAlways);

  if (nrrdFlip(nout, nin, axis)) {
    airMopAdd(mop, err = biffGetDone(NRRD), airFree, airMopAlways);
    fprintf(stderr, "%s: error flipping nrrd:\n%s", me, err);
    airMopError(mop);
    return 1;
  }

  SAVE(out, nout, NULL);

  airMopOkay(mop);
  return 0;
}

/******************************************************************************
 *                                                                            *
 * Comments: Translate device name to the one used internally by kernel. The  *
 *           translation is done based on minor and major device numbers      *
 *           listed in INFO_FILE_NAME . If the names differ it is usually an  *
 *           LVM device which is listed in kernel device mapper.              *
 *                                                                            *
 ******************************************************************************/
static int	get_kernel_devname(const char *devname, char *kernel_devname, size_t max_kernel_devname_len)
{
	FILE		*f;
	char		tmp[MAX_STRING_LEN], name[MAX_STRING_LEN], dev_path[MAX_STRING_LEN];
	int		ret = FAIL;
	zbx_uint64_t	ds[ZBX_DSTAT_MAX], rdev_major, rdev_minor;
	zbx_stat_t	dev_st;

	if ('\0' == *devname)
		return ret;

	*dev_path = '\0';
	if (0 != strncmp(devname, ZBX_DEV_PFX, ZBX_CONST_STRLEN(ZBX_DEV_PFX)))
		strscpy(dev_path, ZBX_DEV_PFX);
	strscat(dev_path, devname);

	if (zbx_stat(dev_path, &dev_st) < 0 || NULL == (f = fopen(INFO_FILE_NAME, "r")))
		return ret;

	while (NULL != fgets(tmp, sizeof(tmp), f))
	{
		PARSE(tmp);
		if (major(dev_st.st_rdev) != rdev_major || minor(dev_st.st_rdev) != rdev_minor)
			continue;

		zbx_strlcpy(kernel_devname, name, max_kernel_devname_len);
		ret = SUCCEED;
		break;
	}
	zbx_fclose(f);

	return ret;
}

static int get_disk_stat(const char *interface, struct disk_stat_s *result)
{
	int ret = SYSINFO_RET_FAIL;
	char line[MAX_STRING_LEN];

	char name[MAX_STRING_LEN];
	
	FILE *f;

	assert(result);

	if(NULL != (f = fopen(INFO_FILE_NAME,"r") ))
	{
		while(fgets(line,MAX_STRING_LEN,f) != NULL)
		{
			PARSE(line);
		
			if(strncmp(name, interface, MAX_STRING_LEN) == 0)
			{
				ret = SYSINFO_RET_OK;
				break;
			}
		}
		zbx_fclose(f);
	}

	if(ret != SYSINFO_RET_OK)
	{
		memset(result, 0, sizeof(struct disk_stat_s));
	}
	
	return ret;
}

int
unrrdu_axdeleteMain(int argc, char **argv, char *me, hestParm *hparm) {
  hestOpt *opt = NULL;
  char *out, *err;
  Nrrd *nin, *nout, *ntmp;
  int pret, _axis;
  unsigned axis;
  airArray *mop;

  hestOptAdd(&opt, "a,axis", "axis", airTypeInt, 1, 1, &_axis, NULL, 
             "dimension (axis index) of the axis to remove");
  OPT_ADD_NIN(nin, "input nrrd");
  OPT_ADD_NOUT(out, "output nrrd");

  mop = airMopNew();
  airMopAdd(mop, opt, (airMopper)hestOptFree, airMopAlways);

  USAGE(_unrrdu_axdeleteInfoL);
  PARSE();
  airMopAdd(mop, opt, (airMopper)hestParseFree, airMopAlways);

  nout = nrrdNew();
  airMopAdd(mop, nout, (airMopper)nrrdNuke, airMopAlways);

  if (-1 == _axis) {
    ntmp = nrrdNew();
    airMopAdd(mop, ntmp, (airMopper)nrrdNuke, airMopAlways);
    if (nrrdCopy(nout, nin)) {
      airMopAdd(mop, err = biffGetDone(NRRD), airFree, airMopAlways);
      fprintf(stderr, "%s: error copying axis:\n%s", me, err);
      airMopError(mop); return 1;
    }
    for (axis=0;
         axis<nout->dim && nout->axis[axis].size > 1;
         axis++);
    while (axis<nout->dim) {
      if (nrrdAxesDelete(ntmp, nout, axis)
          || nrrdCopy(nout, ntmp)) {
        airMopAdd(mop, err = biffGetDone(NRRD), airFree, airMopAlways);
        fprintf(stderr, "%s: error deleting axis:\n%s", me, err);
        airMopError(mop); return 1;
      }
      for (axis=0;
           axis<nout->dim && nout->axis[axis].size > 1;
           axis++);
    }
  } else {
    if (nrrdAxesDelete(nout, nin, _axis)) {
      airMopAdd(mop, err = biffGetDone(NRRD), airFree, airMopAlways);
      fprintf(stderr, "%s: error deleting axis:\n%s", me, err);
      airMopError(mop); return 1;
    }
  }

  SAVE(out, nout, NULL);

  airMopOkay(mop);
  return 0;
}

int
unrrdu_shuffleMain(int argc, const char **argv, const char *me,
                   hestParm *hparm) {
  hestOpt *opt = NULL;
  char *out, *err;
  Nrrd *nin, *nout;
  unsigned int di, axis, permLen, *perm, *iperm, *whichperm;
  size_t *realperm;
  int inverse, pret;
  airArray *mop;

  /* so that long permutations can be read from file */
  hparm->respFileEnable = AIR_TRUE;

  hestOptAdd(&opt, "p,permute", "slc0 slc1", airTypeUInt, 1, -1, &perm, NULL,
             "new slice ordering", &permLen);
  hestOptAdd(&opt, "inv,inverse", NULL, airTypeInt, 0, 0, &inverse, NULL,
             "use inverse of given permutation");
  OPT_ADD_AXIS(axis, "axis to shuffle along");
  OPT_ADD_NIN(nin, "input nrrd");
  OPT_ADD_NOUT(out, "output nrrd");

  mop = airMopNew();
  airMopAdd(mop, opt, (airMopper)hestOptFree, airMopAlways);

  USAGE(_unrrdu_shuffleInfoL);
  PARSE();
  airMopAdd(mop, opt, (airMopper)hestParseFree, airMopAlways);

  nout = nrrdNew();
  airMopAdd(mop, nout, (airMopper)nrrdNuke, airMopAlways);

  /* we have to do error checking on axis in order to do error
     checking on length of permutation */
  if (!( axis < nin->dim )) {
    fprintf(stderr, "%s: axis %d not in valid range [0,%d]\n",
            me, axis, nin->dim-1);
    airMopError(mop);
    return 1;
  }
  if (!( permLen == nin->axis[axis].size )) {
    char stmp[AIR_STRLEN_SMALL];
    fprintf(stderr, "%s: permutation length (%u) != axis %d's size (%s)\n",
            me, permLen, axis,
            airSprintSize_t(stmp, nin->axis[axis].size));
    airMopError(mop);
    return 1;
  }
  if (inverse) {
    iperm = AIR_CALLOC(permLen, unsigned int);
    airMopAdd(mop, iperm, airFree, airMopAlways);
    if (nrrdInvertPerm(iperm, perm, permLen)) {
      fprintf(stderr,
              "%s: couldn't compute inverse of given permutation\n", me);
      airMopError(mop);
      return 1;
    }
    whichperm = iperm;
  } else {

int
tend_evecrgbMain(int argc, char **argv, char *me, hestParm *hparm) {
  int pret;
  hestOpt *hopt = NULL;
  char *perr, *err;
  airArray *mop;

  tenEvecRGBParm *rgbp;
  Nrrd *nin, *nout;
  char *outS;

  mop = airMopNew();
  airMopAdd(mop, hopt, (airMopper)hestOptFree, airMopAlways);

  rgbp = tenEvecRGBParmNew();
  airMopAdd(mop, rgbp, AIR_CAST(airMopper, tenEvecRGBParmNix), airMopAlways);
  
  hestOptAdd(&hopt, "c", "evec index", airTypeUInt, 1, 1, &(rgbp->which), NULL,
             "which eigenvector will be colored. \"0\" for the "
             "principal, \"1\" for the middle, \"2\" for the minor");
  hestOptAdd(&hopt, "a", "aniso", airTypeEnum, 1, 1, &(rgbp->aniso), NULL,
             "Which anisotropy to use for modulating the saturation "
             "of the colors.  " TEN_ANISO_DESC,
             NULL, tenAniso);
  hestOptAdd(&hopt, "t", "thresh", airTypeDouble, 1, 1, &(rgbp->confThresh),
             "0.5", "confidence threshold");
  hestOptAdd(&hopt, "bg", "background", airTypeDouble, 1, 1, &(rgbp->bgGray),
             "0", "gray level to use for voxels who's confidence is zero ");
  hestOptAdd(&hopt, "gr", "gray", airTypeDouble, 1, 1, &(rgbp->isoGray), "0",
             "the gray level to desaturate towards as anisotropy "
             "decreases (while confidence remains 1.0)");
  hestOptAdd(&hopt, "gam", "gamma", airTypeDouble, 1, 1, &(rgbp->gamma), "1",
             "gamma to use on color components");
  hestOptAdd(&hopt, "i", "nin", airTypeOther, 1, 1, &nin, "-",
             "input diffusion tensor volume", NULL, NULL, nrrdHestNrrd);
  hestOptAdd(&hopt, "o", "nout", airTypeString, 1, 1, &outS, "-",
             "output image (floating point)");

  USAGE(_tend_evecrgbInfoL);
  PARSE();
  airMopAdd(mop, hopt, (airMopper)hestParseFree, airMopAlways);

  nout = nrrdNew();
  airMopAdd(mop, nout, (airMopper)nrrdNuke, airMopAlways);
  if (tenEvecRGB(nout, nin, rgbp)) {
    airMopAdd(mop, err=biffGetDone(TEN), airFree, airMopAlways);
    fprintf(stderr, "%s: trouble doing colormapping:\n%s\n", me, err);
    airMopError(mop); return 1;
  }
  if (nrrdSave(outS, nout, NULL)) {
    airMopAdd(mop, err=biffGetDone(NRRD), airFree, airMopAlways);
    fprintf(stderr, "%s: trouble writing:\n%s\n", me, err);
    airMopError(mop); return 1;
  }

  airMopOkay(mop);
  return 0;
}

int
unrrdu_ccfindMain(int argc, char **argv, char *me, hestParm *hparm) {
  hestOpt *opt = NULL;
  char *out, *err, *valS;
  Nrrd *nin, *nout, *nval=NULL;
  airArray *mop;
  int type, pret;
  unsigned int conny;

  hestOptAdd(&opt, "v,values", "filename", airTypeString, 1, 1, &valS, "",
             "Giving a filename here allows you to save out the values "
             "associated with each connect component.  This can be used "
             "later with \"ccmerge -d\".  By default, no record of the "
             "original CC values is kept.");
  hestOptAdd(&opt, "t,type", "type", airTypeOther, 1, 1, &type, "default",
             "type to use for output, to store the CC ID values.  By default "
             "(not using this option), the type used will be the smallest of "
             "uchar, ushort, or int, that can represent all the CC ID values. "
             "Using this option allows one to specify the integral type to "
             "be used.",
             NULL, NULL, &unrrduHestMaybeTypeCB);
  hestOptAdd(&opt, "c,connect", "connectivity", airTypeUInt, 1, 1,
             &conny, NULL,
             "what kind of connectivity to use: the number of coordinates "
             "that vary in order to traverse the neighborhood of a given "
             "sample.  In 2D: \"1\": 4-connected, \"2\": 8-connected");
  OPT_ADD_NIN(nin, "input nrrd");
  OPT_ADD_NOUT(out, "output nrrd");

  mop = airMopNew();
  airMopAdd(mop, opt, (airMopper)hestOptFree, airMopAlways);

  USAGE(_unrrdu_ccfindInfoL);
  PARSE();
  airMopAdd(mop, opt, (airMopper)hestParseFree, airMopAlways);

  nout = nrrdNew();
  airMopAdd(mop, nout, (airMopper)nrrdNuke, airMopAlways);

  if (nrrdCCFind(nout, airStrlen(valS) ? &nval : NULL, nin, type, conny)) {
    airMopAdd(mop, err = biffGetDone(NRRD), airFree, airMopAlways);
    fprintf(stderr, "%s: error doing connected components:\n%s", me, err);
    airMopError(mop);
    return 1;
  }
  if (nval) {
    airMopAdd(mop, nval, (airMopper)nrrdNuke, airMopAlways);
  }

  if (airStrlen(valS)) {
    SAVE(valS, nval, NULL);
  }
  SAVE(out, nout, NULL);

  airMopOkay(mop);
  return 0;
}

int
unrrdu_dataMain(int argc, const char **argv, const char *me,
                hestParm *hparm) {
  hestOpt *opt = NULL;
  char *err, *inS=NULL;
  Nrrd *nin;
  NrrdIoState *nio;
  airArray *mop;
  int car, pret;

  mop = airMopNew();
  hestOptAdd(&opt, NULL, "nin", airTypeString, 1, 1, &inS, NULL,
             "input nrrd");
  airMopAdd(mop, opt, (airMopper)hestOptFree, airMopAlways);

  USAGE(_unrrdu_dataInfoL);
  PARSE();
  airMopAdd(mop, opt, (airMopper)hestParseFree, airMopAlways);

  nio = nrrdIoStateNew();
  airMopAdd(mop, nio, (airMopper)nrrdIoStateNix, airMopAlways);
  nio->skipData = AIR_TRUE;
  nio->keepNrrdDataFileOpen = AIR_TRUE;
  nin = nrrdNew();
  airMopAdd(mop, nin, (airMopper)nrrdNuke, airMopAlways);

  if (nrrdLoad(nin, inS, nio)) {
    airMopAdd(mop, err = biffGetDone(NRRD), airFree, airMopAlways);
    fprintf(stderr, "%s: error reading header:\n%s", me, err);
    airMopError(mop);
    return 1;
  }
  if (_nrrdDataFNNumber(nio) > 1) {
    fprintf(stderr, "%s: sorry, currently can't operate with multiple "
            "detached datafiles\n", me);
    airMopError(mop);
    return 1;
  }
  if (!( nrrdFormatNRRD == nio->format )) {
    fprintf(stderr, "%s: can only print data of NRRD format files\n", me);
    airMopError(mop); return 1;
  }
  car = fgetc(nio->dataFile);
#ifdef _MSC_VER
  /* needed because otherwise printing a carraige return will
     automatically also produce a newline */
  _setmode(_fileno(stdout), _O_BINARY);
#endif
  while (EOF != car) {
    fputc(car, stdout);
    car = fgetc(nio->dataFile);
  }
  airFclose(nio->dataFile);

  airMopOkay(mop);
  return 0;
}

int	get_diskstat(const char *devname, zbx_uint64_t *dstat)
{
	FILE		*f;
	char		tmp[MAX_STRING_LEN], name[MAX_STRING_LEN], dev_path[MAX_STRING_LEN];
	int		i, ret = FAIL, dev_exists = FAIL;
	zbx_uint64_t	ds[ZBX_DSTAT_MAX], rdev_major, rdev_minor;
	zbx_stat_t 	dev_st;
	int		found = 0;

	for (i = 0; i < ZBX_DSTAT_MAX; i++)
		dstat[i] = (zbx_uint64_t)__UINT64_C(0);

	if (NULL != devname && '\0' != *devname && 0 != strcmp(devname, "all"))
	{
		*dev_path = '\0';
		if (0 != strncmp(devname, ZBX_DEV_PFX, ZBX_CONST_STRLEN(ZBX_DEV_PFX)))
			strscpy(dev_path, ZBX_DEV_PFX);
		strscat(dev_path, devname);

		if (zbx_stat(dev_path, &dev_st) == 0)
			dev_exists = SUCCEED;
	}

	if (NULL == (f = fopen(INFO_FILE_NAME, "r")))
		return FAIL;

	while (NULL != fgets(tmp, sizeof(tmp), f))
	{
		PARSE(tmp);

		if (NULL != devname && '\0' != *devname && 0 != strcmp(devname, "all"))
		{
			if (0 != strcmp(name, devname))
			{
				if (SUCCEED != dev_exists
					|| major(dev_st.st_rdev) != rdev_major
					|| minor(dev_st.st_rdev) != rdev_minor)
					continue;
			}
			else
				found = 1;
		}

		dstat[ZBX_DSTAT_R_OPER] += ds[ZBX_DSTAT_R_OPER];
		dstat[ZBX_DSTAT_R_SECT] += ds[ZBX_DSTAT_R_SECT];
		dstat[ZBX_DSTAT_W_OPER] += ds[ZBX_DSTAT_W_OPER];
		dstat[ZBX_DSTAT_W_SECT] += ds[ZBX_DSTAT_W_SECT];

		ret = SUCCEED;

		if (1 == found)
			break;
	}
	zbx_fclose(f);

	return ret;
}

int
unrrdu_spliceMain(int argc, char **argv, char *me, hestParm *hparm) {
  hestOpt *opt = NULL;
  char *out, *err;
  Nrrd *nin, *nout, *nslice;
  unsigned int axis;
  int pret;
  long int _pos[2];
  size_t pos;
  airArray *mop;

  OPT_ADD_AXIS(axis, "axis to splice along");
  hestOptAdd(&opt, "p,position", "pos", airTypeOther, 1, 1, _pos, NULL,
             "position to splice at:\n "
             "\b\bo <int> gives 0-based index\n "
             "\b\bo M-<int> give index relative "
             "to the last sample on the axis (M == #samples-1).",
             NULL, NULL, &unrrduHestPosCB);
  hestOptAdd(&opt, "s,slice", "nslice", airTypeOther, 1, 1, &(nslice), NULL,
             "slice nrrd.  This the slice to be inserted in \"nin\"",
             NULL, NULL, nrrdHestNrrd);
  OPT_ADD_NIN(nin, "input nrrd.  This the nrrd into which the slice will "
              "be inserted");
  OPT_ADD_NOUT(out, "output nrrd");

  mop = airMopNew();
  airMopAdd(mop, opt, (airMopper)hestOptFree, airMopAlways);

  USAGE(_unrrdu_spliceInfoL);
  PARSE();
  airMopAdd(mop, opt, (airMopper)hestParseFree, airMopAlways);
  if (!( axis < nin->dim )) {
    fprintf(stderr, "%s: axis %u not in range [0,%u]\n", me, axis, nin->dim-1);
    return 1;
  }
  if (_pos[0] == -1) {
    fprintf(stderr, "%s: m+<int> specification format meaningless here\n", me);
    return 1;
  }
  pos = _pos[0]*(nin->axis[axis].size-1) + _pos[1];

  nout = nrrdNew();
  airMopAdd(mop, nout, (airMopper)nrrdNuke, airMopAlways);

  if (nrrdSplice(nout, nin, nslice, axis, pos)) {
    airMopAdd(mop, err = biffGetDone(NRRD), airFree, airMopAlways);
    fprintf(stderr, "%s: error splicing nrrd:\n%s", me, err);
    airMopError(mop);
    return 1;
  }

  SAVE(out, nout, NULL);

  airMopOkay(mop);
  return 0;
}

int
limnpu_rastMain(int argc, const char **argv, const char *me,
                hestParm *hparm) {
  hestOpt *hopt = NULL;
  char *err, *perr;
  airArray *mop;
  int pret;

  limnPolyData *pld;
  double min[3], max[3];
  Nrrd *nout;
  char *out;
  int type;
  size_t size[NRRD_DIM_MAX];

  hestOptAdd(&hopt, "min", "min", airTypeDouble, 3, 3, min, NULL,
             "bottom corner");
  hestOptAdd(&hopt, "max", "max", airTypeDouble, 3, 3, max, NULL,
             "top corner");
  hestOptAdd(&hopt, "s", "size", airTypeSize_t, 3, 3, size, NULL,
             "number of samples along each axis");
  hestOptAdd(&hopt, "t", "type", airTypeEnum, 1, 1, &type, "uchar",
             "type of output nrrd",
             NULL, nrrdType);
  hestOptAdd(&hopt, NULL, "input", airTypeOther, 1, 1, &pld, NULL,
             "input polydata filename",
             NULL, NULL, limnHestPolyDataLMPD);
  hestOptAdd(&hopt, NULL, "output", airTypeString, 1, 1, &out, NULL,
             "output nrrd filename");

  mop = airMopNew();
  airMopAdd(mop, hopt, (airMopper)hestOptFree, airMopAlways);

  USAGE(myinfo);
  PARSE();
  airMopAdd(mop, hopt, (airMopper)hestParseFree, airMopAlways);

  nout = nrrdNew();
  airMopAdd(mop, nout, (airMopper)nrrdNuke, airMopAlways);

  if (limnPolyDataRasterize(nout, pld, min, max, size, type)) {
    airMopAdd(mop, err = biffGetDone(LIMN), airFree, airMopAlways);
    fprintf(stderr, "%s: trouble:%s", me, err);
    airMopError(mop);
    return 1;
  }
  if (nrrdSave(out, nout, NULL)) {
    airMopAdd(mop, err = biffGetDone(NRRD), airFree, airMopAlways);
    fprintf(stderr, "%s: trouble:%s", me, err);
    airMopError(mop);
    return 1;
  }

  airMopOkay(mop);
  return 0;
}

int
unrrdu_histaxMain(int argc, char **argv, char *me, hestParm *hparm) {
  hestOpt *opt = NULL;
  char *out, *err;
  Nrrd *nin, *nout;
  int type, bins, pret, blind8BitRange;
  unsigned int axis;
  double min, max;
  airArray *mop;
  NrrdRange *range;

  OPT_ADD_AXIS(axis, "axis to histogram along");
  hestOptAdd(&opt, "b,bin", "bins", airTypeInt, 1, 1, &bins, NULL,
             "# of bins in histogram");
  OPT_ADD_TYPE(type, "output type", "uchar");
  hestOptAdd(&opt, "min,minimum", "value", airTypeDouble, 1, 1, &min, "nan",
             "Value at low end of histogram. Defaults to lowest value "
             "found in input nrrd.");
  hestOptAdd(&opt, "max,maximum", "value", airTypeDouble, 1, 1, &max, "nan",
             "Value at high end of histogram. Defaults to highest value "
             "found in input nrrd.");
  hestOptAdd(&opt, "blind8", "bool", airTypeBool, 1, 1, &blind8BitRange,
             nrrdStateBlind8BitRange ? "true" : "false",
             "Whether to know the range of 8-bit data blindly "
             "(uchar is always [0,255], signed char is [-128,127]).");
  OPT_ADD_NIN(nin, "input nrrd");
  OPT_ADD_NOUT(out, "output nrrd");

  mop = airMopNew();
  airMopAdd(mop, opt, (airMopper)hestOptFree, airMopAlways);

  USAGE(_unrrdu_histaxInfoL);
  PARSE();
  airMopAdd(mop, opt, (airMopper)hestParseFree, airMopAlways);

  nout = nrrdNew();
  airMopAdd(mop, nout, (airMopper)nrrdNuke, airMopAlways);

  range = nrrdRangeNew(min, max);
  airMopAdd(mop, range, (airMopper)nrrdRangeNix, airMopAlways);
  nrrdRangeSafeSet(range, nin, blind8BitRange);
  if (nrrdHistoAxis(nout, nin, range, axis, bins, type)) {
    airMopAdd(mop, err = biffGetDone(NRRD), airFree, airMopAlways);
    fprintf(stderr, "%s: error doing axis histogramming:\n%s", me, err);
    airMopError(mop);
    return 1;
  }
  

  SAVE(out, nout, NULL);

  airMopOkay(mop);
  return 0;
}

int
unrrdu_dhistoMain(int argc, const char **argv, const char *me,
                  hestParm *hparm) {
  hestOpt *opt = NULL;
  char *out, *err;
  Nrrd *nin, *nout;
  int pret, nolog, notick;
  unsigned int size;
  airArray *mop;
  double max;

  hestOptAdd(&opt, "h,height", "height", airTypeUInt, 1, 1, &size, NULL,
             "height of output image (horizontal size is determined by "
             "number of bins in input histogram).");
  hestOptAdd(&opt, "nolog", NULL, airTypeInt, 0, 0, &nolog, NULL,
             "do not show the log-scaled histogram with decade tick-marks");
  hestOptAdd(&opt, "notick", NULL, airTypeInt, 0, 0, &notick, NULL,
             "do not draw the log decade tick marks");
  hestOptAdd(&opt, "max,maximum", "max # hits", airTypeDouble, 1, 1,
             &max, "nan",
             "constrain the top of the drawn histogram to be at this "
             "number of hits.  This will either scale the drawn histogram "
             "downward or clip its top, depending on whether the given max "
             "is higher or lower than the actual maximum bin count.  By "
             "not using this option (the default), the actual maximum bin "
             "count is used");
  OPT_ADD_NIN(nin, "input nrrd");
  OPT_ADD_NOUT(out, "output nrrd");

  mop = airMopNew();
  airMopAdd(mop, opt, (airMopper)hestOptFree, airMopAlways);

  USAGE(_unrrdu_dhistoInfoL);
  PARSE();
  airMopAdd(mop, opt, (airMopper)hestParseFree, airMopAlways);

  nout = nrrdNew();
  airMopAdd(mop, nout, (airMopper)nrrdNuke, airMopAlways);

  if (nrrdHistoDraw(nout, nin, size,
                    nolog ? AIR_FALSE : (notick ? 2 : AIR_TRUE),
                    max)) {
    airMopAdd(mop, err = biffGetDone(NRRD), airFree, airMopAlways);
    fprintf(stderr, "%s: error drawing histogram nrrd:\n%s", me, err);
    airMopError(mop);
    return 1;
  }

  SAVE(out, nout, NULL);

  airMopOkay(mop);
  return 0;
}

int
unrrdu_axmergeMain(int argc, const char **argv, const char *me,
                   hestParm *hparm) {
  hestOpt *opt = NULL;
  char *out, *err;
  Nrrd *nin, *nout[2];
  int *axes, pret, ni;
  unsigned int ii, jj, axesLen;
  airArray *mop;

  hestOptAdd(&opt, "a,axis", "ax0", airTypeInt, 1, -1, &axes, NULL,
             "axis (or axes) to merge.  Each axis index identified is the "
             "lower of the pair of axes that will be merged.  Saying \"-a 2\" "
             "means to merge axis 2 and axis 3 into axis 2.  If multiple "
             "merges are to be done, the indices listed here are for "
             "the axes prior to any merging.", &axesLen);
  OPT_ADD_NIN(nin, "input nrrd");
  OPT_ADD_NOUT(out, "output nrrd");

  mop = airMopNew();
  airMopAdd(mop, opt, (airMopper)hestOptFree, airMopAlways);

  USAGE(_unrrdu_axmergeInfoL);
  PARSE();
  airMopAdd(mop, opt, (airMopper)hestParseFree, airMopAlways);

  airMopAdd(mop, nout[0]=nrrdNew(), (airMopper)nrrdNuke, airMopAlways);
  airMopAdd(mop, nout[1]=nrrdNew(), (airMopper)nrrdNuke, airMopAlways);

  if (axesLen > 1) {
    /* sort merge axes into ascending order */
    qsort(axes, axesLen, sizeof(*axes), nrrdValCompare[nrrdTypeInt]);
  }

  ni = 0;
  for (ii=0; ii<axesLen; ii++) {
    if (nrrdAxesMerge(nout[ni], !ii ? nin : nout[1-ni], axes[ii])) {
      airMopAdd(mop, err = biffGetDone(NRRD), airFree, airMopAlways);
      fprintf(stderr, "%s: error merging axes:\n%s", me, err);
      airMopError(mop);
      return 1;
    }
    for (jj=ii+1; jj<axesLen; jj++) {
      axes[jj] -= 1;
    }
    ni = 1-ni;
  }

  SAVE(out, nout[1-ni], NULL);

  airMopOkay(mop);
  return 0;
}

int
tend_anhistMain(int argc, char **argv, char *me, hestParm *hparm) {
  int pret;
  hestOpt *hopt = NULL;
  char *perr, *err;
  airArray *mop;

  int version, res, right;
  Nrrd *nin, *nout, *nwght;
  char *outS;

  hestOptAdd(&hopt, "v", "westin version", airTypeInt, 1, 1, &version, "1",
             "Which version of Westin's anisotropy metric triple "
             "to use, either \"1\" or \"2\"");
  hestOptAdd(&hopt, "w", "nweight", airTypeOther, 1, 1, &nwght, "",
             "how to weigh contributions to histogram.  By default "
             "(not using this option), the increment is one bin count per "
             "sample, but by giving a nrrd, the value in the nrrd at the "
             "corresponding location will be the bin count increment ",
             NULL, NULL, nrrdHestNrrd);
  hestOptAdd(&hopt, "r", "res", airTypeInt, 1, 1, &res, NULL,
             "resolution of anisotropy plot");
  hestOptAdd(&hopt, "right", NULL, airTypeInt, 0, 0, &right, NULL,
             "sample a right-triangle-shaped region, instead of "
             "a roughly equilateral triangle. ");
  hestOptAdd(&hopt, "i", "nin", airTypeOther, 1, 1, &nin, "-",
             "input diffusion tensor volume", NULL, NULL, nrrdHestNrrd);
  hestOptAdd(&hopt, "o", "nout", airTypeString, 1, 1, &outS, "-",
             "output image (floating point)");

  mop = airMopNew();
  airMopAdd(mop, hopt, (airMopper)hestOptFree, airMopAlways);
  USAGE(_tend_anhistInfoL);
  PARSE();
  airMopAdd(mop, hopt, (airMopper)hestParseFree, airMopAlways);

  nout = nrrdNew();
  airMopAdd(mop, nout, (airMopper)nrrdNuke, airMopAlways);

  if (tenAnisoHistogram(nout, nin, nwght, right, version, res)) {
    airMopAdd(mop, err=biffGetDone(TEN), airFree, airMopAlways);
    fprintf(stderr, "%s: trouble making histogram:\n%s\n", me, err);
    airMopError(mop); return 1;
  }
  if (nrrdSave(outS, nout, NULL)) {
    airMopAdd(mop, err=biffGetDone(NRRD), airFree, airMopAlways);
    fprintf(stderr, "%s: trouble writing:\n%s\n", me, err);
    airMopError(mop); return 1;
  }

  airMopOkay(mop);
  return 0;
}

int
tend_evqMain(int argc, char **argv, char *me, hestParm *hparm) {
  int pret;
  hestOpt *hopt = NULL;
  char *perr, *err;
  airArray *mop;

  int which, aniso, dontScaleByAniso;
  Nrrd *nin, *nout;
  char *outS;

  hestOptAdd(&hopt, "c", "evec index", airTypeInt, 1, 1, &which, "0",
             "Which eigenvector should be quantized: \"0\" for the "
             "direction of fastest diffusion (eigenvector associated "
             "with largest eigenvalue), \"1\" or \"2\" for other two "
             "eigenvectors (associated with middle and smallest eigenvalue)");
  hestOptAdd(&hopt, "a", "aniso", airTypeEnum, 1, 1, &aniso, NULL,
             "Which anisotropy metric to scale the eigenvector "
             "with.  " TEN_ANISO_DESC,
             NULL, tenAniso);
  hestOptAdd(&hopt, "ns", NULL, airTypeInt, 0, 0, &dontScaleByAniso, NULL,
             "Don't attenuate the color by anisotropy.  By default (not "
             "using this option), regions with low or no anisotropy are "
             "very dark colors or black");
  hestOptAdd(&hopt, "i", "nin", airTypeOther, 1, 1, &nin, "-",
             "input diffusion tensor volume", NULL, NULL, nrrdHestNrrd);
  hestOptAdd(&hopt, "o", "nout", airTypeString, 1, 1, &outS, "-",
             "output image (floating point)");

  mop = airMopNew();
  airMopAdd(mop, hopt, (airMopper)hestOptFree, airMopAlways);
  USAGE(_tend_evqInfoL);
  PARSE();
  airMopAdd(mop, hopt, (airMopper)hestParseFree, airMopAlways);

  nout = nrrdNew();
  airMopAdd(mop, nout, (airMopper)nrrdNuke, airMopAlways);
  if (tenEvqVolume(nout, nin, which, aniso, !dontScaleByAniso)) {
    airMopAdd(mop, err=biffGetDone(TEN), airFree, airMopAlways);
    fprintf(stderr, "%s: trouble quantizing eigenvectors:\n%s\n", me, err);
    airMopError(mop); return 1;
  }

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

  airMopOkay(mop);
  return 0;
}

int
unrrdu_gammaMain(int argc, const char **argv, const char *me,
                 hestParm *hparm) {
  hestOpt *opt = NULL;
  char *out, *err;
  Nrrd *nin, *nout;
  double min, max, gamma;
  airArray *mop;
  int pret, blind8BitRange;
  NrrdRange *range;

  hestOptAdd(&opt, "g,gamma", "gamma", airTypeDouble, 1, 1, &gamma, NULL,
             "gamma > 1.0 brightens; gamma < 1.0 darkens. "
             "Negative gammas invert values (like in xv). ");
  hestOptAdd(&opt, "min,minimum", "value", airTypeDouble, 1, 1, &min, "nan",
             "Value to implicitly map to 0.0 prior to calling pow(). "
             "Defaults to lowest value found in input nrrd.");
  hestOptAdd(&opt, "max,maximum", "value", airTypeDouble, 1, 1, &max, "nan",
             "Value to implicitly map to 1.0 prior to calling pow(). "
             "Defaults to highest value found in input nrrd.");
  hestOptAdd(&opt, "blind8", "bool", airTypeBool, 1, 1, &blind8BitRange,
             nrrdStateBlind8BitRange ? "true" : "false",
             "Whether to know the range of 8-bit data blindly "
             "(uchar is always [0,255], signed char is [-128,127]).");
  OPT_ADD_NIN(nin, "input nrrd");
  OPT_ADD_NOUT(out, "output nrrd");

  mop = airMopNew();
  airMopAdd(mop, opt, (airMopper)hestOptFree, airMopAlways);

  USAGE(_unrrdu_gammaInfoL);
  PARSE();
  airMopAdd(mop, opt, (airMopper)hestParseFree, airMopAlways);

  nout = nrrdNew();
  airMopAdd(mop, nout, (airMopper)nrrdNuke, airMopAlways);

  range = nrrdRangeNew(min, max);
  airMopAdd(mop, range, (airMopper)nrrdRangeNix, airMopAlways);
  nrrdRangeSafeSet(range, nin, blind8BitRange);
  if (nrrdArithGamma(nout, nin, range, gamma)) {
    airMopAdd(mop, err = biffGetDone(NRRD), airFree, airMopAlways);
    fprintf(stderr, "%s: error doing gamma:\n%s", me, err);
    airMopError(mop);
    return 1;
  }

  SAVE(out, nout, NULL);

  airMopOkay(mop);
  return 0;
}

int
unrrdu_distMain(int argc, char **argv, char *me, hestParm *hparm) {
  hestOpt *opt = NULL;
  char *out, *err;
  Nrrd *nin, *nout;
  int pret;

  int E, typeOut, invert, sign;
  double thresh;
  airArray *mop;

  hestOptAdd(&opt, "th,thresh", "val", airTypeDouble, 1, 1, &thresh, NULL,
             "threshold value to separate inside from outside");
  hestOptAdd(&opt, "t,type", "type", airTypeEnum, 1, 1, &typeOut, "float",
             "type to save output in", NULL, nrrdType);
  hestOptAdd(&opt, "sgn", NULL, airTypeInt, 0, 0, &sign, NULL,
             "also compute signed (negative) distances inside objects, "
             "instead of leaving them as zero");
  hestOptAdd(&opt, "inv", NULL, airTypeInt, 0, 0, &invert, NULL,
             "values *below* threshold are considered interior to object. "
             "By default (not using this option), values above threshold "
             "are considered interior. ");
  OPT_ADD_NIN(nin, "input nrrd");
  OPT_ADD_NOUT(out, "output nrrd");

  mop = airMopNew();
  airMopAdd(mop, opt, (airMopper)hestOptFree, airMopAlways);

  USAGE(_unrrdu_distInfoL);
  PARSE();
  airMopAdd(mop, opt, (airMopper)hestParseFree, airMopAlways);

  nout = nrrdNew();
  airMopAdd(mop, nout, (airMopper)nrrdNuke, airMopAlways);

  if (sign) {
    E = nrrdDistanceL2Signed(nout, nin, typeOut, NULL, thresh, !invert);
  } else {
    E = nrrdDistanceL2(nout, nin, typeOut, NULL, thresh, !invert);
  }
  if (E) {
    airMopAdd(mop, err = biffGetDone(NRRD), airFree, airMopAlways);
    fprintf(stderr, "%s: error doing distance transform:\n%s", me, err);
    airMopError(mop);
    return 1;
  }

  SAVE(out, nout, NULL);

  airMopOkay(mop);
  return 0;
}

void Mesh::readPlyFile(std::istream& is) {
  uint i, k, nVertices, nFaces;
  MT::String str;
  is >>PARSE("ply") >>PARSE("format") >>str;
  if(str=="ascii") {
    is >>PARSE("1.0");
    is >>PARSE("element vertex") >>nVertices;
    is >>PARSE("property float32 x") >>PARSE("property float32 y") >>PARSE("property float32 z");
    is >>PARSE("property float32 nx") >>PARSE("property float32 ny") >>PARSE("property float32 nz");
    is >>PARSE("element face") >>nFaces;
    is >>PARSE("property list uint8 int32 vertex_indices") >>PARSE("end_header");
    V.resize(nVertices, 3);
    T.resize(nFaces   , 3);
    double nx, ny, nz;
    for(i=0; i<V.d0; i++) {
      is >>V(i, 0) >>V(i, 1) >>V(i, 2) >>nx >>ny >>nz;
    }
    fo