Panda3D
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 
void triangulate()
Does the work of triangulating the specified polygon.
int get_triangle_v0(int n) const
Returns vertex 0 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:424
int get_triangle_v2(int n) const
Returns vertex 2 of the nth triangle 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
int get_num_triangles() const
Returns the number of triangles 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.
A handy class to return random numbers.
Definition: randomizer.h:28
int get_triangle_v1(int n) const
Returns vertex 1 of the nth triangle generated by the previous call to triangulate().
int add_vertex(const LPoint2d &point)
Adds a new vertex to the vertex pool.