home *** CD-ROM | disk | FTP | other *** search
/ RISC DISC 3 / RISC_DISC_3.iso / resources / etexts / gems / gemsv / ch7_5 / construct.c < prev    next >
Encoding:
C/C++ Source or Header  |  1995-04-04  |  24.2 KB  |  1,060 lines
  1. #include "basic.h"
  2. #include <math.h>
  3. #include <string.h>    /* for memset() */
  4.  
  5. static int q_idx;
  6. static int tr_idx;
  7.  
  8. /* Return a new node to be added into the query tree */
  9. static int newnode()
  10. {
  11.   if (q_idx < QSIZE)
  12.     return q_idx++;
  13.   else
  14.     {
  15.       fprintf(stderr, "newnode: Query-table overflow\n");
  16.       return -1;
  17.     }
  18. }
  19.  
  20. /* Return a free trapezoid */
  21. static int newtrap()
  22. {
  23.   if (tr_idx < TRSIZE)
  24.     {
  25.       tr[tr_idx].lseg = -1;
  26.       tr[tr_idx].rseg = -1;
  27.       tr[tr_idx].state = ST_VALID;
  28.       return tr_idx++;
  29.     }
  30.   else
  31.     {
  32.       fprintf(stderr, "newtrap: Trapezoid-table overflow\n");
  33.       return -1;
  34.     }
  35. }
  36.  
  37.  
  38. /* Return the maximum of the two points into the yval structure */
  39. static int _max(yval, v0, v1)
  40.      point_t *yval;
  41.      point_t *v0;
  42.      point_t *v1;
  43. {
  44.   if (v0->y > v1->y + C_EPS)
  45.     *yval = *v0;
  46.   else if (FP_EQUAL(v0->y, v1->y))
  47.     {
  48.       if (v0->x > v1->x + C_EPS)
  49.     *yval = *v0;
  50.       else
  51.     *yval = *v1;
  52.     }
  53.   else
  54.     *yval = *v1;
  55.   
  56.   return 0;
  57. }
  58.  
  59.  
  60. /* Return the minimum of the two points into the yval structure */
  61. static int _min(yval, v0, v1)
  62.      point_t *yval;
  63.      point_t *v0;
  64.      point_t *v1;
  65. {
  66.   if (v0->y < v1->y - C_EPS)
  67.     *yval = *v0;
  68.   else if (FP_EQUAL(v0->y, v1->y))
  69.     {
  70.       if (v0->x < v1->x)
  71.     *yval = *v0;
  72.       else
  73.     *yval = *v1;
  74.     }
  75.   else
  76.     *yval = *v1;
  77.   
  78.   return 0;
  79. }
  80.  
  81.  
  82. int _greater_than(v0, v1)
  83.      point_t *v0;
  84.      point_t *v1;
  85. {
  86.   if (v0->y > v1->y + C_EPS)
  87.     return TRUE;
  88.   else if (v0->y < v1->y - C_EPS)
  89.     return FALSE;
  90.   else
  91.     return (v0->x > v1->x);
  92. }
  93.  
  94.  
  95. int _equal_to(v0, v1)
  96.      point_t *v0;
  97.      point_t *v1;
  98. {
  99.   return (FP_EQUAL(v0->y, v1->y) && FP_EQUAL(v0->x, v1->x));
  100. }
  101.  
  102. int _greater_than_equal_to(v0, v1)
  103.      point_t *v0;
  104.      point_t *v1;
  105. {
  106.   if (v0->y > v1->y + C_EPS)
  107.     return TRUE;
  108.   else if (v0->y < v1->y - C_EPS)
  109.     return FALSE;
  110.   else
  111.     return (v0->x >= v1->x);
  112. }
  113.  
  114. int _less_than(v0, v1)
  115.      point_t *v0;
  116.      point_t *v1;
  117. {
  118.   if (v0->y < v1->y - C_EPS)
  119.     return TRUE;
  120.   else if (v0->y > v1->y + C_EPS)
  121.     return FALSE;
  122.   else
  123.     return (v0->x < v1->x);
  124. }
  125.  
  126.  
  127. /* Initilialise the query structure (Q) and the trapezoid table (T) 
  128.  * when the first segment is added to start the trapezoidation
  129.  */
  130. int init_query_structure(segnum)
  131.      int segnum;
  132. {
  133.   int i1, i2, i3, i4, i5, i6, i7, root;
  134.   int t1, t2, t3, t4;
  135.   segment_t *s = &seg[segnum];
  136.  
  137.   memset((void *)tr, 0, sizeof(tr));
  138.   memset((void *)qs, 0, sizeof(qs));
  139.  
  140.   i1 = newnode();
  141.   qs[i1].nodetype = T_Y;
  142.   _max(&qs[i1].yval, &s->v0, &s->v1); /* root */
  143.   root = i1;
  144.  
  145.   qs[i1].right = i2 = newnode();
  146.   qs[i2].nodetype = T_SINK;
  147.   qs[i2].parent = i1;
  148.  
  149.   qs[i1].left = i3 = newnode();
  150.   qs[i3].nodetype = T_Y;
  151.   _min(&qs[i3].yval, &s->v0, &s->v1); /* root */
  152.   qs[i3].parent = i1;
  153.   
  154.   qs[i3].left = i4 = newnode();
  155.   qs[i4].nodetype = T_SINK;
  156.   qs[i4].parent = i3;
  157.   
  158.   qs[i3].right = i5 = newnode();
  159.   qs[i5].nodetype = T_X;
  160.   qs[i5].segnum = segnum;
  161.   qs[i5].parent = i3;
  162.   
  163.   qs[i5].left = i6 = newnode();
  164.   qs[i6].nodetype = T_SINK;
  165.   qs[i6].parent = i5;
  166.  
  167.   qs[i5].right = i7 = newnode();
  168.   qs[i7].nodetype = T_SINK;
  169.   qs[i7].parent = i5;
  170.  
  171.   t1 = newtrap();        /* middle left */
  172.   t2 = newtrap();        /* middle right */
  173.   t3 = newtrap();        /* bottom-most */
  174.   t4 = newtrap();        /* topmost */
  175.  
  176.   tr[t1].hi = tr[t2].hi = tr[t4].lo = qs[i1].yval;
  177.   tr[t1].lo = tr[t2].lo = tr[t3].hi = qs[i3].yval;
  178.   tr[t4].hi.y = (double) (INFINITY);
  179.   tr[t4].hi.x = (double) (INFINITY);
  180.   tr[t3].lo.y = (double) -1* (INFINITY);
  181.   tr[t3].lo.x = (double) -1* (INFINITY);
  182.   tr[t1].rseg = tr[t2].lseg = segnum;
  183.   tr[t1].u0 = tr[t2].u0 = t4;
  184.   tr[t1].d0 = tr[t2].d0 = t3;
  185.   tr[t4].d0 = tr[t3].u0 = t1;
  186.   tr[t4].d1 = tr[t3].u1 = t2;
  187.   
  188.   tr[t1].sink = i6;
  189.   tr[t2].sink = i7;
  190.   tr[t3].sink = i4;
  191.   tr[t4].sink = i2;
  192.  
  193.   tr[t1].state = tr[t2].state = ST_VALID;
  194.   tr[t3].state = tr[t4].state = ST_VALID;
  195.  
  196.   qs[i2].trnum = t4;
  197.   qs[i4].trnum = t3;
  198.   qs[i6].trnum = t1;
  199.   qs[i7].trnum = t2;
  200.  
  201.   s->is_inserted = TRUE;
  202.   return root;
  203. }
  204.  
  205.  
  206. /* Retun TRUE if the vertex v is to the left of line segment no.
  207.    * segnum
  208.    */
  209.  
  210. static int is_left_of(segnum, v)
  211.      int segnum;
  212.      point_t *v;
  213. {
  214.   segment_t *s = &seg[segnum];
  215.   double area;
  216.   
  217.   if (_greater_than(&s->v1, &s->v0)) /* seg. going upwards */
  218.     {
  219.       if (FP_EQUAL(s->v1.y, v->y))
  220.     {
  221.       if (v->x < s->v1.x)
  222.         area = 1.0;
  223.       else
  224.         area = -1.0;
  225.     }
  226.       else if (FP_EQUAL(s->v0.y, v->y))
  227.     {
  228.       if (v->x < s->v0.x)
  229.         area = 1.0;
  230.       else
  231.         area = -1.0;
  232.     }
  233.       else
  234.     area = CROSS(s->v0, s->v1, (*v));
  235.     }
  236.   else                /* v0 > v1 */
  237.     {
  238.       if (FP_EQUAL(s->v1.y, v->y))
  239.     {
  240.       if (v->x < s->v1.x)
  241.         area = 1.0;
  242.       else
  243.         area = -1.0;
  244.     }
  245.       else if (FP_EQUAL(s->v0.y, v->y))
  246.     {
  247.       if (v->x < s->v0.x)
  248.         area = 1.0;
  249.       else
  250.         area = -1.0;
  251.     }
  252.       else
  253.     area = CROSS(s->v1, s->v0, (*v));
  254.     }
  255.   
  256.   if (area > 0.0)
  257.     return TRUE;
  258.   else 
  259.     return FALSE;
  260. }
  261.  
  262.  
  263. int is_collinear(segnum, v, is_swapped)
  264.      int segnum;
  265.      point_t *v;
  266.      int is_swapped;
  267. {
  268.   int n;
  269.   
  270.   /* First check if the endpoint is already inserted */
  271.   if (!is_swapped)
  272.     n = MODULO_NEXT(segnum + 1, global.nseg);
  273.   else
  274.     if ((n = segnum - 1) == 0)
  275.       n = 1;
  276.   
  277.   return seg[n].is_inserted;
  278. }
  279.  
  280.  
  281. /* This is query routine which determines which trapezoid does the 
  282.  * point v lie in. The return value is the trapezoid number 
  283.  */
  284.  
  285. int locate_endpoint(v, vo, r)
  286.      point_t *v;
  287.      point_t *vo;
  288.      int r;
  289. {
  290.   node_t *rptr = &qs[r];
  291.   
  292.   switch (rptr->nodetype)
  293.     {
  294.     case T_SINK:
  295.       return rptr->trnum;
  296.       
  297.     case T_Y:
  298.       if (_greater_than(v, &rptr->yval)) /* above */
  299.     return locate_endpoint(v, vo, rptr->right);
  300.       else if (_equal_to(v, &rptr->yval)) /* the point is already */
  301.     {                      /* inserted. */
  302.       if (_greater_than(vo, &rptr->yval)) /* above */
  303.         return locate_endpoint(v, vo, rptr->right);
  304.       else 
  305.         return locate_endpoint(v, vo, rptr->left); /* below */        
  306.     }
  307.       else
  308.     return locate_endpoint(v, vo, rptr->left); /* below */
  309.  
  310.     case T_X:
  311.       if (_equal_to(v, &seg[rptr->segnum].v0) || 
  312.            _equal_to(v, &seg[rptr->segnum].v1))
  313.     {
  314.       if (FP_EQUAL(v->y, vo->y)) /* horizontal segment */
  315.         {
  316.           if (vo->x < v->x)
  317.         return locate_endpoint(v, vo, rptr->left); /* left */
  318.           else
  319.         return locate_endpoint(v, vo, rptr->right); /* right */
  320.         }
  321.  
  322.       else if (is_left_of(rptr->segnum, vo))
  323.         return locate_endpoint(v, vo, rptr->left); /* left */
  324.       else
  325.         return locate_endpoint(v, vo, rptr->right); /* right */
  326.     }
  327.       else if (is_left_of(rptr->segnum, v))
  328.     return locate_endpoint(v, vo, rptr->left); /* left */
  329.       else
  330.     return locate_endpoint(v, vo, rptr->right); /* right */    
  331.  
  332.     default:
  333.       fprintf(stderr, "Haggu !!!!!\n");
  334.       break;
  335.     }
  336. }
  337.  
  338.  
  339. /* Thread in the segment into the existing trapezoidation. The 
  340.  * limiting trapezoids are given by tfirst and tlast (which are the
  341.  * trapezoids containing the two endpoints of the segment 
  342.  */
  343.  
  344. int merge_trapezoids(segnum, tfirst, tlast, side)
  345.      int segnum;
  346.      int tfirst;
  347.      int tlast;
  348.      int side;
  349. {
  350.   int t, tnext, cond;
  351.   int ptnext;
  352.  
  353.   /* First merge polys on the LHS */
  354.   t = tfirst;
  355.   while ((t > 0) && _greater_than_equal_to(&tr[t].lo, &tr[tlast].lo))
  356.     {
  357.       if (side == S_LEFT)
  358.     cond = ((((tnext = tr[t].d0) > 0) && (tr[tnext].rseg == segnum)) ||
  359.         (((tnext = tr[t].d1) > 0) && (tr[tnext].rseg == segnum)));
  360.       else
  361.     cond = ((((tnext = tr[t].d0) > 0) && (tr[tnext].lseg == segnum)) ||
  362.         (((tnext = tr[t].d1) > 0) && (tr[tnext].lseg == segnum)));
  363.       
  364.       if (cond)
  365.     {
  366.       if ((tr[t].lseg == tr[tnext].lseg) &&
  367.           (tr[t].rseg == tr[tnext].rseg)) /* good neighbours */
  368.         {                          /* merge them */
  369.           /* Use the upper node as the new node i.e. t */
  370.           
  371.           ptnext = qs[tr[tnext].sink].parent;
  372.           
  373.           if (qs[ptnext].left == tr[tnext].sink)
  374.         qs[ptnext].left = tr[t].sink;
  375.           else
  376.         qs[ptnext].right = tr[t].sink;    /* redirect parent */
  377.           
  378.           
  379.           /* Change the upper neighbours of the lower trapezoids */
  380.           
  381.           if ((tr[t].d0 = tr[tnext].d0) > 0)
  382.         if (tr[tr[t].d0].u0 == tnext)
  383.           tr[tr[t].d0].u0 = t;
  384.         else if (tr[tr[t].d0].u1 == tnext)
  385.           tr[tr[t].d0].u1 = t;
  386.           
  387.           if ((tr[t].d1 = tr[tnext].d1) > 0)
  388.         if (tr[tr[t].d1].u0 == tnext)
  389.           tr[tr[t].d1].u0 = t;
  390.         else if (tr[tr[t].d1].u1 == tnext)
  391.           tr[tr[t].d1].u1 = t;
  392.           
  393.           tr[t].lo = tr[tnext].lo;
  394.           tr[tnext].state = ST_INVALID; /* invalidate the lower */
  395.                             /* trapezium */
  396.         }
  397.       else            /* not good neighbours */
  398.         t = tnext;
  399.     }
  400.       else            /* do not satisfy the outer if */
  401.     t = tnext;
  402.       
  403.     } /* end-while */
  404.        
  405.   return 0;
  406. }
  407.  
  408.  
  409. /* Add in the new segment into the trapezoidation and update Q and T
  410.  * structures 
  411.  */
  412. int add_segment(segnum)
  413.      int segnum;
  414. {
  415.   segment_t s;
  416.   int tu, tl, sk, tfirst, tlast, tnext;
  417.   int tfirstr, tlastr, tfirstl, tlastl;
  418.   int i1, i2, t, tn;
  419.   point_t vper, tpt;
  420.   int tritop = 0, tribot = 0, is_swapped = 0;
  421.   int tmptriseg;
  422.  
  423.   s = seg[segnum];
  424.   if (_greater_than(&s.v1, &s.v0)) /* Get higher vertex in v0 */
  425.     {
  426.       int tmp;
  427.       tpt = s.v0;
  428.       s.v0 = s.v1;
  429.       s.v1 = tpt;
  430.       tmp = s.root0;
  431.       s.root0 = s.root1;
  432.       s.root1 = tmp;
  433.       is_swapped = TRUE;
  434.     }
  435.  
  436.   if ((is_swapped) ? !inserted(segnum, LASTPT) :
  437.        !inserted(segnum, FIRSTPT))     /* insert v0 in the tree */
  438.     {
  439.       int tmp_d;
  440.  
  441.       tu = locate_endpoint(&s.v0, &s.v1, s.root0);
  442.       tl = newtrap();        /* tl is the new lower trapezoid */
  443.       tr[tl].state = ST_VALID;
  444.       tr[tl] = tr[tu];
  445.       tr[tu].lo.y = tr[tl].hi.y = s.v0.y;
  446.       tr[tu].lo.x = tr[tl].hi.x = s.v0.x;
  447.       tr[tu].d0 = tl;      
  448.       tr[tu].d1 = 0;
  449.       tr[tl].u0 = tu;
  450.       tr[tl].u1 = 0;
  451.  
  452.       if (((tmp_d = tr[tl].d0) > 0) && (tr[tmp_d].u0 == tu))
  453.     tr[tmp_d].u0 = tl;
  454.       if (((tmp_d = tr[tl].d0) > 0) && (tr[tmp_d].u1 == tu))
  455.     tr[tmp_d].u1 = tl;
  456.  
  457.       if (((tmp_d = tr[tl].d1) > 0) && (tr[tmp_d].u0 == tu))
  458.     tr[tmp_d].u0 = tl;
  459.       if (((tmp_d = tr[tl].d1) > 0) && (tr[tmp_d].u1 == tu))
  460.     tr[tmp_d].u1 = tl;
  461.  
  462.       /* Now update the query structure and obtain the sinks for the */
  463.       /* two trapezoids */ 
  464.       
  465.       i1 = newnode();        /* Upper trapezoid sink */
  466.       i2 = newnode();        /* Lower trapezoid sink */
  467.       sk = tr[tu].sink;
  468.       
  469.       qs[sk].nodetype = T_Y;
  470.       qs[sk].yval = s.v0;
  471.       qs[sk].segnum = segnum;    /* not really reqd ... maybe later */
  472.       qs[sk].left = i2;
  473.       qs[sk].right = i1;
  474.  
  475.       qs[i1].nodetype = T_SINK;
  476.       qs[i1].trnum = tu;
  477.       qs[i1].parent = sk;
  478.  
  479.       qs[i2].nodetype = T_SINK;
  480.       qs[i2].trnum = tl;
  481.       qs[i2].parent = sk;
  482.  
  483.       tr[tu].sink = i1;
  484.       tr[tl].sink = i2;
  485.       tfirst = tl;
  486.     }
  487.   else                /* v0 already present */
  488.     {       /* Get the topmost intersecting trapezoid */
  489.       vper.x = s.v0.x + EPS * (s.v1.x - s.v0.x);
  490.       vper.y = s.v0.y + EPS * (s.v1.y - s.v0.y);
  491.       tfirst = locate_endpoint(&s.v0, &s.v1, s.root0);
  492.       tritop = 1;
  493.     }
  494.  
  495.  
  496.   if ((is_swapped) ? !inserted(segnum, FIRSTPT) :
  497.        !inserted(segnum, LASTPT))     /* insert v1 in the tree */
  498.     {
  499.       int tmp_d;
  500.  
  501.       tu = locate_endpoint(&s.v1, &s.v0, s.root1);
  502.  
  503.       tl = newtrap();        /* tl is the new lower trapezoid */
  504.       tr[tl].state = ST_VALID;
  505.       tr[tl] = tr[tu];
  506.       tr[tu].lo.y = tr[tl].hi.y = s.v1.y;
  507.       tr[tu].lo.x = tr[tl].hi.x = s.v1.x;
  508.       tr[tu].d0 = tl;      
  509.       tr[tu].d1 = 0;
  510.       tr[tl].u0 = tu;
  511.       tr[tl].u1 = 0;
  512.  
  513.       if (((tmp_d = tr[tl].d0) > 0) && (tr[tmp_d].u0 == tu))
  514.     tr[tmp_d].u0 = tl;
  515.       if (((tmp_d = tr[tl].d0) > 0) && (tr[tmp_d].u1 == tu))
  516.     tr[tmp_d].u1 = tl;
  517.  
  518.       if (((tmp_d = tr[tl].d1) > 0) && (tr[tmp_d].u0 == tu))
  519.     tr[tmp_d].u0 = tl;
  520.       if (((tmp_d = tr[tl].d1) > 0) && (tr[tmp_d].u1 == tu))
  521.     tr[tmp_d].u1 = tl;
  522.       
  523.       /* Now update the query structure and obtain the sinks for the */
  524.       /* two trapezoids */ 
  525.       
  526.       i1 = newnode();        /* Upper trapezoid sink */
  527.       i2 = newnode();        /* Lower trapezoid sink */
  528.       sk = tr[tu].sink;
  529.       
  530.       qs[sk].nodetype = T_Y;
  531.       qs[sk].yval = s.v1;
  532.       qs[sk].segnum = segnum;    /* not really reqd ... maybe later */
  533.       qs[sk].left = i2;
  534.       qs[sk].right = i1;
  535.  
  536.       qs[i1].nodetype = T_SINK;
  537.       qs[i1].trnum = tu;
  538.       qs[i1].parent = sk;
  539.  
  540.       qs[i2].nodetype = T_SINK;
  541.       qs[i2].trnum = tl;
  542.       qs[i2].parent = sk;
  543.  
  544.       tr[tu].sink = i1;
  545.       tr[tl].sink = i2;
  546.       tlast = tu;
  547.     }
  548.   else                /* v1 already present */
  549.     {       /* Get the lowermost intersecting trapezoid */
  550.       vper.x = s.v1.x + EPS * (s.v0.x - s.v1.x);
  551.       vper.y = s.v1.y + EPS * (s.v0.y - s.v1.y);
  552.       tlast = locate_endpoint(&s.v1, &s.v0, s.root1);
  553.       tribot = 1;
  554.     }
  555.   
  556.   /* Thread the segment into the query tree creating a new X-node */
  557.   /* First, split all the trapezoids which are intersected by s into */
  558.   /* two */
  559.  
  560.   t = tfirst;            /* topmost trapezoid */
  561.   
  562.   while ((t > 0) && 
  563.      _greater_than_equal_to(&tr[t].lo, &tr[tlast].lo))
  564.                 /* traverse from top to bot */
  565.     {
  566.       int t_sav, tn_sav;
  567.       sk = tr[t].sink;
  568.       i1 = newnode();        /* left trapezoid sink */
  569.       i2 = newnode();        /* right trapezoid sink */
  570.       
  571.       qs[sk].nodetype = T_X;
  572.       qs[sk].segnum = segnum;
  573.       qs[sk].left = i1;
  574.       qs[sk].right = i2;
  575.  
  576.       qs[i1].nodetype = T_SINK;    /* left trapezoid (use existing one) */
  577.       qs[i1].trnum = t;
  578.       qs[i1].parent = sk;
  579.  
  580.       qs[i2].nodetype = T_SINK;    /* right trapezoid (allocate new) */
  581.       qs[i2].trnum = tn = newtrap();
  582.       tr[tn].state = ST_VALID;
  583.       qs[i2].parent = sk;
  584.  
  585.       if (t == tfirst)
  586.     tfirstr = tn;
  587.       if (_equal_to(&tr[t].lo, &tr[tlast].lo))
  588.     tlastr = tn;
  589.  
  590.       tr[tn] = tr[t];
  591.       tr[t].sink = i1;
  592.       tr[tn].sink = i2;
  593.       t_sav = t;
  594.       tn_sav = tn;
  595.  
  596.       /* error */
  597.  
  598.       if ((tr[t].d0 <= 0) && (tr[t].d1 <= 0)) /* case cannot arise */
  599.     {
  600.       fprintf(stderr, "add_segment: error\n");
  601.       break;
  602.     }
  603.       
  604.       /* only one trapezoid below. partition t into two and make the */
  605.       /* two resulting trapezoids t and tn as the upper neighbours of */
  606.       /* the sole lower trapezoid */
  607.       
  608.       else if ((tr[t].d0 > 0) && (tr[t].d1 <= 0))
  609.     {            /* Only one trapezoid below */
  610.       if ((tr[t].u0 > 0) && (tr[t].u1 > 0))
  611.         {            /* continuation of a chain from abv. */
  612.           if (tr[t].usave > 0) /* three upper neighbours */
  613.         {
  614.           if (tr[t].uside == S_LEFT)
  615.             {
  616.               tr[tn].u0 = tr[t].u1;
  617.               tr[t].u1 = -1;
  618.               tr[tn].u1 = tr[t].usave;
  619.               
  620.               tr[tr[t].u0].d0 = t;
  621.               tr[tr[tn].u0].d0 = tn;
  622.               tr[tr[tn].u1].d0 = tn;
  623.             }
  624.           else        /* intersects in the right */
  625.             {
  626.               tr[tn].u1 = -1;
  627.               tr[tn].u0 = tr[t].u1;
  628.               tr[t].u1 = tr[t].u0;
  629.               tr[t].u0 = tr[t].usave;
  630.  
  631.               tr[tr[t].u0].d0 = t;
  632.               tr[tr[t].u1].d0 = t;
  633.               tr[tr[tn].u0].d0 = tn;              
  634.             }
  635.           
  636.           tr[t].usave = tr[tn].usave = 0;
  637.         }
  638.           else        /* No usave.... simple case */
  639.         {
  640.           tr[tn].u0 = tr[t].u1;
  641.           tr[t].u1 = tr[tn].u1 = -1;
  642.           tr[tr[tn].u0].d0 = tn;
  643.         }
  644.         }
  645.       else 
  646.         {            /* fresh seg. or upward cusp */
  647.           int tmp_u = tr[t].u0;
  648.           int td0, td1;
  649.           if (((td0 = tr[tmp_u].d0) > 0) && 
  650.           ((td1 = tr[tmp_u].d1) > 0))
  651.         {        /* upward cusp */
  652.           if ((tr[td0].rseg > 0) &&
  653.               !is_left_of(tr[td0].rseg, &s.v1))
  654.             {
  655.               tr[t].u0 = tr[t].u1 = tr[tn].u1 = -1;
  656.               tr[tr[tn].u0].d1 = tn;
  657.             }
  658.           else        /* cusp going leftwards */
  659.             { 
  660.               tr[tn].u0 = tr[tn].u1 = tr[t].u1 = -1;
  661.               tr[tr[t].u0].d0 = t;
  662.             }
  663.         }
  664.           else        /* fresh segment */
  665.         {
  666.           tr[tr[t].u0].d0 = t;
  667.           tr[tr[t].u0].d1 = tn;
  668.         }          
  669.         }
  670.       
  671.       if (FP_EQUAL(tr[t].lo.y, tr[tlast].lo.y) && 
  672.           FP_EQUAL(tr[t].lo.x, tr[tlast].lo.x) && tribot)
  673.         {        /* bottom forms a triangle */
  674.  
  675.           if (is_swapped)
  676.         {
  677.           tmptriseg = segnum - 1;
  678.           if (tmptriseg == 0)
  679.             tmptriseg = global.nseg;
  680.         }
  681.           else
  682.         tmptriseg = MODULO_NEXT(segnum + 1, global.nseg);
  683.           
  684.           if ((tmptriseg > 0) && is_left_of(tmptriseg, &s.v0))
  685.         {
  686.                 /* L-R downward cusp */
  687.           tr[tr[t].d0].u0 = t;
  688.           tr[tn].d0 = tr[tn].d1 = -1;
  689.         }
  690.           else
  691.         {
  692.                 /* R-L downward cusp */
  693.           tr[tr[tn].d0].u1 = tn;
  694.           tr[t].d0 = tr[t].d1 = -1;
  695.         }
  696.         }
  697.       else
  698.         {
  699.           if ((tr[tr[t].d0].u0 > 0) && (tr[tr[t].d0].u1 > 0))
  700.         {
  701.           if (tr[tr[t].d0].u0 == t) /* passes thru LHS */
  702.             {
  703.               tr[tr[t].d0].usave = tr[tr[t].d0].u1;
  704.               tr[tr[t].d0].uside = S_LEFT;
  705.             }
  706.           else
  707.             {
  708.               tr[tr[t].d0].usave = tr[tr[t].d0].u0;
  709.               tr[tr[t].d0].uside = S_RIGHT;
  710.             }            
  711.         }
  712.           tr[tr[t].d0].u0 = t;
  713.           tr[tr[t].d0].u1 = tn;
  714.         }
  715.       
  716.       t = tr[t].d0;
  717.     }
  718.  
  719.  
  720.       else if ((tr[t].d0 <= 0) && (tr[t].d1 > 0))
  721.     {            /* Only one trapezoid below */
  722.       if ((tr[t].u0 > 0) && (tr[t].u1 > 0))
  723.         {            /* continuation of a chain from abv. */
  724.           if (tr[t].usave > 0) /* three upper neighbours */
  725.         {
  726.           if (tr[t].uside == S_LEFT)
  727.             {
  728.               tr[tn].u0 = tr[t].u1;
  729.               tr[t].u1 = -1;
  730.               tr[tn].u1 = tr[t].usave;
  731.               
  732.               tr[tr[t].u0].d0 = t;
  733.               tr[tr[tn].u0].d0 = tn;
  734.               tr[tr[tn].u1].d0 = tn;
  735.             }
  736.           else        /* intersects in the right */
  737.             {
  738.               tr[tn].u1 = -1;
  739.               tr[tn].u0 = tr[t].u1;
  740.               tr[t].u1 = tr[t].u0;
  741.               tr[t].u0 = tr[t].usave;
  742.  
  743.               tr[tr[t].u0].d0 = t;
  744.               tr[tr[t].u1].d0 = t;
  745.               tr[tr[tn].u0].d0 = tn;              
  746.             }
  747.           
  748.           tr[t].usave = tr[tn].usave = 0;
  749.         }
  750.           else        /* No usave.... simple case */
  751.         {
  752.           tr[tn].u0 = tr[t].u1;
  753.           tr[t].u1 = tr[tn].u1 = -1;
  754.           tr[tr[tn].u0].d0 = tn;
  755.         }
  756.         }
  757.       else 
  758.         {            /* fresh seg. or upward cusp */
  759.           int tmp_u = tr[t].u0;
  760.           int td0, td1;
  761.           if (((td0 = tr[tmp_u].d0) > 0) && 
  762.           ((td1 = tr[tmp_u].d1) > 0))
  763.         {        /* upward cusp */
  764.           if ((tr[td0].rseg > 0) &&
  765.               !is_left_of(tr[td0].rseg, &s.v1))
  766.             {
  767.               tr[t].u0 = tr[t].u1 = tr[tn].u1 = -1;
  768.               tr[tr[tn].u0].d1 = tn;
  769.             }
  770.           else 
  771.             {
  772.               tr[tn].u0 = tr[tn].u1 = tr[t].u1 = -1;
  773.               tr[tr[t].u0].d0 = t;
  774.             }
  775.         }
  776.           else        /* fresh segment */
  777.         {
  778.           tr[tr[t].u0].d0 = t;
  779.           tr[tr[t].u0].d1 = tn;
  780.         }
  781.         }
  782.       
  783.       if (FP_EQUAL(tr[t].lo.y, tr[tlast].lo.y) && 
  784.           FP_EQUAL(tr[t].lo.x, tr[tlast].lo.x) && tribot)
  785.         {        /* bottom forms a triangle */
  786.           int tmpseg;
  787.           if (is_swapped)
  788.         {
  789.           tmpseg = segnum - 1;
  790.           if (tmpseg == 0)
  791.             tmpseg = global.nseg;
  792.         }
  793.           else
  794.         tmpseg = MODULO_NEXT(segnum + 1, global.nseg);
  795.           
  796.           if ((tmpseg > 0) && is_left_of(tmpseg, &s.v0))
  797.         {
  798.           /* L-R downward cusp */
  799.           tr[tr[t].d1].u0 = t;
  800.           tr[tn].d0 = tr[tn].d1 = -1;
  801.         }
  802.           else
  803.         {
  804.           /* R-L downward cusp */
  805.           tr[tr[tn].d1].u1 = tn;
  806.           tr[t].d0 = tr[t].d1 = -1;
  807.         }
  808.         }        
  809.       else
  810.         {
  811.           if ((tr[tr[t].d1].u0 > 0) && (tr[tr[t].d1].u1 > 0))
  812.         {
  813.           if (tr[tr[t].d1].u0 == t) /* passes thru LHS */
  814.             {
  815.               tr[tr[t].d1].usave = tr[tr[t].d1].u1;
  816.               tr[tr[t].d1].uside = S_LEFT;
  817.             }
  818.           else
  819.             {
  820.               tr[tr[t].d1].usave = tr[tr[t].d1].u0;
  821.               tr[tr[t].d1].uside = S_RIGHT;
  822.             }            
  823.         }
  824.           tr[tr[t].d1].u0 = t;
  825.           tr[tr[t].d1].u1 = tn;
  826.         }
  827.       
  828.       t = tr[t].d1;
  829.     }
  830.  
  831.       /* two trapezoids below. Find out which one is intersected by */
  832.       /* this segment and proceed down that one */
  833.       
  834.       else
  835.     {
  836.       double y0, yt;
  837.       point_t tmppt;
  838.       int i_d0;
  839.  
  840.       i_d0 = FALSE;
  841.       if (FP_EQUAL(tr[t].lo.y, s.v0.y))
  842.         {
  843.           if (tr[t].lo.x > s.v0.x)
  844.         i_d0 = TRUE;
  845.         }
  846.       else
  847.         {
  848.           tmppt.y = y0 = tr[t].lo.y;
  849.           yt = (y0 - s.v0.y)/(s.v1.y - s.v0.y);
  850.           tmppt.x = s.v0.x + yt * (s.v1.x - s.v0.x);
  851.           
  852.           if (_less_than(&tmppt, &tr[t].lo))
  853.         i_d0 = TRUE;
  854.         }
  855.       
  856.       /* check continuity from the top so that the lower-neighbour */
  857.       /* values are properly filled for the upper trapezoid */
  858.  
  859.       if ((tr[t].u0 > 0) && (tr[t].u1 > 0))
  860.         {            /* continuation of a chain from abv. */
  861.           if (tr[t].usave > 0) /* three upper neighbours */
  862.         {
  863.           if (tr[t].uside == S_LEFT)
  864.             {
  865.               tr[tn].u0 = tr[t].u1;
  866.               tr[t].u1 = -1;
  867.               tr[tn].u1 = tr[t].usave;
  868.               
  869.               tr[tr[t].u0].d0 = t;
  870.               tr[tr[tn].u0].d0 = tn;
  871.               tr[tr[tn].u1].d0 = tn;
  872.             }
  873.           else        /* intersects in the right */
  874.             {
  875.               tr[tn].u1 = -1;
  876.               tr[tn].u0 = tr[t].u1;
  877.               tr[t].u1 = tr[t].u0;
  878.               tr[t].u0 = tr[t].usave;
  879.  
  880.               tr[tr[t].u0].d0 = t;
  881.               tr[tr[t].u1].d0 = t;
  882.               tr[tr[tn].u0].d0 = tn;              
  883.             }
  884.           
  885.           tr[t].usave = tr[tn].usave = 0;
  886.         }
  887.           else        /* No usave.... simple case */
  888.         {
  889.           tr[tn].u0 = tr[t].u1;
  890.           tr[tn].u1 = -1;
  891.           tr[t].u1 = -1;
  892.           tr[tr[tn].u0].d0 = tn;
  893.         }
  894.         }
  895.       else 
  896.         {            /* fresh seg. or upward cusp */
  897.           int tmp_u = tr[t].u0;
  898.           int td0, td1;
  899.           if (((td0 = tr[tmp_u].d0) > 0) && 
  900.           ((td1 = tr[tmp_u].d1) > 0))
  901.         {        /* upward cusp */
  902.           if ((tr[td0].rseg > 0) &&
  903.               !is_left_of(tr[td0].rseg, &s.v1))
  904.             {
  905.               tr[t].u0 = tr[t].u1 = tr[tn].u1 = -1;
  906.               tr[tr[tn].u0].d1 = tn;
  907.             }
  908.           else 
  909.             {
  910.               tr[tn].u0 = tr[tn].u1 = tr[t].u1 = -1;
  911.               tr[tr[t].u0].d0 = t;
  912.             }
  913.         }
  914.           else        /* fresh segment */
  915.         {
  916.           tr[tr[t].u0].d0 = t;
  917.           tr[tr[t].u0].d1 = tn;
  918.         }
  919.         }
  920.       
  921.       if (FP_EQUAL(tr[t].lo.y, tr[tlast].lo.y) && 
  922.           FP_EQUAL(tr[t].lo.x, tr[tlast].lo.x) && tribot)
  923.         {
  924.           /* this case arises only at the lowest trapezoid.. i.e.
  925.          tlast, if the lower endpoint of the segment is
  926.          already inserted in the structure */
  927.           
  928.           tr[tr[t].d0].u0 = t;
  929.           tr[tr[t].d0].u1 = -1;
  930.           tr[tr[t].d1].u0 = tn;
  931.           tr[tr[t].d1].u1 = -1;
  932.  
  933.           tr[tn].d0 = tr[t].d1;
  934.           tr[t].d1 = tr[tn].d1 = -1;
  935.           
  936.           tnext = tr[t].d1;          
  937.         }
  938.       else if (i_d0)
  939.                 /* intersecting d0 */
  940.         {
  941.           tr[tr[t].d0].u0 = t;
  942.           tr[tr[t].d0].u1 = tn;
  943.           tr[tr[t].d1].u0 = tn;
  944.           tr[tr[t].d1].u1 = -1;
  945.           
  946.           /* new code to determine the bottom neighbours of the */
  947.           /* newly partitioned trapezoid */
  948.           
  949.           tr[t].d1 = -1;
  950.  
  951.           tnext = tr[t].d0;
  952.         }
  953.       else            /* intersecting d1 */
  954.         {
  955.           tr[tr[t].d0].u0 = t;
  956.           tr[tr[t].d0].u1 = -1;
  957.           tr[tr[t].d1].u0 = t;
  958.           tr[tr[t].d1].u1 = tn;
  959.  
  960.           /* new code to determine the bottom neighbours of the */
  961.           /* newly partitioned trapezoid */
  962.           
  963.           tr[tn].d0 = tr[t].d1;
  964.           tr[tn].d1 = -1;
  965.           
  966.           tnext = tr[t].d1;
  967.         }        
  968.       
  969.       t = tnext;
  970.     }
  971.       
  972.       tr[t_sav].rseg = tr[tn_sav].lseg  = segnum;
  973.     } /* end-while */
  974.   
  975.   /* Now combine those trapezoids which share common segments. We can */
  976.   /* use the pointers to the parent to connect these together. This */
  977.   /* works only because all these new trapezoids have been formed */
  978.   /* due to splitting by the segment, and hence have only one parent */
  979.  
  980.   tfirstl = tfirst; 
  981.   tlastl = tlast;
  982.   merge_trapezoids(segnum, tfirstl, tlastl, S_LEFT);
  983.   merge_trapezoids(segnum, tfirstr, tlastr, S_RIGHT);
  984.  
  985.   seg[segnum].is_inserted = TRUE;
  986.   return 0;
  987. }
  988.  
  989.  
  990. /* Update the roots stored for each of the endpoints of the segment.
  991.  * This is done to speed up the location-query for the endpoint when
  992.  * the segment is inserted into the trapezoidation subsequently
  993.  */
  994. static int find_new_roots(segnum)
  995.      int segnum;
  996. {
  997.   segment_t *s = &seg[segnum];
  998.   point_t vper;
  999.   
  1000.   if (s->is_inserted)
  1001.     return 0;
  1002.  
  1003.   vper.x = s->v0.x + EPS * (s->v1.x - s->v0.x);
  1004.   vper.y = s->v0.y + EPS * (s->v1.y - s->v0.y);
  1005.   s->root0 = locate_endpoint(&s->v0, &s->v1, s->root0);
  1006.   s->root0 = tr[s->root0].sink;
  1007.  
  1008.   vper.x = s->v1.x + EPS * (s->v0.x - s->v1.x);
  1009.   vper.y = s->v1.y + EPS * (s->v0.y - s->v1.y);
  1010.   s->root1 = locate_endpoint(&s->v1, &s->v0, s->root1);
  1011.   s->root1 = tr[s->root1].sink;  
  1012.   return 0;
  1013. }
  1014.  
  1015.  
  1016. /* Main routine to perform trapezoidation */
  1017. int construct_trapezoids(nseg, seg)
  1018.      int nseg;
  1019.      segment_t *seg;
  1020. {
  1021.   register int i;
  1022.   int root, h;
  1023.   
  1024.   q_idx = tr_idx = 1;
  1025.   /* Add the first segment and get the query structure and trapezoid */
  1026.   /* list initialised */
  1027.   root = init_query_structure(choose_segment());
  1028.  
  1029. #ifdef SIMPLE            /* no randomization */
  1030.  
  1031.   for (i = 1; i <= nseg; i++)
  1032.     seg[i].root0 = seg[i].root1 = root;
  1033.   
  1034.   for (i = 2; i <= nseg; i++)
  1035.     add_segment(choose_segment());
  1036.  
  1037. #else
  1038.   
  1039.   for (i = 1; i <= nseg; i++)
  1040.     seg[i].root0 = seg[i].root1 = root;
  1041.   
  1042.   for (h = 1; h <= math_logstar_n(nseg); h++)
  1043.     {
  1044.       for (i = math_N(nseg, h -1) + 1; i <= math_N(nseg, h); i++)
  1045.     add_segment(choose_segment());
  1046.       
  1047.       /* Find a new root for each of the segment endpoints */
  1048.       for (i = 1; i <= nseg; i++)
  1049.     find_new_roots(i);
  1050.     }
  1051.   
  1052.   for (i = math_N(nseg, math_logstar_n(nseg)) + 1; i <= nseg; i++)
  1053.     add_segment(choose_segment());
  1054.   
  1055. #endif
  1056.  
  1057. }
  1058.  
  1059.  
  1060.