Add a special case for line-cylinder intersection, solving in

closed form. This is a fairly good speedup, and handles tangency
well.

But that shows that tangency has other problems; need to classify
edges correctly (whether they point to a coincident surface) in
curved surfaces too. I need to tweak SShell::ClassifyPoint().

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

dsc.h

@@ -103,6 +103,7 @@ public:
     Point2d Plus(Point2d b);
     Point2d Minus(Point2d b);
     Point2d ScaledBy(double s);
+    double DivPivoting(Point2d delta);
     double Dot(Point2d p);
     double DistanceTo(Point2d p);
     double DistanceToLine(Point2d p0, Point2d dp, bool segment);

solvespace.h

@@ -26,6 +26,12 @@ inline double WRAP_NOT_0(double v, double n) {
     while(v <= 0) v += n;
     return v;
 }
+inline double WRAP_SYMMETRIC(double v, double n) {
+    // Clamp it to the range (-n/2, n/2]
+    while(v >   n/2) v -= n;
+    while(v <= -n/2) v += n;
+    return v;
+}
 
 #define SWAP(T, a, b) do { T temp = (a); (a) = (b); (b) = temp; } while(0)
 #define ZERO(v) memset((v), 0, sizeof(*(v)))

srf/boolean.cpp

@@ -47,8 +47,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, true, true);
+        if(agnstA) agnstA->AllPointsIntersecting(prev, *p, &il, true,true,true);
-        if(agnstB) agnstB->AllPointsIntersecting(prev, *p, &il, true, true);
+        if(agnstB) agnstB->AllPointsIntersecting(prev, *p, &il, true,true,true);
 
         if(il.n > 0) {
             // The intersections were generated by intersecting the pwl
@@ -377,12 +377,12 @@ SSurface SSurface::MakeCopyTrimAgainst(SShell *agnst, SShell *parent,
         Point2d am = (auv.Plus(buv)).ScaledBy(0.5);
         Vector pt = ret.PointAt(am.x, am.y);
         // and the outer normal from the trim polygon (within the surface)
-        Vector n = ret.NormalAt(am.x, am.y);
+        Vector surf_n = ret.NormalAt(am.x, am.y);
         Vector ea = ret.PointAt(auv.x, auv.y),
                eb = ret.PointAt(buv.x, buv.y);
-        Vector ptout = n.Cross((eb.Minus(ea)));
+        Vector edge_n = surf_n.Cross((eb.Minus(ea)));
 
-        int c_shell = agnst->ClassifyPoint(pt, ptout);
+        int c_shell = agnst->ClassifyPoint(pt, edge_n, surf_n);
 
         int tag = 0;
         tag |= INSIDE(ORIG, INDIR);
@@ -448,7 +448,7 @@ SSurface SSurface::MakeCopyTrimAgainst(SShell *agnst, SShell *parent,
     ZERO(&poly);
     final.l.ClearTags();
     if(!final.AssemblePolygon(&poly, NULL, true)) {
-        DEBUGEDGELIST(&inter, &ret);
+        DEBUGEDGELIST(&final, &ret);
     }
     poly.Clear();
 

srf/ratpoly.cpp

@@ -479,7 +479,10 @@ bool SSurface::IsCylinder(Vector *center, Vector *axis, double *r,
 
     double rd0 = center->Minus(sb.ctrl[0]).Magnitude(),
            rd2 = center->Minus(sb.ctrl[2]).Magnitude();
-    if(fabs(rd0 - rd2) > LENGTH_EPS) return false;
+    if(fabs(rd0 - rd2) > LENGTH_EPS) {
+        return false;
+    }
+    *r = rd0;
 
     Vector u = r0.WithMagnitude(1),
            v = (axis->Cross(u)).WithMagnitude(1);
@@ -490,11 +493,18 @@ bool SSurface::IsCylinder(Vector *center, Vector *axis, double *r,
     double thetaa = atan2(pa2.y, pa2.x), // in fact always zero due to csys
            thetab = atan2(pb2.y, pb2.x),
            dtheta = WRAP_NOT_0(thetab - thetaa, 2*PI);
+    if(dtheta > PI) {
+        // Not possible with a second order Bezier arc; so we must have
+        // the points backwards.
+        dtheta = 2*PI - dtheta;
+    }
 
-    if(fabs(sb.weight[1] - cos(dtheta/2)) > LENGTH_EPS) return false;
+    if(fabs(sb.weight[1] - cos(dtheta/2)) > LENGTH_EPS) {
+        return false;
+    }
 
-    *start  = (sb.ctrl[0]).Minus(*center);
+    *start  = sb.ctrl[0];
-    *finish = (sb.ctrl[2]).Minus(*center);
+    *finish = sb.ctrl[2];
 
     return true;
 }
@@ -754,6 +764,19 @@ void SSurface::GetAxisAlignedBounding(Vector *ptMax, Vector *ptMin) {
     }
 }
 
+bool SSurface::LineEntirelyOutsideBbox(Vector a, Vector b, bool segment) {
+    Vector amax, amin;
+    GetAxisAlignedBounding(&amax, &amin);
+    if(!Vector::BoundingBoxIntersectsLine(amax, amin, a, b, segment)) {
+        // The line segment could fail to intersect the bbox, but lie entirely
+        // within it and intersect the surface.
+        if(a.OutsideAndNotOn(amax, amin) && b.OutsideAndNotOn(amax, amin)) {
+            return true;
+        }
+    }
+    return false;
+}
+
 void SSurface::MakeEdgesInto(SShell *shell, SEdgeList *sel, bool asUv,
                              SShell *useCurvesFrom) {
     STrimBy *stb;

srf/surface.h

@@ -198,6 +198,7 @@ public:
                           SShell *into);
     void AddExactIntersectionCurve(SBezier *sb, SSurface *srfB,
                           SShell *agnstA, SShell *agnstB, SShell *into);
+
     typedef struct {
         int     tag;
         Point2d p;
@@ -210,7 +211,8 @@ public:
     double DepartureFromCoplanar(void);
     void SplitInHalf(bool byU, SSurface *sa, SSurface *sb);
     void AllPointsIntersecting(Vector a, Vector b,
-                                    List<SInter> *l, bool seg, bool trimmed);
+                                    List<SInter> *l,
+                                    bool seg, bool trimmed, bool inclTangent);
     void AllPointsIntersectingUntrimmed(Vector a, Vector b,
                                             int *cnt, int *level,
                                             List<Inter> *l, bool segment,
@@ -222,6 +224,7 @@ public:
     Vector PointAt(double u, double v);
     void TangentsAt(double u, double v, Vector *tu, Vector *tv);
     Vector NormalAt(double u, double v);
+    bool LineEntirelyOutsideBbox(Vector a, Vector b, bool segment);
     void GetAxisAlignedBounding(Vector *ptMax, Vector *ptMin);
     bool CoincidentWithPlane(Vector n, double d);
     bool CoincidentWith(SSurface *ss, bool sameNormal);
@@ -258,20 +261,20 @@ public:
     void MakeIntersectionCurvesAgainst(SShell *against, SShell *into);
     void MakeClassifyingBsps(void);
     void AllPointsIntersecting(Vector a, Vector b, List<SInter> *il,
-                                bool seg, bool trimmed);
+                                bool seg, bool trimmed, bool inclTangent);
     void MakeCoincidentEdgesInto(SSurface *proto, bool sameNormal,
                                  SEdgeList *el, SShell *useCurvesFrom);
     void CleanupAfterBoolean(void);
 
-    static const int INSIDE            = 100;
+    static const int INSIDE                 = 100;
-    static const int OUTSIDE           = 200;
+    static const int OUTSIDE                = 200;
-    static const int SURF_PARALLEL     = 300;
+    static const int SURF_PARALLEL          = 300;
-    static const int SURF_ANTIPARALLEL = 400;
+    static const int SURF_ANTIPARALLEL      = 400;
-    static const int EDGE_PARALLEL     = 500;
+    static const int EDGE_PARALLEL          = 500;
-    static const int EDGE_ANTIPARALLEL = 600;
+    static const int EDGE_ANTIPARALLEL      = 600;
-    static const int EDGE_TANGENT      = 700;
+    static const int EDGE_TANGENT           = 700;
 
-    int ClassifyPoint(Vector p, Vector out);
+    int ClassifyPoint(Vector p, Vector edge_n, Vector surf_n);
 
 
     void MakeFromCopyOf(SShell *a);

srf/surfinter.cpp

@@ -175,7 +175,8 @@ void SSurface::IntersectAgainst(SSurface *b, SShell *agnstA, SShell *agnstB,
 
             List<SInter> inters;
             ZERO(&inters);
-            sext->AllPointsIntersecting(p0, p0.Plus(dp), &inters, false, false);
+            sext->AllPointsIntersecting(
+                p0, p0.Plus(dp), &inters, false, false, true);
     
             SInter *si;
             for(si = inters.First(); si; si = inters.NextAfter(si)) {
@@ -387,15 +388,7 @@ void SSurface::AllPointsIntersectingUntrimmed(Vector a, Vector b,
 {
     // Test if the line intersects our axis-aligned bounding box; if no, then
     // no possibility of an intersection
-    Vector amax, amin;
+    if(LineEntirelyOutsideBbox(a, b, segment)) return;
-    GetAxisAlignedBounding(&amax, &amin);
-    if(!Vector::BoundingBoxIntersectsLine(amax, amin, a, b, segment)) {
-        // The line segment could fail to intersect the bbox, but lie entirely
-        // within it and intersect the surface.
-        if(a.OutsideAndNotOn(amax, amin) && b.OutsideAndNotOn(amax, amin)) {
-            return;
-        }
-    } 
  
     if(*cnt > 2000) {
         dbp("!!! too many subdivisions (level=%d)!", *level);
@@ -407,7 +400,10 @@ void SSurface::AllPointsIntersectingUntrimmed(Vector a, Vector b,
     // If we might intersect, and the surface is small, then switch to Newton
     // iterations.
     if(DepartureFromCoplanar() < 0.2*SS.ChordTolMm()) {
-        Vector p = (amax.Plus(amin)).ScaledBy(0.5);
+        Vector p = (ctrl[0   ][0   ]).Plus(
+                    ctrl[0   ][degn]).Plus(
+                    ctrl[degm][0   ]).Plus(
+                    ctrl[degm][degn]).ScaledBy(0.25);
         Inter inter;
         sorig->ClosestPointTo(p, &(inter.p.x), &(inter.p.y), false);
         if(sorig->PointIntersectingLine(a, b, &(inter.p.x), &(inter.p.y))) {
@@ -438,11 +434,16 @@ void SSurface::AllPointsIntersectingUntrimmed(Vector a, Vector b,
 // Find all points where a line through a and b intersects our surface, and
 // add them to the list. If seg is true then report only intersections that
 // lie within the finite line segment (not including the endpoints); otherwise
-// we work along the infinite line.
+// we work along the infinite line. And we report either just intersections
+// inside the trim curve, or any intersection with u, v in [0, 1]. And we
+// either disregard or report tangent points.
 //-----------------------------------------------------------------------------
 void SSurface::AllPointsIntersecting(Vector a, Vector b,
-                                     List<SInter> *l, bool seg, bool trimmed)
+                                     List<SInter> *l,
+                                     bool seg, bool trimmed, bool inclTangent)
 {
+    if(LineEntirelyOutsideBbox(a, b, seg)) return;
+
     Vector ba = b.Minus(a);
     double bam = ba.Magnitude();
 
@@ -468,8 +469,62 @@ void SSurface::AllPointsIntersecting(Vector a, Vector b,
             ClosestPointTo(p, &(inter.p.x), &(inter.p.y));
             inters.Add(&inter);
         }
-    } else if(IsCylinder(&center, &axis, &radius, &start, &finish) && 0) {
+    } else if(IsCylinder(&center, &axis, &radius, &start, &finish)) {
-        // XXX, cylinder is easy in closed form
+        // This one can be solved in closed form too.
+        Vector ab = b.Minus(a);
+        if(axis.Cross(ab).Magnitude() < LENGTH_EPS) {
+            // edge is parallel to axis of cylinder, no intersection points
+            return;
+        }
+        // A coordinate system centered at the center of the circle, with
+        // the edge under test horizontal
+        Vector u, v, n = axis.WithMagnitude(1);
+        u = (ab.Minus(n.ScaledBy(ab.Dot(n)))).WithMagnitude(1);
+        v = n.Cross(u);
+        Point2d ap = (a.Minus(center)).DotInToCsys(u, v, n).ProjectXy(),
+                bp = (b.Minus(center)).DotInToCsys(u, v, n).ProjectXy(),
+                sp = (start. Minus(center)).DotInToCsys(u, v, n).ProjectXy(),
+                fp = (finish.Minus(center)).DotInToCsys(u, v, n).ProjectXy();
+
+        double thetas = atan2(sp.y, sp.x), thetaf = atan2(fp.y, fp.x);
+
+        Point2d ip[2];
+        int ip_n = 0;
+        if(fabs(fabs(ap.y) - radius) < LENGTH_EPS) {
+            // tangent
+            if(inclTangent) {
+                ip[0] = Point2d::From(0, ap.y);
+                ip_n = 1;
+            }
+        } else if(fabs(ap.y) < radius) {
+            // two intersections
+            double xint = sqrt(radius*radius - ap.y*ap.y);
+            ip[0] = Point2d::From(-xint, ap.y);
+            ip[1] = Point2d::From( xint, ap.y);
+            ip_n = 2;
+        }
+        int i;
+        for(i = 0; i < ip_n; i++) {
+            double t = (ip[i].Minus(ap)).DivPivoting(bp.Minus(ap));
+            // This is a point on the circle; but is it on the arc?
+            Point2d pp = ap.Plus((bp.Minus(ap)).ScaledBy(t));
+            double theta = atan2(pp.y, pp.x);
+            double dp = WRAP_SYMMETRIC(theta  - thetas, 2*PI),
+                   df = WRAP_SYMMETRIC(thetaf - thetas, 2*PI);
+            double tol = LENGTH_EPS/radius;
+
+            if((df > 0 && ((dp < -tol) || (dp > df + tol))) ||
+               (df < 0 && ((dp >  tol) || (dp < df - tol))))
+            {
+                continue;
+            }
+
+            Vector p = a.Plus((b.Minus(a)).ScaledBy(t));
+
+            Inter inter;
+            ClosestPointTo(p, &(inter.p.x), &(inter.p.y));
+            inters.Add(&inter);
+        }
     } else {
         // General numerical solution by subdivision, fallback
         int cnt = 0, level = 0;
@@ -519,24 +574,25 @@ void SSurface::AllPointsIntersecting(Vector a, Vector b,
 }
 
 void SShell::AllPointsIntersecting(Vector a, Vector b,
-                                   List<SInter> *il, bool seg, bool trimmed)
+                                   List<SInter> *il,
+                                   bool seg, bool trimmed, bool inclTangent)
 {
     SSurface *ss;
     for(ss = surface.First(); ss; ss = surface.NextAfter(ss)) {
-        ss->AllPointsIntersecting(a, b, il, seg, trimmed);
+        ss->AllPointsIntersecting(a, b, il, seg, trimmed, inclTangent);
     }
 }
 
 //-----------------------------------------------------------------------------
-// Does the given point lie on our shell? It might be inside or outside, or
+// Does the given point lie on our shell? There are many cases; inside and
-// it might be on the surface with pout parallel or anti-parallel to the
+// outside are obvious, but then there's all the edge-on-edge and edge-on-face
-// intersecting surface's normal.
+// possibilities.
 //
 // 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) {
+int SShell::ClassifyPoint(Vector p, Vector edge_n, Vector surf_n) {
     List<SInter> l;
     ZERO(&l);
 
@@ -550,7 +606,8 @@ 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));
-        AllPointsIntersecting(p.Minus(ray), p.Plus(ray), &l, false, true);
+        AllPointsIntersecting(
+            p.Minus(ray), p.Plus(ray), &l, false, true, false);
         
         double dmin = VERY_POSITIVE;
         ret = OUTSIDE; // no intersections means it's outside
@@ -569,7 +626,7 @@ int SShell::ClassifyPoint(Vector p, Vector pout) {
 
             // Handle edge-on-edge
             if(d < LENGTH_EPS && si->onEdge && edge_inters < 2) {
-                edge_dotp[edge_inters] = (si->surfNormal).Dot(pout);
+                edge_dotp[edge_inters] = (si->surfNormal).Dot(edge_n);
                 edge_inters++;
             }
 
@@ -577,7 +634,10 @@ int SShell::ClassifyPoint(Vector p, Vector pout) {
                 dmin = d;
                 if(d < LENGTH_EPS) {
                     // Edge-on-face (unless edge-on-edge above supercedes)
-                    if((si->surfNormal).Dot(pout) > 0) {
+                    double dot = (si->surfNormal).Dot(edge_n);
+                    if(fabs(dot) < LENGTH_EPS && 0) {
+                        // TODO revamp this
+                    } else if(dot > 0) {
                         ret = SURF_PARALLEL;
                     } else {
                         ret = SURF_ANTIPARALLEL;

util.cpp

@@ -758,6 +758,14 @@ Point2d Point2d::ScaledBy(double s) {
     return r;
 }
 
+double Point2d::DivPivoting(Point2d delta) {
+    if(fabs(delta.x) > fabs(delta.y)) {
+        return x/delta.x;
+    } else {
+        return y/delta.y;
+    }
+}
+
 double Point2d::MagSquared(void) {
     return x*x + y*y;
 }