Panda3D

triangulator.cxx

00001 // Filename: triangulator.cxx
00002 // Created by:  drose (18Jan07)
00003 //
00004 ////////////////////////////////////////////////////////////////////
00005 //
00006 // PANDA 3D SOFTWARE
00007 // Copyright (c) Carnegie Mellon University.  All rights reserved.
00008 //
00009 // All use of this software is subject to the terms of the revised BSD
00010 // license.  You should have received a copy of this license along
00011 // with this source code in a file named "LICENSE."
00012 //
00013 ////////////////////////////////////////////////////////////////////
00014 
00015 #include "triangulator.h"
00016 #include "randomizer.h"
00017 
00018 ////////////////////////////////////////////////////////////////////
00019 //     Function: Triangulator::Constructor
00020 //       Access: Published
00021 //  Description: 
00022 ////////////////////////////////////////////////////////////////////
00023 Triangulator::
00024 Triangulator() {
00025 }
00026 
00027 ////////////////////////////////////////////////////////////////////
00028 //     Function: Triangulator::clear
00029 //       Access: Published
00030 //  Description: Removes all vertices and polygon specifications from
00031 //               the Triangulator, and prepares it to start over.
00032 ////////////////////////////////////////////////////////////////////
00033 void Triangulator::
00034 clear() {
00035   _vertices.clear();
00036   clear_polygon();
00037 }
00038 
00039 ////////////////////////////////////////////////////////////////////
00040 //     Function: Triangulator::add_vertex
00041 //       Access: Published
00042 //  Description: Adds a new vertex to the vertex pool.  Returns the
00043 //               vertex index number.
00044 ////////////////////////////////////////////////////////////////////
00045 int Triangulator::
00046 add_vertex(const LPoint2d &point) {
00047   int index = (int)_vertices.size();
00048   _vertices.push_back(point);
00049   return index;
00050 }
00051 
00052 ////////////////////////////////////////////////////////////////////
00053 //     Function: Triangulator::clear_polygon
00054 //       Access: Published
00055 //  Description: Removes the current polygon definition (and its set
00056 //               of holes), but does not clear the vertex pool.
00057 ////////////////////////////////////////////////////////////////////
00058 void Triangulator::
00059 clear_polygon() {
00060   _polygon.clear();
00061   _holes.clear();
00062 }
00063 
00064 ////////////////////////////////////////////////////////////////////
00065 //     Function: Triangulator::add_polygon_vertex
00066 //       Access: Published
00067 //  Description: Adds the next consecutive vertex of the polygon.
00068 //               This vertex should index into the vertex pool
00069 //               established by repeated calls to add_vertex().
00070 //
00071 //               The vertices may be listed in either clockwise or
00072 //               counterclockwise order.  Vertices should not be
00073 //               repeated.  In particular, do not repeat the first
00074 //               vertex at the end.
00075 ////////////////////////////////////////////////////////////////////
00076 void Triangulator::
00077 add_polygon_vertex(int index) {
00078   _polygon.push_back(index);
00079 }
00080 
00081 ////////////////////////////////////////////////////////////////////
00082 //     Function: Triangulator::begin_hole
00083 //       Access: Published
00084 //  Description: Finishes the previous hole, if any, and prepares to
00085 //               add a new hole.
00086 ////////////////////////////////////////////////////////////////////
00087 void Triangulator::
00088 begin_hole() {
00089   _holes.push_back(vector_int());
00090 }
00091 
00092 ////////////////////////////////////////////////////////////////////
00093 //     Function: Triangulator::add_hole_vertex
00094 //       Access: Published
00095 //  Description: Adds the next consecutive vertex of the current hole.
00096 //               This vertex should index into the vertex pool
00097 //               established by repeated calls to add_vertex().
00098 //
00099 //               The vertices may be listed in either clockwise or
00100 //               counterclockwise order.  Vertices should not be
00101 //               repeated.
00102 ////////////////////////////////////////////////////////////////////
00103 void Triangulator::
00104 add_hole_vertex(int index) {
00105   nassertv(!_holes.empty());
00106   _holes.back().push_back(index);
00107 }
00108 
00109 ////////////////////////////////////////////////////////////////////
00110 //     Function: Triangulator::triangulate
00111 //       Access: Published
00112 //  Description: Does the work of triangulating the specified polygon.
00113 //               After this call, you may retrieve the new triangles
00114 //               one at a time by iterating through
00115 //               get_triangle_v0/1/2().
00116 ////////////////////////////////////////////////////////////////////
00117 void Triangulator::
00118 triangulate() {
00119   _result.clear();
00120 
00121   // Set up the list of segments.
00122   seg.clear();
00123   seg.push_back(segment_t());  // we don't use the first entry.
00124   make_segment(_polygon, true);
00125 
00126   Holes::const_iterator hi;
00127   for (hi = _holes.begin(); hi != _holes.end(); ++hi) {
00128     make_segment(*hi, false);
00129   }
00130 
00131   // Shuffle the segment index.
00132   int num_segments = (int)seg.size() - 1;
00133   permute.reserve(num_segments);
00134   int i;
00135   for (i = 0; i < num_segments; ++i) {
00136     permute.push_back(i + 1);
00137   }
00138 
00139   // Actually, I'm not sure why we should shuffle the index.  That
00140   // makes the result non-deterministic, and isn't one order--for
00141   // instance, the initial order--as good as any other?
00142   /*
00143   Randomizer randomizer;
00144   for (i = 0; i < num_segments; ++i) {
00145     int j = randomizer.random_int(num_segments);
00146     nassertv(j >= 0 && j < num_segments);
00147     int t = permute[i];
00148     permute[i] = permute[j];
00149     permute[j] = t;
00150   }
00151   */
00152   choose_idx = 0;
00153 
00154   //  cerr << "got " << num_segments << " segments\n";
00155   /*
00156   for (i = 1; i < (int)seg.size(); ++i) {
00157     segment_t &s = seg[i];
00158     printf("  %d. (%g %g), (%g %g)\n", i, s.v0.x, s.v0.y, s.v1.x, s.v1.y);
00159     printf("    root0 = %d, root1 = %d\n", s.root0, s.root1);
00160     printf("    next = %d, prev = %d\n", s.next, s.prev);
00161   }
00162   */
00163 
00164   construct_trapezoids(num_segments);
00165 
00166   //  cerr << "got " << tr.size() - 1 << " trapezoids\n";
00167   /*
00168   for (i = 1; i < (int)tr.size(); ++i) {
00169     trap_t &t = tr[i];
00170     cerr << "  " << i << ". state = " << t.state << "\n";
00171     cerr << "    lseg = " << t.lseg << " rseg = " << t.rseg
00172          << "\n";
00173     cerr << "    hi = " << t.hi.x << " " << t.hi.y << " lo = "
00174          << t.lo.x << " " << t.lo.y << "\n";
00175   }
00176   */
00177 
00178   int nmonpoly = monotonate_trapezoids(num_segments);
00179 
00180   //  cerr << "got " << nmonpoly << " monotone polygons\n";
00181 
00182   triangulate_monotone_polygons(num_segments, nmonpoly);
00183 
00184   /*
00185   Result::iterator ri;
00186   for (ri = _result.begin(); ri != _result.end(); ++ri) {
00187     cerr << "tri: " << (*ri)._v0 << " " << (*ri)._v1 << " "
00188          << (*ri)._v2 << "\n";
00189   }
00190   */
00191 }
00192 
00193 ////////////////////////////////////////////////////////////////////
00194 //     Function: Triangulator::get_num_triangles
00195 //       Access: Published
00196 //  Description: Returns the number of triangles generated by the
00197 //               previous call to triangulate().
00198 ////////////////////////////////////////////////////////////////////
00199 int Triangulator::
00200 get_num_triangles() const {
00201   return _result.size();
00202 }
00203 
00204 ////////////////////////////////////////////////////////////////////
00205 //     Function: Triangulator::get_triangle_v0
00206 //       Access: Published
00207 //  Description: Returns vertex 0 of the nth triangle generated by the
00208 //               previous call to triangulate().
00209 //
00210 //               This is a zero-based index into the vertices added by
00211 //               repeated calls to add_vertex().
00212 ////////////////////////////////////////////////////////////////////
00213 int Triangulator::
00214 get_triangle_v0(int n) const {
00215   nassertr(n >= 0 && n < (int)_result.size(), -1);
00216   return _result[n]._v0;
00217 }
00218 
00219 ////////////////////////////////////////////////////////////////////
00220 //     Function: Triangulator::get_triangle_v1
00221 //       Access: Published
00222 //  Description: Returns vertex 1 of the nth triangle generated by the
00223 //               previous call to triangulate().
00224 //
00225 //               This is a zero-based index into the vertices added by
00226 //               repeated calls to add_vertex().
00227 ////////////////////////////////////////////////////////////////////
00228 int Triangulator::
00229 get_triangle_v1(int n) const {
00230   nassertr(n >= 0 && n < (int)_result.size(), -1);
00231   return _result[n]._v1;
00232 }
00233 
00234 ////////////////////////////////////////////////////////////////////
00235 //     Function: Triangulator::get_triangle_v2
00236 //       Access: Published
00237 //  Description: Returns vertex 2 of the nth triangle generated by the
00238 //               previous call to triangulate().
00239 //
00240 //               This is a zero-based index into the vertices added by
00241 //               repeated calls to add_vertex().
00242 ////////////////////////////////////////////////////////////////////
00243 int Triangulator::
00244 get_triangle_v2(int n) const {
00245   nassertr(n >= 0 && n < (int)_result.size(), -1);
00246   return _result[n]._v2;
00247 }
00248 
00249 
00250 // The remainder of the code in this file is adapted more or less from
00251 // the C code published with the referenced paper.
00252 
00253 #define T_X     1
00254 #define T_Y     2
00255 #define T_SINK  3
00256 
00257 
00258 #define FIRSTPT 1               /* checking whether pt. is inserted */ 
00259 #define LASTPT  2
00260 
00261 
00262 #define REALLY_BIG 1<<30
00263 #define C_EPS 1.0e-7            /* tolerance value: Used for making */
00264                                 /* all decisions about collinearity or */
00265                                 /* left/right of segment. Decrease */
00266                                 /* this value if the input points are */
00267                                 /* spaced very close together */
00268 
00269 
00270 #define S_LEFT 1                /* for merge-direction */
00271 #define S_RIGHT 2
00272 
00273 
00274 #define ST_VALID 1              /* for trapezium state */
00275 #define ST_INVALID 2
00276 
00277 
00278 #define SP_SIMPLE_LRUP 1        /* for splitting trapezoids */
00279 #define SP_SIMPLE_LRDN 2
00280 #define SP_2UP_2DN     3
00281 #define SP_2UP_LEFT    4
00282 #define SP_2UP_RIGHT   5
00283 #define SP_2DN_LEFT    6
00284 #define SP_2DN_RIGHT   7
00285 #define SP_NOSPLIT    -1
00286 
00287 #define TR_FROM_UP 1            /* for traverse-direction */
00288 #define TR_FROM_DN 2
00289 
00290 #define TRI_LHS 1
00291 #define TRI_RHS 2
00292 
00293 
00294 #define CROSS(v0, v1, v2) (((v1).x - (v0).x)*((v2).y - (v0).y) - \
00295                            ((v1).y - (v0).y)*((v2).x - (v0).x))
00296 
00297 #define DOT(v0, v1) ((v0).x * (v1).x + (v0).y * (v1).y)
00298 
00299 #define FP_EQUAL(s, t) (fabs(s - t) <= C_EPS)
00300 
00301 
00302 #define CROSS_SINE(v0, v1) ((v0).x * (v1).y - (v1).x * (v0).y)
00303 #define LENGTH(v0) (sqrt((v0).x * (v0).x + (v0).y * (v0).y))
00304 
00305 
00306 ////////////////////////////////////////////////////////////////////
00307 //     Function: Triangulator::check_left_winding
00308 //       Access: Private
00309 //  Description: Returns true if the list of vertices is
00310 //               counter-clockwise, false if it is clockwise.
00311 ////////////////////////////////////////////////////////////////////
00312 bool Triangulator::
00313 check_left_winding(const vector_int &range) const {
00314   // We do this by computing the polygon's signed area.  If it comes
00315   // out negative, the polygon is right-winding.
00316 
00317   double area = 0.0;
00318   size_t j = range.size() - 1;
00319   for (size_t i = 0; i < range.size(); ++i) {
00320     const LPoint2d &p0 = _vertices[range[j]];
00321     const LPoint2d &p1 = _vertices[range[i]];
00322     area += p0[0] * p1[1] - p0[1] * p1[0];
00323     j = i;
00324   }
00325 
00326   return area >= 0.0;
00327 }
00328 
00329 ////////////////////////////////////////////////////////////////////
00330 //     Function: Triangulator::make_segment
00331 //       Access: Private
00332 //  Description: Converts a linear list of integer vertices to a list
00333 //               of segment_t.  If want_left_winding is true, the list
00334 //               is reversed if necessary to make it left-winding;
00335 //               otherwise, it is reversed to make it right-winding.
00336 ////////////////////////////////////////////////////////////////////
00337 void Triangulator::
00338 make_segment(const vector_int &range, bool want_left_winding) {
00339   int num_points = (int)range.size();
00340   nassertv(num_points >= 2);
00341 
00342   int first = (int)seg.size();
00343   int last = first + num_points - 1;
00344 
00345   if (want_left_winding == check_left_winding(range)) {
00346     // Keep it in its natural order.
00347     int first = (int)seg.size();
00348     int last = first + num_points - 1;
00349     
00350     seg.push_back(segment_t(this, range[0], range[1],
00351                             last, first + 1));
00352     
00353     for (int i = 1; i < num_points - 1; ++i) {
00354       seg.push_back(segment_t(this, range[i], range[i + 1],
00355                               first + i - 1, first + i + 1));
00356     }
00357     
00358     seg.push_back(segment_t(this, range[num_points - 1], range[0],
00359                             last - 1, first));
00360 
00361   } else {
00362     // Reverse it.
00363     seg.push_back(segment_t(this, range[0], range[num_points - 1],
00364                             last, first + 1));
00365     
00366     for (int i = 1; i < num_points - 1; ++i) {
00367       seg.push_back(segment_t(this, range[num_points - i], range[num_points - i - 1],
00368                               first + i - 1, first + i + 1));
00369     }
00370     
00371     seg.push_back(segment_t(this, range[1], range[0],
00372                             last - 1, first));
00373 
00374   }
00375 }
00376   
00377 /* Return the next segment in the generated random ordering of all the */
00378 /* segments in S */
00379 int Triangulator::
00380 choose_segment() {
00381   nassertr(choose_idx < (int)permute.size(), 0);
00382   //  segment_t &s = seg[permute[choose_idx]];
00383   //  cerr << "choose_segment " << permute[choose_idx] << ": " << s.v0.x << ", " << s.v0.y << " to " << s.v1.x << ", " << s.v1.y << "\n";
00384   return permute[choose_idx++];
00385 }
00386 
00387 double Triangulator::
00388 math_log2(double v) {
00389   static const double log2 = log(2.0);
00390   return log(v) / log2;
00391 }
00392 
00393 /* Get log*n for given n */
00394 int Triangulator::
00395 math_logstar_n(int n) {
00396   int i;
00397   double v;
00398   
00399   for (i = 0, v = (double) n; v >= 1; i++)
00400     v = math_log2(v);
00401   
00402   return (i - 1);
00403 }
00404   
00405 
00406 int Triangulator::
00407 math_N(int n, int h) {
00408   int i;
00409   double v;
00410 
00411   for (i = 0, v = (int) n; i < h; i++)
00412     v = math_log2(v);
00413   
00414   return (int) ceil((double) 1.0*n/v);
00415 }
00416 
00417 
00418 /* Return a new node to be added into the query tree */
00419 int Triangulator::newnode() {
00420   int index = (int)qs.size();
00421   qs.push_back(node_t());
00422   //  cerr << "creating new node " << index << "\n";
00423   return index;
00424 }
00425 
00426 /* Return a free trapezoid */
00427 int Triangulator::newtrap() {
00428   int tr_idx = (int)tr.size();
00429   tr.push_back(trap_t());
00430   tr[tr_idx].lseg = -1;
00431   tr[tr_idx].rseg = -1;
00432   tr[tr_idx].state = ST_VALID;
00433   //  cerr << "creating new trapezoid " << tr_idx << "\n";
00434   return tr_idx;
00435 }
00436 
00437 
00438 /* Return the maximum of the two points into the yval structure */
00439 int Triangulator::_max(point_t *yval, point_t *v0, point_t *v1) {
00440   if (v0->y > v1->y + C_EPS)
00441     *yval = *v0;
00442   else if (FP_EQUAL(v0->y, v1->y))
00443     {
00444       if (v0->x > v1->x + C_EPS)
00445         *yval = *v0;
00446       else
00447         *yval = *v1;
00448     }
00449   else
00450     *yval = *v1;
00451   
00452   return 0;
00453 }
00454 
00455 
00456 /* Return the minimum of the two points into the yval structure */
00457 int Triangulator::_min(point_t *yval, point_t *v0, point_t *v1) {
00458   if (v0->y < v1->y - C_EPS)
00459     *yval = *v0;
00460   else if (FP_EQUAL(v0->y, v1->y))
00461     {
00462       if (v0->x < v1->x)
00463         *yval = *v0;
00464       else
00465         *yval = *v1;
00466     }
00467   else
00468     *yval = *v1;
00469   
00470   return 0;
00471 }
00472 
00473 
00474 int Triangulator::
00475 _greater_than(point_t *v0, point_t *v1) {
00476   if (v0->y > v1->y + C_EPS)
00477     return true;
00478   else if (v0->y < v1->y - C_EPS)
00479     return false;
00480   else
00481     return (v0->x > v1->x);
00482 }
00483 
00484 
00485 int Triangulator::
00486 _equal_to(point_t *v0, point_t *v1) {
00487   return (FP_EQUAL(v0->y, v1->y) && FP_EQUAL(v0->x, v1->x));
00488 }
00489 
00490 int Triangulator::
00491 _greater_than_equal_to(point_t *v0, point_t *v1) {
00492   if (v0->y > v1->y + C_EPS)
00493     return true;
00494   else if (v0->y < v1->y - C_EPS)
00495     return false;
00496   else
00497     return (v0->x >= v1->x);
00498 }
00499 
00500 int Triangulator::
00501 _less_than(point_t *v0, point_t *v1) {
00502   if (v0->y < v1->y - C_EPS)
00503     return true;
00504   else if (v0->y > v1->y + C_EPS)
00505     return false;
00506   else
00507     return (v0->x < v1->x);
00508 }
00509 
00510 
00511 /* Initilialise the query structure (Q) and the trapezoid table (T) 
00512  * when the first segment is added to start the trapezoidation. The
00513  * query-tree starts out with 4 trapezoids, one S-node and 2 Y-nodes
00514  *    
00515  *                4
00516  *   -----------------------------------
00517  *                \
00518  *      1          \        2
00519  *                  \
00520  *   -----------------------------------
00521  *                3
00522  */
00523 
00524 int Triangulator::
00525 init_query_structure(int segnum) {
00526   int i1, i2, i3, i4, i5, i6, i7, root;
00527   int t1, t2, t3, t4;
00528   segment_t *s = &seg[segnum];
00529 
00530   tr.clear();
00531   qs.clear();
00532 
00533   // We don't use the first elements.
00534   tr.push_back(trap_t());
00535   qs.push_back(node_t());
00536 
00537   i1 = newnode();
00538   qs[i1].nodetype = T_Y;
00539   _max(&qs[i1].yval, &s->v0, &s->v1); /* root */
00540   root = i1;
00541 
00542   i2 = newnode();
00543   qs[i1].right = i2;
00544   qs[i2].nodetype = T_SINK;
00545   qs[i2].parent = i1;
00546 
00547   i3 = newnode();
00548   qs[i1].left = i3;
00549 
00550   qs[i3].nodetype = T_Y;
00551   _min(&qs[i3].yval, &s->v0, &s->v1); /* root */
00552   qs[i3].parent = i1;
00553   
00554   i4 = newnode();
00555   qs[i3].left = i4;
00556   qs[i4].nodetype = T_SINK;
00557   qs[i4].parent = i3;
00558   
00559   i5 = newnode();
00560   qs[i3].right = i5;
00561   qs[i5].nodetype = T_X;
00562   qs[i5].segnum = segnum;
00563   qs[i5].parent = i3;
00564   
00565   i6 = newnode();
00566   qs[i5].left = i6;
00567   qs[i6].nodetype = T_SINK;
00568   qs[i6].parent = i5;
00569 
00570   i7 = newnode();
00571   qs[i5].right = i7;
00572   qs[i7].nodetype = T_SINK;
00573   qs[i7].parent = i5;
00574 
00575   t1 = newtrap();               /* middle left */
00576   t2 = newtrap();               /* middle right */
00577   t3 = newtrap();               /* bottom-most */
00578   t4 = newtrap();               /* topmost */
00579 
00580   tr[t4].lo = qs[i1].yval;
00581   tr[t2].hi = qs[i1].yval;
00582   tr[t1].hi = qs[i1].yval;
00583   tr[t3].hi = qs[i3].yval;
00584   tr[t2].lo = qs[i3].yval;
00585   tr[t1].lo = qs[i3].yval;
00586   tr[t4].hi.y = (double) (REALLY_BIG);
00587   tr[t4].hi.x = (double) (REALLY_BIG);
00588   tr[t3].lo.y = (double) -1* (REALLY_BIG);
00589   tr[t3].lo.x = (double) -1* (REALLY_BIG);
00590   tr[t2].lseg = segnum;
00591   tr[t1].rseg = segnum;
00592   tr[t2].u0 = t4;
00593   tr[t1].u0 = t4;
00594   tr[t2].d0 = t3;
00595   tr[t1].d0 = t3;
00596   tr[t3].u0 = t1;
00597   tr[t4].d0 = t1;
00598   tr[t3].u1 = t2;
00599   tr[t4].d1 = t2;
00600   
00601   tr[t1].sink = i6;
00602   tr[t2].sink = i7;
00603   tr[t3].sink = i4;
00604   tr[t4].sink = i2;
00605 
00606   tr[t2].state = ST_VALID;
00607   tr[t1].state = ST_VALID;
00608   tr[t4].state = ST_VALID;
00609   tr[t3].state = ST_VALID;
00610 
00611   qs[i2].trnum = t4;
00612   qs[i4].trnum = t3;
00613   qs[i6].trnum = t1;
00614   qs[i7].trnum = t2;
00615 
00616   s->is_inserted = true;
00617   return root;
00618 }
00619 
00620 
00621 /* Retun true if the vertex v is to the left of line segment no.
00622  * segnum. Takes care of the degenerate cases when both the vertices
00623  * have the same y--cood, etc.
00624  */
00625 
00626 int Triangulator::
00627 is_left_of(int segnum, point_t *v) {
00628   segment_t *s = &seg[segnum];
00629   double area;
00630   
00631   if (_greater_than(&s->v1, &s->v0)) /* seg. going upwards */
00632     {
00633       if (FP_EQUAL(s->v1.y, v->y))
00634         {
00635           if (v->x < s->v1.x)
00636             area = 1.0;
00637           else
00638             area = -1.0;
00639         }
00640       else if (FP_EQUAL(s->v0.y, v->y))
00641         {
00642           if (v->x < s->v0.x)
00643             area = 1.0;
00644           else
00645             area = -1.0;
00646         }
00647       else
00648         area = CROSS(s->v0, s->v1, (*v));
00649     }
00650   else                          /* v0 > v1 */
00651     {
00652       if (FP_EQUAL(s->v1.y, v->y))
00653         {
00654           if (v->x < s->v1.x)
00655             area = 1.0;
00656           else
00657             area = -1.0;
00658         }
00659       else if (FP_EQUAL(s->v0.y, v->y))
00660         {
00661           if (v->x < s->v0.x)
00662             area = 1.0;
00663           else
00664             area = -1.0;
00665         }
00666       else
00667         area = CROSS(s->v1, s->v0, (*v));
00668     }
00669   
00670   if (area > 0.0)
00671     return true;
00672   else 
00673     return false;
00674 }
00675 
00676 
00677 
00678 /* Returns true if the corresponding endpoint of the given segment is */
00679 /* already inserted into the segment tree. Use the simple test of */
00680 /* whether the segment which shares this endpoint is already inserted */
00681 
00682 int Triangulator::
00683 inserted(int segnum, int whichpt) {
00684   if (whichpt == FIRSTPT)
00685     return seg[seg[segnum].prev].is_inserted;
00686   else
00687     return seg[seg[segnum].next].is_inserted;
00688 }
00689 
00690 /* This is query routine which determines which trapezoid does the 
00691  * point v lie in. The return value is the trapezoid number. 
00692  */
00693 
00694 int Triangulator::
00695 locate_endpoint(point_t *v, point_t *vo, int r) {
00696   //  cerr << "locate_endpoint(" << v->x << " " << v->y << ", " << vo->x << " " << vo->y << ", " << r << ")\n";
00697   node_t *rptr = &qs[r];
00698   
00699   switch (rptr->nodetype)
00700     {
00701     case T_SINK:
00702       return rptr->trnum;
00703       
00704     case T_Y:
00705       if (_greater_than(v, &rptr->yval)) /* above */
00706         return locate_endpoint(v, vo, rptr->right);
00707       else if (_equal_to(v, &rptr->yval)) /* the point is already */
00708         {                                 /* inserted. */
00709           if (_greater_than(vo, &rptr->yval)) /* above */
00710             return locate_endpoint(v, vo, rptr->right);
00711           else 
00712             return locate_endpoint(v, vo, rptr->left); /* below */
00713         }
00714       else
00715         return locate_endpoint(v, vo, rptr->left); /* below */
00716 
00717     case T_X:
00718       if (_equal_to(v, &seg[rptr->segnum].v0) || 
00719                _equal_to(v, &seg[rptr->segnum].v1))
00720         {
00721           if (FP_EQUAL(v->y, vo->y)) /* horizontal segment */
00722             {
00723               if (vo->x < v->x)
00724                 return locate_endpoint(v, vo, rptr->left); /* left */
00725               else
00726                 return locate_endpoint(v, vo, rptr->right); /* right */
00727             }
00728 
00729           else if (is_left_of(rptr->segnum, vo))
00730             return locate_endpoint(v, vo, rptr->left); /* left */
00731           else
00732             return locate_endpoint(v, vo, rptr->right); /* right */
00733         }
00734       else if (is_left_of(rptr->segnum, v))
00735         return locate_endpoint(v, vo, rptr->left); /* left */
00736       else
00737         return locate_endpoint(v, vo, rptr->right); /* right */
00738 
00739     default:
00740       fprintf(stderr, "Haggu !!!!!\n");
00741       nassertr(false, -1);
00742       return -1;
00743     }
00744 }
00745 
00746 
00747 /* Thread in the segment into the existing trapezoidation. The 
00748  * limiting trapezoids are given by tfirst and tlast (which are the
00749  * trapezoids containing the two endpoints of the segment. Merges all
00750  * possible trapezoids which flank this segment and have been recently
00751  * divided because of its insertion
00752  */
00753 
00754 int Triangulator::
00755 merge_trapezoids(int segnum, int tfirst, int tlast, int side) {
00756   int t, tnext, cond;
00757   int ptnext;
00758 
00759   //  cerr << "merge_trapezoids(" << segnum << ", " << tfirst << ", " << tlast << ", " << side << ")\n";
00760 
00761   /* First merge polys on the LHS */
00762   t = tfirst;
00763   while ((t > 0) && _greater_than_equal_to(&tr[t].lo, &tr[tlast].lo))
00764     {
00765       if (side == S_LEFT)
00766         cond = ((((tnext = tr[t].d0) > 0) && (tr[tnext].rseg == segnum)) ||
00767                 (((tnext = tr[t].d1) > 0) && (tr[tnext].rseg == segnum)));
00768       else
00769         cond = ((((tnext = tr[t].d0) > 0) && (tr[tnext].lseg == segnum)) ||
00770                 (((tnext = tr[t].d1) > 0) && (tr[tnext].lseg == segnum)));
00771       
00772       if (cond)
00773         {
00774           if ((tr[t].lseg == tr[tnext].lseg) &&
00775               (tr[t].rseg == tr[tnext].rseg)) /* good neighbours */
00776             {                                 /* merge them */
00777               /* Use the upper node as the new node i.e. t */
00778 
00779               ptnext = qs[tr[tnext].sink].parent;
00780 
00781               if (qs[ptnext].left == tr[tnext].sink)
00782                 qs[ptnext].left = tr[t].sink;
00783               else
00784                 qs[ptnext].right = tr[t].sink;  /* redirect parent */
00785 
00786 
00787               /* Change the upper neighbours of the lower trapezoids */
00788 
00789               if ((tr[t].d0 = tr[tnext].d0) > 0)
00790                 if (tr[tr[t].d0].u0 == tnext)
00791                   tr[tr[t].d0].u0 = t;
00792                 else if (tr[tr[t].d0].u1 == tnext)
00793                   tr[tr[t].d0].u1 = t;
00794 
00795               if ((tr[t].d1 = tr[tnext].d1) > 0)
00796                 if (tr[tr[t].d1].u0 == tnext)
00797                   tr[tr[t].d1].u0 = t;
00798                 else if (tr[tr[t].d1].u1 == tnext)
00799                   tr[tr[t].d1].u1 = t;
00800 
00801               tr[t].lo = tr[tnext].lo;
00802               tr[tnext].state = ST_INVALID; /* invalidate the lower */
00803                                             /* trapezium */
00804             }
00805           else              /* not good neighbours */
00806             t = tnext;
00807         }
00808       else                  /* do not satisfy the outer if */
00809         t = tnext;
00810       
00811     } /* end-while */
00812        
00813   return 0;
00814 }
00815 
00816 
00817 /* Add in the new segment into the trapezoidation and update Q and T
00818  * structures. First locate the two endpoints of the segment in the
00819  * Q-structure. Then start from the topmost trapezoid and go down to
00820  * the  lower trapezoid dividing all the trapezoids in between .
00821  */
00822 
00823 int Triangulator::
00824 add_segment(int segnum) {
00825   //  cerr << "add_segment(" << segnum << ")\n";
00826 
00827   segment_t s;
00828   //  segment_t *so = &seg[segnum];
00829   int tu, tl, sk, tfirst, tlast; //, tnext;
00830   int tfirstr = 0, tlastr = 0, tfirstl = 0, tlastl = 0;
00831   int i1, i2, t, tn; // t1, t2,
00832   point_t tpt;
00833   int tritop = 0, tribot = 0, is_swapped = 0;
00834   int tmptriseg;
00835 
00836   s = seg[segnum];
00837   if (_greater_than(&s.v1, &s.v0)) /* Get higher vertex in v0 */
00838     {
00839       int tmp;
00840       tpt = s.v0;
00841       s.v0 = s.v1;
00842       s.v1 = tpt;
00843       tmp = s.root0;
00844       s.root0 = s.root1;
00845       s.root1 = tmp;
00846       is_swapped = true;
00847     }
00848 
00849   if ((is_swapped) ? !inserted(segnum, LASTPT) :
00850        !inserted(segnum, FIRSTPT))     /* insert v0 in the tree */
00851     {
00852       int tmp_d;
00853 
00854       tu = locate_endpoint(&s.v0, &s.v1, s.root0);
00855       tl = newtrap();           /* tl is the new lower trapezoid */
00856       tr[tl].state = ST_VALID;
00857       tr[tl] = tr[tu];
00858       tr[tl].hi.y = s.v0.y;
00859       tr[tu].lo.y = s.v0.y;
00860       tr[tl].hi.x = s.v0.x;
00861       tr[tu].lo.x = s.v0.x;
00862       tr[tu].d0 = tl;      
00863       tr[tu].d1 = 0;
00864       tr[tl].u0 = tu;
00865       tr[tl].u1 = 0;
00866 
00867       if (((tmp_d = tr[tl].d0) > 0) && (tr[tmp_d].u0 == tu))
00868         tr[tmp_d].u0 = tl;
00869       if (((tmp_d = tr[tl].d0) > 0) && (tr[tmp_d].u1 == tu))
00870         tr[tmp_d].u1 = tl;
00871 
00872       if (((tmp_d = tr[tl].d1) > 0) && (tr[tmp_d].u0 == tu))
00873         tr[tmp_d].u0 = tl;
00874       if (((tmp_d = tr[tl].d1) > 0) && (tr[tmp_d].u1 == tu))
00875         tr[tmp_d].u1 = tl;
00876 
00877       /* Now update the query structure and obtain the sinks for the */
00878       /* two trapezoids */ 
00879       
00880       i1 = newnode();           /* Upper trapezoid sink */
00881       i2 = newnode();           /* Lower trapezoid sink */
00882       sk = tr[tu].sink;
00883       
00884       qs[sk].nodetype = T_Y;
00885       qs[sk].yval = s.v0;
00886       qs[sk].segnum = segnum;   /* not really reqd ... maybe later */
00887       qs[sk].left = i2;
00888       qs[sk].right = i1;
00889 
00890       qs[i1].nodetype = T_SINK;
00891       qs[i1].trnum = tu;
00892       qs[i1].parent = sk;
00893 
00894       qs[i2].nodetype = T_SINK;
00895       qs[i2].trnum = tl;
00896       qs[i2].parent = sk;
00897 
00898       tr[tu].sink = i1;
00899       tr[tl].sink = i2;
00900       tfirst = tl;
00901     }
00902   else                          /* v0 already present */
00903     {       /* Get the topmost intersecting trapezoid */
00904       tfirst = locate_endpoint(&s.v0, &s.v1, s.root0);
00905       tritop = 1;
00906     }
00907 
00908 
00909   if ((is_swapped) ? !inserted(segnum, FIRSTPT) :
00910        !inserted(segnum, LASTPT))     /* insert v1 in the tree */
00911     {
00912       int tmp_d;
00913 
00914       tu = locate_endpoint(&s.v1, &s.v0, s.root1);
00915 
00916       tl = newtrap();           /* tl is the new lower trapezoid */
00917       tr[tl].state = ST_VALID;
00918       tr[tl] = tr[tu];
00919       tr[tl].hi.y = s.v1.y;
00920       tr[tu].lo.y = s.v1.y;
00921       tr[tl].hi.x = s.v1.x;
00922       tr[tu].lo.x = s.v1.x;
00923       tr[tu].d0 = tl;      
00924       tr[tu].d1 = 0;
00925       tr[tl].u0 = tu;
00926       tr[tl].u1 = 0;
00927 
00928       if (((tmp_d = tr[tl].d0) > 0) && (tr[tmp_d].u0 == tu))
00929         tr[tmp_d].u0 = tl;
00930       if (((tmp_d = tr[tl].d0) > 0) && (tr[tmp_d].u1 == tu))
00931         tr[tmp_d].u1 = tl;
00932 
00933       if (((tmp_d = tr[tl].d1) > 0) && (tr[tmp_d].u0 == tu))
00934         tr[tmp_d].u0 = tl;
00935       if (((tmp_d = tr[tl].d1) > 0) && (tr[tmp_d].u1 == tu))
00936         tr[tmp_d].u1 = tl;
00937       
00938       /* Now update the query structure and obtain the sinks for the */
00939       /* two trapezoids */ 
00940       
00941       i1 = newnode();           /* Upper trapezoid sink */
00942       i2 = newnode();           /* Lower trapezoid sink */
00943       sk = tr[tu].sink;
00944       
00945       qs[sk].nodetype = T_Y;
00946       qs[sk].yval = s.v1;
00947       qs[sk].segnum = segnum;   /* not really reqd ... maybe later */
00948       qs[sk].left = i2;
00949       qs[sk].right = i1;
00950 
00951       qs[i1].nodetype = T_SINK;
00952       qs[i1].trnum = tu;
00953       qs[i1].parent = sk;
00954 
00955       qs[i2].nodetype = T_SINK;
00956       qs[i2].trnum = tl;
00957       qs[i2].parent = sk;
00958 
00959       tr[tu].sink = i1;
00960       tr[tl].sink = i2;
00961       tlast = tu;
00962     }
00963   else                          /* v1 already present */
00964     {       /* Get the lowermost intersecting trapezoid */
00965       tlast = locate_endpoint(&s.v1, &s.v0, s.root1);
00966       tribot = 1;
00967     }
00968   
00969   /* Thread the segment into the query tree creating a new X-node */
00970   /* First, split all the trapezoids which are intersected by s into */
00971   /* two */
00972 
00973   t = tfirst;                   /* topmost trapezoid */
00974   
00975   while ((t > 0) && 
00976          _greater_than_equal_to(&tr[t].lo, &tr[tlast].lo))
00977                                 /* traverse from top to bot */
00978     {
00979       int t_sav, tn_sav;
00980       sk = tr[t].sink;
00981       i1 = newnode();           /* left trapezoid sink */
00982       i2 = newnode();           /* right trapezoid sink */
00983       
00984       qs[sk].nodetype = T_X;
00985       qs[sk].segnum = segnum;
00986       qs[sk].left = i1;
00987       qs[sk].right = i2;
00988 
00989       qs[i1].nodetype = T_SINK; /* left trapezoid (use existing one) */
00990       qs[i1].trnum = t;
00991       qs[i1].parent = sk;
00992 
00993       qs[i2].nodetype = T_SINK; /* right trapezoid (allocate new) */
00994       tn = newtrap();
00995       qs[i2].trnum = tn;
00996       tr[tn].state = ST_VALID;
00997       qs[i2].parent = sk;
00998 
00999       if (t == tfirst)
01000         tfirstr = tn;
01001       if (_equal_to(&tr[t].lo, &tr[tlast].lo))
01002         tlastr = tn;
01003 
01004       tr[tn] = tr[t];
01005       tr[t].sink = i1;
01006       tr[tn].sink = i2;
01007       t_sav = t;
01008       tn_sav = tn;
01009 
01010       /* error */
01011 
01012       if ((tr[t].d0 <= 0) && (tr[t].d1 <= 0)) /* case cannot arise */
01013         {
01014           fprintf(stderr, "add_segment: error\n");
01015           break;
01016         }
01017       
01018       /* only one trapezoid below. partition t into two and make the */
01019       /* two resulting trapezoids t and tn as the upper neighbours of */
01020       /* the sole lower trapezoid */
01021       
01022       else if ((tr[t].d0 > 0) && (tr[t].d1 <= 0))
01023         {                       /* Only one trapezoid below */
01024           if ((tr[t].u0 > 0) && (tr[t].u1 > 0))
01025             {                   /* continuation of a chain from abv. */
01026               if (tr[t].usave > 0) /* three upper neighbours */
01027                 {
01028                   if (tr[t].uside == S_LEFT)
01029                     {
01030                       tr[tn].u0 = tr[t].u1;
01031                       tr[t].u1 = -1;
01032                       tr[tn].u1 = tr[t].usave;
01033 
01034                       tr[tr[t].u0].d0 = t;
01035                       tr[tr[tn].u0].d0 = tn;
01036                       tr[tr[tn].u1].d0 = tn;
01037                     }
01038                   else          /* intersects in the right */
01039                     {
01040                       tr[tn].u1 = -1;
01041                       tr[tn].u0 = tr[t].u1;
01042                       tr[t].u1 = tr[t].u0;
01043                       tr[t].u0 = tr[t].usave;
01044 
01045                       tr[tr[t].u0].d0 = t;
01046                       tr[tr[t].u1].d0 = t;
01047                       tr[tr[tn].u0].d0 = tn;
01048                     }
01049                   tr[tn].usave = 0;
01050                   tr[t].usave = 0;
01051                 }
01052               else              /* No usave.... simple case */
01053                 {
01054                   tr[tn].u0 = tr[t].u1;
01055                   tr[tn].u1 = -1;
01056                   tr[t].u1 = -1;
01057                   tr[tr[tn].u0].d0 = tn;
01058                 }
01059             }
01060           else 
01061             {                   /* fresh seg. or upward cusp */
01062               int tmp_u = tr[t].u0;
01063               int td0, td1;
01064               if (((td0 = tr[tmp_u].d0) > 0) && 
01065                   ((td1 = tr[tmp_u].d1) > 0))
01066                 {               /* upward cusp */
01067                   if ((tr[td0].rseg > 0) &&
01068                       !is_left_of(tr[td0].rseg, &s.v1))
01069                     {
01070                       tr[tn].u1 = -1;
01071                       tr[t].u1 = -1;
01072                       tr[t].u0 = -1;
01073                       tr[tr[tn].u0].d1 = tn;
01074                     }
01075                   else          /* cusp going leftwards */
01076                     { 
01077                       tr[t].u1 = -1;
01078                       tr[tn].u1 = -1;
01079                       tr[tn].u0 = -1;
01080                       tr[tr[t].u0].d0 = t;
01081                     }
01082                 }
01083               else              /* fresh segment */
01084                 {
01085                   tr[tr[t].u0].d0 = t;
01086                   tr[tr[t].u0].d1 = tn;
01087                 }
01088             }
01089 
01090           if (FP_EQUAL(tr[t].lo.y, tr[tlast].lo.y) && 
01091               FP_EQUAL(tr[t].lo.x, tr[tlast].lo.x) && tribot)
01092             {           /* bottom forms a triangle */
01093 
01094               if (is_swapped)
01095                 tmptriseg = seg[segnum].prev;
01096               else
01097                 tmptriseg = seg[segnum].next;
01098 
01099               if ((tmptriseg > 0) && is_left_of(tmptriseg, &s.v0))
01100                 {
01101                                 /* L-R downward cusp */
01102                   tr[tr[t].d0].u0 = t;
01103                   tr[tn].d1 = -1;
01104                   tr[tn].d0 = -1;
01105                 }
01106               else
01107                 {
01108                                 /* R-L downward cusp */
01109                   tr[tr[tn].d0].u1 = tn;
01110                   tr[t].d1 = -1;
01111                   tr[t].d0 = -1;
01112                 }
01113             }
01114           else
01115             {
01116               if ((tr[tr[t].d0].u0 > 0) && (tr[tr[t].d0].u1 > 0))
01117                 {
01118                   if (tr[tr[t].d0].u0 == t) /* passes thru LHS */
01119                     {
01120                       tr[tr[t].d0].usave = tr[tr[t].d0].u1;
01121                       tr[tr[t].d0].uside = S_LEFT;
01122                     }
01123                   else
01124                     {
01125                       tr[tr[t].d0].usave = tr[tr[t].d0].u0;
01126                       tr[tr[t].d0].uside = S_RIGHT;
01127                     }
01128                 }
01129               tr[tr[t].d0].u0 = t;
01130               tr[tr[t].d0].u1 = tn;
01131             }
01132 
01133           t = tr[t].d0;
01134         }
01135 
01136 
01137       else if ((tr[t].d0 <= 0) && (tr[t].d1 > 0))
01138         {                       /* Only one trapezoid below */
01139           if ((tr[t].u0 > 0) && (tr[t].u1 > 0))
01140             {                   /* continuation of a chain from abv. */
01141               if (tr[t].usave > 0) /* three upper neighbours */
01142                 {
01143                   if (tr[t].uside == S_LEFT)
01144                     {
01145                       tr[tn].u0 = tr[t].u1;
01146                       tr[t].u1 = -1;
01147                       tr[tn].u1 = tr[t].usave;
01148 
01149                       tr[tr[t].u0].d0 = t;
01150                       tr[tr[tn].u0].d0 = tn;
01151                       tr[tr[tn].u1].d0 = tn;
01152                     }
01153                   else          /* intersects in the right */
01154                     {
01155                       tr[tn].u1 = -1;
01156                       tr[tn].u0 = tr[t].u1;
01157                       tr[t].u1 = tr[t].u0;
01158                       tr[t].u0 = tr[t].usave;
01159 
01160                       tr[tr[t].u0].d0 = t;
01161                       tr[tr[t].u1].d0 = t;
01162                       tr[tr[tn].u0].d0 = tn;
01163                     }
01164 
01165                   tr[tn].usave = 0;
01166                   tr[t].usave = 0;
01167                 }
01168               else              /* No usave.... simple case */
01169                 {
01170                   tr[tn].u0 = tr[t].u1;
01171                   tr[tn].u1 = -1;
01172                   tr[t].u1 = -1;
01173                   tr[tr[tn].u0].d0 = tn;
01174                 }
01175             }
01176           else 
01177             {                   /* fresh seg. or upward cusp */
01178               int tmp_u = tr[t].u0;
01179               int td0, td1;
01180               if (((td0 = tr[tmp_u].d0) > 0) && 
01181                   ((td1 = tr[tmp_u].d1) > 0))
01182                 {               /* upward cusp */
01183                   if ((tr[td0].rseg > 0) &&
01184                       !is_left_of(tr[td0].rseg, &s.v1))
01185                     {
01186                       tr[tn].u1 = -1;
01187                       tr[t].u1 = -1;
01188                       tr[t].u0 = -1;
01189                       tr[tr[tn].u0].d1 = tn;
01190                     }
01191                   else 
01192                     {
01193                       tr[t].u1 = -1;
01194                       tr[tn].u1 = -1;
01195                       tr[tn].u0 = -1;
01196                       tr[tr[t].u0].d0 = t;
01197                     }
01198                 }
01199               else              /* fresh segment */
01200                 {
01201                   tr[tr[t].u0].d0 = t;
01202                   tr[tr[t].u0].d1 = tn;
01203                 }
01204             }
01205 
01206           if (FP_EQUAL(tr[t].lo.y, tr[tlast].lo.y) && 
01207               FP_EQUAL(tr[t].lo.x, tr[tlast].lo.x) && tribot)
01208             {           /* bottom forms a triangle */
01209               if (is_swapped)
01210                 tmptriseg = seg[segnum].prev;
01211               else
01212                 tmptriseg = seg[segnum].next;
01213 
01214               if ((tmptriseg > 0) && is_left_of(tmptriseg, &s.v0))
01215                 {
01216                   /* L-R downward cusp */
01217                   tr[tr[t].d1].u0 = t;
01218                   tr[tn].d1 = -1;
01219                   tr[tn].d0 = -1;
01220                 }
01221               else
01222                 {
01223                   /* R-L downward cusp */
01224                   tr[tr[tn].d1].u1 = tn;
01225                   tr[t].d1 = -1;
01226                   tr[t].d0 = -1;
01227                 }
01228             }
01229           else
01230             {
01231               if ((tr[tr[t].d1].u0 > 0) && (tr[tr[t].d1].u1 > 0))
01232                 {
01233                   if (tr[tr[t].d1].u0 == t) /* passes thru LHS */
01234                     {
01235                       tr[tr[t].d1].usave = tr[tr[t].d1].u1;
01236                       tr[tr[t].d1].uside = S_LEFT;
01237                     }
01238                   else
01239                     {
01240                       tr[tr[t].d1].usave = tr[tr[t].d1].u0;
01241                       tr[tr[t].d1].uside = S_RIGHT;
01242                     }
01243                 }
01244               tr[tr[t].d1].u0 = t;
01245               tr[tr[t].d1].u1 = tn;
01246             }
01247 
01248           t = tr[t].d1;
01249         }
01250 
01251       /* two trapezoids below. Find out which one is intersected by */
01252       /* this segment and proceed down that one */
01253       
01254       else
01255         {
01256           //      int tmpseg = tr[tr[t].d0].rseg;
01257           double y0, yt;
01258           point_t tmppt;
01259           int tnext, i_d0, i_d1;
01260 
01261           i_d1 = false;
01262           i_d0 = false;
01263           if (FP_EQUAL(tr[t].lo.y, s.v0.y))
01264             {
01265               if (tr[t].lo.x > s.v0.x)
01266                 i_d0 = true;
01267               else
01268                 i_d1 = true;
01269             }
01270           else
01271             {
01272               y0 = tr[t].lo.y;
01273               tmppt.y = y0;
01274               yt = (y0 - s.v0.y)/(s.v1.y - s.v0.y);
01275               tmppt.x = s.v0.x + yt * (s.v1.x - s.v0.x);
01276 
01277               if (_less_than(&tmppt, &tr[t].lo))
01278                 i_d0 = true;
01279               else
01280                 i_d1 = true;
01281             }
01282 
01283           /* check continuity from the top so that the lower-neighbour */
01284           /* values are properly filled for the upper trapezoid */
01285 
01286           if ((tr[t].u0 > 0) && (tr[t].u1 > 0))
01287             {                   /* continuation of a chain from abv. */
01288               if (tr[t].usave > 0) /* three upper neighbours */
01289                 {
01290                   if (tr[t].uside == S_LEFT)
01291                     {
01292                       tr[tn].u0 = tr[t].u1;
01293                       tr[t].u1 = -1;
01294                       tr[tn].u1 = tr[t].usave;
01295 
01296                       tr[tr[t].u0].d0 = t;
01297                       tr[tr[tn].u0].d0 = tn;
01298                       tr[tr[tn].u1].d0 = tn;
01299                     }
01300                   else          /* intersects in the right */
01301                     {
01302                       tr[tn].u1 = -1;
01303                       tr[tn].u0 = tr[t].u1;
01304                       tr[t].u1 = tr[t].u0;
01305                       tr[t].u0 = tr[t].usave;
01306 
01307                       tr[tr[t].u0].d0 = t;
01308                       tr[tr[t].u1].d0 = t;
01309                       tr[tr[tn].u0].d0 = tn;
01310                     }
01311 
01312                   tr[tn].usave = 0;
01313                   tr[t].usave = 0;
01314                 }
01315               else              /* No usave.... simple case */
01316                 {
01317                   tr[tn].u0 = tr[t].u1;
01318                   tr[tn].u1 = -1;
01319                   tr[t].u1 = -1;
01320                   tr[tr[tn].u0].d0 = tn;
01321                 }
01322             }
01323           else 
01324             {                   /* fresh seg. or upward cusp */
01325               int tmp_u = tr[t].u0;
01326               int td0, td1;
01327               if (((td0 = tr[tmp_u].d0) > 0) && 
01328                   ((td1 = tr[tmp_u].d1) > 0))
01329                 {               /* upward cusp */
01330                   if ((tr[td0].rseg > 0) &&
01331                       !is_left_of(tr[td0].rseg, &s.v1))
01332                     {
01333                       tr[tn].u1 = -1;
01334                       tr[t].u1 = -1;
01335                       tr[t].u0 = -1;
01336                       tr[tr[tn].u0].d1 = tn;
01337                     }
01338                   else 
01339                     {
01340                       tr[t].u1 = -1;
01341                       tr[tn].u1 = -1;
01342                       tr[tn].u0 = -1;
01343                       tr[tr[t].u0].d0 = t;
01344                     }
01345                 }
01346               else              /* fresh segment */
01347                 {
01348                   tr[tr[t].u0].d0 = t;
01349                   tr[tr[t].u0].d1 = tn;
01350                 }
01351             }
01352 
01353           if (FP_EQUAL(tr[t].lo.y, tr[tlast].lo.y) && 
01354               FP_EQUAL(tr[t].lo.x, tr[tlast].lo.x) && tribot)
01355             {
01356               /* this case arises only at the lowest trapezoid.. i.e.
01357                  tlast, if the lower endpoint of the segment is
01358                  already inserted in the structure */
01359 
01360               tr[tr[t].d0].u0 = t;
01361               tr[tr[t].d0].u1 = -1;
01362               tr[tr[t].d1].u0 = tn;
01363               tr[tr[t].d1].u1 = -1;
01364 
01365               tr[tn].d0 = tr[t].d1;
01366               tr[tn].d1 = -1;
01367               tr[t].d1 = -1;
01368 
01369               tnext = tr[t].d1;
01370             }
01371           else if (i_d0)
01372                                 /* intersecting d0 */
01373             {
01374               tr[tr[t].d0].u0 = t;
01375               tr[tr[t].d0].u1 = tn;
01376               tr[tr[t].d1].u0 = tn;
01377               tr[tr[t].d1].u1 = -1;
01378 
01379               /* new code to determine the bottom neighbours of the */
01380               /* newly partitioned trapezoid */
01381 
01382               tr[t].d1 = -1;
01383 
01384               tnext = tr[t].d0;
01385             }
01386           else                  /* intersecting d1 */
01387             {
01388               tr[tr[t].d0].u0 = t;
01389               tr[tr[t].d0].u1 = -1;
01390               tr[tr[t].d1].u0 = t;
01391               tr[tr[t].d1].u1 = tn;
01392 
01393               /* new code to determine the bottom neighbours of the */
01394               /* newly partitioned trapezoid */
01395 
01396               tr[tn].d0 = tr[t].d1;
01397               tr[tn].d1 = -1;
01398 
01399               tnext = tr[t].d1;
01400             }
01401 
01402           t = tnext;
01403         }
01404       
01405       tr[tn_sav].lseg  = segnum;
01406       tr[t_sav].rseg = segnum;
01407     } /* end-while */
01408   
01409   /* Now combine those trapezoids which share common segments. We can */
01410   /* use the pointers to the parent to connect these together. This */
01411   /* works only because all these new trapezoids have been formed */
01412   /* due to splitting by the segment, and hence have only one parent */
01413 
01414   tfirstl = tfirst; 
01415   tlastl = tlast;
01416   merge_trapezoids(segnum, tfirstl, tlastl, S_LEFT);
01417   merge_trapezoids(segnum, tfirstr, tlastr, S_RIGHT);
01418 
01419   seg[segnum].is_inserted = true;
01420   return 0;
01421 }
01422 
01423 
01424 /* Update the roots stored for each of the endpoints of the segment.
01425  * This is done to speed up the location-query for the endpoint when
01426  * the segment is inserted into the trapezoidation subsequently
01427  */
01428 int Triangulator::
01429 find_new_roots(int segnum) {
01430   //  cerr << "find_new_roots(" << segnum << ")\n";
01431   segment_t *s = &seg[segnum];
01432   
01433   if (s->is_inserted)
01434     return 0;
01435 
01436   s->root0 = locate_endpoint(&s->v0, &s->v1, s->root0);
01437   s->root0 = tr[s->root0].sink;
01438 
01439   s->root1 = locate_endpoint(&s->v1, &s->v0, s->root1);
01440   s->root1 = tr[s->root1].sink;  
01441   return 0;
01442 }
01443 
01444 
01445 /* Main routine to perform trapezoidation */
01446 int Triangulator::
01447 construct_trapezoids(int nseg) {
01448   //  cerr << "construct_trapezoids(" << nseg << ")\n";
01449   int i;
01450   int root, h;
01451   
01452   /* Add the first segment and get the query structure and trapezoid */
01453   /* list initialised */
01454 
01455   root = init_query_structure(choose_segment());
01456 
01457   for (i = 1; i <= nseg; i++) {
01458     seg[i].root1 = root;
01459     seg[i].root0 = root;
01460   }
01461   
01462   for (h = 1; h <= math_logstar_n(nseg); h++)
01463     {
01464       for (i = math_N(nseg, h -1) + 1; i <= math_N(nseg, h); i++) {
01465         add_segment(choose_segment());
01466       }
01467       
01468       /* Find a new root for each of the segment endpoints */
01469       for (i = 1; i <= nseg; i++)
01470         find_new_roots(i);
01471     }
01472   
01473   for (i = math_N(nseg, math_logstar_n(nseg)) + 1; i <= nseg; i++)
01474     add_segment(choose_segment());
01475 
01476   return 0;
01477 }
01478 
01479 
01480 
01481 /* Function returns true if the trapezoid lies inside the polygon */
01482 int Triangulator::
01483 inside_polygon(trap_t *t) {
01484   int rseg = t->rseg;
01485 
01486   if (t->state == ST_INVALID)
01487     return 0;
01488 
01489   if ((t->lseg <= 0) || (t->rseg <= 0))
01490     return 0;
01491   
01492   if (((t->u0 <= 0) && (t->u1 <= 0)) || 
01493       ((t->d0 <= 0) && (t->d1 <= 0))) /* triangle */
01494     return (_greater_than(&seg[rseg].v1, &seg[rseg].v0));
01495   
01496   return 0;
01497 }
01498 
01499 
01500 /* return a new mon structure from the table */
01501 int Triangulator::
01502 newmon() {
01503   int index = (int)mon.size();
01504   mon.push_back(0);
01505   //  cerr << "newmon " << index << "\n";
01506   return index;
01507 }
01508 
01509 
01510 /* return a new chain element from the table */
01511 int Triangulator::
01512 new_chain_element() {
01513   int index = (int)mchain.size();
01514   mchain.push_back(monchain_t());
01515   //  cerr << "new_chain_element " << index << "\n";
01516   return index;
01517 }
01518 
01519 
01520 double Triangulator::
01521 get_angle(point_t *vp0, point_t *vpnext, point_t *vp1) {
01522   point_t v0, v1;
01523   
01524   v0.x = vpnext->x - vp0->x;
01525   v0.y = vpnext->y - vp0->y;
01526 
01527   v1.x = vp1->x - vp0->x;
01528   v1.y = vp1->y - vp0->y;
01529 
01530   if (CROSS_SINE(v0, v1) >= 0)  /* sine is positive */
01531     return DOT(v0, v1)/LENGTH(v0)/LENGTH(v1);
01532   else
01533     return (-1.0 * DOT(v0, v1)/LENGTH(v0)/LENGTH(v1) - 2);
01534 }
01535 
01536 
01537 /* (v0, v1) is the new diagonal to be added to the polygon. Find which */
01538 /* chain to use and return the positions of v0 and v1 in p and q */ 
01539 int Triangulator::
01540 get_vertex_positions(int v0, int v1, int *ip, int *iq) {
01541   vertexchain_t *vp0, *vp1;
01542   int i;
01543   double angle, temp;
01544   int tp = 0, tq = 0;
01545 
01546   vp0 = &vert[v0];
01547   vp1 = &vert[v1];
01548   
01549   /* p is identified as follows. Scan from (v0, v1) rightwards till */
01550   /* you hit the first segment starting from v0. That chain is the */
01551   /* chain of our interest */
01552   
01553   angle = -4.0;
01554   for (i = 0; i < 4; i++)
01555     {
01556       if (vp0->vnext[i] <= 0)
01557         continue;
01558       if ((temp = get_angle(&vp0->pt, &(vert[vp0->vnext[i]].pt), 
01559                             &vp1->pt)) > angle)
01560         {
01561           angle = temp;
01562           tp = i;
01563         }
01564     }
01565 
01566   *ip = tp;
01567 
01568   /* Do similar actions for q */
01569 
01570   angle = -4.0;
01571   for (i = 0; i < 4; i++)
01572     {
01573       if (vp1->vnext[i] <= 0)
01574         continue;      
01575       if ((temp = get_angle(&vp1->pt, &(vert[vp1->vnext[i]].pt), 
01576                             &vp0->pt)) > angle)
01577         {
01578           angle = temp;
01579           tq = i;
01580         }
01581     }
01582 
01583   *iq = tq;
01584 
01585   return 0;
01586 }
01587 
01588   
01589 /* v0 and v1 are specified in anti-clockwise order with respect to 
01590  * the current monotone polygon mcur. Split the current polygon into 
01591  * two polygons using the diagonal (v0, v1) 
01592  */
01593 int Triangulator::
01594 make_new_monotone_poly(int mcur, int v0, int v1) {
01595   int p, q, ip, iq;
01596   int mnew = newmon();
01597   int i, j, nf0, nf1;
01598   vertexchain_t *vp0, *vp1;
01599 
01600   if (v0 <= 0 || v1 <= 0) {
01601     return -1;
01602   }
01603   
01604   vp0 = &vert[v0];
01605   vp1 = &vert[v1];
01606 
01607   get_vertex_positions(v0, v1, &ip, &iq);
01608 
01609   p = vp0->vpos[ip];
01610   q = vp1->vpos[iq];
01611 
01612   /* At this stage, we have got the positions of v0 and v1 in the */
01613   /* desired chain. Now modify the linked lists */
01614 
01615   i = new_chain_element();      /* for the new list */
01616   j = new_chain_element();
01617 
01618   mchain[i].vnum = v0;
01619   mchain[j].vnum = v1;
01620 
01621   mchain[i].next = mchain[p].next;
01622   mchain[mchain[p].next].prev = i;
01623   mchain[i].prev = j;
01624   mchain[j].next = i;
01625   mchain[j].prev = mchain[q].prev;
01626   mchain[mchain[q].prev].next = j;
01627 
01628   mchain[p].next = q;
01629   mchain[q].prev = p;
01630 
01631   nf0 = vp0->nextfree;
01632   nf1 = vp1->nextfree;
01633 
01634   vp0->vnext[ip] = v1;
01635 
01636   vp0->vpos[nf0] = i;
01637   vp0->vnext[nf0] = mchain[mchain[i].next].vnum;
01638   vp1->vpos[nf1] = j;
01639   vp1->vnext[nf1] = v0;
01640 
01641   vp0->nextfree++;
01642   vp1->nextfree++;
01643 
01644 #ifdef DEBUG
01645   fprintf(stderr, "make_poly: mcur = %d, (v0, v1) = (%d, %d)\n", 
01646           mcur, v0, v1);
01647   fprintf(stderr, "next posns = (p, q) = (%d, %d)\n", p, q);
01648 #endif
01649 
01650   mon[mcur] = p;
01651   mon[mnew] = i;
01652   return mnew;
01653 }
01654 
01655 /* Main routine to get monotone polygons from the trapezoidation of 
01656  * the polygon.
01657  */
01658 
01659 int Triangulator::
01660 monotonate_trapezoids(int n) {
01661   int i;
01662   int tr_start;
01663 
01664   vert.clear();
01665   visited.clear();
01666   mchain.clear();
01667   mon.clear();
01668 
01669   vert.insert(vert.begin(), n + 1, vertexchain_t());
01670   mchain.insert(mchain.begin(), n + 1, monchain_t());
01671 
01672   visited.insert(visited.begin(), tr.size(), 0);
01673   
01674   /* First locate a trapezoid which lies inside the polygon */
01675   /* and which is triangular */
01676   for (i = 1; i < (int)tr.size(); i++)
01677     if (inside_polygon(&tr[i]))
01678       break;
01679   if (i >= (int)tr.size()) {
01680     // No valid trapezoids.
01681     return 0;
01682   }
01683   //  printf("start = %d\n", i);
01684   tr_start = i;
01685   
01686   /* Initialise the mon data-structure and start spanning all the */
01687   /* trapezoids within the polygon */
01688 
01689 #if 0
01690   for (i = 1; i <= n; i++)
01691     {
01692       mchain[i].prev = i - 1;
01693       mchain[i].next = i + 1;
01694       mchain[i].vnum = i;
01695       vert[i].pt = seg[i].v0;
01696       vert[i].vnext[0] = i + 1; /* next vertex */
01697       vert[i].vpos[0] = i;      /* locn. of next vertex */
01698       vert[i].nextfree = 1;
01699       vert[i].user_i = seg[i].v0_i;
01700     }
01701   mchain[1].prev = n;
01702   mchain[n].next = 1;
01703   vert[n].vnext[0] = 1;
01704   vert[n].vpos[0] = n;
01705   mon.push_back(1); /* position of any vertex in the first */
01706                     /* chain  */
01707 
01708 #else
01709 
01710   for (i = 1; i <= n; i++)
01711     {
01712       mchain[i].prev = seg[i].prev;
01713       mchain[i].next = seg[i].next;
01714       mchain[i].vnum = i;
01715       vert[i].pt = seg[i].v0;
01716       vert[i].vnext[0] = seg[i].next; /* next vertex */
01717       vert[i].vpos[0] = i;      /* locn. of next vertex */
01718       vert[i].nextfree = 1;
01719       vert[i].user_i = seg[i].v0_i;
01720     }
01721 
01722   mon.push_back(1); /* position of any vertex in the first */
01723                     /* chain  */
01724 
01725 #endif
01726   
01727   /* traverse the polygon */
01728   if (tr[tr_start].u0 > 0)
01729     traverse_polygon(0, tr_start, tr[tr_start].u0, TR_FROM_UP);
01730   else if (tr[tr_start].d0 > 0)
01731     traverse_polygon(0, tr_start, tr[tr_start].d0, TR_FROM_DN);
01732   
01733   /* return the number of polygons created */
01734   return newmon();
01735 }
01736 
01737 
01738 /* recursively visit all the trapezoids */
01739 int Triangulator::
01740 traverse_polygon(int mcur, int trnum, int from, int dir) {
01741   //  printf("traverse_polygon(%d, %d, %d, %d)\n", mcur, trnum, from, dir);
01742 
01743   if (mcur < 0 || trnum <= 0)
01744     return 0;
01745 
01746   if (visited[trnum])
01747     return 0;
01748 
01749   trap_t *t = &tr[trnum];
01750   //  int howsplit;
01751   int mnew;
01752   int v0, v1;  //, v0next, v1next;
01753   int retval = 0;  //, tmp;
01754   int do_switch = false;
01755 
01756   //  printf("visited size = %d, visited[trnum] = %d\n", visited.size(), visited[trnum]);
01757 
01758   visited[trnum] = true;
01759   
01760   /* We have much more information available here. */
01761   /* rseg: goes upwards   */
01762   /* lseg: goes downwards */
01763 
01764   /* Initially assume that dir = TR_FROM_DN (from the left) */
01765   /* Switch v0 and v1 if necessary afterwards */
01766 
01767 
01768   /* special cases for triangles with cusps at the opposite ends. */
01769   /* take care of this first */
01770   if ((t->u0 <= 0) && (t->u1 <= 0))
01771     {
01772       if ((t->d0 > 0) && (t->d1 > 0)) /* downward opening triangle */
01773         {
01774           v0 = tr[t->d1].lseg;
01775           v1 = t->lseg;
01776           if (from == t->d1)
01777             {
01778               do_switch = true;
01779               mnew = make_new_monotone_poly(mcur, v1, v0);
01780               traverse_polygon(mcur, t->d1, trnum, TR_FROM_UP);
01781               traverse_polygon(mnew, t->d0, trnum, TR_FROM_UP);
01782             }
01783           else
01784             {
01785               mnew = make_new_monotone_poly(mcur, v0, v1);
01786               traverse_polygon(mcur, t->d0, trnum, TR_FROM_UP);
01787               traverse_polygon(mnew, t->d1, trnum, TR_FROM_UP);
01788             }
01789         }
01790       else
01791         {
01792           retval = SP_NOSPLIT;  /* Just traverse all neighbours */
01793           traverse_polygon(mcur, t->u0, trnum, TR_FROM_DN);
01794           traverse_polygon(mcur, t->u1, trnum, TR_FROM_DN);
01795           traverse_polygon(mcur, t->d0, trnum, TR_FROM_UP);
01796           traverse_polygon(mcur, t->d1, trnum, TR_FROM_UP);
01797         }
01798     }
01799   
01800   else if ((t->d0 <= 0) && (t->d1 <= 0))
01801     {
01802       if ((t->u0 > 0) && (t->u1 > 0)) /* upward opening triangle */
01803         {
01804           v0 = t->rseg;
01805           v1 = tr[t->u0].rseg;
01806           if (from == t->u1)
01807             {
01808               do_switch = true;
01809               mnew = make_new_monotone_poly(mcur, v1, v0);
01810               traverse_polygon(mcur, t->u1, trnum, TR_FROM_DN);
01811               traverse_polygon(mnew, t->u0, trnum, TR_FROM_DN);
01812             }
01813           else
01814             {
01815               mnew = make_new_monotone_poly(mcur, v0, v1);
01816               traverse_polygon(mcur, t->u0, trnum, TR_FROM_DN);
01817               traverse_polygon(mnew, t->u1, trnum, TR_FROM_DN);
01818             }
01819         }
01820       else
01821         {
01822           retval = SP_NOSPLIT;  /* Just traverse all neighbours */
01823           traverse_polygon(mcur, t->u0, trnum, TR_FROM_DN);
01824           traverse_polygon(mcur, t->u1, trnum, TR_FROM_DN);
01825           traverse_polygon(mcur, t->d0, trnum, TR_FROM_UP);
01826           traverse_polygon(mcur, t->d1, trnum, TR_FROM_UP);
01827         }
01828     }
01829   
01830   else if ((t->u0 > 0) && (t->u1 > 0)) 
01831     {
01832       if ((t->d0 > 0) && (t->d1 > 0)) /* downward + upward cusps */
01833         {
01834           v0 = tr[t->d1].lseg;
01835           v1 = tr[t->u0].rseg;
01836           retval = SP_2UP_2DN;
01837           if (((dir == TR_FROM_DN) && (t->d1 == from)) ||
01838               ((dir == TR_FROM_UP) && (t->u1 == from)))
01839             {
01840               do_switch = true;
01841               mnew = make_new_monotone_poly(mcur, v1, v0);
01842               traverse_polygon(mcur, t->u1, trnum, TR_FROM_DN);
01843               traverse_polygon(mcur, t->d1, trnum, TR_FROM_UP);
01844               traverse_polygon(mnew, t->u0, trnum, TR_FROM_DN);
01845               traverse_polygon(mnew, t->d0, trnum, TR_FROM_UP);
01846             }
01847           else
01848             {
01849               mnew = make_new_monotone_poly(mcur, v0, v1);
01850               traverse_polygon(mcur, t->u0, trnum, TR_FROM_DN);
01851               traverse_polygon(mcur, t->d0, trnum, TR_FROM_UP);
01852               traverse_polygon(mnew, t->u1, trnum, TR_FROM_DN);
01853               traverse_polygon(mnew, t->d1, trnum, TR_FROM_UP);
01854             }
01855         }
01856       else                      /* only downward cusp */
01857         {
01858           if (_equal_to(&t->lo, &seg[t->lseg].v1))
01859             {
01860               v0 = tr[t->u0].rseg;
01861               v1 = seg[t->lseg].next;
01862 
01863               retval = SP_2UP_LEFT;
01864               if ((dir == TR_FROM_UP) && (t->u0 == from))
01865                 {
01866                   do_switch = true;
01867                   mnew = make_new_monotone_poly(mcur, v1, v0);
01868                   traverse_polygon(mcur, t->u0, trnum, TR_FROM_DN);
01869                   traverse_polygon(mnew, t->d0, trnum, TR_FROM_UP);
01870                   traverse_polygon(mnew, t->u1, trnum, TR_FROM_DN);
01871                   traverse_polygon(mnew, t->d1, trnum, TR_FROM_UP);
01872                 }
01873               else
01874                 {
01875                   mnew = make_new_monotone_poly(mcur, v0, v1);
01876                   traverse_polygon(mcur, t->u1, trnum, TR_FROM_DN);
01877                   traverse_polygon(mcur, t->d0, trnum, TR_FROM_UP);
01878                   traverse_polygon(mcur, t->d1, trnum, TR_FROM_UP);
01879                   traverse_polygon(mnew, t->u0, trnum, TR_FROM_DN);
01880                 }
01881             }
01882           else
01883             {
01884               v0 = t->rseg;
01885               v1 = tr[t->u0].rseg;
01886               retval = SP_2UP_RIGHT;
01887               if ((dir == TR_FROM_UP) && (t->u1 == from))
01888                 {
01889                   do_switch = true;
01890                   mnew = make_new_monotone_poly(mcur, v1, v0);
01891                   traverse_polygon(mcur, t->u1, trnum, TR_FROM_DN);
01892                   traverse_polygon(mnew, t->d1, trnum, TR_FROM_UP);
01893                   traverse_polygon(mnew, t->d0, trnum, TR_FROM_UP);
01894                   traverse_polygon(mnew, t->u0, trnum, TR_FROM_DN);
01895                 }
01896               else
01897                 {
01898                   mnew = make_new_monotone_poly(mcur, v0, v1);
01899                   traverse_polygon(mcur, t->u0, trnum, TR_FROM_DN);
01900                   traverse_polygon(mcur, t->d0, trnum, TR_FROM_UP);
01901                   traverse_polygon(mcur, t->d1, trnum, TR_FROM_UP);
01902                   traverse_polygon(mnew, t->u1, trnum, TR_FROM_DN);
01903                 }
01904             }
01905         }
01906     }
01907   else if ((t->u0 > 0) || (t->u1 > 0)) /* no downward cusp */
01908     {
01909       if ((t->d0 > 0) && (t->d1 > 0)) /* only upward cusp */
01910         {
01911           if (_equal_to(&t->hi, &seg[t->lseg].v0))
01912             {
01913               v0 = tr[t->d1].lseg;
01914               v1 = t->lseg;
01915               retval = SP_2DN_LEFT;
01916               if (!((dir == TR_FROM_DN) && (t->d0 == from)))
01917                 {
01918                   do_switch = true;
01919                   mnew = make_new_monotone_poly(mcur, v1, v0);
01920                   traverse_polygon(mcur, t->u1, trnum, TR_FROM_DN);
01921                   traverse_polygon(mcur, t->d1, trnum, TR_FROM_UP);
01922                   traverse_polygon(mcur, t->u0, trnum, TR_FROM_DN);
01923                   traverse_polygon(mnew, t->d0, trnum, TR_FROM_UP);
01924                 }
01925               else
01926                 {
01927                   mnew = make_new_monotone_poly(mcur, v0, v1);
01928                   traverse_polygon(mcur, t->d0, trnum, TR_FROM_UP);
01929                   traverse_polygon(mnew, t->u0, trnum, TR_FROM_DN);
01930                   traverse_polygon(mnew, t->u1, trnum, TR_FROM_DN);
01931                   traverse_polygon(mnew, t->d1, trnum, TR_FROM_UP);
01932                 }
01933             }
01934           else
01935             {
01936               v0 = tr[t->d1].lseg;
01937               v1 = seg[t->rseg].next;
01938 
01939               retval = SP_2DN_RIGHT;
01940               if ((dir == TR_FROM_DN) && (t->d1 == from))
01941                 {
01942                   do_switch = true;
01943                   mnew = make_new_monotone_poly(mcur, v1, v0);
01944                   traverse_polygon(mcur, t->d1, trnum, TR_FROM_UP);
01945                   traverse_polygon(mnew, t->u1, trnum, TR_FROM_DN);
01946                   traverse_polygon(mnew, t->u0, trnum, TR_FROM_DN);
01947                   traverse_polygon(mnew, t->d0, trnum, TR_FROM_UP);
01948                 }
01949               else
01950                 {
01951                   mnew = make_new_monotone_poly(mcur, v0, v1);
01952                   traverse_polygon(mcur, t->u0, trnum, TR_FROM_DN);
01953                   traverse_polygon(mcur, t->d0, trnum, TR_FROM_UP);
01954                   traverse_polygon(mcur, t->u1, trnum, TR_FROM_DN);
01955                   traverse_polygon(mnew, t->d1, trnum, TR_FROM_UP);
01956                 }
01957             }
01958         }
01959       else                      /* no cusp */
01960         {
01961           if (_equal_to(&t->hi, &seg[t->lseg].v0) &&
01962               _equal_to(&t->lo, &seg[t->rseg].v0))
01963             {
01964               v0 = t->rseg;
01965               v1 = t->lseg;
01966               retval = SP_SIMPLE_LRDN;
01967               if (dir == TR_FROM_UP)
01968                 {
01969                   do_switch = true;
01970                   mnew = make_new_monotone_poly(mcur, v1, v0);
01971                   traverse_polygon(mcur, t->u0, trnum, TR_FROM_DN);
01972                   traverse_polygon(mcur, t->u1, trnum, TR_FROM_DN);
01973                   traverse_polygon(mnew, t->d1, trnum, TR_FROM_UP);
01974                   traverse_polygon(mnew, t->d0, trnum, TR_FROM_UP);
01975                 }
01976               else
01977                 {
01978                   mnew = make_new_monotone_poly(mcur, v0, v1);
01979                   traverse_polygon(mcur, t->d1, trnum, TR_FROM_UP);
01980                   traverse_polygon(mcur, t->d0, trnum, TR_FROM_UP);
01981                   traverse_polygon(mnew, t->u0, trnum, TR_FROM_DN);
01982                   traverse_polygon(mnew, t->u1, trnum, TR_FROM_DN);
01983                 }
01984             }
01985           else if (_equal_to(&t->hi, &seg[t->rseg].v1) &&
01986                    _equal_to(&t->lo, &seg[t->lseg].v1))
01987             {
01988               v0 = seg[t->rseg].next;
01989               v1 = seg[t->lseg].next;
01990 
01991               retval = SP_SIMPLE_LRUP;
01992               if (dir == TR_FROM_UP)
01993                 {
01994                   do_switch = true;
01995                   mnew = make_new_monotone_poly(mcur, v1, v0);
01996                   traverse_polygon(mcur, t->u0, trnum, TR_FROM_DN);
01997                   traverse_polygon(mcur, t->u1, trnum, TR_FROM_DN);
01998                   traverse_polygon(mnew, t->d1, trnum, TR_FROM_UP);
01999                   traverse_polygon(mnew, t->d0, trnum, TR_FROM_UP);
02000                 }
02001               else
02002                 {
02003                   mnew = make_new_monotone_poly(mcur, v0, v1);
02004                   traverse_polygon(mcur, t->d1, trnum, TR_FROM_UP);
02005                   traverse_polygon(mcur, t->d0, trnum, TR_FROM_UP);
02006                   traverse_polygon(mnew, t->u0, trnum, TR_FROM_DN);
02007                   traverse_polygon(mnew, t->u1, trnum, TR_FROM_DN);
02008                 }
02009             }
02010           else                  /* no split possible */
02011             {
02012               retval = SP_NOSPLIT;
02013               traverse_polygon(mcur, t->u0, trnum, TR_FROM_DN);
02014               traverse_polygon(mcur, t->d0, trnum, TR_FROM_UP);
02015               traverse_polygon(mcur, t->u1, trnum, TR_FROM_DN);
02016               traverse_polygon(mcur, t->d1, trnum, TR_FROM_UP);
02017             }
02018         }
02019     }
02020 
02021   return retval;
02022 }
02023 
02024 
02025 /* For each monotone polygon, find the ymax and ymin (to determine the */
02026 /* two y-monotone chains) and pass on this monotone polygon for greedy */
02027 /* triangulation. */
02028 /* Take care not to triangulate duplicate monotone polygons */
02029 
02030 void Triangulator::
02031 triangulate_monotone_polygons(int nvert, int nmonpoly) {
02032   int i;
02033   point_t ymax, ymin;
02034   int p, vfirst, posmax, posmin, v;
02035   int vcount, processed;
02036 
02037   #ifdef DEBUG
02038   for (i = 0; i < nmonpoly; i++)
02039     {
02040       fprintf(stderr, "\n\nPolygon %d: ", i);
02041       vfirst = mchain[mon[i]].vnum;
02042       p = mchain[mon[i]].next;
02043       fprintf (stderr, "%d ", mchain[mon[i]].vnum);
02044       while (mchain[p].vnum != vfirst)
02045         {
02046           fprintf(stderr, "%d ", mchain[p].vnum);
02047           if (mchain[p].vnum == -1) {
02048             fprintf(stderr, " xx");
02049             break;
02050           }
02051           p = mchain[p].next;
02052         }
02053     }
02054   fprintf(stderr, "\n");
02055   #endif
02056 
02057   for (i = 0; i < nmonpoly; i++)
02058     {
02059       vcount = 1;
02060       processed = false;
02061       vfirst = mchain[mon[i]].vnum;
02062       if (vfirst <= 0) {
02063         return;
02064       }
02065       ymin = vert[vfirst].pt;
02066       ymax = ymin;
02067       posmin = mon[i];
02068       posmax = posmin;
02069       mchain[mon[i]].marked = true;
02070       p = mchain[mon[i]].next;
02071       while ((v = mchain[p].vnum) != vfirst)
02072         {
02073           if (v <= 0) {
02074             return;
02075           }
02076           if (mchain[p].marked)
02077            {
02078              processed = true;
02079              break;             /* break from while */
02080            }
02081           else
02082             mchain[p].marked = true;
02083           
02084           if (_greater_than(&vert[v].pt, &ymax))
02085             {
02086               ymax = vert[v].pt;
02087               posmax = p;
02088             }
02089           if (_less_than(&vert[v].pt, &ymin))
02090             {
02091               ymin = vert[v].pt;
02092               posmin = p;
02093             }
02094           p = mchain[p].next;
02095           vcount++;
02096        }
02097 
02098       if (processed)            /* Go to next polygon */
02099         continue;
02100       
02101       if (vcount == 3)          /* already a triangle */
02102         {
02103           _result.push_back(Triangle(this, mchain[p].vnum,
02104                                      mchain[mchain[p].next].vnum,
02105                                      mchain[mchain[p].prev].vnum));
02106         }
02107       else                      /* triangulate the polygon */
02108         {
02109           v = mchain[mchain[posmax].next].vnum;
02110           if (_equal_to(&vert[v].pt, &ymin))
02111             {                   /* LHS is a single line */
02112               triangulate_single_polygon(nvert, posmax, TRI_LHS);
02113             }
02114           else
02115             triangulate_single_polygon(nvert, posmax, TRI_RHS);
02116         }
02117     }
02118 }
02119 
02120 
02121 /* A greedy corner-cutting algorithm to triangulate a y-monotone 
02122  * polygon in O(n) time.
02123  * Joseph O-Rourke, Computational Geometry in C.
02124  */
02125 void Triangulator::
02126 triangulate_single_polygon(int nvert, int posmax, int side) {
02127   int v;
02128   vector_int rc;        /* reflex chain */
02129   int ri;
02130   int endv, tmp, vpos;
02131 
02132   //  cerr << "triangulate_single_polygon(" << nvert << ", " << posmax << ", " << side << ")\n";
02133 
02134   if (side == TRI_RHS)          /* RHS segment is a single segment */
02135     {
02136       rc.push_back(mchain[posmax].vnum);
02137       tmp = mchain[posmax].next;
02138       rc.push_back(mchain[tmp].vnum);
02139       ri = 1;
02140       
02141       vpos = mchain[tmp].next;
02142       v = mchain[vpos].vnum;
02143       
02144       if ((endv = mchain[mchain[posmax].prev].vnum) == 0)
02145         endv = nvert;
02146     }
02147   else                          /* LHS is a single segment */
02148     {
02149       tmp = mchain[posmax].next;
02150       rc.push_back(mchain[tmp].vnum);
02151       tmp = mchain[tmp].next;
02152       rc.push_back(mchain[tmp].vnum);
02153       ri = 1;
02154 
02155       vpos = mchain[tmp].next;
02156       v = mchain[vpos].vnum;
02157 
02158       endv = mchain[posmax].vnum;
02159     }
02160   
02161   while ((v != endv) || (ri > 1))
02162     {
02163       if (v <= 0) {
02164         // Something went wrong.
02165         return;
02166       }
02167       if (ri > 0)               /* reflex chain is non-empty */
02168         {
02169           float crossResult = CROSS(vert[v].pt, vert[rc[ri - 1]].pt,
02170                                     vert[rc[ri]].pt);
02171           if ( crossResult >= 0 )  /* could be convex corner or straight */
02172             {
02173               if ( crossResult > 0)  /* convex corner: cut it off */
02174                 _result.push_back(Triangle(this, rc[ri - 1], rc[ri], v));
02175               /* else : perfectly straight, will be abandoned anyway */
02176               ri--;
02177               rc.pop_back();
02178               nassertv(ri + 1 == (int)rc.size());
02179             }
02180           else /* concave, add v to the chain */
02181             {
02182               ri++;
02183               rc.push_back(v);
02184               nassertv(ri + 1 == (int)rc.size());
02185               vpos = mchain[vpos].next;
02186               v = mchain[vpos].vnum;
02187             }
02188         }
02189       else                      /* reflex-chain empty: add v to the */
02190         {                       /* reflex chain and advance it  */
02191           ri++;
02192           rc.push_back(v);
02193           nassertv(ri + 1 == (int)rc.size());
02194           vpos = mchain[vpos].next;
02195           v = mchain[vpos].vnum;
02196         }
02197     } /* end-while */
02198   
02199   /* reached the bottom vertex. Add in the triangle formed */
02200   _result.push_back(Triangle(this, rc[ri - 1], rc[ri], v));  
02201   ri--;
02202 }
02203 
02204 
 All Classes Functions Variables Enumerations