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