1 //-----------------------------------------------------------------------------
2 // Operations on polygons (planar, of line segment edges).
3 //
4 // Copyright 2008-2013 Jonathan Westhues.
5 //-----------------------------------------------------------------------------
6 #include "solvespace.h"
7
Normal(void)8 Vector STriangle::Normal(void) {
9 Vector ab = b.Minus(a), bc = c.Minus(b);
10 return ab.Cross(bc);
11 }
12
MinAltitude(void)13 double STriangle::MinAltitude(void) {
14 double altA = a.DistanceToLine(b, c.Minus(b)),
15 altB = b.DistanceToLine(c, a.Minus(c)),
16 altC = c.DistanceToLine(a, b.Minus(a));
17
18 return min(altA, min(altB, altC));
19 }
20
ContainsPoint(Vector p)21 bool STriangle::ContainsPoint(Vector p) {
22 Vector n = Normal();
23 if(MinAltitude() < LENGTH_EPS) {
24 // shouldn't happen; zero-area triangle
25 return false;
26 }
27 return ContainsPointProjd(n.WithMagnitude(1), p);
28 }
29
ContainsPointProjd(Vector n,Vector p)30 bool STriangle::ContainsPointProjd(Vector n, Vector p) {
31 Vector ab = b.Minus(a), bc = c.Minus(b), ca = a.Minus(c);
32
33 Vector no_ab = n.Cross(ab);
34 if(no_ab.Dot(p) < no_ab.Dot(a) - LENGTH_EPS) return false;
35
36 Vector no_bc = n.Cross(bc);
37 if(no_bc.Dot(p) < no_bc.Dot(b) - LENGTH_EPS) return false;
38
39 Vector no_ca = n.Cross(ca);
40 if(no_ca.Dot(p) < no_ca.Dot(c) - LENGTH_EPS) return false;
41
42 return true;
43 }
44
FlipNormal(void)45 void STriangle::FlipNormal(void) {
46 swap(a, b);
47 swap(an, bn);
48 }
49
From(STriMeta meta,Vector a,Vector b,Vector c)50 STriangle STriangle::From(STriMeta meta, Vector a, Vector b, Vector c) {
51 STriangle tr = {};
52 tr.meta = meta;
53 tr.a = a;
54 tr.b = b;
55 tr.c = c;
56 return tr;
57 }
58
From(Vector a,Vector b)59 SEdge SEdge::From(Vector a, Vector b) {
60 SEdge se = {};
61 se.a = a;
62 se.b = b;
63 return se;
64 }
65
EdgeCrosses(Vector ea,Vector eb,Vector * ppi,SPointList * spl)66 bool SEdge::EdgeCrosses(Vector ea, Vector eb, Vector *ppi, SPointList *spl) {
67 Vector d = eb.Minus(ea);
68 double t_eps = LENGTH_EPS/d.Magnitude();
69
70 double t, tthis;
71 bool skew;
72 Vector pi;
73 bool inOrEdge0, inOrEdge1;
74
75 Vector dthis = b.Minus(a);
76 double tthis_eps = LENGTH_EPS/dthis.Magnitude();
77
78 if((ea.Equals(a) && eb.Equals(b)) ||
79 (eb.Equals(a) && ea.Equals(b)))
80 {
81 if(ppi) *ppi = a;
82 if(spl) spl->Add(a);
83 return true;
84 }
85
86 // Can't just test if distance between d and a equals distance between d and b;
87 // they could be on opposite sides, since we don't have the sign.
88 double m = sqrt(d.Magnitude()*dthis.Magnitude());
89 if(sqrt(fabs(d.Dot(dthis))) > (m - LENGTH_EPS)) {
90 // The edges are parallel.
91 if(fabs(a.DistanceToLine(ea, d)) > LENGTH_EPS) {
92 // and not coincident, so can't be interesecting
93 return false;
94 }
95 // The edges are coincident. Make sure that neither endpoint lies
96 // on the other
97 bool inters = false;
98 double t;
99 t = a.Minus(ea).DivPivoting(d);
100 if(t > t_eps && t < (1 - t_eps)) inters = true;
101 t = b.Minus(ea).DivPivoting(d);
102 if(t > t_eps && t < (1 - t_eps)) inters = true;
103 t = ea.Minus(a).DivPivoting(dthis);
104 if(t > tthis_eps && t < (1 - tthis_eps)) inters = true;
105 t = eb.Minus(a).DivPivoting(dthis);
106 if(t > tthis_eps && t < (1 - tthis_eps)) inters = true;
107
108 if(inters) {
109 if(ppi) *ppi = a;
110 if(spl) spl->Add(a);
111 return true;
112 } else {
113 // So coincident but disjoint, okay.
114 return false;
115 }
116 }
117
118 // Lines are not parallel, so look for an intersection.
119 pi = Vector::AtIntersectionOfLines(ea, eb, a, b,
120 &skew,
121 &t, &tthis);
122 if(skew) return false;
123
124 inOrEdge0 = (t > -t_eps) && (t < (1 + t_eps));
125 inOrEdge1 = (tthis > -tthis_eps) && (tthis < (1 + tthis_eps));
126
127 if(inOrEdge0 && inOrEdge1) {
128 if(a.Equals(ea) || b.Equals(ea) ||
129 a.Equals(eb) || b.Equals(eb))
130 {
131 // Not an intersection if we share an endpoint with an edge
132 return false;
133 }
134 // But it's an intersection if a vertex of one edge lies on the
135 // inside of the other (or if they cross away from either's
136 // vertex).
137 if(ppi) *ppi = pi;
138 if(spl) spl->Add(pi);
139 return true;
140 }
141 return false;
142 }
143
Clear(void)144 void SEdgeList::Clear(void) {
145 l.Clear();
146 }
147
AddEdge(Vector a,Vector b,int auxA,int auxB)148 void SEdgeList::AddEdge(Vector a, Vector b, int auxA, int auxB) {
149 SEdge e = {};
150 e.a = a;
151 e.b = b;
152 e.auxA = auxA;
153 e.auxB = auxB;
154 l.Add(&e);
155 }
156
AssembleContour(Vector first,Vector last,SContour * dest,SEdge * errorAt,bool keepDir)157 bool SEdgeList::AssembleContour(Vector first, Vector last, SContour *dest,
158 SEdge *errorAt, bool keepDir)
159 {
160 int i;
161
162 dest->AddPoint(first);
163 dest->AddPoint(last);
164
165 do {
166 for(i = 0; i < l.n; i++) {
167 SEdge *se = &(l.elem[i]);
168 if(se->tag) continue;
169
170 if(se->a.Equals(last)) {
171 dest->AddPoint(se->b);
172 last = se->b;
173 se->tag = 1;
174 break;
175 }
176 // Don't allow backwards edges if keepDir is true.
177 if(!keepDir && se->b.Equals(last)) {
178 dest->AddPoint(se->a);
179 last = se->a;
180 se->tag = 1;
181 break;
182 }
183 }
184 if(i >= l.n) {
185 // Couldn't assemble a closed contour; mark where.
186 if(errorAt) {
187 errorAt->a = first;
188 errorAt->b = last;
189 }
190 return false;
191 }
192
193 } while(!last.Equals(first));
194
195 return true;
196 }
197
AssemblePolygon(SPolygon * dest,SEdge * errorAt,bool keepDir)198 bool SEdgeList::AssemblePolygon(SPolygon *dest, SEdge *errorAt, bool keepDir) {
199 dest->Clear();
200
201 bool allClosed = true;
202 for(;;) {
203 Vector first = Vector::From(0, 0, 0);
204 Vector last = Vector::From(0, 0, 0);
205 int i;
206 for(i = 0; i < l.n; i++) {
207 if(!l.elem[i].tag) {
208 first = l.elem[i].a;
209 last = l.elem[i].b;
210 l.elem[i].tag = 1;
211 break;
212 }
213 }
214 if(i >= l.n) {
215 return allClosed;
216 }
217
218 // Create a new empty contour in our polygon, and finish assembling
219 // into that contour.
220 dest->AddEmptyContour();
221 if(!AssembleContour(first, last, &(dest->l.elem[dest->l.n-1]),
222 errorAt, keepDir))
223 {
224 allClosed = false;
225 }
226 // But continue assembling, even if some of the contours are open
227 }
228 }
229
230 //-----------------------------------------------------------------------------
231 // Test if the specified edge crosses any of the edges in our list. Two edges
232 // are not considered to cross if they share an endpoint (within LENGTH_EPS),
233 // but they are considered to cross if they are coincident and overlapping.
234 // If pi is not NULL, then a crossing is returned in that.
235 //-----------------------------------------------------------------------------
AnyEdgeCrossings(Vector a,Vector b,Vector * ppi,SPointList * spl)236 int SEdgeList::AnyEdgeCrossings(Vector a, Vector b,
237 Vector *ppi, SPointList *spl)
238 {
239 int cnt = 0;
240 SEdge *se;
241 for(se = l.First(); se; se = l.NextAfter(se)) {
242 if(se->EdgeCrosses(a, b, ppi, spl)) {
243 cnt++;
244 }
245 }
246 return cnt;
247 }
248
249 //-----------------------------------------------------------------------------
250 // Returns true if the intersecting edge list contains an edge that shares
251 // an endpoint with one of our edges.
252 //-----------------------------------------------------------------------------
ContainsEdgeFrom(SEdgeList * sel)253 bool SEdgeList::ContainsEdgeFrom(SEdgeList *sel) {
254 SEdge *se;
255 for(se = l.First(); se; se = l.NextAfter(se)) {
256 if(sel->ContainsEdge(se)) return true;
257 }
258 return false;
259 }
ContainsEdge(SEdge * set)260 bool SEdgeList::ContainsEdge(SEdge *set) {
261 SEdge *se;
262 for(se = l.First(); se; se = l.NextAfter(se)) {
263 if((se->a).Equals(set->a) && (se->b).Equals(set->b)) return true;
264 if((se->b).Equals(set->a) && (se->a).Equals(set->b)) return true;
265 }
266 return false;
267 }
268
269 //-----------------------------------------------------------------------------
270 // Remove unnecessary edges: if two are anti-parallel then remove both, and if
271 // two are parallel then remove one.
272 //-----------------------------------------------------------------------------
CullExtraneousEdges(void)273 void SEdgeList::CullExtraneousEdges(void) {
274 l.ClearTags();
275 int i, j;
276 for(i = 0; i < l.n; i++) {
277 SEdge *se = &(l.elem[i]);
278 for(j = i+1; j < l.n; j++) {
279 SEdge *set = &(l.elem[j]);
280 if((set->a).Equals(se->a) && (set->b).Equals(se->b)) {
281 // Two parallel edges exist; so keep only the first one.
282 set->tag = 1;
283 }
284 if((set->a).Equals(se->b) && (set->b).Equals(se->a)) {
285 // Two anti-parallel edges exist; so keep neither.
286 se->tag = 1;
287 set->tag = 1;
288 }
289 }
290 }
291 l.RemoveTagged();
292 }
293
294 //-----------------------------------------------------------------------------
295 // Make a kd-tree of edges. This is used for O(log(n)) implementations of stuff
296 // that would naively be O(n).
297 //-----------------------------------------------------------------------------
Alloc(void)298 SKdNodeEdges *SKdNodeEdges::Alloc(void) {
299 SKdNodeEdges *ne = (SKdNodeEdges *)AllocTemporary(sizeof(SKdNodeEdges));
300 *ne = {};
301 return ne;
302 }
Alloc(void)303 SEdgeLl *SEdgeLl::Alloc(void) {
304 SEdgeLl *sell = (SEdgeLl *)AllocTemporary(sizeof(SEdgeLl));
305 *sell = {};
306 return sell;
307 }
From(SEdgeList * sel)308 SKdNodeEdges *SKdNodeEdges::From(SEdgeList *sel) {
309 SEdgeLl *sell = NULL;
310 SEdge *se;
311 for(se = sel->l.First(); se; se = sel->l.NextAfter(se)) {
312 SEdgeLl *n = SEdgeLl::Alloc();
313 n->se = se;
314 n->next = sell;
315 sell = n;
316 }
317 return SKdNodeEdges::From(sell);
318 }
From(SEdgeLl * sell)319 SKdNodeEdges *SKdNodeEdges::From(SEdgeLl *sell) {
320 SKdNodeEdges *n = SKdNodeEdges::Alloc();
321
322 // Compute the midpoints (just mean, though median would be better) of
323 // each component.
324 Vector ptAve = Vector::From(0, 0, 0);
325 SEdgeLl *flip;
326 int totaln = 0;
327 for(flip = sell; flip; flip = flip->next) {
328 ptAve = ptAve.Plus(flip->se->a);
329 ptAve = ptAve.Plus(flip->se->b);
330 totaln++;
331 }
332 ptAve = ptAve.ScaledBy(1.0 / (2*totaln));
333
334 // For each component, see how it splits.
335 int ltln[3] = { 0, 0, 0 }, gtln[3] = { 0, 0, 0 };
336 double badness[3];
337 for(flip = sell; flip; flip = flip->next) {
338 for(int i = 0; i < 3; i++) {
339 if(flip->se->a.Element(i) < ptAve.Element(i) + KDTREE_EPS ||
340 flip->se->b.Element(i) < ptAve.Element(i) + KDTREE_EPS)
341 {
342 ltln[i]++;
343 }
344 if(flip->se->a.Element(i) > ptAve.Element(i) - KDTREE_EPS ||
345 flip->se->b.Element(i) > ptAve.Element(i) - KDTREE_EPS)
346 {
347 gtln[i]++;
348 }
349 }
350 }
351 for(int i = 0; i < 3; i++) {
352 badness[i] = pow((double)ltln[i], 4) + pow((double)gtln[i], 4);
353 }
354
355 // Choose the least bad coordinate to split along.
356 if(badness[0] < badness[1] && badness[0] < badness[2]) {
357 n->which = 0;
358 } else if(badness[1] < badness[2]) {
359 n->which = 1;
360 } else {
361 n->which = 2;
362 }
363 n->c = ptAve.Element(n->which);
364
365 if(totaln < 3 || totaln == gtln[n->which] || totaln == ltln[n->which]) {
366 n->edges = sell;
367 // and we're a leaf node
368 return n;
369 }
370
371 // Sort the edges according to which side(s) of the split plane they're on.
372 SEdgeLl *gtl = NULL, *ltl = NULL;
373 for(flip = sell; flip; flip = flip->next) {
374 if(flip->se->a.Element(n->which) < n->c + KDTREE_EPS ||
375 flip->se->b.Element(n->which) < n->c + KDTREE_EPS)
376 {
377 SEdgeLl *selln = SEdgeLl::Alloc();
378 selln->se = flip->se;
379 selln->next = ltl;
380 ltl = selln;
381 }
382 if(flip->se->a.Element(n->which) > n->c - KDTREE_EPS ||
383 flip->se->b.Element(n->which) > n->c - KDTREE_EPS)
384 {
385 SEdgeLl *selln = SEdgeLl::Alloc();
386 selln->se = flip->se;
387 selln->next = gtl;
388 gtl = selln;
389 }
390 }
391
392 n->lt = SKdNodeEdges::From(ltl);
393 n->gt = SKdNodeEdges::From(gtl);
394 return n;
395 }
396
AnyEdgeCrossings(Vector a,Vector b,int cnt,Vector * pi,SPointList * spl)397 int SKdNodeEdges::AnyEdgeCrossings(Vector a, Vector b, int cnt,
398 Vector *pi, SPointList *spl)
399 {
400 int inters = 0;
401 if(gt && lt) {
402 if(a.Element(which) < c + KDTREE_EPS ||
403 b.Element(which) < c + KDTREE_EPS)
404 {
405 inters += lt->AnyEdgeCrossings(a, b, cnt, pi, spl);
406 }
407 if(a.Element(which) > c - KDTREE_EPS ||
408 b.Element(which) > c - KDTREE_EPS)
409 {
410 inters += gt->AnyEdgeCrossings(a, b, cnt, pi, spl);
411 }
412 } else {
413 SEdgeLl *sell;
414 for(sell = edges; sell; sell = sell->next) {
415 SEdge *se = sell->se;
416 if(se->tag == cnt) continue;
417 if(se->EdgeCrosses(a, b, pi, spl)) {
418 inters++;
419 }
420 se->tag = cnt;
421 }
422 }
423 return inters;
424 }
425
426 //-----------------------------------------------------------------------------
427 // We have an edge list that contains only collinear edges, maybe with more
428 // splits than necessary. Merge any collinear segments that join.
429 //-----------------------------------------------------------------------------
430 static Vector LineStart, LineDirection;
ByTAlongLine(const void * av,const void * bv)431 static int ByTAlongLine(const void *av, const void *bv)
432 {
433 SEdge *a = (SEdge *)av,
434 *b = (SEdge *)bv;
435
436 double ta = (a->a.Minus(LineStart)).DivPivoting(LineDirection),
437 tb = (b->a.Minus(LineStart)).DivPivoting(LineDirection);
438
439 return (ta > tb) ? 1 : -1;
440 }
MergeCollinearSegments(Vector a,Vector b)441 void SEdgeList::MergeCollinearSegments(Vector a, Vector b) {
442 LineStart = a;
443 LineDirection = b.Minus(a);
444 qsort(l.elem, l.n, sizeof(l.elem[0]), ByTAlongLine);
445
446 l.ClearTags();
447 int i;
448 for(i = 1; i < l.n; i++) {
449 SEdge *prev = &(l.elem[i-1]),
450 *now = &(l.elem[i]);
451
452 if((prev->b).Equals(now->a) && prev->auxA == now->auxA) {
453 // The previous segment joins up to us; so merge it into us.
454 prev->tag = 1;
455 now->a = prev->a;
456 }
457 }
458 l.RemoveTagged();
459 }
460
Clear(void)461 void SPointList::Clear(void) {
462 l.Clear();
463 }
464
ContainsPoint(Vector pt)465 bool SPointList::ContainsPoint(Vector pt) {
466 return (IndexForPoint(pt) >= 0);
467 }
468
IndexForPoint(Vector pt)469 int SPointList::IndexForPoint(Vector pt) {
470 int i;
471 for(i = 0; i < l.n; i++) {
472 SPoint *p = &(l.elem[i]);
473 if(pt.Equals(p->p)) {
474 return i;
475 }
476 }
477 // Not found, so return negative to indicate that.
478 return -1;
479 }
480
IncrementTagFor(Vector pt)481 void SPointList::IncrementTagFor(Vector pt) {
482 SPoint *p;
483 for(p = l.First(); p; p = l.NextAfter(p)) {
484 if(pt.Equals(p->p)) {
485 (p->tag)++;
486 return;
487 }
488 }
489 SPoint pa;
490 pa.p = pt;
491 pa.tag = 1;
492 l.Add(&pa);
493 }
494
Add(Vector pt)495 void SPointList::Add(Vector pt) {
496 SPoint p = {};
497 p.p = pt;
498 l.Add(&p);
499 }
500
AddPoint(Vector p)501 void SContour::AddPoint(Vector p) {
502 SPoint sp;
503 sp.tag = 0;
504 sp.p = p;
505
506 l.Add(&sp);
507 }
508
MakeEdgesInto(SEdgeList * el)509 void SContour::MakeEdgesInto(SEdgeList *el) {
510 int i;
511 for(i = 0; i < (l.n - 1); i++) {
512 el->AddEdge(l.elem[i].p, l.elem[i+1].p);
513 }
514 }
515
CopyInto(SContour * dest)516 void SContour::CopyInto(SContour *dest) {
517 SPoint *sp;
518 for(sp = l.First(); sp; sp = l.NextAfter(sp)) {
519 dest->AddPoint(sp->p);
520 }
521 }
522
FindPointWithMinX(void)523 void SContour::FindPointWithMinX(void) {
524 SPoint *sp;
525 xminPt = Vector::From(1e10, 1e10, 1e10);
526 for(sp = l.First(); sp; sp = l.NextAfter(sp)) {
527 if(sp->p.x < xminPt.x) {
528 xminPt = sp->p;
529 }
530 }
531 }
532
ComputeNormal(void)533 Vector SContour::ComputeNormal(void) {
534 Vector n = Vector::From(0, 0, 0);
535
536 for(int i = 0; i < l.n - 2; i++) {
537 Vector u = (l.elem[i+1].p).Minus(l.elem[i+0].p).WithMagnitude(1);
538 Vector v = (l.elem[i+2].p).Minus(l.elem[i+1].p).WithMagnitude(1);
539 Vector nt = u.Cross(v);
540 if(nt.Magnitude() > n.Magnitude()) {
541 n = nt;
542 }
543 }
544 return n.WithMagnitude(1);
545 }
546
AnyEdgeMidpoint(void)547 Vector SContour::AnyEdgeMidpoint(void) {
548 if(l.n < 2) oops();
549 return ((l.elem[0].p).Plus(l.elem[1].p)).ScaledBy(0.5);
550 }
551
IsClockwiseProjdToNormal(Vector n)552 bool SContour::IsClockwiseProjdToNormal(Vector n) {
553 // Degenerate things might happen as we draw; doesn't really matter
554 // what we do then.
555 if(n.Magnitude() < 0.01) return true;
556
557 return (SignedAreaProjdToNormal(n) < 0);
558 }
559
SignedAreaProjdToNormal(Vector n)560 double SContour::SignedAreaProjdToNormal(Vector n) {
561 // An arbitrary 2d coordinate system that has n as its normal
562 Vector u = n.Normal(0);
563 Vector v = n.Normal(1);
564
565 double area = 0;
566 for(int i = 0; i < (l.n - 1); i++) {
567 double u0 = (l.elem[i ].p).Dot(u);
568 double v0 = (l.elem[i ].p).Dot(v);
569 double u1 = (l.elem[i+1].p).Dot(u);
570 double v1 = (l.elem[i+1].p).Dot(v);
571
572 area += ((v0 + v1)/2)*(u1 - u0);
573 }
574 return area;
575 }
576
ContainsPointProjdToNormal(Vector n,Vector p)577 bool SContour::ContainsPointProjdToNormal(Vector n, Vector p) {
578 Vector u = n.Normal(0);
579 Vector v = n.Normal(1);
580
581 double up = p.Dot(u);
582 double vp = p.Dot(v);
583
584 bool inside = false;
585 for(int i = 0; i < (l.n - 1); i++) {
586 double ua = (l.elem[i ].p).Dot(u);
587 double va = (l.elem[i ].p).Dot(v);
588 // The curve needs to be exactly closed; approximation is death.
589 double ub = (l.elem[(i+1)%(l.n-1)].p).Dot(u);
590 double vb = (l.elem[(i+1)%(l.n-1)].p).Dot(v);
591
592 if ((((va <= vp) && (vp < vb)) ||
593 ((vb <= vp) && (vp < va))) &&
594 (up < (ub - ua) * (vp - va) / (vb - va) + ua))
595 {
596 inside = !inside;
597 }
598 }
599
600 return inside;
601 }
602
Reverse(void)603 void SContour::Reverse(void) {
604 l.Reverse();
605 }
606
607
Clear(void)608 void SPolygon::Clear(void) {
609 int i;
610 for(i = 0; i < l.n; i++) {
611 (l.elem[i]).l.Clear();
612 }
613 l.Clear();
614 }
615
AddEmptyContour(void)616 void SPolygon::AddEmptyContour(void) {
617 SContour c = {};
618 l.Add(&c);
619 }
620
MakeEdgesInto(SEdgeList * el)621 void SPolygon::MakeEdgesInto(SEdgeList *el) {
622 int i;
623 for(i = 0; i < l.n; i++) {
624 (l.elem[i]).MakeEdgesInto(el);
625 }
626 }
627
ComputeNormal(void)628 Vector SPolygon::ComputeNormal(void) {
629 if(l.n < 1) return Vector::From(0, 0, 0);
630 return (l.elem[0]).ComputeNormal();
631 }
632
SignedArea(void)633 double SPolygon::SignedArea(void) {
634 SContour *sc;
635 double area = 0;
636 // This returns the true area only if the contours are all oriented
637 // correctly, with the holes backwards from the outer contours.
638 for(sc = l.First(); sc; sc = l.NextAfter(sc)) {
639 area += sc->SignedAreaProjdToNormal(normal);
640 }
641 return area;
642 }
643
ContainsPoint(Vector p)644 bool SPolygon::ContainsPoint(Vector p) {
645 return (WindingNumberForPoint(p) % 2) == 1;
646 }
647
WindingNumberForPoint(Vector p)648 int SPolygon::WindingNumberForPoint(Vector p) {
649 int winding = 0;
650 int i;
651 for(i = 0; i < l.n; i++) {
652 SContour *sc = &(l.elem[i]);
653 if(sc->ContainsPointProjdToNormal(normal, p)) {
654 winding++;
655 }
656 }
657 return winding;
658 }
659
FixContourDirections(void)660 void SPolygon::FixContourDirections(void) {
661 // At output, the contour's tag will be 1 if we reversed it, else 0.
662 l.ClearTags();
663
664 // Outside curve looks counterclockwise, projected against our normal.
665 int i, j;
666 for(i = 0; i < l.n; i++) {
667 SContour *sc = &(l.elem[i]);
668 if(sc->l.n < 2) continue;
669 // The contours may not intersect, but they may share vertices; so
670 // testing a vertex for point-in-polygon may fail, but the midpoint
671 // of an edge is okay.
672 Vector pt = (((sc->l.elem[0]).p).Plus(sc->l.elem[1].p)).ScaledBy(0.5);
673
674 sc->timesEnclosed = 0;
675 bool outer = true;
676 for(j = 0; j < l.n; j++) {
677 if(i == j) continue;
678 SContour *sct = &(l.elem[j]);
679 if(sct->ContainsPointProjdToNormal(normal, pt)) {
680 outer = !outer;
681 (sc->timesEnclosed)++;
682 }
683 }
684
685 bool clockwise = sc->IsClockwiseProjdToNormal(normal);
686 if((clockwise && outer) || (!clockwise && !outer)) {
687 sc->Reverse();
688 sc->tag = 1;
689 }
690 }
691 }
692
IsEmpty(void)693 bool SPolygon::IsEmpty(void) {
694 if(l.n == 0 || l.elem[0].l.n == 0) return true;
695 return false;
696 }
697
AnyPoint(void)698 Vector SPolygon::AnyPoint(void) {
699 if(IsEmpty()) oops();
700 return l.elem[0].l.elem[0].p;
701 }
702
SelfIntersecting(Vector * intersectsAt)703 bool SPolygon::SelfIntersecting(Vector *intersectsAt) {
704 SEdgeList el = {};
705 MakeEdgesInto(&el);
706 SKdNodeEdges *kdtree = SKdNodeEdges::From(&el);
707
708 int cnt = 1;
709 el.l.ClearTags();
710
711 bool ret = false;
712 SEdge *se;
713 for(se = el.l.First(); se; se = el.l.NextAfter(se)) {
714 int inters = kdtree->AnyEdgeCrossings(se->a, se->b, cnt, intersectsAt);
715 if(inters != 1) {
716 ret = true;
717 break;
718 }
719 cnt++;
720 }
721
722 el.Clear();
723 return ret;
724 }
725
726 //-----------------------------------------------------------------------------
727 // Low-quality routines to cutter radius compensate a polygon. Assumes the
728 // polygon is in the xy plane, and the contours all go in the right direction
729 // with respect to normal (0, 0, -1).
730 //-----------------------------------------------------------------------------
OffsetInto(SPolygon * dest,double r)731 void SPolygon::OffsetInto(SPolygon *dest, double r) {
732 int i;
733 dest->Clear();
734 for(i = 0; i < l.n; i++) {
735 SContour *sc = &(l.elem[i]);
736 dest->AddEmptyContour();
737 sc->OffsetInto(&(dest->l.elem[dest->l.n-1]), r);
738 }
739 }
740 //-----------------------------------------------------------------------------
741 // Calculate the intersection point of two lines. Returns true for success,
742 // false if they're parallel.
743 //-----------------------------------------------------------------------------
IntersectionOfLines(double x0A,double y0A,double dxA,double dyA,double x0B,double y0B,double dxB,double dyB,double * xi,double * yi)744 static bool IntersectionOfLines(double x0A, double y0A, double dxA, double dyA,
745 double x0B, double y0B, double dxB, double dyB,
746 double *xi, double *yi)
747 {
748 double A[2][2];
749 double b[2];
750
751 // The line is given to us in the form:
752 // (x(t), y(t)) = (x0, y0) + t*(dx, dy)
753 // so first rewrite it as
754 // (x - x0, y - y0) dot (dy, -dx) = 0
755 // x*dy - x0*dy - y*dx + y0*dx = 0
756 // x*dy - y*dx = (x0*dy - y0*dx)
757
758 // So write the matrix, pre-pivoted.
759 if(fabs(dyA) > fabs(dyB)) {
760 A[0][0] = dyA; A[0][1] = -dxA; b[0] = x0A*dyA - y0A*dxA;
761 A[1][0] = dyB; A[1][1] = -dxB; b[1] = x0B*dyB - y0B*dxB;
762 } else {
763 A[1][0] = dyA; A[1][1] = -dxA; b[1] = x0A*dyA - y0A*dxA;
764 A[0][0] = dyB; A[0][1] = -dxB; b[0] = x0B*dyB - y0B*dxB;
765 }
766
767 // Check the determinant; if it's zero then no solution.
768 if(fabs(A[0][0]*A[1][1] - A[0][1]*A[1][0]) < LENGTH_EPS) {
769 return false;
770 }
771
772 // Solve
773 double v = A[1][0] / A[0][0];
774 A[1][0] -= A[0][0]*v;
775 A[1][1] -= A[0][1]*v;
776 b[1] -= b[0]*v;
777
778 // Back-substitute.
779 *yi = b[1] / A[1][1];
780 *xi = (b[0] - A[0][1]*(*yi)) / A[0][0];
781
782 return true;
783 }
OffsetInto(SContour * dest,double r)784 void SContour::OffsetInto(SContour *dest, double r) {
785 int i;
786
787 for(i = 0; i < l.n; i++) {
788 Vector a, b, c;
789 Vector dp, dn;
790 double thetan, thetap;
791
792 a = l.elem[WRAP(i-1, (l.n-1))].p;
793 b = l.elem[WRAP(i, (l.n-1))].p;
794 c = l.elem[WRAP(i+1, (l.n-1))].p;
795
796 dp = a.Minus(b);
797 thetap = atan2(dp.y, dp.x);
798
799 dn = b.Minus(c);
800 thetan = atan2(dn.y, dn.x);
801
802 // A short line segment in a badly-generated polygon might look
803 // okay but screw up our sense of direction.
804 if(dp.Magnitude() < LENGTH_EPS || dn.Magnitude() < LENGTH_EPS) {
805 continue;
806 }
807
808 if(thetan > thetap && (thetan - thetap) > PI) {
809 thetap += 2*PI;
810 }
811 if(thetan < thetap && (thetap - thetan) > PI) {
812 thetan += 2*PI;
813 }
814
815 if(fabs(thetan - thetap) < (1*PI)/180) {
816 Vector p = { b.x - r*sin(thetap), b.y + r*cos(thetap), 0 };
817 dest->AddPoint(p);
818 } else if(thetan < thetap) {
819 // This is an inside corner. We have two edges, Ep and En. Move
820 // out from their intersection by radius, normal to En, and
821 // then draw a line parallel to En. Move out from their
822 // intersection by radius, normal to Ep, and then draw a second
823 // line parallel to Ep. The point that we want to generate is
824 // the intersection of these two lines--it removes as much
825 // material as we can without removing any that we shouldn't.
826 double px0, py0, pdx, pdy;
827 double nx0, ny0, ndx, ndy;
828 double x = 0.0, y = 0.0;
829
830 px0 = b.x - r*sin(thetap);
831 py0 = b.y + r*cos(thetap);
832 pdx = cos(thetap);
833 pdy = sin(thetap);
834
835 nx0 = b.x - r*sin(thetan);
836 ny0 = b.y + r*cos(thetan);
837 ndx = cos(thetan);
838 ndy = sin(thetan);
839
840 IntersectionOfLines(px0, py0, pdx, pdy,
841 nx0, ny0, ndx, ndy,
842 &x, &y);
843
844 dest->AddPoint(Vector::From(x, y, 0));
845 } else {
846 if(fabs(thetap - thetan) < (6*PI)/180) {
847 Vector pp = { b.x - r*sin(thetap),
848 b.y + r*cos(thetap), 0 };
849 dest->AddPoint(pp);
850
851 Vector pn = { b.x - r*sin(thetan),
852 b.y + r*cos(thetan), 0 };
853 dest->AddPoint(pn);
854 } else {
855 double theta;
856 for(theta = thetap; theta <= thetan; theta += (6*PI)/180) {
857 Vector p = { b.x - r*sin(theta),
858 b.y + r*cos(theta), 0 };
859 dest->AddPoint(p);
860 }
861 }
862 }
863 }
864 }
865
866