From adc910185c353aafc2729f635ab4a6efcb650260 Mon Sep 17 00:00:00 2001 From: Jonathan Westhues Date: Sat, 14 Mar 2009 12:01:20 -0800 Subject: [PATCH] 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 | 2 +- draw.cpp | 4 +- dsc.h | 2 + srf/boolean.cpp | 1 - srf/ratpoly.cpp | 72 +++++++++++++++++++- srf/surface.h | 2 + srf/surfinter.cpp | 167 +++++++++++++++++++++++++++++++++++++++++----- textscreens.cpp | 5 ++ util.cpp | 6 ++ win32/w32main.cpp | 9 +-- wishlist.txt | 4 +- 11 files changed, 246 insertions(+), 28 deletions(-) diff --git a/Makefile b/Makefile index 14c5eeb..62a102d 100644 --- a/Makefile +++ b/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/* diff --git a/draw.cpp b/draw.cpp index 40bfaf5..077c224 100644 --- a/draw.cpp +++ b/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); diff --git a/dsc.h b/dsc.h index dc64f41..e138f7c 100644 --- a/dsc.h +++ b/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); diff --git a/srf/boolean.cpp b/srf/boolean.cpp index f298ec2..2cb4f7c 100644 --- a/srf/boolean.cpp +++ b/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)) { diff --git a/srf/ratpoly.cpp b/srf/ratpoly.cpp index 4faa43c..dcd05bd 100644 --- a/srf/ratpoly.cpp +++ b/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); diff --git a/srf/surface.h b/srf/surface.h index aa1e042..8544bd3 100644 --- a/srf/surface.h +++ b/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 *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); diff --git a/srf/surfinter.cpp b/srf/surfinter.cpp index 3352987..9c3507d 100644 --- a/srf/surfinter.cpp +++ b/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(sb->deg == 1) { + if(0 && 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)) || - (b->degm == 1 && b->degn == 1 && this->IsExtrusion(NULL, NULL))) + if(degm == 1 && degn == 1 && b->degm == 1 && b->degn == 1) { + // 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, - max(amax.y - amin.y, - amax.z - amin.z)); - if(fabs(h) < SS.ChordTolMm()) { + if(DepartureFromCoplanar() < 0.2*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 inters; ZERO(&inters); - // First, get all the intersections between the infinite ray and the - // untrimmed surface. - int cnt = 0, level = 0; - AllPointsIntersectingUntrimmed(a, b, &cnt, &level, &inters, seg, this); + // All the intersections between the line and the surface; either special + // cases that we can quickly solve in closed form, or general numerical. + Vector center, axis; + 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(¢er, &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(); diff --git a/textscreens.cpp b/textscreens.cpp index 82bb781..61eab6f 100644 --- a/textscreens.cpp +++ b/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)); } //----------------------------------------------------------------------------- diff --git a/util.cpp b/util.cpp index c5d3328..c5bd80f 100644 --- a/util.cpp +++ b/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; diff --git a/win32/w32main.cpp b/win32/w32main.cpp index 509fe2e..eaeb94a 100644 --- a/win32/w32main.cpp +++ b/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); - OutputDebugString(buf); + _vsnprintf(buf, sizeof(buf), str, f); 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); diff --git a/wishlist.txt b/wishlist.txt index e68b78b..4bbb8b9 100644 --- a/wishlist.txt +++ b/wishlist.txt @@ -1,10 +1,10 @@ marching algorithm for surface intersection surfaces of revolution (lathed) -fix surface-line intersection -cylinder-line and plane-line special cases +cylinder-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