From 666ea1c047ab2faa9a72759ff5e4ec1572cca924 Mon Sep 17 00:00:00 2001 From: Jonathan Westhues Date: Thu, 18 Jun 2009 23:56:33 -0800 Subject: [PATCH] Add beginnings of marching surface intersection; I can find all the boundary points, at least. That required some changes to what gets passed around (for example because to project a point onto this inexact curve, we need to know which two surfaces it trims so that we can do a Newton's method on them). And fix stupidity in the way that I calculated edge normals; I just did normal in uv space, and there's no particular reason why that would be normal in xyz. So edges in long skinny surfaces failed, for example. [git-p4: depot-paths = "//depot/solvespace/": change = 1990] --- polygon.cpp | 21 ++++++++++++ polygon.h | 9 +++++ srf/boolean.cpp | 86 ++++++++++++++++++++++++----------------------- srf/curve.cpp | 20 +++++++++++ srf/ratpoly.cpp | 4 ++- srf/surface.h | 12 ++++--- srf/surfinter.cpp | 73 ++++++++++++++++++++++++++++++++++++++-- wishlist.txt | 6 +++- 8 files changed, 180 insertions(+), 51 deletions(-) diff --git a/polygon.cpp b/polygon.cpp index cc0ee539..dd6523df 100644 --- a/polygon.cpp +++ b/polygon.cpp @@ -285,6 +285,27 @@ void SEdgeList::MergeCollinearSegments(Vector a, Vector b) { l.RemoveTagged(); } +void SPointList::Clear(void) { + l.Clear(); +} + +bool SPointList::ContainsPoint(Vector pt) { + SPoint *p; + for(p = l.First(); p; p = l.NextAfter(p)) { + if(pt.Equals(p->p)) { + return true; + } + } + return false; +} + +void SPointList::Add(Vector pt) { + SPoint p; + ZERO(&p); + p.p = pt; + l.Add(&p); +} + void SContour::AddPoint(Vector p) { SPoint sp; sp.tag = 0; diff --git a/polygon.h b/polygon.h index e61ba189..34c50d6a 100644 --- a/polygon.h +++ b/polygon.h @@ -42,6 +42,15 @@ public: Vector p; }; +class SPointList { +public: + List l; + + void Clear(void); + bool ContainsPoint(Vector pt); + void Add(Vector pt); +}; + class SContour { public: int tag; diff --git a/srf/boolean.cpp b/srf/boolean.cpp index b8a368fc..56797de8 100644 --- a/srf/boolean.cpp +++ b/srf/boolean.cpp @@ -266,14 +266,13 @@ void SSurface::EdgeNormalsWithinSurface(Point2d auv, Point2d buv, Vector *pt, Vector *enin, Vector *enout, Vector *surfn, - DWORD auxA, SShell *shell) + DWORD auxA, + SShell *shell, SShell *sha, SShell *shb) { // 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); @@ -288,15 +287,35 @@ void SSurface::EdgeNormalsWithinSurface(Point2d auv, Point2d buv, sc->exact.ClosestPointTo(*pt, &t, false); *pt = sc->exact.PointAt(t); ClosestPointTo(*pt, &muv); + } else if(!sc->isExact) { + SSurface *trimmedA = sc->GetSurfaceA(sha, shb), + *trimmedB = sc->GetSurfaceB(sha, shb); + *pt = trimmedA->ClosestPointOnThisAndSurface(trimmedB, *pt); + ClosestPointTo(*pt, &muv); } *surfn = NormalAt(muv.x, muv.y); + // Compute the edge's inner normal in xyz space. + Vector ab = (PointAt(auv)).Minus(PointAt(buv)), + enxyz = (ab.Cross(*surfn)).WithMagnitude(SS.ChordTolMm()); + // And based on that, compute the edge's inner normal in uv space. This + // vector is perpendicular to the edge in xyz, but not necessarily in uv. + Vector tu, tv; + TangentsAt(muv.x, muv.y, &tu, &tv); + Point2d enuv; + enuv.x = enxyz.Dot(tu) / tu.MagSquared(); + enuv.y = enxyz.Dot(tv) / tv.MagSquared(); + // Don't let the magnitude get too tiny at small chord tolerances; we + // will otherwise have numerical problems subtracting nearly-equal + // numbers. + enuv = enuv.WithMagnitude(0.01); + // 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)); + Vector pin = PointAt(muv.Minus(enuv)), + pout = PointAt(muv.Plus(enuv)); *enin = pin.Minus(*pt), *enout = pout.Minus(*pt); } @@ -307,10 +326,14 @@ void SSurface::EdgeNormalsWithinSurface(Point2d auv, Point2d buv, // also need a pointer to the shell that contains our own surface, since that // contains our original trim curves. //----------------------------------------------------------------------------- -SSurface SSurface::MakeCopyTrimAgainst(SShell *agnst, SShell *parent, +SSurface SSurface::MakeCopyTrimAgainst(SShell *parent, + SShell *sha, SShell *shb, SShell *into, - int type, bool opA) + int type) { + bool opA = (parent == sha); + SShell *agnst = opA ? shb : sha; + SSurface ret; // The returned surface is identical, just the trim curves change ret = *this; @@ -400,7 +423,7 @@ SSurface SSurface::MakeCopyTrimAgainst(SShell *agnst, SShell *parent, Vector pt, enin, enout, surfn; ret.EdgeNormalsWithinSurface(auv, buv, &pt, &enin, &enout, &surfn, - se->auxA, into); + se->auxA, into, sha, shb); int indir_shell, outdir_shell, indir_orig, outdir_orig; @@ -424,7 +447,7 @@ SSurface SSurface::MakeCopyTrimAgainst(SShell *agnst, SShell *parent, Vector pt, enin, enout, surfn; ret.EdgeNormalsWithinSurface(auv, buv, &pt, &enin, &enout, &surfn, - se->auxA, into); + se->auxA, into, sha, shb); int indir_shell, outdir_shell, indir_orig, outdir_orig; @@ -467,13 +490,13 @@ SSurface SSurface::MakeCopyTrimAgainst(SShell *agnst, SShell *parent, return ret; } -void SShell::CopySurfacesTrimAgainst(SShell *against, SShell *into, - int type, bool opA) +void SShell::CopySurfacesTrimAgainst(SShell *sha, SShell *shb, SShell *into, + int type) { SSurface *ss; for(ss = surface.First(); ss; ss = surface.NextAfter(ss)) { SSurface ssn; - ssn = ss->MakeCopyTrimAgainst(against, this, into, type, opA); + ssn = ss->MakeCopyTrimAgainst(this, sha, shb, into, type); ss->newH = into->surface.AddAndAssignId(&ssn); I++; } @@ -506,20 +529,9 @@ void SShell::CleanupAfterBoolean(void) { //----------------------------------------------------------------------------- void SShell::RewriteSurfaceHandlesForCurves(SShell *a, SShell *b) { SCurve *sc; - for(sc = curve.First(); sc; sc = curve.NextAfter(sc)) { - if(sc->source == SCurve::FROM_A) { - sc->surfA = a->surface.FindById(sc->surfA)->newH; - sc->surfB = a->surface.FindById(sc->surfB)->newH; - } else if(sc->source == SCurve::FROM_B) { - sc->surfA = b->surface.FindById(sc->surfA)->newH; - sc->surfB = b->surface.FindById(sc->surfB)->newH; - } else if(sc->source == SCurve::FROM_INTERSECTION) { - sc->surfA = a->surface.FindById(sc->surfA)->newH; - sc->surfB = b->surface.FindById(sc->surfB)->newH; - } else { - oops(); - } + sc->surfA = sc->GetSurfaceA(a, b)->newH, + sc->surfB = sc->GetSurfaceB(a, b)->newH; } } @@ -591,17 +603,9 @@ void SShell::MakeFromBoolean(SShell *a, SShell *b, int type) { SCurve *sc; for(sc = curve.First(); sc; sc = curve.NextAfter(sc)) { - SSurface *srfA, *srfB; - if(sc->source == SCurve::FROM_A) { - srfA = a->surface.FindById(sc->surfA); - srfB = a->surface.FindById(sc->surfB); - } else if(sc->source == SCurve::FROM_B) { - srfA = b->surface.FindById(sc->surfA); - srfB = b->surface.FindById(sc->surfB); - } else if(sc->source == SCurve::FROM_INTERSECTION) { - srfA = a->surface.FindById(sc->surfA); - srfB = b->surface.FindById(sc->surfB); - } + SSurface *srfA = sc->GetSurfaceA(a, b), + *srfB = sc->GetSurfaceB(a, b); + sc->RemoveShortSegments(srfA, srfB); } @@ -613,16 +617,14 @@ void SShell::MakeFromBoolean(SShell *a, SShell *b, int type) { a->MakeClassifyingBsps(this); b->MakeClassifyingBsps(this); - // Then trim and copy the surfaces if(b->surface.n == 0 || a->surface.n == 0) { - I = 1023123; - a->CopySurfacesTrimAgainst(b, this, type, true); - b->CopySurfacesTrimAgainst(a, this, type, false); + I = 1000000; } else { I = 0; - a->CopySurfacesTrimAgainst(b, this, type, true); - b->CopySurfacesTrimAgainst(a, this, type, false); } + // Then trim and copy the surfaces + a->CopySurfacesTrimAgainst(a, b, this, type); + b->CopySurfacesTrimAgainst(a, b, this, type); // Now that we've copied the surfaces, we know their new hSurfaces, so // rewrite the curves to refer to the surfaces by their handles in the diff --git a/srf/curve.cpp b/srf/curve.cpp index b1dd282a..b5d50475 100644 --- a/srf/curve.cpp +++ b/srf/curve.cpp @@ -401,6 +401,26 @@ void SCurve::Clear(void) { pts.Clear(); } +SSurface *SCurve::GetSurfaceA(SShell *a, SShell *b) { + if(source == FROM_A) { + return a->surface.FindById(surfA); + } else if(source == FROM_B) { + return b->surface.FindById(surfA); + } else if(source == FROM_INTERSECTION) { + return a->surface.FindById(surfA); + } else oops(); +} + +SSurface *SCurve::GetSurfaceB(SShell *a, SShell *b) { + if(source == FROM_A) { + return a->surface.FindById(surfB); + } else if(source == FROM_B) { + return b->surface.FindById(surfB); + } else if(source == FROM_INTERSECTION) { + return b->surface.FindById(surfB); + } else oops(); +} + //----------------------------------------------------------------------------- // When we split line segments wherever they intersect a surface, we introduce // extra pwl points. This may create very short edges that could be removed diff --git a/srf/ratpoly.cpp b/srf/ratpoly.cpp index 17d57462..b46b5e8d 100644 --- a/srf/ratpoly.cpp +++ b/srf/ratpoly.cpp @@ -415,7 +415,7 @@ Vector SSurface::ClosestPointOnThisAndSurface(SSurface *srf2, Vector p) { (srf[j])->ClosestPointTo(p, &(puv[j]), false); } - for(i = 0; i < 4; i++) { + for(i = 0; i < 10; i++) { Vector tu[2], tv[2], cp[2], n[2]; double d[2]; @@ -428,6 +428,8 @@ Vector SSurface::ClosestPointOnThisAndSurface(SSurface *srf2, Vector p) { d[j] = (n[j]).Dot(cp[j]); } + if((cp[0]).Equals(cp[1], RATPOLY_EPS)) break; + Vector p0 = Vector::AtIntersectionOfPlanes(n[0], d[0], n[1], d[1]), dp = (n[0]).Cross(n[1]); diff --git a/srf/surface.h b/srf/surface.h index 5d56c5f9..5fd08805 100644 --- a/srf/surface.h +++ b/srf/surface.h @@ -159,6 +159,8 @@ public: SCurve MakeCopySplitAgainst(SShell *agnstA, SShell *agnstB, SSurface *srfA, SSurface *srfB); void RemoveShortSegments(SSurface *srfA, SSurface *srfB); + SSurface *GetSurfaceA(SShell *a, SShell *b); + SSurface *GetSurfaceB(SShell *a, SShell *b); void Clear(void); }; @@ -225,9 +227,10 @@ public: 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); + DWORD auxA, + SShell *shell, SShell *sha, SShell *shb); + SSurface MakeCopyTrimAgainst(SShell *parent, SShell *a, SShell *b, + SShell *into, int type); void TrimFromEdgeList(SEdgeList *el, bool asUv); void IntersectAgainst(SSurface *b, SShell *agnstA, SShell *agnstB, SShell *into); @@ -305,7 +308,8 @@ public: static const int AS_INTERSECT = 12; void MakeFromBoolean(SShell *a, SShell *b, int type); void CopyCurvesSplitAgainst(bool opA, SShell *agnst, SShell *into); - void CopySurfacesTrimAgainst(SShell *against, SShell *into, int t, bool a); + void CopySurfacesTrimAgainst(SShell *sha, SShell *shb, SShell *into, + int type); void MakeIntersectionCurvesAgainst(SShell *against, SShell *into); void MakeClassifyingBsps(SShell *useCurvesFrom); void AllPointsIntersecting(Vector a, Vector b, List *il, diff --git a/srf/surfinter.cpp b/srf/surfinter.cpp index 49a4309d..68359bc3 100644 --- a/srf/surfinter.cpp +++ b/srf/surfinter.cpp @@ -276,10 +276,77 @@ void SSurface::IntersectAgainst(SSurface *b, SShell *agnstA, SShell *agnstB, inters.Clear(); lv.Clear(); - } + } else { + // Try intersecting the surfaces numerically, by a marching algorithm. + // First, we find all the intersections between a surface and the + // boundary of the other surface. + SPointList spl; + ZERO(&spl); + int a; + for(a = 0; a < 2; a++) { + SShell *shA = (a == 0) ? agnstA : agnstB, + *shB = (a == 0) ? agnstB : agnstA; + SSurface *srfA = (a == 0) ? this : b, + *srfB = (a == 0) ? b : this; - // need to implement general numerical surface intersection for tough - // cases, just giving up for now + SEdgeList el; + ZERO(&el); + srfA->MakeEdgesInto(shA, &el, false, NULL); + + SEdge *se; + for(se = el.l.First(); se; se = el.l.NextAfter(se)) { + List lsi; + ZERO(&lsi); + + srfB->AllPointsIntersecting(se->a, se->b, &lsi, + true, true, false); + if(lsi.n == 0) continue; + + // Find the other surface that this curve trims. + hSCurve hsc = { se->auxA }; + SCurve *sc = shA->curve.FindById(hsc); + hSSurface hother = (sc->surfA.v == srfA->h.v) ? + sc->surfB : sc->surfA; + SSurface *other = shA->surface.FindById(hother); + + SInter *si; + 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); + } + lsi.Clear(); + } + + el.Clear(); + } + + SPoint *sp; + if(spl.l.n == 2) { + + SCurve sc; + ZERO(&sc); + sc.surfA = h; + sc.surfB = b->h; + 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); + + SCurve split = sc.MakeCopySplitAgainst(agnstA, agnstB, this, b); + sc.Clear(); + + into->curve.AddAndAssignId(&split); + } + spl.Clear(); + } } diff --git a/wishlist.txt b/wishlist.txt index 5e01eef3..655a89f7 100644 --- a/wishlist.txt +++ b/wishlist.txt @@ -1,8 +1,12 @@ +marching algorithm for surface intersection +grid +classification of edges only as necessary (all same for span between intersections) +associative entities from solid model, as a special group +rounding, as a special group ----- line styles (color, thickness) -marching algorithm for surface intersection loop detection IGES export incremental regen of entities