/*
 * Copyright (C) 1997-2009, R3vis Corporation.
 *
 * This library is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Library General Public
 * License as published by the Free Software Foundation; either
 * version 2 of the License, or (at your option) any later version.
 *
 * 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
 * Library General Public License for more details.
 *
 * You should have received a copy of the GNU Library General Public
 * License along with this library; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
 * USA, or visit http://www.gnu.org/copyleft/lgpl.html.
 *
 * Original Contributor:
 *   Wes Bethel, R3vis Corporation, Marin County, California
 *   http://www.r3vis.com/
 * Additional Contributor(s):
 *
 * The OpenRM project is located at http://openrm.sourceforge.net/.
 */
/*
 * $Id: rmvcontour.c,v 1.8 2007/07/22 17:37:17 wes Exp $
 * Version: $Name: v180-alpha-02 $
 * $Revision: 1.8 $
 * $Log: rmvcontour.c,v $
 * Revision 1.8  2007/07/22 17:37:17  wes
 * Adding usrData and usrCallback parameters to all rmv routines
 *
 * Revision 1.7  2005/06/08 18:33:18  wes
 * Code cleanup to eliminate compiler warnings.
 *
 * Revision 1.6  2005/02/19 17:56:32  wes
 * Removed unusued variables, beautification.
 *
 * Revision 1.5  2005/02/19 16:09:13  wes
 * Distro sync and consolidation.
 *
 * Revision 1.4  2005/01/23 17:11:02  wes
 * Copyright updated to 2005.
 *
 * Revision 1.3  2004/01/17 04:09:26  wes
 * Updated copyright line for 2004.
 *
 * Revision 1.2  2003/02/02 02:07:23  wes
 * Updated copyright to 2003.
 *
 * Revision 1.1.1.1  2003/01/28 02:15:23  wes
 * Manual rebuild of rm150 repository.
 *
 * Revision 1.5  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.4  2002/04/30 19:40:00  wes
 * Updated copyright dates.
 *
 * Revision 1.3  2001/03/31 17:10:08  wes
 * v1.4.0-alpha-2 checkin.
 *
 * Revision 1.2  2000/04/20 16:17:45  wes
 * JDB modifications: code rearrangement, additional docs.
 *
 * Revision 1.1  2000/04/17 00:05:23  wes
 * Lots of documentation updates, courtesy of jdb.
 *
 * 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.
 *
 */


#include <rm/rm.h>
#include <rm/rmv.h>

/* PRIVATE declarations */
int threshcross (float d0, float *lev);
void do_contour (int i1, int j1, int i2, int j2, int i3, int j3, int i4, int j4, RMvertex3D (*appgridfunc)(int i, int j, void *usrData), void *appGridUsrData, float (*appdatafunc)(int i, int j, void *usrData), void *appDataUsrData,int usize, int vsize, int offset_mask, float *lev, RMvertex3D *varray, int *nverts, RMcolor4D *this_color, RMcolor4D *carray);

/*
 * ----------------------------------------------------
 * @Name rmvJ3SliceContour
 @pstart
 void rmvJ3SliceContour (RMvertex3D (*appgridfunc)(int i, int j, void *usrData),
		         float (*appdatafunc)(int i, int j, void *usrData),
			 int axis_offset_enum,
			 int iusize,
			 int ivsize,
			 int nlevels,
			 float *levels,
			 RMcolor4D *color_list,
			 RMenum linewidth_enum,
			 RMenum linestyle_enum,
			 RMnode *node)
 @pend

 @astart 
 RMvertex3D (*appgridfunc)(int i, int j, void *usrData) - a handle to a
    caller-supplied function that returns an RMvertex3D (x, y, z)
    corresponding to the grid point (i, j) (input).
		         
 float (*appdatafunc)(int i, int j, void *usrData) - a handle to a caller-supplied
    function that returns a float which is the scalar value at the
    grid point (i, j) (input).
			 
 int axis_offset_enum - an RMenum specifying which axis to offset
    from in terms of the 2D grid.  Must be one of RMV_XAXIS_OFFSET,
    RMV_YAXIS_OFFSET, or RMV_ZAXIS_OFFSET (input).
			 
 int iusize, ivsize - int specifying the dimensions of the 2D grid of
    data (input).

 int nlevels - an int specifying the number of contour levels to
    generate (input).
			 
 float *levels - a handle to an array of scalar values to generate
    contours for (input).
	       	 
 RMcolor4D *color_list - a handle to an array of RMcolor4D colors for
    the respective contours levels (input).
			 
 RMenum linewidth_enum - an RMenum specifying the line width.  Must be
    one of RM_LINEWIDTH_NARROW, RM_LINEWIDTH_MEDIUM,
    RM_LINEWIDTH_HEAVY, or RM_LINEWIDTH_[1..8] (input).

 RMenum linestyle_enum - an RMenum specifying the lline style.  Must
    be one of RM_LINES_SOLID, RM_LINES_DASHED, RM_LINES_DOTTED,
    RM_LINES_DOT_DASH, or RM_LINES_DASH_DASH_DOT (input).
			 
 RMnode *n - a handle to an RMnode to contain the generated contours
    (modified).
 @aend

 @dstart

 Computes contour lnes in a structured 2D dataset.  An arbitrary
 number of contour levels may be specified and each contour level may
 have an associated color.  The algorithm used will find all contours,
 even if disconnected and sparse.  The resulting geometry is added to
 the RMnode.

 No status is currently returned to the caller.

 @dend
 * ----------------------------------------------------
 */
void
rmvJ3SliceContour (RMvertex3D (*appgridfunc)(int i, int j, void *usrData),
		   void *appGridUsrData,
		   float (*appdatafunc)(int i, int j, void *usrData),
		   void *appDataUsrData,
		   int axis_offset_enum,
		   int iusize,
		   int ivsize,
		   int nlevels,
		   float *levels,
		   RMcolor4D *color_list,
		   RMenum linewidth_enum,
		   RMenum linestyle_enum,
		   RMnode *n)
{
    int          i, j, k;
    int          maxverts, nverts = 0;
    int          usize, vsize;
    int          opcode;
    float        d00, d10, d01, d11;
    RMvertex3D  *v;
    RMcolor4D   *c;
    RMprimitive *q;

    /* need assertions */
    q = rmPrimitiveNew(RM_LINES);

    usize = iusize;
    vsize = ivsize;
    maxverts = usize * vsize * 4 * nlevels;
    v = rmVertex3DNew(maxverts);

    /* NOTE: we have the limitation that we support only RGB colors (no alpha) */
    if (color_list != NULL)
	c = rmColor4DNew(maxverts);
    else
	c = NULL;
    
    /* for each quad in the mesh, evaluate whether or not there is a threshold crossing for each of the nlevels */
    for (j = 0; j < (vsize - 1); j++)
    {
	for (i = 0; i < (usize - 1); i++)
	{
	    d00 = (*appdatafunc)(i, j, appDataUsrData);
	    d10 = (*appdatafunc)((i + 1), j, appDataUsrData);
	    d01 = (*appdatafunc)(i, (j + 1), appDataUsrData);
	    d11 = (*appdatafunc)((i + 1), (j + 1), appDataUsrData);
	    
	    for (k = 0; k < nlevels; k++)
	    {
		opcode = 0;
		if (threshcross(d00, (levels + k)))
		    opcode |= 1;
		if (threshcross(d10, (levels + k)))
		    opcode |= 2;
		if (threshcross(d11, (levels + k)))
		    opcode |= 4;
		if (threshcross(d01, (levels + k)))
		    opcode |= 8;
		switch (opcode)
		   {
		   case 0:
		   case 0xF:
		      break;  /* no threshold crossings here */
		      
		   case 1:
		      do_contour(i, (j + 1), i, j, (i + 1), j, i, j,
				 appgridfunc, appGridUsrData,
				 appdatafunc, appDataUsrData,
				 usize, vsize, axis_offset_enum, (levels + k),
				 v, &nverts,
				 (c == NULL) ? NULL : (color_list + k), c);
		      break;
		      
		   case 2:
		      do_contour(i, j, (i + 1), j, (i + 1), (j + 1), (i + 1), j,
				 appgridfunc, appGridUsrData,
				 appdatafunc, appDataUsrData,
				 usize, vsize, axis_offset_enum, (levels + k),
				 v, &nverts,
				 (c == NULL) ? NULL : (color_list + k), c);
		      break;
		      
		   case 3:
		      do_contour(i, (j + 1), i, j, (i + 1), (j + 1), (i + 1), j,
				 appgridfunc, appGridUsrData,
				 appdatafunc, appDataUsrData,
				 usize, vsize, axis_offset_enum, (levels + k),
				 v, &nverts,
				 (c == NULL) ? NULL : (color_list + k), c);
		      break;
		      
		   case 4:
		      do_contour(i, (j + 1), (i + 1), (j + 1), (i + 1), j, (i + 1), (j + 1),
				 appgridfunc, appGridUsrData,
				 appdatafunc, appDataUsrData,
				 usize, vsize, axis_offset_enum, (levels + k),
				 v, &nverts,
				 (c == NULL) ? NULL : (color_list + k), c);
		      break;
		      
		   case 5:
		      
		      do_contour(i, (j + 1), (i + 1), (j + 1), (i + 1), j, (i + 1), (j + 1),
				 appgridfunc, appGridUsrData,
				 appdatafunc, appDataUsrData,
				 usize, vsize, axis_offset_enum, (levels + k),
				 v, &nverts,
				 (c == NULL) ? NULL : (color_list + k), c);
		      
		      do_contour(i, (j + 1), i, j, (i + 1), j, i, j,
				 appgridfunc, appGridUsrData,
				 appdatafunc, appDataUsrData,
				 usize, vsize, axis_offset_enum, (levels + k),
				 v, &nverts,
				 (c == NULL) ? NULL : (color_list + k), c);
		      break;
		      
		   case 6:
		      do_contour(i, j, (i + 1), j, i, (j + 1), (i + 1), (j + 1),
				 appgridfunc, appGridUsrData,
				 appdatafunc, appDataUsrData,
				 usize, vsize, axis_offset_enum, (levels + k),
				 v, &nverts,
				 (c==NULL) ? NULL : (color_list + k), c);
		      break;
		      
		   case 7:
		      do_contour(i, (j + 1), (i + 1), (j + 1), i, (j + 1), i, j,
				 appgridfunc, appGridUsrData,
				 appdatafunc, appDataUsrData,
				 usize, vsize, axis_offset_enum, (levels + k),
				 v, &nverts,
				 (c == NULL) ? NULL : (color_list + k), c);
		      break;
		      
		   case 8:
		      do_contour(i, j, i, (j + 1), (i + 1), (j + 1), i, (j + 1),
				 appgridfunc, appGridUsrData,
				 appdatafunc, appDataUsrData,
				 usize, vsize, axis_offset_enum, (levels + k),
				 v, &nverts,
				 (c == NULL) ? NULL : (color_list + k), c);
		      break;
		      
		   case 9:
		      do_contour((i + 1), (j + 1), i, (j + 1), (i + 1), j, i, j,
				 appgridfunc, appGridUsrData,
				 appdatafunc, appDataUsrData,
				 usize, vsize, axis_offset_enum, (levels + k),
				 v, &nverts,
				 (c == NULL) ? NULL : (color_list + k), c);
		      break;
		      
		   case 0xA:
		      do_contour((i + 1), (j + 1), i, (j + 1), (i + 1), (j + 1), (i + 1), j,
				 appgridfunc, appGridUsrData,
				 appdatafunc, appDataUsrData,
				 usize, vsize, axis_offset_enum, (levels + k),
				 v, &nverts,
				 (c == NULL) ? NULL : (color_list + k), c);
		      
		      do_contour(i, j, i, (j + 1), i, j, (i + 1), j,
				 appgridfunc, appGridUsrData,
				 appdatafunc, appDataUsrData,
				 usize, vsize, axis_offset_enum, (levels + k),
				 v, &nverts,
				 (c == NULL) ? NULL : (color_list + k), c);
		      break;
		      
		   case 0xB:
		      do_contour((i + 1), (j + 1), i, (j + 1), (i + 1), (j + 1), (i + 1), j,
				 appgridfunc, appGridUsrData,
				 appdatafunc, appDataUsrData,
				 usize, vsize, axis_offset_enum, (levels + k),
				 v, &nverts,
				 (c == NULL) ? NULL : (color_list + k),c);
		      break;
		      
		   case 0xC:
		      do_contour(i, j, i, (j + 1), (i + 1), j, (i + 1), (j + 1),
				 appgridfunc, appGridUsrData,
				 appdatafunc, appDataUsrData,
				 usize, vsize, axis_offset_enum, (levels + k),
				 v, &nverts,
				 (c == NULL) ? NULL : (color_list + k), c);
		      break;
		      
		   case 0xD:
		      do_contour((i + 1), j, (i + 1), (j + 1), (i + 1), j, i, j,
				 appgridfunc, appGridUsrData,
				 appdatafunc, appDataUsrData,
				 usize, vsize, axis_offset_enum, (levels + k),
				 v, &nverts,
				 (c == NULL) ? NULL : (color_list + k), c);
		      break;
		      
		   case 0xE:
		      do_contour(i, j, (i + 1), j, i, j, i, (j + 1),
				 appgridfunc, appGridUsrData,
				 appdatafunc, appDataUsrData,
				 usize, vsize, axis_offset_enum, (levels + k),
				 v, &nverts,
				 (c == NULL) ? NULL : (color_list + k),c);
		      
		      break;
		      
		   default: /* bogus case */
		      break;
		      
		   } /* end switch */
	    } /* end loop over k (levels) */
	} /* end loop over i (usize) */
    } /* end loop over j (vsize) */

    rmNodeSetLineWidth(n, linewidth_enum);
    rmNodeSetLineStyle(n, linestyle_enum);
	
    rmPrimitiveSetVertex3D(q, nverts, v, RM_COPY_DATA, NULL);

    if (c)
    {
	rmPrimitiveSetColor4D(q, nverts, c, RM_COPY_DATA, NULL);
	rmColor4DDelete(c);
    }
    
    /* now, add the new primitive onto the node */
    rmNodeAddPrimitive(n, q);

    /* compute the bounding box and center point for the node */
    rmNodeComputeBoundingBox(n);
    {
        RMvertex3D bmin, bmax, center;

	rmPointMinMax((float *)v, nverts, 3, sizeof(RMvertex3D), &bmin,&bmax);
	rmNodeSetBoundingBox(n, &bmin, &bmax);
	center.x = bmin.x + (0.5 * (bmax.x - bmin.x));
	center.y = bmin.y + (0.5 * (bmax.y - bmin.y));
	center.z = bmin.z + (0.5 * (bmax.z - bmin.z));
	rmNodeSetCenter(n, &center);
    }
    rmVertex3DDelete (v);
}


/* PRIVATE */
int
threshcross (float d0,
	     float *lev)
{
    register float rlev;

    rlev = *lev;

    if (rlev >= d0)
	return(1);
    else
	return(0);
}


/* PRIVATE */
void
do_contour (int i1,
	    int j1,		/* below thresh index*/
	    int i2,
	    int j2,		/* above thresh index */
	    int i3,
	    int j3,		/* below thresh index */
	    int i4,
	    int j4,		/* above thresh index*/
	    RMvertex3D (*appgridfunc)(int i, int j, void *usrData),
	    void *appGridUsrData,
	    float (*appdatafunc)(int i, int j, void *usrData),
	    void *appDataUsrData,
	    int usize,
	    int vsize,
	    int offset_mask,
	    float *lev,
	    RMvertex3D *varray,
	    int *nverts,
	    RMcolor4D *this_color,
	    RMcolor4D *carray)
{
    float      d1, d2, d3, d4;
    double     u;
    RMvertex3D v, p1, p2, p3, p4;

    d1 = (*appdatafunc)(i1, j1, appDataUsrData);
    d2 = (*appdatafunc)(i2, j2, appDataUsrData);
    d3 = (*appdatafunc)(i3, j3, appDataUsrData);
    d4 = (*appdatafunc)(i4, j4, appDataUsrData);

    p1 = (*appgridfunc)(i1, j1, appGridUsrData);
    p2 = (*appgridfunc)(i2, j2, appGridUsrData);
    p3 = (*appgridfunc)(i3, j3, appGridUsrData);
    p4 = (*appgridfunc)(i4, j4, appGridUsrData);

#if 0
    /* old way - read from an array */
    d1 = *(data + i1 + (j1 * usize));
    d2 = *(data + i2 + (j2 * usize));
    d3 = *(data + i3 + (j3 * usize));
    d4 = *(data + i4 + (j4 * usize));

    p1.x = xcoords[i1 + (j1 * usize)];
    p1.y = ycoords[i1 + (j1 * usize)];
    p1.z = zcoords[i1 + (j1 * usize)];

    p2.x = xcoords[i2 + (j2 * usize)];
    p2.y = ycoords[i2 + (j2 * usize)];
    p2.z = zcoords[i2 + (j2 * usize)];

    p3.x = xcoords[i3 + (j3 * usize)];
    p3.y = ycoords[i3 + (j3 * usize)];
    p3.z = zcoords[i3 + (j3 * usize)];

    p4.x = xcoords[i4 + (j4 * usize)];
    p4.y = ycoords[i4 + (j4 * usize)];
    p4.z = zcoords[i4 + (j4 * usize)];
#else
    /* foil compiler warning */
    usize = vsize = 0;
#endif

    /* do 1st vertex of contour line */
    u = (*lev - d1) / (d2 - d1);

    v.x = p1.x + (u * (p2.x - p1.x));
    if (offset_mask & RMV_XAXIS_OFFSET)
	v.x += *lev;
    
    v.y = p1.y + (u * (p2.y - p1.y));
    if (offset_mask & RMV_YAXIS_OFFSET)
	v.y += *lev;
    
    v.z = p1.z + (u * (p2.z - p1.z));
    if (offset_mask & RMV_ZAXIS_OFFSET)
	v.z += *lev;

    VCOPY(&v, varray + (*nverts));

    if (carray != NULL)
	VCOPY(this_color, carray + (*nverts));
    *nverts += 1;
    
    /* do 2nd vertex of contour line */
    u = (*lev - d3) / (d4 - d3);

    v.x = p3.x + (u * (p4.x - p3.x));
    if (offset_mask & RMV_XAXIS_OFFSET)
	v.x += *lev;
    
    v.y = p3.y + (u * (p4.y - p3.y));
    if (offset_mask & RMV_YAXIS_OFFSET)
	v.y += *lev;
    
    v.z = p3.z + (u * (p4.z - p3.z));
    if (offset_mask & RMV_ZAXIS_OFFSET)
	v.z += *lev;

    VCOPY(&v, varray + (*nverts));

    if (carray != NULL)
	VCOPY(this_color, carray + (*nverts));
    *nverts += 1;
}
/* EOF */
