NURBS boolean improvement.

The 3-plane intersection code fails to converge when a curve joins two tangent NURBS patches. This adds a new function for intersecting exact curves with a surface to avoid those failures. Fixes simplified test model for issue #315.

src/srf/boolean.cpp

@@ -88,8 +88,15 @@ SCurve SCurve::MakeCopySplitAgainst(SShell *agnstA, SShell *agnstB,
                     }
                 }
 
-                // We're keeping the intersection, so actually refine it.
-                (pi->srf)->PointOnSurfaces(srfA, srfB, &(puv.x), &(puv.y));
+                // We're keeping the intersection, so actually refine it. Finding the intersection
+                // to within EPS is important to match the ends of different chopped trim curves.
+                // The general 3-surface intersection fails to refine for trims where surfaces
+                // are tangent at the curve, but those trims are usually exact, so…
+                if(isExact) {
+                    (pi->srf)->PointOnCurve(&exact, &(puv.x), &(puv.y));
+                } else {
+                    (pi->srf)->PointOnSurfaces(srfA, srfB, &(puv.x), &(puv.y));
+                }
                 pi->p = (pi->srf)->PointAt(puv);
             }
             il.RemoveTagged();

src/srf/ratpoly.cpp

@@ -634,3 +634,58 @@ void SSurface::PointOnSurfaces(SSurface *s1, SSurface *s2, double *up, double *v
     dbp("didn't converge (three surfaces intersecting)");
 }
 
+void SSurface::PointOnCurve(const SBezier *curve, double *up, double *vp)
+{
+    Vector tu,tv,n;
+    double u = *up, v = *vp;
+    Vector ps = PointAt(u, v);
+    // Get initial guesses for t on the curve
+    double tCurve = 0.5;
+    curve->ClosestPointTo(ps, &tCurve, /*mustConverge=*/false);
+    if(tCurve < 0.0) tCurve = 0.0;
+    if(tCurve > 1.0) tCurve = 1.0;
+
+    for(int i = 0; i < 30; i++) {
+        // Approximate the surface by a plane
+        Vector ps = PointAt(u, v);
+        TangentsAt(u, v, &tu, &tv);
+        n = tu.Cross(tv).WithMagnitude(1);
+
+        // point on curve and tangent line direction
+        Vector pc = curve->PointAt(tCurve);
+        Vector tc = curve->TangentAt(tCurve);
+
+        if(ps.Equals(pc, RATPOLY_EPS)) {
+            *up = u;
+            *vp = v;
+            return;
+        }
+
+        //pi is where the curve tangent line intersects the surface tangent plane
+        Vector pi;
+        double d = tc.Dot(n);
+        if (fabs(d) < 1e-10) { // parallel line and plane, guess the average rather than fail
+            pi = pc.Plus(ps).ScaledBy(0.5);
+        } else {
+            pi = pc.Minus(tc.ScaledBy(pc.Minus(ps).Dot(n)/d));
+        }
+
+        // project the point onto the tangent plane and line
+        {
+            Vector n = tu.Cross(tv);
+            Vector ty = n.Cross(tu).ScaledBy(1.0/tu.MagSquared());
+            Vector tx = tv.Cross(n).ScaledBy(1.0/tv.MagSquared());
+
+            Vector dp = pi.Minus(ps);
+            double du = dp.Dot(tx), dv = dp.Dot(ty);
+
+            u += du / tx.MagSquared();
+            v += dv / ty.MagSquared();
+        }
+        tCurve += pi.Minus(pc).Dot(tc) / tc.MagSquared();
+        if(tCurve < 0.0) tCurve = 0.0;
+        if(tCurve > 1.0) tCurve = 1.0;
+    }
+    dbp("didn't converge (surface and curve intersecting)");
+}
+

src/srf/surface.h

@@ -332,6 +332,7 @@ public:
     bool PointIntersectingLine(Vector p0, Vector p1, double *u, double *v) const;
     Vector ClosestPointOnThisAndSurface(SSurface *srf2, Vector p);
     void PointOnSurfaces(SSurface *s1, SSurface *s2, double *u, double *v);
+    void PointOnCurve(const SBezier *curve, double *up, double *vp);
     Vector PointAt(double u, double v) const;
     Vector PointAt(Point2d puv) const;
     void TangentsAt(double u, double v, Vector *tu, Vector *tv) const;

src/srf/surfinter.cpp

@@ -341,7 +341,11 @@ void SSurface::IntersectAgainst(SSurface *b, SShell *agnstA, SShell *agnstB,
                     Vector p = si->p;
                     double u, v;
                     srfB->ClosestPointTo(p, &u, &v);
-                    srfB->PointOnSurfaces(srfA, other, &u, &v);
+                    if(sc->isExact) {
+                        srfB->PointOnCurve(&(sc->exact), &u, &v);
+                    } else {
+                        srfB->PointOnSurfaces(srfA, other, &u, &v);
+                    }
                     p = srfB->PointAt(u, v);
                     if(!spl.ContainsPoint(p)) {
                         SPoint sp;