summaryrefslogtreecommitdiff
path: root/graphics/sketch/bsp.c
diff options
context:
space:
mode:
authorNorbert Preining <norbert@preining.info>2019-09-02 13:46:59 +0900
committerNorbert Preining <norbert@preining.info>2019-09-02 13:46:59 +0900
commite0c6872cf40896c7be36b11dcc744620f10adf1d (patch)
tree60335e10d2f4354b0674ec22d7b53f0f8abee672 /graphics/sketch/bsp.c
Initial commit
Diffstat (limited to 'graphics/sketch/bsp.c')
-rw-r--r--graphics/sketch/bsp.c1237
1 files changed, 1237 insertions, 0 deletions
diff --git a/graphics/sketch/bsp.c b/graphics/sketch/bsp.c
new file mode 100644
index 0000000000..866a6b42d6
--- /dev/null
+++ b/graphics/sketch/bsp.c
@@ -0,0 +1,1237 @@
+/* bsp.c
+ Copyright (C) 2005,2006,2007 Eugene K. Ressler, Jr.
+
+This file is part of Sketch, a small, simple system for making
+3d drawings with LaTeX and the PSTricks or TikZ package.
+
+Sketch is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 3, or (at your option)
+any later version.
+
+Sketch is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with Sketch; see the file COPYING.txt. If not, see
+http://www.gnu.org/copyleft */
+
+#include <stdio.h>
+#include "bsp.h"
+#include "geomio.h"
+
+DECLARE_DYNAMIC_ARRAY_FUNCS (BSP_POLYGON_ATTR, BSP_VERTEX_ATTR,
+ polygon_attr, elt, n_elts, NO_OTHER_INIT)
+// ---- polylines --------------------------------------------------------------
+ static BSP_POLYLINE_NODE *new_bsp_polyline_node (void *attr)
+{
+ BSP_POLYLINE_NODE *n = safe_malloc (sizeof *n);
+ n->tag = BSP_POLYLINE;
+ n->attr = attr;
+ n->prev = n->next = n->mark = n->in = n->on = n->out = NULL;
+ n->first_p = n->last_p = 0;
+ init_box_3d (n->extent);
+ init_polyline_3d (n->polyline);
+ return n;
+}
+
+static void
+delete_bsp_polyline_node (BSP_POLYLINE_NODE * node)
+{
+ clear_polyline_3d (node->polyline);
+ safe_free (node);
+}
+
+static void
+set_bsp_polyline_extent (BSP_POLYLINE_NODE * node)
+{
+ // set extent
+ init_box_3d (node->extent);
+ fold_min_max_polyline_3d (node->extent, node->polyline);
+}
+
+static void
+init_bsp_polyline (BSP_POLYLINE_NODE * node, POLYLINE_3D * polyline)
+{
+ // initial polyline contains first and last points; split ones don't
+ node->first_p = node->last_p = 1;
+ copy_polyline_3d (node->polyline, polyline);
+ set_bsp_polyline_extent (node);
+}
+
+static void
+split_polyline_with_plane (BSP_POLYLINE_NODE * node,
+ PLANE * plane,
+ BSP_POLYLINE_NODE ** in_nodes,
+ BSP_POLYLINE_NODE ** on_nodes,
+ BSP_POLYLINE_NODE ** out_nodes)
+{
+ int i, j, i_side, current_side;
+ BSP_POLYLINE_NODE *current, **list;
+ VECTOR_3D v, dp;
+ POINT_3D isect;
+ FLOAT t;
+
+ // initialize all the output lists to empty
+ *in_nodes = *on_nodes = *out_nodes = NULL;
+
+ // initialize the current polyline that's being "built", copying attributes
+ // of the original
+ current = new_bsp_polyline_node (node->attr);
+
+ // copy source polygon's first_ attribute
+ current->first_p = node->first_p;
+
+ // j is the trail index as we step through vertices with head i
+ j = 0;
+
+ // copy first vertex of polyline onto current output polyline
+ copy_pt_3d (pushed_polyline_3d_v (current->polyline), node->polyline->v[j]);
+
+ // find side of cutting plane that first vertex is on.
+ current_side =
+ (S_IN | S_ON | S_OUT) & pt_side_of_plane (plane, node->polyline->v[j]);
+
+ // loop through all directed segments j->i
+ for (i = 1; i < node->polyline->n_vertices; j = i++)
+ {
+ i_side =
+ (S_IN | S_ON | S_OUT) & pt_side_of_plane (plane,
+ node->polyline->v[i]);
+
+#define HASH(Current, I) ((Current << 3) | I)
+
+ switch (HASH (current_side, i_side))
+ {
+
+ case HASH (S_OUT, S_IN):
+ case HASH (S_IN, S_OUT):
+
+ // vertices straddle the plane; compute the intersection
+ sub_pts_3d (v, node->polyline->v[i], node->polyline->v[j]); // direction vector
+ sub_pts_3d (dp, plane->p, node->polyline->v[j]); // P - L
+ t = dot_3d (dp, plane->n) / dot_3d (v, plane->n); // parameter of intersect
+ add_scaled_vec_to_pt_3d (isect, node->polyline->v[j], v, t); // intersect
+
+ // finish current polyline and add to current list
+ copy_pt_3d (pushed_polyline_3d_v (current->polyline), isect);
+ set_bsp_polyline_extent (current);
+ list = current_side == S_IN ? in_nodes : out_nodes;
+ current->in = (BSP_NODE *) * list;
+ *list = current;
+
+ // start a new one by adding the isect and new point
+ current = new_bsp_polyline_node (node->attr);
+ copy_pt_3d (pushed_polyline_3d_v (current->polyline), isect);
+ copy_pt_3d (pushed_polyline_3d_v (current->polyline),
+ node->polyline->v[i]);
+ current_side = i_side;
+ break;
+
+ case HASH (S_OUT, S_ON):
+ case HASH (S_IN, S_ON):
+
+ // line that was away from the plane joins it;
+ // finish current polyline and add to current list
+ copy_pt_3d (pushed_polyline_3d_v (current->polyline),
+ node->polyline->v[i]);
+ set_bsp_polyline_extent (current);
+ list = current_side == S_IN ? in_nodes : out_nodes;
+ current->in = (BSP_NODE *) * list;
+ *list = current;
+
+ // start a new one if there are still vertices left to process
+ if (i < node->polyline->n_vertices - 1)
+ {
+ current = new_bsp_polyline_node (node->attr);
+ copy_pt_3d (pushed_polyline_3d_v (current->polyline),
+ node->polyline->v[i]);
+ current_side = S_ON;
+ }
+ else
+ {
+ // copy last_p attribute from source
+ current->last_p = node->last_p;
+ current = NULL;
+ current_side = 0;
+ }
+ break;
+
+ case HASH (S_ON, S_OUT):
+ case HASH (S_ON, S_IN):
+
+ // line that was on the plane springs away from it;
+ // if there is more than one point in the current polyline,
+ // complete it and start a new one
+ if (current->polyline->n_vertices > 1)
+ {
+ FLOAT *last_vertex = current->polyline->v[current->polyline->n_vertices - 1]; //remember last vertex
+ set_bsp_polyline_extent (current);
+ current->in = (BSP_NODE *) * on_nodes;
+ *on_nodes = current;
+ current = new_bsp_polyline_node (node->attr);
+ copy_pt_3d (pushed_polyline_3d_v (current->polyline),
+ last_vertex);
+ }
+ // add the new vertex to the current polyline
+ copy_pt_3d (pushed_polyline_3d_v (current->polyline),
+ node->polyline->v[i]);
+ current_side = i_side; // now either in or out
+ break;
+
+ default:
+
+ // nothing has changed, so just add the new point
+ // to the ccurrent polygon
+ copy_pt_3d (pushed_polyline_3d_v (current->polyline),
+ node->polyline->v[i]);
+ break;
+ }
+ }
+ // finish the final polyline
+ if (current)
+ {
+ // copy last_p attribute from source
+ current->last_p = node->last_p;
+
+ set_bsp_polyline_extent (current);
+ list = (current_side & S_IN) ? in_nodes :
+ (current_side & S_OUT) ? out_nodes : on_nodes;
+ current->in = (BSP_NODE *) * list;
+ *list = current;
+ }
+}
+
+static void
+insert_polyline (BSP_TREE * bsp, BSP_POLYLINE_NODE * node)
+{
+ if (*bsp == NULL)
+ {
+ *bsp = (BSP_NODE *) node;
+ }
+ else if ((*bsp)->tag == BSP_POLYLINE)
+ {
+ BSP_POLYLINE_NODE *bsp_node = (BSP_POLYLINE_NODE *) * bsp;
+#ifdef EXPERIMENTAL_OPTIMIZATION
+ node->in = bsp_node;
+ *bsp = (BSP_NODE *) node;
+#else
+ insert_polyline ((BSP_TREE *) & bsp_node->in, node);
+#endif
+ }
+ else
+ { // BSP_POLYGON
+ BSP_POLYGON_NODE *bsp_node = (BSP_POLYGON_NODE *) * bsp;
+ int side = polyline_side_of_plane (node->polyline, bsp_node->plane);
+ if (side == S_IN)
+ insert_polyline (&bsp_node->in, node);
+ else if (side == S_ON)
+ insert_polyline (&bsp_node->on, node);
+ else if (side == S_OUT)
+ insert_polyline (&bsp_node->out, node);
+ else
+ { // S_SPLIT
+ BSP_POLYLINE_NODE *in_polylines, *on_polylines, *out_polylines,
+ *p, *p_next;
+ split_polyline_with_plane (node, bsp_node->plane, &in_polylines,
+ &on_polylines, &out_polylines);
+ // remove polylines from lists and add to respective bsp subtrees recursively
+ for (p = in_polylines; p; p = p_next)
+ {
+ p_next = (BSP_POLYLINE_NODE *) p->in;
+ p->in = NULL;
+ insert_polyline (&bsp_node->in, p);
+ }
+ for (p = on_polylines; p; p = p_next)
+ {
+ p_next = (BSP_POLYLINE_NODE *) p->in;
+ p->in = NULL;
+ insert_polyline (&bsp_node->on, p);
+ }
+ for (p = out_polylines; p; p = p_next)
+ {
+ p_next = (BSP_POLYLINE_NODE *) p->in;
+ p->in = NULL;
+ insert_polyline (&bsp_node->out, p);
+ }
+ delete_bsp_polyline_node (node);
+ }
+ }
+}
+
+void
+add_polyline_to_bsp (BSP_TREE * bsp, POLYLINE_3D * polyline, void *attr)
+{
+ BSP_POLYLINE_NODE *node = new_bsp_polyline_node (attr);
+ init_bsp_polyline (node, polyline);
+ insert_polyline (bsp, node);
+}
+
+// ---- polygons ---------------------------------------------------------------
+
+static BSP_POLYGON_NODE *
+new_bsp_polygon_node (void *attr)
+{
+ BSP_POLYGON_NODE *n = safe_malloc (sizeof *n);
+ n->tag = BSP_POLYGON;
+ n->attr = attr;
+ n->prev = n->next = n->mark = n->in = n->on = n->out = NULL;
+ init_box_3d (n->extent);
+ init_polygon_3d (n->polygon);
+ init_polygon_attr (n->polygon_attr);
+ return n;
+}
+
+static void
+set_bsp_polygon_extent (BSP_POLYGON_NODE * node)
+{
+ // set extent
+ init_box_3d (node->extent);
+ fold_min_max_polygon_3d (node->extent, node->polygon);
+}
+
+static void
+init_bsp_polygon_node (BSP_POLYGON_NODE * node, POLYGON_3D * polygon)
+{
+ int i;
+
+ // fill in the polygon verticies
+ copy_polygon_3d (node->polygon, polygon);
+
+ // fill in the plane
+ find_polygon_plane (node->plane, polygon);
+
+ // mark all edges as border edges
+ setup_polygon_attr (node->polygon_attr, polygon->n_sides);
+ for (i = 0; i < polygon->n_sides; i++)
+ node->polygon_attr->elt[i].border_p = 1;
+ node->polygon_attr->n_elts = polygon->n_sides;
+
+ // fill in the extent
+ set_bsp_polygon_extent (node);
+}
+
+static void
+delete_bsp_polygon_node (BSP_POLYGON_NODE * node)
+{
+ clear_polygon_3d (node->polygon);
+ clear_polygon_attr (node->polygon_attr);
+ safe_free (node);
+}
+
+// decide whether a j->i edge of the the new polygon that has
+// been split from parent is part of the border.
+static int
+is_new_border_p (BSP_VERTEX_ATTR * i_attr,
+ BSP_VERTEX_ATTR * j_attr,
+ BSP_POLYGON_ATTR * parent_attr, int parent_n_sides)
+{
+ int parent_is_border_p;
+
+ if (parent_attr->n_elts != parent_n_sides)
+ die (no_line, "failed assumption on attribute size");
+ parent_is_border_p = parent_attr->elt[j_attr->parent_vtx].border_p;
+ if (!parent_is_border_p)
+ return 0;
+
+ if (i_attr->cut_p)
+ {
+ if (j_attr->cut_p)
+ {
+ return 0;
+ }
+ else
+ {
+ return i_attr->parent_vtx == j_attr->parent_vtx;
+ }
+ }
+ return (i_attr->parent_vtx - j_attr->parent_vtx +
+ parent_n_sides) % parent_n_sides == 1;
+}
+
+static void
+decide_boundaries (BSP_POLYGON_NODE * new_node, BSP_POLYGON_NODE * node)
+{
+ int i, j, last_border_p;
+
+ i = 0;
+ j = new_node->polygon->n_sides - 1;
+ last_border_p =
+ is_new_border_p (&new_node->polygon_attr->elt[i],
+ &new_node->polygon_attr->elt[j],
+ node->polygon_attr, node->polygon->n_sides);
+ for (j = i++; i < new_node->polygon->n_sides; j = i++)
+ {
+ new_node->polygon_attr->elt[j].border_p =
+ is_new_border_p (&new_node->polygon_attr->elt[i],
+ &new_node->polygon_attr->elt[j],
+ node->polygon_attr, node->polygon->n_sides);
+ }
+ new_node->polygon_attr->elt[j].border_p = last_border_p;
+}
+
+static void
+push_polygon_attr (BSP_POLYGON_NODE * node, int parent_vtx, int cut_p)
+{
+ BSP_VERTEX_ATTR *vertex_attr = pushed_polygon_attr_elt (node->polygon_attr);
+ vertex_attr->border_p = 0;
+ vertex_attr->parent_vtx = parent_vtx;
+ vertex_attr->cut_p = cut_p;
+}
+
+static void
+split_polygon_with_plane (BSP_POLYGON_NODE * node,
+ PLANE * plane,
+ BSP_POLYGON_NODE * in_node,
+ BSP_POLYGON_NODE * out_node)
+{
+ int i, j, i_side, j_side;
+ VECTOR_3D v, dp;
+ POINT_3D isect;
+ FLOAT t;
+
+ // reset fill pointers
+ in_node->polygon->n_sides = out_node->polygon->n_sides = 0;
+ in_node->polygon_attr->n_elts = out_node->polygon_attr->n_elts = 0;
+
+ // process all pairs of vertices
+ j = node->polygon->n_sides - 1;
+ i_side = pt_side_of_plane (plane, node->polygon->v[j]);
+ for (i = 0; i < node->polygon->n_sides; j = i++)
+ {
+ j_side = i_side;
+ i_side = pt_side_of_plane (plane, node->polygon->v[i]);
+
+ if ((i_side | j_side) == (S_IN | S_OUT))
+ {
+
+ // the two most recent points straddle the the plane
+ // compute the intersection
+ sub_pts_3d (v, node->polygon->v[i], node->polygon->v[j]); // direction vector
+ sub_pts_3d (dp, plane->p, node->polygon->v[j]); // P - L
+ t = dot_3d (dp, plane->n) / dot_3d (v, plane->n); // parameter of intersect
+ add_scaled_vec_to_pt_3d (isect, node->polygon->v[j], v, t); // intersect
+
+ // put intersect in both polygons
+ copy_pt_3d (pushed_polygon_3d_v (in_node->polygon), isect);
+ copy_pt_3d (pushed_polygon_3d_v (out_node->polygon), isect);
+
+ if (i_side == S_IN)
+ {
+ // edge traversed from outside to in
+ push_polygon_attr (out_node, j, 1);
+ push_polygon_attr (in_node, j, 1);
+ copy_pt_3d (pushed_polygon_3d_v (in_node->polygon),
+ node->polygon->v[i]);
+ push_polygon_attr (in_node, i, 0);
+ }
+ else
+ {
+ // edge traversed from inside to out
+ push_polygon_attr (out_node, j, 1);
+ push_polygon_attr (in_node, j, 1);
+ copy_pt_3d (pushed_polygon_3d_v (out_node->polygon),
+ node->polygon->v[i]);
+ push_polygon_attr (out_node, i, 0);;
+ }
+ }
+ else if (i_side & S_ON)
+ {
+
+ // the current vertex is on the plane
+ // put it in both polygons
+ copy_pt_3d (pushed_polygon_3d_v (in_node->polygon),
+ node->polygon->v[i]);
+ copy_pt_3d (pushed_polygon_3d_v (out_node->polygon),
+ node->polygon->v[i]);
+ push_polygon_attr (in_node, i, 0);
+ push_polygon_attr (out_node, i, 0);
+ }
+ else
+ {
+
+ // the new vertex is not straddling nor on the plane
+ // so we output it to the correct polygon
+ if (i_side == S_IN)
+ {
+ copy_pt_3d (pushed_polygon_3d_v (in_node->polygon),
+ node->polygon->v[i]);
+ push_polygon_attr (in_node, i, 0);
+ }
+ else
+ {
+ copy_pt_3d (pushed_polygon_3d_v (out_node->polygon),
+ node->polygon->v[i]);
+ push_polygon_attr (out_node, i, 0);
+ }
+ }
+ }
+ // fill in the planes
+ find_polygon_plane (in_node->plane, in_node->polygon);
+ find_polygon_plane (out_node->plane, out_node->polygon);
+
+ // set up extents
+ set_bsp_polygon_extent (in_node);
+ set_bsp_polygon_extent (out_node);
+
+ // make a pass around the in and out polygons to fill in boundardy_p
+ decide_boundaries (in_node, node);
+ decide_boundaries (out_node, node);
+}
+
+static void
+insert_polygon (BSP_TREE * bsp, BSP_POLYGON_NODE * node)
+{
+ if (*bsp == NULL)
+ *bsp = (BSP_NODE *) node;
+ else
+ {
+ BSP_POLYGON_NODE *bsp_node = (BSP_POLYGON_NODE *) * bsp;
+ int side = polygon_side_of_plane (node->polygon, bsp_node->plane);
+ if (side & (S_IN | S_ON))
+ insert_polygon (&bsp_node->in, node);
+ else if (side == S_OUT)
+ insert_polygon (&bsp_node->out, node);
+ else
+ {
+ BSP_POLYGON_NODE *in_node = new_bsp_polygon_node (node->attr);
+ BSP_POLYGON_NODE *out_node = new_bsp_polygon_node (node->attr);
+ split_polygon_with_plane (node, bsp_node->plane, in_node, out_node);
+ insert_polygon (&bsp_node->in, in_node);
+ insert_polygon (&bsp_node->out, out_node);
+ delete_bsp_polygon_node (node);
+ }
+ }
+}
+
+void
+add_polygon_to_bsp (BSP_TREE * bsp, POLYGON_3D * polygon, void *attr)
+{
+ BSP_POLYGON_NODE *node = new_bsp_polygon_node (attr);
+ init_bsp_polygon_node (node, polygon);
+ insert_polygon (bsp, node);
+}
+
+static struct
+{
+ BSP_NODE_FUNC func;
+ void *env;
+}
+traverse_closure[1];
+
+static void
+walk_bsp (BSP_NODE * bsp)
+{
+ if (bsp == NULL)
+ return;
+ if (bsp->tag == BSP_POLYGON)
+ {
+ BSP_POLYGON_NODE *node = (BSP_POLYGON_NODE *) bsp;
+ if (node->plane->n[Z] < 0)
+ {
+ walk_bsp (node->out);
+ traverse_closure->func (bsp, traverse_closure->env);
+ walk_bsp (node->on);
+ walk_bsp (node->in);
+ }
+ else
+ {
+ walk_bsp (node->in);
+ traverse_closure->func (bsp, traverse_closure->env);
+ walk_bsp (node->on);
+ walk_bsp (node->out);
+ }
+ }
+ else
+ { // BSP_POLYLINE
+ BSP_POLYLINE_NODE *node = (BSP_POLYLINE_NODE *) bsp;
+ traverse_closure->func (bsp, traverse_closure->env);
+ walk_bsp ((BSP_NODE *) node->in);
+ }
+}
+
+void
+traverse_bsp (BSP_NODE * bsp, BSP_NODE_FUNC func, void *env)
+{
+ traverse_closure->func = func;
+ traverse_closure->env = env;
+ walk_bsp (bsp);
+}
+
+void
+traverse_depth_sort (BSP_NODE * bsp, BSP_NODE_FUNC func, void *env)
+{
+ traverse_closure->func = func;
+ traverse_closure->env = env;
+ for (; bsp; bsp = bsp->next)
+ walk_bsp (bsp);
+}
+
+static void
+indent (FILE * f, int n)
+{
+ while (n-- > 0)
+ fprintf (f, " ");
+}
+
+static void
+print_bsp_internal (FILE * f, BSP_NODE * bsp, int n)
+{
+ if (bsp == NULL)
+ return;
+
+ indent (f, 2 * n);
+ if (bsp->tag == BSP_POLYGON)
+ {
+
+ BSP_POLYGON_NODE *node = (BSP_POLYGON_NODE *) bsp;
+ fprintf (f, "BSPpolygon\n");
+
+ indent (f, 2 * n + 1);
+ print_plane (f, node->plane);
+
+ indent (f, 2 * n + 1);
+ print_polygon_3d (f, node->polygon);
+
+ indent (f, 2 * n + 1);
+ fprintf (f, "in\n");
+ print_bsp_internal (f, node->in, n + 1);
+
+ indent (f, 2 * n + 1);
+ fprintf (f, "on\n");
+ print_bsp_internal (f, node->on, n + 1);
+
+ indent (f, 2 * n + 1);
+ fprintf (f, "out\n");
+ print_bsp_internal (f, node->out, n + 1);
+ }
+ else
+ { // BSP_POLYLINE
+ BSP_POLYLINE_NODE *node = (BSP_POLYLINE_NODE *) bsp;
+ fprintf (f, "BSPpolyline ");
+ print_polyline_3d (f, node->polyline);
+ print_bsp_internal (f, (BSP_NODE *) node->in, n);
+ }
+}
+
+void
+print_bsp (FILE * f, BSP_NODE * bsp)
+{
+ print_bsp_internal (f, bsp, 0);
+}
+
+void
+add_polygon_to_sort (BSP_TREE * bsp, POLYGON_3D * polygon, void *attr)
+{
+ BSP_POLYGON_NODE *node = new_bsp_polygon_node (attr);
+ init_bsp_polygon_node (node, polygon);
+ node->next = *bsp;
+ *bsp = (BSP_NODE *) node;
+}
+
+void
+add_polyline_to_sort (BSP_TREE * bsp, POLYLINE_3D * polyline, void *attr)
+{
+ BSP_POLYLINE_NODE *node = new_bsp_polyline_node (attr);
+ init_bsp_polyline (node, polyline);
+ node->next = *bsp;
+ *bsp = (BSP_NODE *) node;
+}
+
+// quicksort of linked list
+#define FAR_DEPTH(Node) (-(Node)->extent->min[Z])
+#define NEAR_DEPTH(Node) (-(Node)->extent->max[Z])
+
+// leq provides ascending sort, so reverse to get max depth at head of list
+#define LEQ(A,B) (FAR_DEPTH(A) >= FAR_DEPTH(B))
+
+static void
+qs (BSP_NODE * hd, BSP_NODE * tl, BSP_NODE ** rtn)
+{
+ int nlo, nhi;
+ BSP_NODE *lo, *hi, *q, *p;
+
+ /* Invariant: Return head sorted with `tl' appended. */
+ while (hd != NULL)
+ {
+
+ nlo = nhi = 0;
+ lo = hi = NULL;
+ q = hd;
+ p = hd->next;
+
+ /* Reverse ascending prefix onto hi. This gives
+ O(n) behavior on sorted and reverse-sorted inputs. */
+ while (p != NULL && LEQ (p, hd))
+ {
+ hd->next = hi;
+ hi = hd;
+ ++nhi;
+ hd = p;
+ p = p->next;
+ }
+
+ /* If entire list was ascending, we're done. */
+ if (p == NULL)
+ {
+ *rtn = hd;
+ hd->next = hi;
+ q->next = tl;
+ return;
+ }
+
+ /* Partition and count sizes. */
+ while (p != NULL)
+ {
+ q = p->next;
+ if (LEQ (p, hd))
+ {
+ p->next = lo;
+ lo = p;
+ ++nlo;
+ }
+ else
+ {
+ p->next = hi;
+ hi = p;
+ ++nhi;
+ }
+ p = q;
+ }
+
+ /* Recur to establish invariant for sublists of hd,
+ choosing shortest list first to limit stack. */
+ if (nlo < nhi)
+ {
+ qs (lo, hd, rtn);
+ rtn = &hd->next;
+ hd = hi; /* Eliminated tail-recursive call. */
+ }
+ else
+ {
+ qs (hi, tl, &hd->next);
+ tl = hd;
+ hd = lo; /* Eliminated tail-recursive call. */
+ }
+ }
+ /* Base case of recurrence. Invariant is easy here. */
+ *rtn = tl;
+}
+
+static int
+xy_intersect_p (BOX_3D * a, BOX_3D * b)
+{
+ if (a->max[X] < b->min[X]) // a left of b
+ return 0;
+ if (a->min[X] > b->max[X]) // a right of b
+ return 0;
+ if (a->max[Y] < b->min[Y]) // a below b
+ return 0;
+ if (a->min[Y] > b->max[Y]) // a above b
+ return 0;
+ return 1;
+}
+
+#define SHORTEST_ALLOWABLE_SIDE 1e-4
+
+static void
+make_polygon_projection (POLYGON_2D * projection,
+ BSP_POLYGON_NODE * polygon_node)
+{
+ int i, j;
+
+ clear_polygon_2d (projection);
+ if (polygon_node->plane->n[Z] >= 0)
+ {
+ for (i = 0, j = polygon_node->polygon->n_sides - 1;
+ i < polygon_node->polygon->n_sides; j = i++)
+ {
+ if (dist_2d
+ (polygon_node->polygon->v[i],
+ polygon_node->polygon->v[j]) > SHORTEST_ALLOWABLE_SIDE)
+ copy_pt_2d (pushed_polygon_2d_v (projection),
+ polygon_node->polygon->v[i]);
+ }
+ }
+ else
+ {
+ for (i = polygon_node->polygon->n_sides - 1, j = 0; i >= 0; j = i--)
+ {
+ if (dist_2d
+ (polygon_node->polygon->v[i],
+ polygon_node->polygon->v[j]) > SHORTEST_ALLOWABLE_SIDE)
+ copy_pt_2d (pushed_polygon_2d_v (projection),
+ polygon_node->polygon->v[i]);
+ }
+ }
+}
+
+#define OVERLAP_EPS 1e-3
+
+int
+projections_overlap_p (BSP_POLYGON_NODE * p, BSP_POLYGON_NODE * q)
+{
+ int r;
+ POLYGON_2D p_projection[1], q_projection[1], cso[1];
+
+ init_polygon_2d (p_projection);
+ init_polygon_2d (q_projection);
+ init_polygon_2d (cso);
+
+ make_polygon_projection (p_projection, p);
+ make_polygon_projection (q_projection, q);
+ if (p_projection->n_sides < 2 && q_projection->n_sides < 2)
+ {
+ r = 0;
+ }
+ else if (p_projection->n_sides < 2)
+ {
+ r = point_near_convex_polygon_2d_p (p->polygon->v[0], q_projection,
+ OVERLAP_EPS);
+ }
+ else if (q_projection->n_sides < 2)
+ {
+ r = point_near_convex_polygon_2d_p (q->polygon->v[0], p_projection,
+ OVERLAP_EPS);
+ }
+ else
+ {
+ make_cso_polygon_2d (cso, p_projection, origin_2d, q_projection);
+ r = point_near_convex_polygon_2d_p (origin_2d, cso, OVERLAP_EPS);
+ }
+
+ clear_polygon_2d (p_projection);
+ clear_polygon_2d (q_projection);
+ clear_polygon_2d (cso);
+ return r;
+}
+
+#define OVERLAP_TOLERANCE 1e-3
+
+int
+polyline_projection_overlaps_polygon (BSP_POLYLINE_NODE * polyline_node,
+ BSP_POLYGON_NODE * polygon_node)
+{
+ int i, r;
+ POLYGON_2D polygon_projection[1], line_seg_projection[1], cso[1];
+
+ init_polygon_2d (polygon_projection);
+ init_polygon_2d (line_seg_projection);
+ init_polygon_2d (cso);
+
+ make_polygon_projection (polygon_projection, polygon_node);
+ if (polygon_projection->n_sides < 2)
+ {
+ // a point can't obscure a line
+ r = 0;
+ }
+ else if (polyline_node->polyline->n_vertices == 1)
+ {
+ // polyline is single vertex; just check to see if it lies in projection
+ r = point_near_convex_polygon_2d_p (polyline_node->polyline->v[0],
+ polygon_projection,
+ OVERLAP_TOLERANCE);
+ }
+ else
+ {
+
+ // use a two-sided polygon to model each edge;
+ setup_polygon_2d (line_seg_projection, 2);
+ line_seg_projection->n_sides = 2;
+ r = 0;
+ for (i = 0; i < polyline_node->polyline->n_vertices - 1; i++)
+ {
+ copy_pt_2d (line_seg_projection->v[0],
+ polyline_node->polyline->v[i]);
+ copy_pt_2d (line_seg_projection->v[1],
+ polyline_node->polyline->v[i + 1]);
+ make_cso_polygon_2d (cso, line_seg_projection, origin_2d,
+ polygon_projection);
+ r |= point_near_convex_polygon_2d_p (origin_2d, cso,
+ OVERLAP_TOLERANCE);
+ if (r)
+ break;
+ }
+ }
+ clear_polygon_2d (polygon_projection);
+ clear_polygon_2d (line_seg_projection);
+ clear_polygon_2d (cso);
+ return r;
+}
+
+static void
+debug_print (BSP_NODE * p)
+{
+ fprintf (stderr, "\nlist:\n");
+ while (p)
+ {
+ fprintf (stderr, " %p:%sprev=%p near=%.4g far=%.4g\n",
+ p,
+ p->mark ? "*" : " ", p->prev, NEAR_DEPTH (p), FAR_DEPTH (p));
+ p = p->next;
+ }
+}
+
+typedef struct make_list_of_bsp_env_t
+{
+ BSP_NODE *head, *tail;
+ int n;
+}
+MAKE_LIST_OF_BSP_ENV;
+
+static void
+make_list_of_bsp (BSP_NODE * bsp, void *v_env)
+{
+ MAKE_LIST_OF_BSP_ENV *env = (MAKE_LIST_OF_BSP_ENV *) v_env;
+ if (env->tail)
+ {
+ env->tail->next = bsp;
+ bsp->prev = env->tail;
+ bsp->next = NULL;
+ env->tail = bsp;
+ }
+ else
+ {
+ env->head = env->tail = bsp;
+ }
+ ++env->n;
+}
+
+// check invariants in the depth sort list
+static void
+check_consistency (BSP_TREE hd)
+{
+ int n_marks, n_other;
+ BSP_NODE *p, *q;
+
+ n_marks = 0;
+ for (q = NULL, p = hd; p && p->mark; q = p, p = p->next)
+ {
+ n_marks++;
+ if (p->prev != q)
+ {
+ debug_print (hd);
+ die (no_line, "broken prev pointer @ %d (%p != %p)", n_marks,
+ p->prev, q);
+ }
+ if (p->extent->min[X] == 0 && p->extent->max[X] == 0 &&
+ p->extent->min[Y] == 0 && p->extent->max[Y] == 0 &&
+ p->extent->min[Z] == 0 && p->extent->max[Z] == 0)
+ die (no_line, "unset extent");
+ }
+
+ n_other = 0;
+ for (; p; q = p, p = p->next)
+ {
+ n_other++;
+ if (p->mark)
+ die (no_line, "unexpected mark");
+ if (p->prev != q)
+ {
+ debug_print (hd);
+ die (no_line, "broken prev pointer @ %d (%p != %p)",
+ n_marks + n_other, p->prev, q);
+ }
+ if (p->extent->min[X] > p->extent->max[X])
+ die (no_line, "unset extent");
+ if (q && !q->mark && FAR_DEPTH (p) > FAR_DEPTH (q))
+ {
+ debug_print (hd);
+ die (no_line, "far depth out of order @ %d", n_marks + n_other);
+ }
+ }
+}
+
+void
+insert_by_depth (BSP_TREE * hd, BSP_NODE * node)
+{
+ BSP_NODE *p, *q;
+
+ // place p after insert point and q before
+ for (q = NULL, p = *hd;
+ p && (p->mark || FAR_DEPTH (p) > FAR_DEPTH (node)); q = p, p = p->next)
+ /* skip */ ;
+
+ // insert
+ node->prev = q;
+ node->next = p;
+ if (q)
+ q->next = node;
+ else
+ *hd = node;
+ if (p)
+ p->prev = node;
+}
+
+// this is taken almost directly from Newell's 1972 paper except that
+// a BSP is used to resolve intersections and cyclic overlaps and it
+// incorporates polyline objects
+void
+sort_by_depth (BSP_TREE * bsp)
+{
+ int side,
+ n_probes, n_swaps, n_nodes,
+ n_bsps, n_in, n_out, n_ppos, n_plos, n_bsp_in_nodes, n_bsp_out_nodes;
+ BSP_NODE *p, *p_next, *q, *prev, *t, *t_next, *r;
+ BSP_POLYGON_NODE *polygon_node;
+ BSP_POLYLINE_NODE *polyline_node;
+ PLANE *plane;
+ BSP_TREE sub_bsp;
+ MAKE_LIST_OF_BSP_ENV env[1];
+
+ // quicksort on deepest vertex, back-to-front
+ qs (*bsp, NULL, &p);
+
+ // hook up prev pointers and make sure marks are clear
+ n_nodes = 0;
+ for (prev = NULL, q = p; q; prev = q, q = q->next)
+ {
+ q->prev = prev;
+ q->mark = NULL;
+ ++n_nodes;
+ }
+
+ // keep some stats just for fun
+ n_probes = n_swaps = n_bsps = n_bsp_in_nodes = n_bsp_out_nodes =
+ n_ppos = n_plos = 0;
+
+ // debug_print(p);
+
+ // this is now output list
+ r = NULL;
+
+ // goto here whenever the current check fails
+ // for "p", the hopeful deepest polygon
+restart_overlap_check:
+
+ while (p)
+ {
+
+ if (n_nodes < 1000)
+ check_consistency (p);
+
+ // check overlapping objects for necessary swaps.
+ for (q = p->next;
+ q && (q->mark || FAR_DEPTH (q) > NEAR_DEPTH (p)); q = q->next)
+ {
+
+ ++n_probes;
+
+ // rectangular x-y extents don't overlap, so a moo point (utterly meaningless)
+ if (!xy_intersect_p (p->extent, q->extent))
+ continue;
+
+ // two polylines don't matter
+ // DEBUG: it really does if they're different colors
+ if (p->tag == BSP_POLYLINE && q->tag == BSP_POLYLINE)
+ continue;
+
+ // two polygons
+ if (p->tag == BSP_POLYGON && q->tag == BSP_POLYGON)
+ {
+
+ // p is contained wholly in the back halfspace of q
+ polygon_node = (BSP_POLYGON_NODE *) p;
+ plane = ((BSP_POLYGON_NODE *) q)->plane;
+ side = polygon_side_of_plane (polygon_node->polygon, plane);
+ if (side == S_ON ||
+ (plane->n[Z] >= 0 && side == S_IN) ||
+ (plane->n[Z] <= 0 && side == S_OUT))
+ continue;
+
+ // q is contained wholly in the front halfspace of p
+ polygon_node = (BSP_POLYGON_NODE *) q;
+ plane = ((BSP_POLYGON_NODE *) p)->plane;
+ side = polygon_side_of_plane (polygon_node->polygon, plane);
+ if (side == S_ON ||
+ (plane->n[Z] >= 0 && side == S_OUT) ||
+ (plane->n[Z] <= 0 && side == S_IN))
+ continue;
+
+ // projections do not overlap
+ ++n_ppos;
+ if (!projections_overlap_p
+ ((BSP_POLYGON_NODE *) p, (BSP_POLYGON_NODE *) q))
+ continue;
+ }
+
+ if (p->tag == BSP_POLYLINE)
+ { // and q is a polygon
+ polyline_node = (BSP_POLYLINE_NODE *) p;
+ plane = ((BSP_POLYGON_NODE *) q)->plane;
+ side = polyline_side_of_plane (polyline_node->polyline, plane);
+
+ // line is contained wholly in the back halfspace of plane
+ // lines lying on plane should be swapped so plane is drawn first
+ if ((plane->n[Z] >= 0 && side == S_IN) ||
+ (plane->n[Z] <= 0 && side == S_OUT))
+ continue;
+
+ // projections do not overlap
+ ++n_plos;
+ if (!polyline_projection_overlaps_polygon
+ (polyline_node, (BSP_POLYGON_NODE *) q))
+ continue;
+ }
+
+ if (q->tag == BSP_POLYLINE)
+ { // and p is a polygon
+ polyline_node = (BSP_POLYLINE_NODE *) q;
+ plane = ((BSP_POLYGON_NODE *) p)->plane;
+ side = polyline_side_of_plane (polyline_node->polyline, plane);
+
+ // line is on or contained wholly in the front halfspace of the plane
+ // a line lying on the plane can stay where it is
+ if (side == S_ON ||
+ (plane->n[Z] >= 0 && side == S_OUT) ||
+ (plane->n[Z] <= 0 && side == S_IN))
+ continue;
+
+ // projections do not overlap
+ ++n_plos;
+ if (!polyline_projection_overlaps_polygon
+ (polyline_node, (BSP_POLYGON_NODE *) p))
+ continue;
+ }
+
+ if (q->mark)
+ {
+
+ // we've discovered an intersection or cyclic overlap; break it by
+ // putting all the marked nodes in a bsp, then pulling them
+ // out and inserting them back on the list; remember our bsps
+ // need all polygons inserted before all polylines
+
+ ++n_bsps;
+ sub_bsp = NULL;
+ n_in = 0;
+ t = NULL; // use t to hold polylines temporarily
+ while (p && p->mark)
+ {
+ p_next = p->next;
+
+ if (p->tag == BSP_POLYGON)
+ {
+ p->next = p->prev = NULL;
+ insert_polygon (&sub_bsp, (BSP_POLYGON_NODE *) p);
+ }
+ else
+ { // save polyline temporarily
+ p->next = t;
+ t = p;
+ }
+ ++n_in;
+ p = p_next;
+ if (p)
+ p->prev = NULL;
+ }
+ // insert the polylines now that all polygons are complete
+ while (t)
+ {
+ t_next = t->next;
+ t->next = t->prev = NULL;
+ insert_polyline (&sub_bsp, (BSP_POLYLINE_NODE *) t);
+ t = t_next;
+ }
+
+ // now traverse the bsp to get the objects back out, including split ones
+ env->n = 0;
+ env->head = env->tail = NULL;
+ traverse_bsp (sub_bsp, make_list_of_bsp, env);
+ n_out = env->n;
+
+ // splitting should always increase the number of primitives, but
+ // polygons very close in depth can cause split to fail; just shovel
+ // the result polygons to the output with a warning.
+ if (n_out <= n_in)
+ {
+ warn (no_line, "split failed in=%d, out=%d", n_in, n_out);
+ remark (no_line,
+ "if hidden surfaces are wrong, try -b option");
+ for (t = env->head; t; t = t_next)
+ {
+ t_next = t->next;
+ t->next = r;
+ t->in = t->out = t->on = t->prev = t->mark = NULL;
+ r = t;
+ }
+ goto restart_overlap_check;
+ }
+ // re-insert in the sorted list
+ for (t = env->head; t; t = t_next)
+ {
+ t_next = t->next;
+ t->in = t->out = t->on = t->prev = t->next = t->mark = NULL;
+ insert_by_depth (&p, t);
+ }
+
+ n_bsp_in_nodes += n_in;
+ n_bsp_out_nodes += n_out;
+
+ goto restart_overlap_check;
+ }
+ else
+ {
+
+ // no overlap, so pull q forward to head of list
+
+ // unlink q
+ if (q->next)
+ q->next->prev = q->prev;
+ q->prev->next = q->next;
+
+ // mark
+ q->mark = p;
+
+ // push
+ q->prev = NULL;
+
+ q->next = p;
+ p->prev = q;
+ p = q;
+
+ ++n_swaps;
+
+ goto restart_overlap_check;
+ }
+ }
+
+ // overlap check saw no problems; pop head onto return list
+ p_next = p->next;
+ if (p_next)
+ p_next->prev = NULL;
+
+ // push on output
+ p->next = r;
+ p->prev = NULL;
+ r = p;
+
+ // move to next item
+ p = p_next;
+ }
+ // pop from q and push onto q until q is empty
+ q = r;
+ r = NULL;
+ while (q)
+ {
+ t = q;
+ q = q->next; // pop
+ t->next = r;
+ r = t; // push
+ }
+
+ {
+ int n_probes_possible = n_nodes + n_bsp_out_nodes - n_bsp_in_nodes;
+
+ remark (no_line,
+ "node=%d probe=%.1lf swap=%d split=%d (in=%d out=%d) ols=%d/%d",
+ n_nodes,
+ (n_probes_possible >=
+ 0) ? (double) n_probes / n_probes_possible : 0.0, n_swaps,
+ n_bsps, n_bsp_in_nodes, n_bsp_out_nodes, n_ppos, n_plos);
+ }
+
+ *bsp = r;
+}