1 //-----------------------------------------------------------------------------
2 // Top-level functions to compute the Boolean union or difference between
3 // two shells of rational polynomial surfaces.
4 //
5 // Copyright 2008-2013 Jonathan Westhues.
6 //-----------------------------------------------------------------------------
7 #include "solvespace.h"
8
9 static int I;
10
MakeFromUnionOf(SShell * a,SShell * b)11 void SShell::MakeFromUnionOf(SShell *a, SShell *b) {
12 MakeFromBoolean(a, b, AS_UNION);
13 }
14
MakeFromDifferenceOf(SShell * a,SShell * b)15 void SShell::MakeFromDifferenceOf(SShell *a, SShell *b) {
16 MakeFromBoolean(a, b, AS_DIFFERENCE);
17 }
18
19 //-----------------------------------------------------------------------------
20 // Take our original pwl curve. Wherever an edge intersects a surface within
21 // either agnstA or agnstB, split the piecewise linear element. Then refine
22 // the intersection so that it lies on all three relevant surfaces: the
23 // intersecting surface, srfA, and srfB. (So the pwl curve should lie at
24 // the intersection of srfA and srfB.) Return a new pwl curve with everything
25 // split.
26 //-----------------------------------------------------------------------------
27 static Vector LineStart, LineDirection;
ByTAlongLine(const void * av,const void * bv)28 static int ByTAlongLine(const void *av, const void *bv)
29 {
30 SInter *a = (SInter *)av,
31 *b = (SInter *)bv;
32
33 double ta = (a->p.Minus(LineStart)).DivPivoting(LineDirection),
34 tb = (b->p.Minus(LineStart)).DivPivoting(LineDirection);
35
36 return (ta > tb) ? 1 : -1;
37 }
MakeCopySplitAgainst(SShell * agnstA,SShell * agnstB,SSurface * srfA,SSurface * srfB)38 SCurve SCurve::MakeCopySplitAgainst(SShell *agnstA, SShell *agnstB,
39 SSurface *srfA, SSurface *srfB)
40 {
41 SCurve ret;
42 ret = *this;
43 ret.pts = {};
44
45 SCurvePt *p = pts.First();
46 if(!p) oops();
47 SCurvePt prev = *p;
48 ret.pts.Add(p);
49 p = pts.NextAfter(p);
50
51 for(; p; p = pts.NextAfter(p)) {
52 List<SInter> il = {};
53
54 // Find all the intersections with the two passed shells
55 if(agnstA)
56 agnstA->AllPointsIntersecting(prev.p, p->p, &il, true, false, true);
57 if(agnstB)
58 agnstB->AllPointsIntersecting(prev.p, p->p, &il, true, false, true);
59
60 if(il.n > 0) {
61 // The intersections were generated by intersecting the pwl
62 // edge against a surface; so they must be refined to lie
63 // exactly on the original curve.
64 il.ClearTags();
65 SInter *pi;
66 for(pi = il.First(); pi; pi = il.NextAfter(pi)) {
67 if(pi->srf == srfA || pi->srf == srfB) {
68 // The edge certainly intersects the surfaces that it
69 // trims (at its endpoints), but those ones don't count.
70 // They are culled later, but no sense calculating them
71 // and they will cause numerical problems (since two
72 // of the three surfaces they're refined to lie on will
73 // be identical, so the matrix will be singular).
74 pi->tag = 1;
75 continue;
76 }
77
78 Point2d puv;
79 (pi->srf)->ClosestPointTo(pi->p, &puv, false);
80
81 // Split the edge if the intersection lies within the surface's
82 // trim curves, or within the chord tol of the trim curve; want
83 // some slop if points are close to edge and pwl is too coarse,
84 // and it doesn't hurt to split unnecessarily.
85 Point2d dummy = { 0, 0 };
86 int c = (pi->srf->bsp) ? pi->srf->bsp->ClassifyPoint(puv, dummy, pi->srf) : SBspUv::OUTSIDE;
87 if(c == SBspUv::OUTSIDE) {
88 double d = VERY_POSITIVE;
89 if(pi->srf->bsp) d = pi->srf->bsp->MinimumDistanceToEdge(puv, pi->srf);
90 if(d > SS.ChordTolMm()) {
91 pi->tag = 1;
92 continue;
93 }
94 }
95
96 // We're keeping the intersection, so actually refine it.
97 (pi->srf)->PointOnSurfaces(srfA, srfB, &(puv.x), &(puv.y));
98 pi->p = (pi->srf)->PointAt(puv);
99 }
100 il.RemoveTagged();
101
102 // And now sort them in order along the line. Note that we must
103 // do that after refining, in case the refining would make two
104 // points switch places.
105 LineStart = prev.p;
106 LineDirection = (p->p).Minus(prev.p);
107 qsort(il.elem, il.n, sizeof(il.elem[0]), ByTAlongLine);
108
109 // And now uses the intersections to generate our split pwl edge(s)
110 Vector prev = Vector::From(VERY_POSITIVE, 0, 0);
111 for(pi = il.First(); pi; pi = il.NextAfter(pi)) {
112 // On-edge intersection will generate same split point for
113 // both surfaces, so don't create zero-length edge.
114 if(!prev.Equals(pi->p)) {
115 SCurvePt scpt;
116 scpt.tag = 0;
117 scpt.p = pi->p;
118 scpt.vertex = true;
119 ret.pts.Add(&scpt);
120 }
121 prev = pi->p;
122 }
123 }
124
125 il.Clear();
126 ret.pts.Add(p);
127 prev = *p;
128 }
129 return ret;
130 }
131
CopyCurvesSplitAgainst(bool opA,SShell * agnst,SShell * into)132 void SShell::CopyCurvesSplitAgainst(bool opA, SShell *agnst, SShell *into) {
133 SCurve *sc;
134 for(sc = curve.First(); sc; sc = curve.NextAfter(sc)) {
135 SCurve scn = sc->MakeCopySplitAgainst(agnst, NULL,
136 surface.FindById(sc->surfA),
137 surface.FindById(sc->surfB));
138 scn.source = opA ? SCurve::FROM_A : SCurve::FROM_B;
139
140 hSCurve hsc = into->curve.AddAndAssignId(&scn);
141 // And note the new ID so that we can rewrite the trims appropriately
142 sc->newH = hsc;
143 }
144 }
145
TrimFromEdgeList(SEdgeList * el,bool asUv)146 void SSurface::TrimFromEdgeList(SEdgeList *el, bool asUv) {
147 el->l.ClearTags();
148
149 STrimBy stb = {};
150 for(;;) {
151 // Find an edge, any edge; we'll start from there.
152 SEdge *se;
153 for(se = el->l.First(); se; se = el->l.NextAfter(se)) {
154 if(se->tag) continue;
155 break;
156 }
157 if(!se) break;
158 se->tag = 1;
159 stb.start = se->a;
160 stb.finish = se->b;
161 stb.curve.v = se->auxA;
162 stb.backwards = se->auxB ? true : false;
163
164 // Find adjoining edges from the same curve; those should be
165 // merged into a single trim.
166 bool merged;
167 do {
168 merged = false;
169 for(se = el->l.First(); se; se = el->l.NextAfter(se)) {
170 if(se->tag) continue;
171 if(se->auxA != (int)stb.curve.v) continue;
172 if(( se->auxB && !stb.backwards) ||
173 (!se->auxB && stb.backwards)) continue;
174
175 if((se->a).Equals(stb.finish)) {
176 stb.finish = se->b;
177 se->tag = 1;
178 merged = true;
179 } else if((se->b).Equals(stb.start)) {
180 stb.start = se->a;
181 se->tag = 1;
182 merged = true;
183 }
184 }
185 } while(merged);
186
187 if(asUv) {
188 stb.start = PointAt(stb.start.x, stb.start.y);
189 stb.finish = PointAt(stb.finish.x, stb.finish.y);
190 }
191
192 // And add the merged trim, with xyz (not uv like the polygon) pts
193 trim.Add(&stb);
194 }
195 }
196
KeepRegion(int type,bool opA,int shell,int orig)197 static bool KeepRegion(int type, bool opA, int shell, int orig)
198 {
199 bool inShell = (shell == SShell::INSIDE),
200 inSame = (shell == SShell::COINC_SAME),
201 inOpp = (shell == SShell::COINC_OPP),
202 inOrig = (orig == SShell::INSIDE);
203
204 bool inFace = inSame || inOpp;
205
206 // If these are correct, then they should be independent of inShell
207 // if inFace is true.
208 if(!inOrig) return false;
209 switch(type) {
210 case SShell::AS_UNION:
211 if(opA) {
212 return (!inShell && !inFace);
213 } else {
214 return (!inShell && !inFace) || inSame;
215 }
216
217 case SShell::AS_DIFFERENCE:
218 if(opA) {
219 return (!inShell && !inFace);
220 } else {
221 return (inShell && !inFace) || inSame;
222 }
223
224 default: oops();
225 }
226 }
KeepEdge(int type,bool opA,int indir_shell,int outdir_shell,int indir_orig,int outdir_orig)227 static bool KeepEdge(int type, bool opA,
228 int indir_shell, int outdir_shell,
229 int indir_orig, int outdir_orig)
230 {
231 bool keepIn = KeepRegion(type, opA, indir_shell, indir_orig),
232 keepOut = KeepRegion(type, opA, outdir_shell, outdir_orig);
233
234 // If the regions to the left and right of this edge are both in or both
235 // out, then this edge is not useful and should be discarded.
236 if(keepIn && !keepOut) return true;
237 return false;
238 }
239
TagByClassifiedEdge(int bspclass,int * indir,int * outdir)240 static void TagByClassifiedEdge(int bspclass, int *indir, int *outdir)
241 {
242 switch(bspclass) {
243 case SBspUv::INSIDE:
244 *indir = SShell::INSIDE;
245 *outdir = SShell::INSIDE;
246 break;
247
248 case SBspUv::OUTSIDE:
249 *indir = SShell::OUTSIDE;
250 *outdir = SShell::OUTSIDE;
251 break;
252
253 case SBspUv::EDGE_PARALLEL:
254 *indir = SShell::INSIDE;
255 *outdir = SShell::OUTSIDE;
256 break;
257
258 case SBspUv::EDGE_ANTIPARALLEL:
259 *indir = SShell::OUTSIDE;
260 *outdir = SShell::INSIDE;
261 break;
262
263 default:
264 dbp("TagByClassifiedEdge: fail!");
265 *indir = SShell::OUTSIDE;
266 *outdir = SShell::OUTSIDE;
267 break;
268 }
269 }
270
DEBUGEDGELIST(SEdgeList * sel,SSurface * surf)271 static void DEBUGEDGELIST(SEdgeList *sel, SSurface *surf) {
272 dbp("print %d edges", sel->l.n);
273 SEdge *se;
274 for(se = sel->l.First(); se; se = sel->l.NextAfter(se)) {
275 Vector mid = (se->a).Plus(se->b).ScaledBy(0.5);
276 Vector arrow = (se->b).Minus(se->a);
277 swap(arrow.x, arrow.y);
278 arrow.x *= -1;
279 arrow = arrow.WithMagnitude(0.01);
280 arrow = arrow.Plus(mid);
281
282 SS.nakedEdges.AddEdge(surf->PointAt(se->a.x, se->a.y),
283 surf->PointAt(se->b.x, se->b.y));
284 SS.nakedEdges.AddEdge(surf->PointAt(mid.x, mid.y),
285 surf->PointAt(arrow.x, arrow.y));
286 }
287 }
288
289 //-----------------------------------------------------------------------------
290 // We are given src, with at least one edge, and avoid, a list of points to
291 // avoid. We return a chain of edges (that share endpoints), such that no
292 // point within the avoid list ever occurs in the middle of a chain. And we
293 // delete the edges in that chain from our source list.
294 //-----------------------------------------------------------------------------
FindChainAvoiding(SEdgeList * src,SEdgeList * dest,SPointList * avoid)295 void SSurface::FindChainAvoiding(SEdgeList *src, SEdgeList *dest,
296 SPointList *avoid)
297 {
298 if(src->l.n < 1) oops();
299 // Start with an arbitrary edge.
300 dest->l.Add(&(src->l.elem[0]));
301 src->l.ClearTags();
302 src->l.elem[0].tag = 1;
303
304 bool added;
305 do {
306 added = false;
307 // The start and finish of the current edge chain
308 Vector s = dest->l.elem[0].a,
309 f = dest->l.elem[dest->l.n - 1].b;
310
311 // We can attach a new edge at the start or finish, as long as that
312 // start or finish point isn't in the list of points to avoid.
313 bool startOkay = !avoid->ContainsPoint(s),
314 finishOkay = !avoid->ContainsPoint(f);
315
316 // Now look for an unused edge that joins at the start or finish of
317 // our chain (if permitted by the avoid list).
318 SEdge *se;
319 for(se = src->l.First(); se; se = src->l.NextAfter(se)) {
320 if(se->tag) continue;
321 if(startOkay && s.Equals(se->b)) {
322 dest->l.AddToBeginning(se);
323 s = se->a;
324 se->tag = 1;
325 startOkay = !avoid->ContainsPoint(s);
326 } else if(finishOkay && f.Equals(se->a)) {
327 dest->l.Add(se);
328 f = se->b;
329 se->tag = 1;
330 finishOkay = !avoid->ContainsPoint(f);
331 } else {
332 continue;
333 }
334
335 added = true;
336 }
337 } while(added);
338
339 src->l.RemoveTagged();
340 }
341
EdgeNormalsWithinSurface(Point2d auv,Point2d buv,Vector * pt,Vector * enin,Vector * enout,Vector * surfn,uint32_t auxA,SShell * shell,SShell * sha,SShell * shb)342 void SSurface::EdgeNormalsWithinSurface(Point2d auv, Point2d buv,
343 Vector *pt,
344 Vector *enin, Vector *enout,
345 Vector *surfn,
346 uint32_t auxA,
347 SShell *shell, SShell *sha, SShell *shb)
348 {
349 // the midpoint of the edge
350 Point2d muv = (auv.Plus(buv)).ScaledBy(0.5);
351
352 *pt = PointAt(muv);
353
354 // If this edge just approximates a curve, then refine our midpoint so
355 // so that it actually lies on that curve too. Otherwise stuff like
356 // point-on-face tests will fail, since the point won't actually lie
357 // on the other face.
358 hSCurve hc = { auxA };
359 SCurve *sc = shell->curve.FindById(hc);
360 if(sc->isExact && sc->exact.deg != 1) {
361 double t;
362 sc->exact.ClosestPointTo(*pt, &t, false);
363 *pt = sc->exact.PointAt(t);
364 ClosestPointTo(*pt, &muv);
365 } else if(!sc->isExact) {
366 SSurface *trimmedA = sc->GetSurfaceA(sha, shb),
367 *trimmedB = sc->GetSurfaceB(sha, shb);
368 *pt = trimmedA->ClosestPointOnThisAndSurface(trimmedB, *pt);
369 ClosestPointTo(*pt, &muv);
370 }
371
372 *surfn = NormalAt(muv.x, muv.y);
373
374 // Compute the edge's inner normal in xyz space.
375 Vector ab = (PointAt(auv)).Minus(PointAt(buv)),
376 enxyz = (ab.Cross(*surfn)).WithMagnitude(SS.ChordTolMm());
377 // And based on that, compute the edge's inner normal in uv space. This
378 // vector is perpendicular to the edge in xyz, but not necessarily in uv.
379 Vector tu, tv;
380 TangentsAt(muv.x, muv.y, &tu, &tv);
381 Point2d enuv;
382 enuv.x = enxyz.Dot(tu) / tu.MagSquared();
383 enuv.y = enxyz.Dot(tv) / tv.MagSquared();
384
385 // Compute the inner and outer normals of this edge (within the srf),
386 // in xyz space. These are not necessarily antiparallel, if the
387 // surface is curved.
388 Vector pin = PointAt(muv.Minus(enuv)),
389 pout = PointAt(muv.Plus(enuv));
390 *enin = pin.Minus(*pt),
391 *enout = pout.Minus(*pt);
392 }
393
394 //-----------------------------------------------------------------------------
395 // Trim this surface against the specified shell, in the way that's appropriate
396 // for the specified Boolean operation type (and which operand we are). We
397 // also need a pointer to the shell that contains our own surface, since that
398 // contains our original trim curves.
399 //-----------------------------------------------------------------------------
MakeCopyTrimAgainst(SShell * parent,SShell * sha,SShell * shb,SShell * into,int type)400 SSurface SSurface::MakeCopyTrimAgainst(SShell *parent,
401 SShell *sha, SShell *shb,
402 SShell *into,
403 int type)
404 {
405 bool opA = (parent == sha);
406 SShell *agnst = opA ? shb : sha;
407
408 SSurface ret;
409 // The returned surface is identical, just the trim curves change
410 ret = *this;
411 ret.trim = {};
412
413 // First, build a list of the existing trim curves; update them to use
414 // the split curves.
415 STrimBy *stb;
416 for(stb = trim.First(); stb; stb = trim.NextAfter(stb)) {
417 STrimBy stn = *stb;
418 stn.curve = (parent->curve.FindById(stn.curve))->newH;
419 ret.trim.Add(&stn);
420 }
421
422 if(type == SShell::AS_DIFFERENCE && !opA) {
423 // The second operand of a Boolean difference gets turned inside out
424 ret.Reverse();
425 }
426
427 // Build up our original trim polygon; remember the coordinates could
428 // be changed if we just flipped the surface normal, and we are using
429 // the split curves (not the original curves).
430 SEdgeList orig = {};
431 ret.MakeEdgesInto(into, &orig, AS_UV);
432 ret.trim.Clear();
433 // which means that we can't necessarily use the old BSP...
434 SBspUv *origBsp = SBspUv::From(&orig, &ret);
435
436 // And now intersect the other shell against us
437 SEdgeList inter = {};
438
439 SSurface *ss;
440 for(ss = agnst->surface.First(); ss; ss = agnst->surface.NextAfter(ss)) {
441 SCurve *sc;
442 for(sc = into->curve.First(); sc; sc = into->curve.NextAfter(sc)) {
443 if(sc->source != SCurve::FROM_INTERSECTION) continue;
444 if(opA) {
445 if(sc->surfA.v != h.v || sc->surfB.v != ss->h.v) continue;
446 } else {
447 if(sc->surfB.v != h.v || sc->surfA.v != ss->h.v) continue;
448 }
449
450 int i;
451 for(i = 1; i < sc->pts.n; i++) {
452 Vector a = sc->pts.elem[i-1].p,
453 b = sc->pts.elem[i].p;
454
455 Point2d auv, buv;
456 ss->ClosestPointTo(a, &(auv.x), &(auv.y));
457 ss->ClosestPointTo(b, &(buv.x), &(buv.y));
458
459 int c = (ss->bsp) ? ss->bsp->ClassifyEdge(auv, buv, ss) : SBspUv::OUTSIDE;
460 if(c != SBspUv::OUTSIDE) {
461 Vector ta = Vector::From(0, 0, 0);
462 Vector tb = Vector::From(0, 0, 0);
463 ret.ClosestPointTo(a, &(ta.x), &(ta.y));
464 ret.ClosestPointTo(b, &(tb.x), &(tb.y));
465
466 Vector tn = ret.NormalAt(ta.x, ta.y);
467 Vector sn = ss->NormalAt(auv.x, auv.y);
468
469 // We are subtracting the portion of our surface that
470 // lies in the shell, so the in-plane edge normal should
471 // point opposite to the surface normal.
472 bool bkwds = true;
473 if((tn.Cross(b.Minus(a))).Dot(sn) < 0) bkwds = !bkwds;
474 if(type == SShell::AS_DIFFERENCE && !opA) bkwds = !bkwds;
475 if(bkwds) {
476 inter.AddEdge(tb, ta, sc->h.v, 1);
477 } else {
478 inter.AddEdge(ta, tb, sc->h.v, 0);
479 }
480 }
481 }
482 }
483 }
484
485 // Record all the points where more than two edges join, which I will call
486 // the choosing points. If two edges join at a non-choosing point, then
487 // they must either both be kept or both be discarded (since that would
488 // otherwise create an open contour).
489 SPointList choosing = {};
490 SEdge *se;
491 for(se = orig.l.First(); se; se = orig.l.NextAfter(se)) {
492 choosing.IncrementTagFor(se->a);
493 choosing.IncrementTagFor(se->b);
494 }
495 for(se = inter.l.First(); se; se = inter.l.NextAfter(se)) {
496 choosing.IncrementTagFor(se->a);
497 choosing.IncrementTagFor(se->b);
498 }
499 SPoint *sp;
500 for(sp = choosing.l.First(); sp; sp = choosing.l.NextAfter(sp)) {
501 if(sp->tag == 2) {
502 sp->tag = 1;
503 } else {
504 sp->tag = 0;
505 }
506 }
507 choosing.l.RemoveTagged();
508
509 // The list of edges to trim our new surface, a combination of edges from
510 // our original and intersecting edge lists.
511 SEdgeList final = {};
512
513 while(orig.l.n > 0) {
514 SEdgeList chain = {};
515 FindChainAvoiding(&orig, &chain, &choosing);
516
517 // Arbitrarily choose an edge within the chain to classify; they
518 // should all be the same, though.
519 se = &(chain.l.elem[chain.l.n/2]);
520
521 Point2d auv = (se->a).ProjectXy(),
522 buv = (se->b).ProjectXy();
523
524 Vector pt, enin, enout, surfn;
525 ret.EdgeNormalsWithinSurface(auv, buv, &pt, &enin, &enout, &surfn,
526 se->auxA, into, sha, shb);
527
528 int indir_shell, outdir_shell, indir_orig, outdir_orig;
529
530 indir_orig = SShell::INSIDE;
531 outdir_orig = SShell::OUTSIDE;
532
533 agnst->ClassifyEdge(&indir_shell, &outdir_shell,
534 ret.PointAt(auv), ret.PointAt(buv), pt,
535 enin, enout, surfn);
536
537 if(KeepEdge(type, opA, indir_shell, outdir_shell,
538 indir_orig, outdir_orig))
539 {
540 for(se = chain.l.First(); se; se = chain.l.NextAfter(se)) {
541 final.AddEdge(se->a, se->b, se->auxA, se->auxB);
542 }
543 }
544 chain.Clear();
545 }
546
547 while(inter.l.n > 0) {
548 SEdgeList chain = {};
549 FindChainAvoiding(&inter, &chain, &choosing);
550
551 // Any edge in the chain, same as above.
552 se = &(chain.l.elem[chain.l.n/2]);
553
554 Point2d auv = (se->a).ProjectXy(),
555 buv = (se->b).ProjectXy();
556
557 Vector pt, enin, enout, surfn;
558 ret.EdgeNormalsWithinSurface(auv, buv, &pt, &enin, &enout, &surfn,
559 se->auxA, into, sha, shb);
560
561 int indir_shell, outdir_shell, indir_orig, outdir_orig;
562
563 int c_this = (origBsp) ? origBsp->ClassifyEdge(auv, buv, &ret) : SBspUv::OUTSIDE;
564 TagByClassifiedEdge(c_this, &indir_orig, &outdir_orig);
565
566 agnst->ClassifyEdge(&indir_shell, &outdir_shell,
567 ret.PointAt(auv), ret.PointAt(buv), pt,
568 enin, enout, surfn);
569
570 if(KeepEdge(type, opA, indir_shell, outdir_shell,
571 indir_orig, outdir_orig))
572 {
573 for(se = chain.l.First(); se; se = chain.l.NextAfter(se)) {
574 final.AddEdge(se->a, se->b, se->auxA, se->auxB);
575 }
576 }
577 chain.Clear();
578 }
579
580 // Cull extraneous edges; duplicates or anti-parallel pairs. In particular,
581 // we can get duplicate edges if our surface intersects the other shell
582 // at an edge, so that both surfaces intersect coincident (and both
583 // generate an intersection edge).
584 final.CullExtraneousEdges();
585
586 // Use our reassembled edges to trim the new surface.
587 ret.TrimFromEdgeList(&final, true);
588
589 SPolygon poly = {};
590 final.l.ClearTags();
591 if(!final.AssemblePolygon(&poly, NULL, true)) {
592 into->booleanFailed = true;
593 dbp("failed: I=%d, avoid=%d", I, choosing.l.n);
594 DEBUGEDGELIST(&final, &ret);
595 }
596 poly.Clear();
597
598 choosing.Clear();
599 final.Clear();
600 inter.Clear();
601 orig.Clear();
602 return ret;
603 }
604
CopySurfacesTrimAgainst(SShell * sha,SShell * shb,SShell * into,int type)605 void SShell::CopySurfacesTrimAgainst(SShell *sha, SShell *shb, SShell *into,
606 int type)
607 {
608 SSurface *ss;
609 for(ss = surface.First(); ss; ss = surface.NextAfter(ss)) {
610 SSurface ssn;
611 ssn = ss->MakeCopyTrimAgainst(this, sha, shb, into, type);
612 ss->newH = into->surface.AddAndAssignId(&ssn);
613 I++;
614 }
615 }
616
MakeIntersectionCurvesAgainst(SShell * agnst,SShell * into)617 void SShell::MakeIntersectionCurvesAgainst(SShell *agnst, SShell *into) {
618 SSurface *sa;
619 for(sa = surface.First(); sa; sa = surface.NextAfter(sa)) {
620 SSurface *sb;
621 for(sb = agnst->surface.First(); sb; sb = agnst->surface.NextAfter(sb)){
622 // Intersect every surface from our shell against every surface
623 // from agnst; this will add zero or more curves to the curve
624 // list for into.
625 sa->IntersectAgainst(sb, this, agnst, into);
626 }
627 }
628 }
629
CleanupAfterBoolean(void)630 void SShell::CleanupAfterBoolean(void) {
631 SSurface *ss;
632 for(ss = surface.First(); ss; ss = surface.NextAfter(ss)) {
633 ss->edges.Clear();
634 }
635 }
636
637 //-----------------------------------------------------------------------------
638 // All curves contain handles to the two surfaces that they trim. After a
639 // Boolean or assembly, we must rewrite those handles to refer to the curves
640 // by their new IDs.
641 //-----------------------------------------------------------------------------
RewriteSurfaceHandlesForCurves(SShell * a,SShell * b)642 void SShell::RewriteSurfaceHandlesForCurves(SShell *a, SShell *b) {
643 SCurve *sc;
644 for(sc = curve.First(); sc; sc = curve.NextAfter(sc)) {
645 sc->surfA = sc->GetSurfaceA(a, b)->newH,
646 sc->surfB = sc->GetSurfaceB(a, b)->newH;
647 }
648 }
649
650 //-----------------------------------------------------------------------------
651 // Copy all the surfaces and curves from two different shells into a single
652 // shell. The only difficulty is to rewrite all of their handles; we don't
653 // look for any surface intersections, so if two objects interfere then the
654 // result is just self-intersecting. This is used for assembly, since it's
655 // much faster than merging as union.
656 //-----------------------------------------------------------------------------
MakeFromAssemblyOf(SShell * a,SShell * b)657 void SShell::MakeFromAssemblyOf(SShell *a, SShell *b) {
658 booleanFailed = false;
659
660 Vector t = Vector::From(0, 0, 0);
661 Quaternion q = Quaternion::IDENTITY;
662 int i = 0;
663 SShell *ab;
664
665 // First, copy over all the curves. Note which shell (a or b) each curve
666 // came from, but assign it a new ID.
667 SCurve *c, cn;
668 for(i = 0; i < 2; i++) {
669 ab = (i == 0) ? a : b;
670 for(c = ab->curve.First(); c; c = ab->curve.NextAfter(c)) {
671 cn = SCurve::FromTransformationOf(c, t, q, 1.0);
672 cn.source = (i == 0) ? SCurve::FROM_A : SCurve::FROM_B;
673 // surfA and surfB are wrong now, and we can't fix them until
674 // we've assigned IDs to the surfaces. So we'll get that later.
675 c->newH = curve.AddAndAssignId(&cn);
676 }
677 }
678
679 // Likewise copy over all the surfaces.
680 SSurface *s, sn;
681 for(i = 0; i < 2; i++) {
682 ab = (i == 0) ? a : b;
683 for(s = ab->surface.First(); s; s = ab->surface.NextAfter(s)) {
684 sn = SSurface::FromTransformationOf(s, t, q, 1.0, true);
685 // All the trim curve IDs get rewritten; we know the new handles
686 // to the curves since we recorded them in the previous step.
687 STrimBy *stb;
688 for(stb = sn.trim.First(); stb; stb = sn.trim.NextAfter(stb)) {
689 stb->curve = ab->curve.FindById(stb->curve)->newH;
690 }
691 s->newH = surface.AddAndAssignId(&sn);
692 }
693 }
694
695 // Finally, rewrite the surfaces associated with each curve to use the
696 // new handles.
697 RewriteSurfaceHandlesForCurves(a, b);
698 }
699
MakeFromBoolean(SShell * a,SShell * b,int type)700 void SShell::MakeFromBoolean(SShell *a, SShell *b, int type) {
701 booleanFailed = false;
702
703 a->MakeClassifyingBsps(NULL);
704 b->MakeClassifyingBsps(NULL);
705
706 // Copy over all the original curves, splitting them so that a
707 // piecwise linear segment never crosses a surface from the other
708 // shell.
709 a->CopyCurvesSplitAgainst(true, b, this);
710 b->CopyCurvesSplitAgainst(false, a, this);
711
712 // Generate the intersection curves for each surface in A against all
713 // the surfaces in B (which is all of the intersection curves).
714 a->MakeIntersectionCurvesAgainst(b, this);
715
716 SCurve *sc;
717 for(sc = curve.First(); sc; sc = curve.NextAfter(sc)) {
718 SSurface *srfA = sc->GetSurfaceA(a, b),
719 *srfB = sc->GetSurfaceB(a, b);
720
721 sc->RemoveShortSegments(srfA, srfB);
722 }
723
724 // And clean up the piecewise linear things we made as a calculation aid
725 a->CleanupAfterBoolean();
726 b->CleanupAfterBoolean();
727 // Remake the classifying BSPs with the split (and short-segment-removed)
728 // curves
729 a->MakeClassifyingBsps(this);
730 b->MakeClassifyingBsps(this);
731
732 if(b->surface.n == 0 || a->surface.n == 0) {
733 I = 1000000;
734 } else {
735 I = 0;
736 }
737 // Then trim and copy the surfaces
738 a->CopySurfacesTrimAgainst(a, b, this, type);
739 b->CopySurfacesTrimAgainst(a, b, this, type);
740
741 // Now that we've copied the surfaces, we know their new hSurfaces, so
742 // rewrite the curves to refer to the surfaces by their handles in the
743 // result.
744 RewriteSurfaceHandlesForCurves(a, b);
745
746 // And clean up the piecewise linear things we made as a calculation aid
747 a->CleanupAfterBoolean();
748 b->CleanupAfterBoolean();
749 }
750
751 //-----------------------------------------------------------------------------
752 // All of the BSP routines that we use to perform and accelerate polygon ops.
753 //-----------------------------------------------------------------------------
MakeClassifyingBsps(SShell * useCurvesFrom)754 void SShell::MakeClassifyingBsps(SShell *useCurvesFrom) {
755 SSurface *ss;
756 for(ss = surface.First(); ss; ss = surface.NextAfter(ss)) {
757 ss->MakeClassifyingBsp(this, useCurvesFrom);
758 }
759 }
760
MakeClassifyingBsp(SShell * shell,SShell * useCurvesFrom)761 void SSurface::MakeClassifyingBsp(SShell *shell, SShell *useCurvesFrom) {
762 SEdgeList el = {};
763
764 MakeEdgesInto(shell, &el, AS_UV, useCurvesFrom);
765 bsp = SBspUv::From(&el, this);
766 el.Clear();
767
768 edges = {};
769 MakeEdgesInto(shell, &edges, AS_XYZ, useCurvesFrom);
770 }
771
Alloc(void)772 SBspUv *SBspUv::Alloc(void) {
773 return (SBspUv *)AllocTemporary(sizeof(SBspUv));
774 }
775
ByLength(const void * av,const void * bv)776 static int ByLength(const void *av, const void *bv)
777 {
778 SEdge *a = (SEdge *)av,
779 *b = (SEdge *)bv;
780
781 double la = (a->a).Minus(a->b).Magnitude(),
782 lb = (b->a).Minus(b->b).Magnitude();
783
784 // Sort in descending order, longest first. This improves numerical
785 // stability for the normals.
786 return (la < lb) ? 1 : -1;
787 }
From(SEdgeList * el,SSurface * srf)788 SBspUv *SBspUv::From(SEdgeList *el, SSurface *srf) {
789 SEdgeList work = {};
790
791 SEdge *se;
792 for(se = el->l.First(); se; se = el->l.NextAfter(se)) {
793 work.AddEdge(se->a, se->b, se->auxA, se->auxB);
794 }
795 qsort(work.l.elem, work.l.n, sizeof(work.l.elem[0]), ByLength);
796
797 SBspUv *bsp = NULL;
798 for(se = work.l.First(); se; se = work.l.NextAfter(se)) {
799 bsp = InsertOrCreateEdge(bsp, (se->a).ProjectXy(), (se->b).ProjectXy(), srf);
800 }
801
802 work.Clear();
803 return bsp;
804 }
805
806 //-----------------------------------------------------------------------------
807 // The points in this BSP are in uv space, but we want to apply our tolerances
808 // consistently in xyz (i.e., we want to say a point is on-edge if its xyz
809 // distance to that edge is less than LENGTH_EPS, irrespective of its distance
810 // in uv). So we linearize the surface about the point we're considering and
811 // then do the test. That preserves point-on-line relationships, and the only
812 // time we care about exact correctness is when we're very close to the line,
813 // which is when the linearization is accurate.
814 //-----------------------------------------------------------------------------
ScalePoints(Point2d * pt,Point2d * a,Point2d * b,SSurface * srf)815 void SBspUv::ScalePoints(Point2d *pt, Point2d *a, Point2d *b, SSurface *srf) {
816 Vector tu, tv;
817 srf->TangentsAt(pt->x, pt->y, &tu, &tv);
818 double mu = tu.Magnitude(), mv = tv.Magnitude();
819
820 pt->x *= mu; pt->y *= mv;
821 a ->x *= mu; a ->y *= mv;
822 b ->x *= mu; b ->y *= mv;
823 }
ScaledSignedDistanceToLine(Point2d pt,Point2d a,Point2d b,SSurface * srf)824 double SBspUv::ScaledSignedDistanceToLine(Point2d pt, Point2d a, Point2d b,
825 SSurface *srf)
826 {
827 ScalePoints(&pt, &a, &b, srf);
828
829 Point2d n = ((b.Minus(a)).Normal()).WithMagnitude(1);
830 double d = a.Dot(n);
831
832 return pt.Dot(n) - d;
833 }
ScaledDistanceToLine(Point2d pt,Point2d a,Point2d b,bool seg,SSurface * srf)834 double SBspUv::ScaledDistanceToLine(Point2d pt, Point2d a, Point2d b, bool seg,
835 SSurface *srf)
836 {
837 ScalePoints(&pt, &a, &b, srf);
838
839 return pt.DistanceToLine(a, b, seg);
840 }
841
InsertOrCreateEdge(SBspUv * where,const Point2d & ea,const Point2d & eb,SSurface * srf)842 SBspUv *SBspUv::InsertOrCreateEdge(SBspUv *where, const Point2d &ea, const Point2d &eb, SSurface *srf) {
843 if(where == NULL) {
844 SBspUv *ret = Alloc();
845 ret->a = ea;
846 ret->b = eb;
847 return ret;
848 }
849 where->InsertEdge(ea, eb, srf);
850 return where;
851 }
852
InsertEdge(Point2d ea,Point2d eb,SSurface * srf)853 void SBspUv::InsertEdge(Point2d ea, Point2d eb, SSurface *srf) {
854 double dea = ScaledSignedDistanceToLine(ea, a, b, srf),
855 deb = ScaledSignedDistanceToLine(eb, a, b, srf);
856
857 if(fabs(dea) < LENGTH_EPS && fabs(deb) < LENGTH_EPS) {
858 // Line segment is coincident with this one, store in same node
859 SBspUv *m = Alloc();
860 m->a = ea;
861 m->b = eb;
862 m->more = more;
863 more = m;
864 } else if(fabs(dea) < LENGTH_EPS) {
865 // Point A lies on this lie, but point B does not
866 if(deb > 0) {
867 pos = InsertOrCreateEdge(pos, ea, eb, srf);
868 } else {
869 neg = InsertOrCreateEdge(neg, ea, eb, srf);
870 }
871 } else if(fabs(deb) < LENGTH_EPS) {
872 // Point B lies on this lie, but point A does not
873 if(dea > 0) {
874 pos = InsertOrCreateEdge(pos, ea, eb, srf);
875 } else {
876 neg = InsertOrCreateEdge(neg, ea, eb, srf);
877 }
878 } else if(dea > 0 && deb > 0) {
879 pos = InsertOrCreateEdge(pos, ea, eb, srf);
880 } else if(dea < 0 && deb < 0) {
881 neg = InsertOrCreateEdge(neg, ea, eb, srf);
882 } else {
883 // New edge crosses this one; we need to split.
884 Point2d n = ((b.Minus(a)).Normal()).WithMagnitude(1);
885 double d = a.Dot(n);
886 double t = (d - n.Dot(ea)) / (n.Dot(eb.Minus(ea)));
887 Point2d pi = ea.Plus((eb.Minus(ea)).ScaledBy(t));
888 if(dea > 0) {
889 pos = InsertOrCreateEdge(pos, ea, pi, srf);
890 neg = InsertOrCreateEdge(neg, pi, eb, srf);
891 } else {
892 neg = InsertOrCreateEdge(neg, ea, pi, srf);
893 pos = InsertOrCreateEdge(pos, pi, eb, srf);
894 }
895 }
896 return;
897 }
898
ClassifyPoint(Point2d p,Point2d eb,SSurface * srf)899 int SBspUv::ClassifyPoint(Point2d p, Point2d eb, SSurface *srf) {
900
901 double dp = ScaledSignedDistanceToLine(p, a, b, srf);
902
903 if(fabs(dp) < LENGTH_EPS) {
904 SBspUv *f = this;
905 while(f) {
906 Point2d ba = (f->b).Minus(f->a);
907 if(ScaledDistanceToLine(p, f->a, ba, true, srf) < LENGTH_EPS) {
908 if(ScaledDistanceToLine(eb, f->a, ba, false, srf) < LENGTH_EPS){
909 if(ba.Dot(eb.Minus(p)) > 0) {
910 return EDGE_PARALLEL;
911 } else {
912 return EDGE_ANTIPARALLEL;
913 }
914 } else {
915 return EDGE_OTHER;
916 }
917 }
918 f = f->more;
919 }
920 // Pick arbitrarily which side to send it down, doesn't matter
921 int c1 = neg ? neg->ClassifyPoint(p, eb, srf) : OUTSIDE;
922 int c2 = pos ? pos->ClassifyPoint(p, eb, srf) : INSIDE;
923 if(c1 != c2) {
924 dbp("MISMATCH: %d %d %08x %08x", c1, c2, neg, pos);
925 }
926 return c1;
927 } else if(dp > 0) {
928 return pos ? pos->ClassifyPoint(p, eb, srf) : INSIDE;
929 } else {
930 return neg ? neg->ClassifyPoint(p, eb, srf) : OUTSIDE;
931 }
932 }
933
ClassifyEdge(Point2d ea,Point2d eb,SSurface * srf)934 int SBspUv::ClassifyEdge(Point2d ea, Point2d eb, SSurface *srf) {
935 int ret = ClassifyPoint((ea.Plus(eb)).ScaledBy(0.5), eb, srf);
936 if(ret == EDGE_OTHER) {
937 // Perhaps the edge is tangent at its midpoint (and we screwed up
938 // somewhere earlier and failed to split it); try a different
939 // point on the edge.
940 ret = ClassifyPoint(ea.Plus((eb.Minus(ea)).ScaledBy(0.294)), eb, srf);
941 }
942 return ret;
943 }
944
MinimumDistanceToEdge(Point2d p,SSurface * srf)945 double SBspUv::MinimumDistanceToEdge(Point2d p, SSurface *srf) {
946
947 double dn = (neg) ? neg->MinimumDistanceToEdge(p, srf) : VERY_POSITIVE;
948 double dp = (pos) ? pos->MinimumDistanceToEdge(p, srf) : VERY_POSITIVE;
949
950 Point2d as = a, bs = b;
951 ScalePoints(&p, &as, &bs, srf);
952 double d = p.DistanceToLine(as, bs.Minus(as), true);
953
954 return min(d, min(dn, dp));
955 }
956
957