More work on Booleans. This works only for planes, and only for

non-coincident faces. There's also a problem when I don't generate
the full intersection polygon of shell B against a given surface in
shell A; I need to modify the code to not require that.

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

dsc.h

@@ -79,6 +79,7 @@ public:
                                       Vector bmax, Vector bmin);
     bool OutsideAndNotOn(Vector maxv, Vector minv);
     Point2d Project2d(Vector u, Vector v);
+    Point2d ProjectXy(void);
 };
 
 class Point2d {
@@ -88,10 +89,12 @@ public:
     Point2d Plus(Point2d b);
     Point2d Minus(Point2d b);
     Point2d ScaledBy(double s);
+    double Dot(Point2d p);
     double DistanceTo(Point2d p);
     double DistanceToLine(Point2d p0, Point2d dp, bool segment);
     double Magnitude(void);
     Point2d WithMagnitude(double v);
+    Point2d Normal(void);
 };
 
 // A simple list

srf/boolean.cpp

@@ -1,5 +1,7 @@
 #include "solvespace.h"
 
+static int I, N;
+
 void SShell::MakeFromUnionOf(SShell *a, SShell *b) {
     MakeFromBoolean(a, b, AS_UNION);
 }
@@ -8,38 +10,66 @@ void SShell::MakeFromDifferenceOf(SShell *a, SShell *b) {
     MakeFromBoolean(a, b, AS_DIFFERENCE);
 }
 
-SCurve SCurve::MakeCopySplitAgainst(SShell *against) {
+static Vector LineStart, LineDirection;
+static int ByTAlongLine(const void *av, const void *bv)
+{
+    Vector *a = (Vector *)av,
+           *b = (Vector *)bv;
+    
+    double ta = (a->Minus(LineStart)).DivPivoting(LineDirection),
+           tb = (b->Minus(LineStart)).DivPivoting(LineDirection);
+
+    return (ta > tb) ? 1 : -1;
+}
+SCurve SCurve::MakeCopySplitAgainst(SShell *agnstA, SShell *agnstB) {
     SCurve ret;
     ret = *this;
+    ret.interCurve = false;
     ZERO(&(ret.pts));
 
-    Vector *p;
+    Vector *p = pts.First();
-    for(p = pts.First(); p; p = pts.NextAfter(p)) {
+    if(!p) oops();
+    Vector prev = *p;
+    ret.pts.Add(p);
+    p = pts.NextAfter(p);
+
+    for(; p; p = pts.NextAfter(p)) {
+        List<Vector> il;
+        ZERO(&il);
+
+        // Find all the intersections with the two passed shells
+        if(agnstA) agnstA->AllPointsIntersecting(prev, *p, &il);
+        if(agnstB) agnstB->AllPointsIntersecting(prev, *p, &il);
+
+        // If any intersections exist, sort them in order along the
+        // line and add them to the curve.
+        if(il.n > 0) {
+            LineStart = prev;
+            LineDirection = p->Minus(prev);
+            qsort(il.elem, il.n, sizeof(il.elem[0]), ByTAlongLine);
+
+            Vector *pi;
+            for(pi = il.First(); pi; pi = il.NextAfter(pi)) {
+                ret.pts.Add(pi);
+            }
+        }
+
         ret.pts.Add(p);
+        prev = *p;
     }
     return ret;
 }
 
-void SShell::CopyCurvesSplitAgainst(SShell *against, SShell *into) {
+void SShell::CopyCurvesSplitAgainst(SShell *aga, SShell *agb, SShell *into) {
     SCurve *sc;
     for(sc = curve.First(); sc; sc = curve.NextAfter(sc)) {
-        SCurve scn = sc->MakeCopySplitAgainst(against);
+        SCurve scn = sc->MakeCopySplitAgainst(aga, agb);
         hSCurve hsc = into->curve.AddAndAssignId(&scn);
         // And note the new ID so that we can rewrite the trims appropriately
         sc->newH = hsc;
     }
 }
 
-void SShell::MakeEdgeListUseNewCurveIds(SEdgeList *el) {
-    SEdge *se;
-    for(se = el->l.First(); se; se = el->l.NextAfter(se)) {
-        hSCurve oldh = { se->auxA };
-        SCurve *osc = curve.FindById(oldh);
-        se->auxA = osc->newH.v;
-        // auxB is the direction, which is unchanged
-    }
-}
-
 void SSurface::TrimFromEdgeList(SEdgeList *el) {
     el->l.ClearTags();
    
@@ -95,7 +125,8 @@ void SSurface::TrimFromEdgeList(SEdgeList *el) {
 // also need a pointer to the shell that contains our own surface, since that
 // contains our original trim curves.
 //-----------------------------------------------------------------------------
-SSurface SSurface::MakeCopyTrimAgainst(SShell *against, SShell *shell,
+SSurface SSurface::MakeCopyTrimAgainst(SShell *agnst, SShell *parent,
+                                       SShell *into,
                                        int type, bool opA)
 {
     SSurface ret;
@@ -103,14 +134,121 @@ SSurface SSurface::MakeCopyTrimAgainst(SShell *against, SShell *shell,
     ret = *this;
     ZERO(&(ret.trim));
 
-    SEdgeList el;
+    // First, build a list of the existing trim curves; update them to use
-    ZERO(&el);
+    // the split curves.
-    MakeEdgesInto(shell, &el, true);
+    STrimBy *stb;
-    shell->MakeEdgeListUseNewCurveIds(&el);
+    for(stb = trim.First(); stb; stb = trim.NextAfter(stb)) {
+        STrimBy stn = *stb;
+        stn.curve = (parent->curve.FindById(stn.curve))->newH;
+        ret.trim.Add(&stn);
+    }
 
-    ret.TrimFromEdgeList(&el);
+    // Build up our original trim polygon
+    SEdgeList orig;
+    ZERO(&orig);
+    ret.MakeEdgesInto(into, &orig, true);
+    ret.trim.Clear();
 
-    el.Clear();
+
+    // And now intersect the other shell against us
+    SEdgeList inter;
+    ZERO(&inter);
+
+    SSurface *ss;
+    for(ss = agnst->surface.First(); ss; ss = agnst->surface.NextAfter(ss)) {
+        SCurve *sc;
+        for(sc = into->curve.First(); sc; sc = into->curve.NextAfter(sc)) {
+            if(!(sc->interCurve)) continue;
+            if(opA) {
+                if(sc->surfB.v != h.v || sc->surfA.v != ss->h.v) continue;
+            } else {
+                if(sc->surfA.v != h.v || sc->surfB.v != ss->h.v) continue;
+            }
+           
+            int i;
+            for(i = 1; i < sc->pts.n; i++) {
+                Vector a = sc->pts.elem[i-1],
+                       b = sc->pts.elem[i];
+
+                Point2d auv, buv;
+                ss->ClosestPointTo(a, &(auv.x), &(auv.y));
+                ss->ClosestPointTo(b, &(buv.x), &(buv.y));
+
+                int c = ss->bsp->ClassifyEdge(auv, buv);
+                if(c == SBspUv::INSIDE) {
+                    Vector ta = Vector::From(0, 0, 0);
+                    Vector tb = Vector::From(0, 0, 0);
+                    ClosestPointTo(a, &(ta.x), &(ta.y));
+                    ClosestPointTo(b, &(tb.x), &(tb.y));
+
+                    Vector tn = NormalAt(ta.x, ta.y);
+                    Vector sn = ss->NormalAt(auv.x, auv.y);
+
+                    if((tn.Cross(b.Minus(a))).Dot(sn) > 0) {
+                        inter.AddEdge(ta, tb, sc->h.v, 0);
+                    } else { 
+                        inter.AddEdge(tb, ta, sc->h.v, 1);
+                    }
+                }
+            }
+        }
+    }
+
+    SEdgeList final;
+    ZERO(&final);
+
+    if(I == 2) dbp("INTERBSP: %d", inter.l.n);
+    SBspUv *interbsp = SBspUv::From(&inter);
+    if(I == 2) dbp("INTEROVER");
+
+
+    SEdge *se;
+    N = 0;
+    for(se = orig.l.First(); se; se = orig.l.NextAfter(se)) {
+        int c = interbsp->ClassifyEdge(se->a.ProjectXy(), se->b.ProjectXy());
+
+        if(I == 2) dbp("edge from %.3f %.3f %.3f to %.3f %.3f %.3f",
+            CO(PointAt(se->a.x, se->a.y)), CO(PointAt(se->b.x, se->b.y)));
+
+        if(c == SBspUv::OUTSIDE) {
+            if(I == 2) dbp("   keep");
+            final.AddEdge(se->a, se->b, se->auxA, se->auxB);
+        } else {
+            if(I == 2) dbp("    don't keep, %d", c);
+        }
+        N++;
+    }
+
+    for(se = inter.l.First(); se; se = inter.l.NextAfter(se)) {
+
+        if(I == 2) {
+            Vector mid = (se->a).Plus(se->b).ScaledBy(0.5);
+            Vector arrow = (se->b).Minus(se->a);
+            SWAP(double, arrow.x, arrow.y);
+            arrow.x *= -1;
+            arrow = arrow.WithMagnitude(0.03);
+            arrow = arrow.Plus(mid);
+
+            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));
+        }
+
+        int c = bsp->ClassifyEdge(se->a.ProjectXy(), se->b.ProjectXy());
+        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);
+
+    final.Clear();
+    inter.Clear();
+    orig.Clear();
     return ret;
 }
 
@@ -120,8 +258,9 @@ void SShell::CopySurfacesTrimAgainst(SShell *against, SShell *into,
     SSurface *ss;
     for(ss = surface.First(); ss; ss = surface.NextAfter(ss)) {
         SSurface ssn;
-        ssn = ss->MakeCopyTrimAgainst(against, this, type, opA);
+        ssn = ss->MakeCopyTrimAgainst(against, this, into, type, opA);
         into->surface.AddAndAssignId(&ssn);
+        I++;
     }
 }
 
@@ -133,7 +272,7 @@ void SShell::MakeIntersectionCurvesAgainst(SShell *agnst, SShell *into) {
             // Intersect every surface from our shell against every surface
             // from agnst; this will add zero or more curves to the curve
             // list for into.
-            sa->IntersectAgainst(sb, into);
+            sa->IntersectAgainst(sb, agnst, this, into);
         }
     }
 }
@@ -141,30 +280,194 @@ void SShell::MakeIntersectionCurvesAgainst(SShell *agnst, SShell *into) {
 void SShell::CleanupAfterBoolean(void) {
     SSurface *ss;
     for(ss = surface.First(); ss; ss = surface.NextAfter(ss)) {
-        (ss->orig).Clear();
-        (ss->inside).Clear();
-        (ss->onSameNormal).Clear();
-        (ss->onFlipNormal).Clear();
     }
 }
 
 void SShell::MakeFromBoolean(SShell *a, SShell *b, int type) {
+    a->MakeClassifyingBsps();
+    b->MakeClassifyingBsps();
+
     // Copy over all the original curves, splitting them so that a
     // piecwise linear segment never crosses a surface from the other
     // shell.
-    a->CopyCurvesSplitAgainst(b, this);
+    a->CopyCurvesSplitAgainst(b, NULL, this);
-    b->CopyCurvesSplitAgainst(a, this);
+    b->CopyCurvesSplitAgainst(a, NULL, this);
 
     // Generate the intersection curves for each surface in A against all
-    // the surfaces in B
+    // the surfaces in B (which is all of the intersection curves).
     a->MakeIntersectionCurvesAgainst(b, this);
 
-    // Then trim and copy the surfaces
+    if(a->surface.n == 0 || b->surface.n == 0) {
-    a->CopySurfacesTrimAgainst(b, this, type, true);
+        // Then trim and copy the surfaces
-    b->CopySurfacesTrimAgainst(a, this, type, false);
+        I = 100;
+        a->CopySurfacesTrimAgainst(b, this, type, true);
+        b->CopySurfacesTrimAgainst(a, this, type, false);
+    } else {
+        I = -1;
+        a->CopySurfacesTrimAgainst(b, this, type, true);
+        b->CopySurfacesTrimAgainst(a, this, type, false);
+    }
 
     // And clean up the piecewise linear things we made as a calculation aid
     a->CleanupAfterBoolean();
     b->CleanupAfterBoolean();
 }
 
+//-----------------------------------------------------------------------------
+// All of the BSP routines that we use to perform and accelerate polygon ops.
+//-----------------------------------------------------------------------------
+void SShell::MakeClassifyingBsps(void) {
+    SSurface *ss;
+    for(ss = surface.First(); ss; ss = surface.NextAfter(ss)) {
+        ss->MakeClassifyingBsp(this);
+    }
+}
+
+void SSurface::MakeClassifyingBsp(SShell *shell) {
+    SEdgeList el;
+    ZERO(&el);
+
+    MakeEdgesInto(shell, &el, true);
+    bsp = SBspUv::From(&el);
+
+    el.Clear();
+}
+
+SBspUv *SBspUv::Alloc(void) {
+    return (SBspUv *)AllocTemporary(sizeof(SBspUv));
+}
+
+static int ByLength(const void *av, const void *bv)
+{
+    SEdge *a = (SEdge *)av,
+          *b = (SEdge *)bv;
+    
+    double la = (a->a).Minus(a->b).Magnitude(),
+           lb = (b->a).Minus(b->b).Magnitude();
+
+    // Sort in descending order, longest first. This improves numerical
+    // stability for the normals.
+    return (la < lb) ? 1 : -1;
+}
+SBspUv *SBspUv::From(SEdgeList *el) {
+    SEdgeList work;
+    ZERO(&work);
+
+    SEdge *se;
+    for(se = el->l.First(); se; se = el->l.NextAfter(se)) {
+        work.AddEdge(se->a, se->b, se->auxA, se->auxB);
+    }
+    qsort(work.l.elem, work.l.n, sizeof(work.l.elem[0]), ByLength);
+
+    SBspUv *bsp = NULL;
+    for(se = work.l.First(); se; se = work.l.NextAfter(se)) {
+        bsp = bsp->InsertEdge((se->a).ProjectXy(), (se->b).ProjectXy());
+    }
+
+    work.Clear();
+    return bsp;
+}
+
+SBspUv *SBspUv::InsertEdge(Point2d ea, Point2d eb) {
+    if(I == 2) {
+        dbp("insert edge %.3f %.3f to %.3f %.3f", ea.x, ea.y, eb.x, eb.y);
+    }
+
+    if(!this) {
+        SBspUv *ret = Alloc();
+        ret->a = ea;
+        ret->b = eb;
+        return ret;
+    }
+
+    Point2d n = ((b.Minus(a)).Normal()).WithMagnitude(1);
+    double d = a.Dot(n);
+
+    double dea = ea.Dot(n) - d,
+           deb = eb.Dot(n) - d;
+
+    if(fabs(dea) < LENGTH_EPS && fabs(deb) < LENGTH_EPS) {
+        // Line segment is coincident with this one, store in same node
+        SBspUv *m = Alloc();
+        m->a = ea;
+        m->b = eb;
+        m->more = more;
+        more = m;
+    } else if(fabs(dea) < LENGTH_EPS) {
+        // Point A lies on this lie, but point B does not
+        if(deb > 0) {
+            pos = pos->InsertEdge(ea, eb);
+        } else {
+            neg = neg->InsertEdge(ea, eb);
+        }
+    } else if(fabs(deb) < LENGTH_EPS) {
+        // Point B lies on this lie, but point A does not
+        if(dea > 0) {
+            pos = pos->InsertEdge(ea, eb);
+        } else {
+            neg = neg->InsertEdge(ea, eb);
+        }
+    } else if(dea > 0 && deb > 0) {
+        pos = pos->InsertEdge(ea, eb);
+    } else if(dea < 0 && deb < 0) {
+        neg = neg->InsertEdge(ea, eb);
+    } else {
+        // New edge crosses this one; we need to split.
+        double t = (d - n.Dot(ea)) / (n.Dot(eb.Minus(ea)));
+        Point2d pi = ea.Plus((eb.Minus(ea)).ScaledBy(t));
+        if(dea > 0) {
+            pos = pos->InsertEdge(ea, pi);
+            neg = neg->InsertEdge(pi, eb);
+        } else {
+            neg = neg->InsertEdge(ea, pi);
+            pos = pos->InsertEdge(pi, eb);
+        }
+    }
+    return this;
+}
+
+int SBspUv::ClassifyPoint(Point2d p, Point2d eb) {
+    if(!this) return OUTSIDE;
+
+    Point2d n = ((b.Minus(a)).Normal()).WithMagnitude(1);
+    double d = a.Dot(n);
+
+    double dp = p.Dot(n) - d;
+
+    if(I == 2 && N == 5) {
+        dbp("point %.3f %.3f has d=%.3f", p.x, p.y, dp);
+    }
+
+    if(fabs(dp) < LENGTH_EPS) {
+        if(I == 2 && N == 5) dbp("   on line");
+        SBspUv *f = this;
+        while(f) {
+            Point2d ba = (f->b).Minus(f->a);
+            if(p.DistanceToLine(f->a, ba, true) < LENGTH_EPS) {
+                if(eb.DistanceToLine(f->a, ba, false) < LENGTH_EPS) {
+                    if(ba.Dot(eb.Minus(p)) > 0) {
+                        return EDGE_PARALLEL;
+                    } else {
+                        return EDGE_ANTIPARALLEL;
+                    }
+                } else {
+                    return EDGE_OTHER;
+                }
+            }
+            f = f->more;
+        }
+        // Pick arbitrarily which side to send it down, doesn't matter
+        return neg ? neg->ClassifyPoint(p, eb) : OUTSIDE;
+    } else if(dp > 0) {
+        if(I == 2 && N == 5) dbp("   pos");
+        return pos ? pos->ClassifyPoint(p, eb) : INSIDE;
+    } else {
+        if(I == 2 && N == 5) dbp("   neg");
+        return neg ? neg->ClassifyPoint(p, eb) : OUTSIDE;
+    }
+}
+
+int SBspUv::ClassifyEdge(Point2d ea, Point2d eb) {
+    return ClassifyPoint((ea.Plus(eb)).ScaledBy(0.5), eb);
+}
+

srf/ratpoly.cpp

@@ -522,7 +522,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/100)) {
+        if(p0.Equals(p, LENGTH_EPS/1e3)) {
             return;
         }
 

srf/surface.h

@@ -6,6 +6,35 @@
 double Bernstein(int k, int deg, double t);
 double BernsteinDerivative(int k, int deg, double t);
 
+// Utility data structure, a two-dimensional BSP to accelerate polygon
+// operations.
+class SBspUv {
+public:
+    Point2d  a, b;
+
+    SBspUv  *pos;
+    SBspUv  *neg;
+
+    SBspUv  *more;
+
+    static const int INSIDE            = 100;
+    static const int OUTSIDE           = 200;
+    static const int EDGE_PARALLEL     = 300;
+    static const int EDGE_ANTIPARALLEL = 400;
+    static const int EDGE_OTHER        = 500;
+
+    static SBspUv *Alloc(void);
+    static SBspUv *From(SEdgeList *el);
+
+    Point2d IntersectionWith(Point2d a, Point2d b);
+    SBspUv *InsertEdge(Point2d a, Point2d b);
+    int ClassifyPoint(Point2d p, Point2d eb);
+    int ClassifyEdge(Point2d ea, Point2d eb);
+};
+
+// Now the data structures to represent a shell of trimmed rational polynomial
+// surfaces.
+
 class SShell;
 
 class hSSurface {
@@ -94,7 +123,7 @@ public:
     hSSurface       surfB;
 
     static SCurve FromTransformationOf(SCurve *a, Vector t, Quaternion q);
-    SCurve MakeCopySplitAgainst(SShell *against);
+    SCurve MakeCopySplitAgainst(SShell *agnstA, SShell *agnstB);
 
     void Clear(void);
 };
@@ -130,21 +159,20 @@ public:
 
     List<STrimBy>   trim;
 
-    // The trims broken down into piecewise linear segments.
+    // For testing whether a point (u, v) on the surface lies inside the trim
-    SEdgeList       orig;
+    SBspUv          *bsp;
-    SEdgeList       inside;
-    SEdgeList       onSameNormal;
-    SEdgeList       onFlipNormal;
 
     static SSurface FromExtrusionOf(SBezier *spc, Vector t0, Vector t1);
     static SSurface FromPlane(Vector pt, Vector u, Vector v);
     static SSurface FromTransformationOf(SSurface *a, Vector t, Quaternion q, 
                                          bool includingTrims);
 
-    SSurface MakeCopyTrimAgainst(SShell *against, SShell *shell,
+    SSurface MakeCopyTrimAgainst(SShell *against, SShell *parent, SShell *into,
                                  int type, bool opA);
     void TrimFromEdgeList(SEdgeList *el);
-    void IntersectAgainst(SSurface *b, SShell *into);
+    void IntersectAgainst(SSurface *b, SShell *agnstA, SShell *agnstB, 
+                          SShell *into);
+    void AllPointsIntersecting(Vector a, Vector b, List<Vector> *l);
 
     void ClosestPointTo(Vector p, double *u, double *v);
     Vector PointAt(double u, double v);
@@ -154,6 +182,7 @@ public:
 
     void TriangulateInto(SShell *shell, SMesh *sm);
     void MakeEdgesInto(SShell *shell, SEdgeList *sel, bool asUv);
+    void MakeClassifyingBsp(SShell *shell);
 
     void Clear(void);
 };
@@ -172,10 +201,11 @@ public:
     static const int AS_DIFFERENCE = 11;
     static const int AS_INTERSECT  = 12;
     void MakeFromBoolean(SShell *a, SShell *b, int type);
-    void CopyCurvesSplitAgainst(SShell *against, SShell *into);
+    void CopyCurvesSplitAgainst(SShell *agnstA, SShell *agnstB, SShell *into);
     void CopySurfacesTrimAgainst(SShell *against, SShell *into, int t, bool a);
     void MakeIntersectionCurvesAgainst(SShell *against, SShell *into);
-    void MakeEdgeListUseNewCurveIds(SEdgeList *el);
+    void MakeClassifyingBsps(void);
+    void AllPointsIntersecting(Vector a, Vector b, List<Vector> *il);
     void CleanupAfterBoolean(void);
 
     void MakeFromCopyOf(SShell *a);

srf/surfinter.cpp

@@ -1,6 +1,8 @@
 #include "solvespace.h"
 
-void SSurface::IntersectAgainst(SSurface *b, SShell *into) {
+void SSurface::IntersectAgainst(SSurface *b, SShell *agnstA, SShell *agnstB, 
+                                SShell *into)
+{
     Vector amax, amin, bmax, bmin;
     GetAxisAlignedBounding(&amax, &amin);
     b->GetAxisAlignedBounding(&bmax, &bmin);
@@ -37,16 +39,16 @@ void SSurface::IntersectAgainst(SSurface *b, SShell *into) {
         ZERO(&sc);
         sc.surfA = h;
         sc.surfB = b->h;
-        v = inter.Minus(d.WithMagnitude(2*maxl));
+        v = inter.Minus(d.WithMagnitude(5*maxl));
         sc.pts.Add(&v);
-        v = inter.Plus(d.WithMagnitude(2*maxl));
+        v = inter.Plus(d.WithMagnitude(5*maxl));
         sc.pts.Add(&v);
-        sc.interCurve = true;
+
-       
         // Now split the line where it intersects our existing surfaces
-        SCurve split = sc.MakeCopySplitAgainst(into);
+        SCurve split = sc.MakeCopySplitAgainst(agnstA, agnstB);
         sc.Clear();
 
+        split.interCurve = true;
         into->curve.AddAndAssignId(&split);
     }
 
@@ -54,3 +56,34 @@ void SSurface::IntersectAgainst(SSurface *b, SShell *into) {
     // cases, just giving up for now
 }
 
+void SSurface::AllPointsIntersecting(Vector a, Vector b, List<Vector> *l) {
+    if(degm == 1 && degn == 1) {
+        // line-plane intersection
+        Vector p = ctrl[0][0];
+        Vector n = NormalAt(0, 0).WithMagnitude(1);
+        double d = n.Dot(p);
+        if((n.Dot(a) - d < -LENGTH_EPS && n.Dot(b) - d > LENGTH_EPS) ||
+           (n.Dot(b) - d < -LENGTH_EPS && n.Dot(a) - d > LENGTH_EPS))
+        {
+            // It crosses the plane, one point of intersection
+            // (a + t*(b - a)) dot n = d
+            // (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 };
+            ClosestPointTo(pi, &(puv.x), &(puv.y));
+            if(bsp->ClassifyPoint(puv, dummy) != SBspUv::OUTSIDE) {
+                l->Add(&pi);
+            }
+        }
+    }
+}
+
+void SShell::AllPointsIntersecting(Vector a, Vector b, List<Vector> *il) {
+    SSurface *ss;
+    for(ss = surface.First(); ss; ss = surface.NextAfter(ss)) {
+        ss->AllPointsIntersecting(a, b, il);
+    }
+}
+

util.cpp

@@ -515,6 +515,13 @@ Point2d Vector::Project2d(Vector u, Vector v) {
     return p;
 }
 
+Point2d Vector::ProjectXy(void) {
+    Point2d p;
+    p.x = x;
+    p.y = y;
+    return p;
+}
+
 double Vector::DivPivoting(Vector delta) {
     double mx = fabs(delta.x), my = fabs(delta.y), mz = fabs(delta.z);
 
@@ -654,6 +661,10 @@ double Point2d::DistanceTo(Point2d p) {
     return sqrt(dx*dx + dy*dy);
 }
 
+double Point2d::Dot(Point2d p) {
+    return x*p.x + y*p.y;
+}
+
 double Point2d::DistanceToLine(Point2d p0, Point2d dp, bool segment) {
     double m = dp.x*dp.x + dp.y*dp.y;
     if(m < 0.05) return 1e12;
@@ -673,4 +684,10 @@ double Point2d::DistanceToLine(Point2d p0, Point2d dp, bool segment) {
     }
 }
 
+Point2d Point2d::Normal(void) {
+    Point2d ret;
+    ret.x = y;
+    ret.y = -x;
+    return ret;
+}