Panda3D
 All Classes Functions Variables Enumerations
triangulator.cxx
1 // Filename: triangulator.cxx
2 // Created by: drose (18Jan07)
3 //
4 ////////////////////////////////////////////////////////////////////
5 //
6 // PANDA 3D SOFTWARE
7 // Copyright (c) Carnegie Mellon University. All rights reserved.
8 //
9 // All use of this software is subject to the terms of the revised BSD
10 // license. You should have received a copy of this license along
11 // with this source code in a file named "LICENSE."
12 //
13 ////////////////////////////////////////////////////////////////////
14 
15 #include "triangulator.h"
16 #include "randomizer.h"
17 
18 ////////////////////////////////////////////////////////////////////
19 // Function: Triangulator::Constructor
20 // Access: Published
21 // Description:
22 ////////////////////////////////////////////////////////////////////
23 Triangulator::
24 Triangulator() {
25 }
26 
27 ////////////////////////////////////////////////////////////////////
28 // Function: Triangulator::clear
29 // Access: Published
30 // Description: Removes all vertices and polygon specifications from
31 // the Triangulator, and prepares it to start over.
32 ////////////////////////////////////////////////////////////////////
33 void Triangulator::
34 clear() {
35  _vertices.clear();
36  clear_polygon();
37 }
38 
39 ////////////////////////////////////////////////////////////////////
40 // Function: Triangulator::add_vertex
41 // Access: Published
42 // Description: Adds a new vertex to the vertex pool. Returns the
43 // vertex index number.
44 ////////////////////////////////////////////////////////////////////
46 add_vertex(const LPoint2d &point) {
47  int index = (int)_vertices.size();
48  _vertices.push_back(point);
49  return index;
50 }
51 
52 ////////////////////////////////////////////////////////////////////
53 // Function: Triangulator::clear_polygon
54 // Access: Published
55 // Description: Removes the current polygon definition (and its set
56 // of holes), but does not clear the vertex pool.
57 ////////////////////////////////////////////////////////////////////
58 void Triangulator::
60  _polygon.clear();
61  _holes.clear();
62 }
63 
64 ////////////////////////////////////////////////////////////////////
65 // Function: Triangulator::add_polygon_vertex
66 // Access: Published
67 // Description: Adds the next consecutive vertex of the polygon.
68 // This vertex should index into the vertex pool
69 // established by repeated calls to add_vertex().
70 //
71 // The vertices may be listed in either clockwise or
72 // counterclockwise order. Vertices should not be
73 // repeated. In particular, do not repeat the first
74 // vertex at the end.
75 ////////////////////////////////////////////////////////////////////
76 void Triangulator::
77 add_polygon_vertex(int index) {
78  _polygon.push_back(index);
79 }
80 
81 ////////////////////////////////////////////////////////////////////
82 // Function: Triangulator::begin_hole
83 // Access: Published
84 // Description: Finishes the previous hole, if any, and prepares to
85 // add a new hole.
86 ////////////////////////////////////////////////////////////////////
87 void Triangulator::
89  _holes.push_back(vector_int());
90 }
91 
92 ////////////////////////////////////////////////////////////////////
93 // Function: Triangulator::add_hole_vertex
94 // Access: Published
95 // Description: Adds the next consecutive vertex of the current hole.
96 // This vertex should index into the vertex pool
97 // established by repeated calls to add_vertex().
98 //
99 // The vertices may be listed in either clockwise or
100 // counterclockwise order. Vertices should not be
101 // repeated.
102 ////////////////////////////////////////////////////////////////////
103 void Triangulator::
104 add_hole_vertex(int index) {
105  nassertv(!_holes.empty());
106  _holes.back().push_back(index);
107 }
108 
109 ////////////////////////////////////////////////////////////////////
110 // Function: Triangulator::triangulate
111 // Access: Published
112 // Description: Does the work of triangulating the specified polygon.
113 // After this call, you may retrieve the new triangles
114 // one at a time by iterating through
115 // get_triangle_v0/1/2().
116 ////////////////////////////////////////////////////////////////////
117 void Triangulator::
119  _result.clear();
120 
121  // Make sure our index numbers are reasonable.
122  cleanup_polygon_indices(_polygon);
123  Holes::iterator hi;
124  for (hi = _holes.begin(); hi != _holes.end(); ++hi) {
125  cleanup_polygon_indices(*hi);
126  }
127 
128  if (_polygon.size() < 3) {
129  // Degenerate case.
130  return;
131  }
132 
133  // Set up the list of segments.
134  seg.clear();
135  seg.push_back(segment_t()); // we don't use the first entry.
136  make_segment(_polygon, true);
137 
138  for (hi = _holes.begin(); hi != _holes.end(); ++hi) {
139  if ((*hi).size() >= 3) {
140  make_segment(*hi, false);
141  }
142  }
143 
144  // Shuffle the segment index.
145  int num_segments = (int)seg.size() - 1;
146  permute.reserve(num_segments);
147  int i;
148  for (i = 0; i < num_segments; ++i) {
149  permute.push_back(i + 1);
150  }
151 
152  // Actually, I'm not sure why we should shuffle the index. That
153  // makes the result non-deterministic, and isn't one order--for
154  // instance, the initial order--as good as any other?
155  /*
156  Randomizer randomizer;
157  for (i = 0; i < num_segments; ++i) {
158  int j = randomizer.random_int(num_segments);
159  nassertv(j >= 0 && j < num_segments);
160  int t = permute[i];
161  permute[i] = permute[j];
162  permute[j] = t;
163  }
164  */
165  choose_idx = 0;
166 
167  /*
168  //cerr << "got " << num_segments << " segments\n";
169  for (i = 1; i < (int)seg.size(); ++i) {
170  segment_t &s = seg[i];
171  printf(" %d. (%g %g), (%g %g)\n", i, s.v0.x, s.v0.y, s.v1.x, s.v1.y);
172  printf(" root0 = %d, root1 = %d\n", s.root0, s.root1);
173  printf(" next = %d, prev = %d\n", s.next, s.prev);
174  }
175  */
176 
177  while (construct_trapezoids(num_segments) != 0) {
178  // If there's an error, re-shuffle the index and try again.
179  Randomizer randomizer;
180  for (i = 0; i < num_segments; ++i) {
181  int j = randomizer.random_int(num_segments);
182  nassertv(j >= 0 && j < num_segments);
183  int t = permute[i];
184  permute[i] = permute[j];
185  permute[j] = t;
186  }
187  choose_idx = 0;
188 
189  /*
190  //cerr << "got " << num_segments << " segments\n";
191  for (i = 1; i < (int)seg.size(); ++i) {
192  segment_t &s = seg[i];
193  printf(" %d. (%g %g), (%g %g)\n", i, s.v0.x, s.v0.y, s.v1.x, s.v1.y);
194  printf(" root0 = %d, root1 = %d\n", s.root0, s.root1);
195  printf(" next = %d, prev = %d\n", s.next, s.prev);
196  }
197  */
198  }
199 
200  /*
201  //cerr << "got " << tr.size() - 1 << " trapezoids\n";
202  for (i = 1; i < (int)tr.size(); ++i) {
203  trap_t &t = tr[i];
204  //cerr << " " << i << ". state = " << t.state << "\n";
205  //cerr << " lseg = " << t.lseg << " rseg = " << t.rseg << "\n";
206  //cerr << " hi = " << t.hi.x << " " << t.hi.y << " lo = " << t.lo.x << " " << t.lo.y << "\n";
207  }
208  */
209 
210  int nmonpoly = monotonate_trapezoids(num_segments);
211 
212  //cerr << "got " << nmonpoly << " monotone polygons\n";
213 
214  triangulate_monotone_polygons(num_segments, nmonpoly);
215 
216  /*
217  Result::iterator ri;
218  for (ri = _result.begin(); ri != _result.end(); ++ri) {
219  //cerr << "tri: " << (*ri)._v0 << " " << (*ri)._v1 << " " << (*ri)._v2 << "\n";
220  }
221  */
222 }
223 
224 ////////////////////////////////////////////////////////////////////
225 // Function: Triangulator::get_num_triangles
226 // Access: Published
227 // Description: Returns the number of triangles generated by the
228 // previous call to triangulate().
229 ////////////////////////////////////////////////////////////////////
230 int Triangulator::
232  return _result.size();
233 }
234 
235 ////////////////////////////////////////////////////////////////////
236 // Function: Triangulator::get_triangle_v0
237 // Access: Published
238 // Description: Returns vertex 0 of the nth triangle generated by the
239 // previous call to triangulate().
240 //
241 // This is a zero-based index into the vertices added by
242 // repeated calls to add_vertex().
243 ////////////////////////////////////////////////////////////////////
244 int Triangulator::
245 get_triangle_v0(int n) const {
246  nassertr(n >= 0 && n < (int)_result.size(), -1);
247  return _result[n]._v0;
248 }
249 
250 ////////////////////////////////////////////////////////////////////
251 // Function: Triangulator::get_triangle_v1
252 // Access: Published
253 // Description: Returns vertex 1 of the nth triangle generated by the
254 // previous call to triangulate().
255 //
256 // This is a zero-based index into the vertices added by
257 // repeated calls to add_vertex().
258 ////////////////////////////////////////////////////////////////////
259 int Triangulator::
260 get_triangle_v1(int n) const {
261  nassertr(n >= 0 && n < (int)_result.size(), -1);
262  return _result[n]._v1;
263 }
264 
265 ////////////////////////////////////////////////////////////////////
266 // Function: Triangulator::get_triangle_v2
267 // Access: Published
268 // Description: Returns vertex 2 of the nth triangle generated by the
269 // previous call to triangulate().
270 //
271 // This is a zero-based index into the vertices added by
272 // repeated calls to add_vertex().
273 ////////////////////////////////////////////////////////////////////
274 int Triangulator::
275 get_triangle_v2(int n) const {
276  nassertr(n >= 0 && n < (int)_result.size(), -1);
277  return _result[n]._v2;
278 }
279 
280 ////////////////////////////////////////////////////////////////////
281 // Function: Triangulator::cleanup_polygon_indices
282 // Access: Protected
283 // Description: Removes any invalid index numbers from the list.
284 ////////////////////////////////////////////////////////////////////
285 void Triangulator::
286 cleanup_polygon_indices(vector_int &polygon) {
287  // First, check for index bounds.
288  size_t pi = 0;
289  while (pi < polygon.size()) {
290  if (polygon[pi] >= 0 && (size_t)polygon[pi] < _vertices.size()) {
291  // This vertex is OK.
292  ++pi;
293  } else {
294  // This index is out-of-bounds; remove it.
295  polygon.erase(_polygon.begin() + pi);
296  }
297  }
298 
299  // Now, remove any consecutive repeated vertices.
300  pi = 1;
301  while (pi < polygon.size()) {
302  if (_vertices[polygon[pi]] != _vertices[polygon[pi - 1]]) {
303  // This vertex is OK.
304  ++pi;
305  } else {
306  // This vertex repeats the previous one; remove it.
307  polygon.erase(_polygon.begin() + pi);
308  }
309  }
310 
311  if (polygon.size() > 1 && _vertices[polygon.back()] == _vertices[_polygon.front()]) {
312  // The last vertex repeats the first one; remove it.
313  polygon.pop_back();
314  }
315 }
316 
317 
318 // The remainder of the code in this file is adapted more or less from
319 // the C code published with the referenced paper.
320 
321 #define T_X 1
322 #define T_Y 2
323 #define T_SINK 3
324 
325 
326 #define FIRSTPT 1 /* checking whether pt. is inserted */
327 #define LASTPT 2
328 
329 
330 #define REALLY_BIG 1<<30
331 #define C_EPS 1.0e-7 /* tolerance value: Used for making */
332  /* all decisions about collinearity or */
333  /* left/right of segment. Decrease */
334  /* this value if the input points are */
335  /* spaced very close together */
336 
337 
338 #define S_LEFT 1 /* for merge-direction */
339 #define S_RIGHT 2
340 
341 
342 #define ST_VALID 1 /* for trapezium state */
343 #define ST_INVALID 2
344 
345 
346 #define SP_SIMPLE_LRUP 1 /* for splitting trapezoids */
347 #define SP_SIMPLE_LRDN 2
348 #define SP_2UP_2DN 3
349 #define SP_2UP_LEFT 4
350 #define SP_2UP_RIGHT 5
351 #define SP_2DN_LEFT 6
352 #define SP_2DN_RIGHT 7
353 #define SP_NOSPLIT -1
354 
355 #define TR_FROM_UP 1 /* for traverse-direction */
356 #define TR_FROM_DN 2
357 
358 #define TRI_LHS 1
359 #define TRI_RHS 2
360 
361 
362 #define CROSS(v0, v1, v2) (((v1).x - (v0).x)*((v2).y - (v0).y) - \
363  ((v1).y - (v0).y)*((v2).x - (v0).x))
364 
365 #define DOT(v0, v1) ((v0).x * (v1).x + (v0).y * (v1).y)
366 
367 #define FP_EQUAL(s, t) (fabs(s - t) <= C_EPS)
368 
369 
370 #define CROSS_SINE(v0, v1) ((v0).x * (v1).y - (v1).x * (v0).y)
371 #define LENGTH(v0) (sqrt((v0).x * (v0).x + (v0).y * (v0).y))
372 
373 
374 ////////////////////////////////////////////////////////////////////
375 // Function: Triangulator::check_left_winding
376 // Access: Private
377 // Description: Returns true if the list of vertices is
378 // counter-clockwise, false if it is clockwise.
379 ////////////////////////////////////////////////////////////////////
380 bool Triangulator::
381 check_left_winding(const vector_int &range) const {
382  // We do this by computing the polygon's signed area. If it comes
383  // out negative, the polygon is right-winding.
384 
385  double area = 0.0;
386  size_t j = range.size() - 1;
387  for (size_t i = 0; i < range.size(); ++i) {
388  const LPoint2d &p0 = _vertices[range[j]];
389  const LPoint2d &p1 = _vertices[range[i]];
390  area += p0[0] * p1[1] - p0[1] * p1[0];
391  j = i;
392  }
393 
394  return area >= 0.0;
395 }
396 
397 ////////////////////////////////////////////////////////////////////
398 // Function: Triangulator::make_segment
399 // Access: Private
400 // Description: Converts a linear list of integer vertices to a list
401 // of segment_t. If want_left_winding is true, the list
402 // is reversed if necessary to make it left-winding;
403 // otherwise, it is reversed to make it right-winding.
404 ////////////////////////////////////////////////////////////////////
405 void Triangulator::
406 make_segment(const vector_int &range, bool want_left_winding) {
407  int num_points = (int)range.size();
408  nassertv(num_points >= 2);
409 
410  int first = (int)seg.size();
411  int last = first + num_points - 1;
412 
413  if (want_left_winding == check_left_winding(range)) {
414  // Keep it in its natural order.
415  int first = (int)seg.size();
416  int last = first + num_points - 1;
417 
418  seg.push_back(segment_t(this, range[0], range[1],
419  last, first + 1));
420 
421  for (int i = 1; i < num_points - 1; ++i) {
422  seg.push_back(segment_t(this, range[i], range[i + 1],
423  first + i - 1, first + i + 1));
424  }
425 
426  seg.push_back(segment_t(this, range[num_points - 1], range[0],
427  last - 1, first));
428 
429  } else {
430  // Reverse it.
431  seg.push_back(segment_t(this, range[0], range[num_points - 1],
432  last, first + 1));
433 
434  for (int i = 1; i < num_points - 1; ++i) {
435  seg.push_back(segment_t(this, range[num_points - i], range[num_points - i - 1],
436  first + i - 1, first + i + 1));
437  }
438 
439  seg.push_back(segment_t(this, range[1], range[0],
440  last - 1, first));
441 
442  }
443 }
444 
445 /* Return the next segment in the generated random ordering of all the */
446 /* segments in S */
447 int Triangulator::
448 choose_segment() {
449  nassertr(choose_idx < (int)permute.size(), 0);
450  // segment_t &s = seg[permute[choose_idx]];
451  // cerr << "choose_segment " << permute[choose_idx] << ": " << s.v0.x << ", " << s.v0.y << " to " << s.v1.x << ", " << s.v1.y << "\n";
452  return permute[choose_idx++];
453 }
454 
455 double Triangulator::
456 math_log2(double v) {
457  static const double log2 = log(2.0);
458  return log(v) / log2;
459 }
460 
461 /* Get log*n for given n */
462 int Triangulator::
463 math_logstar_n(int n) {
464  int i;
465  double v;
466 
467  for (i = 0, v = (double) n; v >= 1; i++)
468  v = math_log2(v);
469 
470  return (i - 1);
471 }
472 
473 
474 int Triangulator::
475 math_N(int n, int h) {
476  int i;
477  double v;
478 
479  for (i = 0, v = (int) n; i < h; i++)
480  v = math_log2(v);
481 
482  return (int) ceil((double) 1.0*n/v);
483 }
484 
485 
486 /* Return a new node to be added into the query tree */
487 int Triangulator::newnode() {
488  int index = (int)qs.size();
489  qs.push_back(node_t());
490  // cerr << "creating new node " << index << "\n";
491  return index;
492 }
493 
494 /* Return a free trapezoid */
495 int Triangulator::newtrap() {
496  int tr_idx = (int)tr.size();
497  tr.push_back(trap_t());
498  tr[tr_idx].lseg = -1;
499  tr[tr_idx].rseg = -1;
500  tr[tr_idx].state = ST_VALID;
501  // cerr << "creating new trapezoid " << tr_idx << "\n";
502  return tr_idx;
503 }
504 
505 
506 /* Return the maximum of the two points into the yval structure */
507 int Triangulator::_max(point_t *yval, point_t *v0, point_t *v1) {
508  if (v0->y > v1->y + C_EPS)
509  *yval = *v0;
510  else if (FP_EQUAL(v0->y, v1->y))
511  {
512  if (v0->x > v1->x + C_EPS)
513  *yval = *v0;
514  else
515  *yval = *v1;
516  }
517  else
518  *yval = *v1;
519 
520  return 0;
521 }
522 
523 
524 /* Return the minimum of the two points into the yval structure */
525 int Triangulator::_min(point_t *yval, point_t *v0, point_t *v1) {
526  if (v0->y < v1->y - C_EPS)
527  *yval = *v0;
528  else if (FP_EQUAL(v0->y, v1->y))
529  {
530  if (v0->x < v1->x)
531  *yval = *v0;
532  else
533  *yval = *v1;
534  }
535  else
536  *yval = *v1;
537 
538  return 0;
539 }
540 
541 
542 int Triangulator::
543 _greater_than(point_t *v0, point_t *v1) {
544  if (v0->y > v1->y + C_EPS)
545  return true;
546  else if (v0->y < v1->y - C_EPS)
547  return false;
548  else
549  return (v0->x > v1->x);
550 }
551 
552 
553 int Triangulator::
554 _equal_to(point_t *v0, point_t *v1) {
555  return (FP_EQUAL(v0->y, v1->y) && FP_EQUAL(v0->x, v1->x));
556 }
557 
558 int Triangulator::
559 _greater_than_equal_to(point_t *v0, point_t *v1) {
560  if (v0->y > v1->y + C_EPS)
561  return true;
562  else if (v0->y < v1->y - C_EPS)
563  return false;
564  else
565  return (v0->x >= v1->x);
566 }
567 
568 int Triangulator::
569 _less_than(point_t *v0, point_t *v1) {
570  if (v0->y < v1->y - C_EPS)
571  return true;
572  else if (v0->y > v1->y + C_EPS)
573  return false;
574  else
575  return (v0->x < v1->x);
576 }
577 
578 
579 /* Initilialise the query structure (Q) and the trapezoid table (T)
580  * when the first segment is added to start the trapezoidation. The
581  * query-tree starts out with 4 trapezoids, one S-node and 2 Y-nodes
582  *
583  * 4
584  * -----------------------------------
585  * \
586  * 1 \ 2
587  * \
588  * -----------------------------------
589  * 3
590  */
591 
592 int Triangulator::
593 init_query_structure(int segnum) {
594  int i1, i2, i3, i4, i5, i6, i7, root;
595  int t1, t2, t3, t4;
596  segment_t *s = &seg[segnum];
597 
598  tr.clear();
599  qs.clear();
600 
601  // We don't use the first elements.
602  tr.push_back(trap_t());
603  qs.push_back(node_t());
604 
605  i1 = newnode();
606  qs[i1].nodetype = T_Y;
607  _max(&qs[i1].yval, &s->v0, &s->v1); /* root */
608  root = i1;
609 
610  i2 = newnode();
611  qs[i1].right = i2;
612  qs[i2].nodetype = T_SINK;
613  qs[i2].parent = i1;
614 
615  i3 = newnode();
616  qs[i1].left = i3;
617 
618  qs[i3].nodetype = T_Y;
619  _min(&qs[i3].yval, &s->v0, &s->v1); /* root */
620  qs[i3].parent = i1;
621 
622  i4 = newnode();
623  qs[i3].left = i4;
624  qs[i4].nodetype = T_SINK;
625  qs[i4].parent = i3;
626 
627  i5 = newnode();
628  qs[i3].right = i5;
629  qs[i5].nodetype = T_X;
630  qs[i5].segnum = segnum;
631  qs[i5].parent = i3;
632 
633  i6 = newnode();
634  qs[i5].left = i6;
635  qs[i6].nodetype = T_SINK;
636  qs[i6].parent = i5;
637 
638  i7 = newnode();
639  qs[i5].right = i7;
640  qs[i7].nodetype = T_SINK;
641  qs[i7].parent = i5;
642 
643  t1 = newtrap(); /* middle left */
644  t2 = newtrap(); /* middle right */
645  t3 = newtrap(); /* bottom-most */
646  t4 = newtrap(); /* topmost */
647 
648  tr[t4].lo = qs[i1].yval;
649  tr[t2].hi = qs[i1].yval;
650  tr[t1].hi = qs[i1].yval;
651  tr[t3].hi = qs[i3].yval;
652  tr[t2].lo = qs[i3].yval;
653  tr[t1].lo = qs[i3].yval;
654  tr[t4].hi.y = (double) (REALLY_BIG);
655  tr[t4].hi.x = (double) (REALLY_BIG);
656  tr[t3].lo.y = (double) -1* (REALLY_BIG);
657  tr[t3].lo.x = (double) -1* (REALLY_BIG);
658  tr[t2].lseg = segnum;
659  tr[t1].rseg = segnum;
660  tr[t2].u0 = t4;
661  tr[t1].u0 = t4;
662  tr[t2].d0 = t3;
663  tr[t1].d0 = t3;
664  tr[t3].u0 = t1;
665  tr[t4].d0 = t1;
666  tr[t3].u1 = t2;
667  tr[t4].d1 = t2;
668 
669  tr[t1].sink = i6;
670  tr[t2].sink = i7;
671  tr[t3].sink = i4;
672  tr[t4].sink = i2;
673 
674  tr[t2].state = ST_VALID;
675  tr[t1].state = ST_VALID;
676  tr[t4].state = ST_VALID;
677  tr[t3].state = ST_VALID;
678 
679  qs[i2].trnum = t4;
680  qs[i4].trnum = t3;
681  qs[i6].trnum = t1;
682  qs[i7].trnum = t2;
683 
684  s->is_inserted = true;
685  return root;
686 }
687 
688 
689 /* Retun true if the vertex v is to the left of line segment no.
690  * segnum. Takes care of the degenerate cases when both the vertices
691  * have the same y--cood, etc.
692  */
693 
694 int Triangulator::
695 is_left_of(int segnum, point_t *v) {
696  segment_t *s = &seg[segnum];
697  double area;
698 
699  if (_greater_than(&s->v1, &s->v0)) /* seg. going upwards */
700  {
701  if (FP_EQUAL(s->v1.y, v->y))
702  {
703  if (v->x < s->v1.x)
704  area = 1.0;
705  else
706  area = -1.0;
707  }
708  else if (FP_EQUAL(s->v0.y, v->y))
709  {
710  if (v->x < s->v0.x)
711  area = 1.0;
712  else
713  area = -1.0;
714  }
715  else
716  area = CROSS(s->v0, s->v1, (*v));
717  }
718  else /* v0 > v1 */
719  {
720  if (FP_EQUAL(s->v1.y, v->y))
721  {
722  if (v->x < s->v1.x)
723  area = 1.0;
724  else
725  area = -1.0;
726  }
727  else if (FP_EQUAL(s->v0.y, v->y))
728  {
729  if (v->x < s->v0.x)
730  area = 1.0;
731  else
732  area = -1.0;
733  }
734  else
735  area = CROSS(s->v1, s->v0, (*v));
736  }
737 
738  if (area > 0.0)
739  return true;
740  else
741  return false;
742 }
743 
744 
745 
746 /* Returns true if the corresponding endpoint of the given segment is */
747 /* already inserted into the segment tree. Use the simple test of */
748 /* whether the segment which shares this endpoint is already inserted */
749 
750 int Triangulator::
751 inserted(int segnum, int whichpt) {
752  if (whichpt == FIRSTPT)
753  return seg[seg[segnum].prev].is_inserted;
754  else
755  return seg[seg[segnum].next].is_inserted;
756 }
757 
758 /* This is query routine which determines which trapezoid does the
759  * point v lie in. The return value is the trapezoid number.
760  */
761 
762 int Triangulator::
763 locate_endpoint(point_t *v, point_t *vo, int r) {
764  // cerr << "locate_endpoint(" << v->x << " " << v->y << ", " << vo->x << " " << vo->y << ", " << r << ")\n";
765  node_t *rptr = &qs[r];
766 
767  switch (rptr->nodetype)
768  {
769  case T_SINK:
770  return rptr->trnum;
771 
772  case T_Y:
773  if (_greater_than(v, &rptr->yval)) /* above */
774  return locate_endpoint(v, vo, rptr->right);
775  else if (_equal_to(v, &rptr->yval)) /* the point is already */
776  { /* inserted. */
777  if (_greater_than(vo, &rptr->yval)) /* above */
778  return locate_endpoint(v, vo, rptr->right);
779  else
780  return locate_endpoint(v, vo, rptr->left); /* below */
781  }
782  else
783  return locate_endpoint(v, vo, rptr->left); /* below */
784 
785  case T_X:
786  if (_equal_to(v, &seg[rptr->segnum].v0) ||
787  _equal_to(v, &seg[rptr->segnum].v1))
788  {
789  if (FP_EQUAL(v->y, vo->y)) /* horizontal segment */
790  {
791  if (vo->x < v->x)
792  return locate_endpoint(v, vo, rptr->left); /* left */
793  else
794  return locate_endpoint(v, vo, rptr->right); /* right */
795  }
796 
797  else if (is_left_of(rptr->segnum, vo))
798  return locate_endpoint(v, vo, rptr->left); /* left */
799  else
800  return locate_endpoint(v, vo, rptr->right); /* right */
801  }
802  else if (is_left_of(rptr->segnum, v))
803  return locate_endpoint(v, vo, rptr->left); /* left */
804  else
805  return locate_endpoint(v, vo, rptr->right); /* right */
806 
807  default:
808  fprintf(stderr, "Haggu !!!!!\n");
809  nassertr(false, -1);
810  return -1;
811  }
812 }
813 
814 
815 /* Thread in the segment into the existing trapezoidation. The
816  * limiting trapezoids are given by tfirst and tlast (which are the
817  * trapezoids containing the two endpoints of the segment. Merges all
818  * possible trapezoids which flank this segment and have been recently
819  * divided because of its insertion
820  */
821 
822 int Triangulator::
823 merge_trapezoids(int segnum, int tfirst, int tlast, int side) {
824  int t, tnext, cond;
825  int ptnext;
826 
827  // cerr << "merge_trapezoids(" << segnum << ", " << tfirst << ", " << tlast << ", " << side << ")\n";
828 
829  /* First merge polys on the LHS */
830  t = tfirst;
831  while ((t > 0) && _greater_than_equal_to(&tr[t].lo, &tr[tlast].lo))
832  {
833  if (side == S_LEFT)
834  cond = ((((tnext = tr[t].d0) > 0) && (tr[tnext].rseg == segnum)) ||
835  (((tnext = tr[t].d1) > 0) && (tr[tnext].rseg == segnum)));
836  else
837  cond = ((((tnext = tr[t].d0) > 0) && (tr[tnext].lseg == segnum)) ||
838  (((tnext = tr[t].d1) > 0) && (tr[tnext].lseg == segnum)));
839 
840  if (cond)
841  {
842  if ((tr[t].lseg == tr[tnext].lseg) &&
843  (tr[t].rseg == tr[tnext].rseg)) /* good neighbours */
844  { /* merge them */
845  /* Use the upper node as the new node i.e. t */
846 
847  ptnext = qs[tr[tnext].sink].parent;
848 
849  if (qs[ptnext].left == tr[tnext].sink)
850  qs[ptnext].left = tr[t].sink;
851  else
852  qs[ptnext].right = tr[t].sink; /* redirect parent */
853 
854 
855  /* Change the upper neighbours of the lower trapezoids */
856 
857  if ((tr[t].d0 = tr[tnext].d0) > 0) {
858  if (tr[tr[t].d0].u0 == tnext) {
859  tr[tr[t].d0].u0 = t;
860  } else if (tr[tr[t].d0].u1 == tnext) {
861  tr[tr[t].d0].u1 = t;
862  }
863  }
864 
865  if ((tr[t].d1 = tr[tnext].d1) > 0) {
866  if (tr[tr[t].d1].u0 == tnext) {
867  tr[tr[t].d1].u0 = t;
868  } else if (tr[tr[t].d1].u1 == tnext) {
869  tr[tr[t].d1].u1 = t;
870  }
871  }
872 
873  tr[t].lo = tr[tnext].lo;
874  tr[tnext].state = ST_INVALID; /* invalidate the lower */
875  /* trapezium */
876  }
877  else /* not good neighbours */
878  t = tnext;
879  }
880  else /* do not satisfy the outer if */
881  t = tnext;
882 
883  } /* end-while */
884 
885  return 0;
886 }
887 
888 
889 /* Add in the new segment into the trapezoidation and update Q and T
890  * structures. First locate the two endpoints of the segment in the
891  * Q-structure. Then start from the topmost trapezoid and go down to
892  * the lower trapezoid dividing all the trapezoids in between .
893  */
894 
895 int Triangulator::
896 add_segment(int segnum) {
897  //cerr << "add_segment(" << segnum << ")\n";
898 
899  segment_t s;
900  // segment_t *so = &seg[segnum];
901  int tu, tl, sk, tfirst, tlast; //, tnext;
902  int tfirstr = 0, tlastr = 0, tfirstl = 0, tlastl = 0;
903  int i1, i2, t, tn; // t1, t2,
904  point_t tpt;
905  int tritop = 0, tribot = 0, is_swapped = 0;
906  int tmptriseg;
907 
908  s = seg[segnum];
909  if (_greater_than(&s.v1, &s.v0)) /* Get higher vertex in v0 */
910  {
911  int tmp;
912  tpt = s.v0;
913  s.v0 = s.v1;
914  s.v1 = tpt;
915  tmp = s.root0;
916  s.root0 = s.root1;
917  s.root1 = tmp;
918  is_swapped = true;
919  }
920 
921  if ((is_swapped) ? !inserted(segnum, LASTPT) :
922  !inserted(segnum, FIRSTPT)) /* insert v0 in the tree */
923  {
924  int tmp_d;
925 
926  tu = locate_endpoint(&s.v0, &s.v1, s.root0);
927  tl = newtrap(); /* tl is the new lower trapezoid */
928  tr[tl].state = ST_VALID;
929  tr[tl] = tr[tu];
930  tr[tl].hi.y = s.v0.y;
931  tr[tu].lo.y = s.v0.y;
932  tr[tl].hi.x = s.v0.x;
933  tr[tu].lo.x = s.v0.x;
934  tr[tu].d0 = tl;
935  tr[tu].d1 = 0;
936  tr[tl].u0 = tu;
937  tr[tl].u1 = 0;
938 
939  if (((tmp_d = tr[tl].d0) > 0) && (tr[tmp_d].u0 == tu))
940  tr[tmp_d].u0 = tl;
941  if (((tmp_d = tr[tl].d0) > 0) && (tr[tmp_d].u1 == tu))
942  tr[tmp_d].u1 = tl;
943 
944  if (((tmp_d = tr[tl].d1) > 0) && (tr[tmp_d].u0 == tu))
945  tr[tmp_d].u0 = tl;
946  if (((tmp_d = tr[tl].d1) > 0) && (tr[tmp_d].u1 == tu))
947  tr[tmp_d].u1 = tl;
948 
949  /* Now update the query structure and obtain the sinks for the */
950  /* two trapezoids */
951 
952  i1 = newnode(); /* Upper trapezoid sink */
953  i2 = newnode(); /* Lower trapezoid sink */
954  sk = tr[tu].sink;
955 
956  qs[sk].nodetype = T_Y;
957  qs[sk].yval = s.v0;
958  qs[sk].segnum = segnum; /* not really reqd ... maybe later */
959  qs[sk].left = i2;
960  qs[sk].right = i1;
961 
962  qs[i1].nodetype = T_SINK;
963  qs[i1].trnum = tu;
964  qs[i1].parent = sk;
965 
966  qs[i2].nodetype = T_SINK;
967  qs[i2].trnum = tl;
968  qs[i2].parent = sk;
969 
970  tr[tu].sink = i1;
971  tr[tl].sink = i2;
972  tfirst = tl;
973  }
974  else /* v0 already present */
975  { /* Get the topmost intersecting trapezoid */
976  tfirst = locate_endpoint(&s.v0, &s.v1, s.root0);
977  tritop = 1;
978  }
979 
980 
981  if ((is_swapped) ? !inserted(segnum, FIRSTPT) :
982  !inserted(segnum, LASTPT)) /* insert v1 in the tree */
983  {
984  int tmp_d;
985 
986  tu = locate_endpoint(&s.v1, &s.v0, s.root1);
987 
988  tl = newtrap(); /* tl is the new lower trapezoid */
989  tr[tl].state = ST_VALID;
990  tr[tl] = tr[tu];
991  tr[tl].hi.y = s.v1.y;
992  tr[tu].lo.y = s.v1.y;
993  tr[tl].hi.x = s.v1.x;
994  tr[tu].lo.x = s.v1.x;
995  tr[tu].d0 = tl;
996  tr[tu].d1 = 0;
997  tr[tl].u0 = tu;
998  tr[tl].u1 = 0;
999 
1000  if (((tmp_d = tr[tl].d0) > 0) && (tr[tmp_d].u0 == tu))
1001  tr[tmp_d].u0 = tl;
1002  if (((tmp_d = tr[tl].d0) > 0) && (tr[tmp_d].u1 == tu))
1003  tr[tmp_d].u1 = tl;
1004 
1005  if (((tmp_d = tr[tl].d1) > 0) && (tr[tmp_d].u0 == tu))
1006  tr[tmp_d].u0 = tl;
1007  if (((tmp_d = tr[tl].d1) > 0) && (tr[tmp_d].u1 == tu))
1008  tr[tmp_d].u1 = tl;
1009 
1010  /* Now update the query structure and obtain the sinks for the */
1011  /* two trapezoids */
1012 
1013  i1 = newnode(); /* Upper trapezoid sink */
1014  i2 = newnode(); /* Lower trapezoid sink */
1015  sk = tr[tu].sink;
1016 
1017  qs[sk].nodetype = T_Y;
1018  qs[sk].yval = s.v1;
1019  qs[sk].segnum = segnum; /* not really reqd ... maybe later */
1020  qs[sk].left = i2;
1021  qs[sk].right = i1;
1022 
1023  qs[i1].nodetype = T_SINK;
1024  qs[i1].trnum = tu;
1025  qs[i1].parent = sk;
1026 
1027  qs[i2].nodetype = T_SINK;
1028  qs[i2].trnum = tl;
1029  qs[i2].parent = sk;
1030 
1031  tr[tu].sink = i1;
1032  tr[tl].sink = i2;
1033  tlast = tu;
1034  }
1035  else /* v1 already present */
1036  { /* Get the lowermost intersecting trapezoid */
1037  tlast = locate_endpoint(&s.v1, &s.v0, s.root1);
1038  tribot = 1;
1039  }
1040 
1041  /* Thread the segment into the query tree creating a new X-node */
1042  /* First, split all the trapezoids which are intersected by s into */
1043  /* two */
1044 
1045  t = tfirst; /* topmost trapezoid */
1046 
1047  while ((t > 0) &&
1048  _greater_than_equal_to(&tr[t].lo, &tr[tlast].lo))
1049  /* traverse from top to bot */
1050  {
1051  int t_sav, tn_sav;
1052  sk = tr[t].sink;
1053  i1 = newnode(); /* left trapezoid sink */
1054  i2 = newnode(); /* right trapezoid sink */
1055 
1056  qs[sk].nodetype = T_X;
1057  qs[sk].segnum = segnum;
1058  qs[sk].left = i1;
1059  qs[sk].right = i2;
1060 
1061  qs[i1].nodetype = T_SINK; /* left trapezoid (use existing one) */
1062  qs[i1].trnum = t;
1063  qs[i1].parent = sk;
1064 
1065  qs[i2].nodetype = T_SINK; /* right trapezoid (allocate new) */
1066  tn = newtrap();
1067  qs[i2].trnum = tn;
1068  tr[tn].state = ST_VALID;
1069  qs[i2].parent = sk;
1070 
1071  if (t == tfirst)
1072  tfirstr = tn;
1073  if (_equal_to(&tr[t].lo, &tr[tlast].lo))
1074  tlastr = tn;
1075 
1076  tr[tn] = tr[t];
1077  tr[t].sink = i1;
1078  tr[tn].sink = i2;
1079  t_sav = t;
1080  tn_sav = tn;
1081 
1082  /* error */
1083 
1084  if ((tr[t].d0 <= 0) && (tr[t].d1 <= 0)) /* case cannot arise */
1085  {
1086  /* Actually, this case does sometimes arise. Huh. */
1087  fprintf(stderr, "add_segment: error\n");
1088  return 1;
1089  }
1090 
1091  /* only one trapezoid below. partition t into two and make the */
1092  /* two resulting trapezoids t and tn as the upper neighbours of */
1093  /* the sole lower trapezoid */
1094 
1095  else if ((tr[t].d0 > 0) && (tr[t].d1 <= 0))
1096  { /* Only one trapezoid below */
1097  if ((tr[t].u0 > 0) && (tr[t].u1 > 0))
1098  { /* continuation of a chain from abv. */
1099  if (tr[t].usave > 0) /* three upper neighbours */
1100  {
1101  if (tr[t].uside == S_LEFT)
1102  {
1103  tr[tn].u0 = tr[t].u1;
1104  tr[t].u1 = -1;
1105  tr[tn].u1 = tr[t].usave;
1106 
1107  tr[tr[t].u0].d0 = t;
1108  tr[tr[tn].u0].d0 = tn;
1109  tr[tr[tn].u1].d0 = tn;
1110  }
1111  else /* intersects in the right */
1112  {
1113  tr[tn].u1 = -1;
1114  tr[tn].u0 = tr[t].u1;
1115  tr[t].u1 = tr[t].u0;
1116  tr[t].u0 = tr[t].usave;
1117 
1118  tr[tr[t].u0].d0 = t;
1119  tr[tr[t].u1].d0 = t;
1120  tr[tr[tn].u0].d0 = tn;
1121  }
1122  tr[tn].usave = 0;
1123  tr[t].usave = 0;
1124  }
1125  else /* No usave.... simple case */
1126  {
1127  tr[tn].u0 = tr[t].u1;
1128  tr[tn].u1 = -1;
1129  tr[t].u1 = -1;
1130  tr[tr[tn].u0].d0 = tn;
1131  }
1132  }
1133  else
1134  { /* fresh seg. or upward cusp */
1135  int tmp_u = tr[t].u0;
1136  int td0, td1;
1137  if (((td0 = tr[tmp_u].d0) > 0) &&
1138  ((td1 = tr[tmp_u].d1) > 0))
1139  { /* upward cusp */
1140  if ((tr[td0].rseg > 0) &&
1141  !is_left_of(tr[td0].rseg, &s.v1))
1142  {
1143  tr[tn].u1 = -1;
1144  tr[t].u1 = -1;
1145  tr[t].u0 = -1;
1146  tr[tr[tn].u0].d1 = tn;
1147  }
1148  else /* cusp going leftwards */
1149  {
1150  tr[t].u1 = -1;
1151  tr[tn].u1 = -1;
1152  tr[tn].u0 = -1;
1153  tr[tr[t].u0].d0 = t;
1154  }
1155  }
1156  else /* fresh segment */
1157  {
1158  tr[tr[t].u0].d0 = t;
1159  tr[tr[t].u0].d1 = tn;
1160  }
1161  }
1162 
1163  if (FP_EQUAL(tr[t].lo.y, tr[tlast].lo.y) &&
1164  FP_EQUAL(tr[t].lo.x, tr[tlast].lo.x) && tribot)
1165  { /* bottom forms a triangle */
1166 
1167  if (is_swapped)
1168  tmptriseg = seg[segnum].prev;
1169  else
1170  tmptriseg = seg[segnum].next;
1171 
1172  if ((tmptriseg > 0) && is_left_of(tmptriseg, &s.v0))
1173  {
1174  /* L-R downward cusp */
1175  tr[tr[t].d0].u0 = t;
1176  tr[tn].d1 = -1;
1177  tr[tn].d0 = -1;
1178  }
1179  else
1180  {
1181  /* R-L downward cusp */
1182  tr[tr[tn].d0].u1 = tn;
1183  tr[t].d1 = -1;
1184  tr[t].d0 = -1;
1185  }
1186  }
1187  else
1188  {
1189  if ((tr[tr[t].d0].u0 > 0) && (tr[tr[t].d0].u1 > 0))
1190  {
1191  if (tr[tr[t].d0].u0 == t) /* passes thru LHS */
1192  {
1193  tr[tr[t].d0].usave = tr[tr[t].d0].u1;
1194  tr[tr[t].d0].uside = S_LEFT;
1195  }
1196  else
1197  {
1198  tr[tr[t].d0].usave = tr[tr[t].d0].u0;
1199  tr[tr[t].d0].uside = S_RIGHT;
1200  }
1201  }
1202  tr[tr[t].d0].u0 = t;
1203  tr[tr[t].d0].u1 = tn;
1204  }
1205 
1206  t = tr[t].d0;
1207  }
1208 
1209 
1210  else if ((tr[t].d0 <= 0) && (tr[t].d1 > 0))
1211  { /* Only one trapezoid below */
1212  if ((tr[t].u0 > 0) && (tr[t].u1 > 0))
1213  { /* continuation of a chain from abv. */
1214  if (tr[t].usave > 0) /* three upper neighbours */
1215  {
1216  if (tr[t].uside == S_LEFT)
1217  {
1218  tr[tn].u0 = tr[t].u1;
1219  tr[t].u1 = -1;
1220  tr[tn].u1 = tr[t].usave;
1221 
1222  tr[tr[t].u0].d0 = t;
1223  tr[tr[tn].u0].d0 = tn;
1224  tr[tr[tn].u1].d0 = tn;
1225  }
1226  else /* intersects in the right */
1227  {
1228  tr[tn].u1 = -1;
1229  tr[tn].u0 = tr[t].u1;
1230  tr[t].u1 = tr[t].u0;
1231  tr[t].u0 = tr[t].usave;
1232 
1233  tr[tr[t].u0].d0 = t;
1234  tr[tr[t].u1].d0 = t;
1235  tr[tr[tn].u0].d0 = tn;
1236  }
1237 
1238  tr[tn].usave = 0;
1239  tr[t].usave = 0;
1240  }
1241  else /* No usave.... simple case */
1242  {
1243  tr[tn].u0 = tr[t].u1;
1244  tr[tn].u1 = -1;
1245  tr[t].u1 = -1;
1246  tr[tr[tn].u0].d0 = tn;
1247  }
1248  }
1249  else
1250  { /* fresh seg. or upward cusp */
1251  int tmp_u = tr[t].u0;
1252  int td0, td1;
1253  if (((td0 = tr[tmp_u].d0) > 0) &&
1254  ((td1 = tr[tmp_u].d1) > 0))
1255  { /* upward cusp */
1256  if ((tr[td0].rseg > 0) &&
1257  !is_left_of(tr[td0].rseg, &s.v1))
1258  {
1259  tr[tn].u1 = -1;
1260  tr[t].u1 = -1;
1261  tr[t].u0 = -1;
1262  tr[tr[tn].u0].d1 = tn;
1263  }
1264  else
1265  {
1266  tr[t].u1 = -1;
1267  tr[tn].u1 = -1;
1268  tr[tn].u0 = -1;
1269  tr[tr[t].u0].d0 = t;
1270  }
1271  }
1272  else /* fresh segment */
1273  {
1274  tr[tr[t].u0].d0 = t;
1275  tr[tr[t].u0].d1 = tn;
1276  }
1277  }
1278 
1279  if (FP_EQUAL(tr[t].lo.y, tr[tlast].lo.y) &&
1280  FP_EQUAL(tr[t].lo.x, tr[tlast].lo.x) && tribot)
1281  { /* bottom forms a triangle */
1282  if (is_swapped)
1283  tmptriseg = seg[segnum].prev;
1284  else
1285  tmptriseg = seg[segnum].next;
1286 
1287  if ((tmptriseg > 0) && is_left_of(tmptriseg, &s.v0))
1288  {
1289  /* L-R downward cusp */
1290  tr[tr[t].d1].u0 = t;
1291  tr[tn].d1 = -1;
1292  tr[tn].d0 = -1;
1293  }
1294  else
1295  {
1296  /* R-L downward cusp */
1297  tr[tr[tn].d1].u1 = tn;
1298  tr[t].d1 = -1;
1299  tr[t].d0 = -1;
1300  }
1301  }
1302  else
1303  {
1304  if ((tr[tr[t].d1].u0 > 0) && (tr[tr[t].d1].u1 > 0))
1305  {
1306  if (tr[tr[t].d1].u0 == t) /* passes thru LHS */
1307  {
1308  tr[tr[t].d1].usave = tr[tr[t].d1].u1;
1309  tr[tr[t].d1].uside = S_LEFT;
1310  }
1311  else
1312  {
1313  tr[tr[t].d1].usave = tr[tr[t].d1].u0;
1314  tr[tr[t].d1].uside = S_RIGHT;
1315  }
1316  }
1317  tr[tr[t].d1].u0 = t;
1318  tr[tr[t].d1].u1 = tn;
1319  }
1320 
1321  t = tr[t].d1;
1322  }
1323 
1324  /* two trapezoids below. Find out which one is intersected by */
1325  /* this segment and proceed down that one */
1326 
1327  else
1328  {
1329  // int tmpseg = tr[tr[t].d0].rseg;
1330  double y0, yt;
1331  point_t tmppt;
1332  int tnext, i_d0, i_d1;
1333 
1334  i_d1 = false;
1335  i_d0 = false;
1336  if (FP_EQUAL(tr[t].lo.y, s.v0.y))
1337  {
1338  if (tr[t].lo.x > s.v0.x)
1339  i_d0 = true;
1340  else
1341  i_d1 = true;
1342  }
1343  else
1344  {
1345  y0 = tr[t].lo.y;
1346  tmppt.y = y0;
1347  yt = (y0 - s.v0.y)/(s.v1.y - s.v0.y);
1348  tmppt.x = s.v0.x + yt * (s.v1.x - s.v0.x);
1349 
1350  if (_less_than(&tmppt, &tr[t].lo))
1351  i_d0 = true;
1352  else
1353  i_d1 = true;
1354  }
1355 
1356  /* check continuity from the top so that the lower-neighbour */
1357  /* values are properly filled for the upper trapezoid */
1358 
1359  if ((tr[t].u0 > 0) && (tr[t].u1 > 0))
1360  { /* continuation of a chain from abv. */
1361  if (tr[t].usave > 0) /* three upper neighbours */
1362  {
1363  if (tr[t].uside == S_LEFT)
1364  {
1365  tr[tn].u0 = tr[t].u1;
1366  tr[t].u1 = -1;
1367  tr[tn].u1 = tr[t].usave;
1368 
1369  tr[tr[t].u0].d0 = t;
1370  tr[tr[tn].u0].d0 = tn;
1371  tr[tr[tn].u1].d0 = tn;
1372  }
1373  else /* intersects in the right */
1374  {
1375  tr[tn].u1 = -1;
1376  tr[tn].u0 = tr[t].u1;
1377  tr[t].u1 = tr[t].u0;
1378  tr[t].u0 = tr[t].usave;
1379 
1380  tr[tr[t].u0].d0 = t;
1381  tr[tr[t].u1].d0 = t;
1382  tr[tr[tn].u0].d0 = tn;
1383  }
1384 
1385  tr[tn].usave = 0;
1386  tr[t].usave = 0;
1387  }
1388  else /* No usave.... simple case */
1389  {
1390  tr[tn].u0 = tr[t].u1;
1391  tr[tn].u1 = -1;
1392  tr[t].u1 = -1;
1393  tr[tr[tn].u0].d0 = tn;
1394  }
1395  }
1396  else
1397  { /* fresh seg. or upward cusp */
1398  int tmp_u = tr[t].u0;
1399  int td0, td1;
1400  if (((td0 = tr[tmp_u].d0) > 0) &&
1401  ((td1 = tr[tmp_u].d1) > 0))
1402  { /* upward cusp */
1403  if ((tr[td0].rseg > 0) &&
1404  !is_left_of(tr[td0].rseg, &s.v1))
1405  {
1406  tr[tn].u1 = -1;
1407  tr[t].u1 = -1;
1408  tr[t].u0 = -1;
1409  tr[tr[tn].u0].d1 = tn;
1410  }
1411  else
1412  {
1413  tr[t].u1 = -1;
1414  tr[tn].u1 = -1;
1415  tr[tn].u0 = -1;
1416  tr[tr[t].u0].d0 = t;
1417  }
1418  }
1419  else /* fresh segment */
1420  {
1421  tr[tr[t].u0].d0 = t;
1422  tr[tr[t].u0].d1 = tn;
1423  }
1424  }
1425 
1426  if (FP_EQUAL(tr[t].lo.y, tr[tlast].lo.y) &&
1427  FP_EQUAL(tr[t].lo.x, tr[tlast].lo.x) && tribot)
1428  {
1429  /* this case arises only at the lowest trapezoid.. i.e.
1430  tlast, if the lower endpoint of the segment is
1431  already inserted in the structure */
1432 
1433  tr[tr[t].d0].u0 = t;
1434  tr[tr[t].d0].u1 = -1;
1435  tr[tr[t].d1].u0 = tn;
1436  tr[tr[t].d1].u1 = -1;
1437 
1438  tr[tn].d0 = tr[t].d1;
1439  tr[tn].d1 = -1;
1440  tr[t].d1 = -1;
1441 
1442  tnext = tr[t].d1;
1443  }
1444  else if (i_d0)
1445  /* intersecting d0 */
1446  {
1447  tr[tr[t].d0].u0 = t;
1448  tr[tr[t].d0].u1 = tn;
1449  tr[tr[t].d1].u0 = tn;
1450  tr[tr[t].d1].u1 = -1;
1451 
1452  /* new code to determine the bottom neighbours of the */
1453  /* newly partitioned trapezoid */
1454 
1455  tr[t].d1 = -1;
1456 
1457  tnext = tr[t].d0;
1458  }
1459  else /* intersecting d1 */
1460  {
1461  tr[tr[t].d0].u0 = t;
1462  tr[tr[t].d0].u1 = -1;
1463  tr[tr[t].d1].u0 = t;
1464  tr[tr[t].d1].u1 = tn;
1465 
1466  /* new code to determine the bottom neighbours of the */
1467  /* newly partitioned trapezoid */
1468 
1469  tr[tn].d0 = tr[t].d1;
1470  tr[tn].d1 = -1;
1471 
1472  tnext = tr[t].d1;
1473  }
1474 
1475  t = tnext;
1476  }
1477 
1478  tr[tn_sav].lseg = segnum;
1479  tr[t_sav].rseg = segnum;
1480  } /* end-while */
1481 
1482  /* Now combine those trapezoids which share common segments. We can */
1483  /* use the pointers to the parent to connect these together. This */
1484  /* works only because all these new trapezoids have been formed */
1485  /* due to splitting by the segment, and hence have only one parent */
1486 
1487  tfirstl = tfirst;
1488  tlastl = tlast;
1489  merge_trapezoids(segnum, tfirstl, tlastl, S_LEFT);
1490  merge_trapezoids(segnum, tfirstr, tlastr, S_RIGHT);
1491 
1492  seg[segnum].is_inserted = true;
1493  return 0;
1494 }
1495 
1496 
1497 /* Update the roots stored for each of the endpoints of the segment.
1498  * This is done to speed up the location-query for the endpoint when
1499  * the segment is inserted into the trapezoidation subsequently
1500  */
1501 int Triangulator::
1502 find_new_roots(int segnum) {
1503  // cerr << "find_new_roots(" << segnum << ")\n";
1504  segment_t *s = &seg[segnum];
1505 
1506  if (s->is_inserted)
1507  return 0;
1508 
1509  s->root0 = locate_endpoint(&s->v0, &s->v1, s->root0);
1510  s->root0 = tr[s->root0].sink;
1511 
1512  s->root1 = locate_endpoint(&s->v1, &s->v0, s->root1);
1513  s->root1 = tr[s->root1].sink;
1514  return 0;
1515 }
1516 
1517 
1518 /* Main routine to perform trapezoidation */
1519 int Triangulator::
1520 construct_trapezoids(int nseg) {
1521  //cerr << "construct_trapezoids(" << nseg << ")\n";
1522  int i;
1523  int root, h;
1524 
1525  /* Add the first segment and get the query structure and trapezoid */
1526  /* list initialised */
1527 
1528  root = init_query_structure(choose_segment());
1529 
1530  for (i = 1; i <= nseg; i++) {
1531  seg[i].root1 = root;
1532  seg[i].root0 = root;
1533  }
1534 
1535  for (h = 1; h <= math_logstar_n(nseg); h++)
1536  {
1537  for (i = math_N(nseg, h -1) + 1; i <= math_N(nseg, h); i++) {
1538  if (add_segment(choose_segment()) != 0) {
1539  // error in add_segment.
1540  return 1;
1541  }
1542  }
1543 
1544  /* Find a new root for each of the segment endpoints */
1545  for (i = 1; i <= nseg; i++)
1546  find_new_roots(i);
1547  }
1548 
1549  for (i = math_N(nseg, math_logstar_n(nseg)) + 1; i <= nseg; i++)
1550  add_segment(choose_segment());
1551 
1552  return 0;
1553 }
1554 
1555 
1556 
1557 /* Function returns true if the trapezoid lies inside the polygon */
1558 int Triangulator::
1559 inside_polygon(trap_t *t) {
1560  int rseg = t->rseg;
1561 
1562  if (t->state == ST_INVALID)
1563  return 0;
1564 
1565  if ((t->lseg <= 0) || (t->rseg <= 0))
1566  return 0;
1567 
1568  if (((t->u0 <= 0) && (t->u1 <= 0)) ||
1569  ((t->d0 <= 0) && (t->d1 <= 0))) /* triangle */
1570  return (_greater_than(&seg[rseg].v1, &seg[rseg].v0));
1571 
1572  return 0;
1573 }
1574 
1575 
1576 /* return a new mon structure from the table */
1577 int Triangulator::
1578 newmon() {
1579  int index = (int)mon.size();
1580  mon.push_back(0);
1581  // cerr << "newmon " << index << "\n";
1582  return index;
1583 }
1584 
1585 
1586 /* return a new chain element from the table */
1587 int Triangulator::
1588 new_chain_element() {
1589  int index = (int)mchain.size();
1590  mchain.push_back(monchain_t());
1591  // cerr << "new_chain_element " << index << "\n";
1592  return index;
1593 }
1594 
1595 
1596 double Triangulator::
1597 get_angle(point_t *vp0, point_t *vpnext, point_t *vp1) {
1598  point_t v0, v1;
1599 
1600  v0.x = vpnext->x - vp0->x;
1601  v0.y = vpnext->y - vp0->y;
1602 
1603  v1.x = vp1->x - vp0->x;
1604  v1.y = vp1->y - vp0->y;
1605 
1606  if (CROSS_SINE(v0, v1) >= 0) /* sine is positive */
1607  return DOT(v0, v1)/LENGTH(v0)/LENGTH(v1);
1608  else
1609  return (-1.0 * DOT(v0, v1)/LENGTH(v0)/LENGTH(v1) - 2);
1610 }
1611 
1612 
1613 /* (v0, v1) is the new diagonal to be added to the polygon. Find which */
1614 /* chain to use and return the positions of v0 and v1 in p and q */
1615 int Triangulator::
1616 get_vertex_positions(int v0, int v1, int *ip, int *iq) {
1617  vertexchain_t *vp0, *vp1;
1618  int i;
1619  double angle, temp;
1620  int tp = 0, tq = 0;
1621 
1622  vp0 = &vert[v0];
1623  vp1 = &vert[v1];
1624 
1625  /* p is identified as follows. Scan from (v0, v1) rightwards till */
1626  /* you hit the first segment starting from v0. That chain is the */
1627  /* chain of our interest */
1628 
1629  angle = -4.0;
1630  for (i = 0; i < 4; i++)
1631  {
1632  if (vp0->vnext[i] <= 0)
1633  continue;
1634  if ((temp = get_angle(&vp0->pt, &(vert[vp0->vnext[i]].pt),
1635  &vp1->pt)) > angle)
1636  {
1637  angle = temp;
1638  tp = i;
1639  }
1640  }
1641 
1642  *ip = tp;
1643 
1644  /* Do similar actions for q */
1645 
1646  angle = -4.0;
1647  for (i = 0; i < 4; i++)
1648  {
1649  if (vp1->vnext[i] <= 0)
1650  continue;
1651  if ((temp = get_angle(&vp1->pt, &(vert[vp1->vnext[i]].pt),
1652  &vp0->pt)) > angle)
1653  {
1654  angle = temp;
1655  tq = i;
1656  }
1657  }
1658 
1659  *iq = tq;
1660 
1661  return 0;
1662 }
1663 
1664 
1665 /* v0 and v1 are specified in anti-clockwise order with respect to
1666  * the current monotone polygon mcur. Split the current polygon into
1667  * two polygons using the diagonal (v0, v1)
1668  */
1669 int Triangulator::
1670 make_new_monotone_poly(int mcur, int v0, int v1) {
1671  int p, q, ip, iq;
1672  int mnew = newmon();
1673  int i, j, nf0, nf1;
1674  vertexchain_t *vp0, *vp1;
1675 
1676  if (v0 <= 0 || v1 <= 0) {
1677  return -1;
1678  }
1679 
1680  vp0 = &vert[v0];
1681  vp1 = &vert[v1];
1682 
1683  get_vertex_positions(v0, v1, &ip, &iq);
1684 
1685  p = vp0->vpos[ip];
1686  q = vp1->vpos[iq];
1687 
1688  /* At this stage, we have got the positions of v0 and v1 in the */
1689  /* desired chain. Now modify the linked lists */
1690 
1691  i = new_chain_element(); /* for the new list */
1692  j = new_chain_element();
1693 
1694  mchain[i].vnum = v0;
1695  mchain[j].vnum = v1;
1696 
1697  mchain[i].next = mchain[p].next;
1698  mchain[mchain[p].next].prev = i;
1699  mchain[i].prev = j;
1700  mchain[j].next = i;
1701  mchain[j].prev = mchain[q].prev;
1702  mchain[mchain[q].prev].next = j;
1703 
1704  mchain[p].next = q;
1705  mchain[q].prev = p;
1706 
1707  nf0 = vp0->nextfree;
1708  nf1 = vp1->nextfree;
1709 
1710  vp0->vnext[ip] = v1;
1711 
1712  vp0->vpos[nf0] = i;
1713  vp0->vnext[nf0] = mchain[mchain[i].next].vnum;
1714  vp1->vpos[nf1] = j;
1715  vp1->vnext[nf1] = v0;
1716 
1717  vp0->nextfree++;
1718  vp1->nextfree++;
1719 
1720 #ifdef DEBUG
1721  fprintf(stderr, "make_poly: mcur = %d, (v0, v1) = (%d, %d)\n",
1722  mcur, v0, v1);
1723  fprintf(stderr, "next posns = (p, q) = (%d, %d)\n", p, q);
1724 #endif
1725 
1726  mon[mcur] = p;
1727  mon[mnew] = i;
1728  return mnew;
1729 }
1730 
1731 /* Main routine to get monotone polygons from the trapezoidation of
1732  * the polygon.
1733  */
1734 
1735 int Triangulator::
1736 monotonate_trapezoids(int n) {
1737  int i;
1738  int tr_start;
1739 
1740  vert.clear();
1741  visited.clear();
1742  mchain.clear();
1743  mon.clear();
1744 
1745  vert.insert(vert.begin(), n + 1, vertexchain_t());
1746  mchain.insert(mchain.begin(), n + 1, monchain_t());
1747 
1748  visited.insert(visited.begin(), tr.size(), 0);
1749 
1750  /* First locate a trapezoid which lies inside the polygon */
1751  /* and which is triangular */
1752  for (i = 1; i < (int)tr.size(); i++)
1753  if (inside_polygon(&tr[i]))
1754  break;
1755  if (i >= (int)tr.size()) {
1756  // No valid trapezoids.
1757  return 0;
1758  }
1759  // printf("start = %d\n", i);
1760  tr_start = i;
1761 
1762  /* Initialise the mon data-structure and start spanning all the */
1763  /* trapezoids within the polygon */
1764 
1765 #if 0
1766  for (i = 1; i <= n; i++)
1767  {
1768  mchain[i].prev = i - 1;
1769  mchain[i].next = i + 1;
1770  mchain[i].vnum = i;
1771  vert[i].pt = seg[i].v0;
1772  vert[i].vnext[0] = i + 1; /* next vertex */
1773  vert[i].vpos[0] = i; /* locn. of next vertex */
1774  vert[i].nextfree = 1;
1775  vert[i].user_i = seg[i].v0_i;
1776  }
1777  mchain[1].prev = n;
1778  mchain[n].next = 1;
1779  vert[n].vnext[0] = 1;
1780  vert[n].vpos[0] = n;
1781  mon.push_back(1); /* position of any vertex in the first */
1782  /* chain */
1783 
1784 #else
1785 
1786  for (i = 1; i <= n; i++)
1787  {
1788  mchain[i].prev = seg[i].prev;
1789  mchain[i].next = seg[i].next;
1790  mchain[i].vnum = i;
1791  vert[i].pt = seg[i].v0;
1792  vert[i].vnext[0] = seg[i].next; /* next vertex */
1793  vert[i].vpos[0] = i; /* locn. of next vertex */
1794  vert[i].nextfree = 1;
1795  vert[i].user_i = seg[i].v0_i;
1796  }
1797 
1798  mon.push_back(1); /* position of any vertex in the first */
1799  /* chain */
1800 
1801 #endif
1802 
1803  /* traverse the polygon */
1804  if (tr[tr_start].u0 > 0)
1805  traverse_polygon(0, tr_start, tr[tr_start].u0, TR_FROM_UP);
1806  else if (tr[tr_start].d0 > 0)
1807  traverse_polygon(0, tr_start, tr[tr_start].d0, TR_FROM_DN);
1808 
1809  /* return the number of polygons created */
1810  return newmon();
1811 }
1812 
1813 
1814 /* recursively visit all the trapezoids */
1815 int Triangulator::
1816 traverse_polygon(int mcur, int trnum, int from, int dir) {
1817  // printf("traverse_polygon(%d, %d, %d, %d)\n", mcur, trnum, from, dir);
1818 
1819  if (mcur < 0 || trnum <= 0)
1820  return 0;
1821 
1822  if (visited[trnum])
1823  return 0;
1824 
1825  trap_t *t = &tr[trnum];
1826  // int howsplit;
1827  int mnew;
1828  int v0, v1; //, v0next, v1next;
1829  int retval = 0; //, tmp;
1830  int do_switch = false;
1831 
1832  // printf("visited size = %d, visited[trnum] = %d\n", visited.size(), visited[trnum]);
1833 
1834  visited[trnum] = true;
1835 
1836  /* We have much more information available here. */
1837  /* rseg: goes upwards */
1838  /* lseg: goes downwards */
1839 
1840  /* Initially assume that dir = TR_FROM_DN (from the left) */
1841  /* Switch v0 and v1 if necessary afterwards */
1842 
1843 
1844  /* special cases for triangles with cusps at the opposite ends. */
1845  /* take care of this first */
1846  if ((t->u0 <= 0) && (t->u1 <= 0))
1847  {
1848  if ((t->d0 > 0) && (t->d1 > 0)) /* downward opening triangle */
1849  {
1850  v0 = tr[t->d1].lseg;
1851  v1 = t->lseg;
1852  if (from == t->d1)
1853  {
1854  do_switch = true;
1855  mnew = make_new_monotone_poly(mcur, v1, v0);
1856  traverse_polygon(mcur, t->d1, trnum, TR_FROM_UP);
1857  traverse_polygon(mnew, t->d0, trnum, TR_FROM_UP);
1858  }
1859  else
1860  {
1861  mnew = make_new_monotone_poly(mcur, v0, v1);
1862  traverse_polygon(mcur, t->d0, trnum, TR_FROM_UP);
1863  traverse_polygon(mnew, t->d1, trnum, TR_FROM_UP);
1864  }
1865  }
1866  else
1867  {
1868  retval = SP_NOSPLIT; /* Just traverse all neighbours */
1869  traverse_polygon(mcur, t->u0, trnum, TR_FROM_DN);
1870  traverse_polygon(mcur, t->u1, trnum, TR_FROM_DN);
1871  traverse_polygon(mcur, t->d0, trnum, TR_FROM_UP);
1872  traverse_polygon(mcur, t->d1, trnum, TR_FROM_UP);
1873  }
1874  }
1875 
1876  else if ((t->d0 <= 0) && (t->d1 <= 0))
1877  {
1878  if ((t->u0 > 0) && (t->u1 > 0)) /* upward opening triangle */
1879  {
1880  v0 = t->rseg;
1881  v1 = tr[t->u0].rseg;
1882  if (from == t->u1)
1883  {
1884  do_switch = true;
1885  mnew = make_new_monotone_poly(mcur, v1, v0);
1886  traverse_polygon(mcur, t->u1, trnum, TR_FROM_DN);
1887  traverse_polygon(mnew, t->u0, trnum, TR_FROM_DN);
1888  }
1889  else
1890  {
1891  mnew = make_new_monotone_poly(mcur, v0, v1);
1892  traverse_polygon(mcur, t->u0, trnum, TR_FROM_DN);
1893  traverse_polygon(mnew, t->u1, trnum, TR_FROM_DN);
1894  }
1895  }
1896  else
1897  {
1898  retval = SP_NOSPLIT; /* Just traverse all neighbours */
1899  traverse_polygon(mcur, t->u0, trnum, TR_FROM_DN);
1900  traverse_polygon(mcur, t->u1, trnum, TR_FROM_DN);
1901  traverse_polygon(mcur, t->d0, trnum, TR_FROM_UP);
1902  traverse_polygon(mcur, t->d1, trnum, TR_FROM_UP);
1903  }
1904  }
1905 
1906  else if ((t->u0 > 0) && (t->u1 > 0))
1907  {
1908  if ((t->d0 > 0) && (t->d1 > 0)) /* downward + upward cusps */
1909  {
1910  v0 = tr[t->d1].lseg;
1911  v1 = tr[t->u0].rseg;
1912  retval = SP_2UP_2DN;
1913  if (((dir == TR_FROM_DN) && (t->d1 == from)) ||
1914  ((dir == TR_FROM_UP) && (t->u1 == from)))
1915  {
1916  do_switch = true;
1917  mnew = make_new_monotone_poly(mcur, v1, v0);
1918  traverse_polygon(mcur, t->u1, trnum, TR_FROM_DN);
1919  traverse_polygon(mcur, t->d1, trnum, TR_FROM_UP);
1920  traverse_polygon(mnew, t->u0, trnum, TR_FROM_DN);
1921  traverse_polygon(mnew, t->d0, trnum, TR_FROM_UP);
1922  }
1923  else
1924  {
1925  mnew = make_new_monotone_poly(mcur, v0, v1);
1926  traverse_polygon(mcur, t->u0, trnum, TR_FROM_DN);
1927  traverse_polygon(mcur, t->d0, trnum, TR_FROM_UP);
1928  traverse_polygon(mnew, t->u1, trnum, TR_FROM_DN);
1929  traverse_polygon(mnew, t->d1, trnum, TR_FROM_UP);
1930  }
1931  }
1932  else /* only downward cusp */
1933  {
1934  if (_equal_to(&t->lo, &seg[t->lseg].v1))
1935  {
1936  v0 = tr[t->u0].rseg;
1937  v1 = seg[t->lseg].next;
1938 
1939  retval = SP_2UP_LEFT;
1940  if ((dir == TR_FROM_UP) && (t->u0 == from))
1941  {
1942  do_switch = true;
1943  mnew = make_new_monotone_poly(mcur, v1, v0);
1944  traverse_polygon(mcur, t->u0, trnum, TR_FROM_DN);
1945  traverse_polygon(mnew, t->d0, trnum, TR_FROM_UP);
1946  traverse_polygon(mnew, t->u1, trnum, TR_FROM_DN);
1947  traverse_polygon(mnew, t->d1, trnum, TR_FROM_UP);
1948  }
1949  else
1950  {
1951  mnew = make_new_monotone_poly(mcur, v0, v1);
1952  traverse_polygon(mcur, t->u1, trnum, TR_FROM_DN);
1953  traverse_polygon(mcur, t->d0, trnum, TR_FROM_UP);
1954  traverse_polygon(mcur, t->d1, trnum, TR_FROM_UP);
1955  traverse_polygon(mnew, t->u0, trnum, TR_FROM_DN);
1956  }
1957  }
1958  else
1959  {
1960  v0 = t->rseg;
1961  v1 = tr[t->u0].rseg;
1962  retval = SP_2UP_RIGHT;
1963  if ((dir == TR_FROM_UP) && (t->u1 == from))
1964  {
1965  do_switch = true;
1966  mnew = make_new_monotone_poly(mcur, v1, v0);
1967  traverse_polygon(mcur, t->u1, trnum, TR_FROM_DN);
1968  traverse_polygon(mnew, t->d1, trnum, TR_FROM_UP);
1969  traverse_polygon(mnew, t->d0, trnum, TR_FROM_UP);
1970  traverse_polygon(mnew, t->u0, trnum, TR_FROM_DN);
1971  }
1972  else
1973  {
1974  mnew = make_new_monotone_poly(mcur, v0, v1);
1975  traverse_polygon(mcur, t->u0, trnum, TR_FROM_DN);
1976  traverse_polygon(mcur, t->d0, trnum, TR_FROM_UP);
1977  traverse_polygon(mcur, t->d1, trnum, TR_FROM_UP);
1978  traverse_polygon(mnew, t->u1, trnum, TR_FROM_DN);
1979  }
1980  }
1981  }
1982  }
1983  else if ((t->u0 > 0) || (t->u1 > 0)) /* no downward cusp */
1984  {
1985  if ((t->d0 > 0) && (t->d1 > 0)) /* only upward cusp */
1986  {
1987  if (_equal_to(&t->hi, &seg[t->lseg].v0))
1988  {
1989  v0 = tr[t->d1].lseg;
1990  v1 = t->lseg;
1991  retval = SP_2DN_LEFT;
1992  if (!((dir == TR_FROM_DN) && (t->d0 == from)))
1993  {
1994  do_switch = true;
1995  mnew = make_new_monotone_poly(mcur, v1, v0);
1996  traverse_polygon(mcur, t->u1, trnum, TR_FROM_DN);
1997  traverse_polygon(mcur, t->d1, trnum, TR_FROM_UP);
1998  traverse_polygon(mcur, t->u0, trnum, TR_FROM_DN);
1999  traverse_polygon(mnew, t->d0, trnum, TR_FROM_UP);
2000  }
2001  else
2002  {
2003  mnew = make_new_monotone_poly(mcur, v0, v1);
2004  traverse_polygon(mcur, t->d0, trnum, TR_FROM_UP);
2005  traverse_polygon(mnew, t->u0, trnum, TR_FROM_DN);
2006  traverse_polygon(mnew, t->u1, trnum, TR_FROM_DN);
2007  traverse_polygon(mnew, t->d1, trnum, TR_FROM_UP);
2008  }
2009  }
2010  else
2011  {
2012  v0 = tr[t->d1].lseg;
2013  v1 = seg[t->rseg].next;
2014 
2015  retval = SP_2DN_RIGHT;
2016  if ((dir == TR_FROM_DN) && (t->d1 == from))
2017  {
2018  do_switch = true;
2019  mnew = make_new_monotone_poly(mcur, v1, v0);
2020  traverse_polygon(mcur, t->d1, trnum, TR_FROM_UP);
2021  traverse_polygon(mnew, t->u1, trnum, TR_FROM_DN);
2022  traverse_polygon(mnew, t->u0, trnum, TR_FROM_DN);
2023  traverse_polygon(mnew, t->d0, trnum, TR_FROM_UP);
2024  }
2025  else
2026  {
2027  mnew = make_new_monotone_poly(mcur, v0, v1);
2028  traverse_polygon(mcur, t->u0, trnum, TR_FROM_DN);
2029  traverse_polygon(mcur, t->d0, trnum, TR_FROM_UP);
2030  traverse_polygon(mcur, t->u1, trnum, TR_FROM_DN);
2031  traverse_polygon(mnew, t->d1, trnum, TR_FROM_UP);
2032  }
2033  }
2034  }
2035  else /* no cusp */
2036  {
2037  if (_equal_to(&t->hi, &seg[t->lseg].v0) &&
2038  _equal_to(&t->lo, &seg[t->rseg].v0))
2039  {
2040  v0 = t->rseg;
2041  v1 = t->lseg;
2042  retval = SP_SIMPLE_LRDN;
2043  if (dir == TR_FROM_UP)
2044  {
2045  do_switch = true;
2046  mnew = make_new_monotone_poly(mcur, v1, v0);
2047  traverse_polygon(mcur, t->u0, trnum, TR_FROM_DN);
2048  traverse_polygon(mcur, t->u1, trnum, TR_FROM_DN);
2049  traverse_polygon(mnew, t->d1, trnum, TR_FROM_UP);
2050  traverse_polygon(mnew, t->d0, trnum, TR_FROM_UP);
2051  }
2052  else
2053  {
2054  mnew = make_new_monotone_poly(mcur, v0, v1);
2055  traverse_polygon(mcur, t->d1, trnum, TR_FROM_UP);
2056  traverse_polygon(mcur, t->d0, trnum, TR_FROM_UP);
2057  traverse_polygon(mnew, t->u0, trnum, TR_FROM_DN);
2058  traverse_polygon(mnew, t->u1, trnum, TR_FROM_DN);
2059  }
2060  }
2061  else if (_equal_to(&t->hi, &seg[t->rseg].v1) &&
2062  _equal_to(&t->lo, &seg[t->lseg].v1))
2063  {
2064  v0 = seg[t->rseg].next;
2065  v1 = seg[t->lseg].next;
2066 
2067  retval = SP_SIMPLE_LRUP;
2068  if (dir == TR_FROM_UP)
2069  {
2070  do_switch = true;
2071  mnew = make_new_monotone_poly(mcur, v1, v0);
2072  traverse_polygon(mcur, t->u0, trnum, TR_FROM_DN);
2073  traverse_polygon(mcur, t->u1, trnum, TR_FROM_DN);
2074  traverse_polygon(mnew, t->d1, trnum, TR_FROM_UP);
2075  traverse_polygon(mnew, t->d0, trnum, TR_FROM_UP);
2076  }
2077  else
2078  {
2079  mnew = make_new_monotone_poly(mcur, v0, v1);
2080  traverse_polygon(mcur, t->d1, trnum, TR_FROM_UP);
2081  traverse_polygon(mcur, t->d0, trnum, TR_FROM_UP);
2082  traverse_polygon(mnew, t->u0, trnum, TR_FROM_DN);
2083  traverse_polygon(mnew, t->u1, trnum, TR_FROM_DN);
2084  }
2085  }
2086  else /* no split possible */
2087  {
2088  retval = SP_NOSPLIT;
2089  traverse_polygon(mcur, t->u0, trnum, TR_FROM_DN);
2090  traverse_polygon(mcur, t->d0, trnum, TR_FROM_UP);
2091  traverse_polygon(mcur, t->u1, trnum, TR_FROM_DN);
2092  traverse_polygon(mcur, t->d1, trnum, TR_FROM_UP);
2093  }
2094  }
2095  }
2096 
2097  return retval;
2098 }
2099 
2100 
2101 /* For each monotone polygon, find the ymax and ymin (to determine the */
2102 /* two y-monotone chains) and pass on this monotone polygon for greedy */
2103 /* triangulation. */
2104 /* Take care not to triangulate duplicate monotone polygons */
2105 
2106 void Triangulator::
2107 triangulate_monotone_polygons(int nvert, int nmonpoly) {
2108  int i;
2109  point_t ymax, ymin;
2110  int p, vfirst, posmax, posmin, v;
2111  int vcount, processed;
2112 
2113  #ifdef DEBUG
2114  for (i = 0; i < nmonpoly; i++)
2115  {
2116  fprintf(stderr, "\n\nPolygon %d: ", i);
2117  vfirst = mchain[mon[i]].vnum;
2118  p = mchain[mon[i]].next;
2119  fprintf (stderr, "%d ", mchain[mon[i]].vnum);
2120  while (mchain[p].vnum != vfirst)
2121  {
2122  fprintf(stderr, "%d ", mchain[p].vnum);
2123  if (mchain[p].vnum == -1) {
2124  fprintf(stderr, " xx");
2125  break;
2126  }
2127  p = mchain[p].next;
2128  }
2129  }
2130  fprintf(stderr, "\n");
2131  #endif
2132 
2133  for (i = 0; i < nmonpoly; i++)
2134  {
2135  vcount = 1;
2136  processed = false;
2137  vfirst = mchain[mon[i]].vnum;
2138  if (vfirst <= 0) {
2139  return;
2140  }
2141  ymin = vert[vfirst].pt;
2142  ymax = ymin;
2143  posmin = mon[i];
2144  posmax = posmin;
2145  mchain[mon[i]].marked = true;
2146  p = mchain[mon[i]].next;
2147  while ((v = mchain[p].vnum) != vfirst)
2148  {
2149  if (v <= 0) {
2150  return;
2151  }
2152  if (mchain[p].marked)
2153  {
2154  processed = true;
2155  break; /* break from while */
2156  }
2157  else
2158  mchain[p].marked = true;
2159 
2160  if (_greater_than(&vert[v].pt, &ymax))
2161  {
2162  ymax = vert[v].pt;
2163  posmax = p;
2164  }
2165  if (_less_than(&vert[v].pt, &ymin))
2166  {
2167  ymin = vert[v].pt;
2168  posmin = p;
2169  }
2170  p = mchain[p].next;
2171  vcount++;
2172  }
2173 
2174  if (processed) /* Go to next polygon */
2175  continue;
2176 
2177  if (vcount == 3) /* already a triangle */
2178  {
2179  _result.push_back(Triangle(this, mchain[p].vnum,
2180  mchain[mchain[p].next].vnum,
2181  mchain[mchain[p].prev].vnum));
2182  }
2183  else /* triangulate the polygon */
2184  {
2185  v = mchain[mchain[posmax].next].vnum;
2186  if (_equal_to(&vert[v].pt, &ymin))
2187  { /* LHS is a single line */
2188  triangulate_single_polygon(nvert, posmax, TRI_LHS);
2189  }
2190  else
2191  triangulate_single_polygon(nvert, posmax, TRI_RHS);
2192  }
2193  }
2194 }
2195 
2196 
2197 /* A greedy corner-cutting algorithm to triangulate a y-monotone
2198  * polygon in O(n) time.
2199  * Joseph O-Rourke, Computational Geometry in C.
2200  */
2201 void Triangulator::
2202 triangulate_single_polygon(int nvert, int posmax, int side) {
2203  int v;
2204  vector_int rc; /* reflex chain */
2205  int ri;
2206  int endv, tmp, vpos;
2207 
2208  //cerr << "triangulate_single_polygon(" << nvert << ", " << posmax << ", " << side << ")\n";
2209 
2210  if (side == TRI_RHS) /* RHS segment is a single segment */
2211  {
2212  rc.push_back(mchain[posmax].vnum);
2213  tmp = mchain[posmax].next;
2214  rc.push_back(mchain[tmp].vnum);
2215  ri = 1;
2216 
2217  vpos = mchain[tmp].next;
2218  v = mchain[vpos].vnum;
2219 
2220  if ((endv = mchain[mchain[posmax].prev].vnum) == 0)
2221  endv = nvert;
2222  }
2223  else /* LHS is a single segment */
2224  {
2225  tmp = mchain[posmax].next;
2226  rc.push_back(mchain[tmp].vnum);
2227  tmp = mchain[tmp].next;
2228  rc.push_back(mchain[tmp].vnum);
2229  ri = 1;
2230 
2231  vpos = mchain[tmp].next;
2232  v = mchain[vpos].vnum;
2233 
2234  endv = mchain[posmax].vnum;
2235  }
2236 
2237  while ((v != endv) || (ri > 1))
2238  {
2239  //cerr << " v = " << v << " ri = " << ri << " rc = " << rc.size() << " _result = " << _result.size() << "\n";
2240  if (v <= 0) {
2241  // Something went wrong.
2242  return;
2243  }
2244  if (ri > 0) /* reflex chain is non-empty */
2245  {
2246  PN_stdfloat crossResult = CROSS(vert[v].pt, vert[rc[ri - 1]].pt,
2247  vert[rc[ri]].pt);
2248  if ( crossResult >= 0 ) /* could be convex corner or straight */
2249  {
2250  if ( crossResult > 0) /* convex corner: cut it off */
2251  _result.push_back(Triangle(this, rc[ri - 1], rc[ri], v));
2252  /* else : perfectly straight, will be abandoned anyway */
2253  ri--;
2254  rc.pop_back();
2255  nassertv(ri + 1 == (int)rc.size());
2256  }
2257  else /* concave, add v to the chain */
2258  {
2259  ri++;
2260  rc.push_back(v);
2261  nassertv(ri + 1 == (int)rc.size());
2262  vpos = mchain[vpos].next;
2263  v = mchain[vpos].vnum;
2264  }
2265  }
2266  else /* reflex-chain empty: add v to the */
2267  { /* reflex chain and advance it */
2268  ri++;
2269  rc.push_back(v);
2270  nassertv(ri + 1 == (int)rc.size());
2271  vpos = mchain[vpos].next;
2272  v = mchain[vpos].vnum;
2273  }
2274  } /* end-while */
2275 
2276  /* reached the bottom vertex. Add in the triangle formed */
2277  _result.push_back(Triangle(this, rc[ri - 1], rc[ri], v));
2278  ri--;
2279 }
2280 
2281 
int get_triangle_v0(int n) const
Returns vertex 0 of the nth triangle generated by the previous call to triangulate().
void triangulate()
Does the work of triangulating the specified polygon.
int get_triangle_v1(int n) const
Returns vertex 1 of the nth triangle generated by the previous call to triangulate().
void begin_hole()
Finishes the previous hole, if any, and prepares to add a new hole.
This is a two-component point in space.
Definition: lpoint2.h:411
int get_triangle_v2(int n) const
Returns vertex 2 of the nth triangle generated by the previous call to triangulate().
int get_num_triangles() const
Returns the number of triangles generated by the previous call to triangulate().
void add_hole_vertex(int index)
Adds the next consecutive vertex of the current hole.
void clear_polygon()
Removes the current polygon definition (and its set of holes), but does not clear the vertex pool...
int random_int(int range)
Returns a random integer in the range [0, range).
Definition: randomizer.I:55
void clear()
Removes all vertices and polygon specifications from the Triangulator, and prepares it to start over...
void add_polygon_vertex(int index)
Adds the next consecutive vertex of the polygon.
A handy class to return random numbers.
Definition: randomizer.h:28
int add_vertex(const LPoint2d &point)
Adds a new vertex to the vertex pool.