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