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