diff --git a/src/srf/boolean.cpp b/src/srf/boolean.cpp index cb72cd8f..17e9d15c 100644 --- a/src/srf/boolean.cpp +++ b/src/srf/boolean.cpp @@ -88,8 +88,15 @@ SCurve SCurve::MakeCopySplitAgainst(SShell *agnstA, SShell *agnstB, } } - // We're keeping the intersection, so actually refine it. - (pi->srf)->PointOnSurfaces(srfA, srfB, &(puv.x), &(puv.y)); + // We're keeping the intersection, so actually refine it. Finding the intersection + // to within EPS is important to match the ends of different chopped trim curves. + // The general 3-surface intersection fails to refine for trims where surfaces + // are tangent at the curve, but those trims are usually exact, so… + if(isExact) { + (pi->srf)->PointOnCurve(&exact, &(puv.x), &(puv.y)); + } else { + (pi->srf)->PointOnSurfaces(srfA, srfB, &(puv.x), &(puv.y)); + } pi->p = (pi->srf)->PointAt(puv); } il.RemoveTagged(); diff --git a/src/srf/ratpoly.cpp b/src/srf/ratpoly.cpp index 6498cc4a..bd878ab1 100644 --- a/src/srf/ratpoly.cpp +++ b/src/srf/ratpoly.cpp @@ -634,3 +634,58 @@ void SSurface::PointOnSurfaces(SSurface *s1, SSurface *s2, double *up, double *v dbp("didn't converge (three surfaces intersecting)"); } +void SSurface::PointOnCurve(const SBezier *curve, double *up, double *vp) +{ + Vector tu,tv,n; + double u = *up, v = *vp; + Vector ps = PointAt(u, v); + // Get initial guesses for t on the curve + double tCurve = 0.5; + curve->ClosestPointTo(ps, &tCurve, /*mustConverge=*/false); + if(tCurve < 0.0) tCurve = 0.0; + if(tCurve > 1.0) tCurve = 1.0; + + for(int i = 0; i < 30; i++) { + // Approximate the surface by a plane + Vector ps = PointAt(u, v); + TangentsAt(u, v, &tu, &tv); + n = tu.Cross(tv).WithMagnitude(1); + + // point on curve and tangent line direction + Vector pc = curve->PointAt(tCurve); + Vector tc = curve->TangentAt(tCurve); + + if(ps.Equals(pc, RATPOLY_EPS)) { + *up = u; + *vp = v; + return; + } + + //pi is where the curve tangent line intersects the surface tangent plane + Vector pi; + double d = tc.Dot(n); + if (fabs(d) < 1e-10) { // parallel line and plane, guess the average rather than fail + pi = pc.Plus(ps).ScaledBy(0.5); + } else { + pi = pc.Minus(tc.ScaledBy(pc.Minus(ps).Dot(n)/d)); + } + + // project the point onto the tangent plane and line + { + Vector n = tu.Cross(tv); + Vector ty = n.Cross(tu).ScaledBy(1.0/tu.MagSquared()); + Vector tx = tv.Cross(n).ScaledBy(1.0/tv.MagSquared()); + + Vector dp = pi.Minus(ps); + double du = dp.Dot(tx), dv = dp.Dot(ty); + + u += du / tx.MagSquared(); + v += dv / ty.MagSquared(); + } + tCurve += pi.Minus(pc).Dot(tc) / tc.MagSquared(); + if(tCurve < 0.0) tCurve = 0.0; + if(tCurve > 1.0) tCurve = 1.0; + } + dbp("didn't converge (surface and curve intersecting)"); +} + diff --git a/src/srf/surface.h b/src/srf/surface.h index d1918b6c..ba4414a2 100644 --- a/src/srf/surface.h +++ b/src/srf/surface.h @@ -332,6 +332,7 @@ public: bool PointIntersectingLine(Vector p0, Vector p1, double *u, double *v) const; Vector ClosestPointOnThisAndSurface(SSurface *srf2, Vector p); void PointOnSurfaces(SSurface *s1, SSurface *s2, double *u, double *v); + void PointOnCurve(const SBezier *curve, double *up, double *vp); Vector PointAt(double u, double v) const; Vector PointAt(Point2d puv) const; void TangentsAt(double u, double v, Vector *tu, Vector *tv) const; diff --git a/src/srf/surfinter.cpp b/src/srf/surfinter.cpp index 9f6c85cc..0a2f09c0 100644 --- a/src/srf/surfinter.cpp +++ b/src/srf/surfinter.cpp @@ -341,7 +341,11 @@ void SSurface::IntersectAgainst(SSurface *b, SShell *agnstA, SShell *agnstB, Vector p = si->p; double u, v; srfB->ClosestPointTo(p, &u, &v); - srfB->PointOnSurfaces(srfA, other, &u, &v); + if(sc->isExact) { + srfB->PointOnCurve(&(sc->exact), &u, &v); + } else { + srfB->PointOnSurfaces(srfA, other, &u, &v); + } p = srfB->PointAt(u, v); if(!spl.ContainsPoint(p)) { SPoint sp;