Discard intersection curves that lie entirely outside of one

surface's domain of u, v in [0, 1]. Cache the starting guess when
projecting a point into a ratpoly surface, to avoid brute force
searching for a good one every time. Split edges even if they
aren't quite inside the trim curve, since the trim boundaries are
pwl, not exact; unnecessary splits won't hurt, but failure to split
when necessary will. Make the triangulation code use a better (but
not perfect) epsilon, to avoid "can't find ear" failures on very
fine meshes.

And turn on compiler optimization! I had somehow forgotten about
that, and it's a ~2x improvement.

[git-p4: depot-paths = "//depot/solvespace/": change = 2026]
solver
Jonathan Westhues 2009-08-20 20:58:28 -08:00
parent 36182448ce
commit f7f9000c68
11 changed files with 171 additions and 67 deletions

View File

@ -1,7 +1,7 @@
DEFINES = /D_WIN32_WINNT=0x500 /DISOLATION_AWARE_ENABLED /D_WIN32_IE=0x500 /DWIN32_LEAN_AND_MEAN /DWIN32
# Use the multi-threaded static libc because libpng and zlib do; not sure if anything bad
# happens if those mix, but don't want to risk it.
CFLAGS = /W3 /nologo -MT -Iextlib -I..\common\win32 /D_DEBUG /D_CRT_SECURE_NO_WARNINGS /I. /Zi /EHs
CFLAGS = /W3 /nologo -MT -Iextlib -I..\common\win32 /D_DEBUG /D_CRT_SECURE_NO_WARNINGS /I. /Zi /EHs /O2
HEADERS = ..\common\win32\freeze.h ui.h solvespace.h dsc.h sketch.h expr.h polygon.h srf\surface.h
@ -61,7 +61,7 @@ LIBS = user32.lib gdi32.lib comctl32.lib advapi32.lib shell32.lib opengl32.lib g
all: $(OBJDIR)/solvespace.exe
@cp $(OBJDIR)/solvespace.exe .
solvespace t8.slvs
solvespace fail.slvs
clean:
rm -f obj/*

View File

@ -685,8 +685,9 @@ void EntityBase::GenerateEquations(IdList<Equation,hEquation> *l) {
// If the two endpoints of the arc are constrained coincident
// (to make a complete circle), then our distance constraint
// would be redundant and therefore overconstrain things.
Constraint *c;
for(c = SK.constraint.First(); c; c = SK.constraint.NextAfter(c)) {
int i;
for(i = 0; i < SK.constraint.n; i++) {
ConstraintBase *c = &(SK.constraint.elem[i]);
if(c->group.v != group.v) continue;
if(c->type != Constraint::POINTS_COINCIDENT) continue;
@ -696,7 +697,7 @@ void EntityBase::GenerateEquations(IdList<Equation,hEquation> *l) {
break;
}
}
if(c) break;
if(i < SK.constraint.n) break;
Expr *ra = Constraint::Distance(workplane, point[0], point[1]);
Expr *rb = Constraint::Distance(workplane, point[0], point[2]);

View File

@ -1,7 +1,7 @@
DEFINES = /D_WIN32_WINNT=0x500 /DISOLATION_AWARE_ENABLED /D_WIN32_IE=0x500 /DWIN32_LEAN_AND_MEAN /DWIN32 /DLIBRARY
# Use the multi-threaded static libc because libpng and zlib do; not sure if anything bad
# happens if those mix, but don't want to risk it.
CFLAGS = /W3 /nologo -MT -Iextlib -I..\..\common\win32 /D_DEBUG /D_CRT_SECURE_NO_WARNINGS /I. /I.. /Zi /EHs
CFLAGS = /W3 /nologo -MT -Iextlib -I..\..\common\win32 /D_DEBUG /D_CRT_SECURE_NO_WARNINGS /I. /I.. /Zi /EHs /O2
HEADERS = ..\solvespace.h ..\dsc.h ..\sketch.h ..\expr.h slvs.h

View File

@ -76,9 +76,9 @@ public:
void FindPointWithMinX(void);
Vector AnyEdgeMidpoint(void);
bool IsEar(int bp);
bool IsEar(int bp, double scaledEps);
bool BridgeToContour(SContour *sc, SEdgeList *el, List<Vector> *vl);
void ClipEarInto(SMesh *m, int bp);
void ClipEarInto(SMesh *m, int bp, double scaledEps);
void UvTriangulateInto(SMesh *m, SSurface *srf);
};

View File

@ -48,9 +48,9 @@ SCurve SCurve::MakeCopySplitAgainst(SShell *agnstA, SShell *agnstB,
// Find all the intersections with the two passed shells
if(agnstA)
agnstA->AllPointsIntersecting(prev.p, p->p, &il, true, true, true);
agnstA->AllPointsIntersecting(prev.p, p->p, &il, true, false, true);
if(agnstB)
agnstB->AllPointsIntersecting(prev.p, p->p, &il, true, true, true);
agnstB->AllPointsIntersecting(prev.p, p->p, &il, true, false, true);
if(il.n > 0) {
// The intersections were generated by intersecting the pwl
@ -69,10 +69,28 @@ SCurve SCurve::MakeCopySplitAgainst(SShell *agnstA, SShell *agnstB,
pi->tag = 1;
continue;
}
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);
Point2d puv;
(pi->srf)->ClosestPointTo(pi->p, &puv, false);
// Split the edge if the intersection lies within the surface's
// trim curves, or within the chord tol of the trim curve; want
// some slop if points are close to edge and pwl is too coarse,
// and it doesn't hurt to split unnecessarily.
Point2d dummy = { 0, 0 };
int c = pi->srf->bsp->ClassifyPoint(puv, dummy, pi->srf);
if(c == SBspUv::OUTSIDE) {
double d;
d = pi->srf->bsp->MinimumDistanceToEdge(puv, pi->srf);
if(d > SS.ChordTolMm()) {
pi->tag = 1;
continue;
}
}
// We're keeping the intersection, so actually refine it.
(pi->srf)->PointOnSurfaces(srfA, srfB, &(puv.x), &(puv.y));
pi->p = (pi->srf)->PointAt(puv);
}
il.RemoveTagged();
@ -423,7 +441,7 @@ SSurface SSurface::MakeCopyTrimAgainst(SShell *parent,
// the split curves (not the original curves).
SEdgeList orig;
ZERO(&orig);
ret.MakeEdgesInto(into, &orig, true);
ret.MakeEdgesInto(into, &orig, AS_UV);
ret.trim.Clear();
// which means that we can't necessarily use the old BSP...
SBspUv *origBsp = SBspUv::From(&orig, &ret);
@ -763,12 +781,12 @@ void SSurface::MakeClassifyingBsp(SShell *shell, SShell *useCurvesFrom) {
SEdgeList el;
ZERO(&el);
MakeEdgesInto(shell, &el, true, useCurvesFrom);
MakeEdgesInto(shell, &el, AS_UV, useCurvesFrom);
bsp = SBspUv::From(&el, this);
el.Clear();
ZERO(&edges);
MakeEdgesInto(shell, &edges, false, useCurvesFrom);
MakeEdgesInto(shell, &edges, AS_XYZ, useCurvesFrom);
}
SBspUv *SBspUv::Alloc(void) {
@ -942,3 +960,16 @@ int SBspUv::ClassifyEdge(Point2d ea, Point2d eb, SSurface *srf) {
return ret;
}
double SBspUv::MinimumDistanceToEdge(Point2d p, SSurface *srf) {
if(!this) return VERY_POSITIVE;
double dn = neg->MinimumDistanceToEdge(p, srf),
dp = pos->MinimumDistanceToEdge(p, srf);
Point2d as = a, bs = b;
ScalePoints(&p, &as, &bs, srf);
double d = p.DistanceToLine(as, bs.Minus(as), true);
return min(d, min(dn, dp));
}

View File

@ -22,7 +22,7 @@ void SShell::MergeCoincidentSurfaces(void) {
SEdgeList sel;
ZERO(&sel);
si->MakeEdgesInto(this, &sel, false);
si->MakeEdgesInto(this, &sel, SSurface::AS_XYZ);
bool mergedThisTime, merged = false;
do {
@ -42,7 +42,7 @@ void SShell::MergeCoincidentSurfaces(void) {
// less robust.
SEdgeList tel;
ZERO(&tel);
sj->MakeEdgesInto(this, &tel, false);
sj->MakeEdgesInto(this, &tel, SSurface::AS_XYZ);
if(!sel.ContainsEdgeFrom(&tel)) {
tel.Clear();
continue;
@ -52,7 +52,7 @@ void SShell::MergeCoincidentSurfaces(void) {
sj->tag = 1;
merged = true;
mergedThisTime = true;
sj->MakeEdgesInto(this, &sel, false);
sj->MakeEdgesInto(this, &sel, SSurface::AS_XYZ);
sj->trim.Clear();
// All the references to this surface get replaced with the

View File

@ -356,6 +356,19 @@ void SSurface::ClosestPointTo(Vector p, double *u, double *v, bool converge) {
if(p.Equals(ctrl[degm][degn])) { *u = 1; *v = 1; return; }
if(p.Equals(ctrl[0] [degn])) { *u = 0; *v = 1; return; }
// Try whatever the previous guess was. This is likely to do something
// good if we're working our way along a curve or something else where
// we project successive points that are close to each other; something
// like a 20% speedup empirically.
if(converge) {
double ut = cached.x, vt = cached.y;
if(ClosestPointNewton(p, &ut, &vt, converge)) {
cached.x = *u = ut;
cached.y = *v = vt;
return;
}
}
// Search for a reasonable initial guess
int i, j;
if(degm == 1 && degn == 1) {
@ -378,13 +391,27 @@ void SSurface::ClosestPointTo(Vector p, double *u, double *v, bool converge) {
}
}
if(ClosestPointNewton(p, u, v, converge)) {
cached.x = *u;
cached.y = *v;
return;
}
// If we failed to converge, then at least don't return NaN.
if(isnan(*u) || isnan(*v)) {
*u = *v = 0;
}
}
bool SSurface::ClosestPointNewton(Vector p, double *u, double *v, bool converge)
{
// Initial guess is in u, v; refine by Newton iteration.
Vector p0;
for(i = 0; i < (converge ? 25 : 5); i++) {
for(int i = 0; i < (converge ? 25 : 5); i++) {
p0 = PointAt(*u, *v);
if(converge) {
if(p0.Equals(p, RATPOLY_EPS)) {
return;
return true;
}
}
@ -406,10 +433,7 @@ void SSurface::ClosestPointTo(Vector p, double *u, double *v, bool converge) {
dbp("want %.3f %.3f %.3f", CO(p));
dbp("distance = %g", (p.Minus(p0)).Magnitude());
}
if(isnan(*u) || isnan(*v)) {
*u = *v = 0;
}
return false;
}
bool SSurface::PointIntersectingLine(Vector p0, Vector p1, double *u, double *v)

View File

@ -194,10 +194,10 @@ bool SSurface::LineEntirelyOutsideBbox(Vector a, Vector b, bool segment) {
// Generate the piecewise linear approximation of the trim stb, which applies
// to the curve sc.
//-----------------------------------------------------------------------------
void SSurface::MakeTrimEdgesInto(SEdgeList *sel, bool asUv,
void SSurface::MakeTrimEdgesInto(SEdgeList *sel, int flags,
SCurve *sc, STrimBy *stb)
{
Vector prev, prevuv, ptuv;
Vector prev;
bool inCurve = false, empty = true;
double u = 0, v = 0;
@ -212,23 +212,22 @@ void SSurface::MakeTrimEdgesInto(SEdgeList *sel, bool asUv,
increment = 1;
}
for(i = first; i != (last + increment); i += increment) {
Vector *pt = &(sc->pts.elem[i].p);
if(asUv) {
Vector tpt, *pt = &(sc->pts.elem[i].p);
if(flags & AS_UV) {
ClosestPointTo(*pt, &u, &v);
ptuv = Vector::From(u, v, 0);
if(inCurve) {
sel->AddEdge(prevuv, ptuv, sc->h.v, stb->backwards);
empty = false;
}
prevuv = ptuv;
tpt = Vector::From(u, v, 0);
} else {
if(inCurve) {
sel->AddEdge(prev, *pt, sc->h.v, stb->backwards);
empty = false;
}
prev = *pt;
tpt = *pt;
}
if(inCurve) {
sel->AddEdge(prev, tpt, sc->h.v, stb->backwards);
empty = false;
}
prev = tpt; // either uv or xyz, depending on flags
if(pt->Equals(stb->start)) inCurve = true;
if(pt->Equals(stb->finish)) inCurve = false;
}
@ -242,7 +241,7 @@ void SSurface::MakeTrimEdgesInto(SEdgeList *sel, bool asUv,
// the split curves from useCurvesFrom instead of the curves in our own
// shell.
//-----------------------------------------------------------------------------
void SSurface::MakeEdgesInto(SShell *shell, SEdgeList *sel, bool asUv,
void SSurface::MakeEdgesInto(SShell *shell, SEdgeList *sel, int flags,
SShell *useCurvesFrom)
{
STrimBy *stb;
@ -257,7 +256,7 @@ void SSurface::MakeEdgesInto(SShell *shell, SEdgeList *sel, bool asUv,
sc = useCurvesFrom->curve.FindById(sc->newH);
}
MakeTrimEdgesInto(sel, asUv, sc, stb);
MakeTrimEdgesInto(sel, flags, sc, stb);
}
}
@ -384,7 +383,7 @@ void SSurface::MakeSectionEdgesInto(SShell *shell,
sp = fpt;
}
} else {
if(sel) MakeTrimEdgesInto(sel, false, sc, stb);
if(sel) MakeTrimEdgesInto(sel, AS_XYZ, sc, stb);
}
}
}
@ -393,21 +392,21 @@ void SSurface::TriangulateInto(SShell *shell, SMesh *sm) {
SEdgeList el;
ZERO(&el);
MakeEdgesInto(shell, &el, true);
MakeEdgesInto(shell, &el, AS_UV);
SPolygon poly;
ZERO(&poly);
if(el.AssemblePolygon(&poly, NULL, true)) {
int i, start = sm->l.n;
if(degm == 1 && degn == 1) {
// A plane; triangulate any old way
poly.UvTriangulateInto(sm, NULL);
} else if(degm == 1 || degn == 1) {
// A surface with curvature along one direction only; so
// choose the triangulation with chords that lie as much
// as possible within the surface. And since the trim curves
// have been pwl'd to within the desired chord tol, that will
// produce a surface good to within roughly that tol.
//
// If this is just a plane (degree (1, 1)) then the triangulation
// code will notice that, and not bother checking chord tols.
poly.UvTriangulateInto(sm, this);
} else {
// A surface with compound curvature. So we must overlay a
@ -820,7 +819,7 @@ void SShell::MakeFromTransformationOf(SShell *a, Vector t, Quaternion q) {
void SShell::MakeEdgesInto(SEdgeList *sel) {
SSurface *s;
for(s = surface.First(); s; s = surface.NextAfter(s)) {
s->MakeEdgesInto(this, sel, false);
s->MakeEdgesInto(this, sel, SSurface::AS_XYZ);
}
}

View File

@ -38,6 +38,7 @@ public:
SBspUv *InsertEdge(Point2d a, Point2d b, SSurface *srf);
int ClassifyPoint(Point2d p, Point2d eb, SSurface *srf);
int ClassifyEdge(Point2d ea, Point2d eb, SSurface *srf);
double MinimumDistanceToEdge(Point2d p, SSurface *srf);
};
// Now the data structures to represent a shell of trimmed rational polynomial
@ -224,6 +225,10 @@ public:
SBspUv *bsp;
SEdgeList edges;
// For caching our initial (u, v) when doing Newton iterations to project
// a point into our surface.
Point2d cached;
static SSurface FromExtrusionOf(SBezier *spc, Vector t0, Vector t1);
static SSurface FromRevolutionOf(SBezier *sb, Vector pt, Vector axis,
double thetas, double thetaf);
@ -266,6 +271,8 @@ public:
void ClosestPointTo(Vector p, Point2d *puv, bool converge=true);
void ClosestPointTo(Vector p, double *u, double *v, bool converge=true);
bool ClosestPointNewton(Vector p, double *u, double *v, bool converge=true);
bool PointIntersectingLine(Vector p0, Vector p1, double *u, double *v);
Vector ClosestPointOnThisAndSurface(SSurface *srf2, Vector p);
void PointOnSurfaces(SSurface *s1, SSurface *s2, double *u, double *v);
@ -283,9 +290,14 @@ public:
Vector *start, Vector *finish);
void TriangulateInto(SShell *shell, SMesh *sm);
void MakeTrimEdgesInto(SEdgeList *sel, bool asUv, SCurve *sc, STrimBy *stb);
void MakeEdgesInto(SShell *shell, SEdgeList *sel, bool asUv,
// these are intended as bitmasks, even though there's just one now
static const int AS_UV = 0x01;
static const int AS_XYZ = 0x00;
void MakeTrimEdgesInto(SEdgeList *sel, int flags, SCurve *sc, STrimBy *stb);
void MakeEdgesInto(SShell *shell, SEdgeList *sel, int flags,
SShell *useCurvesFrom=NULL);
Vector ExactSurfaceTangentAt(Vector p, SSurface *srfA, SSurface *srfB,
Vector dir);
void MakeSectionEdgesInto(SShell *shell, SEdgeList *sel, SBezierList *sbl);

View File

@ -54,13 +54,41 @@ void SSurface::AddExactIntersectionCurve(SBezier *sb, SSurface *srfB,
sc.Clear();
}
if(0 && sb->deg == 1) {
// Test if the curve lies entirely outside one of the
SCurvePt *scpt;
bool withinA = false, withinB = false;
for(scpt = split.pts.First(); scpt; scpt = split.pts.NextAfter(scpt)) {
double tol = 0.01;
Point2d puv;
ClosestPointTo(scpt->p, &puv);
if(puv.x > -tol && puv.x < 1 + tol &&
puv.y > -tol && puv.y < 1 + tol)
{
withinA = true;
}
srfB->ClosestPointTo(scpt->p, &puv);
if(puv.x > -tol && puv.x < 1 + tol &&
puv.y > -tol && puv.y < 1 + tol)
{
withinB = true;
}
// Break out early, no sense wasting time if we already have the answer.
if(withinA && withinB) break;
}
if(!(withinA && withinB)) {
// Intersection curve lies entirely outside one of the surfaces, so
// it's fake.
split.Clear();
return;
}
if(sb->deg == 2 && 0) {
dbp(" ");
SCurvePt *prev = NULL, *v;
dbp("split.pts.n =%d", split.pts.n);
dbp("split.pts.n = %d", split.pts.n);
for(v = split.pts.First(); v; v = split.pts.NextAfter(v)) {
if(prev) {
Vector e = (prev->p).Minus(v->p).WithMagnitude(-1);
Vector e = (prev->p).Minus(v->p).WithMagnitude(0);
SS.nakedEdges.AddEdge((prev->p).Plus(e), (v->p).Minus(e));
}
prev = v;
@ -292,7 +320,7 @@ void SSurface::IntersectAgainst(SSurface *b, SShell *agnstA, SShell *agnstB,
SEdgeList el;
ZERO(&el);
srfA->MakeEdgesInto(shA, &el, false, NULL);
srfA->MakeEdgesInto(shA, &el, AS_XYZ, NULL);
SEdge *se;
for(se = el.l.First(); se; se = el.l.NextAfter(se)) {
@ -480,7 +508,7 @@ void SShell::MakeCoincidentEdgesInto(SSurface *proto, bool sameNormal,
SSurface *ss;
for(ss = surface.First(); ss; ss = surface.NextAfter(ss)) {
if(proto->CoincidentWith(ss, sameNormal)) {
ss->MakeEdgesInto(this, el, false, useCurvesFrom);
ss->MakeEdgesInto(this, el, SSurface::AS_XYZ, useCurvesFrom);
}
}

View File

@ -202,7 +202,7 @@ haveEdge:
return true;
}
bool SContour::IsEar(int bp) {
bool SContour::IsEar(int bp, double scaledEps) {
int ap = WRAP(bp-1, l.n),
cp = WRAP(bp+1, l.n);
@ -219,7 +219,7 @@ bool SContour::IsEar(int bp) {
}
Vector n = Vector::From(0, 0, -1);
if((tr.Normal()).Dot(n) < LENGTH_EPS) {
if((tr.Normal()).Dot(n) < scaledEps) {
// This vertex is reflex, or between two collinear edges; either way,
// it's not an ear.
return false;
@ -251,7 +251,7 @@ bool SContour::IsEar(int bp) {
return true;
}
void SContour::ClipEarInto(SMesh *m, int bp) {
void SContour::ClipEarInto(SMesh *m, int bp, double scaledEps) {
int ap = WRAP(bp-1, l.n),
cp = WRAP(bp+1, l.n);
@ -260,7 +260,7 @@ void SContour::ClipEarInto(SMesh *m, int bp) {
tr.a = l.elem[ap].p;
tr.b = l.elem[bp].p;
tr.c = l.elem[cp].p;
if(tr.Normal().MagSquared() < LENGTH_EPS*LENGTH_EPS) {
if(tr.Normal().MagSquared() < scaledEps*scaledEps) {
// A vertex with more than two edges will cause us to generate
// zero-area triangles, which must be culled.
} else {
@ -278,8 +278,15 @@ void SContour::ClipEarInto(SMesh *m, int bp) {
}
void SContour::UvTriangulateInto(SMesh *m, SSurface *srf) {
int i;
Vector tu, tv;
srf->TangentsAt(0.5, 0.5, &tu, &tv);
double s = sqrt(tu.MagSquared() + tv.MagSquared());
// We would like to apply our tolerances in xyz; but that would be a lot
// of work, so at least scale the epsilon semi-reasonably. That's
// perfect for square planes, less perfect for anything else.
double scaledEps = LENGTH_EPS / s;
int i;
// Clean the original contour by removing any zero-length edges.
l.ClearTags();
for(i = 1; i < l.n; i++) {
@ -291,7 +298,7 @@ void SContour::UvTriangulateInto(SMesh *m, SSurface *srf) {
// Now calculate the ear-ness of each vertex
for(i = 0; i < l.n; i++) {
(l.elem[i]).ear = IsEar(i) ? SPoint::EAR : SPoint::NOT_EAR;
(l.elem[i]).ear = IsEar(i, scaledEps) ? SPoint::EAR : SPoint::NOT_EAR;
}
bool toggle = false;
@ -299,7 +306,8 @@ void SContour::UvTriangulateInto(SMesh *m, SSurface *srf) {
// Some points may have changed ear-ness, so recalculate
for(i = 0; i < l.n; i++) {
if(l.elem[i].ear == SPoint::UNKNOWN) {
(l.elem[i]).ear = IsEar(i) ? SPoint::EAR : SPoint::NOT_EAR;
(l.elem[i]).ear = IsEar(i, scaledEps) ?
SPoint::EAR : SPoint::NOT_EAR;
}
}
@ -312,7 +320,8 @@ void SContour::UvTriangulateInto(SMesh *m, SSurface *srf) {
for(i = 0; i < l.n; i++) {
int ear = WRAP(i+offset, l.n);
if(l.elem[ear].ear == SPoint::EAR) {
if(!srf) {
if(srf->degm == 1 && srf->degn == 1) {
// This is a plane; any ear is a good ear.
bestEar = ear;
break;
}
@ -322,7 +331,7 @@ void SContour::UvTriangulateInto(SMesh *m, SSurface *srf) {
Vector prev = l.elem[WRAP((i+offset-1), l.n)].p,
next = l.elem[WRAP((i+offset+1), l.n)].p;
double tol = srf->ChordToleranceForEdge(prev, next);
if(tol < bestChordTol - LENGTH_EPS) {
if(tol < bestChordTol - scaledEps) {
bestEar = ear;
bestChordTol = tol;
}
@ -335,10 +344,10 @@ void SContour::UvTriangulateInto(SMesh *m, SSurface *srf) {
dbp("couldn't find an ear! fail");
return;
}
ClipEarInto(m, bestEar);
ClipEarInto(m, bestEar, scaledEps);
}
ClipEarInto(m, 0); // add the last triangle
ClipEarInto(m, 0, scaledEps); // add the last triangle
}
double SSurface::ChordToleranceForEdge(Vector a, Vector b) {