From ae35b3595cf646be78ba6636baa2ad01a5cf118d Mon Sep 17 00:00:00 2001 From: Jonathan Westhues Date: Wed, 3 Jun 2009 19:59:40 -0800 Subject: [PATCH] Revamp the edge classification for Booleans. I no longer make a separate polygon of coincident (with same or opposite normal) faces; I instead test all the edges against the other shell, and have extended the classify-against-shell stuff to handle those cases. And the normals are now perturbed a bit numerically, to either side of the edge, to distinguish tangency from a coincident surface. This seems to work fairly well, although things still tend to fail when the piecewise linear tolerance is too coarse. [git-p4: depot-paths = "//depot/solvespace/": change = 1964] --- graphicswin.cpp | 2 +- solvespace.cpp | 7 ++ srf/boolean.cpp | 251 ++++++++++++++++++++++++---------------------- srf/ratpoly.cpp | 59 ++++++++++- srf/surface.h | 48 +++++---- srf/surfinter.cpp | 197 +++++++++++++++++++++++++++++------- util.cpp | 5 +- wishlist.txt | 3 +- 8 files changed, 392 insertions(+), 180 deletions(-) diff --git a/graphicswin.cpp b/graphicswin.cpp index d937f397..97e32748 100644 --- a/graphicswin.cpp +++ b/graphicswin.cpp @@ -93,7 +93,7 @@ const GraphicsWindow::MenuEntry GraphicsWindow::menu[] = { { 1, "&Vertical\tV", MNU_VERTICAL, 'V', mCon }, { 1, NULL, 0, NULL }, { 1, "&On Point / Curve / Plane\tO", MNU_ON_ENTITY, 'O', mCon }, -{ 1, "E&qual Length / Radius\tQ", MNU_EQUAL, 'Q', mCon }, +{ 1, "E&qual Length / Radius / Angle\tQ", MNU_EQUAL, 'Q', mCon }, { 1, "Length Ra&tio\tZ", MNU_RATIO, 'Z', mCon }, { 1, "At &Midpoint\tM", MNU_AT_MIDPOINT, 'M', mCon }, { 1, "S&ymmetric\tY", MNU_SYMMETRIC, 'Y', mCon }, diff --git a/solvespace.cpp b/solvespace.cpp index 2a0025af..59e2542a 100644 --- a/solvespace.cpp +++ b/solvespace.cpp @@ -199,6 +199,13 @@ void SolveSpace::AfterNewFile(void) { // and the naked edges nakedEdges.Clear(); + // GenerateAll() expects the view to be valid, because it uses that to + // fill in default values for extrusion depths etc. (which won't matter + // here, but just don't let it work on garbage) + SS.GW.offset = Vector::From(0, 0, 0); + SS.GW.projRight = Vector::From(1, 0, 0); + SS.GW.projUp = Vector::From(0, 1, 0); + ReloadAllImported(); GenerateAll(-1, -1); diff --git a/srf/boolean.cpp b/srf/boolean.cpp index 8b720ec8..a487aa4a 100644 --- a/srf/boolean.cpp +++ b/srf/boolean.cpp @@ -158,26 +158,12 @@ void SSurface::TrimFromEdgeList(SEdgeList *el) { } } -// For each edge, we record the membership of the regions to its left and -// right, which we call the "in direction" and "out direction" (wrt its -// outer normal) -#define INDIR (0) -#define OUTDIR (8) -// Regions of interest are the other shell itself, the coincident faces of the -// shell (same or opposite normal) and the original surface. -#define SHELL (0) -#define COINCIDENT_SAME (1) -#define COINCIDENT_OPPOSITE (2) -#define ORIG (3) -// Macro for building bit to test -#define INSIDE(reg, dir) (1 << ((reg)+(dir))) - -static bool KeepRegion(int type, bool opA, int tag, int dir) +static bool KeepRegion(int type, bool opA, int shell, int orig) { - bool inShell = (tag & INSIDE(SHELL, dir)) != 0, - inSame = (tag & INSIDE(COINCIDENT_SAME, dir)) != 0, - inOpp = (tag & INSIDE(COINCIDENT_OPPOSITE, dir)) != 0, - inOrig = (tag & INSIDE(ORIG, dir)) != 0; + bool inShell = (shell == SShell::INSIDE), + inSame = (shell == SShell::COINC_SAME), + inOpp = (shell == SShell::COINC_OPP), + inOrig = (orig == SShell::INSIDE); bool inFace = inSame || inOpp; @@ -204,10 +190,12 @@ static bool KeepRegion(int type, bool opA, int tag, int dir) default: oops(); } } -static bool KeepEdge(int type, bool opA, int tag) +static bool KeepEdge(int type, bool opA, + int indir_shell, int outdir_shell, + int indir_orig, int outdir_orig) { - bool keepIn = KeepRegion(type, opA, tag, INDIR), - keepOut = KeepRegion(type, opA, tag, OUTDIR); + bool keepIn = KeepRegion(type, opA, indir_shell, indir_orig), + keepOut = KeepRegion(type, opA, outdir_shell, outdir_orig); // If the regions to the left and right of this edge are both in or both // out, then this edge is not useful and should be discarded. @@ -215,37 +203,33 @@ static bool KeepEdge(int type, bool opA, int tag) return false; } -static int TagByClassifiedEdge(int bspclass, int reg) +static void TagByClassifiedEdge(int bspclass, int *indir, int *outdir) { switch(bspclass) { case SBspUv::INSIDE: - return INSIDE(reg, OUTDIR) | INSIDE(reg, INDIR); + *indir = SShell::INSIDE; + *outdir = SShell::INSIDE; + break; case SBspUv::OUTSIDE: - return 0; + *indir = SShell::OUTSIDE; + *outdir = SShell::OUTSIDE; + break; case SBspUv::EDGE_PARALLEL: - return INSIDE(reg, INDIR); + *indir = SShell::INSIDE; + *outdir = SShell::OUTSIDE; + break; case SBspUv::EDGE_ANTIPARALLEL: - return INSIDE(reg, OUTDIR); + *indir = SShell::OUTSIDE; + *outdir = SShell::INSIDE; + break; default: oops(); } } -void DBPEDGE(int tag) { - dbp("edge: indir %s %s %s %s ; outdir %s %s %s %s", - (tag & INSIDE(SHELL, INDIR)) ? "shell" : "", - (tag & INSIDE(COINCIDENT_SAME, INDIR)) ? "coinc-same" : "", - (tag & INSIDE(COINCIDENT_OPPOSITE, INDIR)) ? "coinc-opp" : "", - (tag & INSIDE(ORIG, INDIR)) ? "orig" : "", - (tag & INSIDE(SHELL, OUTDIR)) ? "shell" : "", - (tag & INSIDE(COINCIDENT_SAME, OUTDIR)) ? "coinc-same" : "", - (tag & INSIDE(COINCIDENT_OPPOSITE, OUTDIR)) ? "coinc-opp" : "", - (tag & INSIDE(ORIG, OUTDIR)) ? "orig" : ""); -} - void DEBUGEDGELIST(SEdgeList *sel, SSurface *surf) { dbp("print %d edges", sel->l.n); SEdge *se; @@ -264,6 +248,54 @@ void DEBUGEDGELIST(SEdgeList *sel, SSurface *surf) { } } +static char *REGION(int d) { + switch(d) { + case SShell::INSIDE: return "inside"; + case SShell::OUTSIDE: return "outside"; + case SShell::COINC_SAME: return "same"; + case SShell::COINC_OPP: return "opp"; + default: return "xxx"; + } +} + + +void SSurface::EdgeNormalsWithinSurface(Point2d auv, Point2d buv, + Vector *pt, + Vector *enin, Vector *enout, + Vector *surfn, + DWORD auxA, SShell *shell) +{ + // the midpoint of the edge + Point2d muv = (auv.Plus(buv)).ScaledBy(0.5); + // a vector parallel to the edge + Point2d abuv = buv.Minus(auv).WithMagnitude(0.01); + // the edge's inner normal + Point2d enuv = abuv.Normal(); + + *pt = PointAt(muv); + *surfn = NormalAt(muv.x, muv.y); + + // If this edge just approximates a curve, then refine our midpoint so + // so that it actually lies on that curve too. Otherwise stuff like + // point-on-face tests will fail, since the point won't actually lie + // on the other face. + hSCurve hc = { auxA }; + SCurve *sc = shell->curve.FindById(hc); + if(sc->isExact && sc->exact.deg != 1) { + double t; + sc->exact.ClosestPointTo(*pt, &t, false); + *pt = sc->exact.PointAt(t); + } + + // Compute the inner and outer normals of this edge (within the srf), + // in xyz space. These are not necessarily antiparallel, if the + // surface is curved. + Vector pin = PointAt(muv.Plus(enuv)), + pout = PointAt(muv.Minus(enuv)); + *enin = pin.Minus(*pt), + *enout = pout.Minus(*pt); +} + //----------------------------------------------------------------------------- // Trim this surface against the specified shell, in the way that's appropriate // for the specified Boolean operation type (and which operand we are). We @@ -294,7 +326,8 @@ SSurface SSurface::MakeCopyTrimAgainst(SShell *agnst, SShell *parent, } // Build up our original trim polygon; remember the coordinates could - // be changed if we just flipped the surface normal. + // be changed if we just flipped the surface normal, and we are using + // the split curves (not the original curves). SEdgeList orig; ZERO(&orig); ret.MakeEdgesInto(into, &orig, true); @@ -302,22 +335,6 @@ SSurface SSurface::MakeCopyTrimAgainst(SShell *agnst, SShell *parent, // which means that we can't necessarily use the old BSP... SBspUv *origBsp = SBspUv::From(&orig); - // Find any surfaces from the other shell that are coincident with ours, - // and get the intersection polygons, grouping surfaces with the same - // and with opposite normal separately. - SEdgeList sameNormal, oppositeNormal; - ZERO(&sameNormal); - ZERO(&oppositeNormal); - agnst->MakeCoincidentEdgesInto(&ret, true, &sameNormal, into); - agnst->MakeCoincidentEdgesInto(&ret, false, &oppositeNormal, into); - // and cull parallel or anti-parallel pairs, which may occur if multiple - // surfaces are coincident with ours - sameNormal.CullExtraneousEdges(); - oppositeNormal.CullExtraneousEdges(); - // and build the trees for quick in-polygon testing - SBspUv *sameBsp = SBspUv::From(&sameNormal); - SBspUv *oppositeBsp = SBspUv::From(&oppositeNormal); - // And now intersect the other shell against us SEdgeList inter; ZERO(&inter); @@ -373,45 +390,25 @@ SSurface SSurface::MakeCopyTrimAgainst(SShell *agnst, SShell *parent, SEdge *se; for(se = orig.l.First(); se; se = orig.l.NextAfter(se)) { - Point2d auv = (se->a).ProjectXy(), - buv = (se->b).ProjectXy(); + Point2d auv = (se->a).ProjectXy(), + buv = (se->b).ProjectXy(); - int c_same = sameBsp->ClassifyEdge(auv, buv); - int c_opp = oppositeBsp->ClassifyEdge(auv, buv); + Vector pt, enin, enout, surfn; + ret.EdgeNormalsWithinSurface(auv, buv, &pt, &enin, &enout, &surfn, + se->auxA, into); - // Get the midpoint of this edge - Point2d am = (auv.Plus(buv)).ScaledBy(0.5); - Vector pt = ret.PointAt(am.x, am.y); - // and the outer normal from the trim polygon (within the surface) - Vector surf_n = ret.NormalAt(am.x, am.y); - Vector ea = ret.PointAt(auv.x, auv.y), - eb = ret.PointAt(buv.x, buv.y); - Vector edge_n = surf_n.Cross((eb.Minus(ea))); + int indir_shell, outdir_shell, indir_orig, outdir_orig; - int c_shell = agnst->ClassifyPoint(pt, edge_n, surf_n); + indir_orig = SShell::INSIDE; + outdir_orig = SShell::OUTSIDE; - int tag = 0; - tag |= INSIDE(ORIG, INDIR); - tag |= TagByClassifiedEdge(c_same, COINCIDENT_SAME); - tag |= TagByClassifiedEdge(c_opp, COINCIDENT_OPPOSITE); + agnst->ClassifyEdge(&indir_shell, &outdir_shell, + ret.PointAt(auv), ret.PointAt(buv), pt, + enin, enout, surfn); - if(c_shell == SShell::INSIDE) { - tag |= INSIDE(SHELL, INDIR) | INSIDE(SHELL, OUTDIR); - } else if(c_shell == SShell::OUTSIDE) { - tag |= 0; - } else if(c_shell == SShell::SURF_PARALLEL) { - tag |= INSIDE(SHELL, INDIR); - } else if(c_shell == SShell::SURF_ANTIPARALLEL) { - tag |= INSIDE(SHELL, OUTDIR); - } else if(c_shell == SShell::EDGE_PARALLEL) { - tag |= INSIDE(SHELL, INDIR); - } else if(c_shell == SShell::EDGE_ANTIPARALLEL) { - tag |= INSIDE(SHELL, OUTDIR); - } else if(c_shell == SShell::EDGE_TANGENT) { - continue; - } - - if(KeepEdge(type, opA, tag)) { + if(KeepEdge(type, opA, indir_shell, outdir_shell, + indir_orig, outdir_orig)) + { final.AddEdge(se->a, se->b, se->auxA, se->auxB); } } @@ -420,23 +417,22 @@ SSurface SSurface::MakeCopyTrimAgainst(SShell *agnst, SShell *parent, Point2d auv = (se->a).ProjectXy(), buv = (se->b).ProjectXy(); + Vector pt, enin, enout, surfn; + ret.EdgeNormalsWithinSurface(auv, buv, &pt, &enin, &enout, &surfn, + se->auxA, into); + + int indir_shell, outdir_shell, indir_orig, outdir_orig; + int c_this = origBsp->ClassifyEdge(auv, buv); - int c_same = sameBsp->ClassifyEdge(auv, buv); - int c_opp = oppositeBsp->ClassifyEdge(auv, buv); + TagByClassifiedEdge(c_this, &indir_orig, &outdir_orig); - int tag = 0; - tag |= TagByClassifiedEdge(c_this, ORIG); - tag |= TagByClassifiedEdge(c_same, COINCIDENT_SAME); - tag |= TagByClassifiedEdge(c_opp, COINCIDENT_OPPOSITE); + agnst->ClassifyEdge(&indir_shell, &outdir_shell, + ret.PointAt(auv), ret.PointAt(buv), pt, + enin, enout, surfn); - if(type == SShell::AS_DIFFERENCE && !opA) { - // The second operand of a difference gets turned inside out - tag |= INSIDE(SHELL, INDIR); - } else { - tag |= INSIDE(SHELL, OUTDIR); - } - - if(KeepEdge(type, opA, tag)) { + if(KeepEdge(type, opA, indir_shell, outdir_shell, + indir_orig, outdir_orig)) + { final.AddEdge(se->a, se->b, se->auxA, se->auxB); } } @@ -455,12 +451,11 @@ SSurface SSurface::MakeCopyTrimAgainst(SShell *agnst, SShell *parent, final.l.ClearTags(); if(!final.AssemblePolygon(&poly, NULL, true)) { into->booleanFailed = true; + dbp("failed: I=%d", I); DEBUGEDGELIST(&final, &ret); } poly.Clear(); - sameNormal.Clear(); - oppositeNormal.Clear(); final.Clear(); inter.Clear(); orig.Clear(); @@ -496,6 +491,7 @@ void SShell::MakeIntersectionCurvesAgainst(SShell *agnst, SShell *into) { void SShell::CleanupAfterBoolean(void) { SSurface *ss; for(ss = surface.First(); ss; ss = surface.NextAfter(ss)) { + ss->edges.Clear(); } } @@ -576,8 +572,8 @@ void SShell::MakeFromAssemblyOf(SShell *a, SShell *b) { void SShell::MakeFromBoolean(SShell *a, SShell *b, int type) { booleanFailed = false; - a->MakeClassifyingBsps(); - b->MakeClassifyingBsps(); + a->MakeClassifyingBsps(NULL); + b->MakeClassifyingBsps(NULL); // Copy over all the original curves, splitting them so that a // piecwise linear segment never crosses a surface from the other @@ -605,11 +601,21 @@ void SShell::MakeFromBoolean(SShell *a, SShell *b, int type) { sc->RemoveShortSegments(srfA, srfB); } + // And clean up the piecewise linear things we made as a calculation aid + a->CleanupAfterBoolean(); + b->CleanupAfterBoolean(); + // Remake the classifying BSPs with the split (and short-segment-removed) + // curves + a->MakeClassifyingBsps(this); + b->MakeClassifyingBsps(this); + + // Then trim and copy the surfaces if(b->surface.n == 0 || a->surface.n == 0) { - // Then trim and copy the surfaces + I = 1023123; a->CopySurfacesTrimAgainst(b, this, type, true); b->CopySurfacesTrimAgainst(a, this, type, false); } else { + I = 0; a->CopySurfacesTrimAgainst(b, this, type, true); b->CopySurfacesTrimAgainst(a, this, type, false); } @@ -627,21 +633,23 @@ void SShell::MakeFromBoolean(SShell *a, SShell *b, int type) { //----------------------------------------------------------------------------- // All of the BSP routines that we use to perform and accelerate polygon ops. //----------------------------------------------------------------------------- -void SShell::MakeClassifyingBsps(void) { +void SShell::MakeClassifyingBsps(SShell *useCurvesFrom) { SSurface *ss; for(ss = surface.First(); ss; ss = surface.NextAfter(ss)) { - ss->MakeClassifyingBsp(this); + ss->MakeClassifyingBsp(this, useCurvesFrom); } } -void SSurface::MakeClassifyingBsp(SShell *shell) { +void SSurface::MakeClassifyingBsp(SShell *shell, SShell *useCurvesFrom) { SEdgeList el; ZERO(&el); - MakeEdgesInto(shell, &el, true); + MakeEdgesInto(shell, &el, true, useCurvesFrom); bsp = SBspUv::From(&el); - el.Clear(); + + ZERO(&edges); + MakeEdgesInto(shell, &edges, false, useCurvesFrom); } SBspUv *SBspUv::Alloc(void) { @@ -733,7 +741,7 @@ SBspUv *SBspUv::InsertEdge(Point2d ea, Point2d eb) { return this; } -int SBspUv::ClassifyPoint(Point2d p, Point2d eb) { +int SBspUv::ClassifyPoint(Point2d p, Point2d eb, Point2d *ia, Point2d *ib) { if(!this) return OUTSIDE; Point2d n = ((b.Minus(a)).Normal()).WithMagnitude(1); @@ -746,6 +754,11 @@ int SBspUv::ClassifyPoint(Point2d p, Point2d eb) { while(f) { Point2d ba = (f->b).Minus(f->a); if(p.DistanceToLine(f->a, ba, true) < LENGTH_EPS) { + if(ia && ib) { + // 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; @@ -759,16 +772,16 @@ int SBspUv::ClassifyPoint(Point2d p, Point2d eb) { f = f->more; } // Pick arbitrarily which side to send it down, doesn't matter - int c1 = neg ? neg->ClassifyPoint(p, eb) : OUTSIDE; - int c2 = pos ? pos->ClassifyPoint(p, eb) : INSIDE; + int c1 = neg ? neg->ClassifyPoint(p, eb, ia, ib) : OUTSIDE; + int c2 = pos ? pos->ClassifyPoint(p, eb, ia, ib) : 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) : INSIDE; + return pos ? pos->ClassifyPoint(p, eb, ia, ib) : INSIDE; } else { - return neg ? neg->ClassifyPoint(p, eb) : OUTSIDE; + return neg ? neg->ClassifyPoint(p, eb, ia, ib) : OUTSIDE; } } diff --git a/srf/ratpoly.cpp b/srf/ratpoly.cpp index 28e7ff94..17d57462 100644 --- a/srf/ratpoly.cpp +++ b/srf/ratpoly.cpp @@ -130,7 +130,7 @@ Vector SBezier::TangentAt(double t) { return ret; } -void SBezier::ClosestPointTo(Vector p, double *t) { +void SBezier::ClosestPointTo(Vector p, double *t, bool converge) { int i; double minDist = VERY_POSITIVE; *t = 0; @@ -147,7 +147,7 @@ void SBezier::ClosestPointTo(Vector p, double *t) { } Vector p0; - for(i = 0; i < 15; i++) { + for(i = 0; i < (converge ? 15 : 3); i++) { p0 = PointAt(*t); if(p0.Equals(p, RATPOLY_EPS)) { return; @@ -157,7 +157,9 @@ void SBezier::ClosestPointTo(Vector p, double *t) { Vector pc = p.ClosestPointOnLine(p0, dp); *t += (pc.Minus(p0)).DivPivoting(dp); } - dbp("didn't converge (closest point on bezier curve)"); + if(converge) { + dbp("didn't converge (closest point on bezier curve)"); + } } void SBezier::SplitAt(double t, SBezier *bef, SBezier *aft) { @@ -242,6 +244,9 @@ void SBezier::MakePwlWorker(List *l, double ta, double tb) { } } +Vector SSurface::PointAt(Point2d puv) { + return PointAt(puv.x, puv.y); +} Vector SSurface::PointAt(double u, double v) { Vector num = Vector::From(0, 0, 0); double den = 0; @@ -294,12 +299,18 @@ void SSurface::TangentsAt(double u, double v, Vector *tu, Vector *tv) { *tv = tv->ScaledBy(1.0/(den*den)); } +Vector SSurface::NormalAt(Point2d puv) { + return NormalAt(puv.x, puv.y); +} Vector SSurface::NormalAt(double u, double v) { Vector tu, tv; TangentsAt(u, v, &tu, &tv); return tu.Cross(tv); } +void SSurface::ClosestPointTo(Vector p, Point2d *puv, bool converge) { + ClosestPointTo(p, &(puv->x), &(puv->y), converge); +} void SSurface::ClosestPointTo(Vector p, double *u, double *v, bool converge) { // A few special cases first; when control points are coincident the // derivative goes to zero at the conrol points, and would result in @@ -394,6 +405,48 @@ bool SSurface::PointIntersectingLine(Vector p0, Vector p1, double *u, double *v) return false; } +Vector SSurface::ClosestPointOnThisAndSurface(SSurface *srf2, Vector p) { + // This is untested. + int i, j; + Point2d puv[2]; + SSurface *srf[2] = { this, srf2 }; + + for(j = 0; j < 2; j++) { + (srf[j])->ClosestPointTo(p, &(puv[j]), false); + } + + for(i = 0; i < 4; i++) { + Vector tu[2], tv[2], cp[2], n[2]; + double d[2]; + + for(j = 0; j < 2; j++) { + (srf[j])->TangentsAt(puv[j].x, puv[j].y, &(tu[j]), &(tv[j])); + + cp[j] = (srf[j])->PointAt(puv[j]); + + n[j] = ((tu[j]).Cross(tv[j])).WithMagnitude(1); + d[j] = (n[j]).Dot(cp[j]); + } + + Vector p0 = Vector::AtIntersectionOfPlanes(n[0], d[0], n[1], d[1]), + dp = (n[0]).Cross(n[1]); + + Vector pc = p.ClosestPointOnLine(p0, dp); + + // Adjust our guess and iterate + for(j = 0; j < 2; j++) { + Vector dc = pc.Minus(cp[j]); + double du = dc.Dot(tu[j]), dv = dc.Dot(tv[j]); + puv[j].x += du / ((tu[j]).MagSquared()); + puv[j].y += dv / ((tv[j]).MagSquared()); + } + } + + // If this converged, then the two points are actually equal. + return ((srf[0])->PointAt(puv[0])).Plus( + ((srf[1])->PointAt(puv[1]))).ScaledBy(0.5); +} + void SSurface::PointOnSurfaces(SSurface *s1, SSurface *s2, double *up, double *vp) { diff --git a/srf/surface.h b/srf/surface.h index c8bb888b..e5a2555b 100644 --- a/srf/surface.h +++ b/srf/surface.h @@ -31,7 +31,8 @@ public: Point2d IntersectionWith(Point2d a, Point2d b); SBspUv *InsertEdge(Point2d a, Point2d b); - int ClassifyPoint(Point2d p, Point2d eb); + int ClassifyPoint(Point2d p, Point2d eb, + Point2d *ia=NULL, Point2d *ib=NULL); int ClassifyEdge(Point2d ea, Point2d eb); }; @@ -62,7 +63,7 @@ public: Vector PointAt(double t); Vector TangentAt(double t); - void ClosestPointTo(Vector p, double *t); + void ClosestPointTo(Vector p, double *t, bool converge=true); void SplitAt(double t, SBezier *bef, SBezier *aft); Vector Start(void); @@ -184,9 +185,10 @@ public: int tag; Vector p; SSurface *srf; - hSSurface hsrf; - Vector surfNormal; // of the intersecting surface, at pinter - bool onEdge; // pinter is on edge of trim poly + 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. @@ -209,6 +211,7 @@ public: // For testing whether a point (u, v) on the surface lies inside the trim SBspUv *bsp; + SEdgeList edges; static SSurface FromExtrusionOf(SBezier *spc, Vector t0, Vector t1); static SSurface FromRevolutionOf(SBezier *sb, Vector pt, Vector axis, @@ -217,6 +220,10 @@ public: static SSurface FromTransformationOf(SSurface *a, Vector t, Quaternion q, bool includingTrims); + void EdgeNormalsWithinSurface(Point2d auv, Point2d buv, + Vector *pt, Vector *enin, Vector *enout, + Vector *surfn, + DWORD auxA, SShell *shell); SSurface MakeCopyTrimAgainst(SShell *against, SShell *parent, SShell *into, int type, bool opA); void TrimFromEdgeList(SEdgeList *el); @@ -244,11 +251,15 @@ public: List *l, bool segment, SSurface *sorig); + void ClosestPointTo(Vector p, Point2d *puv, bool converge=true); void ClosestPointTo(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); Vector PointAt(double u, double v); + Vector PointAt(Point2d puv); void TangentsAt(double u, double v, Vector *tu, Vector *tv); + Vector NormalAt(Point2d puv); Vector NormalAt(double u, double v); bool LineEntirelyOutsideBbox(Vector a, Vector b, bool segment); void GetAxisAlignedBounding(Vector *ptMax, Vector *ptMin); @@ -263,7 +274,7 @@ public: void MakeEdgesInto(SShell *shell, SEdgeList *sel, bool asUv, SShell *useCurvesFrom=NULL); void MakeSectionEdgesInto(SShell *shell, SEdgeList *sel, SBezierList *sbl); - void MakeClassifyingBsp(SShell *shell); + void MakeClassifyingBsp(SShell *shell, SShell *useCurvesFrom); double ChordToleranceForEdge(Vector a, Vector b); void MakeTriangulationGridInto(List *l, double vs, double vf, bool swapped); @@ -294,7 +305,7 @@ public: void CopyCurvesSplitAgainst(bool opA, SShell *agnst, SShell *into); void CopySurfacesTrimAgainst(SShell *against, SShell *into, int t, bool a); void MakeIntersectionCurvesAgainst(SShell *against, SShell *into); - void MakeClassifyingBsps(void); + void MakeClassifyingBsps(SShell *useCurvesFrom); void AllPointsIntersecting(Vector a, Vector b, List *il, bool seg, bool trimmed, bool inclTangent); void MakeCoincidentEdgesInto(SSurface *proto, bool sameNormal, @@ -302,16 +313,19 @@ public: void RewriteSurfaceHandlesForCurves(SShell *a, SShell *b); void CleanupAfterBoolean(void); - static const int INSIDE = 100; - static const int OUTSIDE = 200; - static const int SURF_PARALLEL = 300; - static const int SURF_ANTIPARALLEL = 400; - static const int EDGE_PARALLEL = 500; - static const int EDGE_ANTIPARALLEL = 600; - static const int EDGE_TANGENT = 700; - - int ClassifyPoint(Vector p, Vector edge_n, Vector surf_n); - + // Definitions when classifying regions of a surface; it is either inside, + // outside, or coincident (with parallel or antiparallel normal) with a + // shell. + static const int INSIDE = 100; + static const int OUTSIDE = 200; + static const int COINC_SAME = 300; + static const int COINC_OPP = 400; + static const double DOTP_TOL; + int ClassifyRegion(Vector edge_n, Vector inter_surf_n, Vector edge_surf_n); + bool ClassifyEdge(int *indir, int *outdir, + Vector ea, Vector eb, + Vector p, + Vector edge_n_in, Vector edge_n_out, Vector surf_n); void MakeFromCopyOf(SShell *a); void MakeFromTransformationOf(SShell *a, Vector trans, Quaternion q); diff --git a/srf/surfinter.cpp b/srf/surfinter.cpp index 961a6d70..b1666595 100644 --- a/srf/surfinter.cpp +++ b/srf/surfinter.cpp @@ -1,5 +1,8 @@ #include "solvespace.h" +// Dot product tolerance for perpendicular. +const double SShell::DOTP_TOL = 1e-3; + extern int FLAG; void SSurface::AddExactIntersectionCurve(SBezier *sb, SSurface *srfB, @@ -220,9 +223,25 @@ void SSurface::IntersectAgainst(SSurface *b, SShell *agnstA, SShell *agnstB, List lv; ZERO(&lv); - Vector axisa = axis.ScaledBy((b->ctrl[0][0]).Dot(axis)), - axisb = axis.ScaledBy((b->ctrl[0][1]).Dot(axis)), - axisc = (axisa.Plus(axisb)).ScaledBy(0.5); + double a_axis0 = ( ctrl[0][0]).Dot(axis), + a_axis1 = ( ctrl[0][1]).Dot(axis), + b_axis0 = (b->ctrl[0][0]).Dot(axis), + b_axis1 = (b->ctrl[0][1]).Dot(axis); + + if(a_axis0 > a_axis1) SWAP(double, a_axis0, a_axis1); + if(b_axis0 > b_axis1) SWAP(double, b_axis0, b_axis1); + + double ab_axis0 = max(a_axis0, b_axis0), + ab_axis1 = min(a_axis1, b_axis1); + + if(fabs(ab_axis0 - ab_axis1) < LENGTH_EPS) { + // The line would be zero-length + return; + } + + Vector axis0 = axis.ScaledBy(ab_axis0), + axis1 = axis.ScaledBy(ab_axis1), + axisc = (axis0.Plus(axis1)).ScaledBy(0.5); oft.MakePwlInto(&lv); @@ -250,7 +269,7 @@ void SSurface::IntersectAgainst(SSurface *b, SShell *agnstA, SShell *agnstB, p = b->PointAt(ub, vb); SBezier bezier; - bezier = SBezier::From(p.Plus(axisa), p.Plus(axisb)); + bezier = SBezier::From(p.Plus(axis0), p.Plus(axis1)); AddExactIntersectionCurve(&bezier, b, agnstA, agnstB, into); } @@ -608,8 +627,8 @@ void SSurface::AllPointsIntersecting(Vector a, Vector b, } // And that it lies inside our trim region - Point2d dummy = { 0, 0 }; - int c = bsp->ClassifyPoint(puv, dummy); + Point2d dummy = { 0, 0 }, ia = { 0, 0 }, ib = { 0, 0 }; + int c = bsp->ClassifyPoint(puv, dummy, &ia, &ib); if(trimmed && c == SBspUv::OUTSIDE) { continue; } @@ -618,9 +637,11 @@ void SSurface::AllPointsIntersecting(Vector a, Vector b, SInter si; si.p = pxyz; si.surfNormal = NormalAt(puv.x, puv.y); - si.hsrf = h; + si.pinter = puv; si.srf = this; si.onEdge = (c != SBspUv::INSIDE); + si.edgeA = ia; + si.edgeB = ib; l->Add(&si); } @@ -637,6 +658,26 @@ void SShell::AllPointsIntersecting(Vector a, Vector b, } } +int SShell::ClassifyRegion(Vector edge_n, Vector inter_surf_n, + Vector edge_surf_n) +{ + double dot = inter_surf_n.Dot(edge_n); + if(fabs(dot) < DOTP_TOL) { + // The edge's surface and the edge-on-face surface + // are coincident. Test the edge's surface normal + // to see if it's with same or opposite normals. + if(inter_surf_n.Dot(edge_surf_n) > 0) { + return COINC_SAME; + } else { + return COINC_OPP; + } + } else if(dot > 0) { + return OUTSIDE; + } else { + return INSIDE; + } +} + //----------------------------------------------------------------------------- // Does the given point lie on our shell? There are many cases; inside and // outside are obvious, but then there's all the edge-on-edge and edge-on-face @@ -646,15 +687,101 @@ void SShell::AllPointsIntersecting(Vector a, Vector b, // using the closest intersection point. If the ray hits a surface on edge, // then just reattempt in a different random direction. //----------------------------------------------------------------------------- -int SShell::ClassifyPoint(Vector p, Vector edge_n, Vector surf_n) { +bool SShell::ClassifyEdge(int *indir, int *outdir, + Vector ea, Vector eb, + Vector p, + Vector edge_n_in, Vector edge_n_out, Vector surf_n) +{ List l; ZERO(&l); srand(0); - int ret, cnt = 0, edge_inters; - double edge_dotp[2]; + // First, check for edge-on-edge + int edge_inters = 0; + Vector inter_surf_n[2], inter_edge_n[2]; + SSurface *srf; + for(srf = surface.First(); srf; srf = surface.NextAfter(srf)) { + if(srf->LineEntirelyOutsideBbox(ea, eb, true)) continue; + SEdgeList *sel = &(srf->edges); + SEdge *se; + for(se = sel->l.First(); se; se = sel->l.NextAfter(se)) { + if((ea.Equals(se->a) && eb.Equals(se->b)) || + (eb.Equals(se->a) && ea.Equals(se->b)) || + p.OnLineSegment(se->a, se->b)) + { + if(edge_inters < 2) { + // Edge-on-edge case + Point2d pm; + srf->ClosestPointTo(p, &pm, false); + // A vector normal to the surface, at the intersection point + inter_surf_n[edge_inters] = srf->NormalAt(pm); + // A vector normal to the intersecting edge (but within the + // intersecting surface) at the intersection point, pointing + // out. + inter_edge_n[edge_inters] = + (inter_surf_n[edge_inters]).Cross((se->b).Minus((se->a))); + } + + edge_inters++; + } + } + } + + if(edge_inters == 2) { + // TODO, make this use the appropriate curved normals + double dotp[2]; + for(int i = 0; i < 2; i++) { + dotp[i] = edge_n_out.Dot(inter_surf_n[i]); + } + + if(fabs(dotp[1]) < DOTP_TOL) { + SWAP(double, dotp[0], dotp[1]); + SWAP(Vector, inter_surf_n[0], inter_surf_n[1]); + SWAP(Vector, inter_edge_n[0], inter_edge_n[1]); + } + + int coinc = (surf_n.Dot(inter_surf_n[0])) > 0 ? COINC_SAME : COINC_OPP; + + if(fabs(dotp[0]) < DOTP_TOL && fabs(dotp[1]) < DOTP_TOL) { + // This is actually an edge on face case, just that the face + // is split into two pieces joining at our edge. + *indir = coinc; + *outdir = coinc; + } else if(fabs(dotp[0]) < DOTP_TOL && dotp[1] > DOTP_TOL) { + if(edge_n_out.Dot(inter_edge_n[0]) > 0) { + *indir = coinc; + *outdir = OUTSIDE; + } else { + *indir = INSIDE; + *outdir = coinc; + } + } else if(fabs(dotp[0]) < DOTP_TOL && dotp[1] < -DOTP_TOL) { + if(edge_n_out.Dot(inter_edge_n[0]) > 0) { + *indir = coinc; + *outdir = INSIDE; + } else { + *indir = OUTSIDE; + *outdir = coinc; + } + } else if(dotp[0] > DOTP_TOL && dotp[1] > DOTP_TOL) { + *indir = INSIDE; + *outdir = OUTSIDE; + } else if(dotp[0] < -DOTP_TOL && dotp[1] < -DOTP_TOL) { + *indir = OUTSIDE; + *outdir = INSIDE; + } else { + // Edge is tangent to the shell at shell's edge, so can't be + // a boundary of the surface. + return false; + } + return true; + } + + if(edge_inters != 0) dbp("bad, edge_inters=%d", edge_inters); + + int cnt = 0; for(;;) { // Cast a ray in a random direction (two-sided so that we test if // the point lies on a surface, but use only one side for in/out @@ -663,8 +790,10 @@ int SShell::ClassifyPoint(Vector p, Vector edge_n, Vector surf_n) { AllPointsIntersecting( p.Minus(ray), p.Plus(ray), &l, false, true, false); + // no intersections means it's outside + *indir = OUTSIDE; + *outdir = OUTSIDE; double dmin = VERY_POSITIVE; - ret = OUTSIDE; // no intersections means it's outside bool onEdge = false; edge_inters = 0; @@ -678,9 +807,9 @@ int SShell::ClassifyPoint(Vector p, Vector edge_n, Vector surf_n) { double d = ((si->p).Minus(p)).Magnitude(); - // Handle edge-on-edge - if(d < LENGTH_EPS && si->onEdge && edge_inters < 2) { - edge_dotp[edge_inters] = (si->surfNormal).Dot(edge_n); + // We actually should never hit this case; it should have been + // handled above. + if(d < LENGTH_EPS && si->onEdge) { edge_inters++; } @@ -689,21 +818,24 @@ int SShell::ClassifyPoint(Vector p, Vector edge_n, Vector surf_n) { if(d < LENGTH_EPS) { // Edge-on-face (unless edge-on-edge above supercedes) - double dot = (si->surfNormal).Dot(edge_n); - if(fabs(dot) < LENGTH_EPS && 0) { - // TODO revamp this - } else if(dot > 0) { - ret = SURF_PARALLEL; - } else { - ret = SURF_ANTIPARALLEL; - } + Point2d pin, pout; + (si->srf)->ClosestPointTo(p.Plus(edge_n_in), &pin, false); + (si->srf)->ClosestPointTo(p.Plus(edge_n_out), &pout, false); + + Vector surf_n_in = (si->srf)->NormalAt(pin), + surf_n_out = (si->srf)->NormalAt(pout); + + *indir = ClassifyRegion(edge_n_in, surf_n_in, surf_n); + *outdir = ClassifyRegion(edge_n_out, surf_n_out, surf_n); } else { // Edge does not lie on surface; either strictly inside // or strictly outside if((si->surfNormal).Dot(ray) > 0) { - ret = INSIDE; + *indir = INSIDE; + *outdir = INSIDE; } else { - ret = OUTSIDE; + *indir = OUTSIDE; + *outdir = OUTSIDE; } } onEdge = si->onEdge; @@ -714,27 +846,16 @@ int SShell::ClassifyPoint(Vector p, Vector edge_n, Vector surf_n) { // If the point being tested lies exactly on an edge of the shell, // then our ray always lies on edge, and that's okay. Otherwise // try again in a different random direction. - if((edge_inters == 2) || !onEdge) break; + if(!onEdge) break; if(cnt++ > 5) { dbp("can't find a ray that doesn't hit on edge!"); dbp("on edge = %d, edge_inters = %d", onEdge, edge_inters); + SS.nakedEdges.AddEdge(ea, eb); break; } } - if(edge_inters == 2) { - double tol = 1e-3; - - if(edge_dotp[0] > -tol && edge_dotp[1] > -tol) { - return EDGE_PARALLEL; - } else if(edge_dotp[0] < tol && edge_dotp[1] < tol) { - return EDGE_ANTIPARALLEL; - } else { - return EDGE_TANGENT; - } - } else { - return ret; - } + return true; } //----------------------------------------------------------------------------- diff --git a/util.cpp b/util.cpp index 7ea23b7a..0b58aed5 100644 --- a/util.cpp +++ b/util.cpp @@ -460,6 +460,8 @@ double Vector::DistanceToLine(Vector p0, Vector dp) { } bool Vector::OnLineSegment(Vector a, Vector b) { + if(this->Equals(a) || this->Equals(b)) return true; + Vector d = b.Minus(a); double m = d.MagSquared(); @@ -468,7 +470,7 @@ bool Vector::OnLineSegment(Vector a, Vector b) { if(distsq >= LENGTH_EPS*LENGTH_EPS) return false; double t = (this->Minus(a)).DivPivoting(d); - // On-endpoint must be tested for separately. + // On-endpoint already tested if(t < 0 || t > 1) return false; return true; } @@ -508,6 +510,7 @@ Vector Vector::ScaledBy(double v) { Vector Vector::WithMagnitude(double v) { double m = Magnitude(); if(m == 0) { + dbp("Vector::WithMagnitude of zero vector!"); return From(v, 0, 0); } else { return ScaledBy(v/m); diff --git a/wishlist.txt b/wishlist.txt index 471b8d73..b403b80d 100644 --- a/wishlist.txt +++ b/wishlist.txt @@ -1,9 +1,10 @@ -marching algorithm for surface intersection tangent intersections +plane detection for solid of revolution ----- line styles (color, thickness) +marching algorithm for surface intersection loop detection IGES and STEP export incremental regen of entities?