--- /dev/null
+#ifndef lint
+static char *RCSid() { return RCSid("$Id: hidden3d.c,v 1.56.2.3 2008/06/02 19:40:10 sfeam Exp $"); }
+#endif
+
+/* 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 "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
+
+/* The code used to die messily if the scale parameters got over-large.
+ * Prevent this from happening due to mousing by locking out the mouse
+ * response. */
+TBOOLEAN disable_mouse_z = FALSE;
+
+/* 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 routine to calculate the cells' position in that array: */
+static int
+coord_to_treecell(coordval x)
+{
+ int index;
+ index = ((((x) / surface_scale) + 1.0) / 2.0) * QUADTREE_GRANULARITY;
+ if (index >= QUADTREE_GRANULARITY)
+ index = QUADTREE_GRANULARITY - 1;
+ else if (index < 0)
+ index = 0;
+
+ return index;
+}
+
+/* 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 __PROTO((struct coordinate GPHUGE *point,
+ lp_style_type *lp_style, TBOOLEAN color_from_column));
+static long int make_edge __PROTO((long int vnum1, long int vnum2,
+ struct lp_style_type *lp,
+ int style, int next));
+static long int store_edge __PROTO((long int vnum1, edge_direction direction,
+ long int crvlen, struct lp_style_type *lp,
+ int style));
+static GP_INLINE double eval_plane_equation __PROTO((t_plane p, p_vertex v));
+static long int store_polygon __PROTO((long int vnum1,
+ polygon_direction direction,
+ long int crvlen));
+static void color_edges __PROTO((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 __PROTO((struct surface_points * plots,
+ int pcount));
+int compare_edges_by_zmin __PROTO((SORTFUNC_ARGS p1, SORTFUNC_ARGS p2));
+int compare_polys_by_zmax __PROTO((SORTFUNC_ARGS p1, SORTFUNC_ARGS p2));
+static void sort_edges_by_z __PROTO((void));
+static void sort_polys_by_z __PROTO((void));
+static TBOOLEAN get_plane __PROTO((p_polygon p, t_plane plane));
+static long split_line_at_ratio __PROTO((long int vnum1, long int vnum2, double w));
+static GP_INLINE double area2D __PROTO((p_vertex v1, p_vertex v2,
+ p_vertex v3));
+static void draw_vertex __PROTO((p_vertex v));
+static GP_INLINE void draw_edge __PROTO((p_edge e, p_vertex v1, p_vertex v2));
+static GP_INLINE void handle_edge_fragment __PROTO((long int edgenum,
+ long int vnum1,
+ long int vnum2,
+ long int firstpoly));
+static int in_front __PROTO((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()
+{
+ struct value t;
+ 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) real(const_express(&t));
+ c_token--;
+ break;
+ case S_HI_NOOFFSET:
+ hiddenBacksideLinetypeOffset = 0;
+ break;
+ case S_HI_TRIANGLEPATTERN:
+ c_token++;
+ hiddenTriangleLinesdrawnPattern = (int) real(const_express(&t));
+ c_token--;
+ break;
+ case S_HI_UNDEFINED:
+ c_token++;
+ tmp = (int) real(const_express(&t));
+ 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 = 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 = 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);
+
+ /* Check if polygon is degenerate */
+ if (V_EQUAL(v1,v2) || V_EQUAL(v2,v3) || V_EQUAL(v3,v1))
+ return (-2);
+
+ /* All else OK, fill in the polygon: */
+
+ p = nextfrom_dynarray(&polygons);
+
+ memcpy (p->vertex, v, sizeof(v));
+#if ! HIDDEN3D_QUADTREE
+ p->next = -1;
+#endif
+
+ /* Some helper macros for repeated 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; \
+ if (min < -surface_scale) disable_mouse_z = TRUE; \
+ } 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; \
+ if (max > surface_scale) disable_mouse_z = TRUE; \
+ } 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 = gp_alloc(2 * max_crvlen * sizeof (long int),
+ "hidden north_polys");
+ these_polygons = gp_alloc(2 * max_crvlen * sizeof (long int),
+ "hidden these_polys");
+ north_edges = gp_alloc(3 * max_crvlen * sizeof (long int),
+ "hidden north_edges");
+ these_edges = 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 = {0};
+
+ 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 = 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 = 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: */
+ {
+ int grid_x, grid_y;
+ 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 = 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
+ int grid_x, grid_y;
+ int grid_x_low, grid_x_high;
+ 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: ***
+ */
projects / gnuplot / blobdiff