From bc70089dd08a447dc25b0133069af8a68738d433 Mon Sep 17 00:00:00 2001 From: Jonathan Westhues Date: Sun, 8 Mar 2009 02:59:57 -0800 Subject: [PATCH] Add code to subdivide (with de Castljau's algorithm) a surface, and use that for surface-line intersections. That has major problems with the heuristic on when to stop and do Newton polishing. There's also an issue with all the Newton stuff when surfaces join tangent. And update the wishlist to reflect current needs. [git-p4: depot-paths = "//depot/solvespace/": change = 1925] --- dsc.h | 2 + solvespace.cpp | 3 + solvespace.h | 3 +- srf/boolean.cpp | 34 +++--- srf/ratpoly.cpp | 8 +- srf/surface.h | 14 +++ srf/surfinter.cpp | 247 ++++++++++++++++++++++++++++++++++---------- srf/triangulate.cpp | 2 +- util.cpp | 34 ++++++ wishlist.txt | 18 +++- 10 files changed, 282 insertions(+), 83 deletions(-) diff --git a/dsc.h b/dsc.h index 3a921ca..dc64f41 100644 --- a/dsc.h +++ b/dsc.h @@ -84,6 +84,8 @@ public: void MakeMaxMin(Vector *maxv, Vector *minv); static bool BoundingBoxesDisjoint(Vector amax, Vector amin, Vector bmax, Vector bmin); + static bool BoundingBoxIntersectsLine(Vector amax, Vector amin, + Vector p0, Vector p1, bool segment); bool OutsideAndNotOn(Vector maxv, Vector minv); Point2d Project2d(Vector u, Vector v); Point2d ProjectXy(void); diff --git a/solvespace.cpp b/solvespace.cpp index 23b9a98..a1f4f72 100644 --- a/solvespace.cpp +++ b/solvespace.cpp @@ -177,6 +177,9 @@ double SolveSpace::StringToMm(char *str) { return atof(str); } } +double SolveSpace::ChordTolMm(void) { + return SS.chordTol / SS.GW.scale; +} void SolveSpace::AfterNewFile(void) { diff --git a/solvespace.h b/solvespace.h index 70b1954..d377c55 100644 --- a/solvespace.h +++ b/solvespace.h @@ -31,7 +31,7 @@ inline double WRAP_NOT_0(double v, double n) { #define ZERO(v) memset((v), 0, sizeof(*(v))) #define CO(v) (v).x, (v).y, (v).z -#define LENGTH_EPS (1e-7) +#define LENGTH_EPS (1e-6) #define VERY_POSITIVE (1e10) #define VERY_NEGATIVE (-1e10) @@ -448,6 +448,7 @@ public: char *MmToString(double v); double ExprToMm(Expr *e); double StringToMm(char *s); + double ChordTolMm(void); // The platform-dependent code calls this before entering the msg loop void Init(char *cmdLine); diff --git a/srf/boolean.cpp b/srf/boolean.cpp index 714c11d..f298ec2 100644 --- a/srf/boolean.cpp +++ b/srf/boolean.cpp @@ -50,24 +50,30 @@ SCurve SCurve::MakeCopySplitAgainst(SShell *agnstA, SShell *agnstB, if(agnstA) agnstA->AllPointsIntersecting(prev, *p, &il, true, true); if(agnstB) agnstB->AllPointsIntersecting(prev, *p, &il, true, true); - // If any intersections exist, sort them in order along the - // line and add them to the curve. if(il.n > 0) { - LineStart = prev; - LineDirection = p->Minus(prev); - qsort(il.elem, il.n, sizeof(il.elem[0]), ByTAlongLine); - - Vector prev = Vector::From(VERY_POSITIVE, 0, 0); + // The intersections were generated by intersecting the pwl + // edge against a surface; so they must be refined to lie + // exactly on the original curve. SInter *pi; for(pi = il.First(); pi; pi = il.NextAfter(pi)) { - // The intersection was generated by intersecting the pwl - // edge against a surface; so it must be refined to lie - // exactly on the original curve. double u, v; (pi->srf)->ClosestPointTo(pi->p, &u, &v, false); (pi->srf)->PointOnSurfaces(srfA, srfB, &u, &v); pi->p = (pi->srf)->PointAt(u, v); + } + // And now sort them in order along the line. Note that we must + // do that after refining, in case the refining would make two + // points switch places. + LineStart = prev; + LineDirection = p->Minus(prev); + qsort(il.elem, il.n, sizeof(il.elem[0]), ByTAlongLine); + + // And now uses the intersections to generate our split pwl edge(s) + Vector prev = Vector::From(VERY_POSITIVE, 0, 0); + for(pi = il.First(); pi; pi = il.NextAfter(pi)) { + double t = (pi->p.Minus(LineStart)).DivPivoting(LineDirection); + dbp("t=%.3f", t); // On-edge intersection will generate same split point for // both surfaces, so don't create zero-length edge. if(!prev.Equals(pi->p)) { @@ -236,6 +242,7 @@ void DBPEDGE(int tag) { } void DEBUGEDGELIST(SEdgeList *sel, SSurface *surf) { + dbp("print %d edges", sel->l.n); SEdge *se; for(se = sel->l.First(); se; se = sel->l.NextAfter(se)) { Vector mid = (se->a).Plus(se->b).ScaledBy(0.5); @@ -372,11 +379,7 @@ SSurface SSurface::MakeCopyTrimAgainst(SShell *agnst, SShell *parent, eb = ret.PointAt(buv.x, buv.y); Vector ptout = n.Cross((eb.Minus(ea))); - if(I == 1 && fabs(((se->b).Minus(se->a).Magnitude()) - 0.166) < 0.01) { - FLAG = 1; - } int c_shell = agnst->ClassifyPoint(pt, ptout); - FLAG = 0; int tag = 0; tag |= INSIDE(ORIG, INDIR); @@ -456,7 +459,7 @@ SSurface SSurface::MakeCopyTrimAgainst(SShell *agnst, SShell *parent, } final.l.RemoveTagged(); -// if(I == 1) DEBUGEDGELIST(&final, &ret); +// if(I == 0) DEBUGEDGELIST(&final, &ret); // Use our reassembled edges to trim the new surface. ret.TrimFromEdgeList(&final); @@ -491,6 +494,7 @@ void SShell::MakeIntersectionCurvesAgainst(SShell *agnst, SShell *into) { // list for into. sa->IntersectAgainst(sb, agnst, this, into); } + FLAG++; } } diff --git a/srf/ratpoly.cpp b/srf/ratpoly.cpp index 1b02b07..4faa43c 100644 --- a/srf/ratpoly.cpp +++ b/srf/ratpoly.cpp @@ -4,7 +4,7 @@ // independently projected into uv and back, to end up equal with the // LENGTH_EPS. Best case that requires LENGTH_EPS/2, but more is better // and convergence should be fast by now. -#define RATPOLY_EPS (LENGTH_EPS/(1e1)) +#define RATPOLY_EPS (LENGTH_EPS/(1e2)) static double Bernstein(int k, int deg, double t) { @@ -131,10 +131,8 @@ void SBezier::MakePwlWorker(List *l, double ta, double tb, Vector off) { double d = max(pm1.DistanceToLine(pa, pb.Minus(pa)), pm2.DistanceToLine(pa, pb.Minus(pa))); - double tol = SS.chordTol/SS.GW.scale; - double step = 1.0/SS.maxSegments; - if((tb - ta) < step || d < tol) { + if((tb - ta) < step || d < SS.ChordTolMm()) { // A previous call has already added the beginning of our interval. pb = pb.Plus(off); l->Add(&pb); @@ -536,7 +534,7 @@ void SSurface::ClosestPointTo(Vector p, double *u, double *v, bool converge) { *u = *v = 0; // a plane, perfect no matter what the initial guess } else { double minDist = VERY_POSITIVE; - double res = 7.0; + double res = (max(degm, degn) == 2) ? 7.0 : 20.0; for(i = 0; i < (int)res; i++) { for(j = 0; j <= (int)res; j++) { double tryu = (i/res), tryv = (j/res); diff --git a/srf/surface.h b/srf/surface.h index 245231e..aa1e042 100644 --- a/srf/surface.h +++ b/srf/surface.h @@ -194,8 +194,22 @@ public: SShell *into); void AddExactIntersectionCurve(SBezier *sb, SSurface *srfB, SShell *agnstA, SShell *agnstB, SShell *into); + typedef struct { + int tag; + Point2d p; + } Inter; + void WeightControlPoints(void); + void UnWeightControlPoints(void); + void CopyRowOrCol(bool row, int this_ij, SSurface *src, int src_ij); + void BlendRowOrCol(bool row, int this_ij, SSurface *a, int a_ij, + SSurface *b, int b_ij); + void SplitInHalf(bool byU, SSurface *sa, SSurface *sb); void AllPointsIntersecting(Vector a, Vector b, List *l, bool seg, bool trimmed); + void AllPointsIntersectingUntrimmed(Vector a, Vector b, + int *cnt, int *level, + List *l, bool segment, + SSurface *sorig); void ClosestPointTo(Vector p, double *u, double *v, bool converge=true); bool PointIntersectingLine(Vector p0, Vector p1, double *u, double *v); diff --git a/srf/surfinter.cpp b/srf/surfinter.cpp index 8bd0c48..3352987 100644 --- a/srf/surfinter.cpp +++ b/srf/surfinter.cpp @@ -11,16 +11,22 @@ void SSurface::AddExactIntersectionCurve(SBezier *sb, SSurface *srfB, sc.surfB = srfB->h; sb->MakePwlInto(&(sc.pts)); -/* - Vector *prev = NULL, *v; - for(v = sc.pts.First(); v; v = sc.pts.NextAfter(v)) { - if(prev) SS.nakedEdges.AddEdge(*prev, *v); - prev = v; - } */ - // Now split the line where it intersects our existing surfaces SCurve split = sc.MakeCopySplitAgainst(agnstA, agnstB, this, srfB); sc.Clear(); + +/* + if(sb->deg == 1) { + dbp(" "); + Vector *prev = NULL, *v; + dbp("split.pts.n =%d", split.pts.n); + for(v = split.pts.First(); v; v = split.pts.NextAfter(v)) { + if(prev) { + SS.nakedEdges.AddEdge(*prev, *v); + } + prev = v; + } + } */ split.source = SCurve::FROM_INTERSECTION; into->curve.AddAndAssignId(&split); @@ -88,7 +94,10 @@ void SSurface::IntersectAgainst(SSurface *b, SShell *agnstA, SShell *agnstB, inters.Clear(); } else { // Direction of extrusion is not parallel to plane; so - // intersection is projection of extruded curve into our plane + // intersection is projection of extruded curve into our plane. + // If both curves are planes, then we could do it either way; + // so choose the one that generates the shorter curve. + // XXX TODO int i; for(i = 0; i <= bezier.deg; i++) { @@ -107,6 +116,177 @@ void SSurface::IntersectAgainst(SSurface *b, SShell *agnstA, SShell *agnstB, // cases, just giving up for now } +void SSurface::WeightControlPoints(void) { + int i, j; + for(i = 0; i <= degm; i++) { + for(j = 0; j <= degn; j++) { + ctrl[i][j] = (ctrl[i][j]).ScaledBy(weight[i][j]); + } + } +} +void SSurface::UnWeightControlPoints(void) { + int i, j; + for(i = 0; i <= degm; i++) { + for(j = 0; j <= degn; j++) { + ctrl[i][j] = (ctrl[i][j]).ScaledBy(1.0/weight[i][j]); + } + } +} +void SSurface::CopyRowOrCol(bool row, int this_ij, SSurface *src, int src_ij) { + if(row) { + int j; + for(j = 0; j <= degn; j++) { + ctrl [this_ij][j] = src->ctrl [src_ij][j]; + weight[this_ij][j] = src->weight[src_ij][j]; + } + } else { + int i; + for(i = 0; i <= degm; i++) { + ctrl [i][this_ij] = src->ctrl [i][src_ij]; + weight[i][this_ij] = src->weight[i][src_ij]; + } + } +} +void SSurface::BlendRowOrCol(bool row, int this_ij, SSurface *a, int a_ij, + SSurface *b, int b_ij) +{ + if(row) { + int j; + for(j = 0; j <= degn; j++) { + Vector c = (a->ctrl [a_ij][j]).Plus(b->ctrl [b_ij][j]); + double w = (a->weight[a_ij][j] + b->weight[b_ij][j]); + ctrl [this_ij][j] = c.ScaledBy(0.5); + weight[this_ij][j] = w / 2; + } + } else { + int i; + for(i = 0; i <= degm; i++) { + Vector c = (a->ctrl [i][a_ij]).Plus(b->ctrl [i][b_ij]); + double w = (a->weight[i][a_ij] + b->weight[i][b_ij]); + ctrl [i][this_ij] = c.ScaledBy(0.5); + weight[i][this_ij] = w / 2; + } + } +} +void SSurface::SplitInHalf(bool byU, SSurface *sa, SSurface *sb) { + sa->degm = sb->degm = degm; + sa->degn = sb->degn = degn; + + // by de Casteljau's algorithm in a projective space; so we must work + // on points (w*x, w*y, w*z, w) + WeightControlPoints(); + + switch(byU ? degm : degn) { + case 1: + sa->CopyRowOrCol (byU, 0, this, 0); + sb->CopyRowOrCol (byU, 1, this, 1); + + sa->BlendRowOrCol(byU, 1, this, 0, this, 1); + sb->BlendRowOrCol(byU, 0, this, 0, this, 1); + break; + + case 2: + sa->CopyRowOrCol (byU, 0, this, 0); + sb->CopyRowOrCol (byU, 2, this, 2); + + sa->BlendRowOrCol(byU, 1, this, 0, this, 1); + sb->BlendRowOrCol(byU, 1, this, 1, this, 2); + + sa->BlendRowOrCol(byU, 2, sa, 1, sb, 1); + sb->BlendRowOrCol(byU, 0, sa, 1, sb, 1); + break; + + case 3: { + SSurface st; + st.degm = degm; st.degn = degn; + + sa->CopyRowOrCol (byU, 0, this, 0); + sb->CopyRowOrCol (byU, 3, this, 3); + + sa->BlendRowOrCol(byU, 1, this, 0, this, 1); + sb->BlendRowOrCol(byU, 2, this, 2, this, 3); + st. BlendRowOrCol(byU, 0, this, 1, this, 2); // scratch var + + sa->BlendRowOrCol(byU, 2, sa, 1, &st, 0); + sb->BlendRowOrCol(byU, 1, sb, 2, &st, 0); + + sa->BlendRowOrCol(byU, 3, sa, 2, sb, 1); + sb->BlendRowOrCol(byU, 0, sa, 2, sb, 1); + break; + } + + default: oops(); + } + + sa->UnWeightControlPoints(); + sb->UnWeightControlPoints(); + UnWeightControlPoints(); +} + +//----------------------------------------------------------------------------- +// Find all points where the indicated finite (if segment) or infinite (if not +// segment) line intersects our surface. Report them in uv space in the list. +// We first do a bounding box check; if the line doesn't intersect, then we're +// done. If it does, then we check how small our surface is. If it's big, +// then we subdivide into quarters and recurse. If it's small, then we refine +// by Newton's method and record the point. +//----------------------------------------------------------------------------- +void SSurface::AllPointsIntersectingUntrimmed(Vector a, Vector b, + int *cnt, int *level, + List *l, bool segment, + SSurface *sorig) +{ + // 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(*cnt > 2000) { + dbp("!!! too many subdivisions (level=%d)!", *level); + return; + } + (*cnt)++; + + // If we might intersect, and the surface is small, then switch to Newton + // iterations. + double h = max(amax.x - amin.x, + max(amax.y - amin.y, + amax.z - amin.z)); + if(fabs(h) < SS.ChordTolMm()) { + Vector p = (amax.Plus(amin)).ScaledBy(0.5); + Inter inter; + sorig->ClosestPointTo(p, &(inter.p.x), &(inter.p.y), false); + if(sorig->PointIntersectingLine(a, b, &(inter.p.x), &(inter.p.y))) { + Vector p = sorig->PointAt(inter.p.x, inter.p.y); + // Debug check, verify that the point lies in both surfaces + // (which it ought to, since the surfaces should be coincident) + double u, v; + ClosestPointTo(p, &u, &v); + l->Add(&inter); + } else { + // Might not converge if line is almost tangent to surface... + } + return; + } + + // But the surface is big, so split it, alternating by u and v + SSurface surf0, surf1; + SplitInHalf((*level & 1) == 0, &surf0, &surf1); + + int nextLevel = (*level) + 1; + (*level) = nextLevel; + surf0.AllPointsIntersectingUntrimmed(a, b, cnt, level, l, segment, sorig); + (*level) = nextLevel; + surf1.AllPointsIntersectingUntrimmed(a, b, cnt, level, l, segment, sorig); +} + //----------------------------------------------------------------------------- // 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 @@ -119,55 +299,17 @@ void SSurface::AllPointsIntersecting(Vector a, Vector b, Vector ba = b.Minus(a); double bam = ba.Magnitude(); - typedef struct { - int tag; - Point2d p; - } Inter; List inters; ZERO(&inters); // First, get all the intersections between the infinite ray and the // untrimmed surface. - int i, j; - for(i = 0; i < degm; i++) { - for(j = 0; j < degn; j++) { - // Intersect the ray with each face in the control polyhedron - Vector p00 = ctrl[i][j], - p01 = ctrl[i][j+1], - p10 = ctrl[i+1][j]; - - Vector u = p01.Minus(p00), v = p10.Minus(p00); - - Vector n = (u.Cross(v)).WithMagnitude(1); - double d = n.Dot(p00); - - bool parallel; - Vector pi = - Vector::AtIntersectionOfPlaneAndLine(n, d, a, b, ¶llel); - if(parallel) continue; - - double ui = (pi.Minus(p00)).Dot(u) / u.MagSquared(), - vi = (pi.Minus(p00)).Dot(v) / v.MagSquared(); - - double tol = 1e-2; - if(ui < -tol || ui > 1 + tol || vi < -tol || vi > 1 + tol) { - continue; - } - - Inter inter; - ClosestPointTo(pi, &inter.p.x, &inter.p.y, false); - if(PointIntersectingLine(a, b, &inter.p.x, &inter.p.y)) { - inters.Add(&inter); - } else { - // No convergence; but this might not be an error, since - // the line can intersect the convex hull more times than - // it intersects the surface itself. - } - } - } + int cnt = 0, level = 0; + AllPointsIntersectingUntrimmed(a, b, &cnt, &level, &inters, seg, this); // Remove duplicate intersection points inters.ClearTags(); + int i, j; for(i = 0; i < inters.n; i++) { for(j = i + 1; j < inters.n; j++) { if(inters.elem[i].p.Equals(inters.elem[j].p)) { @@ -246,11 +388,6 @@ int SShell::ClassifyPoint(Vector p, Vector pout) { bool onEdge = false; edge_inters = 0; - if(FLAG) { - dbp("inters=%d cnt=%d", l.n, cnt); - SS.nakedEdges.AddEdge(p, p.Plus(ray.WithMagnitude(50))); - } - SInter *si; for(si = l.First(); si; si = l.NextAfter(si)) { double t = ((si->p).Minus(p)).DivPivoting(ray); @@ -259,8 +396,6 @@ int SShell::ClassifyPoint(Vector p, Vector pout) { continue; } - if(FLAG) SS.nakedEdges.AddEdge(p, si->p); - double d = ((si->p).Minus(p)).Magnitude(); // Handle edge-on-edge diff --git a/srf/triangulate.cpp b/srf/triangulate.cpp index 823c0ee..b942e00 100644 --- a/srf/triangulate.cpp +++ b/srf/triangulate.cpp @@ -299,7 +299,7 @@ void SContour::UvTriangulateInto(SMesh *m, SSurface *srf) { bestEar = ear; bestChordTol = tol; } - if(bestChordTol < 0.1*(SS.chordTol / SS.GW.scale)) { + if(bestChordTol < 0.1*SS.ChordTolMm()) { break; } } diff --git a/util.cpp b/util.cpp index 8a02d70..c5d3328 100644 --- a/util.cpp +++ b/util.cpp @@ -573,6 +573,40 @@ bool Vector::BoundingBoxesDisjoint(Vector amax, Vector amin, return false; } +bool Vector::BoundingBoxIntersectsLine(Vector amax, Vector amin, + Vector p0, Vector p1, bool segment) +{ + Vector dp = p1.Minus(p0); + double lp = dp.Magnitude(); + dp = dp.ScaledBy(1.0/lp); + + int i, a; + for(i = 0; i < 3; i++) { + int j = WRAP(i+1, 3), k = WRAP(i+2, 3); + if(lp*fabs(dp.Element(i)) < LENGTH_EPS) continue; // parallel to plane + + for(a = 0; a < 2; a++) { + double d = (a == 0) ? amax.Element(i) : amin.Element(i); + // n dot (p0 + t*dp) = d + // (n dot p0) + t * (n dot dp) = d + double t = (d - p0.Element(i)) / dp.Element(i); + Vector p = p0.Plus(dp.ScaledBy(t)); + + if(segment && (t < -LENGTH_EPS || t > (lp+LENGTH_EPS))) continue; + + if(p.Element(j) > amax.Element(j) + LENGTH_EPS) continue; + if(p.Element(k) > amax.Element(k) + LENGTH_EPS) continue; + + if(p.Element(j) < amin.Element(j) - LENGTH_EPS) continue; + if(p.Element(k) < amin.Element(k) - LENGTH_EPS) continue; + + return true; + } + } + + return false; +} + Vector Vector::AtIntersectionOfPlanes(Vector n1, double d1, Vector n2, double d2) { diff --git a/wishlist.txt b/wishlist.txt index b01b280..e68b78b 100644 --- a/wishlist.txt +++ b/wishlist.txt @@ -1,11 +1,19 @@ -leak fixing +marching algorithm for surface intersection +surfaces of revolution (lathed) +fix surface-line intersection +cylinder-line and plane-line special cases +exact boundaries when near pwl trim +step and repeat rotate/translate +short pwl edge avoidance +faces +take consistent pwl with coincident faces +exact curve export (at least for dxf) +hidden line removal from mesh +line styles (color, thickness) - -long term ----- +loop detection incremental regen of entities? -my own plane poly triangulation code -exact curved surfaces