I had been using LENGTH_EPS as the tolerance on both xyz points and

uv points. This is inconsistent, unless the surface happens to be a
plane square with side length one.

So modify the SBspUv tests to take a surface, and measure distance
linearized in that surface. That fixes at least one
mis-classification bug, and doesn't seem to break anything.

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

srf/boolean.cpp

@@ -257,7 +257,7 @@ void DEBUGEDGELIST(SEdgeList *sel, SSurface *surf) {
         Vector arrow = (se->b).Minus(se->a);
         SWAP(double, arrow.x, arrow.y);
         arrow.x *= -1;
-        arrow = arrow.WithMagnitude(0.03);
+        arrow = arrow.WithMagnitude(0.01);
         arrow = arrow.Plus(mid);
 
         SS.nakedEdges.AddEdge(surf->PointAt(se->a.x, se->a.y),
@@ -426,7 +426,7 @@ SSurface SSurface::MakeCopyTrimAgainst(SShell *parent,
     ret.MakeEdgesInto(into, &orig, true);
     ret.trim.Clear();
     // which means that we can't necessarily use the old BSP...
-    SBspUv *origBsp = SBspUv::From(&orig);
+    SBspUv *origBsp = SBspUv::From(&orig, &ret);
 
     // And now intersect the other shell against us
     SEdgeList inter;
@@ -452,7 +452,7 @@ SSurface SSurface::MakeCopyTrimAgainst(SShell *parent,
                 ss->ClosestPointTo(a, &(auv.x), &(auv.y));
                 ss->ClosestPointTo(b, &(buv.x), &(buv.y));
 
-                int c = ss->bsp->ClassifyEdge(auv, buv);
+                int c = ss->bsp->ClassifyEdge(auv, buv, ss);
                 if(c != SBspUv::OUTSIDE) {
                     Vector ta = Vector::From(0, 0, 0);
                     Vector tb = Vector::From(0, 0, 0);
@@ -560,7 +560,7 @@ SSurface SSurface::MakeCopyTrimAgainst(SShell *parent,
 
         int indir_shell, outdir_shell, indir_orig, outdir_orig;
 
-        int c_this = origBsp->ClassifyEdge(auv, buv);
+        int c_this = origBsp->ClassifyEdge(auv, buv, &ret);
         TagByClassifiedEdge(c_this, &indir_orig, &outdir_orig);
 
         agnst->ClassifyEdge(&indir_shell, &outdir_shell,
@@ -764,7 +764,7 @@ void SSurface::MakeClassifyingBsp(SShell *shell, SShell *useCurvesFrom) {
     ZERO(&el);
 
     MakeEdgesInto(shell, &el, true, useCurvesFrom);
-    bsp = SBspUv::From(&el);
+    bsp = SBspUv::From(&el, this);
     el.Clear();
 
     ZERO(&edges);
@@ -787,7 +787,7 @@ static int ByLength(const void *av, const void *bv)
     // stability for the normals.
     return (la < lb) ? 1 : -1;
 }
-SBspUv *SBspUv::From(SEdgeList *el) {
+SBspUv *SBspUv::From(SEdgeList *el, SSurface *srf) {
     SEdgeList work;
     ZERO(&work);
 
@@ -799,14 +799,50 @@ SBspUv *SBspUv::From(SEdgeList *el) {
 
     SBspUv *bsp = NULL;
     for(se = work.l.First(); se; se = work.l.NextAfter(se)) {
-        bsp = bsp->InsertEdge((se->a).ProjectXy(), (se->b).ProjectXy());
+        bsp = bsp->InsertEdge((se->a).ProjectXy(), (se->b).ProjectXy(), srf);
     }
 
     work.Clear();
     return bsp;
 }
 
-SBspUv *SBspUv::InsertEdge(Point2d ea, Point2d eb) {
+//-----------------------------------------------------------------------------
+// The points in this BSP are in uv space, but we want to apply our tolerances
+// consistently in xyz (i.e., we want to say a point is on-edge if its xyz
+// distance to that edge is less than LENGTH_EPS, irrespective of its distance
+// in uv). So we linearize the surface about the point we're considering and
+// then do the test. That preserves point-on-line relationships, and the only
+// time we care about exact correctness is when we're very close to the line,
+// which is when the linearization is accurate.
+//-----------------------------------------------------------------------------
+void SBspUv::ScalePoints(Point2d *pt, Point2d *a, Point2d *b, SSurface *srf) {
+    Vector tu, tv;
+    srf->TangentsAt(pt->x, pt->y, &tu, &tv);
+    double mu = tu.Magnitude(), mv = tv.Magnitude();
+
+    pt->x *= mu; pt->y *= mv;
+    a ->x *= mu; a ->y *= mv;
+    b ->x *= mu; b ->y *= mv;
+}
+double SBspUv::ScaledSignedDistanceToLine(Point2d pt, Point2d a, Point2d b, 
+                                            SSurface *srf)
+{
+    ScalePoints(&pt, &a, &b, srf);
+
+    Point2d n = ((b.Minus(a)).Normal()).WithMagnitude(1);
+    double d = a.Dot(n);
+
+    return pt.Dot(n) - d;
+}
+double SBspUv::ScaledDistanceToLine(Point2d pt, Point2d a, Point2d b, bool seg,
+                                        SSurface *srf)
+{
+    ScalePoints(&pt, &a, &b, srf);
+
+    return pt.DistanceToLine(a, b, seg);
+}
+
+SBspUv *SBspUv::InsertEdge(Point2d ea, Point2d eb, SSurface *srf) {
     if(!this) {
         SBspUv *ret = Alloc();
         ret->a = ea;
@@ -814,11 +850,8 @@ SBspUv *SBspUv::InsertEdge(Point2d ea, Point2d eb) {
         return ret;
     }
 
-    Point2d n = ((b.Minus(a)).Normal()).WithMagnitude(1);
+    double dea = ScaledSignedDistanceToLine(ea, a, b, srf),
-    double d = a.Dot(n);
+           deb = ScaledSignedDistanceToLine(eb, a, b, srf);
-
-    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
@@ -830,55 +863,49 @@ SBspUv *SBspUv::InsertEdge(Point2d ea, Point2d eb) {
     } 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);
+            pos = pos->InsertEdge(ea, eb, srf);
         } else {
-            neg = neg->InsertEdge(ea, eb);
+            neg = neg->InsertEdge(ea, eb, srf);
         }
     } 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);
+            pos = pos->InsertEdge(ea, eb, srf);
         } else {
-            neg = neg->InsertEdge(ea, eb);
+            neg = neg->InsertEdge(ea, eb, srf);
         }
     } else if(dea > 0 && deb > 0) {
-        pos = pos->InsertEdge(ea, eb);
+        pos = pos->InsertEdge(ea, eb, srf);
     } else if(dea < 0 && deb < 0) {
-        neg = neg->InsertEdge(ea, eb);
+        neg = neg->InsertEdge(ea, eb, srf);
     } else {
         // New edge crosses this one; we need to split.
+        Point2d n = ((b.Minus(a)).Normal()).WithMagnitude(1);
+        double d = a.Dot(n);
         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);
+            pos = pos->InsertEdge(ea, pi, srf);
-            neg = neg->InsertEdge(pi, eb);
+            neg = neg->InsertEdge(pi, eb, srf);
         } else {
-            neg = neg->InsertEdge(ea, pi);
+            neg = neg->InsertEdge(ea, pi, srf);
-            pos = pos->InsertEdge(pi, eb);
+            pos = pos->InsertEdge(pi, eb, srf);
         }
     }
     return this;
 }
 
-int SBspUv::ClassifyPoint(Point2d p, Point2d eb, Point2d *ia, Point2d *ib) {
+int SBspUv::ClassifyPoint(Point2d p, Point2d eb, SSurface *srf) {
     if(!this) return OUTSIDE;
 
-    Point2d n = ((b.Minus(a)).Normal()).WithMagnitude(1);
+    double dp = ScaledSignedDistanceToLine(p, a, b, srf);
-    double d = a.Dot(n);
-
-    double dp = p.Dot(n) - d;
 
     if(fabs(dp) < LENGTH_EPS) {
         SBspUv *f = this;
         while(f) {
             Point2d ba = (f->b).Minus(f->a);
-            if(p.DistanceToLine(f->a, ba, true) < LENGTH_EPS) {
+            if(ScaledDistanceToLine(p, f->a, ba, true, srf) < LENGTH_EPS) {
-                if(ia && ib) {
+                if(ScaledDistanceToLine(eb, f->a, ba, false, srf) < LENGTH_EPS){
-                    // Tell the caller which edge it happens to lie on.
-                    *ia = f->a;
-                    *ib = f->b;
-                }
-                if(eb.DistanceToLine(f->a, ba, false) < LENGTH_EPS) {
                     if(ba.Dot(eb.Minus(p)) > 0) {
                         return EDGE_PARALLEL;
                     } else {
@@ -891,26 +918,26 @@ int SBspUv::ClassifyPoint(Point2d p, Point2d eb, Point2d *ia, Point2d *ib) {
             f = f->more;
         }
         // Pick arbitrarily which side to send it down, doesn't matter
-        int c1 =  neg ? neg->ClassifyPoint(p, eb, ia, ib) : OUTSIDE;
+        int c1 =  neg ? neg->ClassifyPoint(p, eb, srf) : OUTSIDE;
-        int c2 =  pos ? pos->ClassifyPoint(p, eb, ia, ib) : INSIDE;
+        int c2 =  pos ? pos->ClassifyPoint(p, eb, srf) : INSIDE;
         if(c1 != c2) {
             dbp("MISMATCH: %d %d %08x %08x", c1, c2, neg, pos);
         }
         return c1;
     } else if(dp > 0) {
-        return pos ? pos->ClassifyPoint(p, eb, ia, ib) : INSIDE;
+        return pos ? pos->ClassifyPoint(p, eb, srf) : INSIDE;
     } else {
-        return neg ? neg->ClassifyPoint(p, eb, ia, ib) : OUTSIDE;
+        return neg ? neg->ClassifyPoint(p, eb, srf) : OUTSIDE;
     }
 }
 
-int SBspUv::ClassifyEdge(Point2d ea, Point2d eb) {
+int SBspUv::ClassifyEdge(Point2d ea, Point2d eb, SSurface *srf) {
-    int ret = ClassifyPoint((ea.Plus(eb)).ScaledBy(0.5), eb);
+    int ret = ClassifyPoint((ea.Plus(eb)).ScaledBy(0.5), eb, srf);
     if(ret == EDGE_OTHER) {
         // Perhaps the edge is tangent at its midpoint (and we screwed up
         // somewhere earlier and failed to split it); try a different
         // point on the edge.
-        ret = ClassifyPoint(ea.Plus((eb.Minus(ea)).ScaledBy(0.294)), eb);
+        ret = ClassifyPoint(ea.Plus((eb.Minus(ea)).ScaledBy(0.294)), eb, srf);
     }
     return ret;
 }

srf/raycast.cpp

@@ -359,8 +359,8 @@ void SSurface::AllPointsIntersecting(Vector a, Vector b,
         }
        
         // And that it lies inside our trim region
-        Point2d dummy = { 0, 0 }, ia = { 0, 0 }, ib = { 0, 0 };
+        Point2d dummy = { 0, 0 };
-        int c = bsp->ClassifyPoint(puv, dummy, &ia, &ib);
+        int c = bsp->ClassifyPoint(puv, dummy, this);
         if(trimmed && c == SBspUv::OUTSIDE) {
             continue;
         }
@@ -372,8 +372,6 @@ void SSurface::AllPointsIntersecting(Vector a, Vector b,
         si.pinter = puv;
         si.srf = this;
         si.onEdge = (c != SBspUv::INSIDE);
-        si.edgeA = ia;
-        si.edgeB = ib;
         l->Add(&si);
     }
 
@@ -528,8 +526,8 @@ bool SShell::ClassifyEdge(int *indir, int *outdir,
         Vector pp = srf->PointAt(puv);
 
         if((pp.Minus(p)).Magnitude() > LENGTH_EPS) continue;
-        Point2d dummy = { 0, 0 }, ia = { 0, 0 }, ib = { 0, 0 };
+        Point2d dummy = { 0, 0 };
-        int c = srf->bsp->ClassifyPoint(puv, dummy, &ia, &ib);
+        int c = srf->bsp->ClassifyPoint(puv, dummy, srf);
         if(c == SBspUv::OUTSIDE) continue;
 
         // Edge-on-face (unless edge-on-edge above superceded)

srf/surface.h

@@ -27,13 +27,17 @@ public:
     static const int EDGE_OTHER        = 500;
 
     static SBspUv *Alloc(void);
-    static SBspUv *From(SEdgeList *el);
+    static SBspUv *From(SEdgeList *el, SSurface *srf);
 
-    Point2d IntersectionWith(Point2d a, Point2d b);
+    void ScalePoints(Point2d *pt, Point2d *a, Point2d *b, SSurface *srf);
-    SBspUv *InsertEdge(Point2d a, Point2d b);
+    double ScaledSignedDistanceToLine(Point2d pt, Point2d a, Point2d b,
-    int ClassifyPoint(Point2d p, Point2d eb,
+        SSurface *srf);
-                                        Point2d *ia=NULL, Point2d *ib=NULL);
+    double ScaledDistanceToLine(Point2d pt, Point2d a, Point2d b, bool seg,
-    int ClassifyEdge(Point2d ea, Point2d eb);
+        SSurface *srf);
+
+    SBspUv *InsertEdge(Point2d a, Point2d b, SSurface *srf);
+    int ClassifyPoint(Point2d p, Point2d eb, SSurface *srf);
+    int ClassifyEdge(Point2d ea, Point2d eb, SSurface *srf);
 };
 
 // Now the data structures to represent a shell of trimmed rational polynomial
@@ -191,7 +195,6 @@ public:
     Point2d     pinter;
     Vector      surfNormal;     // of the intersecting surface, at pinter
     bool        onEdge;         // pinter is on edge of trim poly
-    Point2d     edgeA, edgeB;   // the edge that pinter is on
 };
 
 // A rational polynomial surface in Bezier form.

srf/surfinter.cpp

@@ -395,9 +395,9 @@ void SSurface::IntersectAgainst(SSurface *b, SShell *agnstA, SShell *agnstB,
                     tol = (npc.Minus(np)).Magnitude();
 
                     if(tol > maxtol*0.8) {
-                        step *= 0.95;
+                        step *= 0.90;
                     } else {
-                        step /= 0.95;
+                        step /= 0.90;
                     }
 
                     if((tol < maxtol) && (tol > maxtol/2)) {