/*
 * $Id: rmviso.c,v 1.7 2007/07/22 17:37:17 wes Exp $
 * Version: $Name: v180-alpha-02 $
 * $Revision: 1.7 $
 * $Log: rmviso.c,v $
 * Revision 1.7  2007/07/22 17:37:17  wes
 * Adding usrData and usrCallback parameters to all rmv routines
 *
 * Revision 1.6  2006/11/19 17:16:55  wes
 * Docs update for rmvK3MarchingCubes
 *
 * Revision 1.5  2005/09/12 04:06:02  wes
 * Minor documentation updates.
 *
 * Revision 1.4  2005/06/08 18:33:18  wes
 * Code cleanup to eliminate compiler warnings.
 *
 * Revision 1.3  2005/02/19 16:09:13  wes
 * Distro sync and consolidation.
 *
 * Revision 1.2  2005/02/12 00:34:08  wes
 * Removed archaic PROTO() macros around function declarations.
 *
 * Revision 1.1.1.1  2003/01/28 02:15:23  wes
 * Manual rebuild of rm150 repository.
 *
 * Revision 1.7  2003/01/16 22:21:20  wes
 * Updated all source files to reflect new organization of header files:
 * all header files formerly located in include/rmaux, include/rmi, include/rmv
 * are now located in include/rm.
 *
 * Revision 1.6  2003/01/09 16:40:41  wes
 * Fixed a bug in the isosurface code that showed up when doing block decomp
 * isosurfaces. The bug showed up as an incorrect z-size for the volume inside
 * the application-defined grid function callback.
 *
 * Revision 1.5  2001/03/31 17:10:08  wes
 * v1.4.0-alpha-2 checkin.
 *
 * Revision 1.4  2000/12/02 17:25:44  wes
 * Added parameters to app grid and data functions to support
 * multi-threaded use.
 *
 * Revision 1.3  2000/11/17 03:14:40  wes
 * In the pursuit of thread-safety, all static variables were removed
 * and replaced with dynamically-allocated arrays. This will allow
 * the isosurface generator code to be completely thread-safe (callable
 * from multiple simultaneous application threads).
 *
 * Revision 1.2  2000/08/28 01:34:50  wes
 * Removed unused code (old khoros stuff).
 *
 * Revision 1.1  2000/04/17 00:05:24  wes
 * Lots of documentation updates, courtesy of jdb.
 *
 * Revision 1.2  2000/02/29 23:43:59  wes
 * Compile warning cleanups.
 *
 * Revision 1.1.1.1  2000/02/28 21:29:40  wes
 * OpenRM 1.2 Checkin
 *
 * Revision 1.1.1.1  2000/02/28 17:18:48  wes
 * Initial entry - pre-RM120 release, source base for OpenRM 1.2.
 *
 */



/*
 * Copyright (c) 1996,  The Regents of the University of California.
 * All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions
 * are met:
 * 1. Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer.
 * 2. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the
 *    documentation and/or other materials provided with the distribution.
 * 3. All advertising materials mentioning features or use of this software
 *    must display the following acknowledgement:
 *      This product includes software developed by the Visualization
 *      Group at Lawrence Berkeley National Laboratory.
 * 4. Neither the name of the University of California, Berkeley nor of the 
 *    Lawrence Berkeley National Laboratory may be used to endorse or 
 *    promote products derived from this software without specific prior 
 *    written permission.
 *
 * THIS SOFTWARE IS PROVIDED BY THE CONTRIBUTORS ``AS IS'' AND
 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
 * ARE DISCLAIMED.  IN NO EVENT SHALL THE CONTRIBUTORS BE LIABLE
 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
 * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
 * SUCH DAMAGE.
 *
 *
 * This work is supported by the U. S. Department of Energy under contract 
 * number DE-AC03-76SF00098 between the U. S. Department of Energy and the 
 * University of California.
 *
 */
 

#include <rm/rm.h>
#include <rm/rmv.h>
#include "../rm/rmprivat.h"
#include "rmvprivt.h"
#include "rmvcell.h"

typedef struct slice
{
    int             usize, vsize;
    float         **xpts, **ypts, **zpts;
    double        **data;
    unsigned char **mask;
} slice;

#define BUFSIZE 65536 /* max number of triangles per flushed primitive */

/* PRIVATE declarations */

void        local_doIsosurface (RMvertex3D (*appgridfunc)(int i, int j, int k, int isize, int jsize, int ksize, float *baseX, float *baseY, float *baseZ, void *usrData), void *appGridUsrData, float (*appdatafunc)(int i, int j, int k, int isize, int jsize, int ksize, float *baseData, void *usrData), void *appDataUsrData, float (*appdata2func)(int i, int j, int k, int isize, int jsize, int ksize, float *baseData, void *usrData), void *appData2UsrData, const RMvisMap *vmap, int iusize, int ivsize, int iwsize, float isolevel, RMnode *n, RMvertex3D *bmin, RMvertex3D *bmax, RMvertex3D *vertexBuffer, RMvertex3D *normalBuffer, RMcolor4D *colorBuffer, int *totalTrianglesReturn, int *bufferedTrianglesReturn, float *baseX, float *baseY, float *baseZ, float *primaryData, float *secondaryData);

static void flush_triangles (RMvertex3D *, RMvertex3D *, RMcolor4D *, RMnode *n, int nTriangles);

void        generate_triangles (int index, int ui, int vi, int wi, slice *slice_prev, slice *slice0, slice *slice1, slice *slice_next, double level, int flip_normals, int do_colors, slice *cslice0, slice *cslice1, const RMvisMap *vmap, RMnode *n, RMvertex3D *vertexBuffer, RMvertex3D *normalBuffer, RMcolor4D *colorBuffer, int *totalTrianglesReturn, int *bufferedTrianglesReturn);

static void form_indices (int **, unsigned char **, unsigned char **, int, int);
static void malloc_slice (slice *, int, int);
static void free_slice (slice *);
static void malloc_2d_iarray (int ***, int, int);
static void free_2d_iarray (int ***);
static void load_slice (slice *slicep, int w, int usize, int vsize, int wsize, float isolevel, RMvertex3D (*appgridfunc)(int i, int j, int k, int isize, int jsize, int ksize, float *baseX, float *baseY, float *baseZ, void *usrData), void *appGridUsrData, float (*appdatafunc)(int i, int j, int k, int isize, int jsize, int ksize, float *baseData, void *usrData), void *appDataUsrData, float *baseX, float *baseY, float *baseZ, float *baseData);

/*
 * ----------------------------------------------------
 * @Name rmvK3MarchingCubes
 @pstart
 RMenum rmvK3MarchingCubes (RMvertex3D (*appgridfunc)(int i, int j, int k, int isize, int jsize, int ksize, float *baseX, float *baseY, float *baseZ, *void *appGridUsrData),
                            void *appGridUsrData,
		            float (*appdatafunc)(int i, int j, int k, int isize, int jsize, int ksize, float *baseData, void *appDataUsrData),
			    void *appDataUsrData,
			    float (*appdata2func)(int i, int j, int k, int isize, int jsize, int ksize, float *baseData2, void *appData2UsrData),
			    void *appData2UsrData,
			    const RMvisMap *vmap,
			    int isize,
			    int jsize,
			    int ksize,
			    float isolevel,
			    RMnode *n,
			    float *baseX,
			    float *baseY,
			    float *baseZ,
			    float *primaryData,
			    float *secondaryData)
 @pend

 @astart
 RMvertex3D (*appgridfunc)(int i, int j, int k, int isize, int jsize, int ksize, float *baseX, float *baseY, float *baseZ, void *appGridUsrData) - (input,
    required) -- A handle to a caller-supplied function that returns an
    RMvertex3D (x, y, z) corresponding to the grid point (i, j, k). RMV
    provides several parameters to the application callback: the size of
    the grid (isize, jsize, ksize); pointers to the base coordinate
    arrays (float *baseX, *baseY and *baseZ); the current (i,j,k) index being
    evaluated by RMV; and a void * to optional user data (usrGridUsrData).
  
 void *appGridUsrData - (input, optional) you can pass in a void * to
    user data; this void * will be provided as a parameter to the user-defined
    grid evalution callback. If you don't want to pass in a void * to any
    user data for the grid evaluation callback, pass in a NULL for this
    parameter.

 float (*appdatafunc)(int i, int j, int k, int isize, int jsize, int ksize, float *baseData, void *usrData) -- (input, required) -- A handle to a
    caller-supplied function that returns a float, which is the scalar value
    at the grid point (i, j, k). RMV provides several parameters to the
    application callback: the size of the grid (isize, jsize, ksize);
    pointers to the base data array (float *baseData); the current (i,j,k)
    index being  evaluated by RMV; and a void * to optional user data
    (usrData).

 void *appDataUsrData - (optional, input) you can pass in a void * to user
    data; this void * will be provided as a parameter to the user-defined
    data field evaluation callback. If you don't want to pass in a void *
    to any user data for the primary data field grid evaluation callback,
    pass in a NULL for this parameter.
			    
 float (*appdata2func)(int i, int j, int k, int isize, int jsize, int ksize, float *baseData2, void *appData2UsrData) - (input, optional) A handle to a
    caller-supplied function that returns a float, which is the scalar value
    at the grid point (i, j, k) in the secondary data field. RMV uses the float
    returned by appdata2func in conjunction with the RMvismap to compute
    vertex color.  RMV provides several parameters to the
    application callback: the size of the grid (isize, jsize, ksize);
    pointers to the base data array (float *baseData); the current (i,j,k)
    index being  evaluated by RMV; and a void * to optional user data
    (usrData). If you don't want to specify an app callback for a secondary
    data field, pass in a value of NULL for this parameter.

 void *appData2UsrData - (optional, input) you can pass in a void * to user
    data; this void * will be provided as a parameter to the user-defined
    secondary data field evaluation callback. If you don't want to pass in
    a void * to any user data for the secondary data field evaluation
    callback, pass in a NULL for this parameter.
			    
 const RMvisMap *vmap - a handle to an RMvisMap object (input).
			    
 int isize, jsize, ksize - int specifying the (u, v, w) dimensions
    of the data grid (input).

 float isolevel - a float specifying the isosurface scalar value
    (input).
			    
 RMnode *n - a handle to an RMnode (modified).

 float *baseX, *baseY, *baseZ, *primaryData, *secondaryData (input) - these
     are caller-supplied pointers to arrays corresponding to the x, y, and z
     coordinates, the primary scalar data field being isocontoured, and
     the optional secondary scalar data field used to colorize isosurface
     triangles. These arrays are not modified by rmv3KMarchingCubes, but
     instead are passed to the application-supplied callback routines.
     The objective of these pointers is to help alleviate the burden of
     requiring the application to maintain static pointers to data and
     coordinates.
 @aend

 @dstart

 Computes an isosurface using the Marching Cubes algorithm.  RMV invokes
 the caller-supplied functions at each (i,j,k) index of the structured
 grid to obtain vertex and data values at each grid point, the computes
 the surface intersecting each grid cell at the value isoLevel. If
 there is a secondary scalar data field input, then per-vertex colors
 for the isosurface triangles will be computed by computing an interpolated
 value for the secondary data field at the location intersected by the
 isosurface.

 This routine assumes a valid output RMnode to which the triangles
 generated by Marching Cubes are added.  Upon success, RM_CHILL is
 returned and the RMnode contains the geometry for the isosurface.
 Otherwise, RM_WHACKED is returned.
 
 @dend 
 * ----------------------------------------------------
 */
RMenum
rmvK3MarchingCubes (RMvertex3D (*appgridfunc)(int i, int j, int k, int isize, int jsize, int ksize, float *baseX, float *baseY, float *baseZ, void *usrData),
		    void *appGridUsrData,
		    float (*appdatafunc)(int i, int j, int k, int isize, int jsize, int ksize, float *baseData, void *usrData),
		    void *appDataUsrData,
		    float (*appdata2func)(int i, int j, int k, int isize, int jsize, int ksize, float *baseData, void *usrData),
		    void *appData2UsrData,
 		    const RMvisMap *vmap,
		    int iusize,
		    int ivsize,
		    int iwsize,
		    float isolevel,
		    RMnode *n,
		    float *baseX,
		    float *baseY,
		    float *baseZ,
		    float *primaryData,
		    float *secondaryData)
{
    RMvertex3D bmin, bmax;
    RMvertex3D *vertexBuffer, *normalBuffer;
    RMcolor4D *colorBuffer;
    int nTriangles=0, totalTriangles=0;

    /* assertions */
    if ((RM_ASSERT(appgridfunc, "rmvK3MarchingCubes() error: the grid callback function is NULL.") == RM_WHACKED) ||
	(RM_ASSERT(appdatafunc, "rmvK3MarchingCubes() error: the data callback function is NULL.") == RM_WHACKED) ||
	(RM_ASSERT(n, "rmvK3MarchingCubes() error: the return RMnode pointer is NULL. ") == RM_WHACKED))
	return(RM_WHACKED);

    if ((vmap == NULL && appdata2func != NULL) || (vmap != NULL && appdata2func == NULL))
    {
	rmError("rmvK3MarchingCubes() error: both the visualization colormap and the secondary data callback must both be NULL or both be defined. ");
	return(RM_WHACKED);
    }

    bmin.x = bmin.y = bmin.z = RM_MAXFLOAT;
    bmax.x = bmax.y = bmax.z = RM_MINFLOAT;

    totalTriangles = 0;
    nTriangles = 0;
    
    vertexBuffer = (RMvertex3D *)malloc(sizeof(RMvertex3D)*BUFSIZE*3);
    normalBuffer = (RMvertex3D *)malloc(sizeof(RMvertex3D)*BUFSIZE*3);
    
    if ((vmap != NULL) && (appdata2func != NULL))
	colorBuffer = rmColor4DNew(BUFSIZE * 3);
    else
	colorBuffer = NULL;

    local_doIsosurface(appgridfunc, appGridUsrData,
		       appdatafunc, appDataUsrData,
		       appdata2func, appData2UsrData,
		       vmap, iusize, ivsize, iwsize, isolevel,
		       n, &bmin, &bmax,
		       vertexBuffer, normalBuffer, colorBuffer,
		       &totalTriangles, &nTriangles, baseX, baseY, baseZ,
		       primaryData, secondaryData);


    /* do a final flush */
    flush_triangles(vertexBuffer, normalBuffer, colorBuffer, n, nTriangles);

    free((void *)vertexBuffer);
    free((void *)normalBuffer);
    
    if (colorBuffer != NULL)
	free((void *)colorBuffer);

    return(RM_CHILL);
}

/*
 * ----------------------------------------------------
 * @Name rmvK3MarchingCubesVertsNormals
 @pstart
 RMenum rmvK3MarchingCubesVertsNormals (RMvertex3D (*appgridfunc)(int i, int j, int k, int isize, int jsize, int ksize, float *baseX, float *baseY, float *baseZ, *void *appGridUsrData),
                            void *appGridUsrData,
		            float (*appdatafunc)(int i, int j, int k, int isize, int jsize, int ksize, float *baseData, void *appDataUsrData),
			    void *appDataUsrData,
			    float (*appdata2func)(int i, int j, int k, int isize, int jsize, int ksize, float *baseData2, void *appData2UsrData),
			    void *appData2UsrData,
			    const RMvisMap *vmap,
			    int isize,
			    int jsize,
			    int ksize,
			    float isolevel,
			    RMnode *n,
			    float *baseX,
			    float *baseY,
			    float *baseZ,
			    float *primaryData,
			    float *secondaryData,
			    int *returnNVerts,
			    RMvertex3D **returnVerts,
			    RMvertex3D **returnNormals)
 @pend

 @astart
 RMvertex3D (*appgridfunc)(int i, int j, int k, int isize, int jsize, int ksize, float *baseX, float *baseY, float *baseZ, void *appGridUsrData) - (input,
    required) -- A handle to a caller-supplied function that returns an
    RMvertex3D (x, y, z) corresponding to the grid point (i, j, k). RMV
    provides several parameters to the application callback: the size of
    the grid (isize, jsize, ksize); pointers to the base coordinate
    arrays (float *baseX, *baseY and *baseZ); the current (i,j,k) index being
    evaluated by RMV; and a void * to optional user data (usrGridUsrData).
  
 void *appGridUsrData - (input, optional) you can pass in a void * to
    user data; this void * will be provided as a parameter to the user-defined
    grid evalution callback. If you don't want to pass in a void * to any
    user data for the grid evaluation callback, pass in a NULL for this
    parameter.

 float (*appdatafunc)(int i, int j, int k, int isize, int jsize, int ksize, float *baseData, void *usrData) -- (input, required) -- A handle to a
    caller-supplied function that returns a float, which is the scalar value
    at the grid point (i, j, k). RMV provides several parameters to the
    application callback: the size of the grid (isize, jsize, ksize);
    pointers to the base data array (float *baseData); the current (i,j,k)
    index being  evaluated by RMV; and a void * to optional user data
    (usrData).

 void *appDataUsrData - (optional, input) you can pass in a void * to user
    data; this void * will be provided as a parameter to the user-defined
    data field evaluation callback. If you don't want to pass in a void *
    to any user data for the primary data field grid evaluation callback,
    pass in a NULL for this parameter.
			    
 float (*appdata2func)(int i, int j, int k, int isize, int jsize, int ksize, float *baseData2, void *appData2UsrData) - (input, optional) A handle to a
    caller-supplied function that returns a float, which is the scalar value
    at the grid point (i, j, k) in the secondary data field. RMV uses the float
    returned by appdata2func in conjunction with the RMvismap to compute
    vertex color.  RMV provides several parameters to the
    application callback: the size of the grid (isize, jsize, ksize);
    pointers to the base data array (float *baseData); the current (i,j,k)
    index being  evaluated by RMV; and a void * to optional user data
    (usrData). If you don't want to specify an app callback for a secondary
    data field, pass in a value of NULL for this parameter.

 void *appData2UsrData - (optional, input) you can pass in a void * to user
    data; this void * will be provided as a parameter to the user-defined
    secondary data field evaluation callback. If you don't want to pass in
    a void * to any user data for the secondary data field evaluation
    callback, pass in a NULL for this parameter.
			    
 const RMvisMap *vmap - a handle to an RMvisMap object (input).
			    
 int isize, jsize, ksize - int specifying the (u, v, w) dimensions
    of the data grid (input).

 float isolevel - a float specifying the isosurface scalar value
    (input).
			    
 RMnode *n - a handle to an RMnode (modified).

 float *baseX, *baseY, *baseZ, *primaryData, *secondaryData (input) - these
     are caller-supplied pointers to arrays corresponding to the x, y, and z
     coordinates, the primary scalar data field being isocontoured, and
     the optional secondary scalar data field used to colorize isosurface
     triangles. These arrays are not modified by rmv3KMarchingCubes, but
     instead are passed to the application-supplied callback routines.
     The objective of these pointers is to help alleviate the burden of
     requiring the application to maintain static pointers to data and
     coordinates.
     
 int *returnNVerts - an int pointer (modified). If non-null, will be set
     to the number of triangle vertices generated by the isocontouring
     operation.

 RMvertex3D **returnVerts - a pointer to an RMvertex3D pointer (modified).
     If the RMvertex3D ** is non-null, this routine will create an array
     of RMvertex3D containing the triangle vertices generated by the
     isocontouring operation and pass them back to the caller. The caller
     should rmVertex3DFree() these vertices when no longer needed. 

 RMvertex3D **returnNormals - a pointer to an RMvertex3D pointer (modified).
     If the RMvertex3D ** is non-null, this routine will create an array
     of RMvertex3D containing the per-vertex normals generated by the
     isocontouring operation and pass them back to the caller. The caller
     should rmVertex3DFree() these vertices when no longer needed. 
 @aend

 @dstart

 Computes an isosurface using the Marching Cubes algorithm.  RMV invokes
 the caller-supplied functions at each (i,j,k) index of the structured
 grid to obtain vertex and data values at each grid point, the computes
 the surface intersecting each grid cell at the value isoLevel. If
 there is a secondary scalar data field input, then per-vertex colors
 for the isosurface triangles will be computed by computing an interpolated
 value for the secondary data field at the location intersected by the
 isosurface.


 This routine is just like rmvK3MarchingCubes, but it offers the ability 
 to return a copy of the triangle vertices and normals computed by the
 isocontouring routine. The count of the number of vertices is returned
 through the parameter returnNVerts. The vertices are placed into a new
 array that is then handed back via the returnVerts parameter. Similarly,
 the normals are placed into a new array that is then handed back via the
 returnNormals parameter.
 
 This routine assumes a valid output RMnode to which the triangles
 generated by Marching Cubes are added.  Upon success, RM_CHILL is
 returned and the RMnode contains the geometry for the isosurface.
 Otherwise, RM_WHACKED is returned.

 11/21/06 Implementation Notes. The way we implement returning vertices 
 and normals is to run the isocontouring routine as usual, and then 
 harvest all the vertex/normal data from the RMprimitives at the
 node "isoNode." Since the node isoNode may contain RMprimitives prior
 to this routine being called, internally, we take note of the number
 of RMprimitives in the node before generating the isocontour, and
 harvest vertex/normal information only from the newly created
 RMprimitives.

 @dend 
 * ----------------------------------------------------
 */
RMenum
rmvK3MarchingCubesVertsNormals (RMvertex3D (*appgridfunc)(int i, int j, int k, int isize, int jsize, int ksize, float *baseX, float *baseY, float *baseZ, void *usrData),
				void *appGridUsrData,
				float (*appdatafunc)(int i, int j, int k, int isize, int jsize, int ksize, float *baseData, void *usrData),
				void *appDataUsrData,
				float (*appdata2func)(int i, int j, int k, int isize, int jsize, int ksize, float *baseData, void *usrData),
				void *appData2UsrData,
				const RMvisMap *vmap,
				int iusize,
				int ivsize,
				int iwsize,
				float isolevel,
				RMnode *isoNode,
				float *baseX,
				float *baseY,
				float *baseZ,
				float *primaryData,
				float *secondaryData,
				int *returnNVerts,
				RMvertex3D **returnVerts,
				RMvertex3D **returnNormals)
{
    RMvertex3D bmin, bmax;
    RMvertex3D *vertexBuffer, *normalBuffer;
    RMcolor4D *colorBuffer;
    int nTriangles=0, totalTriangles=0;

    /* 11/21/06 additions for harvesting verts/normals */
    int startNPrims, endNPrims;
    int nVerts=0;
    RMvertex3D *rv=NULL, *rn=NULL;
    int i;

    /* assertions */
    if ((RM_ASSERT(appgridfunc, "rmvK3MarchingCubesVertsNormals() error: the grid callback function is NULL.") == RM_WHACKED) ||
	(RM_ASSERT(appdatafunc, "rmvK3MarchingCubesVertsNormals() error: the data callback function is NULL.") == RM_WHACKED) ||
	(RM_ASSERT(isoNode, "rmvK3MarchingCubesVertsNormals() error: the return RMnode pointer is NULL. ") == RM_WHACKED))
	return(RM_WHACKED);

    if ((vmap == NULL && appdata2func != NULL) || (vmap != NULL && appdata2func == NULL))
    {
	rmError("rmvK3MarchingCubes() error: both the visualization colormap and the secondary data callback must both be NULL or both be defined. ");
	return(RM_WHACKED);
    }

    bmin.x = bmin.y = bmin.z = RM_MAXFLOAT;
    bmax.x = bmax.y = bmax.z = RM_MINFLOAT;

    totalTriangles = 0;
    nTriangles = 0;
    
    vertexBuffer = (RMvertex3D *)malloc(sizeof(RMvertex3D)*BUFSIZE*3);
    normalBuffer = (RMvertex3D *)malloc(sizeof(RMvertex3D)*BUFSIZE*3);
    
    if ((vmap != NULL) && (appdata2func != NULL))
	colorBuffer = rmColor4DNew(BUFSIZE * 3);
    else
	colorBuffer = NULL;

    startNPrims = rmNodeGetNumPrims(isoNode);

    local_doIsosurface(appgridfunc, appGridUsrData,
		       appdatafunc, appDataUsrData,
		       appdata2func, appData2UsrData,
		       vmap, iusize, ivsize, iwsize, isolevel,
		       isoNode, &bmin, &bmax,
		       vertexBuffer, normalBuffer, colorBuffer,
		       &totalTriangles, &nTriangles, baseX, baseY, baseZ,
		       primaryData, secondaryData);

    /* do a final flush */
    flush_triangles(vertexBuffer, normalBuffer, colorBuffer, isoNode, nTriangles);

    free((void *)vertexBuffer);
    free((void *)normalBuffer);
    
    if (colorBuffer != NULL)
	free((void *)colorBuffer);

    endNPrims = rmNodeGetNumPrims(isoNode);

    /* 
     *
     */
    nVerts=0;
    rv=NULL, rn=NULL;

    for (i=startNPrims; i<endNPrims; i++)
    {
	int vstride, nv, vveclen;
	int nstride, nn, nveclen;
	RMvertex3D *v=NULL, *n=NULL;

	/* there should be an API call for this */   
	RMprimitive *p = isoNode->prims[i];

	if (p == NULL)
	    continue; 		/* a degenerate case */

	/* assume we have RMvertex3D normals and verts - a reasonable 
	   assumption since we just created them. */
	private_rmGetBlobData(BLOB_VERTEX_INDEX, p, &vstride, &nv, (void **)&v, &vveclen);
	private_rmGetBlobData(BLOB_NORMAL_INDEX, p, &nstride, &nn, (void **)&n, &nveclen);

	/* could do sanity check on v and n, nverts & nnormals */

	if (returnVerts != NULL)
	{
	    rv = (RMvertex3D *)realloc(rv, (nVerts + nv)*sizeof(RMvertex3D));
	    memcpy((void *)(rv+nVerts), (void *)v, nv*sizeof(RMvertex3D));
	}

	if (returnNormals != NULL)
	{
	    rn = (RMvertex3D *)realloc(rn, (nVerts + nn)*sizeof(RMvertex3D));
	    memcpy((void *)(rn+nVerts), (void *)n, nn*sizeof(RMvertex3D));
	}
	nVerts += nv;
    }

    if (returnNVerts != NULL)
	*returnNVerts = nVerts;
    if (returnVerts != NULL)
	*returnVerts = rv;
    if (returnNormals != NULL)
	*returnNormals = rn;

    return(RM_CHILL);
}


/* PRIVATE
 * 
 * marching cubes workhorse
 */
void
local_doIsosurface (RMvertex3D (*appgridfunc)(int i, int j, int k, int isize, int jsize, int ksize, float *baseX, float *baseY, float *baseZ, void *usrData),
		    void *appGridUsrData,
		    float (*appdatafunc)(int i, int j, int k, int isize, int jsize, int ksize, float *baseData, void *usrData),
		    void *appDataUsrData,
		    float (*appdata2func)(int i, int j, int k, int isize, int jsize, int ksize, float *baseData, void *usrData),
		    void *appData2UsrData,
		    const RMvisMap *vmap,
		    int iusize,
		    int ivsize,
		    int iwsize,
		    float isolevel,
		    RMnode *n,
		    RMvertex3D *bmin,
		    RMvertex3D *bmax,
		    RMvertex3D *vbuf,
		    RMvertex3D *nbuf,
		    RMcolor4D *cbuf,
		    int *totalTrianglesReturn,
		    int *bufferedTrianglesReturn,
		    float *baseX,
		    float *baseY,
		    float *baseZ,
		    float *primaryData,
		    float *secondaryData)

{
    /*
     * allocate space for data, coordinate and normal buffers.  each
     * of these buffers must be large enough to accomodate 4 slices
     * of data, where each slice is u*v in size.
     *
     * the mask buffer is used to indicate whether or not a given
     * vertex is above (mask==1) the threshold, or below (mask==0)
     * the threshold.
     */
    int             i, j, k, do_color, flip = 1; /* flip needs to be promoted to a 1st class parm */
    int           **indices;
    slice          *cslice0 = NULL, *cslice1 = NULL;
    slice *slice_prev, *slice0, *slice1, *slice_next, *slice_temp;

    /* foil compiler warning */
    bmin = bmax = NULL;

    do_color = ((vmap != NULL) && (appdata2func != NULL));

    slice_prev = (slice *)malloc(sizeof(slice));
    slice0 = (slice *)malloc(sizeof(slice));
    slice1 = (slice *)malloc(sizeof(slice));
    slice_next = (slice *)malloc(sizeof(slice));

    if (do_color != 0)
    {
	cslice0 = (slice *)malloc(sizeof(slice));
	cslice1 = (slice *)malloc(sizeof(slice));
	malloc_slice(cslice0, iusize, ivsize);
	malloc_slice(cslice1, iusize, ivsize);
    }

    malloc_slice(slice_prev, iusize, ivsize);
    malloc_slice(slice0, iusize, ivsize);
    malloc_slice(slice1, iusize, ivsize);
    malloc_slice(slice_next, iusize, ivsize);
    malloc_2d_iarray(&indices, iusize, ivsize);

    load_slice(slice_prev, 0, iusize, ivsize, iwsize, isolevel, appgridfunc, appGridUsrData, appdatafunc, appDataUsrData, baseX, baseY, baseZ, primaryData);
    load_slice(slice0, 0, iusize, ivsize, iwsize, isolevel, appgridfunc, appGridUsrData, appdatafunc, appDataUsrData, baseX, baseY, baseZ, primaryData);
    load_slice(slice1, 1, iusize, ivsize, iwsize, isolevel, appgridfunc, appGridUsrData, appdatafunc, appDataUsrData, baseX, baseY, baseZ, primaryData);

    if (do_color != 0)
    {
	load_slice(cslice0, 0, iusize, ivsize, iwsize, isolevel, NULL, NULL, appdata2func, appData2UsrData, baseX, baseY, baseZ, secondaryData);
	load_slice(cslice1, 1, iusize, ivsize, iwsize, isolevel, NULL, NULL, appdata2func, appData2UsrData, baseX, baseY, baseZ, secondaryData);
    }

    if (iwsize == 2)
	load_slice(slice_next, 1, iusize, ivsize, iwsize, isolevel, appgridfunc, appGridUsrData, appdatafunc, appDataUsrData, baseX, baseY, baseZ, primaryData);
    else
        load_slice(slice_next, 2, iusize, ivsize, iwsize, isolevel, appgridfunc, appGridUsrData, appdatafunc, appDataUsrData, baseX, baseY, baseZ, primaryData);

    for (k = 0; k < (iwsize - 1); k++)
    {
	form_indices(indices, slice0->mask, slice1->mask, iusize, ivsize);
	
	for (j = 0; j < (ivsize - 1); j++)
	{
	    for (i = 0; i < (iusize - 1); i++)
	    {
		if (indices[j][i] != 0)
		    generate_triangles(indices[j][i], i, j, k,
				       slice_prev, slice0, slice1, slice_next,
				       isolevel, flip, do_color,
				       cslice0, cslice1, vmap, n,
				       vbuf, nbuf, cbuf,
				       totalTrianglesReturn,
				       bufferedTrianglesReturn);
	    }
	}
	slice_temp = slice_prev;
	slice_prev = slice0;
	slice0 = slice1;
	slice1 = slice_next;
	slice_next = slice_temp;

	slice_temp = cslice0;
	cslice0 = cslice1;
	cslice1 = slice_temp;

	if (k < (iwsize - 3))
	{
	    load_slice(slice_next, (k + 3), iusize, ivsize, iwsize, isolevel, appgridfunc, appGridUsrData, appdatafunc, appDataUsrData, baseX, baseY, baseZ, primaryData);
	    if (cslice1 != NULL)
		load_slice(cslice1, (k + 2), iusize, ivsize, iwsize, isolevel, NULL, NULL, appdata2func, appData2UsrData, baseX, baseY, baseZ, secondaryData);
	}
    }
    free_slice(slice_prev);
    free_slice(slice0);
    free_slice(slice1);
    free_slice(slice_next);
    free_2d_iarray(&indices);

    if (cslice0 != NULL)
	free_slice(cslice0);
    if (cslice1 != NULL)
	free_slice(cslice1);

    free(slice_prev);
    free(slice0);
    free(slice1);
    free(slice_next);

    if (cslice0 != NULL)
	free(cslice0);
	
    if (cslice1 != NULL)
	free(cslice1);
}


/* PRIVATE */
static void
form_indices (int **indices,
              unsigned char **mask0,
	      unsigned char **mask1,
	      int u,
	      int v)
{
    int i, j, r;
        
    for (j = 0; j < (v - 1); j++)
    {
	for (i = 0; i < (u - 1); i++)
	{
	    r = 0;
	    if (mask0[j][i])
		r += 1;
	    if (mask0[j][i + 1])
		r += 2;
	    if (mask0[j + 1][i + 1])
		r += 4;
	    if (mask0[j + 1][i])
		r += 8;
	    if (mask1[j][i])
		r += 16;
	    if (mask1[j][i + 1])
		r += 32;
	    if (mask1[j + 1][i + 1])
		r += 64;
	    if (mask1[j + 1][i])
		r += 128;
	    
	    indices[j][i] = r;
	}
    }
}

/* PRIVATE */
void
add_triangle (RMvertex3D **vlist,
	      RMvertex3D **nlist,
	      RMcolor4D *clist,
	      RMnode *n,
	      RMvertex3D *vbuf,
	      RMvertex3D *nbuf,
	      RMcolor4D *cbuf,
	      int *totalTrianglesReturn,
	      int *bufferedTrianglesReturn)
{
    int    i;
    double norm;
    int indx;

    indx = (*bufferedTrianglesReturn) * 3;

    for (i = 0; i < 3; i++, indx++)
    {
	vbuf[indx] = *vlist[i];
	
	norm = (nlist[i]->x * nlist[i]->x) + (nlist[i]->y * nlist[i]->y) + (nlist[i]->z * nlist[i]->z);

	if (norm != 0.0)
	{
	    norm = 1.0 / sqrt(norm);
	    nlist[i]->x *= norm;
	    nlist[i]->y *= norm;
	    nlist[i]->z *= norm;
	}
	nbuf[indx] = *nlist[i];
	
	if (clist != NULL)
	    cbuf[indx] = clist[i];
	
    }
    *bufferedTrianglesReturn += 1;
    *totalTrianglesReturn += 1;

    if (*bufferedTrianglesReturn == BUFSIZE)
    {
	/*
	 * if we've filled the buffer, dump it and reset the
	 * accumulator to zero.
	 */
	flush_triangles(vbuf, nbuf, cbuf, n,
			*bufferedTrianglesReturn);
	*bufferedTrianglesReturn = 0;
    }
}



/* PRIVATE */
static void
malloc_slice (slice *slicep, int u, int v)
{
    unsigned char *m2;
    int            i, t;
    float         *fx, *fy, *fz;
    double        *d2;

    slicep->usize = u;
    slicep->vsize = v;

    fx = (float *)malloc(sizeof(float) * u * v);
    slicep->xpts = (float **)malloc(sizeof(float *) * v);

    fy = (float *)malloc(sizeof(float) * u * v);
    slicep->ypts = (float **)malloc(sizeof(float *) * v);

    fz = (float *)malloc(sizeof(float) * u * v);
    slicep->zpts = (float **)malloc(sizeof(float *) * v);

    d2 = (double *)malloc(sizeof(double) * u * v);
    slicep->data = (double **)malloc(sizeof(double *) * v);

    m2 = (unsigned char *)malloc(sizeof(unsigned char) * u * v);
    slicep->mask = (unsigned char **)malloc(sizeof(unsigned char *) * v);
    
    for (i = 0, t = 0; i < v; i++, t += u)
    {
	slicep->xpts[i] = fx + t;
	slicep->ypts[i] = fy + t;
	slicep->zpts[i] = fz + t;
	slicep->data[i] = d2 + t;
	slicep->mask[i] = m2 + t;
    }
}


/* PRIVATE */
static void
free_slice (slice *slicep)
{
    free(slicep->xpts[0]);
    free(slicep->ypts[0]);
    free(slicep->zpts[0]);
    free(slicep->data[0]);
    free(slicep->mask[0]);

    free(slicep->xpts);
    free(slicep->ypts);
    free(slicep->zpts);
    free(slicep->data);
    free(slicep->mask);
}


/* PRIVATE */
static void
malloc_2d_iarray (int ***iarray, int u, int v) 
{
    int   i, index;
    int **t;
    int  *t2;

    t2 = (int *)malloc(sizeof(int) * u * v);
    t = (int **)malloc(sizeof(int *) * v);
    
    for (i = 0, index = 0; i < v; i++, index += u)
    {
	t[i] = t2 + index;
    }
    *iarray = t;
}


/* PRIVATE */
static void
free_2d_iarray (int ***iarray)
{
    int **t;

    t = *iarray;
    free(t[0]);
    free(t);
}


/* PRIVATE */
static void
load_slice (slice *slicep,
	    int w,
	    int usize,
	    int vsize,
	    int wsize,
	    float isolevel,
	    RMvertex3D (*appgridfunc)(int i, int j, int k, int isize, int jsize, int ksize, float *baseX, float *baseY, float *baseZ, void *usrData),
	    void *appGridUsrData,
	    float (*appdatafunc)(int i, int j, int k, int isize, int jsize, int ksize, float *baseData, void *usrData),
	    void *appDataUsrData,
	    float *baseX,
	    float *baseY,
	    float *baseZ,
	    float *baseData)
{
    unsigned char *mask_ptr;
    int            i, j;
    double         dlev;
    double        *data_ptr;
    RMvertex3D     gridpoint;

    /* compute the binary threshold mask */
    dlev = isolevel;

    for (j = 0; j < vsize; j++)
    {
	data_ptr = &(slicep->data[j][0]);
	mask_ptr = &(slicep->mask[j][0]);
	
	for (i = 0; i < usize; i++, data_ptr++, mask_ptr++)
	{
	    if (appgridfunc != NULL)
	    {
		gridpoint = (*appgridfunc)(i, j, w, usize, vsize, wsize, baseX, baseY, baseZ, appGridUsrData);
		slicep->xpts[j][i] = gridpoint.x;
		slicep->ypts[j][i] = gridpoint.y;
		slicep->zpts[j][i] = gridpoint.z;
	    }
	    *data_ptr = (*appdatafunc)(i, j, w, usize, vsize, wsize, baseData, appDataUsrData); /* float to double */
	    if (*data_ptr >= dlev)
		*mask_ptr = 1;
	    else
		*mask_ptr = 0;
	}
    }

#if 0
    /* debug: accumulate the source object's bounding box */
    for (j = 0; j < v; j++)
    {
	for (i = 0; i < u; i++)
	{
	    if (slicep->xpts[j][i] < input_loc_min.x)
		input_loc_min.x = slicep->xpts[j][i];
	    
	    if (slicep->xpts[j][i] > input_loc_max.x)
		input_loc_max.x = slicep->xpts[j][i];
	    
	    if (slicep->ypts[j][i] < input_loc_min.y)
		input_loc_min.y = slicep->ypts[j][i];
	    
	    if (slicep->ypts[j][i] > input_loc_max.y)
		input_loc_max.y = slicep->ypts[j][i];
	    
	    if (slicep->zpts[j][i] < input_loc_min.z)
		input_loc_min.z = slicep->zpts[j][i];
	    
	    if (slicep->zpts[j][i] > input_loc_max.z)
		input_loc_max.z = slicep->zpts[j][i];
	}
    }
#endif
}

#include "rmvcell.h"

static void compute_normal (RMvertex3D *, int, int, int, int, int, slice *, slice *, slice *, slice *, int);

/* PRIVATE */
double
GET_T (double d1,
       double d2,
       double t)
{
    double p;

    p = 1.0 / (d2 - d1);
    p = (t - d1) * p;
    if ((p > 1.0) || (p < 0.0))
	fprintf(stderr, "parametric value out of range.\n");
    
    return(p);
}


/* PRIVATE */
void
generate_triangles (int index,
		    int ui,
		    int vi,
		    int wi,
		    slice *slice_prev,
		    slice *slice0,
		    slice *slice1,
		    slice *slice_next,
		    double level,
		    int flip_normals,
		    int do_colors,
		    slice *cslice0,
		    slice *cslice1,
		    const RMvisMap *vmap,
		    RMnode *n,
		    RMvertex3D *vbuf,
		    RMvertex3D *nbuf,
		    RMcolor4D *cbuf,
		    int *totalTrianglesReturn,
		    int *bufferedTrianglesReturn)
{
    /*
     * this routine is rather verbose and needs to be
     * condensed. w.bethel, 4/16/99
     */
    int        i, nedges, current_edge, usize, vsize, npolys,index2;
    int        do_normals = 1;
    double     x1, x2, y1, y2, z1, z2;
    double     parametric_loc = 0.0;
    double     data1, data2, data3, data4, data5, data6, data7, data8;
    double     cdata1 = 0.0, cdata2 = 0.0, cdata3 = 0.0, cdata4 = 0.0, cdata5 = 0.0, cdata6 = 0.0, cdata7 = 0.0, cdata8 = 0.0;
    RMvertex3D norm1, norm2;

    float     crossings[13][3];
    float     grad[13][3];
    RMcolor4D tcolor[13];	/* easier to work just with 4D color */

    /* foil compiler warning */
    wi = 0;

    data1 = slice0->data[vi][ui];
    data2 = slice0->data[vi][ui+  1];
    data3 = slice0->data[vi + 1][ui + 1];
    data4 = slice0->data[vi + 1][ui];

    data5 = slice1->data[vi][ui];
    data6 = slice1->data[vi][ui + 1];
    data7 = slice1->data[vi + 1][ui + 1];
    data8 = slice1->data[vi + 1][ui];

    if (do_colors)
    {
	cdata1 = cslice0->data[vi][ui];
	cdata2 = cslice0->data[vi][ui + 1];
	cdata3 = cslice0->data[vi + 1][ui + 1];
	cdata4 = cslice0->data[vi + 1][ui];

	cdata5 = cslice1->data[vi][ui];
	cdata6 = cslice1->data[vi][ui + 1];
	cdata7 = cslice1->data[vi + 1][ui + 1];
	cdata8 = cslice1->data[vi + 1][ui];
    }
    usize = slice0->usize;
    vsize = slice0->vsize;
    
    nedges = cell_table[index].nedges;

    for (i = 0; i < nedges; i++)
    {
	current_edge = cell_table[index].edges[i];
	switch (current_edge)
	   {
	   case 1:
	      x1 = slice0->xpts[vi][ui];
	      x2 = slice0->xpts[vi][ui + 1];
	      y1 = slice0->ypts[vi][ui];
	      y2 = slice0->ypts[vi][ui + 1];
	      z1 = slice0->zpts[vi][ui];
	      z2 = slice0->zpts[vi][ui + 1];
	      
	      if (data2 == data1)
		 crossings[1][0] = x1;
	      else
		 {
		    parametric_loc = GET_T(data1, data2, level);
		    crossings[1][0] = x1 + (parametric_loc * (x2 - x1));
		 }
	      
	      crossings[1][1] = y1 + (parametric_loc * (y2 - y1)); 
	      crossings[1][2] = z1 + (parametric_loc * (z2 - z1));
	      
	      if (do_normals)
		 {
		    compute_normal(&norm1, ui, vi, 0, usize, vsize, slice_prev, slice0, slice1, slice_next, flip_normals);
		    compute_normal(&norm2, (ui + 1), vi, 0, usize, vsize, slice_prev, slice0, slice1, slice_next, flip_normals);
		    grad[1][0] = norm1.x + (parametric_loc*(norm2.x - norm1.x));
		    grad[1][1] = norm1.y + (parametric_loc*(norm2.y - norm1.y));
		    grad[1][2] = norm1.z + (parametric_loc*(norm2.z - norm1.z));
		 }
	      
	      if (do_colors)
		 {
		    int   indx;
		    float t;
		    
		    t = cdata1 + (parametric_loc * (cdata2 - cdata1));
		    indx = rmVismapIndexFromData(vmap, t);
		    rmVismapGetColor4D(vmap, indx, (tcolor + 1));
		 }
	      break;
	      
	   case 2:
	      x1 = slice0->xpts[vi][ui + 1];
	      x2 = slice0->xpts[vi + 1][ui + 1];
	      y1 = slice0->ypts[vi][ui + 1];
	      y2 = slice0->ypts[vi + 1][ui + 1];
	      z1 = slice0->zpts[vi][ui + 1];
	      z2 = slice0->zpts[vi + 1][ui + 1];
	      
	      if (data3 == data2)
		 crossings[2][1] = y1;
	      else
		 {
		    parametric_loc = GET_T(data2,data3,level);
		    crossings[2][1] = y1 + parametric_loc * (y2 - y1);
		 }
	      crossings[2][0] = x1 + (parametric_loc * (x2 - x1)); /* x2 */
	      crossings[2][2] = z1 + (parametric_loc * (z2 - z1)); 
	      
	      if (do_normals)
		 {
		    compute_normal(&norm1, (ui + 1), vi, 0, usize, vsize, slice_prev, slice0, slice1, slice_next, flip_normals);
		    compute_normal(&norm2, (ui + 1), (vi + 1), 0, usize, vsize, slice_prev, slice0, slice1, slice_next, flip_normals);
		    grad[2][0] = norm1.x + (parametric_loc * (norm2.x - norm1.x));
		    grad[2][1] = norm1.y + (parametric_loc * (norm2.y - norm1.y));
		    grad[2][2] = norm1.z + (parametric_loc * (norm2.z - norm1.z));
		 }
	      
	      if (do_colors)
		 {
		    int   indx;
		    float t;
		    
		    t = cdata2 + (parametric_loc * (cdata3 - cdata2));
		    indx = rmVismapIndexFromData(vmap, t);
		    rmVismapGetColor4D(vmap, indx, (tcolor + 2));
		 }
	      break;
	      
	   case 3:
	      x1 = slice0->xpts[vi + 1][ui + 1];
	      x2 = slice0->xpts[vi + 1][ui];
	      y1 = slice0->ypts[vi + 1][ui + 1];
	      y2 = slice0->ypts[vi + 1][ui];
	      z1 = slice0->zpts[vi + 1][ui + 1];
	      z2 = slice0->zpts[vi + 1][ui];
	      
	      if (data3 == data4)
		 crossings[3][0] = x1;
	      else
		 {
		    parametric_loc = GET_T(data3, data4, level);
		    crossings[3][0] = x1 + (parametric_loc * (x2 - x1));
		 }
	      
	      crossings[3][1] = y1 + (parametric_loc * (y2 - y1));  /* y2 in paper */
	      crossings[3][2] = z1 + (parametric_loc * (z2 - z1));
	      
	      if (do_normals)
		 {
		    compute_normal(&norm1, (ui + 1), (vi + 1), 0, usize, vsize, slice_prev, slice0, slice1, slice_next, flip_normals);
		    compute_normal(&norm2, ui, (vi + 1), 0, usize, vsize, slice_prev, slice0, slice1, slice_next, flip_normals);
		    grad[3][0] = norm1.x + (parametric_loc * (norm2.x - norm1.x));
		    grad[3][1] = norm1.y + (parametric_loc * (norm2.y - norm1.y));
		    grad[3][2] = norm1.z + (parametric_loc * (norm2.z - norm1.z));
		 }
	      if (do_colors)
		 {
		    int   indx;
		    float t;
		    
		    t = cdata3 + (parametric_loc * (cdata4 - cdata3));
		    indx = rmVismapIndexFromData(vmap, t);
		    rmVismapGetColor4D(vmap, indx, (tcolor + 3));
		 }
	      break;
	      
	   case 4:
	      x1 = slice0->xpts[vi + 1][ui];
	      x2 = slice0->xpts[vi][ui];
	      y1 = slice0->ypts[vi + 1][ui];
	      y2 = slice0->ypts[vi][ui];
	      z1 = slice0->zpts[vi + 1][ui];
	      z2 = slice0->zpts[vi][ui];
	      
	      if (data4 == data1)
		 crossings[4][1] = y1;
	      else
		 {
		    parametric_loc = GET_T(data4, data1, level);
		    crossings[4][1] = y1 + (parametric_loc * (y2 - y1));
		 }
	      
	      crossings[4][0] = x1 + (parametric_loc * (x2 - x1));
	      crossings[4][2] = z1 + (parametric_loc * (z2 - z1));
	      
	      if (do_normals)
		 {
		    compute_normal(&norm1, ui,(vi + 1), 0, usize, vsize, slice_prev, slice0, slice1, slice_next, flip_normals);
		    compute_normal(&norm2, ui, vi, 0, usize, vsize, slice_prev, slice0, slice1, slice_next, flip_normals);
		    grad[4][0] = norm1.x + (parametric_loc * (norm2.x - norm1.x));
		    grad[4][1] = norm1.y + (parametric_loc * (norm2.y - norm1.y));
		    grad[4][2] = norm1.z + (parametric_loc * (norm2.z - norm1.z));
		 }
	      
	      if (do_colors)
		 {
		    int   indx;
		    float t;
		    
		    t = cdata4 + (parametric_loc * (cdata1 - cdata4));
		    indx = rmVismapIndexFromData(vmap, t);
		    rmVismapGetColor4D(vmap, indx, (tcolor + 4));
		 }
	      break;
	      
	   case 5:
	      x1 = slice1->xpts[vi][ui];
	      x2 = slice1->xpts[vi][ui + 1];
	      y1 = slice1->ypts[vi][ui];
	      y2 = slice1->ypts[vi][ui + 1];
	      z1 = slice1->zpts[vi][ui];
	      z2 = slice1->zpts[vi][ui + 1];
	      
	      if (data6 == data5)
		 crossings[5][0] = x1;
	      else
		 {
		    parametric_loc = GET_T(data5, data6, level);
		    crossings[5][0] = x1 + (parametric_loc * (x2 - x1));
		 }
	      
	      crossings[5][1] = y1 + (parametric_loc * (y2 - y1)); 
	      crossings[5][2] = z1 + (parametric_loc * (z2 - z1)); /* z2 in paper */
	      
	      if (do_normals)
		 {
		    compute_normal(&norm1, ui, vi, 1, usize, vsize, slice_prev, slice0, slice1, slice_next, flip_normals);
		    compute_normal(&norm2, (ui + 1), vi, 1, usize, vsize, slice_prev, slice0, slice1, slice_next, flip_normals);
		    grad[5][0] = norm1.x + (parametric_loc * (norm2.x - norm1.x));
		    grad[5][1] = norm1.y + (parametric_loc * (norm2.y - norm1.y));
		    grad[5][2] = norm1.z + (parametric_loc * (norm2.z - norm1.z));
		 }
	      
	      if (do_colors)
		 {
		    int   indx;
		    float t;
		    
		    t = cdata5 + (parametric_loc * (cdata6 - cdata5));
		    indx = rmVismapIndexFromData(vmap, t);
		    rmVismapGetColor4D(vmap, indx, (tcolor + 5));
		 }
	      break;
	      
	   case 6:
	      x1 = slice1->xpts[vi][ui + 1];
	      x2 = slice1->xpts[vi + 1][ui + 1];
	      y1 = slice1->ypts[vi][ui + 1];
	      y2 = slice1->ypts[vi + 1][ui + 1];
	      z1 = slice1->zpts[vi][ui + 1];
	      z2 = slice1->zpts[vi + 1][ui + 1];
	      
	      if (data7 == data6)
		 crossings[6][1] = y1;
	      else
		 {
		    parametric_loc = GET_T(data6, data7, level);
		    crossings[6][1] = y1 + (parametric_loc * (y2 - y1));
		 }
	      crossings[6][0] = x1 + (parametric_loc * (x2 - x1)); /* x2 in paper */
	      crossings[6][2] = z1 + (parametric_loc * (z2 - z1)); /* z2 in paper */
	      
	      if (do_normals)
		 {
		    compute_normal(&norm1, (ui + 1), vi, 1, usize, vsize, slice_prev, slice0, slice1, slice_next, flip_normals);
		    compute_normal(&norm2, (ui + 1), (vi + 1), 1, usize, vsize, slice_prev, slice0, slice1, slice_next, flip_normals);
		    grad[6][0] = norm1.x + (parametric_loc * (norm2.x - norm1.x));
		    grad[6][1] = norm1.y + (parametric_loc * (norm2.y - norm1.y));
		    grad[6][2] = norm1.z + (parametric_loc * (norm2.z - norm1.z));
		 }
	      
	      if (do_colors)
		 {
		    int   indx;
		    float t;
		    
		    t = cdata6 + (parametric_loc * (cdata7 - cdata6));
		    indx = rmVismapIndexFromData(vmap, t);
		    rmVismapGetColor4D(vmap, indx, (tcolor + 6));
		 }
	      break;
	      
	   case 7:
	      x1 = slice1->xpts[vi + 1][ui + 1];
	      x2 = slice1->xpts[vi + 1][ui];
	      y1 = slice1->ypts[vi + 1][ui + 1];
	      y2 = slice1->ypts[vi + 1][ui];
	      z1 = slice1->zpts[vi + 1][ui + 1];
	      z2 = slice1->zpts[vi + 1][ui];
	      
	      if (data7 == data8)
		 crossings[7][0] = x1;
	      else
		 {
		    parametric_loc = GET_T(data7, data8, level); 
		    crossings[7][0] = x1 + (parametric_loc * (x2 - x1));
		 }
	      crossings[7][1] = y1 + (parametric_loc * (y2 - y1)); /* y2 in paper */
	      crossings[7][2] = z1 + (parametric_loc * (z2 - z1)); /* z2 in paper */
	      
	      if (do_normals)
		 {
		    compute_normal(&norm1, (ui + 1), (vi + 1), 1, usize, vsize, slice_prev, slice0, slice1, slice_next, flip_normals);
		    compute_normal(&norm2, ui, (vi + 1), 1, usize, vsize, slice_prev, slice0, slice1, slice_next, flip_normals);
		    grad[7][0] = norm1.x + (parametric_loc * (norm2.x - norm1.x));
		    grad[7][1] = norm1.y + (parametric_loc * (norm2.y - norm1.y));
		    grad[7][2] = norm1.z + (parametric_loc * (norm2.z - norm1.z));
		 }
	      
	      if (do_colors)
		 {
		    int   indx;
		    float t;
		    
		    t = cdata7 + (parametric_loc * (cdata8 - cdata7));
		    indx = rmVismapIndexFromData(vmap, t);
		    rmVismapGetColor4D(vmap, indx, (tcolor + 7));
		 }
	      break;
	      
	   case 8:
	      x1 = slice1->xpts[vi + 1][ui];
	      x2 = slice1->xpts[vi][ui];
	      y1 = slice1->ypts[vi + 1][ui];
	      y2 = slice1->ypts[vi][ui];
	      z1 = slice1->zpts[vi + 1][ui];
	      z2 = slice1->zpts[vi][ui];
	      
	      if (data8 == data5)
		 crossings[8][1] = y1;
	      else
		 {
		    parametric_loc = GET_T(data8, data5, level);
		    crossings[8][1] = y1 + (parametric_loc * (y2 - y1));
		 }
	      
	      crossings[8][0] = x1 + (parametric_loc * (x2 - x1)); 
	      crossings[8][2] = z1 + (parametric_loc * (z2 - z1)); /* z2 in paper */
	      
	      if (do_normals)
		 {
		    compute_normal(&norm1, ui, (vi + 1), 1, usize, vsize, slice_prev, slice0, slice1, slice_next, flip_normals);
		    compute_normal(&norm2, ui, vi, 1, usize, vsize, slice_prev, slice0, slice1, slice_next, flip_normals);
		    grad[8][0] = norm1.x + (parametric_loc * (norm2.x - norm1.x));
		    grad[8][1] = norm1.y + (parametric_loc * (norm2.y - norm1.y));
		    grad[8][2] = norm1.z + (parametric_loc * (norm2.z - norm1.z));
		 }
	      
	      if (do_colors)
		 {
		    int   indx;
		    float t;
		    
		    t = cdata8 + (parametric_loc * (cdata5 - cdata8));
		    indx = rmVismapIndexFromData(vmap, t);
		    rmVismapGetColor4D(vmap, indx, (tcolor + 8));
		 }
	      break;
	      
	   case 9:
	      x1 = slice0->xpts[vi][ui];
	      x2 = slice1->xpts[vi][ui];	
	      y1 = slice0->ypts[vi][ui];
	      y2 = slice1->ypts[vi][ui];
	      z1 = slice0->zpts[vi][ui];
	      z2 = slice1->zpts[vi][ui];
	      
	      if (data5 == data1)
		 crossings[9][2] = z1;
	      else
		 {
		    parametric_loc = GET_T(data1, data5, level);
		    crossings[9][2] = z1 + (parametric_loc * (z2 - z1));
		 }
	      
	      crossings[9][1] = y1 + (parametric_loc * (y2 - y1));
	      crossings[9][0] = x1 + (parametric_loc * (x2 - x1));
	      
	      if (do_normals)
		 {
		    compute_normal(&norm1, ui, vi, 0, usize, vsize, slice_prev, slice0, slice1, slice_next, flip_normals);
		    compute_normal(&norm2, ui, vi, 1, usize, vsize, slice_prev, slice0, slice1, slice_next, flip_normals);
		    grad[9][0] = norm1.x + (parametric_loc * (norm2.x - norm1.x));
		    grad[9][1] = norm1.y + (parametric_loc * (norm2.y - norm1.y));
		    grad[9][2] = norm1.z + (parametric_loc * (norm2.z - norm1.z));
		 }
	      
	      if (do_colors)
		 {
		    int   indx;
		    float t;
		    
		    t = cdata1 + (parametric_loc * (cdata5 - cdata1));
		    indx = rmVismapIndexFromData(vmap, t);
		    rmVismapGetColor4D(vmap, indx, (tcolor + 9));
		 }
	      break;
	      
	   case 10:
	      x1 = slice0->xpts[vi][ui + 1];
	      x2 = slice1->xpts[vi][ui + 1];
	      y1 = slice0->ypts[vi][ui + 1];
	      y2 = slice1->ypts[vi][ui + 1];
	      z1 = slice0->zpts[vi][ui + 1];
	      z2 = slice1->zpts[vi][ui + 1];
	      
	      if (data6 == data2)
		 crossings[10][2] = z1;
	      else
		 {
		    parametric_loc = GET_T(data2, data6, level);
		    crossings[10][2] = z1 + (parametric_loc * (z2 - z1));
		 }
	      
	      crossings[10][1] = y1 + (parametric_loc * (y2 - y1)); 
	      crossings[10][0] = x1 + (parametric_loc * (x2 - x1));  /* x2 in paper */
	      
	      if (do_normals)
		 {
		    compute_normal(&norm1, (ui + 1), vi, 0, usize, vsize, slice_prev, slice0, slice1, slice_next, flip_normals);
		    compute_normal(&norm2, (ui + 1), vi, 1, usize, vsize, slice_prev, slice0, slice1, slice_next, flip_normals);
		    grad[10][0] = norm1.x + (parametric_loc * (norm2.x - norm1.x));
		    grad[10][1] = norm1.y + (parametric_loc * (norm2.y - norm1.y));
		    grad[10][2] = norm1.z + (parametric_loc * (norm2.z - norm1.z));
		 }
	      
	      if (do_colors)
		 {
		    int   indx;
		    float t;
		    
		    t = cdata2 + (parametric_loc * (cdata6 - cdata2));
		    indx = rmVismapIndexFromData(vmap, t);
		    rmVismapGetColor4D(vmap, indx, (tcolor + 10));
		 }
	      break;
	      
	   case 11:
	      x1 = slice0->xpts[vi + 1][ui];
	      x2 = slice1->xpts[vi + 1][ui];
	      y1 = slice0->ypts[vi + 1][ui];
	      y2 = slice1->ypts[vi + 1][ui];
	      z1 = slice0->zpts[vi + 1][ui];
	      z2 = slice1->zpts[vi + 1][ui];
	      
	      if (data8 == data4)
		 crossings[11][2] = z1;
	      else
		 {
		    parametric_loc = GET_T(data4, data8, level);
		    crossings[11][2] = z1 + (parametric_loc * (z2 - z1));
		 }
	      
	      crossings[11][1] = y1 + (parametric_loc * (y2 - y1)); /* y2 in paper */
	      crossings[11][0] = x1 + (parametric_loc * (x2 - x1));
	      
	      if (do_normals)
		 {
		    compute_normal(&norm1, ui, (vi + 1), 0, usize, vsize, slice_prev, slice0, slice1, slice_next, flip_normals);
		    compute_normal(&norm2, ui, (vi + 1), 1, usize, vsize, slice_prev, slice0, slice1, slice_next, flip_normals);
		    grad[11][0] = norm1.x + (parametric_loc * (norm2.x - norm1.x));
		    grad[11][1] = norm1.y + (parametric_loc * (norm2.y - norm1.y));
		    grad[11][2] = norm1.z + (parametric_loc * (norm2.z - norm1.z));
		 }
	      if (do_colors)
		 {
		    int   indx;
		    float t;
		    t = cdata4 + (parametric_loc * (cdata8 - cdata4));
		    indx = rmVismapIndexFromData(vmap, t);
		    rmVismapGetColor4D(vmap, indx, (tcolor + 11));
		 }
	      break;
	      
	   case 12:
	      x1 = slice0->xpts[vi + 1][ui + 1];
	      x2 = slice1->xpts[vi + 1][ui + 1];
	      y1 = slice0->ypts[vi + 1][ui + 1];
	      y2 = slice1->ypts[vi + 1][ui + 1];
	      z1 = slice0->zpts[vi + 1][ui + 1];
	      z2 = slice1->zpts[vi + 1][ui + 1];
	      
	      if (data7 == data3)
		 crossings[12][2] = z1;
	      else
		 {
		    parametric_loc = GET_T(data3, data7, level);
		    crossings[12][2] = z1 + (parametric_loc * (z2 - z1));
		 }
	      crossings[12][1] = y1 + (parametric_loc * (y2 - y1));  /* y2 in paper */
	      crossings[12][0] = x1 + (parametric_loc * (x2 - x1)); /* x2 in paper */
	      
	      if (do_normals)
		 {
		    compute_normal(&norm1, (ui + 1), (vi + 1), 0, usize, vsize, slice_prev, slice0, slice1, slice_next, flip_normals);
		    compute_normal(&norm2, (ui + 1), (vi + 1), 1, usize, vsize, slice_prev, slice0, slice1, slice_next, flip_normals);
		    grad[12][0] = norm1.x + (parametric_loc * (norm2.x - norm1.x));
		    grad[12][1] = norm1.y + (parametric_loc * (norm2.y - norm1.y));
		    grad[12][2] = norm1.z + (parametric_loc * (norm2.z - norm1.z));
		 }
	      if (do_colors)
		 {
		    int   indx;
		    float t;
		    
		    t = cdata3 + (parametric_loc * (cdata7 - cdata3));
		    indx = rmVismapIndexFromData(vmap, t);
		    rmVismapGetColor4D(vmap, indx, (tcolor + 12));
		 }
	      break;
	   } /* switch for edge */
    } /* loop over number of edges (i) */

    npolys = cell_table[index].npolys;
    index2 = 0;
    {
	RMvertex3D *l_vlist[3];	 /* pointers to 3 triangle vertices */
	RMvertex3D *l_nlist[3];	 /* pointers to 3 vertex normals */
	RMcolor4D   l_colors[3];

	for (i = 0; i < npolys; i++)
	{
	    l_vlist[0] = (RMvertex3D *)&crossings[cell_table[index].polys[index2]][0];

	    if (do_normals)
		l_nlist[0] = (RMvertex3D *)&grad[cell_table[index].polys[index2]][0];
	    if (do_colors)
		l_colors[0] = tcolor[cell_table[index].polys[index2]];
	    
	    index2++;

	    l_vlist[1] = (RMvertex3D *)&crossings[cell_table[index].polys[index2]][0];
	    
	    if (do_normals)
		l_nlist[1] = (RMvertex3D *)&grad[cell_table[index].polys[index2]][0];

	    if (do_colors)
		l_colors[1] = tcolor[cell_table[index].polys[index2]];
	    
	    index2++;

	    l_vlist[2] = (RMvertex3D *)&crossings[cell_table[index].polys[index2]][0];
	    if (do_normals)
		l_nlist[2] = (RMvertex3D *)&grad[cell_table[index].polys[index2]][0];
	    if (do_colors)
		l_colors[2] = tcolor[cell_table[index].polys[index2]];
	    
	    index2++;
	    add_triangle(l_vlist, l_nlist,
			 ((do_colors != 0) ? l_colors : NULL),
			 n, vbuf, nbuf, cbuf,
			 totalTrianglesReturn, bufferedTrianglesReturn);
	}
    }
}


/* PRIVATE */
static void
compute_normal (RMvertex3D *norm,
	        int ui,
	        int vi,
	        int wi,
	        int usize,
	        int vsize,
	        slice *slice_prev,
	        slice *slice0,
	        slice *slice1,
	        slice *slice_next,
	        int flip_normals)
{
    int        i_minus, i_plus, j_minus, j_plus;
    float      jac[3][3];
    float      det;
    float      inv[3][3];
    double     mag, sign;
    slice     *sp0,* sp1, *sp2;
    RMvertex3D data;
    
    if (wi == 0)
    {
	sp0 = slice_prev;
	sp1 = slice0;
	sp2 = slice1;
    }
    else
    {
	sp0 = slice0;
	sp1 = slice1;
	sp2 = slice_next;
    }
    
    /* get neighborhood stencil */
    i_plus = (ui == (usize - 1)) ? (usize - 1) : (ui + 1);
    i_minus = (ui == 0) ? 0 : (ui - 1);
    j_plus = (vi == (vsize - 1)) ? (vsize - 1) : (vi + 1);
    j_minus = (vi == 0) ? 0 : (vi - 1);

    /* Get the derivatives of the data with respect to (u,v,w) */
    data.x = sp1->data[vi][i_plus] - sp1->data[vi][i_minus];
    data.y = sp1->data[j_plus][ui] - sp1->data[j_minus][ui];
    data.z = sp2->data[vi][ui] - sp0->data[vi][ui];

    /* Here I put the Jacobian (a 3x3 matrix) for (x,y,z) with respect to
     * (u,v,w) into a 4x4 matrix "jacobi" (why I do this is explained below)
     */
    jac[0][0] = sp1->xpts[vi][i_plus] - sp1->xpts[vi][i_minus];
    jac[0][1] = sp1->xpts[j_plus][ui] - sp1->xpts[j_minus][ui];
    jac[0][2] = sp2->xpts[vi][ui] - sp0->xpts[vi][ui];

    jac[1][0] = sp1->ypts[vi][i_plus] - sp1->ypts[vi][i_minus];
    jac[1][1] = sp1->ypts[j_plus][ui] - sp1->ypts[j_minus][ui];
    jac[1][2] = sp2->ypts[vi][ui] - sp0->ypts[vi][ui];

    jac[2][0] = sp1->zpts[vi][i_plus] - sp1->zpts[vi][i_minus];
    jac[2][1] = sp1->zpts[j_plus][ui] - sp1->zpts[j_minus][ui];
    jac[2][2] = sp2->zpts[vi][ui] - sp0->zpts[vi][ui];

    /* By taking the inverse of this Jacobian, I get the transpose of the
     * Jacobian of (u,v,w) with respect to (x,y,z).  This is what I need to
     * convert the derivatives for "data" from (u,v,w) coordinates to (x,y,z)
     * coordinates (i.e. the gradient which can be used for the normal).
     */
    if ((jac[0][1] == 0) && (jac[0][2] == 0) && (jac[1][0] == 0) && (jac[1][2] == 0) && (jac[2][0] == 0) && (jac[2][1] == 0))
     {
	if (jac[0][0] != 0)
	{
	    norm->x = data.x / jac[0][0];
	}
	else
	{
	    fprintf(stderr, "rmv: Coordinate system degenerate in x at index (%d, %d, %d)\n", ui, vi, wi);
	    norm->x = 1.0;
	}

	if (jac[1][1] != 0)
	{
	    norm->y = data.y / jac[1][1];
	}
	else
	{
	    fprintf(stderr, "rmv: Coordinate system degenerate in y at index (%d, %d, %d)\n", ui, vi, wi);
	    norm->y = 1.0;
	}

	if (jac[2][2] != 0)
	{
	    norm->z = data.z / jac[2][2];
	}
	else
	{
	    fprintf(stderr, "rmv: Coordinate system degenerate in z at index (%d, %d, %d)\n", ui, vi, wi);
	    norm->z = 1.0;
	}
    }
    else
    {
        det = (jac[0][0] * ((jac[1][1] * jac[2][2]) - (jac[2][1] * jac[1][2]))) 
	    - (jac[0][1] * ((jac[1][0] * jac[2][2]) - (jac[2][0] * jac[1][2]))) 
	    + (jac[0][2] * ((jac[1][0] * jac[2][1]) - (jac[2][0] * jac[1][1])));

	if (det != 0.0)
	{
	    det = 1.0 / det;

	    inv[0][0] =  ((jac[1][1] * jac[2][2]) - (jac[2][1] * jac[1][2])) * det;
	    inv[0][1] = -((jac[0][1] * jac[2][2]) - (jac[2][1] * jac[0][2])) * det;
	    inv[0][2] =  ((jac[0][1] * jac[1][2]) - (jac[1][1] * jac[0][2])) * det;

	    inv[1][0] = -((jac[1][0] * jac[2][2]) - (jac[2][0] * jac[1][2])) * det;
	    inv[1][1] =  ((jac[0][0] * jac[2][2]) - (jac[2][0] * jac[0][2])) * det;
	    inv[1][2] = -((jac[0][0] * jac[1][2]) - (jac[1][0] * jac[0][2])) * det;

	    inv[2][0] =  ((jac[1][0] * jac[2][1]) - (jac[2][0] * jac[1][1])) * det;
	    inv[2][1] = -((jac[0][0] * jac[2][1]) - (jac[2][0] * jac[0][1])) * det;
	    inv[2][2] =  ((jac[0][0] * jac[1][1]) - (jac[1][0] * jac[0][1])) * det;

	    /* Use the chain rule to get the "normal":
	     *
	     *    f  = f  * u  + f  * v  + f  * w
	     *     x    u    x    v    x    w    x
	     *
	     *    f  = f  * u  + f  * v  + f  * w
	     *     y    u    y    v    y    w    y
	     *
	     *    f  = f  * u  + f  * v  + f  * w
	     *     z    u    z    v    z    w    z
	     */

	    norm->x = (data.x * inv[0][0]) + (data.y * inv[1][0]) + (data.z * inv[2][0]);
	    norm->y = (data.x * inv[0][1]) + (data.y * inv[1][1]) + (data.z * inv[2][1]);
	    norm->z = (data.x * inv[0][2]) + (data.y * inv[1][2]) + (data.z * inv[2][2]);
	}
	else
	{
	    fprintf(stderr, "Coordinate system degenerate at index (%d,%d,%d)\n", ui, vi, wi);
	    norm->x = 1.0;
	    norm->y = 0.0;
	    norm->z = 0.0;
	}
    }

    /* adjust sign and normalize */
    if (flip_normals)
	sign = -1.0;
    else
	sign = 1.0;

    mag = (norm->x * norm->x) + (norm->y * norm->y) + (norm->z * norm->z);
    if (mag != 0.0)
    {
	mag = sqrt(mag);
	mag = 1.0 / mag;
	norm->x *= (mag * sign);
	norm->y *= mag * sign;
	norm->z *= mag * sign;
    }
}

/* PRIVATE */
static void
flush_triangles (RMvertex3D *vlist,
		 RMvertex3D *nlist,
		 RMcolor4D *clist,
		 RMnode *n,
		 int nTriangles)
{
    RMprimitive *p = rmPrimitiveNew(RM_TRIANGLES);

    fprintf(stderr, " writing a buffer of %d triangles. \n", nTriangles);
  
    if (nTriangles == 0)
	return;

    rmPrimitiveSetVertex3D(p, (nTriangles * 3), vlist, RM_COPY_DATA, NULL);
    rmPrimitiveSetNormal3D(p, (nTriangles * 3), nlist, RM_COPY_DATA, NULL);

    if (clist != NULL)
	rmPrimitiveSetColor4D(p, (nTriangles * 3), clist, RM_COPY_DATA, NULL);
    rmNodeAddPrimitive(n, p);
    
#if 0
    /* debug */
    {
       int i;

       for (i = 0; i < (ntriangles * 3); i++)
	  printf("\t%d\t%g,%g,%g\t%g,%g,%g\n", i, vlist[i].x, vlist[i].y, vlist[i].z, nlist[i].x, nlist[i].y, nlist[i].z);
    }
#endif
}

/* EOF */
