diff --git a/dsc.h b/dsc.h index 3a921ca0..dc64f41c 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 23b9a98c..a1f4f721 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 70b1954b..d377c55e 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 714c11d0..f298ec2a 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 1b02b074..4faa43c7 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 245231ef..aa1e042d 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 8bd0c480..33529877 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 823c0ee7..b942e009 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 8a02d70e..c5d3328d 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 b01b2804..e68b78b3 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