/* GNUPLOT - hidden3d.c */

/*[
 * Copyright 1986 - 1993, 1998, 1999, 2004   Thomas Williams, Colin Kelley
 *
 * Permission to use, copy, and distribute this software and its
 * documentation for any purpose with or without fee is hereby granted,
 * provided that the above copyright notice appear in all copies and
 * that both that copyright notice and this permission notice appear
 * in supporting documentation.
 *
 * Permission to modify the software is granted, but not the right to
 * distribute the complete modified source code.  Modifications are to
 * be distributed as patches to the released version.  Permission to
 * distribute binaries produced by compiling modified sources is granted,
 * provided you
 *   1. distribute the corresponding source modifications from the
 *    released version in the form of a patch file along with the binaries,
 *   2. add special version identification to distinguish your version
 *    in addition to the base release version number,
 *   3. provide your name and address as the primary contact for the
 *    support of your modified version, and
 *   4. retain our contact information in regard to use of the base
 *    software.
 * Permission to distribute the released version of the source code along
 * with corresponding source modifications in the form of a patch file is
 * granted with same provisions 2 through 4 for binary distributions.
 *
 * This software is provided "as is" without express or implied warranty
 * to the extent permitted by applicable law.
]*/


/*
 * 1999 Hans-Bernhard Broeker (Broeker@physik.rwth-aachen.de)
 *
 * Major rewrite, affecting just about everything
 *
 */
#include <math.h>
#include "hidden3d.h"
#include "color.h"
#include "pm3d.h"
#include "alloc.h"
#include "axis.h"
#include "command.h"
#include "dynarray.h"
#include "graph3d.h"
#include "tables.h"
#include "term_api.h"
#include "util.h"
#include "util3d.h"


/*************************/
/* Configuration section */
/*************************/

/* If this module is compiled with HIDDEN3D_GRIDBOX = 1 defined, it
 * will store the information about {x|y}{min|max} in an other
 * (additional!) form: a bit mask, with each bit representing one
 * horizontal or vertical strip of the screen. The bits for strips a
 * polygon spans are set to one. This allows to test for xy overlap
 * of an edge with a polygon simply by comparing bit patterns.  */
#ifndef HIDDEN3D_GRIDBOX
#define HIDDEN3D_GRIDBOX 0
#endif

/* HBB 19991204: new code started to finally implement a spatially
 * ordered data structure to store the polygons in. This is meant to
 * speed up the HLR process. Before, the hot spot of hidden3d was the
 * loop in in_front, where by far most of the polygons are rejected by
 * the first test, already. The idea is to _not_ to loop over all
 * those polygons far away from the edge under consideration, in the
 * first place. Instead, store the polygons in an xy grid of lists,
 * so we can select a sample of these lists to test a given edge
 * against. */
#ifndef HIDDEN3D_QUADTREE
#define HIDDEN3D_QUADTREE 0
#endif
#if HIDDEN3D_QUADTREE && HIDDEN3D_GRIDBOX
# warning HIDDEN3D_QUADTREE & HIDDEN3D_GRIDBOX do not work together, sensibly!
#endif

/* If you don't want the color-distinction between the
 * 'top' and 'bottom' sides of the surface, like I do, then just compile
 * with -DBACKSIDE_LINETYPE_OFFSET = 0. */
#ifndef BACKSIDE_LINETYPE_OFFSET
# define BACKSIDE_LINETYPE_OFFSET 1
#endif

/* This #define lets you choose if the diagonals that
 * divide each original quadrangle in two triangles will be drawn
 * visible or not: do draw them, define it to be 7L, otherwise let be
 * 3L */
#ifndef TRIANGLE_LINESDRAWN_PATTERN
# define TRIANGLE_LINESDRAWN_PATTERN 3L
#endif

/* Handle out-of-range or undefined points. Compares the maximum
 * marking (0=inrange, 1=outrange, 2=undefined) of the coordinates of
 * a vertex to the value #defined here. If not less, the vertex is
 * rejected, and all edges that hit it, as well. NB: if this is set to
 * anything above 1, gnuplot may crash with a floating point exception
 * in hidden3d. You get what you asked for ... */
#ifndef HANDLE_UNDEFINED_POINTS
# define HANDLE_UNDEFINED_POINTS 1
#endif
/* Symbolic value for 'do not handle Undefined Points specially' */
#define UNHANDLED (UNDEFINED+1)

/* If both subtriangles of a quad were cancelled, try if using the
 * other diagonal is better. This only makes a difference if exactly
 * one vertex of the quad is unusable, and that one is on the 'usual'
 * tried diagonal. In such a case, the border of the hole in the
 * surface will be less rough than with the previous method, as the
 * border follows the undefined region as close as it can. */
#ifndef SHOW_ALTERNATIVE_DIAGONAL
# define SHOW_ALTERNATIVE_DIAGONAL 1
#endif

/* If the two triangles in a quad are both drawn, and they show
 * different sides to the user (the quad is 'bent over'), then it's
 * preferrable to force the diagonal being visible to avoid other
 * parts of the scene being obscured by a line the user can't
 * see. This avoids unnecessary user surprises. */
#ifndef HANDLE_BENTOVER_QUADRANGLES
# define HANDLE_BENTOVER_QUADRANGLES 1
#endif

/* The actual configuration is stored in these variables, modifiable
 * at runtime through 'set hidden3d' options */
static int hiddenBacksideLinetypeOffset = BACKSIDE_LINETYPE_OFFSET;
static long hiddenTriangleLinesdrawnPattern = TRIANGLE_LINESDRAWN_PATTERN;
static int hiddenHandleUndefinedPoints = HANDLE_UNDEFINED_POINTS;
static int hiddenShowAlternativeDiagonal = SHOW_ALTERNATIVE_DIAGONAL;
static int hiddenHandleBentoverQuadrangles = HANDLE_BENTOVER_QUADRANGLES;


/**************************************************************/
/**************************************************************
 * The 'real' code begins, here.                              *
 *                                                            *
 * first: types and global variables                          *
 **************************************************************/
/**************************************************************/

/* precision of calculations in normalized space. Coordinates closer to
 * each other than an absolute difference of EPSILON are considered
 * equal, by some of the routines in this module. */
#define EPSILON 1e-5

/* Some inexact operations: == , > , >=, sign() */
#define EQ(X,Y)  (fabs( (X) - (Y) ) < EPSILON)	/* X == Y */
#define GR(X,Y)  ((X) >  (Y) + EPSILON)		/* X >  Y */
#define GE(X,Y)  ((X) >= (Y) - EPSILON)		/* X >= Y */
#define SIGN(X)  ( ((X)<-EPSILON) ? -1: ((X)>EPSILON) )

/* A plane equation, stored as a four-element vector. The equation
 * itself is: x*p[0]+y*p[1]+z*p[2]+1*p[3]=0 */
typedef coordval t_plane[4];

/* One edge of the mesh. The edges are (currently) organized into a
 * linked list as a method of traversing them back-to-front. */
typedef struct edge {
    long v1, v2;		/* the vertices at either end */
    int style;			/* linetype index */
    struct lp_style_type *lp;	/* line/point style attributes */
    long next;			/* index of next edge in z-sorted list */
} edge;
typedef edge GPHUGE *p_edge;


/* One polygon (actually: always a triangle) of the surface
 * mesh(es). */
#define POLY_NVERT 3
typedef struct polygon {
    long vertex[POLY_NVERT];    /* The vertices (indices on vlist) */
    /* min/max in all three directions */
    coordval xmin, xmax, ymin, ymax, zmin, zmax;
    t_plane plane;		/* the plane coefficients */
    TBOOLEAN frontfacing;	/* is polygon facing front- or backwards? */
#if ! HIDDEN3D_QUADTREE
    long next;			/* index of next polygon in z-sorted list */
#endif
#if HIDDEN3D_GRIDBOX
    unsigned long xbits;	/* x coverage mask of bounding box */
    unsigned long ybits;	/* y coverage mask of bounding box */
#endif
} polygon;
typedef polygon GPHUGE *p_polygon;

#if HIDDEN3D_GRIDBOX
# define UINT_BITS (CHAR_BIT * sizeof(unsigned int))
# define COORD_TO_BITMASK(x,shift)							\
  (~0U << (unsigned int) ((((x) / surface_scale) + 1.0) / 2.0 * UINT_BITS + (shift)))
# define CALC_BITRANGE(range_min, range_max)				 \
  ((~COORD_TO_BITMASK((range_max), 1)) & COORD_TO_BITMASK(range_min, 0))
#endif

/* Enumeration of possible types of line, for use with the
 * store_edge() function. Influences the position in the grid the
 * second vertex will be put to, relative to the one that is passed
 * in, as another argument to that function. edir_none is for
 * single-pixel 'blind' edges, which exist only to facilitate output
 * of 'points' style splots.
 *
 * Directions are interpreted in a pseudo-geographical coordinate
 * system of the data grid: within the isoline, we count from left to
 * right (west to east), and the isolines themselves are counted from
 * top to bottom, described as north and south. */
typedef enum edge_direction {
    edir_west, edir_north,
    edir_NW, edir_NE,
    edir_impulse, edir_point,
    edir_vector
} edge_direction;

/* direction into which the polygon is facing (the corner with the
 * right angle, inside the mesh, that is). The reference identifiying
 * the whole cell is always the lower right, i.e. southeast one. */
typedef enum polygon_direction {
    pdir_NE, pdir_SE, pdir_SW, pdir_NW
} polygon_direction;

/* Three dynamical arrays that describe what we have to plot: */
static dynarray vertices, edges, polygons;

/* convenience #defines to make the generic vector useable as typed arrays */
#define vlist ((p_vertex) vertices.v)
#define plist ((p_polygon) polygons.v)
#define elist ((p_edge) edges.v)

static long pfirst;		/* first polygon in zsorted chain*/
static long efirst;		/* first edges in zsorted chain */

#if HIDDEN3D_QUADTREE
/* HBB 20000716: spatially oriented hierarchical data structure to
 * store polygons in. For now, it's a simple xy grid of z-sorted
 * lists. A single polygon can appear in several lists, if it spans
 * cell borders */
typedef struct qtreelist {
    long p;			/* the polygon */
    long next;			/* next element in this chain */
} qtreelist;
typedef qtreelist GPHUGE *p_qtreelist;

/* the number of cells in x and y direction: */
# ifndef QUADTREE_GRANULARITY
#  define QUADTREE_GRANULARITY 10
# endif
/* indices of the heads of all the cells' chains: */
static long quadtree[QUADTREE_GRANULARITY][QUADTREE_GRANULARITY];
/* and a macro to calculate the cells' position in that array: */
#define COORD_TO_TREECELL(x)								\
    ((unsigned int)(((((x) / surface_scale) + 1.0) / 2.0) * QUADTREE_GRANULARITY))

/* the dynarray to actually store all that stuff in: */
static dynarray qtree;
#define qlist ((p_qtreelist) qtree.v)
#endif /* HIDDEN3D_QUADTREE*/

/* Prototypes for internal functions of this module. */
static long int store_vertex (struct coordinate GPHUGE *point,
				      lp_style_type *lp_style, TBOOLEAN color_from_column);
static long int make_edge (long int vnum1, long int vnum2,
				   struct lp_style_type *lp,
				   int style, int next);
static long int store_edge (long int vnum1, edge_direction direction,
				    long int crvlen, struct lp_style_type *lp,
				    int style);
static GP_INLINE double eval_plane_equation (t_plane p, p_vertex v);
static long int store_polygon (long int vnum1,
				       polygon_direction direction,
				       long int crvlen);
static void color_edges (long int new_edge, long int old_edge,
				 long int new_poly, long int old_poly,
				 int style_above, int style_below);
static void build_networks (struct surface_points * plots,
				    int pcount);
int compare_edges_by_zmin (SORTFUNC_ARGS p1, SORTFUNC_ARGS p2);
int compare_polys_by_zmax (SORTFUNC_ARGS p1, SORTFUNC_ARGS p2);
static void sort_edges_by_z (void);
static void sort_polys_by_z (void);
static TBOOLEAN get_plane (p_polygon p, t_plane plane);
static long split_line_at_ratio (long int vnum1, long int vnum2, double w);
static GP_INLINE double area2D (p_vertex v1, p_vertex v2,
					p_vertex v3);
static void draw_vertex (p_vertex v);
static GP_INLINE void draw_edge (p_edge e, p_vertex v1, p_vertex v2);
static GP_INLINE void handle_edge_fragment (long int edgenum,
						    long int vnum1,
						    long int vnum2,
						    long int firstpoly);
static int in_front (long int edgenum,
			     long int vnum1, long int vnum2,
			     long int *firstpoly);


/* Set the options for hidden3d. To be called from set.c, when the
 * user has begun a command with 'set hidden3d', to parse the rest of
 * that command */
void
set_hidden3doptions()
{
    int tmp;

    while (!END_OF_COMMAND) {
	switch (lookup_table(&set_hidden3d_tbl[0], c_token)) {
	case S_HI_DEFAULTS:
	    /* reset all parameters to defaults */
	    reset_hidden3doptions();
	    c_token++;
	    if (!END_OF_COMMAND)
		int_error(c_token,
			  "No further options allowed after 'defaults'");
	    return;
	    break;
	case S_HI_OFFSET:
	    c_token++;
	    hiddenBacksideLinetypeOffset = int_expression();
	    c_token--;
	    break;
	case S_HI_NOOFFSET:
	    hiddenBacksideLinetypeOffset = 0;
	    break;
	case S_HI_TRIANGLEPATTERN:
	    c_token++;
	    hiddenTriangleLinesdrawnPattern = int_expression();
	    c_token--;
	    break;
	case S_HI_UNDEFINED:
	    c_token++;
	    tmp = int_expression();
	    if (tmp <= 0 || tmp > UNHANDLED)
		tmp = UNHANDLED;
	    hiddenHandleUndefinedPoints = tmp;
	    c_token--;
	    break;
	case S_HI_NOUNDEFINED:
	    hiddenHandleUndefinedPoints = UNHANDLED;
	    break;
	case S_HI_ALTDIAGONAL:
	    hiddenShowAlternativeDiagonal = 1;
	    break;
	case S_HI_NOALTDIAGONAL:
	    hiddenShowAlternativeDiagonal = 0;
	    break;
	case S_HI_BENTOVER:
	    hiddenHandleBentoverQuadrangles = 1;
	    break;
	case S_HI_NOBENTOVER:
	    hiddenHandleBentoverQuadrangles = 0;
	    break;
	case S_HI_INVALID:
	    int_error(c_token, "No such option to hidden3d (or wrong order)");
	default:
	    break;
	}
	c_token++;
    }
}

void
show_hidden3doptions()
{
    fprintf(stderr,"\
\t  Back side of surfaces has linestyle offset of %d\n\
\t  Bit-Mask of Lines to draw in each triangle is %ld\n\
\t  %d: ",
	    hiddenBacksideLinetypeOffset, hiddenTriangleLinesdrawnPattern,
	    hiddenHandleUndefinedPoints);

    switch (hiddenHandleUndefinedPoints) {
    case OUTRANGE:
	fputs("Outranged and undefined datapoints are omitted from the surface.\n",
	      stderr);
	break;
    case UNDEFINED:
	fputs("Only undefined datapoints are omitted from the surface.\n",
	      stderr);
	break;
    case UNHANDLED:
	fputs("Will not check for undefined datapoints (may cause crashes).\n",
	      stderr);
	break;
    default:
	fputs("Value stored for undefined datapoint handling is illegal!!!\n",
	      stderr);
	break;
    }

    fprintf(stderr,"\
\t  Will %suse other diagonal if it gives a less jaggy outline\n\
\t  Will %sdraw diagonal visibly if quadrangle is 'bent over'\n",
	    hiddenShowAlternativeDiagonal ? "" : "not ",
	    hiddenHandleBentoverQuadrangles ? "" : "not ");
}

/* Implements proper 'save'ing of the new hidden3d options... */
void
save_hidden3doptions(FILE *fp)
{
    if (!hidden3d) {
	fputs("unset hidden3d\n", fp);
	return;
    }
    fprintf(fp, "set hidden3d offset %d trianglepattern %ld undefined %d %saltdiagonal %sbentover\n",
	    hiddenBacksideLinetypeOffset,
	    hiddenTriangleLinesdrawnPattern,
	    hiddenHandleUndefinedPoints,
	    hiddenShowAlternativeDiagonal ? "" : "no",
	    hiddenHandleBentoverQuadrangles ? "" : "no");
}

/* Initialize the necessary steps for hidden line removal and
   initialize global variables. */
void
init_hidden_line_removal()
{
    /* Check for some necessary conditions to be set elsewhere: */
    /* HandleUndefinedPoints mechanism depends on these: */
    assert(OUTRANGE == 1);
    assert(UNDEFINED == 2);

    /* Re-mapping of this value makes the test easier in the critical
     * section */
    if (hiddenHandleUndefinedPoints < OUTRANGE)
	hiddenHandleUndefinedPoints = UNHANDLED;

    init_dynarray(&vertices, sizeof(vertex), 100, 100);
    init_dynarray(&edges, sizeof(edge), 100, 100);
    init_dynarray(&polygons, sizeof(polygon), 100, 100);
#if HIDDEN3D_QUADTREE
    init_dynarray(&qtree, sizeof(qtreelist), 100, 100);
#endif

}

/* Reset the hidden line data to a fresh start. */
void
reset_hidden_line_removal()
{
    vertices.end = 0;
    edges.end = 0;
    polygons.end = 0;
#if HIDDEN3D_QUADTREE
    qtree.end = 0;
#endif
}


/* Terminates the hidden line removal process.                  */
/* Free any memory allocated by init_hidden_line_removal above. */
void
term_hidden_line_removal()
{
    free_dynarray(&polygons);
    free_dynarray(&edges);
    free_dynarray(&vertices);
#if HIDDEN3D_QUADTREE
    free_dynarray(&qtree);
#endif
}


#if 0 /* UNUSED ! */
/* Do we see the top or bottom of the polygon, or is it 'on edge'? */
#define GET_SIDE(vlst,csign)						\
do {									\
    double ctmp =							\
	vlist[vlst[0]].x * (vlist[vlst[1]].y - vlist[vlst[2]].y) +	\
	vlist[vlst[1]].x * (vlist[vlst[2]].y - vlist[vlst[0]].y) +	\
	vlist[vlst[2]].x * (vlist[vlst[0]].y - vlist[vlst[1]].y);	\
    csign = SIGN (ctmp);						\
} while (0)
#endif /* UNUSED */

static long int
store_vertex (
    struct coordinate GPHUGE * point,
    lp_style_type *lp_style,
    TBOOLEAN color_from_column)
{
    p_vertex thisvert = (struct vertex *)nextfrom_dynarray(&vertices);

    thisvert->lp_style = lp_style;
    if ((int) point->type >= hiddenHandleUndefinedPoints) {
	FLAG_VERTEX_AS_UNDEFINED(*thisvert);
	return (-1);
    }
    map3d_xyz(point->x, point->y, point->z, thisvert);
    if (color_from_column) {
	thisvert->real_z = point->CRD_COLOR;
	thisvert->lp_style->pm3d_color.lt = LT_COLORFROMCOLUMN;
    } else
	thisvert->real_z = point->z;

#ifdef HIDDEN3D_VAR_PTSIZE
    /* Store pointer back to original point */
    /* Needed to support variable pointsize */
    thisvert->original = point;
#endif
	
    return (thisvert - vlist);
}

/* A part of store_edge that does the actual storing. Used by
 * in_front(), as well, so I separated it out. */
static long int make_edge(
    long vnum1, long vnum2,
    struct lp_style_type *lp,
    int style,
    int next)
{
    p_edge thisedge = (struct edge *)nextfrom_dynarray(&edges);
    p_vertex v1 = vlist + vnum1;
    p_vertex v2 = vlist + vnum2;

    /* ensure z ordering inside each edge */
    if (v1->z >= v2->z) {
	thisedge->v1 = vnum1;
	thisedge->v2 = vnum2;
    } else {
	thisedge->v1 = vnum2;
	thisedge->v2 = vnum1;
    }

    thisedge->style = style;
    thisedge->lp = lp;
    thisedge->next = next;

    return thisedge - elist;
}

/* store the edge from vnum1 to vnum2 into the edge list. Ensure that
 * the vertex with higher z is stored in v1, to ease sorting by zmax */
static long int
store_edge(
    long int vnum1,
    edge_direction direction,
    long int crvlen,
    struct lp_style_type *lp,
    int style)
{
    p_vertex v1 = vlist + vnum1;
    p_vertex v2 = NULL;		/* just in case: initialize... */
    long int vnum2;
    unsigned int drawbits = (0x1 << direction);

    switch (direction) {
    case edir_vector:
	v2 = v1 + 1;
	drawbits = 0;
	break;
    case edir_west:
	v2 = v1 - 1;
	break;
    case edir_north:
	v2 = v1 - crvlen;
	break;
    case edir_NW:
	v2 = v1 - crvlen - 1;
	break;
    case edir_NE:
	v2 = v1 - crvlen;
	v1 -= 1;
	drawbits >>= 1;		/* altDiag is handled like normal NW one */
	break;
    case edir_impulse:
	v2 = v1 - 1;
	drawbits = 0;		/* don't care about the triangle pattern */
	break;
    case edir_point:
	v2 = v1;
	drawbits = 0;		/* nothing to draw, but disable check */
	break;
    }

    vnum2 = v2 - vlist;

    if (VERTEX_IS_UNDEFINED(*v1)
	|| VERTEX_IS_UNDEFINED(*v2)) {
	return -2;
    }

    if (drawbits &&		/* no bits set: 'blind' edge --> no test! */
	! (hiddenTriangleLinesdrawnPattern & drawbits)
	)
	style = -3;

    return make_edge(vnum1, vnum2, lp, style, -1);
}


/* Calculate the normal equation coefficients of the plane of polygon
 * 'p'. Uses is the 'signed projected area' method. Its benefit is
 * that it doesn't rely on only three of the vertices of 'p', as the
 * naive cross product method does. */
static TBOOLEAN
get_plane(p_polygon poly, t_plane plane)
{
    int i;
    p_vertex v1, v2;
    double x, y, z, s;
    TBOOLEAN frontfacing=TRUE;

    /* calculate the signed areas of the polygon projected onto the
     * planes x=0, y=0 and z=0, respectively. The three areas form
     * the components of the plane's normal vector: */
    v1 = vlist + poly->vertex[POLY_NVERT - 1];
    v2 = vlist + poly->vertex[0];
    plane[0] = (v1->y - v2->y) * (v1->z + v2->z);
    plane[1] = (v1->z - v2->z) * (v1->x + v2->x);
    plane[2] = (v1->x - v2->x) * (v1->y + v2->y);
    for (i = 1; i < POLY_NVERT; i++) {
	v1 = v2;
	v2 = vlist + poly->vertex[i];
	plane[0] += (v1->y - v2->y) * (v1->z + v2->z);
	plane[1] += (v1->z - v2->z) * (v1->x + v2->x);
	plane[2] += (v1->x - v2->x) * (v1->y + v2->y);
    }

    /* Normalize the resulting normal vector */
    s = sqrt(plane[0] * plane[0] + plane[1] * plane[1] + plane[2] * plane[2]);

    if (GE(0.0, s)) {
	/* The normal vanishes, i.e. the polygon is degenerate. We build
	 * another vector that is orthogonal to the line of the polygon */
	v1 = vlist + poly->vertex[0];
	for (i = 1; i < POLY_NVERT; i++) {
	    v2 = vlist + poly->vertex[i];
	    if (!V_EQUAL(v1, v2))
		break;
	}

	/* build (x,y,z) that should be linear-independant from <v1, v2> */
	x = v1->x;
	y = v1->y;
	z = v1->z;
	if (EQ(y, v2->y))
	    y += 1.0;
	else
	    x += 1.0;

	/* Re-do the signed area computations */
	plane[0] = v1->y * (v2->z - z) + v2->y * (z - v1->z) + y * (v1->z - v2->z);
	plane[1] = v1->z * (v2->x - x) + v2->z * (x - v1->x) + z * (v1->x - v2->x);
	plane[2] = v1->x * (v2->y - y) + v2->x * (y - v1->y) + x * (v1->y - v2->y);
	s = sqrt(plane[0] * plane[0] + plane[1] * plane[1] + plane[2] * plane[2]);
    }

    /* ensure that normalized c is > 0 */
    if (plane[2] < 0.0) {
	s *= -1.0;
	frontfacing = FALSE;
    }

    plane[0] /= s;
    plane[1] /= s;
    plane[2] /= s;

    /* Now we have the normalized normal vector, insert one of the
     * vertices into the equation to get 'd'. For an even better result,
     * an average over all the vertices might be used */
    plane[3] = -plane[0] * v1->x - plane[1] * v1->y - plane[2] * v1->z;

    return frontfacing;
}


/* Evaluate the plane equation represented a four-vector for the given
 * vector. For points in the plane, this should result in values ==0.
 * < 0 is 'away' from the polygon, > 0 is infront of it */
static GP_INLINE double
eval_plane_equation(t_plane p, p_vertex v)
{
    return (p[0]*v->x + p[1]*v->y + p[2]*v->z + p[3]);
}


/* Build the data structure for this polygon. The return value is the
 * index of the newly generated polygon. This is memorized for access
 * to polygons in the previous isoline, from the next-following
 * one. */
static long int
store_polygon(long vnum1, polygon_direction direction, long crvlen)
{
    long int v[POLY_NVERT];
    p_vertex v1, v2, v3;
    p_polygon p;

    switch (direction) {
    case pdir_NE:
	v[0] = vnum1;
	v[2] = vnum1 - crvlen;
	v[1] = v[2] - 1;
	break;
    case pdir_SW:
	/* triangle points southwest, here */
	v[0] = vnum1;
	v[1] = vnum1 - 1;
	v[2] = v[1] - crvlen;
	break;
    case pdir_SE:
	/* alt-diagonal, case 1: southeast triangle: */
	v[0] = vnum1;
	v[2] = vnum1 - crvlen;
	v[1] = vnum1 - 1;
	break;
    case pdir_NW:
	v[2] = vnum1 - crvlen;
	v[0] = vnum1 - 1;
	v[1] = v[0] - crvlen;
	break;
    }

    v1 = vlist + v[0];
    v2 = vlist + v[1];
    v3 = vlist + v[2];

    if (VERTEX_IS_UNDEFINED(*v1)
	|| VERTEX_IS_UNDEFINED(*v2)
	|| VERTEX_IS_UNDEFINED(*v3)
	)
	return (-2);

    /* All else OK, fill in the polygon: */

    p = (struct polygon *)nextfrom_dynarray(&polygons);

    memcpy (p->vertex, v, sizeof(v));
#if ! HIDDEN3D_QUADTREE
    p->next = -1;
#endif

    /* Some helper macros for repeted code blocks: */

    /* Gets Minimum 'var' value of polygon 'poly' into variable
     * 'min. C is one of x, y, or z: */
#define GET_MIN(poly, var, min)			\
    do {					\
	int i;					\
	long *v = poly->vertex;			\
						\
	min = vlist[*v++].var;			\
	for (i = 1; i< POLY_NVERT; i++, v++)	\
	    if (vlist[*v].var < min)		\
		min = vlist[*v].var;		\
        assert(min >= -surface_scale);		\
    } while (0)

    /* Gets Maximum 'var' value of polygon 'poly', as with GET_MIN */
#define GET_MAX(poly, var, max)			\
    do {					\
	int i;					\
	long *v = poly->vertex;			\
						\
	max = vlist[*v++].var;			\
	for (i = 1; i< POLY_NVERT; i++, v++)	\
	    if (vlist[*v].var > max)		\
		max = vlist[*v].var;		\
	assert(max <= surface_scale);		\
    } while (0)

    GET_MIN(p, x, p->xmin);
    GET_MIN(p, y, p->ymin);
    GET_MIN(p, z, p->zmin);
    GET_MAX(p, x, p->xmax);
    GET_MAX(p, y, p->ymax);
    GET_MAX(p, z, p->zmax);
#undef GET_MIN
#undef GET_MAX

#if HIDDEN3D_GRIDBOX
    p->xbits = CALC_BITRANGE(p->xmin, p->xmax);
    p->ybits = CALC_BITRANGE(p->ymin, p->ymax);
#endif

    p->frontfacing = get_plane(p, p->plane);

    return (p - plist);
}

/* color edges, based on the orientation of polygon(s). One of the two
 * edges passed in is a new one, meaning there is no other polygon
 * sharing it, yet. The other, 'old' edge is common to the new polygon
 * and another one, which was created earlier on. If these two polygon
 * differ in their orientation (one front-, the other backsided to the
 * viewer), this routine has to resolve that conflict.  Edge colours
 * are changed only if the edge wasn't invisible, before */
static void
color_edges(
    long int new_edge,		/* index of 'new', conflictless edge */
    long int old_edge,		/* index of 'old' edge, may conflict */
    long int new_poly,		/* index of current polygon */
    long int old_poly,		/* index of poly sharing old_edge */
    int above,			/* style number for front of polygons */
    int below)			/* style number for backside of polys */
{
    int casenumber;

    if (new_poly > -2) {
	/* new polygon was built successfully */
	if (old_poly <= -2)
	    /* old polygon doesn't exist. Use new_polygon for both: */
	    old_poly = new_poly;

	casenumber =
	    (plist[new_poly].frontfacing ? 1 : 0)
	    + 2 * (plist[old_poly].frontfacing ? 1 : 0);
	switch (casenumber) {
	case 0:
	    /* both backfacing */
	    if (elist[new_edge].style >= -2)
		elist[new_edge].style	= below;
	    if (elist[old_edge].style >= -2)
		elist[old_edge].style = below;
	    break;
	case 2:
	    if (elist[new_edge].style >= -2)
		elist[new_edge].style = below;
	    /* FALLTHROUGH */
	case 1:
	    /* new front-, old one backfacing, or */
	    /* new back-, old one frontfacing */
	    if (((new_edge == old_edge)
		 && hiddenHandleBentoverQuadrangles) /* a diagonal edge! */
		|| (elist[old_edge].style >= -2)) {
		/* conflict has occured: two polygons meet here, with opposige
		 * sides being shown. What's to do?
		 * 1) find a vertex of one polygon outside this common
		 * edge
		 * 2) check wether it's in front of or behind the
		 * other polygon's plane
		 * 3) if in front, color the edge accoring to the
		 * vertex' polygon, otherwise, color like the other
		 * polygon */
		long int vnum1 = elist[old_edge].v1;
		long int vnum2 = elist[old_edge].v2;
		p_polygon p = plist + new_poly;
		long int pvert = -1;
		double point_to_plane;

		if (p->vertex[0] == vnum1) {
		    if (p->vertex[1] == vnum2) {
			pvert = p->vertex[2];
		    } else if (p->vertex[2] == vnum2) {
			pvert = p->vertex[1];
		    }
		} else if (p->vertex[1] == vnum1) {
		    if (p->vertex[0] == vnum2) {
			pvert = p->vertex[2];
		    } else if (p->vertex[2] == vnum2) {
			pvert = p->vertex[0];
		    }
		} else if (p->vertex[2] == vnum1) {
		    if (p->vertex[0] == vnum2) {
			pvert = p->vertex[1];
		    } else if (p->vertex[1] == vnum2) {
			pvert = p->vertex[0];
		    }
		}
		assert (pvert >= 0);

		point_to_plane =
		    eval_plane_equation(plist[old_poly].plane, vlist + pvert);

		if (point_to_plane > 0) {
		    /* point in new_poly is in front of old_poly plane */
		    elist[old_edge].style = p->frontfacing ? above : below;
		}	else {
		    elist[old_edge].style =
			plist[old_poly].frontfacing ? above : below;
		}
	    }
	    break;
	case 3:
	    /* both frontfacing: nothing to do */
	    break;
	} /* switch */
    } else {
	/* Ooops? build_networks() must have guessed incorrectly that
	 * this polygon should exist. */
	return;
    }
}


/* This somewhat monstrous routine fills the vlist, elist and plist
 * dynamic arrays with values from all those plots. It strives to
 * respect all the topological linkage between vertices, edges and
 * polygons. E.g., it has to find the correct color for each edge,
 * based on the orientation of the two polygons sharing it, WRT both
 * the observer and each other. */
/* NEW FEATURE HBB 20000715: allow non-grid datasets too, by storing
 * only vertices and 'direct' edges, but no polygons or 'cross' edges
 * */
static void
build_networks(struct surface_points *plots, int pcount)
{
    long int i;
    struct surface_points *this_plot;
    int surface;		/* count the surfaces (i.e. sub-plots) */
    long int crv, ncrvs;	/* count isolines */
    long int nverts;		/* count vertices */
    long int max_crvlen;	/* maximal length of isoline in any plot */
    long int nv, ne, np;	/* local poly/edge/vertex counts */
    long int *north_polygons;	/* stores polygons of isoline above */
    long int *these_polygons;	/* same, being built for use by next turn */
    long int *north_edges;	/* stores edges of polyline above */
    long int *these_edges;	/* same, being built for use by next turn */
    struct iso_curve *icrvs;
    int above = -3, below;	/* line type for edges of front/back side*/
    struct lp_style_type *lp;	/* pointer to line and point properties */

    /* Count out the initial sizes needed for the polygon and vertex
     * lists. */
    nv = ne = np = 0;
    max_crvlen = -1;

    for (this_plot = plots, surface = 0;
	surface < pcount;
	this_plot = this_plot->next_sp, surface++) {
	long int crvlen;
	
	/* Quietly skip empty plots */
	if (this_plot->plot_type == NODATA)
	    continue;

	crvlen = this_plot->iso_crvs->p_count;

	/* Allow individual plots to opt out of hidden3d calculations */
	if (this_plot->opt_out_of_hidden3d)
	    continue;

	/* register maximal isocurve length. Only necessary for
	 * grid-topology plots that will create polygons, so I can do
	 * it here, already. */
	if (crvlen > max_crvlen)
	    max_crvlen = crvlen;

	/* count 'curves' (i.e. isolines) and vertices in this plot */
	nverts = 0;
	if(this_plot->plot_type == FUNC3D) {
	    ncrvs = 0;
	    for(icrvs = this_plot->iso_crvs;
		icrvs; icrvs = icrvs->next) {
		ncrvs++;
	    }
	    nverts += ncrvs * crvlen;
	} else if(this_plot->plot_type == DATA3D) {
	    ncrvs = this_plot->num_iso_read;
	    if (this_plot->has_grid_topology)
		nverts += ncrvs * crvlen;
	    else if (this_plot->plot_style == VECTOR)
		nverts += this_plot->iso_crvs->p_count;
	    else {
		/* have to check each isoline separately: */
		for (icrvs = this_plot->iso_crvs; icrvs; icrvs = icrvs->next)
		    nverts += icrvs->p_count;
	    }
	} else {
	    graph_error("Plot type is neither function nor data");
	    return;
	}

	/* To avoid possibly suprising error messages, several 2d-only
	 * plot styles are mapped to others, that are genuinely
	 * available in 3d. */
	switch (this_plot->plot_style) {
	case LINESPOINTS:
	case STEPS:
	case FSTEPS:
	case HISTEPS:
	case LINES:
	    nv += nverts;
	    ne += nverts - ncrvs;
	    if (this_plot->has_grid_topology) {
		ne += 2 * nverts - ncrvs - 2 * crvlen + 1;
		np += 2 * (ncrvs - 1) * (crvlen - 1);
	    }
	    break;
	case BOXES:
	case FILLEDCURVES:
	case IMPULSES:
	case VECTOR:
	    nv += 2 * nverts;
	    ne += nverts;
	    break;
	case POINTSTYLE:
	default:
	    /* treat all remaining ones like 'points' */
	    nv += nverts;
	    ne += nverts; /* a 'phantom edge' per isolated point */
	    break;
	} /* switch */
    } /* for (plots) */

    /* Check for no data at all */
    if (max_crvlen <= 0)
	return;

    /* allocate all the lists to the size we need: */
    resize_dynarray(&vertices, nv);
    resize_dynarray(&edges, ne);
    resize_dynarray(&polygons, np);

    /* allocate the storage for polygons and edges of the isoline just
     * above the current one, to allow easy access to them from the
     * current isoline */
    north_polygons = (long*)gp_alloc(2 * max_crvlen * sizeof (long int),
			      "hidden north_polys");
    these_polygons = (long*)gp_alloc(2 * max_crvlen * sizeof (long int),
			      "hidden these_polys");
    north_edges = (long*)gp_alloc(3 * max_crvlen * sizeof (long int),
			   "hidden north_edges");
    these_edges = (long*)gp_alloc(3 * max_crvlen * sizeof (long int),
			   "hidden these_edges");

    /* initialize the lists, all in one large loop. This is different
     * from the previous approach, which went over the vertices,
     * first, and only then, in new loop, built polygons */
    for (this_plot = plots, surface = 0;
	 surface < pcount;
	 this_plot = this_plot->next_sp, surface++) {
	lp_style_type *lp_style = NULL;
	TBOOLEAN color_from_column = this_plot->pm3d_color_from_column;
	long int crvlen;

	/* Quietly skip empty plots */
	if (this_plot->plot_type == NODATA)
	    continue;

	crvlen = this_plot->iso_crvs->p_count;

	/* Allow individual plots to opt out of hidden3d calculations */
	if (this_plot->opt_out_of_hidden3d)
	    continue;

	lp = &(this_plot->lp_properties);
	above = this_plot->lp_properties.l_type;
	below = this_plot->lp_properties.l_type + hiddenBacksideLinetypeOffset;

	/* calculate the point symbol type: */
	/* Assumest hat upstream functions have made sure this is
	 * initialized sensibly --- thou hast been warned */
	lp_style = &(this_plot->lp_properties);

	/* HBB 20000715: new initialization code block for non-grid
	 * structured datasets. Sufficiently different from the rest
	 * to warrant separate code, I think. */
	if (! this_plot->has_grid_topology) {
	    for (crv = 0, icrvs = this_plot->iso_crvs;
		 icrvs;
		 crv++, icrvs = icrvs->next) {
		struct coordinate GPHUGE *points = icrvs->points;
		long int previousvertex = -1;

#ifdef EAM_DATASTRINGS
		/* To handle labels we must look inside a separate list */
		/* rather than just walking through the points arrays.  */
		if (this_plot->plot_style == LABELPOINTS) {
		    struct text_label *label;
		    long int thisvertex;
		    struct coordinate labelpoint;
		    memset( &labelpoint , 0 , sizeof( struct coordinate ));

		    lp->pointflag = 1; /* Labels can use the code for hidden points */
		    for (label = this_plot->labels->next; label != NULL; label = label->next) {
			labelpoint.x = label->place.x;
			labelpoint.y = label->place.y;
			labelpoint.z = label->place.z;
			thisvertex = store_vertex(&labelpoint, 
				&(this_plot->lp_properties), color_from_column);
			if (thisvertex < 0)
			    continue;
			(vlist+thisvertex)->label = label;
			store_edge(thisvertex, edir_point, crvlen, lp, above);
		    }
		} else
#endif

		for (i = 0; i < icrvs->p_count; i++) {
		    long int thisvertex, basevertex;

		    thisvertex = store_vertex(points + i, lp_style,
					      color_from_column);

		    if (this_plot->plot_style == VECTOR) {
			store_vertex(icrvs->next->points+i, 0, 0);
		    }

		    if (thisvertex < 0) {
			previousvertex = thisvertex;
			continue;
		    }

		    switch (this_plot->plot_style) {
		    case LINESPOINTS:
		    case STEPS:
		    case FSTEPS:
		    case HISTEPS:
		    case LINES:
			if (previousvertex >= 0)
			    store_edge(thisvertex, edir_west, 0, lp, above);
			break;
		    case VECTOR:
			store_edge(thisvertex, edir_vector, 0, lp, above);
			break;
		    case BOXES:
		    case FILLEDCURVES:
		    case IMPULSES:
			/* set second vertex to the low end of zrange */
			{
			    coordval remember_z = points[i].z;

			    points[i].z = axis_array[FIRST_Z_AXIS].min;
			    basevertex = store_vertex(points + i, lp_style,
						      color_from_column);
			    points[i].z = remember_z;
			}
			if (basevertex > 0)
			    store_edge(basevertex, edir_impulse, 0, lp, above);
			break;

		    case POINTSTYLE:
		    default:	/* treat all the others like 'points' */
			store_edge(thisvertex, edir_point, crvlen, lp, above);
			break;
		    } /* switch(plot_style) */

		    previousvertex = thisvertex;
		} /* for(vertex) */
	    } /* for(crv) */

	    continue;		/* done with this plot! */
	}

	/* initialize stored indices of north-of-this-isoline polygons and
	 * edges properly */
	for (i=0; i < this_plot->iso_crvs->p_count; i++) {
	    north_polygons[2 * i]
		= north_polygons[2 * i + 1]
		= north_edges[3 * i]
		= north_edges[3 * i + 1]
		= north_edges[3 * i + 2]
		= -3;
	}

	for (crv = 0, icrvs = this_plot->iso_crvs;
	     icrvs;
	     crv++, icrvs = icrvs->next) {
	    struct coordinate GPHUGE *points = icrvs->points;

	    for (i = 0; i < icrvs->p_count; i++) {
		long int thisvertex, basevertex;
		long int e1, e2, e3;
		long int pnum;

		thisvertex = store_vertex(points + i, lp_style,
					  color_from_column);

		/* Preset the pointers to the polygons and edges
		 * belonging to this isoline */
		these_polygons[2 * i] = these_polygons[2 * i + 1]
		    = these_edges[3 * i] = these_edges[3 * i + 1]
		    = these_edges[3 * i + 2]
		    = -3;

		switch (this_plot->plot_style) {
		case LINESPOINTS:
		case STEPS:
		case FSTEPS:
		case HISTEPS:
		case LINES:
		    if (i > 0) {
			/* not first point, so we might want to set up
			 * the edge(s) to the left of this vertex */
			if (thisvertex < 0) {
			    if ((crv > 0)
				&& (hiddenShowAlternativeDiagonal)
				) {
				/* this vertex is invalid, but the
				 * other three might still form a
				 * valid triangle, facing northwest to
				 * do that, we'll need the 'wrong'
				 * diagonal, which goes from SW to NE:
				 * */
				these_edges[i*3+2] = e3
				    = store_edge(vertices.end - 1, edir_NE, crvlen,
						 lp, above);
				if (e3 > -2) {
				    /* don't store this polygon for
				     * later: it doesn't share edges
				     * with any others to the south or
				     * east, so there's need to */
				    pnum
					= store_polygon(vertices.end - 1, pdir_NW, crvlen);
				    /* The other two edges of this
				     * polygon need to be checked
				     * against the neighboring
				     * polygons' orientations, before
				     * being coloured */
				    color_edges(e3, these_edges[3*(i-1) +1],
						pnum, these_polygons[2*(i-1) + 1],
						above, below);
				    color_edges(e3, north_edges[3*i],
						pnum, north_polygons[2*i], above, below);
				}
			    }
			    break; /* nothing else to do for invalid vertex */
			}

			/* Coming here means that the current vertex
			 * is valid: check the other three of this
			 * cell, by trying to set up the edges from
			 * this one to there */
			these_edges[i*3] = e1
			    = store_edge(thisvertex, edir_west, crvlen, lp, above);

			if (crv > 0) { /* vertices to the north exist */
			    these_edges[i*3 + 1] = e2
				= store_edge(thisvertex, edir_north, crvlen, lp, above);
			    these_edges[i*3 + 2] = e3
				= store_edge(thisvertex, edir_NW, crvlen, lp, above);
			    if (e3 > -2) {
				/* diagonal edge of this cell is OK,
				 * so try to build both the polygons:
				 * */
				if (e1 > -2) {
				    /* one pair of edges is valid: put
				     * first polygon, which points
				     * towards the southwest */
				    these_polygons[2*i] = pnum
					= store_polygon(thisvertex, pdir_SW, crvlen);
				    color_edges(e1, these_edges[3*(i-1)+1],
						pnum, these_polygons[2*(i-1)+ 1], above, below );
				}
				if (e2 > -2) {
				    /* other pair of two is fine, put
				     * the northeast polygon: */
				    these_polygons[2*i + 1] = pnum
					= store_polygon(thisvertex, pdir_NE, crvlen);
				    color_edges(e2, north_edges[3*i],
						pnum, north_polygons[2*i], above, below);
				}
				/* In case these two new polygons
				 * differ in orientation, find good
				 * coloring of the diagonal */
				color_edges(e3, e3, these_polygons[2*i],
					    these_polygons[2*i+1], above, below);
			    } /* if e3 valid */
			    else if ((e1 > -2) && (e2 > -2)
				     && hiddenShowAlternativeDiagonal) {
				/* looks like all but the north-west
				 * vertex are usable, so we set up the
				 * southeast-pointing triangle, using
				 * the 'wrong' diagonal: */
				these_edges[3*i + 2] = e3
				    = store_edge(thisvertex, edir_NE, crvlen, lp, above);
				if (e3 > -2) {
				    /* fill this polygon into *both*
				     * polygon places for this
				     * quadrangle, as this triangle
				     * coincides with both edges that
				     * will be used by later polygons
				     * */
				    these_polygons[2*i] = these_polygons[2*i+1] = pnum
					= store_polygon(thisvertex, pdir_SE, crvlen);
				    /* This case is somewhat special:
				     * all edges are new, so there is
				     * no other polygon orientation to
				     * consider */
				    if (!plist[pnum].frontfacing)
					elist[e1].style = elist[e2].style = elist[e3].style
					    = below;
				}
			    }
			}
		    } else if ((crv > 0)
			       && (thisvertex >= 0)) {
			/* We're at the west border of the grid, but
			 * not on the north one: put vertical end-wall
			 * edge:*/
			these_edges[3*i + 1] =
			    store_edge(thisvertex, edir_north, crvlen, lp, above);
		    }
		    break;

		case BOXES:
		case FILLEDCURVES:
		case IMPULSES:
		    if (thisvertex < 0)
			break;

		    /* set second vertex to the low end of zrange */
		    {
			coordval remember_z = points[i].z;

			points[i].z = axis_array[FIRST_Z_AXIS].min;
			basevertex = store_vertex(points + i, lp_style,
						  color_from_column);
			points[i].z = remember_z;
		    }
		    if (basevertex > 0)
			store_edge(basevertex, edir_impulse, 0, lp, above);
		    break;

		case POINTSTYLE:
		default:	/* treat all the others like 'points' */
		    if (thisvertex < 0) /* Ignore invalid vertex */
			break;
		    store_edge(thisvertex, edir_point, crvlen, lp, above);
		    break;
		} /* switch */
	    } /* for(i) */

	    /* Swap the 'north' lists of polygons and edges with
	     * 'these' ones, which have been filled in the pass
	     * through this isocurve */
	    {
		long int *temp = north_polygons;
		north_polygons = these_polygons;
		these_polygons = temp;

		temp = north_edges;
		north_edges = these_edges;
		these_edges = temp;
	    }
	} /* for(isocrv) */
    } /* for(plot) */

    free (these_polygons);
    free (north_polygons);
    free (these_edges);
    free (north_edges);
}

/* Sort the elist in order of growing zmax. Uses qsort on an array of
 * plist indices, and then fills in the 'next' fields in struct
 * polygon to store the resulting order inside the plist */
/* HBB 20010720: removed 'static' to avoid HP-sUX gcc bug */
int
compare_edges_by_zmin(SORTFUNC_ARGS p1, SORTFUNC_ARGS p2)
{
    return SIGN(vlist[elist[*(const long *) p1].v2].z
		- vlist[elist[*(const long *) p2].v2].z);
}

static void sort_edges_by_z()
{
	long *sortarray, i;
	p_edge this_;
	
	if (!edges.end)
		return;
	
	sortarray = (long *)gp_alloc((unsigned long) sizeof(long) * edges.end,
		"hidden sort edges");
	/* initialize sortarray with an identity mapping */
	for (i = 0; i < edges.end; i++)
		sortarray[i] = i;
	/* sort it */
	qsort(sortarray, (size_t) edges.end, sizeof(long), compare_edges_by_zmin);
	
	/* traverse plist in the order given by sortarray, and set the
	* 'next' pointers */
	this_ = elist + sortarray[0];
	for (i = 1; i < edges.end; i++) 
	{
		this_->next = sortarray[i];
		this_ = elist + sortarray[i];
	}
	this_->next = -1L;
	
	/* 'efirst' is the index of the leading element of plist */
	efirst = sortarray[0];
	
	free(sortarray);
}

/* HBB 20010720: removed 'static' to avoid HP-sUX gcc bug */
int
compare_polys_by_zmax(SORTFUNC_ARGS p1, SORTFUNC_ARGS p2)
{
    return (SIGN(plist[*(const long *) p1].zmax
		 - plist[*(const long *) p2].zmax));
}

static void sort_polys_by_z()
{
	long *sortarray, i;
	p_polygon this_;
	
	if (!polygons.end)
		return;
	
	sortarray = (long *)gp_alloc((unsigned long) sizeof(long) * polygons.end,
		"hidden sortarray");
	
	/* initialize sortarray with an identity mapping */
	for (i = 0; i < polygons.end; i++)
		sortarray[i] = i;
	
	/* sort it */
	qsort(sortarray, (size_t) polygons.end, sizeof(long),
		compare_polys_by_zmax);
	
		/* traverse plist in the order given by sortarray, and set the
	* 'next' pointers */
#if HIDDEN3D_QUADTREE
	/* HBB 20000716: Loop backwards, to ease construction of
	* linked lists from the head: */
	{
		unsigned int grid_x, grid_y;
		unsigned int grid_x_low, grid_x_high, grid_y_low, grid_y_high;
		
		for (grid_x = 0; grid_x < QUADTREE_GRANULARITY; grid_x++)
			for (grid_y = 0; grid_y < QUADTREE_GRANULARITY; grid_y++)
				quadtree[grid_x][grid_y] = -1;
			
			for (i=polygons.end - 1; i >= 0; i--) {
				this = plist + sortarray[i];
				
				grid_x_low = COORD_TO_TREECELL(this_->xmin);
				grid_x_high = COORD_TO_TREECELL(this_->xmax);
				grid_y_low = COORD_TO_TREECELL(this_->ymin);
				grid_y_high = COORD_TO_TREECELL(this_->ymax);
				
				for (grid_x = grid_x_low; grid_x <= grid_x_high; grid_x++) {
					for (grid_y = grid_y_low; grid_y <= grid_y_high; grid_y++) {
						p_qtreelist newhead = nextfrom_dynarray(&qtree);
						
						newhead->next = quadtree[grid_x][grid_y];
						newhead->p = sortarray[i];
						
						quadtree[grid_x][grid_y] = newhead - qlist;
					}
				}
			}
	}
	
#else /* HIDDEN3D_QUADTREE */
	this_ = plist + sortarray[0];
	for (i = 1; i < polygons.end; i++) 
	{
		this_->next = sortarray[i];
		this_ = plist + sortarray[i];
	}
	this_->next = -1L;
	/* 'pfirst' is the index of the leading element of plist */
#endif /* HIDDEN3D_QUADTREE */
	pfirst = sortarray[0];
	
	free(sortarray);
}


#define LASTBITS (USHRT_BITS -1)	/* ????? */

/************************************************/
/*******            Drawing the polygons ********/
/************************************************/

/* draw a single vertex as a point symbol, if requested by the chosen
 * plot style (linespoints, points, or dots...) */
static void
draw_vertex(p_vertex v)
{
    unsigned int x, y;

    TERMCOORD(v, x, y);
    if (v->lp_style && v->lp_style->p_type >= 0 && !clip_point(x,y)) {
	int colortype = v->lp_style->pm3d_color.type;

#ifdef EAM_DATASTRINGS
	if (v->label)  {
	    write_label(x,y, v->label);
	    v->lp_style = NULL;
	    return;
	}
#endif

	/* EAM DEBUG - Check for extra point properties */
	if (colortype == TC_LT)
	    /* Should have been set already! */
	    ;
	else if (colortype == TC_RGB && v->lp_style->pm3d_color.lt == LT_COLORFROMCOLUMN)
	    set_rgbcolor(v->real_z);
	else if (colortype == TC_RGB)
	    set_rgbcolor(v->lp_style->pm3d_color.lt);
	else if (colortype == TC_Z)
	    set_color( cb2gray(z2cb(v->real_z)) );

#ifdef HIDDEN3D_VAR_PTSIZE
	if (v->lp_style->p_size == PTSZ_VARIABLE)
	    (term->pointsize)(pointsize * v->original->CRD_PTSIZE);
#endif

	(term->point)(x,y, v->lp_style->p_type);

	/* vertex has been drawn --> flag it as done */
	v->lp_style = NULL;
    }
}


/* The function that actually does the drawing of the visible portions
 * of lines */
/* HBB 20001108: changed to take the pointers to the end vertices as
 * additional arguments. */
static void
draw_edge(p_edge e, p_vertex v1, p_vertex v2)
{
    assert (e >= elist && e < elist + edges.end);

    draw3d_line_unconditional(v1, v2, e->lp, e->style);
    if (e->lp->pointflag) {
	draw_vertex(v1);
	draw_vertex(v2);
    }
}


/*************************************************************/
/*************************************************************/
/*******   The depth sort algorithm (in_front) and its  ******/
/*******   whole lot of helper functions                ******/
/*************************************************************/
/*************************************************************/

/* Split a given line segment into two at an inner point. The inner
 * point is specified as a fraction of the line-length (0 is V1, 1 is
 * V2) */
/* HBB 20001108: changed to now take two vertex pointers as its
 * arguments, rather than an edge pointer. */
/* HBB 20001204: changed interface again. Now use vertex indices,
 * rather than pointers, to avoid problems with dangling pointers
 * after nextfrom_dynarray() call. */
static long split_line_at_ratio(
				long vnum1, long vnum2, 	/* vertex indices of line to split */
				double w)			/* where to split it */
{
	p_vertex v;
	
	if (EQ(w, 0.0))
		return vnum1;
	if (EQ(w, 1.0))
		return vnum2;
	
	/* Create a new vertex */
	v = (struct vertex *)nextfrom_dynarray(&vertices);
	
	v->x = (vlist[vnum2].x - vlist[vnum1].x) * w + vlist[vnum1].x;
	v->y = (vlist[vnum2].y - vlist[vnum1].y) * w + vlist[vnum1].y;
	v->z = (vlist[vnum2].z - vlist[vnum1].z) * w + vlist[vnum1].z;
	v->real_z = (vlist[vnum2].real_z - vlist[vnum1].real_z) * w
		+ vlist[vnum1].real_z;
	
	/* no point symbol for vertices generated by splitting an edge */
	v->lp_style = NULL;
	
	/* additional checks to prevent adding unnecessary vertices */
	if (V_EQUAL(v, vlist + vnum1)) 
	{
		droplast_dynarray(&vertices);
		return vnum1;
	}
	if (V_EQUAL(v, vlist + vnum2)) {
		droplast_dynarray(&vertices);
		return vnum2;
	}
	
	return (v - vlist);
}


/* Compute the 'signed area' of 3 points in their 2d projection
 * to the x-y plane. Essentially the z component of the crossproduct.
 * Should come out positive if v1, v2, v3 are ordered counter-clockwise */

static GP_INLINE double
area2D(p_vertex v1, p_vertex v2, p_vertex v3)
{
    double
	dx12 = v2->x - v1->x,	/* x/y components of (v2-v1) and (v3-v1) */
	dx13 = v3->x - v1->x,
	dy12 = v2->y - v1->y,
	dy13 = v3->y - v1->y;
    return (dx12 * dy13 - dy12 * dx13);
}

/* Utility routine: takes an edge and makes a new one, which is a fragment
 * of the old one. The fragment inherits the line style and stuff of the
 * given edge; only the two new vertices are different. The new edge
 * is then passed to in_front, for recursive handling */
/* HBB 20001108: Changed from edge pointer to edge index. Don't
 * allocate a fresh anymore, as this is no longer needed after the
 * change to in_front().  What remains of this function may no longer
 * be worth having. I.e. it can be replaced by a direct recursion call
 * of in_front(), sometime soon. */
static GP_INLINE void handle_edge_fragment(long edgenum, long vnum1, long vnum2, long firstpoly)
{
#if !HIDDEN3D_QUADTREE
    /* Avoid checking against the same polygon again. */
    firstpoly = plist[firstpoly].next;
#endif
    in_front(edgenum, vnum1, vnum2, &firstpoly);
}

/*********************************************************************/
/* The actual heart of all this: determines if edge at index 'edgenum'
 * of the elist is in_front of all the polygons, or not. If necessary,
 * it will recursively call itself to isolate more than one visible
 * fragment of the input edge. Wherever possible, recursion is
 * avoided, by in-place modification of the edge.
 *
 * The visible fragments are then drawn by a call to 'draw_edge' from
 * inside this routine. */
/*********************************************************************/
/* HBB 20001108: changed to now take the vertex numbers as additional
 * arguments. The idea is to not overwrite the endpoint stored with
 * the edge, so Test 2 will catch on even after the subject edge has
 * been split up before one of its two polygons is tested against
 * it. FIXME: allocates new vertices when splitting, but never frees
 * them, currently. */

static int in_front(
    long edgenum,		/* number of the edge in elist */
    long vnum1, long vnum2,	/* numbers of its endpoints */
    long *firstpoly)		/* first plist index to consider */
{
    p_polygon p;		/* pointer to current testing polygon */
    long int polynum;		/* ... and its index in the plist */
    p_vertex v1, v2;		/* pointers to vertices of input edge */

    coordval xmin, xmax;	/* all of these are for the edge */
    coordval ymin, ymax;
    coordval zmin, zmax;
#if HIDDEN3D_GRIDBOX
    unsigned int xextent;	/* extent bitmask in x direction */
    unsigned int yextent;	/* same, in y direction */

# define SET_XEXTENT \
  xextent = CALC_BITRANGE(xmin, xmax);
# define SET_YEXTENT \
  yextent = CALC_BITRANGE(ymin, ymax);
#else
# define SET_XEXTENT /* nothing */
# define SET_YEXTENT /* nothing */
#endif
#if HIDDEN3D_QUADTREE
    unsigned int grid_x, grid_y;
    unsigned int grid_x_low, grid_x_high;
    unsigned int grid_y_low, grid_y_high;
    long listhead;
#endif

    /* zmin of the edge, as it started out. This is needed separately to
     * allow modifying '*firstpoly', without moving it too far to the
     * front. */
    coordval first_zmin;

    /* macro for eliminating tail-recursion inside in_front: when the
     * current edge is modified, recompute all function-wide status
     * variables. Note that it guarantees that v1 is always closer to
     * the viewer than v2 (in z direction) */
    /* HBB 20001108: slightly changed so it can be called with vnum1
     * and vnum2 as its arguments, too */
#define setup_edge(vert1, vert2)		\
    do {					\
	if (vlist[vert1].z > vlist[vert2].z) {	\
	    v1 = vlist + (vert1);		\
	    v2 = vlist + (vert2);		\
	} else {				\
	    v1 = vlist + (vert2);		\
	    v2 = vlist + (vert1);		\
	}					\
	vnum1 = v1 - vlist;			\
	vnum2 = v2 - vlist;			\
	zmax = v1->z;	zmin = v2->z;		\
						\
	if (v1->x > v2->x) {			\
	    xmin = v2->x;	xmax = v1->x;	\
	} else {				\
	    xmin = v1->x;	xmax = v2->x;	\
	}					\
	SET_XEXTENT;				\
						\
	if (v1->y > v2->y) {			\
	    ymin = v2->y;	ymax = v1->y;	\
	} else {				\
	    ymin = v1->y;	ymax = v2->y;	\
	}					\
	SET_YEXTENT;				\
    } while (0) /* end macro setup_edge */

    /* use the macro for initial setup, too: */
    setup_edge(vnum1, vnum2);

    first_zmin = zmin;

#if HIDDEN3D_QUADTREE
    grid_x_low = COORD_TO_TREECELL(xmin);
    grid_x_high = COORD_TO_TREECELL(xmax);
    grid_y_low = COORD_TO_TREECELL(ymin);
    grid_y_high = COORD_TO_TREECELL(ymax);

    for (grid_x = grid_x_low; grid_x <= grid_x_high; grid_x ++)
	for (grid_y = grid_y_low; grid_y <= grid_y_high; grid_y ++)
	    for (listhead = quadtree[grid_x][grid_y];
		 listhead >= 0;
		 listhead = qlist[listhead].next)
#else /* HIDDEN3D_QUADTREE */
    /* loop over all the polygons in the sorted list, starting at the
     * currently first (i.e. furthest, from the viewer) polgon. */
    for (polynum = *firstpoly; polynum >=0; polynum = p->next)
#endif /* HIDDEN3D_QUADTREE */
	{
	    /* shortcut variables for the three vertices of 'p':*/
	    p_vertex w1, w2, w3;
	    /* signed areas of each of e's vertices wrt. the edges of
	     * p. I store them explicitly, because they are used more
	     * than once, eventually */
	    double e_side[2][POLY_NVERT];
	    /* signed areas of each of p's vertices wrt. to edge e. Stored
	     * for re-use in intersection calculations, as well */
	    double p_side[POLY_NVERT];
	    /* values got from throwing the edge vertices into the plane
	     * equation of the polygon. Used for testing if the vertices are
	     * in front of or behind the plane, and also to determine the
	     * point where the edge goes through the polygon, by
	     * interpolation. */
	    double v1_rel_pplane, v2_rel_pplane;
	    /* Orientation of polygon wrt. to the eye: front or back side
	     * visible? */
	    coordval polyfacing;
	    int whichside;
	    /* stores classification of cases as 4 2-bit patterns, mainly */
	    unsigned int classification[POLY_NVERT + 1];

#if HIDDEN3D_QUADTREE
	    polynum = qlist[listhead].p;
#endif
	    p = plist + polynum;

	    /* OK, off we go with the real work. This algo is mainly
	     * following the one of 'HLines.java', as described in the
	     * book 'Computer Graphics for Java Programmers', by Dutch
	     * professor Leen Ammeraal, published by J.Wiley & Sons,
	     * ISBN 0 471 98142 7. It happens to the only(!) actual
	     * code or algorithm example I got out of a question to
	     * the Newsgroup comp.graphics.algorithms, asking for a
	     * hidden *line* removal algorithm. */

	    /* Test 1 (2D): minimax tests. Do x/y ranges of polygon
	     * and edge have any overlap? */
	    if (0
#if HIDDEN3D_GRIDBOX
		/* First, check by comparing the extent bit patterns: */
		|| (!(xextent & p->xbits))
		|| (!(yextent & p->ybits))
#endif
		|| (p->xmax < xmin)
		|| (p->xmin > xmax)
		|| (p->ymax < ymin)
		|| (p->ymin > ymax)
		)
		continue;

	    /* Tests 2 and 3 switched... */

	    /* Test 3 (3D): Is edge completely in front of polygon? */
	    if (p->zmax < zmin) {
		/* Polygon completely behind this edge. Move start of
		 * relevant plist to this point, to speed up next
		 * run. This makes use of the fact that elist is also
		 * kept in upwardly sorted order of zmin, i.e. the
		 * condition found here will also hold for all coming
		 * edges in the list */
		if (p->zmax < first_zmin)
		    *firstpoly = polynum;
		continue;	/* this polygon is done with */
	    }

	    /* Test 2 (0D): does edge belong to this very polygon? */
	    /* 20001108: to make this rejector more effective, do keep
	     * the original edge vertices unchanged */
	    if (1
		&& (0
		    || (p->vertex[0] == elist[edgenum].v1)
		    || (p->vertex[1] == elist[edgenum].v1)
		    || (p->vertex[2] == elist[edgenum].v1)
		    )
		&& (0
		    || (p->vertex[0] == elist[edgenum].v2)
		    || (p->vertex[1] == elist[edgenum].v2)
		    || (p->vertex[2] == elist[edgenum].v2)
		    )
		)
		continue;

	    w1 = vlist + p->vertex[0];
	    w2 = vlist + p->vertex[1];
	    w3 = vlist + p->vertex[2];


	    /* Test 4 (2D): Is edge fully on the 'wrong side' of one of the
	     * polygon edges?
	     *
	     * Done by comparing signed areas (cross-product z components of
	     * both edge endpoints) to that of the polygon itself (as given by
	     * the plane normal z component). If both v1 and v2 are found on
	     * the 'outside' side of a polygon edge (interpreting it as
	     * infinite line), the polygon can't obscure this edge. All the
	     * signed areas are remembered, for later tests. */
	    polyfacing = p->plane[2];
	    whichside = SIGN(polyfacing);
	    if (!whichside)
		/* We're looking at the border of this polygon. It
		 * looks like an infitisemally thin line, i.e. it'll
		 * be practically invisible*/
		/* FIXME: should such a polygon have been stored in the plist,
		 * in the first place??? */
		continue;

	    /* p->plane[2] is now forced >=0, so it doesn't tell anything
	     * about which side the polygon is facing */
	    whichside = (p->frontfacing) ? -1 : 1;

#if 0
	    /* Just make sure I got this the right way round: */
	    if (p->frontfacing)
		assert (GE(area2D(w1, w2, w3),0));
	    else
		assert (GE(0,area2D(w1, w2, w3)));
#endif

	    /* classify both endpoints of the test edge, by their
	     * position relative to some plane that defines the volume
	     * shadowed by a triangle. Used both by macro 'checkside',
	     * and for individual checking, below --> put this into a
	     * macro */
#define classify(par1, par2, classnum)					\
	    do {							\
		if (GR(par1, 0)) {					\
		    /* v1 strictly outside --> v2 in-plane is enough,	\
                       already : */					\
		    classification[classnum]				\
			= (1 << 0) | ((GE(par2, 0)) << 1);		\
		} else if (GR(par2, 0)) {				\
		    /* v2 strictly outside --> v1 in-plane is enough,	\
                       already : */					\
		    classification[classnum]				\
			= ((GE(par1, 0)) << 0) | (1 << 1);		\
		} else {						\
		    /* as neither one of them is strictly outside,	\
		     * there's no need for any edge-splitting, at all.	\
		     * */						\
		    classification[classnum] = 0;			\
		}							\
	    } while (0)

	    /* code block used thrice, so put into a macro. The
	     * 'area2D' calls return <0, ==0 or >0 if the three
	     * vertices passed to it are in clockwise order, colinear,
	     * or in conter-clockwise order, respectively. So the sign
	     * of this can be used to determine if a point is inside
	     * or outside the shadow volume. */
#define checkside(a,b,s)					\
	    e_side[0][s] = whichside * area2D((a), (b), v1);	\
	    e_side[1][s] = whichside * area2D((a), (b), v2);	\
	    classify(e_side[0][s], e_side[1][s], s);		\
	    if (classification[s] == 3)				\
		continue;

	    checkside(w1, w2, 0);
	    checkside(w2, w3, 1);
	    checkside(w3, w1, 2);
#undef checkside		/* get rid of that macro, again */


	    /* Test 5 (2D): Does the whole polygon lie on one and the same
	     * side of the tested edge? Again, area2D is used to determine on
	     * which side of a given edge a vertex lies. */
	    /* HBB 20000621: this test only makes sense if this is a
	     * 'real' edge, not just a single point */
	    /* HBB 20051206: but p_side[] may still have to be computed, to be used
	     * by Test 9 */
	    p_side[0] = area2D(v1, v2, w1);
	    p_side[1] = area2D(v1, v2, w2);
	    p_side[2] = area2D(v1, v2, w3);
	    if (v1 != v2) {
		/* HBB 20001104: made this more restrictive: only
		 * reject p if areas are greater than 0, i.e. don't
		 * allow EQ(..,0). Otherwise, edges coincident with a
		 * polygon edge in 2D will not be hidden by it. */
		/* FIXME HBB 20001108: some more code will be needed
		 * here or further down below. It will currently fail
		 * if two edges fall exactly on top of each other
		 * (like the crossover line in the Klein-bottle demo):
		 * both edges will be hidden by the other's polygon,
		 * essentially */
#if 0
		if (0
		    || (GR(p_side[0], 0)
			&& GR(p_side[1], 0)
			&& GR(p_side[2], 0)
			)
		    || (GR(0 , p_side[0])
			&& GR(0 , p_side[1])
			&& GR(0 , p_side[2])
			)
		    )
		    continue;
#else
		/* HBB 20020406: try to repair this */
		if (0
		    || (GE(p_side[0], 0)
			&& GE(p_side[1], 0)
			&& GE(p_side[2], 0)
			)
		    || (GE(0 , p_side[0])
			&& GE(0 , p_side[1])
			&& GE(0 , p_side[2])
			)
		    )
		    continue;
#endif
	    }

	    /* Test 6 (3D): does the whole edge lie on the viewer's
	     * side of the polygon's plane? */
	    v1_rel_pplane = eval_plane_equation(p->plane, v1);
	    v2_rel_pplane = eval_plane_equation(p->plane, v2);

	    /* Condense into classification bit pattern for v1 and v2
	     * wrt.  this plane, and store into entry 'POLY_NVERT' of
	     * the classification array, following the three edges of
	     * the polygon. */
	    classify(v1_rel_pplane, v2_rel_pplane, POLY_NVERT);

	    if (classification[POLY_NVERT] == 3)
		continue;									/* edge completely in front of polygon */

	    /* Test 7 (2D): is edge completely inside the triangle?
	     * The stored results from Test 4 make this rather easy to
	     * check: if a point is oriented equally wrt. all three
	     * edges, and also to the triangle itself, the point is
	     * inside the 2D triangle: */

	    /* search for any one bits (--> they mean 'outside') */
	    if (! (classification[0]
		   | classification[1]
		   | classification[2])) {

		/* Edge is completely inside the polygon's 2D
		 * projection. Unlike Ammeraal, I have to check for
		 * edges penetrating polygons, as well. */

		/* Check if edge really doesn't intersect the triangle
		 * (this happens in roughly one third of the
		 * cases). If both vertices have been classified as
		 * 'behind', the result is true even in the case of
		 * possibly intersecting polygons */
		/* HBB 20020408: softened this check to allow edges
		 * not to be hidden by polygons they're fully inside.
		 * I.e. only declare the edge invisible if at least
		 * one of its endpoints is truly behind, not inside
		 * the polygon's plane. ---> FIXME 20001108 above */
		if (!classification[POLY_NVERT]
		    && (GR(0, v1_rel_pplane) || GR(0, v2_rel_pplane)))
		    return 0;
	    }


	    /* Test 8 (3D): If no edge intersects any polygon, the edge must
	     * lie infront of the poly if the 'inside' vertex in front of it.
	     * This test doesn't work if lines might intersect polygons, so I
	     * took it out of my version of Ammeraal's algorithm. */


	    /* Test 9 (3D): The final 'catch-all' handler. We now know
	     * that the edges passes through the walls of the
	     * semifinite polygonal cylinder of invisible space
	     * stamped out by the polygon, somewhere, so we have to
	     * compute the intersection points with the 2D projection
	     * of the polygon and check whether they are infront of or
	     * behind the polygon, to isolate the visible fragments of
	     * the line. Method is to calc. the [0..1] parameter that
	     * describes the way from V1 to V2, for all intersections
	     * of the edge with any of the planes that define the
	     * 'shadow volume' of this polygon, and deciding what to
	     * do, with all this.
	     *
	     * The remaining work can be classified by the values of
	     * the variables calculated in previous tests:
	     * v{1/2}_rel_pplane, e_side[2][POLY_NVERT],
	     * p_side[POLY_NVERT]. Setting aside cases where one of
	     * them is zero, this leaves a total of 2^11 = 2048
	     * separate cases, in principal. Several of those have
	     * been sorted out, already, by the previous tests. That
	     * leaves us with 3*27 cases, essentially. The '3' cases
	     * differ by which of the three vertices of the polygon is
	     * alone on its side of the test edge.  That already fixes
	     * which two edges might intersect the test edge.  The 27
	     * cases differ by which of the edge points lies on which
	     * of the three planes in space (two sides, and the
	     * polygon plane) that remain interesting for the tests.
	     * That would normally make 2^2 = 4 cases per plane. But
	     * the 'both outside' cases for all three planes have
	     * already been taken care of, earlier, leaving 3 of those
	     * cases to be handled. This is the same for all three
	     * relevant planes, so we get 3^3 = 27 cases. */

	    {
		/* the interception point parameters: */
		double front_hit, this_hit, hit_in_poly;
		long newvert[2];
		int side1, side2;
		double hit1, hit2;
		unsigned int full_class;

		/* find out which of the three edges would be intersected: */
		if ((p_side[0] > 0) == (p_side[1] > 0)) {
		    /* first two vertices on the same side, the third
		     * is on the 'other' side of the test edge */
		    side1 = 1;
		    side2 = 2;
		} else if ((p_side[0] > 0) == (p_side[2] > 0)) {
		    /* vertex '1' on the other side */
		    side1 = 0;
		    side2 = 1;
		} else {
		    /* we now know that it must be the first vertex
		     * that is on the 'other' side, as not all three
		     * can be on different sides. */
		    side1 = 0;
		    side2 = 2;
		}

		/* Carry out the classification into those 27 cases,
		 * based upon classification bits precomputed above,
		 * and calculate the intersection parameters, as
		 * appropriate. Start by regarding the polygon's
		 * plane: '1' means point is in front of plane.
		 *
		 * First, some utility macros: */
#define cleanup_hit(hit)			\
		do {				\
		    if (EQ(hit,0))		\
			hit=0;			\
		    else if (EQ(hit,1))		\
			hit=1;			\
		} while (0)

#define hit_in_edge(hit) ((hit >= 0) && (hit <= 1))

		/* find the intersection point through the front
		 * plane, if any: */
		front_hit = 1e10;
		if (classification[3]) {
		    if (SIGN(v1_rel_pplane - v2_rel_pplane)) {
			front_hit = v1_rel_pplane / (v1_rel_pplane - v2_rel_pplane);
			cleanup_hit(front_hit);
		    }
		    if (!hit_in_edge(front_hit)) {
			front_hit = 1e10;
			classification[3] = 0;	/* not really intersecting */
		    }
		}


		/* Find intersection points with the two relevant side
		 * planes of the shadow volume. A macro contains the
		 * code, to be used for both side1 and side2, to
		 * reduce code overhead: */
#define check_out_hit(hit, side)					\
		do {							\
		    hit = 1e10;	/* default is 'way off' */		\
		    if (classification[side]) {				\
			/* Not both inside, i.e. there has to be an	\
			 * intersection point: */			\
			hit_in_poly = - whichside * p_side[side]	\
			    / (e_side[0][side] - e_side[1][side]);	\
			cleanup_hit(hit_in_poly);			\
			if (hit_in_edge(hit_in_poly))	{		\
			    this_hit = e_side[0][side]			\
				/ (e_side[0][side] - e_side[1][side]);	\
			    cleanup_hit(this_hit);			\
			    if (hit_in_edge(this_hit)) {		\
				hit = this_hit;				\
			    }						\
			}						\
		    }							\
		    if (hit == 1e10)					\
			/* doesn't really intersect */			\
			classification[side] = 0; 			\
		} while (0)

		check_out_hit(hit1, side1);
		check_out_hit(hit2, side2);
#undef check_out_hit

		/* OK, now we know the position of all the hits that
		 * might be needed, given by their positions between
		 * the test edge's vertices. */

		/* Macro: combine all the classifications into a
		 * single classification bitpattern, for easier
		 * listing of cases */
#define makeclass(s2, s1, p) (16 * s2 + 4 * s1 + p)

		full_class = makeclass(classification[side2],
				       classification[side1],
				       classification[3]);

		if (!full_class)
		    /* all suspected intersections cleared, already! */
		    continue;

		/* First sort out all cases where one endpoint is
		 * attributed to more than one intersecting
		 * plane. Only one of those can be 'the real'
		 * intersection of the edge with the surface of the
		 * shadow volume surface. The others are intersections
		 * with the planes defining the volumes, but outside
		 * the actual faces.
		 *
		 * This cuts down the 27 to 13 cases (000, 001, 002,
		 * 012, and their permutations). Start with v1: */
		switch (full_class & makeclass(1, 1, 1)) {
		case makeclass(0,1,1):
		    /* v1 is out 'side1' and infront of 'p' */
		    if (front_hit < hit1)
			classification[3] = 0;
		    else
			classification[side1] = 0;
		    break;

		case makeclass(1,0,1):
		    /* v1 is out 'side2' and infront of 'p' */
		    if (front_hit < hit2)
			classification[3] = 0;
		    else
			classification[side2] = 0;
		    break;

		    /* The following two cases never trigger, even in
		     * a full run through all.dem. Neither do the
		     * corresponding 2 cases for the vertex2
		     * classification. I think this is so because the
		     * determination of 'side1' and 'side2' already
		     * took care of this: if the edge vertex is out
		     * two sides, those two will _not_ be chosen as
		     * 'side1' and 'side2', because one of them is
		     * completely on one side of the edge.  */
		case makeclass(1,1,0):
		    /* v1 out both sides, but behind the front */
		    if (hit1 < hit2)
			/* hit2 is closer to the hidden vertex v2, so
			 * it's the relevant intersection: */
			classification[side1] = 0;
		    else
			classification[side2] = 0;
		    break;

		case makeclass(1,1,1):
		    /* v1 out both sides, and in front of p (--> v2 is
		     * hidden) */
		    if (front_hit < hit1) {
			classification[3] = 0;
			if (hit1 < hit2)
			    classification[side1] = 0;
			else
			    classification[side2] = 0;
		    } else {
			classification[side1] = 0;
			if (front_hit < hit2)
			    classification[3] = 0;
			else
			    classification[side2] = 0;
		    }
		    break;

		default:
		    ; /* do nothing */
		} /* switch(v1 classes) */


		/* And the same, for v2 */
		switch (full_class & makeclass(2, 2, 2)) {
		case makeclass(0,2,2):
		    /* v2 is out 'side1' and infront of 'p' */
		    if (front_hit > hit1)
			classification[3] = 0;
		    else
			classification[side1] = 0;
		    break;

		case makeclass(2,0,2):
		    /* v2 out side 'side2' and infront of 'p' */
		    if (front_hit > hit2)
			classification[3] = 0;
		    else
			classification[side2] = 0;
		    break;

		    /* See note in v1 handling about these two cases
		     * being impossible... */
		case makeclass(2,2,0):
		    /* v2 out both sides, but behind the front */
		    if (hit1 > hit2)
			/* hit2 is closer to the hidden vertex v2, so
			 * it's the relevant intersection: */
			classification[side1] = 0;
		    else
			classification[side2] = 0;
		    break;

		case makeclass(2,2,2):
		    /* v2 out both sides, and in front of p (--> v1 is
		     * hidden) */
		    if (front_hit > hit1) {
			classification[3] = 0;
			if (hit1 > hit2)
			    classification[side1] = 0;
			else
			    classification[side2] = 0;
		    } else {
			classification[side1] = 0;
			if (front_hit > hit2)
			    classification[3] = 0;
			else
			    classification[side2] = 0;
		    }
		    break;

		default:
		    ; /* do nothing */
		} /* switch(v2 classes) */


		/* recalculate full_class after all these possible changes: */
		full_class = makeclass(classification[side2],
				       classification[side1],
				       classification[3]);

		/* This monster switch catches all the 13 remaining cases
		 * individually. */
		switch (full_class) {
		default:
		    printf("in_front: (T9) class %d is illegal. Should never happen.\n",
			   full_class);
		    break;
		case makeclass(0,0,0):
		    /* can happen, by resetting of classification bits
		     * inside this test */
		    break;

/*--------- The simplest cases first: only one intersection point -----*/
		    /* Code is the same for all six cases, or nearly so: */
#define handle_singleplane_hit(vert, hitpar)				    \
		    newvert[0] = split_line_at_ratio(vnum1, vnum2, hitpar); \
		    setup_edge(vert, newvert[0]);			    \
		    break;

		case makeclass(0,0,1): /* v1 is in front of poly: */
		    handle_singleplane_hit(vnum1, front_hit);
		case makeclass(0,0,2): /* v2 is in front of poly: */
		    handle_singleplane_hit(vnum2, front_hit);
		case makeclass(0,1,0): /* v1 is out side1: */
		    handle_singleplane_hit(vnum1,hit1);
		case makeclass(0,2,0): /* v2 is out side1: */
		    handle_singleplane_hit(vnum2,hit1);
		case makeclass(1,0,0): /* v1 is out side2 */
		    handle_singleplane_hit(vnum1,hit2);
		case makeclass(2,0,0): /* v2 is out side2 */
		    handle_singleplane_hit(vnum2,hit2);

/*----------- cases with 2 certain intersections: --------------*/
		case makeclass(1,2,0):
		case makeclass(2,1,0):
		    /* Both behind the front, v1 out one side, v2 out
		     * the other.  These two cases can be handled by
		     * common code: */
		    newvert[0] = split_line_at_ratio(vnum1, vnum2, hit1);
		    newvert[1] = split_line_at_ratio(vnum1, vnum2, hit2);
		    if (hit1 < hit2) {
			/* the fragments are v1--hit1 and hit2--v2 */
			handle_edge_fragment(edgenum, newvert[1], vnum2, polynum);
			setup_edge(vnum1, newvert[0]);
		    } else {
			/* the fragments are v2--hit1 and hit2--v1 */
			handle_edge_fragment(edgenum, vnum1, newvert[1], polynum);
			setup_edge(newvert[0], vnum2);
		    }
		    break;

/*----------- cases with either 2 or no intersections: --------------*/

		    /* Mainly identical code block to be used 4 times
		     * --> macro */
		    /* HBB 20001108: up until today, this part of the
		     * was severely buggy. But I do think I've got it
		     * fixed up, this time */
#define handle_outside_behind_vs_infront(hitvar1, hitvar2)			   \
		    /* vnum1 out plane of 'hitvar1' and in back, vnum2		   \
                     * in front */						   \
		    if ((hitvar1) < (hitvar2)) {				   \
			newvert[0] = split_line_at_ratio(vnum1, vnum2, (hitvar1)); \
			newvert[1] = split_line_at_ratio(vnum1, vnum2, (hitvar2)); \
			handle_edge_fragment(edgenum, newvert[1], vnum2, polynum); \
			setup_edge(vnum1, newvert[0]);				   \
		    }								   \
		    /* otherwise, the edge missed the shadow volume		   \
		     * --> nothing to do */					   \
		    break;

		case makeclass(0,1,2): /* v1 out side1 and in back, v2 in front */
		    handle_outside_behind_vs_infront(hit1, front_hit);

		case makeclass(0,2,1): /* v2 out side1 and in back, v1 in front */
		    handle_outside_behind_vs_infront(front_hit, hit1);

		case makeclass(1,0,2): /* v1 out side2 and in back, v2 in front */
		    handle_outside_behind_vs_infront(hit2, front_hit);

		case makeclass(2,0,1): /* v2 out side2 and in back, v1 in front */
		    handle_outside_behind_vs_infront(front_hit, hit2);

#undef handle_outside_behind_vs_infront

		} /* end of switch through the cases */

	    } /* end of part 'T9' */
	} /* for (polygons in list) */

    /* Came here, so there's something left of this edge, which needs
     * to be drawn.  But the vertices are different, now, so copy our
     * new vertices back into 'e' */

    draw_edge(elist + edgenum, vlist + vnum1, vlist + vnum2);

    return 1;
}


/* HBB 20000617: reimplemented this routine from scratch */
/* Externally callable function to draw a line, but hide it behind the
 * visible surface. */
/* NB: The p_vertex arguments are not allowed to be pointers into the
 * hidden3d 'vlist' structure. If they are, they may become invalid
 * before they're used, because of the nextfrom_dynarray() call. */
void draw_line_hidden( p_vertex v1, p_vertex v2,	/* pointers to the end vertices */
		struct lp_style_type *lp)	/* line and point style to draw in */
{
    long int vstore1, vstore2;
    long int edgenum;
    long int temp_pfirst;

    /* If no polygons have been stored, nothing can be hidden, and we
     * can't use in_front() because the datastructures are partly
     * invalid. So just draw the line and be done with it */
    if (!polygons.end) {
	draw3d_line_unconditional(v1, v2, lp, lp->l_type);
	return;
    }

    /* Copy two vertices into hidden3d arrays: */
    nextfrom_dynarray(&vertices);
    vstore1 = vertices.end - 1;
    vlist[vstore1] = *v1;
    if (v2) {
	vlist[vstore1].lp_style = NULL;
	nextfrom_dynarray(&vertices);
	vstore2 = vertices.end - 1;
	vlist[vstore2] = *v2;
	vlist[vstore2].lp_style = NULL;
    } else {
	/* v2 == NULL --> this is a point symbol to be drawn. Make two
	 * vertex pointers the same, and set up the 'style' field */
	vstore2 = vstore1;
	vlist[vstore2].lp_style = lp;
    }

    /* store the edge into the hidden3d datastructures */
    edgenum = make_edge(vstore1, vstore2, lp, lp->l_type, -1);

    /* remove hidden portions of the line, and draw what remains */
    temp_pfirst = pfirst;
    in_front(edgenum, elist[edgenum].v1, elist[edgenum].v2, &temp_pfirst);

    /* release allocated storage slots: */
    droplast_dynarray(&edges);
    droplast_dynarray(&vertices);
    if (v2)
	droplast_dynarray(&vertices);
}

/***********************************************************************
 * and, finally, the 'mother function' that uses all these lots of tools
 ***********************************************************************/
void
plot3d_hidden(struct surface_points *plots, int pcount)
{
    /* make vertices, edges and polygons out of all the plots */
    build_networks(plots, pcount);

    if (! edges.end) {
	/* No drawable edges found. Free all storage and bail out. */
	term_hidden_line_removal();
	graph_error("*All* edges undefined or out of range, thus no plot.");
    }

    if (! polygons.end) {
	/* No polygons anything could be hidden behind... */

	sort_edges_by_z();
	while (efirst >= 0) {
	    draw_edge(elist+efirst, vlist + elist[efirst].v1, vlist + elist[efirst].v2);
	    efirst = elist[efirst].next;
	}
    } else {
	long int temporary_pfirst;

	/* Presort edges in z order */
	sort_edges_by_z();
	/* Presort polygons in z order */
	sort_polys_by_z();

	temporary_pfirst = pfirst;

	while (efirst >=0) {
	    if (elist[efirst].style >= -2) /* skip invisible edges */
		in_front(efirst, elist[efirst].v1, elist[efirst].v2,
			 &temporary_pfirst);
	    efirst = elist[efirst].next;
	}
    }

    /* Free memory */
    /* FIXME: anything to free? */
}

void
reset_hidden3doptions()
{
    hiddenBacksideLinetypeOffset = BACKSIDE_LINETYPE_OFFSET;
    hiddenTriangleLinesdrawnPattern = TRIANGLE_LINESDRAWN_PATTERN;
    hiddenHandleUndefinedPoints = HANDLE_UNDEFINED_POINTS;
    hiddenShowAlternativeDiagonal = SHOW_ALTERNATIVE_DIAGONAL;
    hiddenHandleBentoverQuadrangles = HANDLE_BENTOVER_QUADRANGLES;
}

/* Emacs editing help for HBB:
 * Local Variables: ***
 * c-basic-offset: 4 ***
 * End: ***
 */
