zcy
/
solvespace
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Revamp the edge classification for Booleans. I no longer make a 

separate polygon of coincident (with same or opposite normal)
faces; I instead test all the edges against the other shell, and
have extended the classify-against-shell stuff to handle those
cases.

And the normals are now perturbed a bit numerically, to either side
of the edge, to distinguish tangency from a coincident surface.

This seems to work fairly well, although things still tend to fail
when the piecewise linear tolerance is too coarse.

[git-p4: depot-paths = "//depot/solvespace/": change = 1964]
solver
Jonathan Westhues 2009-06-03 19:59:40 -08:00
parent 24891c0141
commit ae35b3595c
8 changed files with 392 additions and 180 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
graphicswin.cpp
View File

@ -93,7 +93,7 @@ const GraphicsWindow::MenuEntry GraphicsWindow::menu[] = {
{ 1, "&Vertical\tV", MNU_VERTICAL, 'V', mCon },
{ 1, NULL, 0, NULL },
{ 1, "&On Point / Curve / Plane\tO", MNU_ON_ENTITY, 'O', mCon },
{ 1, "E&qual Length / Radius\tQ", MNU_EQUAL, 'Q', mCon },
{ 1, "E&qual Length / Radius / Angle\tQ", MNU_EQUAL, 'Q', mCon },
{ 1, "Length Ra&tio\tZ", MNU_RATIO, 'Z', mCon },
{ 1, "At &Midpoint\tM", MNU_AT_MIDPOINT, 'M', mCon },
{ 1, "S&ymmetric\tY", MNU_SYMMETRIC, 'Y', mCon },

7
solvespace.cpp
View File

@ -199,6 +199,13 @@ void SolveSpace::AfterNewFile(void) {
// and the naked edges
nakedEdges.Clear();
// GenerateAll() expects the view to be valid, because it uses that to
// fill in default values for extrusion depths etc. (which won't matter
// here, but just don't let it work on garbage)
SS.GW.offset = Vector::From(0, 0, 0);
SS.GW.projRight = Vector::From(1, 0, 0);
SS.GW.projUp = Vector::From(0, 1, 0);
ReloadAllImported();
GenerateAll(-1, -1);

251
srf/boolean.cpp
View File

@ -158,26 +158,12 @@ void SSurface::TrimFromEdgeList(SEdgeList *el) {
}
}
// For each edge, we record the membership of the regions to its left and
// right, which we call the "in direction" and "out direction" (wrt its
// outer normal)
#define INDIR (0)
#define OUTDIR (8)
// Regions of interest are the other shell itself, the coincident faces of the
// shell (same or opposite normal) and the original surface.
#define SHELL (0)
#define COINCIDENT_SAME (1)
#define COINCIDENT_OPPOSITE (2)
#define ORIG (3)
// Macro for building bit to test
#define INSIDE(reg, dir) (1 << ((reg)+(dir)))
static bool KeepRegion(int type, bool opA, int tag, int dir)
static bool KeepRegion(int type, bool opA, int shell, int orig)
{
bool inShell = (tag & INSIDE(SHELL, dir)) != 0,
inSame = (tag & INSIDE(COINCIDENT_SAME, dir)) != 0,
inOpp = (tag & INSIDE(COINCIDENT_OPPOSITE, dir)) != 0,
inOrig = (tag & INSIDE(ORIG, dir)) != 0;
bool inShell = (shell == SShell::INSIDE),
inSame = (shell == SShell::COINC_SAME),
inOpp = (shell == SShell::COINC_OPP),
inOrig = (orig == SShell::INSIDE);
bool inFace = inSame || inOpp;
@ -204,10 +190,12 @@ static bool KeepRegion(int type, bool opA, int tag, int dir)
default: oops();
}
}
static bool KeepEdge(int type, bool opA, int tag)
static bool KeepEdge(int type, bool opA,
int indir_shell, int outdir_shell,
int indir_orig, int outdir_orig)
{
bool keepIn = KeepRegion(type, opA, tag, INDIR),
keepOut = KeepRegion(type, opA, tag, OUTDIR);
bool keepIn = KeepRegion(type, opA, indir_shell, indir_orig),
keepOut = KeepRegion(type, opA, outdir_shell, outdir_orig);
// If the regions to the left and right of this edge are both in or both
// out, then this edge is not useful and should be discarded.
@ -215,37 +203,33 @@ static bool KeepEdge(int type, bool opA, int tag)
return false;
}
static int TagByClassifiedEdge(int bspclass, int reg)
static void TagByClassifiedEdge(int bspclass, int *indir, int *outdir)
{
switch(bspclass) {
case SBspUv::INSIDE:
return INSIDE(reg, OUTDIR) | INSIDE(reg, INDIR);
*indir = SShell::INSIDE;
*outdir = SShell::INSIDE;
break;
case SBspUv::OUTSIDE:
return 0;
*indir = SShell::OUTSIDE;
*outdir = SShell::OUTSIDE;
break;
case SBspUv::EDGE_PARALLEL:
return INSIDE(reg, INDIR);
*indir = SShell::INSIDE;
*outdir = SShell::OUTSIDE;
break;
case SBspUv::EDGE_ANTIPARALLEL:
return INSIDE(reg, OUTDIR);
*indir = SShell::OUTSIDE;
*outdir = SShell::INSIDE;
break;
default: oops();
}
}
void DBPEDGE(int tag) {
dbp("edge: indir %s %s %s %s ; outdir %s %s %s %s",
(tag & INSIDE(SHELL, INDIR)) ? "shell" : "",
(tag & INSIDE(COINCIDENT_SAME, INDIR)) ? "coinc-same" : "",
(tag & INSIDE(COINCIDENT_OPPOSITE, INDIR)) ? "coinc-opp" : "",
(tag & INSIDE(ORIG, INDIR)) ? "orig" : "",
(tag & INSIDE(SHELL, OUTDIR)) ? "shell" : "",
(tag & INSIDE(COINCIDENT_SAME, OUTDIR)) ? "coinc-same" : "",
(tag & INSIDE(COINCIDENT_OPPOSITE, OUTDIR)) ? "coinc-opp" : "",
(tag & INSIDE(ORIG, OUTDIR)) ? "orig" : "");
}
void DEBUGEDGELIST(SEdgeList *sel, SSurface *surf) {
dbp("print %d edges", sel->l.n);
SEdge *se;
@ -264,6 +248,54 @@ void DEBUGEDGELIST(SEdgeList *sel, SSurface *surf) {
}
}
static char *REGION(int d) {
switch(d) {
case SShell::INSIDE: return "inside";
case SShell::OUTSIDE: return "outside";
case SShell::COINC_SAME: return "same";
case SShell::COINC_OPP: return "opp";
default: return "xxx";
}
}
void SSurface::EdgeNormalsWithinSurface(Point2d auv, Point2d buv,
Vector *pt,
Vector *enin, Vector *enout,
Vector *surfn,
DWORD auxA, SShell *shell)
{
// 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);
*surfn = NormalAt(muv.x, muv.y);
// If this edge just approximates a curve, then refine our midpoint so
// so that it actually lies on that curve too. Otherwise stuff like
// point-on-face tests will fail, since the point won't actually lie
// on the other face.
hSCurve hc = { auxA };
SCurve *sc = shell->curve.FindById(hc);
if(sc->isExact && sc->exact.deg != 1) {
double t;
sc->exact.ClosestPointTo(*pt, &t, false);
*pt = sc->exact.PointAt(t);
}
// 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));
*enin = pin.Minus(*pt),
*enout = pout.Minus(*pt);
}
//-----------------------------------------------------------------------------
// Trim this surface against the specified shell, in the way that's appropriate
// for the specified Boolean operation type (and which operand we are). We
@ -294,7 +326,8 @@ SSurface SSurface::MakeCopyTrimAgainst(SShell *agnst, SShell *parent,
}
// Build up our original trim polygon; remember the coordinates could
// be changed if we just flipped the surface normal.
// be changed if we just flipped the surface normal, and we are using
// the split curves (not the original curves).
SEdgeList orig;
ZERO(&orig);
ret.MakeEdgesInto(into, &orig, true);
@ -302,22 +335,6 @@ SSurface SSurface::MakeCopyTrimAgainst(SShell *agnst, SShell *parent,
// which means that we can't necessarily use the old BSP...
SBspUv *origBsp = SBspUv::From(&orig);
// Find any surfaces from the other shell that are coincident with ours,
// and get the intersection polygons, grouping surfaces with the same
// and with opposite normal separately.
SEdgeList sameNormal, oppositeNormal;
ZERO(&sameNormal);
ZERO(&oppositeNormal);
agnst->MakeCoincidentEdgesInto(&ret, true, &sameNormal, into);
agnst->MakeCoincidentEdgesInto(&ret, false, &oppositeNormal, into);
// and cull parallel or anti-parallel pairs, which may occur if multiple
// surfaces are coincident with ours
sameNormal.CullExtraneousEdges();
oppositeNormal.CullExtraneousEdges();
// and build the trees for quick in-polygon testing
SBspUv *sameBsp = SBspUv::From(&sameNormal);
SBspUv *oppositeBsp = SBspUv::From(&oppositeNormal);
// And now intersect the other shell against us
SEdgeList inter;
ZERO(&inter);
@ -373,45 +390,25 @@ SSurface SSurface::MakeCopyTrimAgainst(SShell *agnst, SShell *parent,
SEdge *se;
for(se = orig.l.First(); se; se = orig.l.NextAfter(se)) {
Point2d auv = (se->a).ProjectXy(),
buv = (se->b).ProjectXy();
Point2d auv = (se->a).ProjectXy(),
buv = (se->b).ProjectXy();
int c_same = sameBsp->ClassifyEdge(auv, buv);
int c_opp = oppositeBsp->ClassifyEdge(auv, buv);
Vector pt, enin, enout, surfn;
ret.EdgeNormalsWithinSurface(auv, buv, &pt, &enin, &enout, &surfn,
se->auxA, into);
// Get the midpoint of this edge
Point2d am = (auv.Plus(buv)).ScaledBy(0.5);
Vector pt = ret.PointAt(am.x, am.y);
// and the outer normal from the trim polygon (within the surface)
Vector surf_n = ret.NormalAt(am.x, am.y);
Vector ea = ret.PointAt(auv.x, auv.y),
eb = ret.PointAt(buv.x, buv.y);
Vector edge_n = surf_n.Cross((eb.Minus(ea)));
int indir_shell, outdir_shell, indir_orig, outdir_orig;
int c_shell = agnst->ClassifyPoint(pt, edge_n, surf_n);
indir_orig = SShell::INSIDE;
outdir_orig = SShell::OUTSIDE;
int tag = 0;
tag |= INSIDE(ORIG, INDIR);
tag |= TagByClassifiedEdge(c_same, COINCIDENT_SAME);
tag |= TagByClassifiedEdge(c_opp, COINCIDENT_OPPOSITE);
agnst->ClassifyEdge(&indir_shell, &outdir_shell,
ret.PointAt(auv), ret.PointAt(buv), pt,
enin, enout, surfn);
if(c_shell == SShell::INSIDE) {
tag |= INSIDE(SHELL, INDIR) | INSIDE(SHELL, OUTDIR);
} else if(c_shell == SShell::OUTSIDE) {
tag |= 0;
} else if(c_shell == SShell::SURF_PARALLEL) {
tag |= INSIDE(SHELL, INDIR);
} else if(c_shell == SShell::SURF_ANTIPARALLEL) {
tag |= INSIDE(SHELL, OUTDIR);
} else if(c_shell == SShell::EDGE_PARALLEL) {
tag |= INSIDE(SHELL, INDIR);
} else if(c_shell == SShell::EDGE_ANTIPARALLEL) {
tag |= INSIDE(SHELL, OUTDIR);
} else if(c_shell == SShell::EDGE_TANGENT) {
continue;
}
if(KeepEdge(type, opA, tag)) {
if(KeepEdge(type, opA, indir_shell, outdir_shell,
indir_orig, outdir_orig))
{
final.AddEdge(se->a, se->b, se->auxA, se->auxB);
}
}
@ -420,23 +417,22 @@ SSurface SSurface::MakeCopyTrimAgainst(SShell *agnst, SShell *parent,
Point2d auv = (se->a).ProjectXy(),
buv = (se->b).ProjectXy();
Vector pt, enin, enout, surfn;
ret.EdgeNormalsWithinSurface(auv, buv, &pt, &enin, &enout, &surfn,
se->auxA, into);
int indir_shell, outdir_shell, indir_orig, outdir_orig;
int c_this = origBsp->ClassifyEdge(auv, buv);
int c_same = sameBsp->ClassifyEdge(auv, buv);
int c_opp = oppositeBsp->ClassifyEdge(auv, buv);
TagByClassifiedEdge(c_this, &indir_orig, &outdir_orig);
int tag = 0;
tag |= TagByClassifiedEdge(c_this, ORIG);
tag |= TagByClassifiedEdge(c_same, COINCIDENT_SAME);
tag |= TagByClassifiedEdge(c_opp, COINCIDENT_OPPOSITE);
agnst->ClassifyEdge(&indir_shell, &outdir_shell,
ret.PointAt(auv), ret.PointAt(buv), pt,
enin, enout, surfn);
if(type == SShell::AS_DIFFERENCE && !opA) {
// The second operand of a difference gets turned inside out
tag |= INSIDE(SHELL, INDIR);
} else {
tag |= INSIDE(SHELL, OUTDIR);
}
if(KeepEdge(type, opA, tag)) {
if(KeepEdge(type, opA, indir_shell, outdir_shell,
indir_orig, outdir_orig))
{
final.AddEdge(se->a, se->b, se->auxA, se->auxB);
}
}
@ -455,12 +451,11 @@ SSurface SSurface::MakeCopyTrimAgainst(SShell *agnst, SShell *parent,
final.l.ClearTags();
if(!final.AssemblePolygon(&poly, NULL, true)) {
into->booleanFailed = true;
dbp("failed: I=%d", I);
DEBUGEDGELIST(&final, &ret);
}
poly.Clear();
sameNormal.Clear();
oppositeNormal.Clear();
final.Clear();
inter.Clear();
orig.Clear();
@ -496,6 +491,7 @@ void SShell::MakeIntersectionCurvesAgainst(SShell *agnst, SShell *into) {
void SShell::CleanupAfterBoolean(void) {
SSurface *ss;
for(ss = surface.First(); ss; ss = surface.NextAfter(ss)) {
ss->edges.Clear();
}
}
@ -576,8 +572,8 @@ void SShell::MakeFromAssemblyOf(SShell *a, SShell *b) {
void SShell::MakeFromBoolean(SShell *a, SShell *b, int type) {
booleanFailed = false;
a->MakeClassifyingBsps();
b->MakeClassifyingBsps();
a->MakeClassifyingBsps(NULL);
b->MakeClassifyingBsps(NULL);
// Copy over all the original curves, splitting them so that a
// piecwise linear segment never crosses a surface from the other
@ -605,11 +601,21 @@ void SShell::MakeFromBoolean(SShell *a, SShell *b, int type) {
sc->RemoveShortSegments(srfA, srfB);
}
// And clean up the piecewise linear things we made as a calculation aid
a->CleanupAfterBoolean();
b->CleanupAfterBoolean();
// Remake the classifying BSPs with the split (and short-segment-removed)
// curves
a->MakeClassifyingBsps(this);
b->MakeClassifyingBsps(this);
// Then trim and copy the surfaces
if(b->surface.n == 0 || a->surface.n == 0) {
// Then trim and copy the surfaces
I = 1023123;
a->CopySurfacesTrimAgainst(b, this, type, true);
b->CopySurfacesTrimAgainst(a, this, type, false);
} else {
I = 0;
a->CopySurfacesTrimAgainst(b, this, type, true);
b->CopySurfacesTrimAgainst(a, this, type, false);
}
@ -627,21 +633,23 @@ void SShell::MakeFromBoolean(SShell *a, SShell *b, int type) {
//-----------------------------------------------------------------------------
// All of the BSP routines that we use to perform and accelerate polygon ops.
//-----------------------------------------------------------------------------
void SShell::MakeClassifyingBsps(void) {
void SShell::MakeClassifyingBsps(SShell *useCurvesFrom) {
SSurface *ss;
for(ss = surface.First(); ss; ss = surface.NextAfter(ss)) {
ss->MakeClassifyingBsp(this);
ss->MakeClassifyingBsp(this, useCurvesFrom);
}
}
void SSurface::MakeClassifyingBsp(SShell *shell) {
void SSurface::MakeClassifyingBsp(SShell *shell, SShell *useCurvesFrom) {
SEdgeList el;
ZERO(&el);
MakeEdgesInto(shell, &el, true);
MakeEdgesInto(shell, &el, true, useCurvesFrom);
bsp = SBspUv::From(&el);
el.Clear();
ZERO(&edges);
MakeEdgesInto(shell, &edges, false, useCurvesFrom);
}
SBspUv *SBspUv::Alloc(void) {
@ -733,7 +741,7 @@ SBspUv *SBspUv::InsertEdge(Point2d ea, Point2d eb) {
return this;
}
int SBspUv::ClassifyPoint(Point2d p, Point2d eb) {
int SBspUv::ClassifyPoint(Point2d p, Point2d eb, Point2d *ia, Point2d *ib) {
if(!this) return OUTSIDE;
Point2d n = ((b.Minus(a)).Normal()).WithMagnitude(1);
@ -746,6 +754,11 @@ int SBspUv::ClassifyPoint(Point2d p, Point2d eb) {
while(f) {
Point2d ba = (f->b).Minus(f->a);
if(p.DistanceToLine(f->a, ba, true) < LENGTH_EPS) {
if(ia && ib) {
// Tell the caller which edge it happens to lie on.
*ia = f->a;
*ib = f->b;
}
if(eb.DistanceToLine(f->a, ba, false) < LENGTH_EPS) {
if(ba.Dot(eb.Minus(p)) > 0) {
return EDGE_PARALLEL;
@ -759,16 +772,16 @@ int SBspUv::ClassifyPoint(Point2d p, Point2d eb) {
f = f->more;
}
// Pick arbitrarily which side to send it down, doesn't matter
int c1 = neg ? neg->ClassifyPoint(p, eb) : OUTSIDE;
int c2 = pos ? pos->ClassifyPoint(p, eb) : INSIDE;
int c1 = neg ? neg->ClassifyPoint(p, eb, ia, ib) : OUTSIDE;
int c2 = pos ? pos->ClassifyPoint(p, eb, ia, ib) : INSIDE;
if(c1 != c2) {
dbp("MISMATCH: %d %d %08x %08x", c1, c2, neg, pos);
}
return c1;
} else if(dp > 0) {
return pos ? pos->ClassifyPoint(p, eb) : INSIDE;
return pos ? pos->ClassifyPoint(p, eb, ia, ib) : INSIDE;
} else {
return neg ? neg->ClassifyPoint(p, eb) : OUTSIDE;
return neg ? neg->ClassifyPoint(p, eb, ia, ib) : OUTSIDE;
}
}

59
srf/ratpoly.cpp
View File

@ -130,7 +130,7 @@ Vector SBezier::TangentAt(double t) {
return ret;
}
void SBezier::ClosestPointTo(Vector p, double *t) {
void SBezier::ClosestPointTo(Vector p, double *t, bool converge) {
int i;
double minDist = VERY_POSITIVE;
*t = 0;
@ -147,7 +147,7 @@ void SBezier::ClosestPointTo(Vector p, double *t) {
}
Vector p0;
for(i = 0; i < 15; i++) {
for(i = 0; i < (converge ? 15 : 3); i++) {
p0 = PointAt(*t);
if(p0.Equals(p, RATPOLY_EPS)) {
return;
@ -157,7 +157,9 @@ void SBezier::ClosestPointTo(Vector p, double *t) {
Vector pc = p.ClosestPointOnLine(p0, dp);
*t += (pc.Minus(p0)).DivPivoting(dp);
}
dbp("didn't converge (closest point on bezier curve)");
if(converge) {
dbp("didn't converge (closest point on bezier curve)");
}
}
void SBezier::SplitAt(double t, SBezier *bef, SBezier *aft) {
@ -242,6 +244,9 @@ void SBezier::MakePwlWorker(List<Vector> *l, double ta, double tb) {
}
}
Vector SSurface::PointAt(Point2d puv) {
return PointAt(puv.x, puv.y);
}
Vector SSurface::PointAt(double u, double v) {
Vector num = Vector::From(0, 0, 0);
double den = 0;
@ -294,12 +299,18 @@ void SSurface::TangentsAt(double u, double v, Vector *tu, Vector *tv) {
*tv = tv->ScaledBy(1.0/(den*den));
}
Vector SSurface::NormalAt(Point2d puv) {
return NormalAt(puv.x, puv.y);
}
Vector SSurface::NormalAt(double u, double v) {
Vector tu, tv;
TangentsAt(u, v, &tu, &tv);
return tu.Cross(tv);
}
void SSurface::ClosestPointTo(Vector p, Point2d *puv, bool converge) {
ClosestPointTo(p, &(puv->x), &(puv->y), converge);
}
void SSurface::ClosestPointTo(Vector p, double *u, double *v, bool converge) {
// A few special cases first; when control points are coincident the
// derivative goes to zero at the conrol points, and would result in
@ -394,6 +405,48 @@ bool SSurface::PointIntersectingLine(Vector p0, Vector p1, double *u, double *v)
return false;
}
Vector SSurface::ClosestPointOnThisAndSurface(SSurface *srf2, Vector p) {
// This is untested.
int i, j;
Point2d puv[2];
SSurface *srf[2] = { this, srf2 };
for(j = 0; j < 2; j++) {
(srf[j])->ClosestPointTo(p, &(puv[j]), false);
}
for(i = 0; i < 4; i++) {
Vector tu[2], tv[2], cp[2], n[2];
double d[2];
for(j = 0; j < 2; j++) {
(srf[j])->TangentsAt(puv[j].x, puv[j].y, &(tu[j]), &(tv[j]));
cp[j] = (srf[j])->PointAt(puv[j]);
n[j] = ((tu[j]).Cross(tv[j])).WithMagnitude(1);
d[j] = (n[j]).Dot(cp[j]);
}
Vector p0 = Vector::AtIntersectionOfPlanes(n[0], d[0], n[1], d[1]),
dp = (n[0]).Cross(n[1]);
Vector pc = p.ClosestPointOnLine(p0, dp);
// Adjust our guess and iterate
for(j = 0; j < 2; j++) {
Vector dc = pc.Minus(cp[j]);
double du = dc.Dot(tu[j]), dv = dc.Dot(tv[j]);
puv[j].x += du / ((tu[j]).MagSquared());
puv[j].y += dv / ((tv[j]).MagSquared());
}
}
// If this converged, then the two points are actually equal.
return ((srf[0])->PointAt(puv[0])).Plus(
((srf[1])->PointAt(puv[1]))).ScaledBy(0.5);
}
void SSurface::PointOnSurfaces(SSurface *s1, SSurface *s2,
double *up, double *vp)
{

48
srf/surface.h
View File

@ -31,7 +31,8 @@ public:
Point2d IntersectionWith(Point2d a, Point2d b);
SBspUv *InsertEdge(Point2d a, Point2d b);
int ClassifyPoint(Point2d p, Point2d eb);
int ClassifyPoint(Point2d p, Point2d eb,
Point2d *ia=NULL, Point2d *ib=NULL);
int ClassifyEdge(Point2d ea, Point2d eb);
};
@ -62,7 +63,7 @@ public:
Vector PointAt(double t);
Vector TangentAt(double t);
void ClosestPointTo(Vector p, double *t);
void ClosestPointTo(Vector p, double *t, bool converge=true);
void SplitAt(double t, SBezier *bef, SBezier *aft);
Vector Start(void);
@ -184,9 +185,10 @@ public:
int tag;
Vector p;
SSurface *srf;
hSSurface hsrf;
Vector surfNormal; // of the intersecting surface, at pinter
bool onEdge; // pinter is on edge of trim poly
Point2d pinter;
Vector surfNormal; // of the intersecting surface, at pinter
bool onEdge; // pinter is on edge of trim poly
Point2d edgeA, edgeB; // the edge that pinter is on
};
// A rational polynomial surface in Bezier form.
@ -209,6 +211,7 @@ public:
// For testing whether a point (u, v) on the surface lies inside the trim
SBspUv *bsp;
SEdgeList edges;
static SSurface FromExtrusionOf(SBezier *spc, Vector t0, Vector t1);
static SSurface FromRevolutionOf(SBezier *sb, Vector pt, Vector axis,
@ -217,6 +220,10 @@ public:
static SSurface FromTransformationOf(SSurface *a, Vector t, Quaternion q,
bool includingTrims);
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);
void TrimFromEdgeList(SEdgeList *el);
@ -244,11 +251,15 @@ public:
List<Inter> *l, bool segment,
SSurface *sorig);
void ClosestPointTo(Vector p, Point2d *puv, bool converge=true);
void ClosestPointTo(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);
Vector PointAt(double u, double v);
Vector PointAt(Point2d puv);
void TangentsAt(double u, double v, Vector *tu, Vector *tv);
Vector NormalAt(Point2d puv);
Vector NormalAt(double u, double v);
bool LineEntirelyOutsideBbox(Vector a, Vector b, bool segment);
void GetAxisAlignedBounding(Vector *ptMax, Vector *ptMin);
@ -263,7 +274,7 @@ public:
void MakeEdgesInto(SShell *shell, SEdgeList *sel, bool asUv,
SShell *useCurvesFrom=NULL);
void MakeSectionEdgesInto(SShell *shell, SEdgeList *sel, SBezierList *sbl);
void MakeClassifyingBsp(SShell *shell);
void MakeClassifyingBsp(SShell *shell, SShell *useCurvesFrom);
double ChordToleranceForEdge(Vector a, Vector b);
void MakeTriangulationGridInto(List<double> *l, double vs, double vf,
bool swapped);
@ -294,7 +305,7 @@ public:
void CopyCurvesSplitAgainst(bool opA, SShell *agnst, SShell *into);
void CopySurfacesTrimAgainst(SShell *against, SShell *into, int t, bool a);
void MakeIntersectionCurvesAgainst(SShell *against, SShell *into);
void MakeClassifyingBsps(void);
void MakeClassifyingBsps(SShell *useCurvesFrom);
void AllPointsIntersecting(Vector a, Vector b, List<SInter> *il,
bool seg, bool trimmed, bool inclTangent);
void MakeCoincidentEdgesInto(SSurface *proto, bool sameNormal,
@ -302,16 +313,19 @@ public:
void RewriteSurfaceHandlesForCurves(SShell *a, SShell *b);
void CleanupAfterBoolean(void);
static const int INSIDE = 100;
static const int OUTSIDE = 200;
static const int SURF_PARALLEL = 300;
static const int SURF_ANTIPARALLEL = 400;
static const int EDGE_PARALLEL = 500;
static const int EDGE_ANTIPARALLEL = 600;
static const int EDGE_TANGENT = 700;
int ClassifyPoint(Vector p, Vector edge_n, Vector surf_n);
// Definitions when classifying regions of a surface; it is either inside,
// outside, or coincident (with parallel or antiparallel normal) with a
// shell.
static const int INSIDE = 100;
static const int OUTSIDE = 200;
static const int COINC_SAME = 300;
static const int COINC_OPP = 400;
static const double DOTP_TOL;
int ClassifyRegion(Vector edge_n, Vector inter_surf_n, Vector edge_surf_n);
bool ClassifyEdge(int *indir, int *outdir,
Vector ea, Vector eb,
Vector p,
Vector edge_n_in, Vector edge_n_out, Vector surf_n);
void MakeFromCopyOf(SShell *a);
void MakeFromTransformationOf(SShell *a, Vector trans, Quaternion q);

197
srf/surfinter.cpp
View File

@ -1,5 +1,8 @@
#include "solvespace.h"
// Dot product tolerance for perpendicular.
const double SShell::DOTP_TOL = 1e-3;
extern int FLAG;
void SSurface::AddExactIntersectionCurve(SBezier *sb, SSurface *srfB,
@ -220,9 +223,25 @@ void SSurface::IntersectAgainst(SSurface *b, SShell *agnstA, SShell *agnstB,
List<Vector> lv;
ZERO(&lv);
Vector axisa = axis.ScaledBy((b->ctrl[0][0]).Dot(axis)),
axisb = axis.ScaledBy((b->ctrl[0][1]).Dot(axis)),
axisc = (axisa.Plus(axisb)).ScaledBy(0.5);
double a_axis0 = ( ctrl[0][0]).Dot(axis),
a_axis1 = ( ctrl[0][1]).Dot(axis),
b_axis0 = (b->ctrl[0][0]).Dot(axis),
b_axis1 = (b->ctrl[0][1]).Dot(axis);
if(a_axis0 > a_axis1) SWAP(double, a_axis0, a_axis1);
if(b_axis0 > b_axis1) SWAP(double, b_axis0, b_axis1);
double ab_axis0 = max(a_axis0, b_axis0),
ab_axis1 = min(a_axis1, b_axis1);
if(fabs(ab_axis0 - ab_axis1) < LENGTH_EPS) {
// The line would be zero-length
return;
}
Vector axis0 = axis.ScaledBy(ab_axis0),
axis1 = axis.ScaledBy(ab_axis1),
axisc = (axis0.Plus(axis1)).ScaledBy(0.5);
oft.MakePwlInto(&lv);
@ -250,7 +269,7 @@ void SSurface::IntersectAgainst(SSurface *b, SShell *agnstA, SShell *agnstB,
p = b->PointAt(ub, vb);
SBezier bezier;
bezier = SBezier::From(p.Plus(axisa), p.Plus(axisb));
bezier = SBezier::From(p.Plus(axis0), p.Plus(axis1));
AddExactIntersectionCurve(&bezier, b, agnstA, agnstB, into);
}
@ -608,8 +627,8 @@ void SSurface::AllPointsIntersecting(Vector a, Vector b,
}
// And that it lies inside our trim region
Point2d dummy = { 0, 0 };
int c = bsp->ClassifyPoint(puv, dummy);
Point2d dummy = { 0, 0 }, ia = { 0, 0 }, ib = { 0, 0 };
int c = bsp->ClassifyPoint(puv, dummy, &ia, &ib);
if(trimmed && c == SBspUv::OUTSIDE) {
continue;
}
@ -618,9 +637,11 @@ void SSurface::AllPointsIntersecting(Vector a, Vector b,
SInter si;
si.p = pxyz;
si.surfNormal = NormalAt(puv.x, puv.y);
si.hsrf = h;
si.pinter = puv;
si.srf = this;
si.onEdge = (c != SBspUv::INSIDE);
si.edgeA = ia;
si.edgeB = ib;
l->Add(&si);
}
@ -637,6 +658,26 @@ void SShell::AllPointsIntersecting(Vector a, Vector b,
}
}
int SShell::ClassifyRegion(Vector edge_n, Vector inter_surf_n,
Vector edge_surf_n)
{
double dot = inter_surf_n.Dot(edge_n);
if(fabs(dot) < DOTP_TOL) {
// The edge's surface and the edge-on-face surface
// are coincident. Test the edge's surface normal
// to see if it's with same or opposite normals.
if(inter_surf_n.Dot(edge_surf_n) > 0) {
return COINC_SAME;
} else {
return COINC_OPP;
}
} else if(dot > 0) {
return OUTSIDE;
} else {
return INSIDE;
}
}
//-----------------------------------------------------------------------------
// Does the given point lie on our shell? There are many cases; inside and
// outside are obvious, but then there's all the edge-on-edge and edge-on-face
@ -646,15 +687,101 @@ void SShell::AllPointsIntersecting(Vector a, Vector b,
// using the closest intersection point. If the ray hits a surface on edge,
// then just reattempt in a different random direction.
//-----------------------------------------------------------------------------
int SShell::ClassifyPoint(Vector p, Vector edge_n, Vector surf_n) {
bool SShell::ClassifyEdge(int *indir, int *outdir,
Vector ea, Vector eb,
Vector p,
Vector edge_n_in, Vector edge_n_out, Vector surf_n)
{
List<SInter> l;
ZERO(&l);
srand(0);
int ret, cnt = 0, edge_inters;
double edge_dotp[2];
// First, check for edge-on-edge
int edge_inters = 0;
Vector inter_surf_n[2], inter_edge_n[2];
SSurface *srf;
for(srf = surface.First(); srf; srf = surface.NextAfter(srf)) {
if(srf->LineEntirelyOutsideBbox(ea, eb, true)) continue;
SEdgeList *sel = &(srf->edges);
SEdge *se;
for(se = sel->l.First(); se; se = sel->l.NextAfter(se)) {
if((ea.Equals(se->a) && eb.Equals(se->b)) ||
(eb.Equals(se->a) && ea.Equals(se->b)) ||
p.OnLineSegment(se->a, se->b))
{
if(edge_inters < 2) {
// Edge-on-edge case
Point2d pm;
srf->ClosestPointTo(p, &pm, false);
// A vector normal to the surface, at the intersection point
inter_surf_n[edge_inters] = srf->NormalAt(pm);
// A vector normal to the intersecting edge (but within the
// intersecting surface) at the intersection point, pointing
// out.
inter_edge_n[edge_inters] =
(inter_surf_n[edge_inters]).Cross((se->b).Minus((se->a)));
}
edge_inters++;
}
}
}
if(edge_inters == 2) {
// TODO, make this use the appropriate curved normals
double dotp[2];
for(int i = 0; i < 2; i++) {
dotp[i] = edge_n_out.Dot(inter_surf_n[i]);
}
if(fabs(dotp[1]) < DOTP_TOL) {
SWAP(double, dotp[0], dotp[1]);
SWAP(Vector, inter_surf_n[0], inter_surf_n[1]);
SWAP(Vector, inter_edge_n[0], inter_edge_n[1]);
}
int coinc = (surf_n.Dot(inter_surf_n[0])) > 0 ? COINC_SAME : COINC_OPP;
if(fabs(dotp[0]) < DOTP_TOL && fabs(dotp[1]) < DOTP_TOL) {
// This is actually an edge on face case, just that the face
// is split into two pieces joining at our edge.
*indir = coinc;
*outdir = coinc;
} else if(fabs(dotp[0]) < DOTP_TOL && dotp[1] > DOTP_TOL) {
if(edge_n_out.Dot(inter_edge_n[0]) > 0) {
*indir = coinc;
*outdir = OUTSIDE;
} else {
*indir = INSIDE;
*outdir = coinc;
}
} else if(fabs(dotp[0]) < DOTP_TOL && dotp[1] < -DOTP_TOL) {
if(edge_n_out.Dot(inter_edge_n[0]) > 0) {
*indir = coinc;
*outdir = INSIDE;
} else {
*indir = OUTSIDE;
*outdir = coinc;
}
} else if(dotp[0] > DOTP_TOL && dotp[1] > DOTP_TOL) {
*indir = INSIDE;
*outdir = OUTSIDE;
} else if(dotp[0] < -DOTP_TOL && dotp[1] < -DOTP_TOL) {
*indir = OUTSIDE;
*outdir = INSIDE;
} else {
// Edge is tangent to the shell at shell's edge, so can't be
// a boundary of the surface.
return false;
}
return true;
}
if(edge_inters != 0) dbp("bad, edge_inters=%d", edge_inters);
int cnt = 0;
for(;;) {
// Cast a ray in a random direction (two-sided so that we test if
// the point lies on a surface, but use only one side for in/out
@ -663,8 +790,10 @@ int SShell::ClassifyPoint(Vector p, Vector edge_n, Vector surf_n) {
AllPointsIntersecting(
p.Minus(ray), p.Plus(ray), &l, false, true, false);
// no intersections means it's outside
*indir = OUTSIDE;
*outdir = OUTSIDE;
double dmin = VERY_POSITIVE;
ret = OUTSIDE; // no intersections means it's outside
bool onEdge = false;
edge_inters = 0;
@ -678,9 +807,9 @@ int SShell::ClassifyPoint(Vector p, Vector edge_n, Vector surf_n) {
double d = ((si->p).Minus(p)).Magnitude();
// Handle edge-on-edge
if(d < LENGTH_EPS && si->onEdge && edge_inters < 2) {
edge_dotp[edge_inters] = (si->surfNormal).Dot(edge_n);
// We actually should never hit this case; it should have been
// handled above.
if(d < LENGTH_EPS && si->onEdge) {
edge_inters++;
}
@ -689,21 +818,24 @@ int SShell::ClassifyPoint(Vector p, Vector edge_n, Vector surf_n) {
if(d < LENGTH_EPS) {
// Edge-on-face (unless edge-on-edge above supercedes)
double dot = (si->surfNormal).Dot(edge_n);
if(fabs(dot) < LENGTH_EPS && 0) {
// TODO revamp this
} else if(dot > 0) {
ret = SURF_PARALLEL;
} else {
ret = SURF_ANTIPARALLEL;
}
Point2d pin, pout;
(si->srf)->ClosestPointTo(p.Plus(edge_n_in), &pin, false);
(si->srf)->ClosestPointTo(p.Plus(edge_n_out), &pout, false);
Vector surf_n_in = (si->srf)->NormalAt(pin),
surf_n_out = (si->srf)->NormalAt(pout);
*indir = ClassifyRegion(edge_n_in, surf_n_in, surf_n);
*outdir = ClassifyRegion(edge_n_out, surf_n_out, surf_n);
} else {
// Edge does not lie on surface; either strictly inside
// or strictly outside
if((si->surfNormal).Dot(ray) > 0) {
ret = INSIDE;
*indir = INSIDE;
*outdir = INSIDE;
} else {
ret = OUTSIDE;
*indir = OUTSIDE;
*outdir = OUTSIDE;
}
}
onEdge = si->onEdge;
@ -714,27 +846,16 @@ int SShell::ClassifyPoint(Vector p, Vector edge_n, Vector surf_n) {
// If the point being tested lies exactly on an edge of the shell,
// then our ray always lies on edge, and that's okay. Otherwise
// try again in a different random direction.
if((edge_inters == 2) || !onEdge) break;
if(!onEdge) break;
if(cnt++ > 5) {
dbp("can't find a ray that doesn't hit on edge!");
dbp("on edge = %d, edge_inters = %d", onEdge, edge_inters);
SS.nakedEdges.AddEdge(ea, eb);
break;
}
}
if(edge_inters == 2) {
double tol = 1e-3;
if(edge_dotp[0] > -tol && edge_dotp[1] > -tol) {
return EDGE_PARALLEL;
} else if(edge_dotp[0] < tol && edge_dotp[1] < tol) {
return EDGE_ANTIPARALLEL;
} else {
return EDGE_TANGENT;
}
} else {
return ret;
}
return true;
}
//-----------------------------------------------------------------------------

5
util.cpp
View File

@ -460,6 +460,8 @@ double Vector::DistanceToLine(Vector p0, Vector dp) {
}
bool Vector::OnLineSegment(Vector a, Vector b) {
if(this->Equals(a) || this->Equals(b)) return true;
Vector d = b.Minus(a);
double m = d.MagSquared();
@ -468,7 +470,7 @@ bool Vector::OnLineSegment(Vector a, Vector b) {
if(distsq >= LENGTH_EPS*LENGTH_EPS) return false;
double t = (this->Minus(a)).DivPivoting(d);
// On-endpoint must be tested for separately.
// On-endpoint already tested
if(t < 0 || t > 1) return false;
return true;
}
@ -508,6 +510,7 @@ Vector Vector::ScaledBy(double v) {
Vector Vector::WithMagnitude(double v) {
double m = Magnitude();
if(m == 0) {
dbp("Vector::WithMagnitude of zero vector!");
return From(v, 0, 0);
} else {
return ScaledBy(v/m);

3
wishlist.txt
View File

@ -1,9 +1,10 @@
marching algorithm for surface intersection
tangent intersections
plane detection for solid of revolution
-----
line styles (color, thickness)
marching algorithm for surface intersection
loop detection
IGES and STEP export
incremental regen of entities?