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