From 7f3dd91bd93e1d0e28658d8e1a44d2e3d5d235af Mon Sep 17 00:00:00 2001 From: Jonathan Westhues Date: Thu, 19 Mar 2009 09:40:11 -0800 Subject: [PATCH] Add a special case for line-cylinder intersection, solving in closed form. This is a fairly good speedup, and handles tangency well. But that shows that tangency has other problems; need to classify edges correctly (whether they point to a coincident surface) in curved surfaces too. I need to tweak SShell::ClassifyPoint(). [git-p4: depot-paths = "//depot/solvespace/": change = 1933] --- dsc.h | 1 + solvespace.h | 6 +++ srf/boolean.cpp | 12 ++--- srf/ratpoly.cpp | 31 +++++++++++-- srf/surface.h | 23 +++++----- srf/surfinter.cpp | 110 +++++++++++++++++++++++++++++++++++----------- util.cpp | 8 ++++ 7 files changed, 146 insertions(+), 45 deletions(-) diff --git a/dsc.h b/dsc.h index 84d9222..05e531f 100644 --- a/dsc.h +++ b/dsc.h @@ -103,6 +103,7 @@ public: Point2d Plus(Point2d b); Point2d Minus(Point2d b); Point2d ScaledBy(double s); + double DivPivoting(Point2d delta); double Dot(Point2d p); double DistanceTo(Point2d p); double DistanceToLine(Point2d p0, Point2d dp, bool segment); diff --git a/solvespace.h b/solvespace.h index a182f6d..30d4cad 100644 --- a/solvespace.h +++ b/solvespace.h @@ -26,6 +26,12 @@ inline double WRAP_NOT_0(double v, double n) { while(v <= 0) v += n; return v; } +inline double WRAP_SYMMETRIC(double v, double n) { + // Clamp it to the range (-n/2, n/2] + while(v > n/2) v -= n; + while(v <= -n/2) v += n; + return v; +} #define SWAP(T, a, b) do { T temp = (a); (a) = (b); (b) = temp; } while(0) #define ZERO(v) memset((v), 0, sizeof(*(v))) diff --git a/srf/boolean.cpp b/srf/boolean.cpp index edde3d9..3bee9fd 100644 --- a/srf/boolean.cpp +++ b/srf/boolean.cpp @@ -47,8 +47,8 @@ SCurve SCurve::MakeCopySplitAgainst(SShell *agnstA, SShell *agnstB, ZERO(&il); // Find all the intersections with the two passed shells - if(agnstA) agnstA->AllPointsIntersecting(prev, *p, &il, true, true); - if(agnstB) agnstB->AllPointsIntersecting(prev, *p, &il, true, true); + if(agnstA) agnstA->AllPointsIntersecting(prev, *p, &il, true,true,true); + if(agnstB) agnstB->AllPointsIntersecting(prev, *p, &il, true,true,true); if(il.n > 0) { // The intersections were generated by intersecting the pwl @@ -377,12 +377,12 @@ SSurface SSurface::MakeCopyTrimAgainst(SShell *agnst, SShell *parent, 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 n = ret.NormalAt(am.x, am.y); + 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 ptout = n.Cross((eb.Minus(ea))); + Vector edge_n = surf_n.Cross((eb.Minus(ea))); - int c_shell = agnst->ClassifyPoint(pt, ptout); + int c_shell = agnst->ClassifyPoint(pt, edge_n, surf_n); int tag = 0; tag |= INSIDE(ORIG, INDIR); @@ -448,7 +448,7 @@ SSurface SSurface::MakeCopyTrimAgainst(SShell *agnst, SShell *parent, ZERO(&poly); final.l.ClearTags(); if(!final.AssemblePolygon(&poly, NULL, true)) { - DEBUGEDGELIST(&inter, &ret); + DEBUGEDGELIST(&final, &ret); } poly.Clear(); diff --git a/srf/ratpoly.cpp b/srf/ratpoly.cpp index 8a80cb0..12aa57c 100644 --- a/srf/ratpoly.cpp +++ b/srf/ratpoly.cpp @@ -479,7 +479,10 @@ bool SSurface::IsCylinder(Vector *center, Vector *axis, double *r, double rd0 = center->Minus(sb.ctrl[0]).Magnitude(), rd2 = center->Minus(sb.ctrl[2]).Magnitude(); - if(fabs(rd0 - rd2) > LENGTH_EPS) return false; + if(fabs(rd0 - rd2) > LENGTH_EPS) { + return false; + } + *r = rd0; Vector u = r0.WithMagnitude(1), v = (axis->Cross(u)).WithMagnitude(1); @@ -490,11 +493,18 @@ bool SSurface::IsCylinder(Vector *center, Vector *axis, double *r, 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(dtheta > PI) { + // Not possible with a second order Bezier arc; so we must have + // the points backwards. + dtheta = 2*PI - dtheta; + } - if(fabs(sb.weight[1] - cos(dtheta/2)) > LENGTH_EPS) return false; + if(fabs(sb.weight[1] - cos(dtheta/2)) > LENGTH_EPS) { + return false; + } - *start = (sb.ctrl[0]).Minus(*center); - *finish = (sb.ctrl[2]).Minus(*center); + *start = sb.ctrl[0]; + *finish = sb.ctrl[2]; return true; } @@ -754,6 +764,19 @@ void SSurface::GetAxisAlignedBounding(Vector *ptMax, Vector *ptMin) { } } +bool SSurface::LineEntirelyOutsideBbox(Vector a, Vector b, bool segment) { + Vector amax, amin; + GetAxisAlignedBounding(&amax, &amin); + if(!Vector::BoundingBoxIntersectsLine(amax, amin, a, b, segment)) { + // The line segment could fail to intersect the bbox, but lie entirely + // within it and intersect the surface. + if(a.OutsideAndNotOn(amax, amin) && b.OutsideAndNotOn(amax, amin)) { + return true; + } + } + return false; +} + void SSurface::MakeEdgesInto(SShell *shell, SEdgeList *sel, bool asUv, SShell *useCurvesFrom) { STrimBy *stb; diff --git a/srf/surface.h b/srf/surface.h index a022e14..894ce9c 100644 --- a/srf/surface.h +++ b/srf/surface.h @@ -198,6 +198,7 @@ public: SShell *into); void AddExactIntersectionCurve(SBezier *sb, SSurface *srfB, SShell *agnstA, SShell *agnstB, SShell *into); + typedef struct { int tag; Point2d p; @@ -210,7 +211,8 @@ public: double DepartureFromCoplanar(void); void SplitInHalf(bool byU, SSurface *sa, SSurface *sb); void AllPointsIntersecting(Vector a, Vector b, - List *l, bool seg, bool trimmed); + List *l, + bool seg, bool trimmed, bool inclTangent); void AllPointsIntersectingUntrimmed(Vector a, Vector b, int *cnt, int *level, List *l, bool segment, @@ -222,6 +224,7 @@ public: Vector PointAt(double u, double v); void TangentsAt(double u, double v, Vector *tu, Vector *tv); Vector NormalAt(double u, double v); + bool LineEntirelyOutsideBbox(Vector a, Vector b, bool segment); void GetAxisAlignedBounding(Vector *ptMax, Vector *ptMin); bool CoincidentWithPlane(Vector n, double d); bool CoincidentWith(SSurface *ss, bool sameNormal); @@ -258,20 +261,20 @@ public: void MakeIntersectionCurvesAgainst(SShell *against, SShell *into); void MakeClassifyingBsps(void); void AllPointsIntersecting(Vector a, Vector b, List *il, - bool seg, bool trimmed); + bool seg, bool trimmed, bool inclTangent); void MakeCoincidentEdgesInto(SSurface *proto, bool sameNormal, SEdgeList *el, SShell *useCurvesFrom); 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; + 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 out); + int ClassifyPoint(Vector p, Vector edge_n, Vector surf_n); void MakeFromCopyOf(SShell *a); diff --git a/srf/surfinter.cpp b/srf/surfinter.cpp index 658b60e..751f7ba 100644 --- a/srf/surfinter.cpp +++ b/srf/surfinter.cpp @@ -175,7 +175,8 @@ void SSurface::IntersectAgainst(SSurface *b, SShell *agnstA, SShell *agnstB, List inters; ZERO(&inters); - sext->AllPointsIntersecting(p0, p0.Plus(dp), &inters, false, false); + sext->AllPointsIntersecting( + p0, p0.Plus(dp), &inters, false, false, true); SInter *si; for(si = inters.First(); si; si = inters.NextAfter(si)) { @@ -387,16 +388,8 @@ void SSurface::AllPointsIntersectingUntrimmed(Vector a, Vector b, { // Test if the line intersects our axis-aligned bounding box; if no, then // no possibility of an intersection - Vector amax, amin; - GetAxisAlignedBounding(&amax, &amin); - if(!Vector::BoundingBoxIntersectsLine(amax, amin, a, b, segment)) { - // The line segment could fail to intersect the bbox, but lie entirely - // within it and intersect the surface. - if(a.OutsideAndNotOn(amax, amin) && b.OutsideAndNotOn(amax, amin)) { - return; - } - } - + if(LineEntirelyOutsideBbox(a, b, segment)) return; + if(*cnt > 2000) { dbp("!!! too many subdivisions (level=%d)!", *level); dbp("degm = %d degn = %d", degm, degn); @@ -407,7 +400,10 @@ void SSurface::AllPointsIntersectingUntrimmed(Vector a, Vector b, // If we might intersect, and the surface is small, then switch to Newton // iterations. if(DepartureFromCoplanar() < 0.2*SS.ChordTolMm()) { - Vector p = (amax.Plus(amin)).ScaledBy(0.5); + Vector p = (ctrl[0 ][0 ]).Plus( + ctrl[0 ][degn]).Plus( + ctrl[degm][0 ]).Plus( + ctrl[degm][degn]).ScaledBy(0.25); Inter inter; sorig->ClosestPointTo(p, &(inter.p.x), &(inter.p.y), false); if(sorig->PointIntersectingLine(a, b, &(inter.p.x), &(inter.p.y))) { @@ -438,11 +434,16 @@ void SSurface::AllPointsIntersectingUntrimmed(Vector a, Vector b, // Find all points where a line through a and b intersects our surface, and // add them to the list. If seg is true then report only intersections that // lie within the finite line segment (not including the endpoints); otherwise -// we work along the infinite line. +// we work along the infinite line. And we report either just intersections +// inside the trim curve, or any intersection with u, v in [0, 1]. And we +// either disregard or report tangent points. //----------------------------------------------------------------------------- void SSurface::AllPointsIntersecting(Vector a, Vector b, - List *l, bool seg, bool trimmed) + List *l, + bool seg, bool trimmed, bool inclTangent) { + if(LineEntirelyOutsideBbox(a, b, seg)) return; + Vector ba = b.Minus(a); double bam = ba.Magnitude(); @@ -468,8 +469,62 @@ void SSurface::AllPointsIntersecting(Vector a, Vector b, ClosestPointTo(p, &(inter.p.x), &(inter.p.y)); inters.Add(&inter); } - } else if(IsCylinder(¢er, &axis, &radius, &start, &finish) && 0) { - // XXX, cylinder is easy in closed form + } else if(IsCylinder(¢er, &axis, &radius, &start, &finish)) { + // This one can be solved in closed form too. + Vector ab = b.Minus(a); + if(axis.Cross(ab).Magnitude() < LENGTH_EPS) { + // edge is parallel to axis of cylinder, no intersection points + return; + } + // A coordinate system centered at the center of the circle, with + // the edge under test horizontal + Vector u, v, n = axis.WithMagnitude(1); + u = (ab.Minus(n.ScaledBy(ab.Dot(n)))).WithMagnitude(1); + v = n.Cross(u); + Point2d ap = (a.Minus(center)).DotInToCsys(u, v, n).ProjectXy(), + bp = (b.Minus(center)).DotInToCsys(u, v, n).ProjectXy(), + sp = (start. Minus(center)).DotInToCsys(u, v, n).ProjectXy(), + fp = (finish.Minus(center)).DotInToCsys(u, v, n).ProjectXy(); + + double thetas = atan2(sp.y, sp.x), thetaf = atan2(fp.y, fp.x); + + Point2d ip[2]; + int ip_n = 0; + if(fabs(fabs(ap.y) - radius) < LENGTH_EPS) { + // tangent + if(inclTangent) { + ip[0] = Point2d::From(0, ap.y); + ip_n = 1; + } + } else if(fabs(ap.y) < radius) { + // two intersections + double xint = sqrt(radius*radius - ap.y*ap.y); + ip[0] = Point2d::From(-xint, ap.y); + ip[1] = Point2d::From( xint, ap.y); + ip_n = 2; + } + int i; + for(i = 0; i < ip_n; i++) { + double t = (ip[i].Minus(ap)).DivPivoting(bp.Minus(ap)); + // This is a point on the circle; but is it on the arc? + Point2d pp = ap.Plus((bp.Minus(ap)).ScaledBy(t)); + double theta = atan2(pp.y, pp.x); + double dp = WRAP_SYMMETRIC(theta - thetas, 2*PI), + df = WRAP_SYMMETRIC(thetaf - thetas, 2*PI); + double tol = LENGTH_EPS/radius; + + if((df > 0 && ((dp < -tol) || (dp > df + tol))) || + (df < 0 && ((dp > tol) || (dp < df - tol)))) + { + continue; + } + + Vector p = a.Plus((b.Minus(a)).ScaledBy(t)); + + Inter inter; + ClosestPointTo(p, &(inter.p.x), &(inter.p.y)); + inters.Add(&inter); + } } else { // General numerical solution by subdivision, fallback int cnt = 0, level = 0; @@ -519,24 +574,25 @@ void SSurface::AllPointsIntersecting(Vector a, Vector b, } void SShell::AllPointsIntersecting(Vector a, Vector b, - List *il, bool seg, bool trimmed) + List *il, + bool seg, bool trimmed, bool inclTangent) { SSurface *ss; for(ss = surface.First(); ss; ss = surface.NextAfter(ss)) { - ss->AllPointsIntersecting(a, b, il, seg, trimmed); + ss->AllPointsIntersecting(a, b, il, seg, trimmed, inclTangent); } } //----------------------------------------------------------------------------- -// Does the given point lie on our shell? It might be inside or outside, or -// it might be on the surface with pout parallel or anti-parallel to the -// intersecting surface's normal. +// 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 +// possibilities. // // To calculate, we intersect a ray through p with our shell, and classify // 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 pout) { +int SShell::ClassifyPoint(Vector p, Vector edge_n, Vector surf_n) { List l; ZERO(&l); @@ -550,7 +606,8 @@ int SShell::ClassifyPoint(Vector p, Vector pout) { // 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); + AllPointsIntersecting( + p.Minus(ray), p.Plus(ray), &l, false, true, false); double dmin = VERY_POSITIVE; ret = OUTSIDE; // no intersections means it's outside @@ -569,7 +626,7 @@ int SShell::ClassifyPoint(Vector p, Vector pout) { // Handle edge-on-edge if(d < LENGTH_EPS && si->onEdge && edge_inters < 2) { - edge_dotp[edge_inters] = (si->surfNormal).Dot(pout); + edge_dotp[edge_inters] = (si->surfNormal).Dot(edge_n); edge_inters++; } @@ -577,7 +634,10 @@ int SShell::ClassifyPoint(Vector p, Vector pout) { dmin = d; if(d < LENGTH_EPS) { // Edge-on-face (unless edge-on-edge above supercedes) - if((si->surfNormal).Dot(pout) > 0) { + 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; diff --git a/util.cpp b/util.cpp index 82f9c0e..5efafb7 100644 --- a/util.cpp +++ b/util.cpp @@ -758,6 +758,14 @@ Point2d Point2d::ScaledBy(double s) { return r; } +double Point2d::DivPivoting(Point2d delta) { + if(fabs(delta.x) > fabs(delta.y)) { + return x/delta.x; + } else { + return y/delta.y; + } +} + double Point2d::MagSquared(void) { return x*x + y*y; }