diff --git a/Makefile b/Makefile index 09a502b7..f3dbebe7 100644 --- a/Makefile +++ b/Makefile @@ -1,7 +1,7 @@ DEFINES = /D_WIN32_WINNT=0x500 /DISOLATION_AWARE_ENABLED /D_WIN32_IE=0x500 /DWIN32_LEAN_AND_MEAN /DWIN32 # Use the multi-threaded static libc because libpng and zlib do; not sure if anything bad # happens if those mix, but don't want to risk it. -CFLAGS = /W3 /nologo -MT -Iextlib -I..\common\win32 /D_DEBUG /D_CRT_SECURE_NO_WARNINGS /I. /Zi /EHs +CFLAGS = /W3 /nologo -MT -Iextlib -I..\common\win32 /D_DEBUG /D_CRT_SECURE_NO_WARNINGS /I. /Zi /EHs /O2 HEADERS = ..\common\win32\freeze.h ui.h solvespace.h dsc.h sketch.h expr.h polygon.h srf\surface.h @@ -61,7 +61,7 @@ LIBS = user32.lib gdi32.lib comctl32.lib advapi32.lib shell32.lib opengl32.lib g all: $(OBJDIR)/solvespace.exe @cp $(OBJDIR)/solvespace.exe . - solvespace t8.slvs + solvespace fail.slvs clean: rm -f obj/* diff --git a/entity.cpp b/entity.cpp index b03247d7..b2aeae23 100644 --- a/entity.cpp +++ b/entity.cpp @@ -685,8 +685,9 @@ void EntityBase::GenerateEquations(IdList *l) { // If the two endpoints of the arc are constrained coincident // (to make a complete circle), then our distance constraint // would be redundant and therefore overconstrain things. - Constraint *c; - for(c = SK.constraint.First(); c; c = SK.constraint.NextAfter(c)) { + int i; + for(i = 0; i < SK.constraint.n; i++) { + ConstraintBase *c = &(SK.constraint.elem[i]); if(c->group.v != group.v) continue; if(c->type != Constraint::POINTS_COINCIDENT) continue; @@ -696,7 +697,7 @@ void EntityBase::GenerateEquations(IdList *l) { break; } } - if(c) break; + if(i < SK.constraint.n) break; Expr *ra = Constraint::Distance(workplane, point[0], point[1]); Expr *rb = Constraint::Distance(workplane, point[0], point[2]); diff --git a/exposed/Makefile b/exposed/Makefile index 90b0c19a..67b1eec7 100644 --- a/exposed/Makefile +++ b/exposed/Makefile @@ -1,7 +1,7 @@ DEFINES = /D_WIN32_WINNT=0x500 /DISOLATION_AWARE_ENABLED /D_WIN32_IE=0x500 /DWIN32_LEAN_AND_MEAN /DWIN32 /DLIBRARY # Use the multi-threaded static libc because libpng and zlib do; not sure if anything bad # happens if those mix, but don't want to risk it. -CFLAGS = /W3 /nologo -MT -Iextlib -I..\..\common\win32 /D_DEBUG /D_CRT_SECURE_NO_WARNINGS /I. /I.. /Zi /EHs +CFLAGS = /W3 /nologo -MT -Iextlib -I..\..\common\win32 /D_DEBUG /D_CRT_SECURE_NO_WARNINGS /I. /I.. /Zi /EHs /O2 HEADERS = ..\solvespace.h ..\dsc.h ..\sketch.h ..\expr.h slvs.h diff --git a/polygon.h b/polygon.h index ac3b6bf0..6e811086 100644 --- a/polygon.h +++ b/polygon.h @@ -76,9 +76,9 @@ public: void FindPointWithMinX(void); Vector AnyEdgeMidpoint(void); - bool IsEar(int bp); + bool IsEar(int bp, double scaledEps); bool BridgeToContour(SContour *sc, SEdgeList *el, List *vl); - void ClipEarInto(SMesh *m, int bp); + void ClipEarInto(SMesh *m, int bp, double scaledEps); void UvTriangulateInto(SMesh *m, SSurface *srf); }; diff --git a/srf/boolean.cpp b/srf/boolean.cpp index f6ab4808..d87b32a1 100644 --- a/srf/boolean.cpp +++ b/srf/boolean.cpp @@ -48,9 +48,9 @@ SCurve SCurve::MakeCopySplitAgainst(SShell *agnstA, SShell *agnstB, // Find all the intersections with the two passed shells if(agnstA) - agnstA->AllPointsIntersecting(prev.p, p->p, &il, true, true, true); + agnstA->AllPointsIntersecting(prev.p, p->p, &il, true, false, true); if(agnstB) - agnstB->AllPointsIntersecting(prev.p, p->p, &il, true, true, true); + agnstB->AllPointsIntersecting(prev.p, p->p, &il, true, false, true); if(il.n > 0) { // The intersections were generated by intersecting the pwl @@ -69,10 +69,28 @@ SCurve SCurve::MakeCopySplitAgainst(SShell *agnstA, SShell *agnstB, pi->tag = 1; continue; } - 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); + + Point2d puv; + (pi->srf)->ClosestPointTo(pi->p, &puv, false); + + // Split the edge if the intersection lies within the surface's + // trim curves, or within the chord tol of the trim curve; want + // some slop if points are close to edge and pwl is too coarse, + // and it doesn't hurt to split unnecessarily. + Point2d dummy = { 0, 0 }; + int c = pi->srf->bsp->ClassifyPoint(puv, dummy, pi->srf); + if(c == SBspUv::OUTSIDE) { + double d; + d = pi->srf->bsp->MinimumDistanceToEdge(puv, pi->srf); + if(d > SS.ChordTolMm()) { + pi->tag = 1; + continue; + } + } + + // We're keeping the intersection, so actually refine it. + (pi->srf)->PointOnSurfaces(srfA, srfB, &(puv.x), &(puv.y)); + pi->p = (pi->srf)->PointAt(puv); } il.RemoveTagged(); @@ -423,7 +441,7 @@ SSurface SSurface::MakeCopyTrimAgainst(SShell *parent, // the split curves (not the original curves). SEdgeList orig; ZERO(&orig); - ret.MakeEdgesInto(into, &orig, true); + ret.MakeEdgesInto(into, &orig, AS_UV); ret.trim.Clear(); // which means that we can't necessarily use the old BSP... SBspUv *origBsp = SBspUv::From(&orig, &ret); @@ -763,12 +781,12 @@ void SSurface::MakeClassifyingBsp(SShell *shell, SShell *useCurvesFrom) { SEdgeList el; ZERO(&el); - MakeEdgesInto(shell, &el, true, useCurvesFrom); + MakeEdgesInto(shell, &el, AS_UV, useCurvesFrom); bsp = SBspUv::From(&el, this); el.Clear(); ZERO(&edges); - MakeEdgesInto(shell, &edges, false, useCurvesFrom); + MakeEdgesInto(shell, &edges, AS_XYZ, useCurvesFrom); } SBspUv *SBspUv::Alloc(void) { @@ -942,3 +960,16 @@ int SBspUv::ClassifyEdge(Point2d ea, Point2d eb, SSurface *srf) { return ret; } +double SBspUv::MinimumDistanceToEdge(Point2d p, SSurface *srf) { + if(!this) return VERY_POSITIVE; + + double dn = neg->MinimumDistanceToEdge(p, srf), + dp = pos->MinimumDistanceToEdge(p, srf); + + Point2d as = a, bs = b; + ScalePoints(&p, &as, &bs, srf); + double d = p.DistanceToLine(as, bs.Minus(as), true); + + return min(d, min(dn, dp)); +} + diff --git a/srf/merge.cpp b/srf/merge.cpp index b9f9bbb7..0a37e893 100644 --- a/srf/merge.cpp +++ b/srf/merge.cpp @@ -22,7 +22,7 @@ void SShell::MergeCoincidentSurfaces(void) { SEdgeList sel; ZERO(&sel); - si->MakeEdgesInto(this, &sel, false); + si->MakeEdgesInto(this, &sel, SSurface::AS_XYZ); bool mergedThisTime, merged = false; do { @@ -42,7 +42,7 @@ void SShell::MergeCoincidentSurfaces(void) { // less robust. SEdgeList tel; ZERO(&tel); - sj->MakeEdgesInto(this, &tel, false); + sj->MakeEdgesInto(this, &tel, SSurface::AS_XYZ); if(!sel.ContainsEdgeFrom(&tel)) { tel.Clear(); continue; @@ -52,7 +52,7 @@ void SShell::MergeCoincidentSurfaces(void) { sj->tag = 1; merged = true; mergedThisTime = true; - sj->MakeEdgesInto(this, &sel, false); + sj->MakeEdgesInto(this, &sel, SSurface::AS_XYZ); sj->trim.Clear(); // All the references to this surface get replaced with the diff --git a/srf/ratpoly.cpp b/srf/ratpoly.cpp index af9aafac..a7a32d72 100644 --- a/srf/ratpoly.cpp +++ b/srf/ratpoly.cpp @@ -356,6 +356,19 @@ void SSurface::ClosestPointTo(Vector p, double *u, double *v, bool converge) { if(p.Equals(ctrl[degm][degn])) { *u = 1; *v = 1; return; } if(p.Equals(ctrl[0] [degn])) { *u = 0; *v = 1; return; } + // Try whatever the previous guess was. This is likely to do something + // good if we're working our way along a curve or something else where + // we project successive points that are close to each other; something + // like a 20% speedup empirically. + if(converge) { + double ut = cached.x, vt = cached.y; + if(ClosestPointNewton(p, &ut, &vt, converge)) { + cached.x = *u = ut; + cached.y = *v = vt; + return; + } + } + // Search for a reasonable initial guess int i, j; if(degm == 1 && degn == 1) { @@ -378,13 +391,27 @@ void SSurface::ClosestPointTo(Vector p, double *u, double *v, bool converge) { } } + if(ClosestPointNewton(p, u, v, converge)) { + cached.x = *u; + cached.y = *v; + return; + } + + // If we failed to converge, then at least don't return NaN. + if(isnan(*u) || isnan(*v)) { + *u = *v = 0; + } +} + +bool SSurface::ClosestPointNewton(Vector p, double *u, double *v, bool converge) +{ // Initial guess is in u, v; refine by Newton iteration. Vector p0; - for(i = 0; i < (converge ? 25 : 5); i++) { + for(int i = 0; i < (converge ? 25 : 5); i++) { p0 = PointAt(*u, *v); if(converge) { if(p0.Equals(p, RATPOLY_EPS)) { - return; + return true; } } @@ -406,10 +433,7 @@ void SSurface::ClosestPointTo(Vector p, double *u, double *v, bool converge) { dbp("want %.3f %.3f %.3f", CO(p)); dbp("distance = %g", (p.Minus(p0)).Magnitude()); } - - if(isnan(*u) || isnan(*v)) { - *u = *v = 0; - } + return false; } bool SSurface::PointIntersectingLine(Vector p0, Vector p1, double *u, double *v) diff --git a/srf/surface.cpp b/srf/surface.cpp index 34433003..699553c0 100644 --- a/srf/surface.cpp +++ b/srf/surface.cpp @@ -194,10 +194,10 @@ bool SSurface::LineEntirelyOutsideBbox(Vector a, Vector b, bool segment) { // Generate the piecewise linear approximation of the trim stb, which applies // to the curve sc. //----------------------------------------------------------------------------- -void SSurface::MakeTrimEdgesInto(SEdgeList *sel, bool asUv, +void SSurface::MakeTrimEdgesInto(SEdgeList *sel, int flags, SCurve *sc, STrimBy *stb) { - Vector prev, prevuv, ptuv; + Vector prev; bool inCurve = false, empty = true; double u = 0, v = 0; @@ -212,23 +212,22 @@ void SSurface::MakeTrimEdgesInto(SEdgeList *sel, bool asUv, increment = 1; } for(i = first; i != (last + increment); i += increment) { - Vector *pt = &(sc->pts.elem[i].p); - if(asUv) { + Vector tpt, *pt = &(sc->pts.elem[i].p); + + if(flags & AS_UV) { ClosestPointTo(*pt, &u, &v); - ptuv = Vector::From(u, v, 0); - if(inCurve) { - sel->AddEdge(prevuv, ptuv, sc->h.v, stb->backwards); - empty = false; - } - prevuv = ptuv; + tpt = Vector::From(u, v, 0); } else { - if(inCurve) { - sel->AddEdge(prev, *pt, sc->h.v, stb->backwards); - empty = false; - } - prev = *pt; + tpt = *pt; } + if(inCurve) { + sel->AddEdge(prev, tpt, sc->h.v, stb->backwards); + empty = false; + } + + prev = tpt; // either uv or xyz, depending on flags + if(pt->Equals(stb->start)) inCurve = true; if(pt->Equals(stb->finish)) inCurve = false; } @@ -242,7 +241,7 @@ void SSurface::MakeTrimEdgesInto(SEdgeList *sel, bool asUv, // the split curves from useCurvesFrom instead of the curves in our own // shell. //----------------------------------------------------------------------------- -void SSurface::MakeEdgesInto(SShell *shell, SEdgeList *sel, bool asUv, +void SSurface::MakeEdgesInto(SShell *shell, SEdgeList *sel, int flags, SShell *useCurvesFrom) { STrimBy *stb; @@ -257,7 +256,7 @@ void SSurface::MakeEdgesInto(SShell *shell, SEdgeList *sel, bool asUv, sc = useCurvesFrom->curve.FindById(sc->newH); } - MakeTrimEdgesInto(sel, asUv, sc, stb); + MakeTrimEdgesInto(sel, flags, sc, stb); } } @@ -384,7 +383,7 @@ void SSurface::MakeSectionEdgesInto(SShell *shell, sp = fpt; } } else { - if(sel) MakeTrimEdgesInto(sel, false, sc, stb); + if(sel) MakeTrimEdgesInto(sel, AS_XYZ, sc, stb); } } } @@ -393,21 +392,21 @@ void SSurface::TriangulateInto(SShell *shell, SMesh *sm) { SEdgeList el; ZERO(&el); - MakeEdgesInto(shell, &el, true); + MakeEdgesInto(shell, &el, AS_UV); SPolygon poly; ZERO(&poly); if(el.AssemblePolygon(&poly, NULL, true)) { int i, start = sm->l.n; if(degm == 1 && degn == 1) { - // A plane; triangulate any old way - poly.UvTriangulateInto(sm, NULL); - } else if(degm == 1 || degn == 1) { // A surface with curvature along one direction only; so // choose the triangulation with chords that lie as much // as possible within the surface. And since the trim curves // have been pwl'd to within the desired chord tol, that will // produce a surface good to within roughly that tol. + // + // If this is just a plane (degree (1, 1)) then the triangulation + // code will notice that, and not bother checking chord tols. poly.UvTriangulateInto(sm, this); } else { // A surface with compound curvature. So we must overlay a @@ -820,7 +819,7 @@ void SShell::MakeFromTransformationOf(SShell *a, Vector t, Quaternion q) { void SShell::MakeEdgesInto(SEdgeList *sel) { SSurface *s; for(s = surface.First(); s; s = surface.NextAfter(s)) { - s->MakeEdgesInto(this, sel, false); + s->MakeEdgesInto(this, sel, SSurface::AS_XYZ); } } diff --git a/srf/surface.h b/srf/surface.h index 01d52887..99f252c2 100644 --- a/srf/surface.h +++ b/srf/surface.h @@ -38,6 +38,7 @@ public: SBspUv *InsertEdge(Point2d a, Point2d b, SSurface *srf); int ClassifyPoint(Point2d p, Point2d eb, SSurface *srf); int ClassifyEdge(Point2d ea, Point2d eb, SSurface *srf); + double MinimumDistanceToEdge(Point2d p, SSurface *srf); }; // Now the data structures to represent a shell of trimmed rational polynomial @@ -224,6 +225,10 @@ public: SBspUv *bsp; SEdgeList edges; + // For caching our initial (u, v) when doing Newton iterations to project + // a point into our surface. + Point2d cached; + static SSurface FromExtrusionOf(SBezier *spc, Vector t0, Vector t1); static SSurface FromRevolutionOf(SBezier *sb, Vector pt, Vector axis, double thetas, double thetaf); @@ -266,6 +271,8 @@ public: void ClosestPointTo(Vector p, Point2d *puv, bool converge=true); void ClosestPointTo(Vector p, double *u, double *v, bool converge=true); + bool ClosestPointNewton(Vector p, double *u, double *v, bool converge=true); + bool PointIntersectingLine(Vector p0, Vector p1, double *u, double *v); Vector ClosestPointOnThisAndSurface(SSurface *srf2, Vector p); void PointOnSurfaces(SSurface *s1, SSurface *s2, double *u, double *v); @@ -283,9 +290,14 @@ public: Vector *start, Vector *finish); void TriangulateInto(SShell *shell, SMesh *sm); - void MakeTrimEdgesInto(SEdgeList *sel, bool asUv, SCurve *sc, STrimBy *stb); - void MakeEdgesInto(SShell *shell, SEdgeList *sel, bool asUv, + + // these are intended as bitmasks, even though there's just one now + static const int AS_UV = 0x01; + static const int AS_XYZ = 0x00; + void MakeTrimEdgesInto(SEdgeList *sel, int flags, SCurve *sc, STrimBy *stb); + void MakeEdgesInto(SShell *shell, SEdgeList *sel, int flags, SShell *useCurvesFrom=NULL); + Vector ExactSurfaceTangentAt(Vector p, SSurface *srfA, SSurface *srfB, Vector dir); void MakeSectionEdgesInto(SShell *shell, SEdgeList *sel, SBezierList *sbl); diff --git a/srf/surfinter.cpp b/srf/surfinter.cpp index 7f9326bf..7fd49d7a 100644 --- a/srf/surfinter.cpp +++ b/srf/surfinter.cpp @@ -54,13 +54,41 @@ void SSurface::AddExactIntersectionCurve(SBezier *sb, SSurface *srfB, sc.Clear(); } - if(0 && sb->deg == 1) { + // Test if the curve lies entirely outside one of the + SCurvePt *scpt; + bool withinA = false, withinB = false; + for(scpt = split.pts.First(); scpt; scpt = split.pts.NextAfter(scpt)) { + double tol = 0.01; + Point2d puv; + ClosestPointTo(scpt->p, &puv); + if(puv.x > -tol && puv.x < 1 + tol && + puv.y > -tol && puv.y < 1 + tol) + { + withinA = true; + } + srfB->ClosestPointTo(scpt->p, &puv); + if(puv.x > -tol && puv.x < 1 + tol && + puv.y > -tol && puv.y < 1 + tol) + { + withinB = true; + } + // Break out early, no sense wasting time if we already have the answer. + if(withinA && withinB) break; + } + if(!(withinA && withinB)) { + // Intersection curve lies entirely outside one of the surfaces, so + // it's fake. + split.Clear(); + return; + } + + if(sb->deg == 2 && 0) { dbp(" "); SCurvePt *prev = NULL, *v; - dbp("split.pts.n =%d", split.pts.n); + dbp("split.pts.n = %d", split.pts.n); for(v = split.pts.First(); v; v = split.pts.NextAfter(v)) { if(prev) { - Vector e = (prev->p).Minus(v->p).WithMagnitude(-1); + Vector e = (prev->p).Minus(v->p).WithMagnitude(0); SS.nakedEdges.AddEdge((prev->p).Plus(e), (v->p).Minus(e)); } prev = v; @@ -292,7 +320,7 @@ void SSurface::IntersectAgainst(SSurface *b, SShell *agnstA, SShell *agnstB, SEdgeList el; ZERO(&el); - srfA->MakeEdgesInto(shA, &el, false, NULL); + srfA->MakeEdgesInto(shA, &el, AS_XYZ, NULL); SEdge *se; for(se = el.l.First(); se; se = el.l.NextAfter(se)) { @@ -480,7 +508,7 @@ void SShell::MakeCoincidentEdgesInto(SSurface *proto, bool sameNormal, SSurface *ss; for(ss = surface.First(); ss; ss = surface.NextAfter(ss)) { if(proto->CoincidentWith(ss, sameNormal)) { - ss->MakeEdgesInto(this, el, false, useCurvesFrom); + ss->MakeEdgesInto(this, el, SSurface::AS_XYZ, useCurvesFrom); } } diff --git a/srf/triangulate.cpp b/srf/triangulate.cpp index 22917aac..8813292e 100644 --- a/srf/triangulate.cpp +++ b/srf/triangulate.cpp @@ -202,7 +202,7 @@ haveEdge: return true; } -bool SContour::IsEar(int bp) { +bool SContour::IsEar(int bp, double scaledEps) { int ap = WRAP(bp-1, l.n), cp = WRAP(bp+1, l.n); @@ -219,7 +219,7 @@ bool SContour::IsEar(int bp) { } Vector n = Vector::From(0, 0, -1); - if((tr.Normal()).Dot(n) < LENGTH_EPS) { + if((tr.Normal()).Dot(n) < scaledEps) { // This vertex is reflex, or between two collinear edges; either way, // it's not an ear. return false; @@ -251,7 +251,7 @@ bool SContour::IsEar(int bp) { return true; } -void SContour::ClipEarInto(SMesh *m, int bp) { +void SContour::ClipEarInto(SMesh *m, int bp, double scaledEps) { int ap = WRAP(bp-1, l.n), cp = WRAP(bp+1, l.n); @@ -260,7 +260,7 @@ void SContour::ClipEarInto(SMesh *m, int bp) { tr.a = l.elem[ap].p; tr.b = l.elem[bp].p; tr.c = l.elem[cp].p; - if(tr.Normal().MagSquared() < LENGTH_EPS*LENGTH_EPS) { + if(tr.Normal().MagSquared() < scaledEps*scaledEps) { // A vertex with more than two edges will cause us to generate // zero-area triangles, which must be culled. } else { @@ -278,8 +278,15 @@ void SContour::ClipEarInto(SMesh *m, int bp) { } void SContour::UvTriangulateInto(SMesh *m, SSurface *srf) { - int i; + Vector tu, tv; + srf->TangentsAt(0.5, 0.5, &tu, &tv); + double s = sqrt(tu.MagSquared() + tv.MagSquared()); + // We would like to apply our tolerances in xyz; but that would be a lot + // of work, so at least scale the epsilon semi-reasonably. That's + // perfect for square planes, less perfect for anything else. + double scaledEps = LENGTH_EPS / s; + int i; // Clean the original contour by removing any zero-length edges. l.ClearTags(); for(i = 1; i < l.n; i++) { @@ -291,7 +298,7 @@ void SContour::UvTriangulateInto(SMesh *m, SSurface *srf) { // Now calculate the ear-ness of each vertex for(i = 0; i < l.n; i++) { - (l.elem[i]).ear = IsEar(i) ? SPoint::EAR : SPoint::NOT_EAR; + (l.elem[i]).ear = IsEar(i, scaledEps) ? SPoint::EAR : SPoint::NOT_EAR; } bool toggle = false; @@ -299,7 +306,8 @@ void SContour::UvTriangulateInto(SMesh *m, SSurface *srf) { // Some points may have changed ear-ness, so recalculate for(i = 0; i < l.n; i++) { if(l.elem[i].ear == SPoint::UNKNOWN) { - (l.elem[i]).ear = IsEar(i) ? SPoint::EAR : SPoint::NOT_EAR; + (l.elem[i]).ear = IsEar(i, scaledEps) ? + SPoint::EAR : SPoint::NOT_EAR; } } @@ -312,7 +320,8 @@ void SContour::UvTriangulateInto(SMesh *m, SSurface *srf) { for(i = 0; i < l.n; i++) { int ear = WRAP(i+offset, l.n); if(l.elem[ear].ear == SPoint::EAR) { - if(!srf) { + if(srf->degm == 1 && srf->degn == 1) { + // This is a plane; any ear is a good ear. bestEar = ear; break; } @@ -322,7 +331,7 @@ void SContour::UvTriangulateInto(SMesh *m, SSurface *srf) { Vector prev = l.elem[WRAP((i+offset-1), l.n)].p, next = l.elem[WRAP((i+offset+1), l.n)].p; double tol = srf->ChordToleranceForEdge(prev, next); - if(tol < bestChordTol - LENGTH_EPS) { + if(tol < bestChordTol - scaledEps) { bestEar = ear; bestChordTol = tol; } @@ -335,10 +344,10 @@ void SContour::UvTriangulateInto(SMesh *m, SSurface *srf) { dbp("couldn't find an ear! fail"); return; } - ClipEarInto(m, bestEar); + ClipEarInto(m, bestEar, scaledEps); } - ClipEarInto(m, 0); // add the last triangle + ClipEarInto(m, 0, scaledEps); // add the last triangle } double SSurface::ChordToleranceForEdge(Vector a, Vector b) {