Compute surface intersections in a way that is closer to what I

will do for real; now handling the special cases of plane against a
surface of extrusion. Still need to fix up line-surface
intersection to work for curved things, but then some simple curved
cases should work (as well as plane-plane).

[git-p4: depot-paths = "//depot/solvespace/": change = 1919]

dsc.h

@@ -49,6 +49,9 @@ public:
                                         Vector b0, Vector b1,
                                         bool *skew,
                                         double *pa=NULL, double *pb=NULL);
+    static Vector AtIntersectionOfPlaneAndLine(Vector n, double d,
+                                               Vector p0, Vector p1,
+                                               bool *parallel);
 
     double Element(int i);
     bool Equals(Vector v, double tol=LENGTH_EPS);
@@ -93,8 +96,10 @@ public:
     double DistanceTo(Point2d p);
     double DistanceToLine(Point2d p0, Point2d dp, bool segment);
     double Magnitude(void);
+    double MagSquared(void);
     Point2d WithMagnitude(double v);
     Point2d Normal(void);
+    bool Equals(Point2d v, double tol=LENGTH_EPS);
 };
 
 // A simple list

srf/boolean.cpp

@@ -38,8 +38,8 @@ SCurve SCurve::MakeCopySplitAgainst(SShell *agnstA, SShell *agnstB) {
         ZERO(&il);
 
         // Find all the intersections with the two passed shells
-        if(agnstA) agnstA->AllPointsIntersecting(prev, *p, &il);
+        if(agnstA) agnstA->AllPointsIntersecting(prev, *p, &il, true, true);
-        if(agnstB) agnstB->AllPointsIntersecting(prev, *p, &il);
+        if(agnstB) agnstB->AllPointsIntersecting(prev, *p, &il, true, true);
 
         // If any intersections exist, sort them in order along the
         // line and add them to the curve.
@@ -426,7 +426,7 @@ SSurface SSurface::MakeCopyTrimAgainst(SShell *agnst, SShell *parent,
     }
     final.l.RemoveTagged();
 
-//    if(I == 10) DEBUGEDGELIST(&final, &ret);
+//    if(I == 0) DEBUGEDGELIST(&final, &ret);
 
     // Use our reassembled edges to trim the new surface.
     ret.TrimFromEdgeList(&final);
@@ -485,13 +485,13 @@ void SShell::MakeFromBoolean(SShell *a, SShell *b, int type) {
     a->MakeIntersectionCurvesAgainst(b, this);
     
     I = 100;
-    if(b->surface.n == 0 || a->surface.n == 0 || a->surface.n == 6) {
+    if(b->surface.n == 0 || a->surface.n == 0) {
         // Then trim and copy the surfaces
         a->CopySurfacesTrimAgainst(b, this, type, true);
         b->CopySurfacesTrimAgainst(a, this, type, false);
     } else {
-        I = 0;
         a->CopySurfacesTrimAgainst(b, this, type, true);
+        I = -1;
         b->CopySurfacesTrimAgainst(a, this, type, false);
     }
 

srf/ratpoly.cpp

@@ -387,6 +387,33 @@ SSurface SSurface::FromExtrusionOf(SBezier *sb, Vector t0, Vector t1) {
     return ret;
 }
 
+bool SSurface::IsExtrusion(SBezier *of, Vector *alongp) {
+    int i;
+
+    if(degn != 1) return false;
+
+    Vector along = (ctrl[0][1]).Minus(ctrl[0][0]);
+    for(i = 0; i <= degm; i++) {
+        if((fabs(weight[i][1] - weight[i][0]) < LENGTH_EPS) &&
+           ((ctrl[i][1]).Minus(ctrl[i][0])).Equals(along))
+        {
+            continue;
+        }
+        return false;
+    }
+
+    // yes, we are a surface of extrusion; copy the original curve and return
+    if(of) {
+        for(i = 0; i <= degm; i++) {
+            of->weight[i] = weight[i][0];
+            of->ctrl[i] = ctrl[i][0];
+        }
+        of->deg = degm;
+        *alongp = along;
+    }
+    return true;
+}
+
 SSurface SSurface::FromPlane(Vector pt, Vector u, Vector v) {
     SSurface ret;
     ZERO(&ret);

srf/surface.h

@@ -181,7 +181,10 @@ public:
     void TrimFromEdgeList(SEdgeList *el);
     void IntersectAgainst(SSurface *b, SShell *agnstA, SShell *agnstB, 
                           SShell *into);
-    void AllPointsIntersecting(Vector a, Vector b, List<SInter> *l);
+    void AddExactIntersectionCurve(SBezier *sb, hSSurface hsb,
+                          SShell *agnstA, SShell *agnstB, SShell *into);
+    void AllPointsIntersecting(Vector a, Vector b,
+                                    List<SInter> *l, bool seg, bool trimmed);
 
     void ClosestPointTo(Vector p, double *u, double *v);
     Vector PointAt(double u, double v);
@@ -190,6 +193,7 @@ public:
     void GetAxisAlignedBounding(Vector *ptMax, Vector *ptMin);
     bool CoincidentWithPlane(Vector n, double d);
     bool CoincidentWith(SSurface *ss, bool sameNormal);
+    bool IsExtrusion(SBezier *of, Vector *along);
 
     void TriangulateInto(SShell *shell, SMesh *sm);
     void MakeEdgesInto(SShell *shell, SEdgeList *sel, bool asUv);
@@ -217,7 +221,8 @@ public:
     void CopySurfacesTrimAgainst(SShell *against, SShell *into, int t, bool a);
     void MakeIntersectionCurvesAgainst(SShell *against, SShell *into);
     void MakeClassifyingBsps(void);
-    void AllPointsIntersecting(Vector a, Vector b, List<SInter> *il);
+    void AllPointsIntersecting(Vector a, Vector b, List<SInter> *il,
+                                bool seg, bool trimmed);
     void MakeCoincidentEdgesInto(SSurface *proto, bool sameNormal,
                                  SEdgeList *el);
     void CleanupAfterBoolean(void);

srf/surfinter.cpp

@@ -1,5 +1,28 @@
 #include "solvespace.h"
 
+void SSurface::AddExactIntersectionCurve(SBezier *sb, hSSurface hsb,
+                            SShell *agnstA, SShell *agnstB, SShell *into)
+{
+    SCurve sc;
+    ZERO(&sc);
+    sc.surfA = h;
+    sc.surfB = hsb;
+    sb->MakePwlInto(&(sc.pts));
+
+    Vector *prev = NULL, *v;
+    for(v = sc.pts.First(); v; v = sc.pts.NextAfter(v)) {
+      if(prev) SS.nakedEdges.AddEdge(*prev, *v);
+        prev = v;
+    }
+
+    // Now split the line where it intersects our existing surfaces
+    SCurve split = sc.MakeCopySplitAgainst(agnstA, agnstB);
+    sc.Clear();
+
+    split.interCurve = true;
+    into->curve.AddAndAssignId(&split);
+}
+
 void SSurface::IntersectAgainst(SSurface *b, SShell *agnstA, SShell *agnstB, 
                                 SShell *into)
 {
@@ -12,86 +35,174 @@ void SSurface::IntersectAgainst(SSurface *b, SShell *agnstA, SShell *agnstB,
         return;
     }
 
-    if(degm == 1 && degn == 1 && b->degm == 1 && b->degn == 1) {
+    if((degm == 1 && degn == 1 && b->IsExtrusion(NULL, NULL)) ||
-        // Plane-plane intersection, easy; result is a line
+       (b->degm == 1 && b->degn == 1 && this->IsExtrusion(NULL, NULL)))
-        Vector pta = ctrl[0][0], ptb = b->ctrl[0][0];
+    {
-        Vector na = NormalAt(0, 0), nb = b->NormalAt(0, 0);
+        // The intersection between a plane and a surface of extrusion
-        na = na.WithMagnitude(1);
+        SSurface *splane, *sext;
-        nb = nb.WithMagnitude(1);
+        if(degm == 1 && degn == 1) {
-
+            splane = this;
-        Vector d = (na.Cross(nb));
+            sext = b;
-
+        } else {
-        if(d.Magnitude() < LENGTH_EPS) {
+            splane = b;
-            // parallel planes, no intersection
+            sext = this;
-            return;
         }
 
-        Vector inter = Vector::AtIntersectionOfPlanes(na, na.Dot(pta),
+        Vector n = splane->NormalAt(0, 0).WithMagnitude(1), along;
-                                                      nb, nb.Dot(ptb));
+        double d = n.Dot(splane->PointAt(0, 0));
+        SBezier bezier;
+        (void)sext->IsExtrusion(&bezier, &along);
 
-        // The intersection curve can't be longer than the longest curve
+        if(fabs(n.Dot(along)) < LENGTH_EPS) {
-        // that lies in both planes, which is the diagonal of the shorter;
+            // Direction of extrusion is parallel to plane; so intersection
-        // so just pick one, and then give some slop, not critical.
+            // is zero or more lines. Build a line within the plane, and
-        double maxl = ((ctrl[0][0]).Minus(ctrl[1][1])).Magnitude();
+            // normal to the direction of extrusion, and intersect that line
+            // against the surface; each intersection point corresponds to
+            // a line.
+            Vector pm, alu, p0, dp;
+            // a point halfway along the extrusion
+            pm = ((sext->ctrl[0][0]).Plus(sext->ctrl[0][1])).ScaledBy(0.5);
+            alu = along.WithMagnitude(1);
+            dp = (n.Cross(along)).WithMagnitude(1);
+            // n, alu, and dp form an orthogonal csys; set n component to
+            // place it on the plane, alu component to lie halfway along
+            // extrusion, and dp component doesn't matter so zero
+            p0 = n.ScaledBy(d).Plus(alu.ScaledBy(pm.Dot(alu)));
 
-        Vector v;
+            List<SInter> inters;
-        SCurve sc;
+            ZERO(&inters);
-        ZERO(&sc);
+            sext->AllPointsIntersecting(p0, p0.Plus(dp), &inters, false, false);
-        sc.surfA = h;
-        sc.surfB = b->h;
-        v = inter.Minus(d.WithMagnitude(5*maxl));
-        sc.pts.Add(&v);
-        v = inter.Plus(d.WithMagnitude(5*maxl));
-        sc.pts.Add(&v);
     
-        // Now split the line where it intersects our existing surfaces
+            SInter *si;
-        SCurve split = sc.MakeCopySplitAgainst(agnstA, agnstB);
+            for(si = inters.First(); si; si = inters.NextAfter(si)) {
-        sc.Clear();
+                Vector al = along.ScaledBy(0.5);
+                SBezier bezier;
+                bezier = SBezier::From((si->p).Minus(al), (si->p).Plus(al));
 
-        split.interCurve = true;
+                AddExactIntersectionCurve(&bezier, b->h, agnstA, agnstB, into);
-        into->curve.AddAndAssignId(&split);
+            }
+
+            inters.Clear();
+        } else {
+            // Direction of extrusion is not parallel to plane; so
+            // intersection is projection of extruded curve into our plane
+
+            int i;
+            for(i = 0; i <= bezier.deg; i++) {
+                Vector p0 = bezier.ctrl[i],
+                       p1 = p0.Plus(along);
+
+                bezier.ctrl[i] =
+                    Vector::AtIntersectionOfPlaneAndLine(n, d, p0, p1, NULL);
+            }
+
+            AddExactIntersectionCurve(&bezier, b->h, agnstA, agnstB, into);
+        }
     }
 
     // need to implement general numerical surface intersection for tough
     // cases, just giving up for now
 }
 
-void SSurface::AllPointsIntersecting(Vector a, Vector b, List<SInter> *l) {
+//-----------------------------------------------------------------------------
-    if(degm == 1 && degn == 1) {
+// Find all points where a line through a and b intersects our surface, and
-        // line-plane intersection
+// add them to the list. If seg is true then report only intersections that
-        Vector p = ctrl[0][0];
+// lie within the finite line segment (not including the endpoints); otherwise
-        Vector n = NormalAt(0, 0).WithMagnitude(1);
+// we work along the infinite line.
-        double d = n.Dot(p);
+//-----------------------------------------------------------------------------
-        if((n.Dot(a) - d < -LENGTH_EPS && n.Dot(b) - d > LENGTH_EPS) ||
+void SSurface::AllPointsIntersecting(Vector a, Vector b,
-           (n.Dot(b) - d < -LENGTH_EPS && n.Dot(a) - d > LENGTH_EPS))
+                                     List<SInter> *l, bool seg, bool trimmed)
-        {
+{
-            // It crosses the plane, one point of intersection
+    Vector ba = b.Minus(a);
-            // (a + t*(b - a)) dot n = d
+    double bam = ba.Magnitude();
-            // (a dot n) + t*((b - a) dot n) = d
-            // 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 };
+    typedef struct {
-            ClosestPointTo(pi, &(puv.x), &(puv.y));
+        int     tag;
-            int c = bsp->ClassifyPoint(puv, dummy);
+        Point2d p;
+    } Inter;
+    List<Inter> inters;
+    ZERO(&inters);
 
-            if(c != SBspUv::OUTSIDE) {
+    // First, get all the intersections between the infinite ray and the
-                SInter si;
+    // untrimmed surface.
-                si.p = pi;
+    int i, j;
-                si.surfNormal = NormalAt(puv.x, puv.y);
+    for(i = 0; i < degm; i++) {
-                si.surface = h;
+        for(j = 0; j < degn; j++) {
-                si.onEdge = (c != SBspUv::INSIDE);
+            // Intersect the ray with each face in the control polyhedron
-                l->Add(&si);
+            Vector p00 = ctrl[i][j],
+                   p01 = ctrl[i][j+1],
+                   p10 = ctrl[i+1][j];
+
+            Vector u = p01.Minus(p00), v = p10.Minus(p00);
+
+            Vector n = (u.Cross(v)).WithMagnitude(1);
+            double d = n.Dot(p00);
+
+            bool parallel;
+            Vector pi =
+                Vector::AtIntersectionOfPlaneAndLine(n, d, a, b, &parallel);
+            if(parallel) continue;
+
+            double ui = (pi.Minus(p00)).Dot(u) / u.MagSquared(),
+                   vi = (pi.Minus(p00)).Dot(v) / v.MagSquared();
+
+            double tol = 1e-2;
+            if(ui < -tol || ui > 1 + tol || vi < -tol || vi > 1 + tol) {
+                continue;
+            }
+
+            Inter inter;
+            ClosestPointTo(pi, &inter.p.x, &inter.p.y);
+            inters.Add(&inter);
+        }
+    }
+
+    // Remove duplicate intersection points
+    inters.ClearTags();
+    for(i = 0; i < inters.n; i++) {
+        for(j = i + 1; j < inters.n; j++) {
+            if(inters.elem[i].p.Equals(inters.elem[j].p)) {
+                inters.elem[j].tag = 1;
             }
         }
     }
+    inters.RemoveTagged();
+
+    for(i = 0; i < inters.n; i++) {
+        Point2d puv = inters.elem[i].p;
+
+        // Make sure the point lies within the finite line segment
+        Vector pxyz = PointAt(puv.x, puv.y);
+        double t = (pxyz.Minus(a)).DivPivoting(ba);
+        if(seg && (t > 1 - LENGTH_EPS/bam || t < LENGTH_EPS/bam)) {
+            continue;
+        }
+       
+        // And that it lies inside our trim region
+        Point2d dummy = { 0, 0 };
+        int c = bsp->ClassifyPoint(puv, dummy);
+        if(trimmed && c == SBspUv::OUTSIDE) {
+            continue;
+        }
+
+        // It does, so generate the intersection
+        SInter si;
+        si.p = pxyz;
+        si.surfNormal = NormalAt(puv.x, puv.y);
+        si.surface = h;
+        si.onEdge = (c != SBspUv::INSIDE);
+        l->Add(&si);
+    }
+
+    inters.Clear();
 }
 
-void SShell::AllPointsIntersecting(Vector a, Vector b, List<SInter> *il) {
+void SShell::AllPointsIntersecting(Vector a, Vector b,
+                                   List<SInter> *il, bool seg, bool trimmed)
+{
     SSurface *ss;
     for(ss = surface.First(); ss; ss = surface.NextAfter(ss)) {
-        ss->AllPointsIntersecting(a, b, il);
+        ss->AllPointsIntersecting(a, b, il, seg, trimmed);
     }
 }
 
@@ -118,8 +229,7 @@ int SShell::ClassifyPoint(Vector p, Vector pout) {
         // 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, false, true);
-        AllPointsIntersecting(p.Minus(ray), p.Plus(ray), &l);
         
         double dmin = VERY_POSITIVE;
         ret = OUTSIDE; // no intersections means it's outside

util.cpp

@@ -620,6 +620,26 @@ Vector Vector::AtIntersectionOfLines(Vector a0, Vector a1,
     return pi;
 }
 
+Vector Vector::AtIntersectionOfPlaneAndLine(Vector n, double d,
+                                            Vector p0, Vector p1,
+                                            bool *parallel)
+{
+    Vector dp = p1.Minus(p0);
+
+    if(fabs(n.Dot(dp)) < LENGTH_EPS) {
+        if(parallel) *parallel = true;
+        return Vector::From(0, 0, 0);
+    }
+
+    if(parallel) *parallel = false;
+
+    // n dot (p0 + t*dp) = d
+    // (n dot p0) + t * (n dot dp) = d
+    double t = (d - n.Dot(p0)) / (n.Dot(dp));
+
+    return p0.Plus(dp.ScaledBy(t));
+}
+
 Point2d Point2d::Plus(Point2d b) {
     Point2d r;
     r.x = x + b.x;
@@ -641,6 +661,10 @@ Point2d Point2d::ScaledBy(double s) {
     return r;
 }
 
+double Point2d::MagSquared(void) {
+    return x*x + y*y;
+}
+
 double Point2d::Magnitude(void) {
     return sqrt(x*x + y*y);
 }
@@ -691,3 +715,10 @@ Point2d Point2d::Normal(void) {
     return ret;
 }
 
+bool Point2d::Equals(Point2d v, double tol) {
+    double dx = v.x - x; if(dx < -tol || dx > tol) return false;
+    double dy = v.y - y; if(dy < -tol || dy > tol) return false;
+
+    return (this->Minus(v)).MagSquared() < tol*tol;
+}
+