From 90842131ff167e71252ac3abc694884b34331584 Mon Sep 17 00:00:00 2001 From: Jonathan Westhues Date: Mon, 9 Feb 2009 04:40:48 -0800 Subject: [PATCH] Make Boolean union work when the shells have coincident plane faces. Still on planes only, no curved surface intersections. [git-p4: depot-paths = "//depot/solvespace/": change = 1912] --- solvespace.h | 5 +- srf/boolean.cpp | 113 ++++++++++++++++++++++++++------- srf/ratpoly.cpp | 32 ++++++---- srf/surface.h | 12 +++- srf/surfinter.cpp | 155 +++++++++++++++++++++++++++++++++++++++------- 5 files changed, 255 insertions(+), 62 deletions(-) diff --git a/solvespace.h b/solvespace.h index acab5e6..2414e14 100644 --- a/solvespace.h +++ b/solvespace.h @@ -20,7 +20,6 @@ inline int WRAP(int v, int n) { while(v < 0) v += n; return v; } - inline double WRAP_NOT_0(double v, double n) { // Clamp it to the range (0, n] while(v > n) v -= n; @@ -50,6 +49,10 @@ typedef signed short SWORD; #include #include +inline double Random(double vmax) { + return (vmax*rand()) / RAND_MAX; +} + class Expr; class ExprVector; class ExprQuaternion; diff --git a/srf/boolean.cpp b/srf/boolean.cpp index 8975f5c..4e4cb79 100644 --- a/srf/boolean.cpp +++ b/srf/boolean.cpp @@ -47,10 +47,16 @@ SCurve SCurve::MakeCopySplitAgainst(SShell *agnstA, SShell *agnstB) { LineStart = prev; LineDirection = p->Minus(prev); qsort(il.elem, il.n, sizeof(il.elem[0]), ByTAlongLine); - + + Vector prev = Vector::From(VERY_POSITIVE, 0, 0); SInter *pi; for(pi = il.First(); pi; pi = il.NextAfter(pi)) { - ret.pts.Add(&(pi->p)); + // On-edge intersection will generate same split point for + // both surfaces, so don't create zero-length edge. + if(!prev.Equals(pi->p)) { + ret.pts.Add(&(pi->p)); + } + prev = pi->p; } } @@ -150,6 +156,17 @@ SSurface SSurface::MakeCopyTrimAgainst(SShell *agnst, SShell *parent, ret.MakeEdgesInto(into, &orig, true); ret.trim.Clear(); + // Find any surfaces from the other shell that are coincident with ours, + // and get the intersection polygons, grouping surfaces with the same + // and with opposite normal separately. + SEdgeList sameNormal, oppositeNormal; + ZERO(&sameNormal); + ZERO(&oppositeNormal); + agnst->MakeCoincidentEdgesInto(this, true, &sameNormal); + agnst->MakeCoincidentEdgesInto(this, false, &oppositeNormal); + // and build the trees for quick in-polygon testing + SBspUv *sameBsp = SBspUv::From(&sameNormal); + SBspUv *oppositeBsp = SBspUv::From(&oppositeNormal); // And now intersect the other shell against us SEdgeList inter; @@ -176,7 +193,7 @@ SSurface SSurface::MakeCopyTrimAgainst(SShell *agnst, SShell *parent, ss->ClosestPointTo(b, &(buv.x), &(buv.y)); int c = ss->bsp->ClassifyEdge(auv, buv); - if(c == SBspUv::INSIDE) { + if(c != SBspUv::OUTSIDE) { Vector ta = Vector::From(0, 0, 0); Vector tb = Vector::From(0, 0, 0); ClosestPointTo(a, &(ta.x), &(ta.y)); @@ -198,23 +215,78 @@ SSurface SSurface::MakeCopyTrimAgainst(SShell *agnst, SShell *parent, SEdgeList final; ZERO(&final); - SBspUv *interbsp = SBspUv::From(&inter); - SEdge *se; for(se = orig.l.First(); se; se = orig.l.NextAfter(se)) { - Vector ea = PointAt(se->a.x, se->a.y), - eb = PointAt(se->b.x, se->b.y); - int c = agnst->ClassifyPoint(ea.Plus(eb).ScaledBy(0.5)); + Point2d auv = (se->a).ProjectXy(), + buv = (se->b).ProjectXy(); - if(c == SShell::OUTSIDE) { + int c_same = sameBsp->ClassifyEdge(auv, buv); + int c_opp = oppositeBsp->ClassifyEdge(auv, buv); + + // Get the midpoint of this edge + Point2d am = (auv.Plus(buv)).ScaledBy(0.5); + Vector pt = PointAt(am.x, am.y); + // and the outer normal from the trim polygon (within the surface) + Vector n = NormalAt(am.x, am.y); + Vector ea = PointAt(auv.x, auv.y), + eb = PointAt(buv.x, buv.y); + Vector ptout = n.Cross((eb.Minus(ea))); + + int c_shell = agnst->ClassifyPoint(pt, ptout); + + bool keep; + if(c_opp != SBspUv::OUTSIDE) { + // Edge lies on coincident (opposite normals) surface of agnst + keep = (c_opp == SBspUv::OUTSIDE ) || + (c_opp == SBspUv::EDGE_ANTIPARALLEL ); + + } else if(c_same != SBspUv::OUTSIDE) { + // Edge lies on coincident (same normals) surface of agnst + if(opA) { + keep = true; + } else { + keep = false; + } + + } else { + // Edge does not lie on a coincident surface + keep = (c_shell == SShell::OUTSIDE ) || + (c_shell == SShell::ON_ANTIPARALLEL ); + } + + if(keep) { final.AddEdge(se->a, se->b, se->auxA, se->auxB); } } for(se = inter.l.First(); se; se = inter.l.NextAfter(se)) { - int c = bsp->ClassifyEdge(se->a.ProjectXy(), se->b.ProjectXy()); + Point2d auv = (se->a).ProjectXy(), + buv = (se->b).ProjectXy(); - if(I == 4) { + int c_this = bsp->ClassifyEdge(auv, buv); + int c_same = sameBsp->ClassifyEdge(auv, buv); + int c_opp = oppositeBsp->ClassifyEdge(auv, buv); + + bool keep; + if(c_opp != SBspUv::OUTSIDE) { + keep = (c_this == SBspUv::INSIDE); + } else if(c_same != SBspUv::OUTSIDE) { + if(opA) { + keep = false; + } else { + keep = (c_this == SBspUv::INSIDE); + } + } else { + keep = (c_this == SBspUv::INSIDE); + } + + if(keep) { + final.AddEdge(se->b, se->a, se->auxA, !se->auxB); + } + } + + for(se = final.l.First(); se; se = final.l.NextAfter(se)) { + if(I == 0) { Vector mid = (se->a).Plus(se->b).ScaledBy(0.5); Vector arrow = (se->b).Minus(se->a); SWAP(double, arrow.x, arrow.y); @@ -222,24 +294,17 @@ SSurface SSurface::MakeCopyTrimAgainst(SShell *agnst, SShell *parent, arrow = arrow.WithMagnitude(0.03); arrow = arrow.Plus(mid); - if(c == SBspUv::INSIDE) { - SS.nakedEdges.AddEdge(PointAt(se->a.x, se->a.y), - PointAt(se->b.x, se->b.y)); - SS.nakedEdges.AddEdge(PointAt(mid.x, mid.y), - PointAt(arrow.x, arrow.y)); - } + SS.nakedEdges.AddEdge(PointAt(se->a.x, se->a.y), + PointAt(se->b.x, se->b.y)); + SS.nakedEdges.AddEdge(PointAt(mid.x, mid.y), + PointAt(arrow.x, arrow.y)); } - - if(c == SBspUv::INSIDE) { - final.AddEdge(se->b, se->a, se->auxA, !se->auxB); - } - } - - for(se = final.l.First(); se; se = final.l.NextAfter(se)) { } ret.TrimFromEdgeList(&final); + sameNormal.Clear(); + oppositeNormal.Clear(); final.Clear(); inter.Clear(); orig.Clear(); diff --git a/srf/ratpoly.cpp b/srf/ratpoly.cpp index b5c9447..bdc6dfd 100644 --- a/srf/ratpoly.cpp +++ b/srf/ratpoly.cpp @@ -498,18 +498,23 @@ Vector SSurface::NormalAt(double u, double v) { void SSurface::ClosestPointTo(Vector p, double *u, double *v) { int i, j; - double minDist = 1e10; - double res = 7.0; - for(i = 0; i < (int)res; i++) { - for(j = 0; j <= (int)res; j++) { - double tryu = (i/res), tryv = (j/res); - - Vector tryp = PointAt(tryu, tryv); - double d = (tryp.Minus(p)).Magnitude(); - if(d < minDist) { - *u = tryu; - *v = tryv; - minDist = d; + + if(degm == 1 && degn == 1) { + *u = *v = 0; // a plane, perfect no matter what the initial guess + } else { + double minDist = 1e10; + double res = 7.0; + for(i = 0; i < (int)res; i++) { + for(j = 0; j <= (int)res; j++) { + double tryu = (i/res), tryv = (j/res); + + Vector tryp = PointAt(tryu, tryv); + double d = (tryp.Minus(p)).Magnitude(); + if(d < minDist) { + *u = tryu; + *v = tryv; + minDist = d; + } } } } @@ -522,7 +527,7 @@ void SSurface::ClosestPointTo(Vector p, double *u, double *v) { // independently projected into uv and back, to end up equal with // the LENGTH_EPS. Best case that requires LENGTH_EPS/2, but more // is better and convergence should be fast by now. - if(p0.Equals(p, LENGTH_EPS/1e3)) { + if(p0.Equals(p, LENGTH_EPS/1e2)) { return; } @@ -540,6 +545,7 @@ void SSurface::ClosestPointTo(Vector p, double *u, double *v) { dbp("didn't converge"); dbp("have %.3f %.3f %.3f", CO(p0)); dbp("want %.3f %.3f %.3f", CO(p)); + dbp("distance = %g", (p.Minus(p0)).Magnitude()); if(isnan(*u) || isnan(*v)) { *u = *v = 0; } diff --git a/srf/surface.h b/srf/surface.h index cf5a59b..3c9776a 100644 --- a/srf/surface.h +++ b/srf/surface.h @@ -149,8 +149,9 @@ public: class SInter { public: Vector p; - double dot; // between line and surface's normal hSSurface surface; + Vector surfNormal; // of the intersecting surface, at pinter + bool onEdge; // pinter is on edge of trim poly }; // A rational polynomial surface in Bezier form. @@ -187,6 +188,8 @@ public: void TangentsAt(double u, double v, Vector *tu, Vector *tv); Vector NormalAt(double u, double v); void GetAxisAlignedBounding(Vector *ptMax, Vector *ptMin); + bool CoincidentWithPlane(Vector n, double d); + bool CoincidentWith(SSurface *ss, bool sameNormal); void TriangulateInto(SShell *shell, SMesh *sm); void MakeEdgesInto(SShell *shell, SEdgeList *sel, bool asUv); @@ -214,12 +217,15 @@ public: void MakeIntersectionCurvesAgainst(SShell *against, SShell *into); void MakeClassifyingBsps(void); void AllPointsIntersecting(Vector a, Vector b, List *il); + void MakeCoincidentEdgesInto(SSurface *proto, bool sameNormal, + SEdgeList *el); void CleanupAfterBoolean(void); static const int INSIDE = 100; static const int OUTSIDE = 200; - static const int ON_SURFACE = 300; - int ClassifyPoint(Vector p); + static const int ON_PARALLEL = 300; + static const int ON_ANTIPARALLEL = 400; + int ClassifyPoint(Vector p, Vector out); void MakeFromCopyOf(SShell *a); diff --git a/srf/surfinter.cpp b/srf/surfinter.cpp index a433423..c63903b 100644 --- a/srf/surfinter.cpp +++ b/srf/surfinter.cpp @@ -71,13 +71,17 @@ void SSurface::AllPointsIntersecting(Vector a, Vector b, List *l) { // t = (d - (a dot n))/((b - a) dot n) double t = (d - a.Dot(n)) / ((b.Minus(a)).Dot(n)); Vector pi = a.Plus((b.Minus(a)).ScaledBy(t)); + Point2d puv, dummy = { 0, 0 }; ClosestPointTo(pi, &(puv.x), &(puv.y)); - if(bsp->ClassifyPoint(puv, dummy) != SBspUv::OUTSIDE) { + int c = bsp->ClassifyPoint(puv, dummy); + + if(c != SBspUv::OUTSIDE) { SInter si; si.p = pi; - si.dot = NormalAt(puv.x, puv.y).Dot(b.Minus(a)); + si.surfNormal = NormalAt(puv.x, puv.y); si.surface = h; + si.onEdge = (c != SBspUv::INSIDE); l->Add(&si); } } @@ -91,34 +95,143 @@ void SShell::AllPointsIntersecting(Vector a, Vector b, List *il) { } } -int SShell::ClassifyPoint(Vector p) { +//----------------------------------------------------------------------------- +// Does the given point lie on our shell? It might be inside or outside, or +// it might be on the surface with pout parallel or anti-parallel to the +// intersecting surface's normal. +// +// To calculate, we intersect a ray through p with our shell, and classify +// using the closest intersection point. If the ray hits a surface on edge, +// then just reattempt in a different random direction. +//----------------------------------------------------------------------------- +int SShell::ClassifyPoint(Vector p, Vector pout) { List l; ZERO(&l); - Vector d = Vector::From(1e5, 0, 0); // direction is arbitrary - // but it does need to be a one-sided ray - AllPointsIntersecting(p, p.Plus(d), &l); - - double dmin = VERY_POSITIVE; - int ret = OUTSIDE; // no intersections means it's outside + srand(0); - SInter *si; - for(si = l.First(); si; si = l.NextAfter(si)) { - double d = ((si->p).Minus(p)).Magnitude(); + int ret, cnt = 0; + for(;;) { + // Cast a ray in a random direction (two-sided so that we test if + // the point lies on a surface, but use only one side for in/out + // testing) + Vector ray = Vector::From(Random(1), Random(1), Random(1)); + ray = ray.WithMagnitude(1e4); + AllPointsIntersecting(p.Minus(ray), p.Plus(ray), &l); + + double dmin = VERY_POSITIVE; + ret = OUTSIDE; // no intersections means it's outside + bool onEdge = false; - if(d < dmin) { - dmin = d; - if(d < LENGTH_EPS) { - ret = ON_SURFACE; - } else if(si->dot > 0) { - ret = INSIDE; - } else { - ret = OUTSIDE; + SInter *si; + for(si = l.First(); si; si = l.NextAfter(si)) { + double t = ((si->p).Minus(p)).DivPivoting(ray); + if(t*ray.Magnitude() < -LENGTH_EPS) { + // wrong side, doesn't count + continue; } + + double d = ((si->p).Minus(p)).Magnitude(); + + if(d < dmin) { + dmin = d; + if(d < LENGTH_EPS) { + // Lies on the surface + if((si->surfNormal).Dot(pout) > 0) { + ret = ON_PARALLEL; + } else { + ret = ON_ANTIPARALLEL; + } + } else { + // Does not lie on this surface; inside or out, depending + // on the normal + if((si->surfNormal).Dot(ray) > 0) { + ret = INSIDE; + } else { + ret = OUTSIDE; + } + } + onEdge = si->onEdge; + } + } + l.Clear(); + + // If the point being tested lies exactly on an edge of the shell, + // then our ray always lies on edge, and that's okay. + if(ret == ON_PARALLEL || ret == ON_ANTIPARALLEL || !onEdge) break; + if(cnt++ > 10) { + dbp("can't find a ray that doesn't hit on edge!"); + break; } } - l.Clear(); return ret; } +//----------------------------------------------------------------------------- +// Are two surfaces coincident, with the same (or with opposite) normals? +// Currently handles planes only. +//----------------------------------------------------------------------------- +bool SSurface::CoincidentWith(SSurface *ss, bool sameNormal) { + if(degm != 1 || degn != 1) return false; + if(ss->degm != 1 || ss->degn != 1) return false; + + Vector p = ctrl[0][0]; + Vector n = NormalAt(0, 0).WithMagnitude(1); + double d = n.Dot(p); + + if(!ss->CoincidentWithPlane(n, d)) return false; + + Vector n2 = ss->NormalAt(0, 0); + if(sameNormal) { + if(n2.Dot(n) < 0) return false; + } else { + if(n2.Dot(n) > 0) return false; + } + + return true; +} + +bool SSurface::CoincidentWithPlane(Vector n, double d) { + if(degm != 1 || degn != 1) return false; + if(fabs(n.Dot(ctrl[0][0]) - d) > LENGTH_EPS) return false; + if(fabs(n.Dot(ctrl[0][1]) - d) > LENGTH_EPS) return false; + if(fabs(n.Dot(ctrl[1][0]) - d) > LENGTH_EPS) return false; + if(fabs(n.Dot(ctrl[1][1]) - d) > LENGTH_EPS) return false; + + return true; +} + +//----------------------------------------------------------------------------- +// In our shell, find all surfaces that are coincident with the prototype +// surface (with same or opposite normal, as specified), and copy all of +// their trim polygons into el. The edges are returned in uv coordinates for +// the prototype surface. +//----------------------------------------------------------------------------- +void SShell::MakeCoincidentEdgesInto(SSurface *proto, bool sameNormal, + SEdgeList *el) +{ + SSurface *ss; + for(ss = surface.First(); ss; ss = surface.NextAfter(ss)) { + if(proto->CoincidentWith(ss, sameNormal)) { + ss->MakeEdgesInto(this, el, false); + } + } + + SEdge *se; + for(se = el->l.First(); se; se = el->l.NextAfter(se)) { + double ua, va, ub, vb; + proto->ClosestPointTo(se->a, &ua, &va); + proto->ClosestPointTo(se->b, &ub, &vb); + + if(sameNormal) { + se->a = Vector::From(ua, va, 0); + se->b = Vector::From(ub, vb, 0); + } else { + // Flip normal, so flip all edge directions + se->b = Vector::From(ua, va, 0); + se->a = Vector::From(ub, vb, 0); + } + } +} +