diff --git a/drawentity.cpp b/drawentity.cpp index cc8f9be..9a4af1d 100644 --- a/drawentity.cpp +++ b/drawentity.cpp @@ -240,11 +240,11 @@ void Entity::GenerateBezierCurves(SBezierList *sbl) { ArcGetAngles(&thetaa, &thetab, &dtheta); } int i, n; - if(dtheta > 3*PI/2) { + if(dtheta > (3*PI/2 + 0.01)) { n = 4; - } else if(dtheta > PI) { + } else if(dtheta > (PI + 0.01)) { n = 3; - } else if(dtheta > PI/2) { + } else if(dtheta > (PI/2 + 0.01)) { n = 2; } else { n = 1; diff --git a/dsc.h b/dsc.h index 1ba704d..55a6d02 100644 --- a/dsc.h +++ b/dsc.h @@ -74,7 +74,7 @@ public: Vector DotInToCsys(Vector u, Vector v, Vector n); Vector ScaleOutOfCsys(Vector u, Vector v, Vector n); double DistanceToLine(Vector p0, Vector dp); - bool OnLineSegment(Vector a, Vector b); + bool OnLineSegment(Vector a, Vector b, double tol=LENGTH_EPS); Vector ClosestPointOnLine(Vector p0, Vector dp); double Magnitude(void); double MagSquared(void); diff --git a/polygon.h b/polygon.h index 34c50d6..f29efb6 100644 --- a/polygon.h +++ b/polygon.h @@ -40,6 +40,7 @@ public: int ear; Vector p; + Vector auxv; }; class SPointList { diff --git a/srf/ratpoly.cpp b/srf/ratpoly.cpp index b46b5e8..ac89d6c 100644 --- a/srf/ratpoly.cpp +++ b/srf/ratpoly.cpp @@ -147,7 +147,7 @@ void SBezier::ClosestPointTo(Vector p, double *t, bool converge) { } Vector p0; - for(i = 0; i < (converge ? 15 : 3); i++) { + for(i = 0; i < (converge ? 15 : 5); i++) { p0 = PointAt(*t); if(p0.Equals(p, RATPOLY_EPS)) { return; @@ -345,7 +345,7 @@ void SSurface::ClosestPointTo(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 ? 15 : 3); i++) { + for(i = 0; i < (converge ? 25 : 5); i++) { p0 = PointAt(*u, *v); if(converge) { if(p0.Equals(p, RATPOLY_EPS)) { @@ -443,6 +443,10 @@ Vector SSurface::ClosestPointOnThisAndSurface(SSurface *srf2, Vector p) { puv[j].y += dv / ((tv[j]).MagSquared()); } } + if(i >= 10) { + dbp("this and srf, didn't converge, d=%g", + (puv[0].Minus(puv[1])).Magnitude()); + } // If this converged, then the two points are actually equal. return ((srf[0])->PointAt(puv[0])).Plus( @@ -461,7 +465,7 @@ void SSurface::PointOnSurfaces(SSurface *s1, SSurface *s2, (srf[2])->ClosestPointTo(p, &(u[2]), &(v[2]), false); int i, j; - for(i = 0; i < 15; i++) { + for(i = 0; i < 20; i++) { // Approximate each surface by a plane Vector p[3], tu[3], tv[3], n[3]; double d[3]; diff --git a/srf/surfinter.cpp b/srf/surfinter.cpp index 68359bc..0fbe031 100644 --- a/srf/surfinter.cpp +++ b/srf/surfinter.cpp @@ -51,7 +51,7 @@ void SSurface::AddExactIntersectionCurve(SBezier *sb, SSurface *srfB, split = sc.MakeCopySplitAgainst(agnstA, agnstB, this, srfB); sc.Clear(); } - + if(0 && sb->deg == 1) { dbp(" "); SCurvePt *prev = NULL, *v; @@ -192,7 +192,6 @@ void SSurface::IntersectAgainst(SSurface *b, SShell *agnstA, SShell *agnstB, Vector al = along.ScaledBy(0.5); SBezier bezier; bezier = SBezier::From((si->p).Minus(al), (si->p).Plus(al)); - AddExactIntersectionCurve(&bezier, b, agnstA, agnstB, into); } @@ -313,10 +312,21 @@ void SSurface::IntersectAgainst(SSurface *b, SShell *agnstA, SShell *agnstB, for(si = lsi.First(); si; si = lsi.NextAfter(si)) { Vector p = si->p; double u, v; - srfA->ClosestPointTo(p, &u, &v); - srfA->PointOnSurfaces(srfB, other, &u, &v); - p = srfA->PointAt(u, v); - if(!spl.ContainsPoint(p)) spl.Add(p); + srfB->ClosestPointTo(p, &u, &v); + srfB->PointOnSurfaces(srfA, other, &u, &v); + p = srfB->PointAt(u, v); + if(!spl.ContainsPoint(p)) { + SPoint sp; + sp.p = p; + // We also need the edge normal, so that we know in + // which direction to march. + srfA->ClosestPointTo(p, &u, &v); + Vector n = srfA->NormalAt(u, v); + sp.auxv = n.Cross((se->b).Minus(se->a)); + sp.auxv = (sp.auxv).WithMagnitude(1); + + spl.l.Add(&sp); + } } lsi.Clear(); } @@ -324,9 +334,7 @@ void SSurface::IntersectAgainst(SSurface *b, SShell *agnstA, SShell *agnstB, el.Clear(); } - SPoint *sp; - if(spl.l.n == 2) { - + while(spl.l.n >= 2) { SCurve sc; ZERO(&sc); sc.surfA = h; @@ -334,16 +342,82 @@ void SSurface::IntersectAgainst(SSurface *b, SShell *agnstA, SShell *agnstB, sc.isExact = false; sc.source = SCurve::FROM_INTERSECTION; - SCurvePt scpt; - scpt.p = (spl.l.elem[0].p); - sc.pts.Add(&scpt); - scpt.p = (spl.l.elem[1].p); - sc.pts.Add(&scpt); + Vector start = spl.l.elem[0].p, + startv = spl.l.elem[0].auxv; + spl.l.ClearTags(); + spl.l.elem[0].tag = 1; + spl.l.RemoveTagged(); + SCurvePt padd; + ZERO(&padd); + padd.vertex = true; + padd.p = start; + sc.pts.Add(&padd); + Point2d pa, pb; + Vector np, npc; + // Better to start with a too-small step, so that we don't miss + // features of the curve entirely. + bool fwd; + double tol, step = SS.ChordTolMm(); + for(a = 0; a < 100; a++) { + ClosestPointTo(start, &pa); + b->ClosestPointTo(start, &pb); + + Vector na = NormalAt(pa).WithMagnitude(1), + nb = b->NormalAt(pb).WithMagnitude(1); + + if(a == 0) { + Vector dp = nb.Cross(na); + if(dp.Dot(startv) < 0) { + // We want to march in the more inward direction. + fwd = true; + } else { + fwd = false; + } + } + + int i = 0; + do { + if(i++ > 20) break; + + Vector dp = nb.Cross(na); + if(!fwd) dp = dp.ScaledBy(-1); + dp = dp.WithMagnitude(step); + + np = start.Plus(dp); + npc = ClosestPointOnThisAndSurface(b, np); + tol = (npc.Minus(np)).Magnitude(); + + if(tol > SS.ChordTolMm()*0.8) { + step *= 0.95; + } else { + step /= 0.95; + } + } while(tol > SS.ChordTolMm() || tol < SS.ChordTolMm()/2); + + SPoint *sp; + for(sp = spl.l.First(); sp; sp = spl.l.NextAfter(sp)) { + if((sp->p).OnLineSegment(start, npc, 2*SS.ChordTolMm())) { + sp->tag = 1; + a = 1000; + npc = sp->p; + } + } + + padd.p = npc; + padd.vertex = (a == 1000); + sc.pts.Add(&padd); + + start = npc; + } + + spl.l.RemoveTagged(); + + // And now we split and insert the curve SCurve split = sc.MakeCopySplitAgainst(agnstA, agnstB, this, b); sc.Clear(); - into->curve.AddAndAssignId(&split); + break; } spl.Clear(); } @@ -848,13 +922,46 @@ bool SShell::ClassifyEdge(int *indir, int *outdir, } if(edge_inters != 0) dbp("bad, edge_inters=%d", edge_inters); - + + // Next, check for edge-on-surface. The ray-casting for edge-inside-shell + // would catch this too, but test separately, for speed (since many edges + // are on surface) and for numerical stability, so we don't pick up + // the additional error from the line intersection. + + for(srf = surface.First(); srf; srf = surface.NextAfter(srf)) { + if(srf->LineEntirelyOutsideBbox(ea, eb, true)) continue; + + Point2d puv; + srf->ClosestPointTo(p, &(puv.x), &(puv.y), false); + Vector pp = srf->PointAt(puv); + + if((pp.Minus(p)).Magnitude() > LENGTH_EPS) continue; + Point2d dummy = { 0, 0 }, ia = { 0, 0 }, ib = { 0, 0 }; + int c = srf->bsp->ClassifyPoint(puv, dummy, &ia, &ib); + if(c == SBspUv::OUTSIDE) continue; + + // Edge-on-face (unless edge-on-edge above superceded) + Point2d pin, pout; + srf->ClosestPointTo(p.Plus(edge_n_in), &pin, false); + srf->ClosestPointTo(p.Plus(edge_n_out), &pout, false); + + Vector surf_n_in = srf->NormalAt(pin), + surf_n_out = srf->NormalAt(pout); + + *indir = ClassifyRegion(edge_n_in, surf_n_in, surf_n); + *outdir = ClassifyRegion(edge_n_out, surf_n_out, surf_n); + return true; + } + + // Edge is not on face or on edge; so it's either inside or outside + // the shell, and we'll determine which by raycasting. 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 // testing) Vector ray = Vector::From(Random(1), Random(1), Random(1)); + AllPointsIntersecting( p.Minus(ray), p.Plus(ray), &l, false, true, false); @@ -883,28 +990,14 @@ bool SShell::ClassifyEdge(int *indir, int *outdir, if(d < dmin) { dmin = d; - - if(d < LENGTH_EPS) { - // Edge-on-face (unless edge-on-edge above supercedes) - 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); + // Edge does not lie on surface; either strictly inside + // or strictly outside + if((si->surfNormal).Dot(ray) > 0) { + *indir = INSIDE; + *outdir = INSIDE; } else { - // Edge does not lie on surface; either strictly inside - // or strictly outside - if((si->surfNormal).Dot(ray) > 0) { - *indir = INSIDE; - *outdir = INSIDE; - } else { - *indir = OUTSIDE; - *outdir = OUTSIDE; - } + *indir = OUTSIDE; + *outdir = OUTSIDE; } onEdge = si->onEdge; } diff --git a/util.cpp b/util.cpp index e2a8006..0df2d76 100644 --- a/util.cpp +++ b/util.cpp @@ -459,15 +459,15 @@ double Vector::DistanceToLine(Vector p0, Vector dp) { return ((this->Minus(p0)).Cross(dp)).Magnitude() / m; } -bool Vector::OnLineSegment(Vector a, Vector b) { - if(this->Equals(a) || this->Equals(b)) return true; +bool Vector::OnLineSegment(Vector a, Vector b, double tol) { + if(this->Equals(a, tol) || this->Equals(b, tol)) return true; Vector d = b.Minus(a); double m = d.MagSquared(); double distsq = ((this->Minus(a)).Cross(d)).MagSquared() / m; - if(distsq >= LENGTH_EPS*LENGTH_EPS) return false; + if(distsq >= tol*tol) return false; double t = (this->Minus(a)).DivPivoting(d); // On-endpoint already tested