Add plane-plane intersection as a special case (to generate the

trimmed line), and plane-line intersection. Terminate the Bezier
surface subdivision on a chord tolerance, and that seems okay now.
And print info about the graphics adapter in the text window, could
be useful.

Also have a cylinder-detection routine that works; should special
case those surfaces in closed form since they are common, but not
doing it yet.

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

Makefile

@@ -53,7 +53,7 @@ LIBS = user32.lib gdi32.lib comctl32.lib advapi32.lib shell32.lib opengl32.lib g
 
 all: $(OBJDIR)/solvespace.exe
     @cp $(OBJDIR)/solvespace.exe .
-    solvespace t.slvs
+    solvespace tttt.slvs
 
 clean:
 	rm -f obj/*

draw.cpp

@@ -1024,10 +1024,10 @@ void GraphicsWindow::Paint(int w, int h) {
     // And the naked edges, if the user did Analyze -> Show Naked Edges.
     glLineWidth(7);
     glEnable(GL_LINE_STIPPLE);
-    glLineStipple(4, 0x5555);
+    glLineStipple(1, 0x5555);
     glColor3d(1, 0, 0);
     glxDrawEdges(&(SS.nakedEdges));
-    glLineStipple(4, 0xaaaa);
+    glLineStipple(1, 0xaaaa);
     glColor3d(0, 0, 0);
     glxDrawEdges(&(SS.nakedEdges));
     glDisable(GL_LINE_STIPPLE);

dsc.h

@@ -95,6 +95,8 @@ class Point2d {
 public:
     double x, y;
 
+    static Point2d From(double x, double y);
+
     Point2d Plus(Point2d b);
     Point2d Minus(Point2d b);
     Point2d ScaledBy(double s);

srf/boolean.cpp

@@ -73,7 +73,6 @@ SCurve SCurve::MakeCopySplitAgainst(SShell *agnstA, SShell *agnstB,
             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)) {

srf/ratpoly.cpp

@@ -50,7 +50,42 @@ static double Bernstein(int k, int deg, double t)
 
 double BernsteinDerivative(int k, int deg, double t)
 {
-    return deg*(Bernstein(k-1, deg-1, t) - Bernstein(k, deg-1, t));
+    switch(deg) {
+        case 0:
+            return 0;
+            break;
+
+        case 1:
+            if(k == 0) {
+                return -1;
+            } else if(k = 1) {
+                return 1;
+            }
+            break;
+
+        case 2:
+            if(k == 0) {
+                return -2 + 2*t;
+            } else if(k == 1) {
+                return 2 - 4*t;
+            } else if(k == 2) {
+                return 2*t;
+            }
+            break;
+
+        case 3:
+            if(k == 0) {
+                return -3 + 6*t - 3*t*t;
+            } else if(k == 1) {
+                return 3 - 12*t + 9*t*t;
+            } else if(k == 2) {
+                return 6*t - 9*t*t;
+            } else if(k == 3) {
+                return 3*t*t;
+            }
+            break;
+    }
+    oops();
 }
 
 SBezier SBezier::From(Vector p0, Vector p1) {
@@ -418,6 +453,41 @@ bool SSurface::IsExtrusion(SBezier *of, Vector *alongp) {
     return true;
 }
 
+bool SSurface::IsCylinder(Vector *center, Vector *axis, double *r,
+                             double *dtheta)
+{
+    SBezier sb;
+    if(!IsExtrusion(&sb, axis)) return false;
+    if(sb.deg != 2) return false;
+
+    Vector t0 = (sb.ctrl[0]).Minus(sb.ctrl[1]),
+           t2 = (sb.ctrl[2]).Minus(sb.ctrl[1]),
+           r0 = axis->Cross(t0),
+           r2 = axis->Cross(t2);
+
+    *center = Vector::AtIntersectionOfLines(sb.ctrl[0], (sb.ctrl[0]).Plus(r0),
+                                            sb.ctrl[2], (sb.ctrl[2]).Plus(r2),
+                                            NULL, NULL, NULL);
+
+    double rd0 = center->Minus(sb.ctrl[0]).Magnitude(),
+           rd2 = center->Minus(sb.ctrl[2]).Magnitude();
+    if(fabs(rd0 - rd2) > LENGTH_EPS) return false;
+
+    Vector u = r0.WithMagnitude(1),
+           v = (axis->Cross(u)).WithMagnitude(1);
+    Point2d c2  = center->Project2d(u, v),
+            pa2 = (sb.ctrl[0]).Project2d(u, v).Minus(c2),
+            pb2 = (sb.ctrl[2]).Project2d(u, v).Minus(c2);
+    
+    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(fabs(sb.weight[1] - cos(*dtheta/2)) > LENGTH_EPS) return false;
+
+    return true;
+}
+
 SSurface SSurface::FromPlane(Vector pt, Vector u, Vector v) {
     SSurface ret;
     ZERO(&ret);

srf/surface.h

@@ -203,6 +203,7 @@ public:
     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);
+    double DepartureFromCoplanar(void);
     void SplitInHalf(bool byU, SSurface *sa, SSurface *sb);
     void AllPointsIntersecting(Vector a, Vector b,
                                     List<SInter> *l, bool seg, bool trimmed);
@@ -221,6 +222,7 @@ public:
     bool CoincidentWithPlane(Vector n, double d);
     bool CoincidentWith(SSurface *ss, bool sameNormal);
     bool IsExtrusion(SBezier *of, Vector *along);
+    bool IsCylinder(Vector *center, Vector *axis, double *r, double *dtheta);
 
     void TriangulateInto(SShell *shell, SMesh *sm);
     void MakeEdgesInto(SShell *shell, SEdgeList *sel, bool asUv);

srf/surfinter.cpp

@@ -15,8 +15,7 @@ void SSurface::AddExactIntersectionCurve(SBezier *sb, SSurface *srfB,
     SCurve split = sc.MakeCopySplitAgainst(agnstA, agnstB, this, srfB);
     sc.Clear();
    
-/*
+    if(0 && sb->deg == 1) {
-    if(sb->deg == 1) {
         dbp(" ");
         Vector *prev = NULL, *v;
         dbp("split.pts.n =%d", split.pts.n);
@@ -26,7 +25,7 @@ void SSurface::AddExactIntersectionCurve(SBezier *sb, SSurface *srfB,
             }
             prev = v;
         }
-    } */
+    }
 
     split.source = SCurve::FROM_INTERSECTION;
     into->curve.AddAndAssignId(&split);
@@ -44,8 +43,68 @@ void SSurface::IntersectAgainst(SSurface *b, SShell *agnstA, SShell *agnstB,
         return;
     }
 
-    if((degm == 1 && degn == 1 && b->IsExtrusion(NULL, NULL)) ||
+    if(degm == 1 && degn == 1 && b->degm == 1 && b->degn == 1) {
-       (b->degm == 1 && b->degn == 1 && this->IsExtrusion(NULL, NULL)))
+        // Line-line intersection; it's a plane or nothing.
+        Vector na = NormalAt(0, 0).WithMagnitude(1),
+               nb = b->NormalAt(0, 0).WithMagnitude(1);
+        double da = na.Dot(PointAt(0, 0)),
+               db = nb.Dot(b->PointAt(0, 0));
+
+        Vector dl = na.Cross(nb);
+        if(dl.Magnitude() < LENGTH_EPS) return; // parallel planes
+        dl = dl.WithMagnitude(1);
+        Vector p = Vector::AtIntersectionOfPlanes(na, da, nb, db);
+
+        // Trim it to the region 0 <= {u,v} <= 1 for each plane; not strictly
+        // necessary, since line will be split and excess edges culled, but
+        // this improves speed and robustness.
+        int i;
+        double tmax = VERY_POSITIVE, tmin = VERY_NEGATIVE;
+        for(i = 0; i < 2; i++) {
+            SSurface *s = (i == 0) ? this : b;
+            Vector tu, tv;
+            s->TangentsAt(0, 0, &tu, &tv);
+
+            double up, vp, ud, vd;
+            s->ClosestPointTo(p, &up, &vp);
+            ud = (dl.Dot(tu)) / tu.MagSquared();
+            vd = (dl.Dot(tv)) / tv.MagSquared();
+
+            // so u = up + t*ud
+            //    v = vp + t*vd
+            if(ud > LENGTH_EPS) {
+                tmin = max(tmin, -up/ud);
+                tmax = min(tmax, (1 - up)/ud);
+            } else if(ud < -LENGTH_EPS) {
+                tmax = min(tmax, -up/ud);
+                tmin = max(tmin, (1 - up)/ud);
+            } else {
+                if(up < -LENGTH_EPS || up > 1 + LENGTH_EPS) {
+                    // u is constant, and outside [0, 1]
+                    tmax = VERY_NEGATIVE;
+                }
+            }
+            if(vd > LENGTH_EPS) {
+                tmin = max(tmin, -vp/vd);
+                tmax = min(tmax, (1 - vp)/vd);
+            } else if(vd < -LENGTH_EPS) {
+                tmax = min(tmax, -vp/vd);
+                tmin = max(tmin, (1 - vp)/vd);
+            } else {
+                if(vp < -LENGTH_EPS || vp > 1 + LENGTH_EPS) {
+                    // v is constant, and outside [0, 1]
+                    tmax = VERY_NEGATIVE;
+                }
+            }
+        }
+
+        if(tmax > tmin) {
+            SBezier bezier = SBezier::From(p.Plus(dl.ScaledBy(tmin)),
+                                           p.Plus(dl.ScaledBy(tmax)));
+            AddExactIntersectionCurve(&bezier, b, agnstA, agnstB, into);
+        }
+    } else if((degm == 1 && degn == 1 && b->IsExtrusion(NULL, NULL)) ||
+              (b->degm == 1 && b->degn == 1 && this->IsExtrusion(NULL, NULL)))
     {
         // The intersection between a plane and a surface of extrusion
         SSurface *splane, *sext;
@@ -95,10 +154,6 @@ void SSurface::IntersectAgainst(SSurface *b, SShell *agnstA, SShell *agnstB,
         } else {
             // Direction of extrusion is not parallel to plane; so
             // 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++) {
                 Vector p0 = bezier.ctrl[i],
@@ -116,6 +171,64 @@ void SSurface::IntersectAgainst(SSurface *b, SShell *agnstA, SShell *agnstB,
     // cases, just giving up for now
 }
 
+
+double SSurface::DepartureFromCoplanar(void) {
+    int i, j;
+    int ia, ja, ib, jb, ic, jc;
+    double best;
+
+    // Grab three points to define a plane; first choose (0, 0) arbitrarily.
+    ia = ja = 0;
+    // Then the point farthest from pt a.
+    best = VERY_NEGATIVE;
+    for(i = 0; i <= degm; i++) {
+        for(j = 0; j <= degn; j++) {
+            if(i == ia && j == ja) continue;
+
+            double dist = (ctrl[i][j]).Minus(ctrl[ia][ja]).Magnitude();
+            if(dist > best) {
+                best = dist;
+                ib = i;
+                jb = j;
+            }
+        }
+    }
+    // Then biggest magnitude of ab cross ac.
+    best = VERY_NEGATIVE;
+    for(i = 0; i <= degm; i++) {
+        for(j = 0; j <= degn; j++) {
+            if(i == ia && j == ja) continue;
+            if(i == ib && j == jb) continue;
+
+            double mag =
+                ((ctrl[ia][ja].Minus(ctrl[ib][jb]))).Cross(
+                 (ctrl[ia][ja].Minus(ctrl[i ][j ]))).Magnitude();
+            if(mag > best) {
+                best = mag;
+                ic = i;
+                jc = j;
+            }
+        }
+    }
+
+    Vector n = ((ctrl[ia][ja].Minus(ctrl[ib][jb]))).Cross(
+                (ctrl[ia][ja].Minus(ctrl[ic][jc])));
+    n = n.WithMagnitude(1);
+    double d = (ctrl[ia][ja]).Dot(n);
+
+    // Finally, calculate the deviation from each point to the plane.
+    double farthest = VERY_NEGATIVE;
+    for(i = 0; i <= degm; i++) {
+        for(j = 0; j <= degn; j++) {
+            double dist = fabs(n.Dot(ctrl[i][j]) - d);
+            if(dist > farthest) {
+                farthest = dist;
+            }
+        }
+    }
+    return farthest;
+}
+
 void SSurface::WeightControlPoints(void) {
     int i, j;
     for(i = 0; i <= degm; i++) {
@@ -250,16 +363,14 @@ void SSurface::AllPointsIntersectingUntrimmed(Vector a, Vector b,
    
     if(*cnt > 2000) {
         dbp("!!! too many subdivisions (level=%d)!", *level);
+        dbp("degm = %d degn = %d", degm, degn);
         return;
     }
     (*cnt)++;
 
     // If we might intersect, and the surface is small, then switch to Newton
     // iterations.
-    double h = max(amax.x - amin.x,
+    if(DepartureFromCoplanar() < 0.2*SS.ChordTolMm()) {
-               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);
@@ -302,10 +413,32 @@ void SSurface::AllPointsIntersecting(Vector a, Vector b,
     List<Inter> inters;
     ZERO(&inters);
 
-    // First, get all the intersections between the infinite ray and the
+    // All the intersections between the line and the surface; either special
-    // untrimmed surface.
+    // cases that we can quickly solve in closed form, or general numerical.
-    int cnt = 0, level = 0;
+    Vector center, axis;
-    AllPointsIntersectingUntrimmed(a, b, &cnt, &level, &inters, seg, this);
+    double radius, dtheta;
+    if(degm == 1 && degn == 1) {
+        // Against a plane, easy.
+        Vector n = NormalAt(0, 0).WithMagnitude(1);
+        double d = n.Dot(PointAt(0, 0));
+        // Trim to line segment now if requested, don't generate points that
+        // would just get discarded later.
+        if(!seg ||
+           (n.Dot(a) > d + LENGTH_EPS && n.Dot(b) < d - LENGTH_EPS) ||
+           (n.Dot(b) > d + LENGTH_EPS && n.Dot(a) < d - LENGTH_EPS))
+        {
+            Vector p = Vector::AtIntersectionOfPlaneAndLine(n, d, a, b, NULL);
+            Inter inter;
+            ClosestPointTo(p, &(inter.p.x), &(inter.p.y));
+            inters.Add(&inter);
+        }
+    } else if(IsCylinder(&center, &axis, &radius, &dtheta) && 0) {
+        // XXX, cylinder is easy in closed form
+    } else {
+        // General numerical solution by subdivision, fallback
+        int cnt = 0, level = 0;
+        AllPointsIntersectingUntrimmed(a, b, &cnt, &level, &inters, seg, this);
+    }
 
     // Remove duplicate intersection points
     inters.ClearTags();

textscreens.cpp

@@ -672,6 +672,11 @@ void TextWindow::ShowConfiguration(void) {
         (SS.drawBackFaces ? "" : "yes"), (SS.drawBackFaces ? "yes" : ""),
         &ScreenChangeBackFaces,
         (!SS.drawBackFaces ? "" : "no"), (!SS.drawBackFaces ? "no" : ""));
+    
+    Printf(false, "");
+    Printf(false, " %Ftgl vendor   %E%s", glGetString(GL_VENDOR));
+    Printf(false, " %Ft   renderer %E%s", glGetString(GL_RENDERER));
+    Printf(false, " %Ft   version  %E%s", glGetString(GL_VERSION));
 }
 
 //-----------------------------------------------------------------------------

util.cpp

@@ -716,6 +716,12 @@ Vector Vector::AtIntersectionOfPlanes(Vector na, double da,
     return Vector::From(detx/det, dety/det, detz/det);
 }
 
+Point2d Point2d::From(double x, double y) {
+    Point2d r;
+    r.x = x; r.y = y;
+    return r;
+}
+
 Point2d Point2d::Plus(Point2d b) {
     Point2d r;
     r.x = x + b.x;

win32/w32main.cpp

@@ -54,9 +54,10 @@ void dbp(char *str, ...)
     va_list f;
     static char buf[1024*50];
     va_start(f, str);
-    vsprintf(buf, str, f);
+    _vsnprintf(buf, sizeof(buf), str, f);
-    OutputDebugString(buf);
     va_end(f);
+
+    OutputDebugString(buf);
 }
 
 
@@ -603,7 +604,7 @@ static void CreateGlContext(void)
     HDC hdc = GetDC(GraphicsWnd);
 
     PIXELFORMATDESCRIPTOR pfd;
-    int pixelFormat; 
+    int pixelFormat;
 
     memset(&pfd, 0, sizeof(pfd));
     pfd.nSize = sizeof(PIXELFORMATDESCRIPTOR); 
@@ -623,7 +624,7 @@ static void CreateGlContext(void)
     if(!SetPixelFormat(hdc, pixelFormat, &pfd)) oops();
 
     GraphicsHpgl = wglCreateContext(hdc); 
-    wglMakeCurrent(hdc, GraphicsHpgl); 
+    wglMakeCurrent(hdc, GraphicsHpgl);
 
     // Create a bitmap font in a display list, for DrawWithBitmapFont().
     SelectObject(hdc, FixedFont);

wishlist.txt

@@ -1,10 +1,10 @@
 
 marching algorithm for surface intersection
 surfaces of revolution (lathed)
-fix surface-line intersection
+cylinder-line special cases
-cylinder-line and plane-line special cases
 exact boundaries when near pwl trim
 step and repeat rotate/translate
+tangent intersections
 short pwl edge avoidance
 faces
 take consistent pwl with coincident faces