initial setup of thesis repository

[cluster_expansion_thesis.git] / little_helpers / tikz / sketch-0.2.161 / bsp.c

/* 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;
}