text: /*                           V O L . C
   * BRL-CAD
   *
   * Copyright (c) 1989-2021 United States Government as represented by
   * the U.S. Army Research Laboratory.
   *
   * This library is free software; you can redistribute it and/or
   * modify it under the terms of the GNU Lesser General Public License
   * version 2.1 as published by the Free Software Foundation.
   *
   * This library is distributed in the hope that it will be useful, but
   * WITHOUT ANY WARRANTY; without even the implied warranty of
   * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
   * Lesser General Public License for more details.
   *
   * You should have received a copy of the GNU Lesser General Public
   * License along with this file; see the file named COPYING for more
   * information.
   */
  /** @addtogroup primitives */
  /** @{ */
  /** @file primitives/vol/vol.c
   *
   * Intersect a ray with a 3-D volume.  The volume is described as a
   * concatenation of bw(5) files.
   *
   */
  /** @} */
  
  #include "common.h"
  
  #include <stddef.h>
  #include <ctype.h>
  #include <math.h>
  #include <string.h>
  #include "bio.h"
  
  #include "bu/parallel.h"
  #include "vmath.h"
  #include "rt/db4.h"
  #include "nmg.h"
  #include "rt/geom.h"
  #include "raytrace.h"
  
  #include "../fixpt.h"
  
  
  /*
    NOTES:
    Changed small to small1 for win32 compatibility
  */
  
  
  struct rt_vol_specific {
      struct rt_vol_internal vol_i;
      vect_t vol_xnorm;	/* local +X norm in model coords */
      vect_t vol_ynorm;
      vect_t vol_znorm;
      mat_t vol_mat;	/* model to ideal space */
      vect_t vol_origin;	/* local coords of grid origin (0, 0, 0) for now */
      vect_t vol_large;	/* local coords of XYZ max */
  };
  #define VOL_NULL ((struct rt_vol_specific *)0)
  
  #define VOL_O(m) bu_offsetof(struct rt_vol_internal, m)
  
  const struct bu_structparse rt_vol_parse[] = {
      {"%s", RT_VOL_NAME_LEN, "file", bu_offsetof(struct rt_vol_internal, name), BU_STRUCTPARSE_FUNC_NULL, NULL, NULL },
      {"%s", RT_VOL_NAME_LEN, "name", bu_offsetof(struct rt_vol_internal, name), BU_STRUCTPARSE_FUNC_NULL, NULL, NULL },
      {"%c", 1, "src",	VOL_O(datasrc),	BU_STRUCTPARSE_FUNC_NULL, NULL, NULL },
      {"%d", 1, "w", VOL_O(xdim), BU_STRUCTPARSE_FUNC_NULL, NULL, NULL },
      {"%d", 1, "n", VOL_O(ydim), BU_STRUCTPARSE_FUNC_NULL, NULL, NULL },
      {"%d", 1, "d", VOL_O(zdim), BU_STRUCTPARSE_FUNC_NULL, NULL, NULL },
      {"%d", 1, "lo", VOL_O(lo), BU_STRUCTPARSE_FUNC_NULL, NULL, NULL },
      {"%d", 1, "hi", VOL_O(hi), BU_STRUCTPARSE_FUNC_NULL, NULL, NULL },
      {"%f", ELEMENTS_PER_VECT, "size", bu_offsetof(struct rt_vol_internal, cellsize), BU_STRUCTPARSE_FUNC_NULL, NULL, NULL },
      {"%f", 16, "mat", bu_offsetof(struct rt_vol_internal, mat), BU_STRUCTPARSE_FUNC_NULL, NULL, NULL },
      {"", 0, (char *)0, 0, BU_STRUCTPARSE_FUNC_NULL, NULL, NULL }
  };
  
  
  extern void rt_vol_plate(point_t a, point_t b, point_t c, point_t d,
  			 mat_t mat, struct bu_list *vlfree, struct bu_list *vhead, struct rt_vol_internal *vip);
  extern int rt_retrieve_binunif(struct rt_db_internal *intern, const struct db_i *dbip, const char *name);
  extern int rt_binunif_describe(struct bu_vls  *str, const struct rt_db_internal *ip, int verbose, double mm2local);
  /*
   * Codes to represent surface normals.
   * In a bitmap, there are only 4 possible normals.
   * With this code, reverse the sign to reverse the direction.
   * As always, the normal is expected to point outwards.
   */
  #define NORM_ZPOS 3
  #define NORM_YPOS 2
  #define NORM_XPOS 1
  #define NORM_XNEG (-1)
  #define NORM_YNEG (-2)
  #define NORM_ZNEG (-3)
  
  /*
   * Regular bit addressing is used:  (0..W-1, 0..N-1),
   * but the bitmap is stored with two cells of zeros all around,
   * so permissible subscripts run (-2..W+1, -2..N+1).
   * This eliminates special-case code for the boundary conditions.
   */
  #define VOL_XWIDEN 2
  #define VOL_YWIDEN 2
  #define VOL_ZWIDEN 2
  #define VOLIDX(_xdim, _ydim, _xx, _yy, _zz) \
  	(((_zz)+VOL_ZWIDEN) \
  	 * ((_ydim) + VOL_YWIDEN*2) \
  	 + ((_yy)+VOL_YWIDEN)) \
  	* ((_xdim) + VOL_XWIDEN*2) \
  	+ (_xx)+VOL_XWIDEN
  #define VOLMAP(_map, _xdim, _ydim, _xx, _yy, _zz) (_map)[ VOLIDX(_xdim, _ydim, _xx, _yy, _zz) ]
  #define VOL(_vip, _xx, _yy, _zz) VOLMAP((_vip)->map, (_vip)->xdim, (_vip)->ydim, _xx, _yy, _zz)
  
  #define OK(_vip, _v)	((_v) >= (_vip)->lo && (_v) <= (_vip)->hi)
  
  static int rt_vol_normtab[3] = { NORM_XPOS, NORM_YPOS, NORM_ZPOS };
  
  /**
   * Transform the ray into local coordinates of the volume ("ideal space").
   * Step through the 3-D array, in local coordinates.
   * Return intersection segments.
   *
   */
  int
  rt_vol_shot(struct soltab *stp, register struct xray *rp, struct application *ap, struct seg *seghead)
  {
      register struct rt_vol_specific *volp =
  	(struct rt_vol_specific *)stp->st_specific;
      vect_t invdir;
      double t0;	/* in point of cell */
      double t1;	/* out point of cell */
      double tmax;	/* out point of entire grid */
      vect_t t;	/* next t value for XYZ cell plane intersect */
      vect_t delta;	/* spacing of XYZ cell planes along ray */
      int igrid[3];/* Grid cell coordinates of cell (integerized) */
      vect_t P;	/* hit point */
      int inside;	/* inside/outside a solid flag */
      int in_axis;
      int axis_set = 0;
      int out_axis;
      int j;
      struct xray ideal_ray;
  
      /* Transform actual ray into ideal space at origin in X-Y plane */
      MAT4X3PNT(ideal_ray.r_pt, volp->vol_mat, rp->r_pt);
      MAT4X3VEC(ideal_ray.r_dir, volp->vol_mat, rp->r_dir);
      rp = &ideal_ray;	/* XXX */
  
      /* Compute the inverse of the direction cosines */
      if (!ZERO(rp->r_dir[X])) {
  	invdir[X] = 1.0/rp->r_dir[X];
      } else {
  	invdir[X] = INFINITY;
  	rp->r_dir[X] = 0.0;
      }
      if (!ZERO(rp->r_dir[Y])) {
  	invdir[Y] = 1.0/rp->r_dir[Y];
      } else {
  	invdir[Y] = INFINITY;
  	rp->r_dir[Y] = 0.0;
      }
      if (!ZERO(rp->r_dir[Z])) {
  	invdir[Z] = 1.0/rp->r_dir[Z];
      } else {
  	invdir[Z] = INFINITY;
  	rp->r_dir[Z] = 0.0;
      }
  
      /* intersect ray with ideal grid rpp */
      VSETALL(P, 0);
      if (! rt_in_rpp(rp, invdir, P, volp->vol_large))
  	return 0;	/* MISS */
      VJOIN1(P, rp->r_pt, rp->r_min, rp->r_dir);	/* P is hit point */
      if (RT_G_DEBUG&RT_DEBUG_VOL)VPRINT("vol_large", volp->vol_large);
      if (RT_G_DEBUG&RT_DEBUG_VOL)VPRINT("vol_origin", volp->vol_origin);
      if (RT_G_DEBUG&RT_DEBUG_VOL)VPRINT("r_pt", rp->r_pt);
      if (RT_G_DEBUG&RT_DEBUG_VOL)VPRINT("P", P);
      if (RT_G_DEBUG&RT_DEBUG_VOL)VPRINT("cellsize", volp->vol_i.cellsize);
      t0 = rp->r_min;
      tmax = rp->r_max;
      if (RT_G_DEBUG&RT_DEBUG_VOL)bu_log("[shoot: r_min=%g, r_max=%g]\n", rp->r_min, rp->r_max);
  
      /* find grid cell where ray first hits ideal space bounding RPP */
      igrid[X] = (P[X] - volp->vol_origin[X]) / volp->vol_i.cellsize[X];
      igrid[Y] = (P[Y] - volp->vol_origin[Y]) / volp->vol_i.cellsize[Y];
      igrid[Z] = (P[Z] - volp->vol_origin[Z]) / volp->vol_i.cellsize[Z];
      if (igrid[X] < 0) {
  	igrid[X] = 0;
      } else if ((size_t)igrid[X] >= volp->vol_i.xdim) {
  	igrid[X] = volp->vol_i.xdim-1;
      }
      if (igrid[Y] < 0) {
  	igrid[Y] = 0;
      } else if ((size_t)igrid[Y] >= volp->vol_i.ydim) {
  	igrid[Y] = volp->vol_i.ydim-1;
      }
      if (igrid[Z] < 0) {
  	igrid[Z] = 0;
      } else if ((size_t)igrid[Z] >= volp->vol_i.zdim) {
  	igrid[Z] = volp->vol_i.zdim-1;
      }
      if (RT_G_DEBUG&RT_DEBUG_VOL)bu_log("igrid=(%d, %d, %d)\n", igrid[X], igrid[Y], igrid[Z]);
  
      /* X setup */
      if (ZERO(rp->r_dir[X])) {
  	t[X] = INFINITY;
  	delta[X] = 0;
      } else {
  	j = igrid[X];
  	if (rp->r_dir[X] < 0) j++;
  	t[X] = (volp->vol_origin[X] + j*volp->vol_i.cellsize[X] -
  		rp->r_pt[X]) * invdir[X];
  	delta[X] = volp->vol_i.cellsize[X] * fabs(invdir[X]);
  
  	/* face of entry into first cell -- max initial t value */
  	t0 = t[X];
  	in_axis = X;
  	axis_set = 1;
      }
      /* Y setup */
      if (ZERO(rp->r_dir[Y])) {
  	t[Y] = INFINITY;
  	delta[Y] = 0;
      } else {
  	j = igrid[Y];
  	if (rp->r_dir[Y] < 0) j++;
  	t[Y] = (volp->vol_origin[Y] + j*volp->vol_i.cellsize[Y] -
  		rp->r_pt[Y]) * invdir[Y];
  	delta[Y] = volp->vol_i.cellsize[Y] * fabs(invdir[Y]);
  
  	/* face of entry into first cell -- max initial t value */
  	if (axis_set) {
  	    if (t[Y] > t0) {
  		t0 = t[Y];
  		in_axis = Y;
  	    }
  	}
  	else {
  	    t0 = t[Y];
  	    in_axis = Y;
  	    axis_set = 1;
  	}
      }
      /* Z setup */
      if (ZERO(rp->r_dir[Z])) {
  	t[Z] = INFINITY;
  	delta[Z] = 0;
      } else {
  	j = igrid[Z];
  	if (rp->r_dir[Z] < 0) j++;
  	t[Z] = (volp->vol_origin[Z] + j*volp->vol_i.cellsize[Z] -
  		rp->r_pt[Z]) * invdir[Z];
  	delta[Z] = volp->vol_i.cellsize[Z] * fabs(invdir[Z]);
  
  	/* face of entry into first cell -- max initial t value */
  	if (axis_set) {
  	    if (t[Z] > t0) {
  		t0 = t[Z];
  		in_axis = Z;
  	    }
  	}
  	else {
  	    t0 = t[Z];
  	    in_axis = Z;
  	    axis_set = 1;
  	}
      }
  
      /* The delta[] elements *must* be positive, as t must increase */
      if (RT_G_DEBUG&RT_DEBUG_VOL)bu_log("t[X] = %g, delta[X] = %g\n", t[X], delta[X]);
      if (RT_G_DEBUG&RT_DEBUG_VOL)bu_log("t[Y] = %g, delta[Y] = %g\n", t[Y], delta[Y]);
      if (RT_G_DEBUG&RT_DEBUG_VOL)bu_log("t[Z] = %g, delta[Z] = %g\n", t[Z], delta[Z]);
  
      if (!axis_set) bu_log("ERROR vol: no valid entry face found\n");
  
      if (RT_G_DEBUG&RT_DEBUG_VOL)bu_log("Entry axis is %s, t0=%g\n", in_axis==X ? "X" : (in_axis==Y?"Y":"Z"), t0);
  
      /* Advance to next exits */
      t[X] += delta[X];
      t[Y] += delta[Y];
      t[Z] += delta[Z];
  
      /* Ensure that next exit is after first entrance */
      if (t[X] < t0) {
  	bu_log("*** advancing t[X]\n");
  	t[X] += delta[X];
      }
      if (t[Y] < t0) {
  	bu_log("*** advancing t[Y]\n");
  	t[Y] += delta[Y];
      }
      if (t[Z] < t0) {
  	bu_log("*** advancing t[Z]\n");
  	t[Z] += delta[Z];
      }
      if (RT_G_DEBUG&RT_DEBUG_VOL) VPRINT("Exit t[]", t);
  
      inside = 0;
  
      while (t0 < tmax) {
  	int val;
  	struct seg *segp;
  
  	/* find minimum exit t value */
  	if (t[X] < t[Y]) {
  	    if (t[Z] < t[X]) {
  		out_axis = Z;
  		t1 = t[Z];
  	    } else {
  		out_axis = X;
  		t1 = t[X];
  	    }
  	} else {
  	    if (t[Z] < t[Y]) {
  		out_axis = Z;
  		t1 = t[Z];
  	    } else {
  		out_axis = Y;
  		t1 = t[Y];
  	    }
  	}
  
  	/* Ray passes through cell igrid[XY] from t0 to t1 */
  	val = VOL(&volp->vol_i, igrid[X], igrid[Y], igrid[Z]);
  	if (RT_G_DEBUG&RT_DEBUG_VOL)bu_log("igrid [%d %d %d] from %g to %g, val=%d\n",
  					igrid[X], igrid[Y], igrid[Z],
  					t0, t1, val);
  	if (RT_G_DEBUG&RT_DEBUG_VOL)bu_log("Exit axis is %s, t[]=(%g, %g, %g)\n",
  					out_axis==X ? "X" : (out_axis==Y?"Y":"Z"),
  					t[X], t[Y], t[Z]);
  
  	if (t1 <= t0) bu_log("ERROR vol t1=%g < t0=%g\n", t1, t0);
  	if (!inside) {
  	    if (OK(&volp->vol_i, (size_t)val)) {
  		/* Handle the transition from vacuum to solid */
  		/* Start of segment (entering a full voxel) */
  		inside = 1;
  
  		RT_GET_SEG(segp, ap->a_resource);
  		segp->seg_stp = stp;
  		segp->seg_in.hit_dist = t0;
  
  		/* Compute entry normal */
  		if (rp->r_dir[in_axis] < 0) {
  		    /* Go left, entry norm goes right */
  		    segp->seg_in.hit_surfno =
  			rt_vol_normtab[in_axis];
  		} else {
  		    /* go right, entry norm goes left */
  		    segp->seg_in.hit_surfno =
  			(-rt_vol_normtab[in_axis]);
  		}
  		BU_LIST_INSERT(&(seghead->l), &(segp->l));
  		if (RT_G_DEBUG&RT_DEBUG_VOL) bu_log("START t=%g, surfno=%d\n",
  						 t0, segp->seg_in.hit_surfno);
  	    } else {
  		/* Do nothing, marching through void */
  	    }
  	} else {
  	    if (OK(&volp->vol_i, (size_t)val)) {
  		/* Do nothing, marching through solid */
  	    } else {
  		register struct seg *tail;
  		/* End of segment (now in an empty voxel) */
  		/* Handle transition from solid to vacuum */
  		inside = 0;
  
  		tail = BU_LIST_LAST(seg, &(seghead->l));
  		tail->seg_out.hit_dist = t0;
  
  		/* Compute exit normal */
  		if (rp->r_dir[in_axis] < 0) {
  		    /* Go left, exit normal goes left */
  		    tail->seg_out.hit_surfno =
  			(-rt_vol_normtab[in_axis]);
  		} else {
  		    /* go right, exit norm goes right */
  		    tail->seg_out.hit_surfno =
  			rt_vol_normtab[in_axis];
  		}
  		if (RT_G_DEBUG&RT_DEBUG_VOL) bu_log("END t=%g, surfno=%d\n",
  						 t0, tail->seg_out.hit_surfno);
  	    }
  	}
  
  	/* Take next step */
  	t0 = t1;
  	in_axis = out_axis;
  	t[out_axis] += delta[out_axis];
  	if (rp->r_dir[out_axis] > 0) {
  	    igrid[out_axis]++;
  	} else {
  	    igrid[out_axis]--;
  	}
      }
  
      if (inside) {
  	register struct seg *tail;
  
  	/* Close off the final segment */
  	tail = BU_LIST_LAST(seg, &(seghead->l));
  	tail->seg_out.hit_dist = tmax;
  
  	/* Compute exit normal.  Previous out_axis is now in_axis */
  	if (rp->r_dir[in_axis] < 0) {
  	    /* Go left, exit normal goes left */
  	    tail->seg_out.hit_surfno = (-rt_vol_normtab[in_axis]);
  	} else {
  	    /* go right, exit norm goes right */
  	    tail->seg_out.hit_surfno = rt_vol_normtab[in_axis];
  	}
  	if (RT_G_DEBUG&RT_DEBUG_VOL) bu_log("closed END t=%g, surfno=%d\n",
  					 tmax, tail->seg_out.hit_surfno);
      }
  
      if (BU_LIST_IS_EMPTY(&(seghead->l)))
  	return 0;
      return 2;
  }
  
  
  /**
   * Read in the information from the string solid record.
   * Then, as a service to the application, read in the bitmap
   * and set up some of the associated internal variables.
   */
  int
  rt_vol_import4(struct rt_db_internal *ip, const struct bu_external *ep, const fastf_t *mat, const struct db_i *dbip)
  {
      union record *rp;
      register struct rt_vol_internal *vip;
      struct bu_vls str = BU_VLS_INIT_ZERO;
      FILE *fp;
      int nbytes;
      size_t y;
      size_t z;
      mat_t tmat;
      size_t ret;
  
      if (dbip) RT_CK_DBI(dbip);
  
      BU_CK_EXTERNAL(ep);
      rp = (union record *)ep->ext_buf;
      if (rp->u_id != DBID_STRSOL) {
  	bu_log("rt_vol_import4: defective strsol record\n");
  	return -1;
      }
  
      RT_CK_DB_INTERNAL(ip);
      ip->idb_major_type = DB5_MAJORTYPE_BRLCAD;
      ip->idb_type = ID_VOL;
      ip->idb_meth = &OBJ[ID_VOL];
      BU_ALLOC(ip->idb_ptr, struct rt_vol_internal);
  
      vip = (struct rt_vol_internal *)ip->idb_ptr;
      vip->magic = RT_VOL_INTERNAL_MAGIC;
  
      /* Establish defaults */
      MAT_IDN(vip->mat);
      vip->lo = 0;
      vip->hi = 255;
  
      /* Default VOL cell size in ideal coordinates is one unit/cell */
      VSETALL(vip->cellsize, 1);
  
      /* default source is file for compatibility with old VOL */
      vip->datasrc = RT_VOL_SRC_FILE;
  
      bu_vls_strcpy(&str, rp->ss.ss_args);
      if (bu_struct_parse(&str, rt_vol_parse, (char *)vip, NULL) < 0) {
  	bu_vls_free(&str);
  	return -2;
      }
      bu_vls_free(&str);
  
      /* Check for reasonable values */
      if (vip->name[0] == '\0' || vip->xdim < 1 ||
  	vip->ydim < 1 || vip->zdim < 1 || vip->mat[15] <= 0.0 ||
  	vip->hi > 255) {
  	bu_struct_print("Unreasonable VOL parameters", rt_vol_parse,
  			(char *)vip);
  	return -1;
      }
  
      /* Apply any modeling transforms to get final matrix */
      if (mat == NULL) mat = bn_mat_identity;
      bn_mat_mul(tmat, mat, vip->mat);
      MAT_COPY(vip->mat, tmat);
  
      /* Get bit map from .bw(5) file */
      nbytes = (vip->xdim+VOL_XWIDEN*2)*
  	(vip->ydim+VOL_YWIDEN*2)*
  	(vip->zdim+VOL_ZWIDEN*2);
      vip->map = (unsigned char *)bu_calloc(1, nbytes, "vol_import4 bitmap");
  
      bu_semaphore_acquire(BU_SEM_SYSCALL);		/* lock */
      if ((fp = fopen(vip->name, "rb")) == NULL) {
  	perror(vip->name);
  	bu_semaphore_release(BU_SEM_SYSCALL);		/* unlock */
  	return -1;
      }
      bu_semaphore_release(BU_SEM_SYSCALL);		/* unlock */
  
      /* Because of in-memory padding, read each scanline separately */
      for (z = 0; z < vip->zdim; z++) {
  	for (y = 0; y < vip->ydim; y++) {
  	    void *data = &VOLMAP(vip->map, vip->xdim, vip->ydim, 0, y, z);
  	    bu_semaphore_acquire(BU_SEM_SYSCALL);		/* lock */
  	    ret = fread(data, vip->xdim, 1, fp); /* res_syscall */
  	    bu_semaphore_release(BU_SEM_SYSCALL);		/* unlock */
  	    if (ret < 1) {
  		bu_log("rt_vol_import4(%s): Unable to read whole VOL, y=%zu, z=%zu\n",
  		       vip->name, y, z);
  		bu_semaphore_acquire(BU_SEM_SYSCALL);		/* lock */
  		fclose(fp);
  		bu_semaphore_release(BU_SEM_SYSCALL);		/* unlock */
  		return -1;
  	    }
  	}
      }
      bu_semaphore_acquire(BU_SEM_SYSCALL);		/* lock */
      fclose(fp);
      bu_semaphore_release(BU_SEM_SYSCALL);		/* unlock */
      return 0;
  }
  
  
  /**
   * The name will be added by the caller.
   */
  int
  rt_vol_export4(struct bu_external *ep, const struct rt_db_internal *ip, double local2mm, const struct db_i *dbip)
  {
      struct rt_vol_internal *vip;
      struct rt_vol_internal vol;	/* scaled version */
      union record *rec;
      struct bu_vls str = BU_VLS_INIT_ZERO;
  
      if (dbip) RT_CK_DBI(dbip);
  
      RT_CK_DB_INTERNAL(ip);
      if (ip->idb_type != ID_VOL) return -1;
      vip = (struct rt_vol_internal *)ip->idb_ptr;
      RT_VOL_CK_MAGIC(vip);
      vol = *vip;			/* struct copy */
  
      /* Apply scale factor */
      vol.mat[15] /= local2mm;
  
      BU_CK_EXTERNAL(ep);
      ep->ext_nbytes = sizeof(union record)*DB_SS_NGRAN;
      ep->ext_buf = (uint8_t *)bu_calloc(1, ep->ext_nbytes, "vol external");
      rec = (union record *)ep->ext_buf;
  
      bu_vls_struct_print(&str, rt_vol_parse, (char *)&vol);
  
      rec->ss.ss_id = DBID_STRSOL;
      bu_strlcpy(rec->ss.ss_keyword, "vol", sizeof(rec->ss.ss_keyword));
      bu_strlcpy(rec->ss.ss_args, bu_vls_addr(&str), DB_SS_LEN);
      bu_vls_free(&str);
  
      return 0;
  }
  
  
  static size_t
  vol_from_file(const char *file, size_t xdim, size_t ydim, size_t zdim, unsigned char **map)
  {
      size_t y;
      size_t z;
      size_t ret = 0;
      size_t nbytes;
      FILE *fp;
  
      /* Get bit map from .bw(5) file */
      nbytes = (xdim+VOL_XWIDEN*2)*
  	(ydim+VOL_YWIDEN*2)*
  	(zdim+VOL_ZWIDEN*2);
      *map = (unsigned char *)bu_calloc(1, nbytes, "vol_import4 bitmap");
  
      bu_semaphore_acquire(BU_SEM_SYSCALL);		/* lock */
      if ((fp = fopen(file, "rb")) == NULL) {
  	perror(file);
  	bu_semaphore_release(BU_SEM_SYSCALL);		/* unlock */
  	return 0;
      }
      bu_semaphore_release(BU_SEM_SYSCALL);		/* unlock */
  
      /* Because of in-memory padding, read each scanline separately */
      for (z = 0; z < zdim; z++) {
  	for (y = 0; y < ydim; y++) {
  	    size_t fret;
  	    void *data = &VOLMAP(*map, xdim, ydim, 0, y, z);
  
  	    bu_semaphore_acquire(BU_SEM_SYSCALL);	/* lock */
  	    fret = fread(data, xdim, 1, fp);		/* res_syscall */
  	    bu_semaphore_release(BU_SEM_SYSCALL);	/* unlock */
  	    if (fret < 1) {
  		bu_log("rt_vol_import4(%s): Unable to read whole VOL, y=%zu, z=%zu\n", file, y, z);
  		bu_semaphore_acquire(BU_SEM_SYSCALL);	/* lock */
  		fclose(fp);
  		bu_semaphore_release(BU_SEM_SYSCALL);	/* unlock */
  		return 0;
  	    }
  	    ret += xdim;
  	}
      }
      bu_semaphore_acquire(BU_SEM_SYSCALL);		/* lock */
      fclose(fp);
      bu_semaphore_release(BU_SEM_SYSCALL);		/* unlock */
  
      return ret;
  }
  
  
  /**
   * Read VOL data from external file
   * Returns :
   * 0 success
   * !0 fail
   */
  static int
  vol_file_data(struct rt_vol_internal *vip)
  {
     size_t nbytes;
  
     size_t bytes = vip->xdim * vip->ydim * vip->zdim;
   	nbytes = vol_from_file(vip->name, vip->xdim, vip->ydim, vip->zdim, &vip->map);
   	if (nbytes != bytes) {
   	    bu_log("WARNING: unexpected VOL bytes (read %zu, expected %zu) in %s\n", nbytes, bytes, vip->name);
   	}
  
     return 0;
  }
  
  
  /**
   * Read VOL data from in database object
   * Returns :
   * 0 success
   * !0 fail
   */
  static int
  get_obj_data(struct rt_vol_internal *vip, const struct db_i *dbip)
  {
      struct rt_binunif_internal *bip;
      int ret;
      int nbytes;
  
      BU_ALLOC(vip->bip, struct rt_db_internal);
  
      ret = rt_retrieve_binunif(vip->bip, dbip, vip->name);
      if (ret)
  	return -1;
  
      if (RT_G_DEBUG & RT_DEBUG_HF) {
  	bu_log("db_internal magic: 0x%08x  major: %d  minor: %d\n",
  		vip->bip->idb_magic,
  		vip->bip->idb_major_type,
  		vip->bip->idb_minor_type);
      }
  
      bip = (struct rt_binunif_internal *)vip->bip->idb_ptr;
  
      if (RT_G_DEBUG & RT_DEBUG_HF)
  	bu_log("binunif magic: 0x%08x  type: %d count:%zu data[0]:%u\n",
  		bip->magic, bip->type, bip->count, bip->u.uint8[0]);
  
      if (bip->type != DB5_MINORTYPE_BINU_8BITINT_U
  	|| (size_t)bip->count != (size_t)(vip->xdim*vip->ydim*vip->zdim))
      {
  	size_t i = 0;
  	size_t size;
  	struct bu_vls binudesc = BU_VLS_INIT_ZERO;
  	rt_binunif_describe(&binudesc, vip->bip, 0, dbip->dbi_base2local);
  
  	/* skip the first title line*/
  	size = bu_vls_strlen(&binudesc);
  	while (size > 0 && i < size && bu_vls_cstr(&binudesc)[0] != '\n') {
  	    bu_vls_nibble(&binudesc, 1);
  	}
  	if (bu_vls_cstr(&binudesc)[0] == '\n')
  	    bu_vls_nibble(&binudesc, 1);
  
  	bu_log("ERROR: Binary object '%s' has invalid data (expected type %d, found %d).\n"
  	       "       Expecting %zu 8-bit unsigned char (nuc) integer data values.\n"
  	       "       Encountered %s\n",
  	       vip->name,
  	       DB5_MINORTYPE_BINU_8BITINT_U,
  	       bip->type,
  	       (size_t)(vip->xdim*vip->ydim*vip->zdim),
  	       bu_vls_cstr(&binudesc));
  	return -2;
      }
  
      nbytes = (vip->xdim+VOL_XWIDEN*2)*(vip->ydim+VOL_YWIDEN*2)*(vip->zdim+VOL_ZWIDEN*2);
  
      size_t y, z;
      if (!vip->map) {
  	unsigned char* cp;
  
    vip->map = (unsigned char *)bu_calloc(1, nbytes, "vol_import4 bitmap");
  	cp = (unsigned char *)bip->u.uint8;
  
    for (z = 0; z < vip->zdim; z++) {
  	   for (y = 0; y < vip->ydim; y++) {
         void *data = &VOLMAP(vip->map, vip->xdim, vip->ydim, 0, y, z);
  
  	     memcpy(data, cp, vip->xdim);
         cp+= vip->xdim;
  	   }
    }
      }
      return 0;
  }
  
  /**
   * Retrieve VOL data from data source
   * Returns :
   * 0 success
   * !0 fail
   */
  static int
  get_vol_data(struct rt_vol_internal *vip, const mat_t mat, const struct db_i *dbip)
  {
    mat_t tmp;
    char *p;
  
    /* Apply Modelling transform */
    bn_mat_mul(tmp, mat, vip->mat);
    MAT_COPY(vip->mat, tmp);
    p = vip->name;
  
    switch (vip->datasrc) {
  case RT_VOL_SRC_FILE:
      /* Retrieve the data from an external file */
      if (RT_G_DEBUG & RT_DEBUG_HF)
    bu_log("getting data from file \"%s\"\n", p);
  
      if(vol_file_data(vip) != 0) {
        return 1;
        p = "file";
      }
      else {
        return 0;
      }
      break;
  case RT_VOL_SRC_OBJ:
      /* Retrieve the data from an internal db object */
      if (RT_G_DEBUG & RT_DEBUG_HF)
    bu_log("getting data from object \"%s\"\n", p);
  
      if (get_obj_data(vip, dbip) != 0) {
    p = "object";
    return 1;
      } else {
    RT_CK_DB_INTERNAL(vip->bip);
    RT_CK_BINUNIF(vip->bip->idb_ptr);
    return 0;
      }
      break;
  default:
  bu_log("%s:%d Odd vol data src '%c' s/b '%c' or '%c'\n",
    __FILE__, __LINE__, vip->datasrc,
    RT_VOL_SRC_FILE, RT_VOL_SRC_OBJ);
    }
  
    bu_log("%s", dbip->dbi_filename);
    return 0; //temporary
  }
  
  
  /**
   * Read in the information from the string solid record.
   * Then, as a service to the application, read in the bitmap
   * and set up some of the associated internal variables.
   */
  int
  rt_vol_import5(struct rt_db_internal *ip, const struct bu_external *ep, const fastf_t *mat, const struct db_i *dbip)
  {
      register struct rt_vol_internal *vip;
      struct bu_vls str = BU_VLS_INIT_ZERO;
      mat_t tmat;
  
      if (dbip) RT_CK_DBI(dbip);
  
      BU_CK_EXTERNAL(ep);
  
      RT_CK_DB_INTERNAL(ip);
      ip->idb_major_type = DB5_MAJORTYPE_BRLCAD;
      ip->idb_type = ID_VOL;
      ip->idb_meth = &OBJ[ID_VOL];
      BU_ALLOC(ip->idb_ptr, struct rt_vol_internal);
  
      vip = (struct rt_vol_internal *)ip->idb_ptr;
      vip->magic = RT_VOL_INTERNAL_MAGIC;
  
      /* Establish defaults */
      MAT_IDN(vip->mat);
      vip->lo = 0;
      vip->hi = 255;
  
      /* default source is file for compatibility with old VOL */
      vip->datasrc = RT_VOL_SRC_FILE;
  
      /* Default VOL cell size in ideal coordinates is one unit/cell */
      VSETALL(vip->cellsize, 1);
  
      bu_vls_strncpy(&str, (const char *)ep->ext_buf, ep->ext_nbytes);
      if (bu_struct_parse(&str, rt_vol_parse, (char *)vip, NULL) < 0) {
  	bu_vls_free(&str);
  	return -2;
      }
      bu_vls_free(&str);
      /* Check for reasonable values */
      if (vip->name[0] == '\0'
  	|| vip->xdim < 1 || vip->ydim < 1 || vip->zdim < 1
  	|| vip->mat[15] <= 0.0 || vip->hi > 255)
      {
  	bu_struct_print("Unreasonable VOL parameters", rt_vol_parse, (char *)vip);
  	return -1;
      }
  
      if (!mat)
  	mat = bn_mat_identity;
  
      /* Apply any modeling transforms to get final matrix */
      bn_mat_mul(tmat, mat, vip->mat);
      MAT_COPY(vip->mat, tmat);
  
      if (get_vol_data(vip, mat, dbip) == 1)
    bu_log("Couldn't find the associated file/object %s",vip->name);
  
      return 0;
  }
  
  
  /**
   * The name will be added by the caller.
   */
  int
  rt_vol_export5(struct bu_external *ep, const struct rt_db_internal *ip, double local2mm, const struct db_i *dbip)
  {
      struct rt_vol_internal *vip;
      struct rt_vol_internal vol;	/* scaled version */
      struct bu_vls str = BU_VLS_INIT_ZERO;
  
      if (dbip) RT_CK_DBI(dbip);
  
      RT_CK_DB_INTERNAL(ip);
      if (ip->idb_type != ID_VOL) return -1;
      vip = (struct rt_vol_internal *)ip->idb_ptr;
      RT_VOL_CK_MAGIC(vip);
      vol = *vip;			/* struct copy */
  
      /* Apply scale factor */
      vol.mat[15] /= local2mm;
  
      BU_CK_EXTERNAL(ep);
  
      bu_vls_struct_print(&str, rt_vol_parse, (char *)&vol);
      ep->ext_nbytes = bu_vls_strlen(&str);
      ep->ext_buf = (uint8_t *)bu_calloc(1, ep->ext_nbytes, "vol external");
  
      bu_strlcpy((char *)ep->ext_buf, bu_vls_addr(&str), ep->ext_nbytes);
      bu_vls_free(&str);
  
      return 0;
  }
  
  
  /**
   * Make human-readable formatted presentation of this solid.
   * First line describes type of solid.
   * Additional lines are indented one tab, and give parameter values.
   */
  int
  rt_vol_describe(struct bu_vls *str, const struct rt_db_internal *ip, int UNUSED(verbose), double mm2local)
  {
      struct rt_vol_internal *vip = (struct rt_vol_internal *)ip->idb_ptr;
      register int i;
      struct bu_vls substr = BU_VLS_INIT_ZERO;
      vect_t local;
  
      RT_VOL_CK_MAGIC(vip);
  
      bu_vls_strcat(str, "thresholded volumetric solid (VOL)\n");
  
      /* pretty-print dimensions in local, not storage (mm) units */
      VSCALE(local, vip->cellsize, mm2local);
      if (vip->datasrc == RT_VOL_SRC_FILE) {
      bu_vls_printf(&substr, "\tfile=\"%s\"\n\tw=%u n=%u d=%u\n\tlo=%u hi=%u\n\tsize=%g,%g,%g\n",
  		  vip->name,
  		  vip->xdim, vip->ydim, vip->zdim, vip->lo, vip->hi,
  		  V3INTCLAMPARGS(local));
      } else {
        bu_vls_printf(&substr, "\tobject name=\"%s\"\n\tw=%u n=%u d=%u\n\tlo=%u hi=%u\n\tsize=%g,%g,%g\n",
    		  vip->name,
    		  vip->xdim, vip->ydim, vip->zdim, vip->lo, vip->hi,
    		  V3INTCLAMPARGS(local));
      }
  
      bu_vls_vlscat(str, &substr);
  
      bu_vls_strcat(str, "\tmat=");
      for (i = 0; i < 15; i++) {
  	bu_vls_trunc(&substr, 0);
  	bu_vls_printf(&substr, "%g, ", INTCLAMP(vip->mat[i]));
  	bu_vls_vlscat(str, &substr);
      }
      bu_vls_trunc(&substr, 0);
      bu_vls_printf(&substr, "%g\n", INTCLAMP(vip->mat[i]));
      bu_vls_vlscat(str, &substr);
  
      bu_vls_free(&substr);
  
      return 0;
  }
  
  
  /**
   * Free the storage associated with the rt_db_internal version of this solid.
   */
  void
  rt_vol_ifree(struct rt_db_internal *ip)
  {
      register struct rt_vol_internal *vip;
  
      RT_CK_DB_INTERNAL(ip);
  
      vip = (struct rt_vol_internal *)ip->idb_ptr;
      RT_VOL_CK_MAGIC(vip);
  
      /* should be stolen by vol_specific, but check just in case */
      if (vip->map)
  	bu_free((char *)vip->map, "vol bitmap");
  
      vip->magic = 0;			/* sanity */
      vip->map = (unsigned char *)0;
      bu_free((char *)vip, "vol ifree");
      ip->idb_ptr = ((void *)0);	/* sanity */
  }
  
  
  /**
   * Calculate bounding RPP for vol
   */
  int
  rt_vol_bbox(struct rt_db_internal *ip, point_t *min, point_t *max, const struct bn_tol *UNUSED(tol))
  {
      register struct rt_vol_internal *vip;
      vect_t v1, localspace;
  
      RT_CK_DB_INTERNAL(ip);
      vip = (struct rt_vol_internal *)ip->idb_ptr;
      RT_VOL_CK_MAGIC(vip);
  
      /* Find bounding RPP of rotated local RPP */
      VSETALL(v1, 0);
      VSET(localspace, vip->xdim*vip->cellsize[0], vip->ydim*vip->cellsize[1], vip->zdim*vip->cellsize[2]);/* type conversion */
      bg_rotate_bbox((*min), (*max), vip->mat, v1, localspace);
      return 0;
  }
  
  
  /**
   * Returns -
   * 0 OK
   * !0 Failure
   *
   * Implicit return -
   * A struct rt_vol_specific is created, and its address is stored
   * in stp->st_specific for use by rt_vol_shot().
   */
  int
  rt_vol_prep(struct soltab *stp, struct rt_db_internal *ip, struct rt_i *rtip)
  {
      struct rt_vol_internal *vip;
      register struct rt_vol_specific *volp;
      vect_t norm;
      vect_t radvec;
      vect_t diam;
  
      RT_CK_SOLTAB(stp);
      RT_CK_DB_INTERNAL(ip);
      if (rtip) RT_CK_RTI(rtip);
  
      vip = (struct rt_vol_internal *)ip->idb_ptr;
      RT_VOL_CK_MAGIC(vip);
  
      BU_GET(volp, struct rt_vol_specific);
      volp->vol_i = *vip;		/* struct copy */
      vip->map = (unsigned char *)0;	/* "steal" the bitmap storage */
  
      /* build Xform matrix from model(world) to ideal(local) space */
      bn_mat_inv(volp->vol_mat, vip->mat);
  
      /* Pre-compute the necessary normals.  Rotate only. */
      VSET(norm, 1, 0, 0);
      MAT3X3VEC(volp->vol_xnorm, vip->mat, norm);
      VSET(norm, 0, 1, 0);
      MAT3X3VEC(volp->vol_ynorm, vip->mat, norm);
      VSET(norm, 0, 0, 1);
      MAT3X3VEC(volp->vol_znorm, vip->mat, norm);
  
      stp->st_specific = (void *)volp;
  
      /* Find bounding RPP of rotated local RPP */
      if (rt_vol_bbox(ip, &(stp->st_min), &(stp->st_max), &rtip->rti_tol)) return 1;
      VSET(volp->vol_large,
  	 volp->vol_i.xdim*vip->cellsize[0], volp->vol_i.ydim*vip->cellsize[1], volp->vol_i.zdim*vip->cellsize[2]);/* type conversion */
  
      /* for now, VOL origin in ideal coordinates is at origin */
      VSETALL(volp->vol_origin, 0);
      VADD2(volp->vol_large, volp->vol_large, volp->vol_origin);
  
      VSUB2(diam, stp->st_max, stp->st_min);
      VADD2SCALE(stp->st_center, stp->st_min, stp->st_max, 0.5);
      VSCALE(radvec, diam, 0.5);
      stp->st_aradius = stp->st_bradius = MAGNITUDE(radvec);
  
      return 0;		/* OK */
  }
  
  
  void
  rt_vol_print(register const struct soltab *stp)
  {
      register const struct rt_vol_specific *volp =
  	(struct rt_vol_specific *)stp->st_specific;
  
      bu_log("vol file = %s\n", volp->vol_i.name);
      bu_log("dimensions = (%u, %u, %u)\n",
  	   volp->vol_i.xdim, volp->vol_i.ydim,
  	   volp->vol_i.zdim);
      VPRINT("model cellsize", volp->vol_i.cellsize);
      VPRINT("model grid origin", volp->vol_origin);
  }
  
  
  /**
   * Given one ray distance, return the normal and
   * entry/exit point.
   * This is mostly a matter of translating the stored
   * code into the proper normal.
   */
  void
  rt_vol_norm(register struct hit *hitp, struct soltab *stp, register struct xray *rp)
  {
      register struct rt_vol_specific *volp =
  	(struct rt_vol_specific *)stp->st_specific;
  
      VJOIN1(hitp->hit_point, rp->r_pt, hitp->hit_dist, rp->r_dir);
  
      switch (hitp->hit_surfno) {
  	case NORM_XPOS:
  	    VMOVE(hitp->hit_normal, volp->vol_xnorm);
  	    break;
  	case NORM_XNEG:
  	    VREVERSE(hitp->hit_normal, volp->vol_xnorm);
  	    break;
  
  	case NORM_YPOS:
  	    VMOVE(hitp->hit_normal, volp->vol_ynorm);
  	    break;
  	case NORM_YNEG:
  	    VREVERSE(hitp->hit_normal, volp->vol_ynorm);
  	    break;
  
  	case NORM_ZPOS:
  	    VMOVE(hitp->hit_normal, volp->vol_znorm);
  	    break;
  	case NORM_ZNEG:
  	    VREVERSE(hitp->hit_normal, volp->vol_znorm);
  	    break;
  
  	default:
  	    bu_log("rt_vol_norm(%s): surfno=%d bad\n",
  		   stp->st_name, hitp->hit_surfno);
  	    VSETALL(hitp->hit_normal, 0);
  	    break;
      }
  }
  
  
  /**
   * Everything has sharp edges.  This makes things easy.
   */
  void
  rt_vol_curve(register struct curvature *cvp, register struct hit *hitp, struct soltab *stp)
  {
      if (!cvp || !hitp)
  	return;
      RT_CK_HIT(hitp);
      if (stp) RT_CK_SOLTAB(stp);
  
      bn_vec_ortho(cvp->crv_pdir, hitp->hit_normal);
      cvp->crv_c1 = cvp->crv_c2 = 0;
  }
  
  
  /**
   * Map the hit point in 2-D into the range 0..1
   * untransformed X becomes U, and Y becomes V.
   */
  void
  rt_vol_uv(struct application *ap, struct soltab *stp, register struct hit *hitp, register struct uvcoord *uvp)
  {
      if (ap) RT_CK_APPLICATION(ap);
      if (stp) RT_CK_SOLTAB(stp);
      if (hitp) RT_CK_SOLTAB(stp);
      if (!uvp)
  	return;
  
      /* XXX uv should be xy in ideal space */
  }
  
  
  void
  rt_vol_free(struct soltab *stp)
  {
      register struct rt_vol_specific *volp =
  	(struct rt_vol_specific *)stp->st_specific;
  
      /* specific steals map from vip, release here */
      if (volp->vol_i.map) {
  	bu_free((char *)volp->vol_i.map, "vol_map");
  	volp->vol_i.map = NULL; /* sanity */
      }
      BU_PUT(volp, struct rt_vol_specific);
  }
  
  
  int
  rt_vol_plot(struct bu_list *vhead, struct rt_db_internal *ip, const struct bg_tess_tol *UNUSED(ttol), const struct bn_tol *UNUSED(tol), const struct bview *UNUSED(info))
  {
      register struct rt_vol_internal *vip;
      int x, y, z;
      register short v1, v2;
      point_t a, b, c, d;
  
      BU_CK_LIST_HEAD(vhead);
      RT_CK_DB_INTERNAL(ip);
      struct bu_list *vlfree = &RTG.rtg_vlfree;
      vip = (struct rt_vol_internal *)ip->idb_ptr;
      RT_VOL_CK_MAGIC(vip);
  
      if (!vip->map) {
  	return 1;
      }
  
      /*
       * Scan across in Z & X.  For each X position, scan down Y,
       * looking for the longest run of edge.
       */
      for (z=-1; z<=(int)vip->zdim; z++) {
  	for (x=-1; x<=(int)vip->xdim; x++) {
  	    for (y=-1; y<=(int)vip->ydim; y++) {
  		v1 = VOL(vip, x, y, z);
  		v2 = VOL(vip, x+1, y, z);
  		if (OK(vip, (size_t)v1) == OK(vip, (size_t)v2)) continue;
  		/* Note start point, continue scan */
  		VSET(a, x+0.5, y-0.5, z-0.5);
  		VSET(b, x+0.5, y-0.5, z+0.5);
  		for (++y; y<=(int)vip->ydim; y++) {
  		    v1 = VOL(vip, x, y, z);
  		    v2 = VOL(vip, x+1, y, z);
  		    if (OK(vip, (size_t)v1) == OK(vip, (size_t)v2))
  			break;
  		}
  		/* End of run of edge.  One cell beyond. */
  		VSET(c, x+0.5, y-0.5, z+0.5);
  		VSET(d, x+0.5, y-0.5, z-0.5);
  		rt_vol_plate(a, b, c, d, vip->mat, vlfree, vhead, vip);
  	    }
  	}
      }
  
      /*
       * Scan up in Z & Y.  For each Y position, scan across X
       */
      for (z=-1; z<=(int)vip->zdim; z++) {
  	for (y=-1; y<=(int)vip->ydim; y++) {
  	    for (x=-1; x<=(int)vip->xdim; x++) {
  		v1 = VOL(vip, x, y, z);
  		v2 = VOL(vip, x, y+1, z);
  		if (OK(vip, (size_t)v1) == OK(vip, (size_t)v2)) continue;
  		/* Note start point, continue scan */
  		VSET(a, x-0.5, y+0.5, z-0.5);
  		VSET(b, x-0.5, y+0.5, z+0.5);
  		for (++x; x<=(int)vip->xdim; x++) {
  		    v1 = VOL(vip, x, y, z);
  		    v2 = VOL(vip, x, y+1, z);
  		    if (OK(vip, (size_t)v1) == OK(vip, (size_t)v2))
  			break;
  		}
  		/* End of run of edge.  One cell beyond */
  		VSET(c, (x-0.5), (y+0.5), (z+0.5));
  		VSET(d, (x-0.5), (y+0.5), (z-0.5));
  		rt_vol_plate(a, b, c, d, vip->mat, vlfree, vhead, vip);
  	    }
  	}
      }
  
      /*
       * Scan across in Y & X.  For each X position pair edge, scan up Z.
       */
      for (x=-1; x<=(int)vip->xdim; x++) {
  	for (z=-1; z<=(int)vip->zdim; z++) {
  	    for (y=-1; y<=(int)vip->ydim; y++) {
  		v1 = VOL(vip, x, y, z);
  		v2 = VOL(vip, x, y, z+1);
  		if (OK(vip, (size_t)v1) == OK(vip, (size_t)v2)) continue;
  		/* Note start point, continue scan */
  		VSET(a, (x-0.5), (y-0.5), (z+0.5));
  		VSET(b, (x+0.5), (y-0.5), (z+0.5));
  		for (++y; y<=(int)vip->ydim; y++) {
  		    v1 = VOL(vip, x, y, z);
  		    v2 = VOL(vip, x, y, z+1);
  		    if (OK(vip, (size_t)v1) == OK(vip, (size_t)v2))
  			break;
  		}
  		/* End of run of edge.  One cell beyond */
  		VSET(c, (x+0.5), (y-0.5), (z+0.5));
  		VSET(d, (x-0.5), (y-0.5), (z+0.5));
  		rt_vol_plate(a, b, c, d, vip->mat, vlfree, vhead, vip);
  	    }
  	}
      }
      return 0;
  }
  
  
  void
  rt_vol_plate(point_t a, point_t b, point_t c, point_t d, register mat_t mat, struct bu_list *vlfree, register struct bu_list *vhead, register struct rt_vol_internal *vip)
  {
      point_t s;		/* scaled original point */
      point_t arot, prot;
  
      BU_CK_LIST_HEAD(vhead);
  
      VELMUL(s, vip->cellsize, a);
      MAT4X3PNT(arot, mat, s);
      BV_ADD_VLIST(vlfree, vhead, arot, BV_VLIST_LINE_MOVE);
  
      VELMUL(s, vip->cellsize, b);
      MAT4X3PNT(prot, mat, s);
      BV_ADD_VLIST(vlfree, vhead, prot, BV_VLIST_LINE_DRAW);
  
      VELMUL(s, vip->cellsize, c);
      MAT4X3PNT(prot, mat, s);
      BV_ADD_VLIST(vlfree, vhead, prot, BV_VLIST_LINE_DRAW);
  
      VELMUL(s, vip->cellsize, d);
      MAT4X3PNT(prot, mat, s);
      BV_ADD_VLIST(vlfree, vhead, prot, BV_VLIST_LINE_DRAW);
  
      BV_ADD_VLIST(vlfree, vhead, arot, BV_VLIST_LINE_DRAW);
  }
  
  
  int
  rt_vol_tess(struct nmgregion **r, struct model *m, struct rt_db_internal *ip, const struct bg_tess_tol *ttol, const struct bn_tol *tol)
  {
      struct rt_vol_internal *vip;
      register size_t x, y, z;
      int i;
      struct shell *s;
      struct vertex *verts[4];
      struct faceuse *fu;
      struct model *m_tmp;
      struct nmgregion *r_tmp;
  
      NMG_CK_MODEL(m);
      BN_CK_TOL(tol);
      BG_CK_TESS_TOL(ttol);
  
      RT_CK_DB_INTERNAL(ip);
      vip = (struct rt_vol_internal *)ip->idb_ptr;
      RT_VOL_CK_MAGIC(vip);
  
      /* make region, shell, vertex */
      m_tmp = nmg_mm();
      r_tmp = nmg_mrsv(m_tmp);
      s = BU_LIST_FIRST(shell, &r_tmp->s_hd);
  
      for (x = 0; x < vip->xdim; x++) {
  	for (y = 0; y < vip->ydim; y++) {
  	    for (z = 0; z < vip->zdim; z++) {
  		point_t pt, pt1;
  
  		/* skip empty cells */
  		if (!OK(vip, VOL(vip, x, y, z)))
  		    continue;
  
  		/* check neighboring cells, make a face where needed */
  
  		/* check z+1 */
  		if (!OK(vip, VOL(vip, x, y, z+1))) {
  		    for (i = 0; i < 4; i++)
  			verts[i] = (struct vertex *)NULL;
  
  		    fu = nmg_cface(s, verts, 4);
  
  		    VSET(pt, x+.5, y-.5, z+.5);
  		    VELMUL(pt1, vip->cellsize, pt);
  		    MAT4X3PNT(pt, vip->mat, pt1);
  		    nmg_vertex_gv(verts[0], pt);
  		    VSET(pt, x+.5, y+.5, z+.5);
  		    VELMUL(pt1, vip->cellsize, pt);
  		    MAT4X3PNT(pt, vip->mat, pt1);
  		    nmg_vertex_gv(verts[1], pt);
  		    VSET(pt, x-.5, y+.5, z+.5);
  		    VELMUL(pt1, vip->cellsize, pt);
  		    MAT4X3PNT(pt, vip->mat, pt1);
  		    nmg_vertex_gv(verts[2], pt);
  		    VSET(pt, x-.5, y-.5, z+.5);
  		    VELMUL(pt1, vip->cellsize, pt);
  		    MAT4X3PNT(pt, vip->mat, pt1);
  		    nmg_vertex_gv(verts[3], pt);
  
  		    if (nmg_fu_planeeqn(fu, tol))
  			goto fail;
  		}
  
  		/* check z-1 */
  		if (!OK(vip, VOL(vip, x, y, z-1))) {
  		    for (i = 0; i < 4; i++)
  			verts[i] = (struct vertex *)NULL;
  
  		    fu = nmg_cface(s, verts, 4);
  
  		    VSET(pt, x+.5, y-.5, z-.5);
  		    VELMUL(pt1, vip->cellsize, pt);
  		    MAT4X3PNT(pt, vip->mat, pt1);
  		    nmg_vertex_gv(verts[3], pt);
  		    VSET(pt, x+.5, y+.5, z-.5);
  		    VELMUL(pt1, vip->cellsize, pt);
  		    MAT4X3PNT(pt, vip->mat, pt1);
  		    nmg_vertex_gv(verts[2], pt);
  		    VSET(pt, x-.5, y+.5, z-.5);
  		    VELMUL(pt1, vip->cellsize, pt);
  		    MAT4X3PNT(pt, vip->mat, pt1);
  		    nmg_vertex_gv(verts[1], pt);
  		    VSET(pt, x-.5, y-.5, z-.5);
  		    VELMUL(pt1, vip->cellsize, pt);
  		    MAT4X3PNT(pt, vip->mat, pt1);
  		    nmg_vertex_gv(verts[0], pt);
  
  		    if (nmg_fu_planeeqn(fu, tol))
  			goto fail;
  		}
  
  		/* check y+1 */
  		if (!OK(vip, VOL(vip, x, y+1, z))) {
  		    for (i = 0; i < 4; i++)
  			verts[i] = (struct vertex *)NULL;
  
  		    fu = nmg_cface(s, verts, 4);
  
  		    VSET(pt, x+.5, y+.5, z+.5);
  		    VELMUL(pt1, vip->cellsize, pt);
  		    MAT4X3PNT(pt, vip->mat, pt1);
  		    nmg_vertex_gv(verts[0], pt);
  		    VSET(pt, x+.5, y+.5, z-.5);
  		    VELMUL(pt1, vip->cellsize, pt);
  		    MAT4X3PNT(pt, vip->mat, pt1);
  		    nmg_vertex_gv(verts[1], pt);
  		    VSET(pt, x-.5, y+.5, z-.5);
  		    VELMUL(pt1, vip->cellsize, pt);
  		    MAT4X3PNT(pt, vip->mat, pt1);
  		    nmg_vertex_gv(verts[2], pt);
  		    VSET(pt, x-.5, y+.5, z+.5);
  		    VELMUL(pt1, vip->cellsize, pt);
  		    MAT4X3PNT(pt, vip->mat, pt1);
  		    nmg_vertex_gv(verts[3], pt);
  
  		    if (nmg_fu_planeeqn(fu, tol))
  			goto fail;
  		}
  
  		/* check y-1 */
  		if (!OK(vip, VOL(vip, x, y-1, z))) {
  		    for (i = 0; i < 4; i++)
  			verts[i] = (struct vertex *)NULL;
  
  		    fu = nmg_cface(s, verts, 4);
  
  		    VSET(pt, x+.5, y-.5, z+.5);
  		    VELMUL(pt1, vip->cellsize, pt);
  		    MAT4X3PNT(pt, vip->mat, pt1);
  		    nmg_vertex_gv(verts[3], pt);
  		    VSET(pt, x+.5, y-.5, z-.5);
  		    VELMUL(pt1, vip->cellsize, pt);
  		    MAT4X3PNT(pt, vip->mat, pt1);
  		    nmg_vertex_gv(verts[2], pt);
  		    VSET(pt, x-.5, y-.5, z-.5);
  		    VELMUL(pt1, vip->cellsize, pt);
  		    MAT4X3PNT(pt, vip->mat, pt1);
  		    nmg_vertex_gv(verts[1], pt);
  		    VSET(pt, x-.5, y-.5, z+.5);
  		    VELMUL(pt1, vip->cellsize, pt);
  		    MAT4X3PNT(pt, vip->mat, pt1);
  		    nmg_vertex_gv(verts[0], pt);
  
  		    if (nmg_fu_planeeqn(fu, tol))
  			goto fail;
  		}
  
  		/* check x+1 */
  		if (!OK(vip, VOL(vip, x+1, y, z))) {
  		    for (i = 0; i < 4; i++)
  			verts[i] = (struct vertex *)NULL;
  
  		    fu = nmg_cface(s, verts, 4);
  
  		    VSET(pt, x+.5, y-.5, z-.5);
  		    VELMUL(pt1, vip->cellsize, pt);
  		    MAT4X3PNT(pt, vip->mat, pt1);
  		    nmg_vertex_gv(verts[0], pt);
  		    VSET(pt, x+.5, y+.5, z-.5);
  		    VELMUL(pt1, vip->cellsize, pt);
  		    MAT4X3PNT(pt, vip->mat, pt1);
  		    nmg_vertex_gv(verts[1], pt);
  		    VSET(pt, x+.5, y+.5, z+.5);
  		    VELMUL(pt1, vip->cellsize, pt);
  		    MAT4X3PNT(pt, vip->mat, pt1);
  		    nmg_vertex_gv(verts[2], pt);
  		    VSET(pt, x+.5, y-.5, z+.5);
  		    VELMUL(pt1, vip->cellsize, pt);
  		    MAT4X3PNT(pt, vip->mat, pt1);
  		    nmg_vertex_gv(verts[3], pt);
  
  		    if (nmg_fu_planeeqn(fu, tol))
  			goto fail;
  		}
  
  		/* check x-1 */
  		if (!OK(vip, VOL(vip, x-1, y, z))) {
  		    for (i = 0; i < 4; i++)
  			verts[i] = (struct vertex *)NULL;
  
  		    fu = nmg_cface(s, verts, 4);
  
  		    VSET(pt, x-.5, y-.5, z-.5);
  		    VELMUL(pt1, vip->cellsize, pt);
  		    MAT4X3PNT(pt, vip->mat, pt1);
  		    nmg_vertex_gv(verts[3], pt);
  		    VSET(pt, x-.5, y+.5, z-.5);
  		    VELMUL(pt1, vip->cellsize, pt);
  		    MAT4X3PNT(pt, vip->mat, pt1);
  		    nmg_vertex_gv(verts[2], pt);
  		    VSET(pt, x-.5, y+.5, z+.5);
  		    VELMUL(pt1, vip->cellsize, pt);
  		    MAT4X3PNT(pt, vip->mat, pt1);
  		    nmg_vertex_gv(verts[1], pt);
  		    VSET(pt, x-.5, y-.5, z+.5);
  		    VELMUL(pt1, vip->cellsize, pt);
  		    MAT4X3PNT(pt, vip->mat, pt1);
  		    nmg_vertex_gv(verts[0], pt);
  
  		    if (nmg_fu_planeeqn(fu, tol))
  			goto fail;
  		}
  	    }
  	}
      }
  
      nmg_region_a(r_tmp, tol);
  
      /* fuse model */
      nmg_model_fuse(m_tmp, &RTG.rtg_vlfree, tol);
  
      /* simplify shell */
      nmg_shell_coplanar_face_merge(s, tol, 1, &RTG.rtg_vlfree);
  
      /* kill snakes */
      for (BU_LIST_FOR(fu, faceuse, &s->fu_hd)) {
  	struct loopuse *lu;
  
  	NMG_CK_FACEUSE(fu);
  
  	if (fu->orientation != OT_SAME)
  	    continue;
  
  	for (BU_LIST_FOR(lu, loopuse, &fu->lu_hd))
  	    (void)nmg_kill_snakes(lu,&RTG.rtg_vlfree);
      }
  
      (void)nmg_unbreak_region_edges((uint32_t *)(&s->l),&RTG.rtg_vlfree);
  
      (void)nmg_mark_edges_real((uint32_t *)&s->l,&RTG.rtg_vlfree);
  
      nmg_merge_models(m, m_tmp);
      *r = r_tmp;
  
      return 0;
  
  fail:
      nmg_km(m_tmp);
      *r = (struct nmgregion *)NULL;
  
      return -1;
  }
  
  
  int
  rt_vol_params(struct pc_pc_set *UNUSED(ps), const struct rt_db_internal *ip)
  {
      if (ip) RT_CK_DB_INTERNAL(ip);
  
      return 0;			/* OK */
  }
  
  
  void
  rt_vol_centroid(point_t *cent, const struct rt_db_internal *ip)
  {
      register struct rt_vol_internal *vip;
      size_t x, y, z;
      size_t cnt;
      fastf_t x_tot, y_tot, z_tot;
      point_t p;
  
      RT_CK_DB_INTERNAL(ip);
      vip = (struct rt_vol_internal *)ip->idb_ptr;
      RT_VOL_CK_MAGIC(vip);
  
      cnt = 0;
      x_tot = y_tot = z_tot = 0.0;
  
      for (x = 0; x < vip->xdim; x++) {
  	for (y = 0; y < vip->ydim; y++) {
  	    for (z = 0; z < vip->zdim; z++) {
  		if (OK(vip, VOL(vip, x, y, z))) {
  		    x_tot += (x+.5) * (vip->cellsize[X]);
  		    y_tot += (y+.5) * (vip->cellsize[Y]);
  		    z_tot += (z+.5) * (vip->cellsize[Z]);
  		    cnt++;
  		}
  	    }
  	}
      }
  
      p[X]=x_tot/cnt;
      p[Y]=y_tot/cnt;
      p[Z]=z_tot/cnt;
  
      MAT4X3PNT(*cent, vip->mat, p);
  }
  
  
  /*
   * Computes the surface area of a volume by transforming each of the vertices by
   * the matrix and then summing the area of the faces of each cell necessary.
   * The vertices are numbered from left to right, front to back, bottom to top.
   */
  void
  rt_vol_surf_area(fastf_t *area, const struct rt_db_internal *ip)
  {
      struct rt_vol_internal *vip;
      unsigned int x, y, z;
      point_t x0, x1, x2, x3, x4, x5, x6, x7;
      point_t _x0, _x1, _x2, _x3, _x4, _x5, _x6, _x7;
      vect_t d3_0, d2_1, d7_4, d6_5, d6_0, d4_2, d4_1, d5_0, d5_3, d7_1, d7_2, d6_3, _cross;
      fastf_t _x, _y, _z, det, _area = 0.0;
  
      if (area == NULL || ip == NULL) return;
      RT_CK_DB_INTERNAL(ip);
      vip = (struct rt_vol_internal *)ip->idb_ptr;
      RT_VOL_CK_MAGIC(vip);
  
      det = fabs(bn_mat_determinant(vip->mat));
      if (EQUAL(det, 0.0)) {
  	*area = -1.0;
  	return;
      }
  
      for (x = 0; x < vip->xdim; x++) {
  	for (y = 0; y < vip->ydim; y++) {
  	    for (z = 0; z < vip->zdim; z++) {
  		if (OK(vip, VOL(vip, x, y, z))) {
  		    _x = (fastf_t)x;
  		    _y = (fastf_t)y;
  		    _z = (fastf_t)z;
  
  		    VSET(_x0, _x - 0.5 * vip->cellsize[X], _y - 0.5 * vip->cellsize[Y], _z - 0.5 * vip->cellsize[Z]);
  		    VSET(_x1, _x + 0.5 * vip->cellsize[X], _y - 0.5 * vip->cellsize[Y], _z - 0.5 * vip->cellsize[Z]);
  		    VSET(_x2, _x - 0.5 * vip->cellsize[X], _y + 0.5 * vip->cellsize[Y], _z - 0.5 * vip->cellsize[Z]);
  		    VSET(_x3, _x + 0.5 * vip->cellsize[X], _y + 0.5 * vip->cellsize[Y], _z - 0.5 * vip->cellsize[Z]);
  		    VSET(_x4, _x - 0.5 * vip->cellsize[X], _y - 0.5 * vip->cellsize[Y], _z + 0.5 * vip->cellsize[Z]);
  		    VSET(_x5, _x + 0.5 * vip->cellsize[X], _y - 0.5 * vip->cellsize[Y], _z + 0.5 * vip->cellsize[Z]);
  		    VSET(_x6, _x - 0.5 * vip->cellsize[X], _y + 0.5 * vip->cellsize[Y], _z + 0.5 * vip->cellsize[Z]);
  		    VSET(_x7, _x + 0.5 * vip->cellsize[X], _y + 0.5 * vip->cellsize[Y], _z + 0.5 * vip->cellsize[Z]);
  
  		    MAT4X3PNT(x0, vip->mat, _x0);
  		    MAT4X3PNT(x1, vip->mat, _x1);
  		    MAT4X3PNT(x2, vip->mat, _x2);
  		    MAT4X3PNT(x3, vip->mat, _x3);
  		    MAT4X3PNT(x4, vip->mat, _x4);
  		    MAT4X3PNT(x5, vip->mat, _x5);
  		    MAT4X3PNT(x6, vip->mat, _x6);
  		    MAT4X3PNT(x7, vip->mat, _x7);
  
  		    if (!OK(vip, VOL(vip, x + 1, y, z))) {
  			VSUB2(d5_3, x5, x3);
  			VSUB2(d7_1, x7, x1);
  			VCROSS(_cross, d5_3, d7_1);
  			_area += 0.5 * MAGNITUDE(_cross);
  		    }
  		    if (!OK(vip, VOL(vip, x - 1, y, z))) {
  			VSUB2(d6_0, x6, x0);
  			VSUB2(d4_2, x4, x2);
  			VCROSS(_cross, d6_0, d4_2);
  			_area += 0.5 * MAGNITUDE(_cross);
  		    }
  
  		    if (!OK(vip, VOL(vip, x, y + 1, z))) {
  			VSUB2(d7_2, x7, x2);
  			VSUB2(d6_3, x6, x3);
  			VCROSS(_cross, d7_2, d6_3);
  			_area += 0.5 * MAGNITUDE(_cross);
  		    }
  		    if (!OK(vip, VOL(vip, x, y - 1, z))) {
  			VSUB2(d4_1, x4, x1);
  			VSUB2(d5_0, x5, x0);
  			VCROSS(_cross, d4_1, d5_0);
  			_area += 0.5 * MAGNITUDE(_cross);
  		    }
  
  		    if (!OK(vip, VOL(vip, x, y, z + 1))) {
  			VSUB2(d3_0, x3, x0);
  			VSUB2(d2_1, x2, x1);
  			VCROSS(_cross, d3_0, d2_1);
  			_area += 0.5 * MAGNITUDE(_cross);
  		    }
  		    if (!OK(vip, VOL(vip, x, y, z - 1))) {
  			VSUB2(d7_4, x7, x4);
  			VSUB2(d6_5, x6, x5);
  			VCROSS(_cross, d7_4, d6_5);
  			_area += 0.5 * MAGNITUDE(_cross);
  		    }
  		}
  	    }
  	}
      }
      *area = _area;
  }
  
  
  /*
   * A paper at http://www.osti.gov/scitech/servlets/purl/632793 gives a method
   * for calculating the volume of hexahedrens from the eight vertices.
   *
   * The eight vertices are calculated, then transformed by the matrix and the
   * volume calculated from that.
   */
  void
  rt_vol_volume(fastf_t *volume, const struct rt_db_internal *ip)
  {
      struct rt_vol_internal *vip;
      unsigned int x, y, z;
      point_t x0, x1, x2, x3, x4, x5, x6, x7;
      point_t _x0, _x1, _x2, _x3, _x4, _x5, _x6, _x7;
      vect_t d7_0, d1_0, d3_5, d4_0, d5_6, d2_0, d6_3, cross1, cross2, cross3;
      fastf_t _x, _y, _z, det, _vol = 0.0;
  
      if (volume == NULL || ip == NULL) return;
      RT_CK_DB_INTERNAL(ip);
      vip = (struct rt_vol_internal *)ip->idb_ptr;
      RT_VOL_CK_MAGIC(vip);
  
      det = fabs(bn_mat_determinant(vip->mat));
      if (EQUAL(det, 0.0)) {
  	*volume = -1.0;
  	return;
      }
  
      for (x = 0; x < vip->xdim; x++) {
  	for (y = 0; y < vip->ydim; y++) {
  	    for (z = 0; z < vip->zdim; z++) {
  		if (OK(vip, VOL(vip, x, y, z))) {
  		    _x = (fastf_t)x;
  		    _y = (fastf_t)y;
  		    _z = (fastf_t)z;
  
  		    VSET(_x0, _x - 0.5 * vip->cellsize[X], _y - 0.5 * vip->cellsize[Y], _z - 0.5 * vip->cellsize[Z]);
  		    VSET(_x1, _x + 0.5 * vip->cellsize[X], _y - 0.5 * vip->cellsize[Y], _z - 0.5 * vip->cellsize[Z]);
  		    VSET(_x2, _x - 0.5 * vip->cellsize[X], _y + 0.5 * vip->cellsize[Y], _z - 0.5 * vip->cellsize[Z]);
  		    VSET(_x3, _x + 0.5 * vip->cellsize[X], _y + 0.5 * vip->cellsize[Y], _z - 0.5 * vip->cellsize[Z]);
  		    VSET(_x4, _x - 0.5 * vip->cellsize[X], _y - 0.5 * vip->cellsize[Y], _z + 0.5 * vip->cellsize[Z]);
  		    VSET(_x5, _x + 0.5 * vip->cellsize[X], _y - 0.5 * vip->cellsize[Y], _z + 0.5 * vip->cellsize[Z]);
  		    VSET(_x6, _x - 0.5 * vip->cellsize[X], _y + 0.5 * vip->cellsize[Y], _z + 0.5 * vip->cellsize[Z]);
  		    VSET(_x7, _x + 0.5 * vip->cellsize[X], _y + 0.5 * vip->cellsize[Y], _z + 0.5 * vip->cellsize[Z]);
  
  		    MAT4X3PNT(x0, vip->mat, _x0);
  		    MAT4X3PNT(x1, vip->mat, _x1);
  		    MAT4X3PNT(x2, vip->mat, _x2);
  		    MAT4X3PNT(x3, vip->mat, _x3);
  		    MAT4X3PNT(x4, vip->mat, _x4);
  		    MAT4X3PNT(x5, vip->mat, _x5);
  		    MAT4X3PNT(x6, vip->mat, _x6);
  		    MAT4X3PNT(x7, vip->mat, _x7);
  
  		    VSUB2(d7_0, x7, x0);
  		    VSUB2(d1_0, x1, x0);
  		    VSUB2(d3_5, x3, x5);
  		    VCROSS(cross1, d1_0, d3_5);
  
  		    VSUB2(d4_0, x4, x0);
  		    VSUB2(d5_6, x5, x6);
  		    VCROSS(cross2, d4_0, d5_6);
  
  		    VSUB2(d2_0, x2, x0);
  		    VSUB2(d6_3, x6, x3);
  		    VCROSS(cross3, d2_0, d6_3);
  
  		    _vol += (1.0/6.0) * (VDOT(d7_0, cross1) + VDOT(d7_0, cross2) + VDOT(d7_0, cross3));
  		}
  	    }
  	}
      }
      *volume = fabs(_vol);
  }
  
  
  /*
   * Local Variables:
   * mode: C
   * tab-width: 8
   * indent-tabs-mode: t
   * c-file-style: "stroustrup"
   * End:
   * ex: shiftwidth=4 tabstop=8
   */
  