Add code to subdivide (with de Castljau's algorithm) a surface, and

use that for surface-line intersections. That has major problems
with the heuristic on when to stop and do Newton polishing.

There's also an issue with all the Newton stuff when surfaces join
tangent.

And update the wishlist to reflect current needs.

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

dsc.h

@@ -84,6 +84,8 @@ public:
     void MakeMaxMin(Vector *maxv, Vector *minv);
     static bool BoundingBoxesDisjoint(Vector amax, Vector amin,
                                       Vector bmax, Vector bmin);
+    static bool BoundingBoxIntersectsLine(Vector amax, Vector amin,
+                                          Vector p0, Vector p1, bool segment);
     bool OutsideAndNotOn(Vector maxv, Vector minv);
     Point2d Project2d(Vector u, Vector v);
     Point2d ProjectXy(void);

solvespace.cpp

@@ -177,6 +177,9 @@ double SolveSpace::StringToMm(char *str) {
         return atof(str);
     }
 }
+double SolveSpace::ChordTolMm(void) {
+    return SS.chordTol / SS.GW.scale;
+}
 
 
 void SolveSpace::AfterNewFile(void) {

solvespace.h

@@ -31,7 +31,7 @@ inline double WRAP_NOT_0(double v, double n) {
 #define ZERO(v) memset((v), 0, sizeof(*(v)))
 #define CO(v) (v).x, (v).y, (v).z
 
-#define LENGTH_EPS      (1e-7)
+#define LENGTH_EPS      (1e-6)
 #define VERY_POSITIVE   (1e10)
 #define VERY_NEGATIVE   (-1e10)
 
@@ -448,6 +448,7 @@ public:
     char *MmToString(double v);
     double ExprToMm(Expr *e);
     double StringToMm(char *s);
+    double ChordTolMm(void);
 
     // The platform-dependent code calls this before entering the msg loop
     void Init(char *cmdLine);

srf/boolean.cpp

@@ -50,24 +50,30 @@ SCurve SCurve::MakeCopySplitAgainst(SShell *agnstA, SShell *agnstB,
         if(agnstA) agnstA->AllPointsIntersecting(prev, *p, &il, true, true);
         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.
         if(il.n > 0) {
-            LineStart = prev;
+            // The intersections were generated by intersecting the pwl
-            LineDirection = p->Minus(prev);
+            // edge against a surface; so they must be refined to lie
-            qsort(il.elem, il.n, sizeof(il.elem[0]), ByTAlongLine);
+            // exactly on the original curve.
-    
-            Vector prev = Vector::From(VERY_POSITIVE, 0, 0);
             SInter *pi;
             for(pi = il.First(); pi; pi = il.NextAfter(pi)) {
-                // The intersection was generated by intersecting the pwl
-                // edge against a surface; so it must be refined to lie
-                // exactly on the original curve.
                 double u, v;
                 (pi->srf)->ClosestPointTo(pi->p, &u, &v, false);
                 (pi->srf)->PointOnSurfaces(srfA, srfB, &u, &v);
                 pi->p = (pi->srf)->PointAt(u, v);
+            }
 
+            // And now sort them in order along the line. Note that we must
+            // do that after refining, in case the refining would make two
+            // points switch places.
+            LineStart = prev;
+            LineDirection = p->Minus(prev);
+            qsort(il.elem, il.n, sizeof(il.elem[0]), ByTAlongLine);
+   
+            // And now uses the intersections to generate our split pwl edge(s)
+            Vector prev = Vector::From(VERY_POSITIVE, 0, 0);
+            for(pi = il.First(); pi; pi = il.NextAfter(pi)) {
+                double t = (pi->p.Minus(LineStart)).DivPivoting(LineDirection);
+                dbp("t=%.3f", t);
                 // On-edge intersection will generate same split point for
                 // both surfaces, so don't create zero-length edge.
                 if(!prev.Equals(pi->p)) {
@@ -236,6 +242,7 @@ void DBPEDGE(int tag) {
 }
 
 void DEBUGEDGELIST(SEdgeList *sel, SSurface *surf) {
+    dbp("print %d edges", sel->l.n);
     SEdge *se;
     for(se = sel->l.First(); se; se = sel->l.NextAfter(se)) {
         Vector mid = (se->a).Plus(se->b).ScaledBy(0.5);
@@ -372,11 +379,7 @@ SSurface SSurface::MakeCopyTrimAgainst(SShell *agnst, SShell *parent,
                eb = ret.PointAt(buv.x, buv.y);
         Vector ptout = n.Cross((eb.Minus(ea)));
 
-        if(I == 1 && fabs(((se->b).Minus(se->a).Magnitude()) - 0.166) < 0.01) {
-            FLAG = 1;
-        }
         int c_shell = agnst->ClassifyPoint(pt, ptout);
-        FLAG = 0;
 
         int tag = 0;
         tag |= INSIDE(ORIG, INDIR);
@@ -456,7 +459,7 @@ SSurface SSurface::MakeCopyTrimAgainst(SShell *agnst, SShell *parent,
     }
     final.l.RemoveTagged();
 
-//    if(I == 1) DEBUGEDGELIST(&final, &ret);
+//    if(I == 0) DEBUGEDGELIST(&final, &ret);
 
     // Use our reassembled edges to trim the new surface.
     ret.TrimFromEdgeList(&final);
@@ -491,6 +494,7 @@ void SShell::MakeIntersectionCurvesAgainst(SShell *agnst, SShell *into) {
             // list for into.
             sa->IntersectAgainst(sb, agnst, this, into);
         }
+        FLAG++;
     }
 }
 

srf/ratpoly.cpp

@@ -4,7 +4,7 @@
 // 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.
-#define RATPOLY_EPS (LENGTH_EPS/(1e1))
+#define RATPOLY_EPS (LENGTH_EPS/(1e2))
 
 static double Bernstein(int k, int deg, double t)
 {
@@ -131,10 +131,8 @@ void SBezier::MakePwlWorker(List<Vector> *l, double ta, double tb, Vector off) {
     double d = max(pm1.DistanceToLine(pa, pb.Minus(pa)),
                    pm2.DistanceToLine(pa, pb.Minus(pa)));
 
-    double tol = SS.chordTol/SS.GW.scale;
-
     double step = 1.0/SS.maxSegments;
-    if((tb - ta) < step || d < tol) {
+    if((tb - ta) < step || d < SS.ChordTolMm()) {
         // A previous call has already added the beginning of our interval.
         pb = pb.Plus(off);
         l->Add(&pb);
@@ -536,7 +534,7 @@ void SSurface::ClosestPointTo(Vector p, double *u, double *v, bool converge) {
         *u = *v = 0; // a plane, perfect no matter what the initial guess
     } else {
         double minDist = VERY_POSITIVE;
-        double res = 7.0;
+        double res = (max(degm, degn) == 2) ? 7.0 : 20.0;
         for(i = 0; i < (int)res; i++) {
             for(j = 0; j <= (int)res; j++) {
                 double tryu = (i/res), tryv = (j/res);

srf/surface.h

@@ -194,8 +194,22 @@ public:
                           SShell *into);
     void AddExactIntersectionCurve(SBezier *sb, SSurface *srfB,
                           SShell *agnstA, SShell *agnstB, SShell *into);
+    typedef struct {
+        int     tag;
+        Point2d p;
+    } Inter;
+    void WeightControlPoints(void);
+    void UnWeightControlPoints(void);
+    void CopyRowOrCol(bool row, int this_ij, SSurface *src, int src_ij);
+    void BlendRowOrCol(bool row, int this_ij, SSurface *a, int a_ij,
+                                              SSurface *b, int b_ij);
+    void SplitInHalf(bool byU, SSurface *sa, SSurface *sb);
     void AllPointsIntersecting(Vector a, Vector b,
                                     List<SInter> *l, bool seg, bool trimmed);
+    void AllPointsIntersectingUntrimmed(Vector a, Vector b,
+                                            int *cnt, int *level,
+                                            List<Inter> *l, bool segment,
+                                            SSurface *sorig);
 
     void ClosestPointTo(Vector p, double *u, double *v, bool converge=true);
     bool PointIntersectingLine(Vector p0, Vector p1, double *u, double *v);

srf/surfinter.cpp

@@ -11,16 +11,22 @@ void SSurface::AddExactIntersectionCurve(SBezier *sb, SSurface *srfB,
     sc.surfB = srfB->h;
     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, this, srfB);
     sc.Clear();
+   
+/*
+    if(sb->deg == 1) {
+        dbp(" ");
+        Vector *prev = NULL, *v;
+        dbp("split.pts.n =%d", split.pts.n);
+        for(v = split.pts.First(); v; v = split.pts.NextAfter(v)) {
+            if(prev) {
+                SS.nakedEdges.AddEdge(*prev, *v);
+            }
+            prev = v;
+        }
+    } */
 
     split.source = SCurve::FROM_INTERSECTION;
     into->curve.AddAndAssignId(&split);
@@ -88,7 +94,10 @@ void SSurface::IntersectAgainst(SSurface *b, SShell *agnstA, SShell *agnstB,
             inters.Clear();
         } else {
             // Direction of extrusion is not parallel to plane; so
-            // intersection is projection of extruded curve into our plane
+            // intersection is projection of extruded curve into our plane.
+            // If both curves are planes, then we could do it either way;
+            // so choose the one that generates the shorter curve.
+            // XXX TODO
 
             int i;
             for(i = 0; i <= bezier.deg; i++) {
@@ -107,6 +116,177 @@ void SSurface::IntersectAgainst(SSurface *b, SShell *agnstA, SShell *agnstB,
     // cases, just giving up for now
 }
 
+void SSurface::WeightControlPoints(void) {
+    int i, j;
+    for(i = 0; i <= degm; i++) {
+        for(j = 0; j <= degn; j++) {
+            ctrl[i][j] = (ctrl[i][j]).ScaledBy(weight[i][j]);
+        }
+    }
+}
+void SSurface::UnWeightControlPoints(void) {
+    int i, j;
+    for(i = 0; i <= degm; i++) {
+        for(j = 0; j <= degn; j++) {
+            ctrl[i][j] = (ctrl[i][j]).ScaledBy(1.0/weight[i][j]);
+        }
+    }
+}
+void SSurface::CopyRowOrCol(bool row, int this_ij, SSurface *src, int src_ij) {
+    if(row) {
+        int j;
+        for(j = 0; j <= degn; j++) {
+            ctrl  [this_ij][j] = src->ctrl  [src_ij][j];
+            weight[this_ij][j] = src->weight[src_ij][j];
+        }
+    } else {
+        int i;
+        for(i = 0; i <= degm; i++) {
+            ctrl  [i][this_ij] = src->ctrl  [i][src_ij];
+            weight[i][this_ij] = src->weight[i][src_ij];
+        }
+    }
+}
+void SSurface::BlendRowOrCol(bool row, int this_ij, SSurface *a, int a_ij,
+                                                    SSurface *b, int b_ij)
+{
+    if(row) {
+        int j;
+        for(j = 0; j <= degn; j++) {
+            Vector c = (a->ctrl  [a_ij][j]).Plus(b->ctrl  [b_ij][j]);
+            double w = (a->weight[a_ij][j]   +   b->weight[b_ij][j]);
+            ctrl  [this_ij][j] = c.ScaledBy(0.5);
+            weight[this_ij][j] = w / 2;
+        }
+    } else {
+        int i;
+        for(i = 0; i <= degm; i++) {
+            Vector c = (a->ctrl  [i][a_ij]).Plus(b->ctrl  [i][b_ij]);
+            double w = (a->weight[i][a_ij]   +   b->weight[i][b_ij]);
+            ctrl  [i][this_ij] = c.ScaledBy(0.5);
+            weight[i][this_ij] = w / 2;
+        }
+    }
+}
+void SSurface::SplitInHalf(bool byU, SSurface *sa, SSurface *sb) {
+    sa->degm = sb->degm = degm;
+    sa->degn = sb->degn = degn;
+
+    // by de Casteljau's algorithm in a projective space; so we must work
+    // on points (w*x, w*y, w*z, w)
+    WeightControlPoints();
+
+    switch(byU ? degm : degn) {
+        case 1:
+            sa->CopyRowOrCol (byU, 0, this, 0);
+            sb->CopyRowOrCol (byU, 1, this, 1);
+
+            sa->BlendRowOrCol(byU, 1, this, 0, this, 1);
+            sb->BlendRowOrCol(byU, 0, this, 0, this, 1);
+            break;
+
+        case 2:
+            sa->CopyRowOrCol (byU, 0, this, 0);
+            sb->CopyRowOrCol (byU, 2, this, 2);
+
+            sa->BlendRowOrCol(byU, 1, this, 0, this, 1);
+            sb->BlendRowOrCol(byU, 1, this, 1, this, 2);
+
+            sa->BlendRowOrCol(byU, 2, sa,   1, sb,   1);
+            sb->BlendRowOrCol(byU, 0, sa,   1, sb,   1);
+            break;
+
+        case 3: {
+            SSurface st;
+            st.degm = degm; st.degn = degn;
+
+            sa->CopyRowOrCol (byU, 0, this, 0);
+            sb->CopyRowOrCol (byU, 3, this, 3);
+
+            sa->BlendRowOrCol(byU, 1, this, 0, this, 1);
+            sb->BlendRowOrCol(byU, 2, this, 2, this, 3);
+            st. BlendRowOrCol(byU, 0, this, 1, this, 2); // scratch var
+
+            sa->BlendRowOrCol(byU, 2, sa,   1, &st,  0);
+            sb->BlendRowOrCol(byU, 1, sb,   2, &st,  0);
+
+            sa->BlendRowOrCol(byU, 3, sa,   2, sb,   1);
+            sb->BlendRowOrCol(byU, 0, sa,   2, sb,   1);
+            break;
+        }
+
+        default: oops();
+    }
+
+    sa->UnWeightControlPoints();
+    sb->UnWeightControlPoints();
+    UnWeightControlPoints();
+}
+
+//-----------------------------------------------------------------------------
+// Find all points where the indicated finite (if segment) or infinite (if not
+// segment) line intersects our surface. Report them in uv space in the list.
+// We first do a bounding box check; if the line doesn't intersect, then we're
+// done. If it does, then we check how small our surface is. If it's big,
+// then we subdivide into quarters and recurse. If it's small, then we refine
+// by Newton's method and record the point.
+//-----------------------------------------------------------------------------
+void SSurface::AllPointsIntersectingUntrimmed(Vector a, Vector b,
+                                              int *cnt, int *level,
+                                              List<Inter> *l, bool segment,
+                                              SSurface *sorig)
+{
+    // Test if the line intersects our axis-aligned bounding box; if no, then
+    // no possibility of an intersection
+    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;
+        }
+    } 
+   
+    if(*cnt > 2000) {
+        dbp("!!! too many subdivisions (level=%d)!", *level);
+        return;
+    }
+    (*cnt)++;
+
+    // If we might intersect, and the surface is small, then switch to Newton
+    // iterations.
+    double h = max(amax.x - amin.x,
+               max(amax.y - amin.y,
+                   amax.z - amin.z));
+    if(fabs(h) < SS.ChordTolMm()) {
+        Vector p = (amax.Plus(amin)).ScaledBy(0.5);
+        Inter inter;
+        sorig->ClosestPointTo(p, &(inter.p.x), &(inter.p.y), false);
+        if(sorig->PointIntersectingLine(a, b, &(inter.p.x), &(inter.p.y))) {
+            Vector p = sorig->PointAt(inter.p.x, inter.p.y);
+            // Debug check, verify that the point lies in both surfaces
+            // (which it ought to, since the surfaces should be coincident)
+            double u, v;
+            ClosestPointTo(p, &u, &v);
+            l->Add(&inter);
+        } else {
+            // Might not converge if line is almost tangent to surface...
+        }
+        return;
+    }
+
+    // But the surface is big, so split it, alternating by u and v
+    SSurface surf0, surf1;
+    SplitInHalf((*level & 1) == 0, &surf0, &surf1);
+
+    int nextLevel = (*level) + 1;
+    (*level) = nextLevel;
+    surf0.AllPointsIntersectingUntrimmed(a, b, cnt, level, l, segment, sorig);
+    (*level) = nextLevel;
+    surf1.AllPointsIntersectingUntrimmed(a, b, cnt, level, l, segment, sorig);
+}
+
 //-----------------------------------------------------------------------------
 // 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
@@ -119,55 +299,17 @@ void SSurface::AllPointsIntersecting(Vector a, Vector b,
     Vector ba = b.Minus(a);
     double bam = ba.Magnitude();
 
-    typedef struct {
-        int     tag;
-        Point2d p;
-    } Inter;
     List<Inter> inters;
     ZERO(&inters);
 
     // First, get all the intersections between the infinite ray and the
     // untrimmed surface.
-    int i, j;
+    int cnt = 0, level = 0;
-    for(i = 0; i < degm; i++) {
+    AllPointsIntersectingUntrimmed(a, b, &cnt, &level, &inters, seg, this);
-        for(j = 0; j < degn; j++) {
-            // Intersect the ray with each face in the control polyhedron
-            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, false);
-            if(PointIntersectingLine(a, b, &inter.p.x, &inter.p.y)) {
-                inters.Add(&inter);
-            } else {
-                // No convergence; but this might not be an error, since
-                // the line can intersect the convex hull more times than
-                // it intersects the surface itself.
-            }
-        }
-    }
 
     // Remove duplicate intersection points
     inters.ClearTags();
+    int i, j;
     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)) {
@@ -246,11 +388,6 @@ int SShell::ClassifyPoint(Vector p, Vector pout) {
         bool onEdge = false;
         edge_inters = 0;
 
-        if(FLAG) {
-            dbp("inters=%d cnt=%d", l.n, cnt);
-            SS.nakedEdges.AddEdge(p, p.Plus(ray.WithMagnitude(50)));
-        }
-
         SInter *si;
         for(si = l.First(); si; si = l.NextAfter(si)) {
             double t = ((si->p).Minus(p)).DivPivoting(ray);
@@ -259,8 +396,6 @@ int SShell::ClassifyPoint(Vector p, Vector pout) {
                 continue;
             }
 
-            if(FLAG) SS.nakedEdges.AddEdge(p, si->p);
-
             double d = ((si->p).Minus(p)).Magnitude();
 
             // Handle edge-on-edge

srf/triangulate.cpp

@@ -299,7 +299,7 @@ void SContour::UvTriangulateInto(SMesh *m, SSurface *srf) {
                     bestEar = ear;
                     bestChordTol = tol;
                 }
-                if(bestChordTol < 0.1*(SS.chordTol / SS.GW.scale)) {
+                if(bestChordTol < 0.1*SS.ChordTolMm()) {
                     break;
                 }
             }

util.cpp

@@ -573,6 +573,40 @@ bool Vector::BoundingBoxesDisjoint(Vector amax, Vector amin,
     return false;
 }
 
+bool Vector::BoundingBoxIntersectsLine(Vector amax, Vector amin,
+                                       Vector p0, Vector p1, bool segment)
+{
+    Vector dp = p1.Minus(p0);
+    double lp = dp.Magnitude();
+    dp = dp.ScaledBy(1.0/lp);
+
+    int i, a;
+    for(i = 0; i < 3; i++) {
+        int j = WRAP(i+1, 3), k = WRAP(i+2, 3);
+        if(lp*fabs(dp.Element(i)) < LENGTH_EPS) continue; // parallel to plane
+
+        for(a = 0; a < 2; a++) {
+            double d = (a == 0) ? amax.Element(i) : amin.Element(i);
+            // n dot (p0 + t*dp) = d
+            // (n dot p0) + t * (n dot dp) = d
+            double t = (d - p0.Element(i)) / dp.Element(i);
+            Vector p = p0.Plus(dp.ScaledBy(t));
+
+            if(segment && (t < -LENGTH_EPS || t > (lp+LENGTH_EPS))) continue;
+
+            if(p.Element(j) > amax.Element(j) + LENGTH_EPS) continue;
+            if(p.Element(k) > amax.Element(k) + LENGTH_EPS) continue;
+
+            if(p.Element(j) < amin.Element(j) - LENGTH_EPS) continue;
+            if(p.Element(k) < amin.Element(k) - LENGTH_EPS) continue;
+
+            return true;
+        }
+    }
+
+    return false;
+}
+
 Vector Vector::AtIntersectionOfPlanes(Vector n1, double d1,
                                       Vector n2, double d2)
 {

wishlist.txt

@@ -1,11 +1,19 @@
 
-leak fixing
+marching algorithm for surface intersection
+surfaces of revolution (lathed)
+fix surface-line intersection
+cylinder-line and plane-line special cases
+exact boundaries when near pwl trim
+step and repeat rotate/translate
+short pwl edge avoidance
+faces
+take consistent pwl with coincident faces
+exact curve export (at least for dxf)
+hidden line removal from mesh
+line styles (color, thickness)
 
-
-long term
 -----
+loop detection
 incremental regen of entities?
-my own plane poly triangulation code
-exact curved surfaces