From d3dcd8fb2369948acca8d7a96bf4ba9c9bfb9c36 Mon Sep 17 00:00:00 2001 From: Jonathan Westhues Date: Sun, 21 Jun 2009 01:02:36 -0800 Subject: [PATCH] Now we are actually marching. There seems to be either a numerical problem or a tendency to generate backwards edges or both, need to debug that. But it generates the curve, and begins to work. And change the edge classification. Now instead of testing for point-on-surface using the results of the raycasting, test for point-on-surface as a separate step. That stops us from picking up the additional numerical error from the surface-line intersection, which may be significant if the ray is parallel or almost parallel to the surface. [git-p4: depot-paths = "//depot/solvespace/": change = 1991] --- drawentity.cpp | 6 +- dsc.h | 2 +- polygon.h | 1 + srf/ratpoly.cpp | 10 ++- srf/surfinter.cpp | 167 ++++++++++++++++++++++++++++++++++++---------- util.cpp | 6 +- 6 files changed, 145 insertions(+), 47 deletions(-) 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