/******************************************************************************
 *
 * 9998sketch/bsp.c
 *
 * (c) Eugene K. Ressler, Jr. 2004, 2005
 *
 * A program related to the book
 *
 * Fundamentals of Virtual World Simulation
 * by Eugene K. Ressler, Jr.
 *
 * No warranty of correctness or capability of any kind is expressed or implied.
 *
 * The author grants free use of this code for any purpose as long any text,
 * executable program, library file, or source code distributed to others 
 * and employing any portion of this code clearly cites the book named above.
 *
 ******************************************************************************/

#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;
}