Logo Search packages:      
Sourcecode: s3d version File versions  Download package

sei_construct.c

/*
 * sei_construct.c
 *
 * Copyright (C) 2006-2008 Simon Wunderlich <dotslash@packetmixer.de>
 * Copyright (C) 1994 A. Narkhede and D .Manocha, who released their code
 * for public domain:
 * <snip>
 *
 * This code is in the public domain. Specifically, we give to the public
 * domain all rights for future licensing of the source code, all resale
 * rights, and all publishing rights.
 *
 * UNC-CH GIVES NO WARRANTY, EXPRESSED OR IMPLIED, FOR THE SOFTWARE
 * AND/OR DOCUMENTATION PROVIDED, INCLUDING, WITHOUT LIMITATION, WARRANTY
 * OF MERCHANTABILITY AND WARRANTY OF FITNESS FOR A PARTICULAR PURPOSE.
 *
 *
 *                                 - Atul Narkhede (narkhede@cs.unc.edu)
 * </snip>
 *
 * This file is part of the s3d API, the API of s3d (the 3d network display server).
 * See http://s3d.berlios.de/ for more updates.
 *
 * The s3d API is free software; you can redistribute it and/or modify
 * it under the terms of the GNU Lesser General Public License as published by
 * the Free Software Foundation; either version 2.1 of the License, or
 * (at your option) any later version.
 *
 * The s3d API 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 Lesser General Public License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public License
 * along with the s3d API; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
 */


#include "sei_triangulate.h"
#include <math.h>
#include <string.h> /* memset() */


node_t qs[QSIZE];  /* Query structure */
trap_t tr[TRSIZE];  /* Trapezoid structure */
segment_t seg[SEGSIZE];  /* Segment table */

static int q_idx;
static int tr_idx;

/* Return a new node to be added into the query tree */
static int newnode(void)
{
      if (q_idx < QSIZE)
            return q_idx++;
      else {
            errs("sei:newnode()", "Query-table overflow");
            return -1;
      }
}

/* Return a free trapezoid */
static int newtrap(void)
{
      if (tr_idx < TRSIZE) {
            tr[tr_idx].lseg = -1;
            tr[tr_idx].rseg = -1;
            tr[tr_idx].state = ST_VALID;
            return tr_idx++;
      } else {
            errs("sei:newtrap()", "Trapezoid-table overflow");
            return -1;
      }
}


/* Return the maximum of the two points into the yval structure */
static int _max(point_t *yval, point_t *v0, point_t *v1)
{
      if (v0->y > v1->y + C_EPS)
            *yval = *v0;
      else if (FP_EQUAL(v0->y, v1->y)) {
            if (v0->x > v1->x + C_EPS)
                  *yval = *v0;
            else
                  *yval = *v1;
      } else
            *yval = *v1;

      return 0;
}


/* Return the minimum of the two points into the yval structure */
static int _min(point_t *yval, point_t *v0, point_t *v1)
{
      if (v0->y < v1->y - C_EPS)
            *yval = *v0;
      else if (FP_EQUAL(v0->y, v1->y)) {
            if (v0->x < v1->x)
                  *yval = *v0;
            else
                  *yval = *v1;
      } else
            *yval = *v1;

      return 0;
}


int _greater_than(point_t *v0, point_t *v1)
{
      if (v0->y > v1->y + C_EPS)
            return TRUE;
      else if (v0->y < v1->y - C_EPS)
            return FALSE;
      else
            return (v0->x > v1->x);
}


int _equal_to(point_t *v0, point_t *v1)
{
      return (FP_EQUAL(v0->y, v1->y) && FP_EQUAL(v0->x, v1->x));
}

int _greater_than_equal_to(point_t *v0, point_t *v1)
{
      if (v0->y > v1->y + C_EPS)
            return TRUE;
      else if (v0->y < v1->y - C_EPS)
            return FALSE;
      else
            return (v0->x >= v1->x);
}

int _less_than(point_t *v0, point_t *v1)
{
      if (v0->y < v1->y - C_EPS)
            return TRUE;
      else if (v0->y > v1->y + C_EPS)
            return FALSE;
      else
            return (v0->x < v1->x);
}


/* Initilialise the query structure (Q) and the trapezoid table (T)
 * when the first segment is added to start the trapezoidation. The
 * query-tree starts out with 4 trapezoids, one S-node and 2 Y-nodes
 *
 *                4
 *   -----------------------------------
 *      \
 *   1    \        2
 *        \
 *   -----------------------------------
 *                3
 */

static int init_query_structure(int segnum)
{
      int i1, i2, i3, i4, i5, i6, i7, root;
      int t1, t2, t3, t4;
      segment_t *s = &seg[segnum];

      q_idx = tr_idx = 1;
      memset((void *)tr, 0, sizeof(tr));
      memset((void *)qs, 0, sizeof(qs));

      i1 = newnode();
      qs[i1].nodetype = T_Y;
      _max(&qs[i1].yval, &s->v0, &s->v1); /* root */
      root = i1;

      qs[i1].right = i2 = newnode();
      qs[i2].nodetype = T_SINK;
      qs[i2].parent = i1;

      qs[i1].left = i3 = newnode();
      qs[i3].nodetype = T_Y;
      _min(&qs[i3].yval, &s->v0, &s->v1); /* root */
      qs[i3].parent = i1;

      qs[i3].left = i4 = newnode();
      qs[i4].nodetype = T_SINK;
      qs[i4].parent = i3;

      qs[i3].right = i5 = newnode();
      qs[i5].nodetype = T_X;
      qs[i5].segnum = segnum;
      qs[i5].parent = i3;

      qs[i5].left = i6 = newnode();
      qs[i6].nodetype = T_SINK;
      qs[i6].parent = i5;

      qs[i5].right = i7 = newnode();
      qs[i7].nodetype = T_SINK;
      qs[i7].parent = i5;

      t1 = newtrap();  /* middle left */
      t2 = newtrap();  /* middle right */
      t3 = newtrap();  /* bottom-most */
      t4 = newtrap();  /* topmost */

      tr[t1].hi = tr[t2].hi = tr[t4].lo = qs[i1].yval;
      tr[t1].lo = tr[t2].lo = tr[t3].hi = qs[i3].yval;
      tr[t4].hi.y = (double)(INFINITY);
      tr[t4].hi.x = (double)(INFINITY);
      tr[t3].lo.y = (double) - 1 * (INFINITY);
      tr[t3].lo.x = (double) - 1 * (INFINITY);
      tr[t1].rseg = tr[t2].lseg = segnum;
      tr[t1].u0 = tr[t2].u0 = t4;
      tr[t1].d0 = tr[t2].d0 = t3;
      tr[t4].d0 = tr[t3].u0 = t1;
      tr[t4].d1 = tr[t3].u1 = t2;

      tr[t1].sink = i6;
      tr[t2].sink = i7;
      tr[t3].sink = i4;
      tr[t4].sink = i2;

      tr[t1].state = tr[t2].state = ST_VALID;
      tr[t3].state = tr[t4].state = ST_VALID;

      qs[i2].trnum = t4;
      qs[i4].trnum = t3;
      qs[i6].trnum = t1;
      qs[i7].trnum = t2;

      s->is_inserted = TRUE;
      return root;
}


/* Retun TRUE if the vertex v is to the left of line segment no.
 * segnum. Takes care of the degenerate cases when both the vertices
 * have the same y--cood, etc.
 */

static int is_left_of(int segnum, point_t *v)
{
      segment_t *s = &seg[segnum];
      double area;

      if (_greater_than(&s->v1, &s->v0)) { /* seg. going upwards */
            if (FP_EQUAL(s->v1.y, v->y)) {
                  if (v->x < s->v1.x)
                        area = 1.0;
                  else
                        area = -1.0;
            } else if (FP_EQUAL(s->v0.y, v->y)) {
                  if (v->x < s->v0.x)
                        area = 1.0;
                  else
                        area = -1.0;
            } else
                  area = CROSS(s->v0, s->v1, (*v));
      } else {   /* v0 > v1 */
            if (FP_EQUAL(s->v1.y, v->y)) {
                  if (v->x < s->v1.x)
                        area = 1.0;
                  else
                        area = -1.0;
            } else if (FP_EQUAL(s->v0.y, v->y)) {
                  if (v->x < s->v0.x)
                        area = 1.0;
                  else
                        area = -1.0;
            } else
                  area = CROSS(s->v1, s->v0, (*v));
      }

      if (area > 0.0)
            return TRUE;
      else
            return FALSE;
}



/* Returns true if the corresponding endpoint of the given segment is */
/* already inserted into the segment tree. Use the simple test of */
/* whether the segment which shares this endpoint is already inserted */

static int inserted(int segnum, int whichpt)
{
      if (whichpt == FIRSTPT)
            return seg[seg[segnum].prev].is_inserted;
      else
            return seg[seg[segnum].next].is_inserted;
}

/* This is query routine which determines which trapezoid does the
 * point v lie in. The return value is the trapezoid number.
 */

int locate_endpoint(point_t *v, point_t *vo, int r)
{
      node_t *rptr = &qs[r];

      switch (rptr->nodetype) {
      case T_SINK:
            return rptr->trnum;

      case T_Y:
            if (_greater_than(v, &rptr->yval)) /* above */
                  return locate_endpoint(v, vo, rptr->right);
            else if (_equal_to(v, &rptr->yval)) /* the point is already */
            {             /* inserted. */
                  if (_greater_than(vo, &rptr->yval)) /* above */
                        return locate_endpoint(v, vo, rptr->right);
                  else
                        return locate_endpoint(v, vo, rptr->left); /* below */
            } else
                  return locate_endpoint(v, vo, rptr->left); /* below */

      case T_X:
            if (_equal_to(v, &seg[rptr->segnum].v0) ||
                            _equal_to(v, &seg[rptr->segnum].v1)) {
                  if (FP_EQUAL(v->y, vo->y)) { /* horizontal segment */
                        if (vo->x < v->x)
                              return locate_endpoint(v, vo, rptr->left); /* left */
                        else
                              return locate_endpoint(v, vo, rptr->right); /* right */
                  }

                  else if (is_left_of(rptr->segnum, vo))
                        return locate_endpoint(v, vo, rptr->left); /* left */
                  else
                        return locate_endpoint(v, vo, rptr->right); /* right */
            } else if (is_left_of(rptr->segnum, v))
                  return locate_endpoint(v, vo, rptr->left); /* left */
            else
                  return locate_endpoint(v, vo, rptr->right); /* right */

      default:
            errs("sei:locate_endpoint()", "Haggu!!!! (whatever)");
            break;
      }
      return(-1);
}


/* Thread in the segment into the existing trapezoidation. The
 * limiting trapezoids are given by tfirst and tlast (which are the
 * trapezoids containing the two endpoints of the segment. Merges all
 * possible trapezoids which flank this segment and have been recently
 * divided because of its insertion
 */

static int merge_trapezoids(int segnum, int tfirst, int tlast, int side)
{
      int t, tnext, cond;
      int ptnext;

      /* First merge polys on the LHS */
      t = tfirst;
      while ((t > 0) && _greater_than_equal_to(&tr[t].lo, &tr[tlast].lo)) {
            if (side == S_LEFT)
                  cond = ((((tnext = tr[t].d0) > 0) && (tr[tnext].rseg == segnum)) ||
                          (((tnext = tr[t].d1) > 0) && (tr[tnext].rseg == segnum)));
            else
                  cond = ((((tnext = tr[t].d0) > 0) && (tr[tnext].lseg == segnum)) ||
                          (((tnext = tr[t].d1) > 0) && (tr[tnext].lseg == segnum)));

            if (cond) {
                  if ((tr[t].lseg == tr[tnext].lseg) &&
                                  (tr[t].rseg == tr[tnext].rseg)) /* good neighbours */
                  {                 /* merge them */
                        /* Use the upper node as the new node i.e. t */

                        ptnext = qs[tr[tnext].sink].parent;

                        if (qs[ptnext].left == tr[tnext].sink)
                              qs[ptnext].left = tr[t].sink;
                        else
                              qs[ptnext].right = tr[t].sink; /* redirect parent */


                        /* Change the upper neighbours of the lower trapezoids */

                        if ((tr[t].d0 = tr[tnext].d0) > 0) {
                              if (tr[tr[t].d0].u0 == tnext)
                                    tr[tr[t].d0].u0 = t;
                              else if (tr[tr[t].d0].u1 == tnext)
                                    tr[tr[t].d0].u1 = t;
                        }

                        if ((tr[t].d1 = tr[tnext].d1) > 0) {
                              if (tr[tr[t].d1].u0 == tnext)
                                    tr[tr[t].d1].u0 = t;
                              else if (tr[tr[t].d1].u1 == tnext)
                                    tr[tr[t].d1].u1 = t;
                        }

                        tr[t].lo = tr[tnext].lo;
                        tr[tnext].state = ST_INVALID; /* invalidate the lower */
                        /* trapezium */
                  } else     /* not good neighbours */
                        t = tnext;
            } else     /* do not satisfy the outer if */
                  t = tnext;

      } /* end-while */

      return 0;
}


/* Add in the new segment into the trapezoidation and update Q and T
 * structures. First locate the two endpoints of the segment in the
 * Q-structure. Then start from the topmost trapezoid and go down to
 * the  lower trapezoid dividing all the trapezoids in between .
 */

static int add_segment(int segnum)
{
      segment_t s;
      int tu, tl, sk, tfirst, tlast;
      int tfirstr = 0, tlastr = 0, tfirstl, tlastl;
      int i1, i2, t, tn;
      point_t tpt;
      int tribot = 0, is_swapped = 0;
      int tmptriseg;
      int tmpseg = 1;

      s = seg[segnum];
      if (_greater_than(&s.v1, &s.v0)) { /* Get higher vertex in v0 */
            int tmp;
            tpt = s.v0;
            s.v0 = s.v1;
            s.v1 = tpt;
            tmp = s.root0;
            s.root0 = s.root1;
            s.root1 = tmp;
            is_swapped = TRUE;
      }

      if ((is_swapped) ? !inserted(segnum, LASTPT) :
                      !inserted(segnum, FIRSTPT)) {   /* insert v0 in the tree */
            int tmp_d;

            tu = locate_endpoint(&s.v0, &s.v1, s.root0);
            tl = newtrap();  /* tl is the new lower trapezoid */
            tr[tl].state = ST_VALID;
            tr[tl] = tr[tu];
            tr[tu].lo.y = tr[tl].hi.y = s.v0.y;
            tr[tu].lo.x = tr[tl].hi.x = s.v0.x;
            tr[tu].d0 = tl;
            tr[tu].d1 = 0;
            tr[tl].u0 = tu;
            tr[tl].u1 = 0;

            if (((tmp_d = tr[tl].d0) > 0) && (tr[tmp_d].u0 == tu))
                  tr[tmp_d].u0 = tl;
            if (((tmp_d = tr[tl].d0) > 0) && (tr[tmp_d].u1 == tu))
                  tr[tmp_d].u1 = tl;

            if (((tmp_d = tr[tl].d1) > 0) && (tr[tmp_d].u0 == tu))
                  tr[tmp_d].u0 = tl;
            if (((tmp_d = tr[tl].d1) > 0) && (tr[tmp_d].u1 == tu))
                  tr[tmp_d].u1 = tl;

            /* Now update the query structure and obtain the sinks for the */
            /* two trapezoids */

            i1 = newnode();  /* Upper trapezoid sink */
            i2 = newnode();  /* Lower trapezoid sink */
            sk = tr[tu].sink;

            qs[sk].nodetype = T_Y;
            qs[sk].yval = s.v0;
            qs[sk].segnum = segnum; /* not really reqd ... maybe later */
            qs[sk].left = i2;
            qs[sk].right = i1;

            qs[i1].nodetype = T_SINK;
            qs[i1].trnum = tu;
            qs[i1].parent = sk;

            qs[i2].nodetype = T_SINK;
            qs[i2].trnum = tl;
            qs[i2].parent = sk;

            tr[tu].sink = i1;
            tr[tl].sink = i2;
            tfirst = tl;
      } else {    /* v0 already present */
            /* Get the topmost intersecting trapezoid */
            tfirst = locate_endpoint(&s.v0, &s.v1, s.root0);
      }


      if ((is_swapped) ? !inserted(segnum, FIRSTPT) :
                      !inserted(segnum, LASTPT)) {   /* insert v1 in the tree */
            int tmp_d;

            tu = locate_endpoint(&s.v1, &s.v0, s.root1);

            tl = newtrap();  /* tl is the new lower trapezoid */
            tr[tl].state = ST_VALID;
            tr[tl] = tr[tu];
            tr[tu].lo.y = tr[tl].hi.y = s.v1.y;
            tr[tu].lo.x = tr[tl].hi.x = s.v1.x;
            tr[tu].d0 = tl;
            tr[tu].d1 = 0;
            tr[tl].u0 = tu;
            tr[tl].u1 = 0;

            if (((tmp_d = tr[tl].d0) > 0) && (tr[tmp_d].u0 == tu))
                  tr[tmp_d].u0 = tl;
            if (((tmp_d = tr[tl].d0) > 0) && (tr[tmp_d].u1 == tu))
                  tr[tmp_d].u1 = tl;

            if (((tmp_d = tr[tl].d1) > 0) && (tr[tmp_d].u0 == tu))
                  tr[tmp_d].u0 = tl;
            if (((tmp_d = tr[tl].d1) > 0) && (tr[tmp_d].u1 == tu))
                  tr[tmp_d].u1 = tl;

            /* Now update the query structure and obtain the sinks for the */
            /* two trapezoids */

            i1 = newnode();  /* Upper trapezoid sink */
            i2 = newnode();  /* Lower trapezoid sink */
            sk = tr[tu].sink;

            qs[sk].nodetype = T_Y;
            qs[sk].yval = s.v1;
            qs[sk].segnum = segnum; /* not really reqd ... maybe later */
            qs[sk].left = i2;
            qs[sk].right = i1;

            qs[i1].nodetype = T_SINK;
            qs[i1].trnum = tu;
            qs[i1].parent = sk;

            qs[i2].nodetype = T_SINK;
            qs[i2].trnum = tl;
            qs[i2].parent = sk;

            tr[tu].sink = i1;
            tr[tl].sink = i2;
            tlast = tu;
      } else    /* v1 already present */
      {       /* Get the lowermost intersecting trapezoid */
            tlast = locate_endpoint(&s.v1, &s.v0, s.root1);
            tribot = 1;
      }

      /* Thread the segment into the query tree creating a new X-node */
      /* First, split all the trapezoids which are intersected by s into */
      /* two */

      t = tfirst;   /* topmost trapezoid */

      while ((t > 0) &&
                      _greater_than_equal_to(&tr[t].lo, &tr[tlast].lo))
            /* traverse from top to bot */
      {
            int t_sav, tn_sav;
            sk = tr[t].sink;
            i1 = newnode();  /* left trapezoid sink */
            i2 = newnode();  /* right trapezoid sink */

            qs[sk].nodetype = T_X;
            qs[sk].segnum = segnum;
            qs[sk].left = i1;
            qs[sk].right = i2;

            qs[i1].nodetype = T_SINK; /* left trapezoid (use existing one) */
            qs[i1].trnum = t;
            qs[i1].parent = sk;

            qs[i2].nodetype = T_SINK; /* right trapezoid (allocate new) */
            qs[i2].trnum = tn = newtrap();
            tr[tn].state = ST_VALID;
            qs[i2].parent = sk;

            if (t == tfirst)
                  tfirstr = tn;
            if (_equal_to(&tr[t].lo, &tr[tlast].lo))
                  tlastr = tn;

            tr[tn] = tr[t];
            tr[t].sink = i1;
            tr[tn].sink = i2;
            t_sav = t;
            tn_sav = tn;

            /* error */

            if ((tr[t].d0 <= 0) && (tr[t].d1 <= 0)) { /* case cannot arise */
                  errs("sei:add_segment()", "error");
                  break;
            }

            /* only one trapezoid below. partition t into two and make the */
            /* two resulting trapezoids t and tn as the upper neighbours of */
            /* the sole lower trapezoid */

            else if ((tr[t].d0 > 0) && (tr[t].d1 <= 0)) {   /* Only one trapezoid below */
                  if ((tr[t].u0 > 0) && (tr[t].u1 > 0)) {   /* continuation of a chain from abv. */
                        if (tr[t].usave > 0) { /* three upper neighbours */
                              if (tr[t].uside == S_LEFT) {
                                    tr[tn].u0 = tr[t].u1;
                                    tr[t].u1 = -1;
                                    tr[tn].u1 = tr[t].usave;

                                    tr[tr[t].u0].d0 = t;
                                    tr[tr[tn].u0].d0 = tn;
                                    tr[tr[tn].u1].d0 = tn;
                              } else { /* intersects in the right */
                                    tr[tn].u1 = -1;
                                    tr[tn].u0 = tr[t].u1;
                                    tr[t].u1 = tr[t].u0;
                                    tr[t].u0 = tr[t].usave;

                                    tr[tr[t].u0].d0 = t;
                                    tr[tr[t].u1].d0 = t;
                                    tr[tr[tn].u0].d0 = tn;
                              }

                              tr[t].usave = tr[tn].usave = 0;
                        } else { /* No usave.... simple case */
                              tr[tn].u0 = tr[t].u1;
                              tr[t].u1 = tr[tn].u1 = -1;
                              tr[tr[tn].u0].d0 = tn;
                        }
                  } else {   /* fresh seg. or upward cusp */
                        int tmp_u = tr[t].u0;
                        int td0;
                        if (((td0 = tr[tmp_u].d0) > 0) &&
                                        (tr[tmp_u].d1 > 0)) {  /* upward cusp */
                              if ((tr[td0].rseg > 0) &&
                                              !is_left_of(tr[td0].rseg, &s.v1)) {
                                    tr[t].u0 = tr[t].u1 = tr[tn].u1 = -1;
                                    tr[tr[tn].u0].d1 = tn;
                              } else { /* cusp going leftwards */
                                    tr[tn].u0 = tr[tn].u1 = tr[t].u1 = -1;
                                    tr[tr[t].u0].d0 = t;
                              }
                        } else { /* fresh segment */
                              tr[tr[t].u0].d0 = t;
                              tr[tr[t].u0].d1 = tn;
                        }
                  }

                  if (FP_EQUAL(tr[t].lo.y, tr[tlast].lo.y) &&
                                  FP_EQUAL(tr[t].lo.x, tr[tlast].lo.x) && tribot) {  /* bottom forms a triangle */

                        if (is_swapped)
                              tmptriseg = seg[segnum].prev;
                        else
                              tmptriseg = seg[segnum].next;

                        if ((tmptriseg > 0) && is_left_of(tmptriseg, &s.v0)) {
                              /* L-R downward cusp */
                              tr[tr[t].d0].u0 = t;
                              tr[tn].d0 = tr[tn].d1 = -1;
                        } else {
                              /* R-L downward cusp */
                              tr[tr[tn].d0].u1 = tn;
                              tr[t].d0 = tr[t].d1 = -1;
                        }
                  } else {
                        if ((tr[tr[t].d0].u0 > 0) && (tr[tr[t].d0].u1 > 0)) {
                              if (tr[tr[t].d0].u0 == t) { /* passes thru LHS */
                                    tr[tr[t].d0].usave = tr[tr[t].d0].u1;
                                    tr[tr[t].d0].uside = S_LEFT;
                              } else {
                                    tr[tr[t].d0].usave = tr[tr[t].d0].u0;
                                    tr[tr[t].d0].uside = S_RIGHT;
                              }
                        }
                        tr[tr[t].d0].u0 = t;
                        tr[tr[t].d0].u1 = tn;
                  }

                  t = tr[t].d0;
            }


            else if ((tr[t].d0 <= 0) && (tr[t].d1 > 0)) {   /* Only one trapezoid below */
                  if ((tr[t].u0 > 0) && (tr[t].u1 > 0)) {   /* continuation of a chain from abv. */
                        if (tr[t].usave > 0) { /* three upper neighbours */
                              if (tr[t].uside == S_LEFT) {
                                    tr[tn].u0 = tr[t].u1;
                                    tr[t].u1 = -1;
                                    tr[tn].u1 = tr[t].usave;

                                    tr[tr[t].u0].d0 = t;
                                    tr[tr[tn].u0].d0 = tn;
                                    tr[tr[tn].u1].d0 = tn;
                              } else { /* intersects in the right */
                                    tr[tn].u1 = -1;
                                    tr[tn].u0 = tr[t].u1;
                                    tr[t].u1 = tr[t].u0;
                                    tr[t].u0 = tr[t].usave;

                                    tr[tr[t].u0].d0 = t;
                                    tr[tr[t].u1].d0 = t;
                                    tr[tr[tn].u0].d0 = tn;
                              }

                              tr[t].usave = tr[tn].usave = 0;
                        } else { /* No usave.... simple case */
                              tr[tn].u0 = tr[t].u1;
                              tr[t].u1 = tr[tn].u1 = -1;
                              tr[tr[tn].u0].d0 = tn;
                        }
                  } else {   /* fresh seg. or upward cusp */
                        int tmp_u = tr[t].u0;
                        int td0;
                        if (((td0 = tr[tmp_u].d0) > 0) &&
                                        (tr[tmp_u].d1 > 0)) {  /* upward cusp */
                              if ((tr[td0].rseg > 0) &&
                                              !is_left_of(tr[td0].rseg, &s.v1)) {
                                    tr[t].u0 = tr[t].u1 = tr[tn].u1 = -1;
                                    tr[tr[tn].u0].d1 = tn;
                              } else {
                                    tr[tn].u0 = tr[tn].u1 = tr[t].u1 = -1;
                                    tr[tr[t].u0].d0 = t;
                              }
                        } else { /* fresh segment */
                              tr[tr[t].u0].d0 = t;
                              tr[tr[t].u0].d1 = tn;
                        }
                  }

                  if (FP_EQUAL(tr[t].lo.y, tr[tlast].lo.y) &&
                                  FP_EQUAL(tr[t].lo.x, tr[tlast].lo.x) && tribot) {  /* bottom forms a triangle */

                        if ((tmpseg > 0) && is_left_of(tmpseg, &s.v0)) {
                              /* L-R downward cusp */
                              tr[tr[t].d1].u0 = t;
                              tr[tn].d0 = tr[tn].d1 = -1;
                        } else {
                              /* R-L downward cusp */
                              tr[tr[tn].d1].u1 = tn;
                              tr[t].d0 = tr[t].d1 = -1;
                        }
                  } else {
                        if ((tr[tr[t].d1].u0 > 0) && (tr[tr[t].d1].u1 > 0)) {
                              if (tr[tr[t].d1].u0 == t) { /* passes thru LHS */
                                    tr[tr[t].d1].usave = tr[tr[t].d1].u1;
                                    tr[tr[t].d1].uside = S_LEFT;
                              } else {
                                    tr[tr[t].d1].usave = tr[tr[t].d1].u0;
                                    tr[tr[t].d1].uside = S_RIGHT;
                              }
                        }
                        tr[tr[t].d1].u0 = t;
                        tr[tr[t].d1].u1 = tn;
                  }

                  t = tr[t].d1;
            }

            /* two trapezoids below. Find out which one is intersected by */
            /* this segment and proceed down that one */

            else {
                  double y0, yt;
                  point_t tmppt;
                  int tnext, i_d0;
                  tmpseg = tr[tr[t].d0].rseg;

                  i_d0 = FALSE;
                  if (FP_EQUAL(tr[t].lo.y, s.v0.y)) {
                        if (tr[t].lo.x > s.v0.x)
                              i_d0 = TRUE;
                  } else {
                        tmppt.y = y0 = tr[t].lo.y;
                        yt = (y0 - s.v0.y) / (s.v1.y - s.v0.y);
                        tmppt.x = s.v0.x + yt * (s.v1.x - s.v0.x);

                        if (_less_than(&tmppt, &tr[t].lo))
                              i_d0 = TRUE;
                  }

                  /* check continuity from the top so that the lower-neighbour */
                  /* values are properly filled for the upper trapezoid */

                  if ((tr[t].u0 > 0) && (tr[t].u1 > 0)) {   /* continuation of a chain from abv. */
                        if (tr[t].usave > 0) { /* three upper neighbours */
                              if (tr[t].uside == S_LEFT) {
                                    tr[tn].u0 = tr[t].u1;
                                    tr[t].u1 = -1;
                                    tr[tn].u1 = tr[t].usave;

                                    tr[tr[t].u0].d0 = t;
                                    tr[tr[tn].u0].d0 = tn;
                                    tr[tr[tn].u1].d0 = tn;
                              } else { /* intersects in the right */
                                    tr[tn].u1 = -1;
                                    tr[tn].u0 = tr[t].u1;
                                    tr[t].u1 = tr[t].u0;
                                    tr[t].u0 = tr[t].usave;

                                    tr[tr[t].u0].d0 = t;
                                    tr[tr[t].u1].d0 = t;
                                    tr[tr[tn].u0].d0 = tn;
                              }

                              tr[t].usave = tr[tn].usave = 0;
                        } else { /* No usave.... simple case */
                              tr[tn].u0 = tr[t].u1;
                              tr[tn].u1 = -1;
                              tr[t].u1 = -1;
                              tr[tr[tn].u0].d0 = tn;
                        }
                  } else {   /* fresh seg. or upward cusp */
                        int tmp_u = tr[t].u0;
                        int td0;
                        if (((td0 = tr[tmp_u].d0) > 0) &&
                                        (tr[tmp_u].d1 > 0)) {  /* upward cusp */
                              if ((tr[td0].rseg > 0) &&
                                              !is_left_of(tr[td0].rseg, &s.v1)) {
                                    tr[t].u0 = tr[t].u1 = tr[tn].u1 = -1;
                                    tr[tr[tn].u0].d1 = tn;
                              } else {
                                    tr[tn].u0 = tr[tn].u1 = tr[t].u1 = -1;
                                    tr[tr[t].u0].d0 = t;
                              }
                        } else { /* fresh segment */
                              tr[tr[t].u0].d0 = t;
                              tr[tr[t].u0].d1 = tn;
                        }
                  }

                  if (FP_EQUAL(tr[t].lo.y, tr[tlast].lo.y) &&
                                  FP_EQUAL(tr[t].lo.x, tr[tlast].lo.x) && tribot) {
                        /* this case arises only at the lowest trapezoid.. i.e.
                           tlast, if the lower endpoint of the segment is
                           already inserted in the structure */

                        tr[tr[t].d0].u0 = t;
                        tr[tr[t].d0].u1 = -1;
                        tr[tr[t].d1].u0 = tn;
                        tr[tr[t].d1].u1 = -1;

                        tr[tn].d0 = tr[t].d1;
                        tr[t].d1 = tr[tn].d1 = -1;

                        tnext = tr[t].d1;
                  } else if (i_d0)
                        /* intersecting d0 */
                  {
                        tr[tr[t].d0].u0 = t;
                        tr[tr[t].d0].u1 = tn;
                        tr[tr[t].d1].u0 = tn;
                        tr[tr[t].d1].u1 = -1;

                        /* new code to determine the bottom neighbours of the */
                        /* newly partitioned trapezoid */

                        tr[t].d1 = -1;

                        tnext = tr[t].d0;
                  } else {  /* intersecting d1 */
                        tr[tr[t].d0].u0 = t;
                        tr[tr[t].d0].u1 = -1;
                        tr[tr[t].d1].u0 = t;
                        tr[tr[t].d1].u1 = tn;

                        /* new code to determine the bottom neighbours of the */
                        /* newly partitioned trapezoid */

                        tr[tn].d0 = tr[t].d1;
                        tr[tn].d1 = -1;

                        tnext = tr[t].d1;
                  }

                  t = tnext;
            }

            tr[t_sav].rseg = tr[tn_sav].lseg  = segnum;
      } /* end-while */

      /* Now combine those trapezoids which share common segments. We can */
      /* use the pointers to the parent to connect these together. This */
      /* works only because all these new trapezoids have been formed */
      /* due to splitting by the segment, and hence have only one parent */

      tfirstl = tfirst;
      tlastl = tlast;
      merge_trapezoids(segnum, tfirstl, tlastl, S_LEFT);
      merge_trapezoids(segnum, tfirstr, tlastr, S_RIGHT);

      seg[segnum].is_inserted = TRUE;
      return 0;
}


/* Update the roots stored for each of the endpoints of the segment.
 * This is done to speed up the location-query for the endpoint when
 * the segment is inserted into the trapezoidation subsequently
 */
static int find_new_roots(int segnum)
{
      segment_t *s = &seg[segnum];

      if (s->is_inserted)
            return 0;

      s->root0 = locate_endpoint(&s->v0, &s->v1, s->root0);
      s->root0 = tr[s->root0].sink;

      s->root1 = locate_endpoint(&s->v1, &s->v0, s->root1);
      s->root1 = tr[s->root1].sink;
      return 0;
}


/* Main routine to perform trapezoidation */
int construct_trapezoids(int nseg)
{
      register int i;
      int root, h;

      /* Add the first segment and get the query structure and trapezoid */
      /* list initialised */

      root = init_query_structure(choose_segment());

      for (i = 1; i <= nseg; i++)
            seg[i].root0 = seg[i].root1 = root;

      for (h = 1; h <= math_logstar_n(nseg); h++) {
            for (i = math_N(nseg, h - 1) + 1; i <= math_N(nseg, h); i++)
                  add_segment(choose_segment());

            /* Find a new root for each of the segment endpoints */
            for (i = 1; i <= nseg; i++)
                  find_new_roots(i);
      }

      for (i = math_N(nseg, math_logstar_n(nseg)) + 1; i <= nseg; i++)
            add_segment(choose_segment());

      return 0;
}



Generated by  Doxygen 1.6.0   Back to index