/* bsp.c
   Copyright (C) 2005,2006,2007,2008 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;
}