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

@@ -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 <gl/gl.h>
 #include <gl/glu.h>
 
+inline double Random(double vmax) {
+    return (vmax*rand()) / RAND_MAX;
+}
+
 class Expr;
 class ExprVector;
 class ExprQuaternion;

srf/boolean.cpp

@@ -48,9 +48,15 @@ SCurve SCurve::MakeCopySplitAgainst(SShell *agnstA, SShell *agnstB) {
             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),
+        Point2d auv = (se->a).ProjectXy(),
-               eb = PointAt(se->b.x, se->b.y);
+                buv = (se->b).ProjectXy();
-        int c = agnst->ClassifyPoint(ea.Plus(eb).ScaledBy(0.5));
 
-        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),
-                SS.nakedEdges.AddEdge(PointAt(se->a.x, se->a.y),
+                                  PointAt(se->b.x, se->b.y));
-                                      PointAt(se->b.x, se->b.y));
+            SS.nakedEdges.AddEdge(PointAt(mid.x, mid.y),
-                SS.nakedEdges.AddEdge(PointAt(mid.x, mid.y),
+                                  PointAt(arrow.x, arrow.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();

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);
+    if(degm == 1 && degn == 1) {
-            double d = (tryp.Minus(p)).Magnitude();
+        *u = *v = 0; // a plane, perfect no matter what the initial guess
-            if(d < minDist) {
+    } else {
-                *u = tryu;
+        double minDist = 1e10;
-                *v = tryv;
+        double res = 7.0;
-                minDist = d;
+        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;
     }

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<SInter> *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;
+    static const int ON_PARALLEL       = 300;
-    int ClassifyPoint(Vector p);
+    static const int ON_ANTIPARALLEL   = 400;
+    int ClassifyPoint(Vector p, Vector out);
 
 
     void MakeFromCopyOf(SShell *a);

srf/surfinter.cpp

@@ -71,13 +71,17 @@ void SSurface::AllPointsIntersecting(Vector a, Vector b, List<SInter> *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<SInter> *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<SInter> l;
     ZERO(&l);
 
-    Vector d = Vector::From(1e5, 0, 0); // direction is arbitrary
+    srand(0);
-    // but it does need to be a one-sided ray
-    AllPointsIntersecting(p, p.Plus(d), &l);
 
-    double dmin = VERY_POSITIVE;
+    int ret, cnt = 0;
-    int ret = OUTSIDE; // no intersections means it's outside
+    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);
         
-    SInter *si;
+        double dmin = VERY_POSITIVE;
-    for(si = l.First(); si; si = l.NextAfter(si)) {
+        ret = OUTSIDE; // no intersections means it's outside
-        double d = ((si->p).Minus(p)).Magnitude();
+        bool onEdge = false;
 
-        if(d < dmin) {
+        SInter *si;
-            dmin = d;
+        for(si = l.First(); si; si = l.NextAfter(si)) {
-            if(d < LENGTH_EPS) {
+            double t = ((si->p).Minus(p)).DivPivoting(ray);
-                ret = ON_SURFACE;
+            if(t*ray.Magnitude() < -LENGTH_EPS) {
-            } else if(si->dot > 0) {
+                // wrong side, doesn't count
-                ret = INSIDE;
+                continue;
-            } else {
-                ret = OUTSIDE;
             }
+
+            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);
+        }
+    }
+}
+