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

Add O(n*log(n)) self-intersection check for polygons. That was 

previously a major bottleneck, and is now at least 10x faster for a
practical section.

And fix a horrible uninit memory bug in the triangle kd-tree stuff;
the split planes were apparently random. This would have slowed
things down, but not caused an incorrect result; except when it
ends up NaN, which is the reason I noticed.

[git-p4: depot-paths = "//depot/solvespace/": change = 2073]
solver
Jonathan Westhues 2009-11-09 03:51:38 -08:00
parent f3eeae673f
commit 61da5969f6
4 changed files with 253 additions and 90 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

25
mesh.cpp
View File

@ -347,8 +347,6 @@ DWORD SMesh::FirstIntersectionWith(Point2d mp) {
return face;
}
#define KDTREE_EPS (20*LENGTH_EPS) // nice and sloppy
STriangleLl *STriangleLl::Alloc(void)
{ return (STriangleLl *)AllocTemporary(sizeof(STriangleLl)); }
SKdNode *SKdNode::Alloc(void)
@ -378,18 +376,21 @@ SKdNode *SKdNode::From(SMesh *m) {
tll = tn;
}
return SKdNode::From(tll, 0);
return SKdNode::From(tll);
}
SKdNode *SKdNode::From(STriangleLl *tll, int which) {
SKdNode *SKdNode::From(STriangleLl *tll) {
int which = 0;
SKdNode *ret = Alloc();
if(!tll) goto leaf;
if(!tll) {
goto leaf;
}
int i;
int gtc[3] = { 0, 0, 0 }, ltc[3] = { 0, 0, 0 }, allc = 0;
double badness[3];
double split[3];
double badness[3] = { 0, 0, 0 };
double split[3] = { 0, 0, 0 };
for(i = 0; i < 3; i++) {
int tcnt = 0;
STriangleLl *ll;
@ -433,8 +434,9 @@ SKdNode *SKdNode::From(STriangleLl *tll, int which) {
which = 2;
}
if(allc < 10) goto leaf;
if(allc == gtc[which] || allc == ltc[which]) goto leaf;
if(allc < 3 || allc == gtc[which] || allc == ltc[which]) {
goto leaf;
}
STriangleLl *ll;
STriangleLl *lgt = NULL, *llt = NULL;
@ -467,12 +469,11 @@ SKdNode *SKdNode::From(STriangleLl *tll, int which) {
ret->which = which;
ret->c = split[which];
ret->gt = SKdNode::From(lgt, (which + 1) % 3);
ret->lt = SKdNode::From(llt, (which + 1) % 3);
ret->gt = SKdNode::From(lgt);
ret->lt = SKdNode::From(llt);
return ret;
leaf:
// dbp("leaf: allc=%d gtc=%d ltc=%d which=%d", allc, gtc[which], ltc[which], which);
ret->tris = tll;
return ret;
}

276
polygon.cpp
View File

@ -57,6 +57,75 @@ SEdge SEdge::From(Vector a, Vector b) {
return se;
}
bool SEdge::EdgeCrosses(Vector ea, Vector eb, Vector *ppi, SPointList *spl) {
Vector d = eb.Minus(ea);
double t_eps = LENGTH_EPS/d.Magnitude();
double dist_a, dist_b;
double t, tthis;
bool skew;
Vector pi;
bool inOrEdge0, inOrEdge1;
Vector dthis = b.Minus(a);
double tthis_eps = LENGTH_EPS/dthis.Magnitude();
if(ea.Equals(a) && eb.Equals(b)) return true;
if(eb.Equals(a) && ea.Equals(b)) return true;
dist_a = a.DistanceToLine(ea, d),
dist_b = b.DistanceToLine(ea, d);
// Can't just test if dist_a equals dist_b; they could be on opposite
// sides, since it's unsigned.
double m = sqrt(d.Magnitude()*dthis.Magnitude());
if(sqrt(fabs(d.Dot(dthis))) > (m - LENGTH_EPS)) {
// The edges are parallel.
if(fabs(dist_a) > LENGTH_EPS) {
// and not coincident, so can't be interesecting
return false;
}
// The edges are coincident. Make sure that neither endpoint lies
// on the other
double t;
t = a.Minus(ea).DivPivoting(d);
if(t > t_eps && t < (1 - t_eps)) return true;
t = b.Minus(ea).DivPivoting(d);
if(t > t_eps && t < (1 - t_eps)) return true;
t = ea.Minus(a).DivPivoting(dthis);
if(t > tthis_eps && t < (1 - tthis_eps)) return true;
t = eb.Minus(a).DivPivoting(dthis);
if(t > tthis_eps && t < (1 - tthis_eps)) return true;
// So coincident but disjoint, okay.
return false;
}
// Lines are not parallel, so look for an intersection.
pi = Vector::AtIntersectionOfLines(ea, eb, a, b,
&skew,
&t, &tthis);
if(skew) return false;
inOrEdge0 = (t > -t_eps) && (t < (1 + t_eps));
inOrEdge1 = (tthis > -tthis_eps) && (tthis < (1 + tthis_eps));
if(inOrEdge0 && inOrEdge1) {
if(a.Equals(ea) || b.Equals(ea) ||
a.Equals(eb) || b.Equals(eb))
{
// Not an intersection if we share an endpoint with an edge
return false;
}
// But it's an intersection if a vertex of one edge lies on the
// inside of the other (or if they cross away from either's
// vertex).
if(ppi) *ppi = pi;
if(spl) spl->Add(pi);
return true;
}
return false;
}
void SEdgeList::Clear(void) {
l.Clear();
}
@ -151,79 +220,12 @@ bool SEdgeList::AssemblePolygon(SPolygon *dest, SEdge *errorAt, bool keepDir) {
int SEdgeList::AnyEdgeCrossings(Vector a, Vector b,
Vector *ppi, SPointList *spl)
{
Vector d = b.Minus(a);
double t_eps = LENGTH_EPS/d.Magnitude();
int cnt = 0;
SEdge *se;
for(se = l.First(); se; se = l.NextAfter(se)) {
double dist_a, dist_b;
double t, tse;
bool skew;
Vector pi;
bool inOrEdge0, inOrEdge1;
Vector dse = (se->b).Minus(se->a);
double tse_eps = LENGTH_EPS/dse.Magnitude();
if(a.Equals(se->a) && b.Equals(se->b)) goto intersects;
if(b.Equals(se->a) && a.Equals(se->b)) goto intersects;
dist_a = (se->a).DistanceToLine(a, d),
dist_b = (se->b).DistanceToLine(a, d);
// Can't just test if dist_a equals dist_b; they could be on opposite
// sides, since it's unsigned.
double m = sqrt(d.Magnitude()*dse.Magnitude());
if(sqrt(fabs(d.Dot(dse))) > (m - LENGTH_EPS)) {
// The edges are parallel.
if(fabs(dist_a) > LENGTH_EPS) {
// and not coincident, so can't be interesecting
continue;
}
// The edges are coincident. Make sure that neither endpoint lies
// on the other
double t;
t = ((se->a).Minus(a)).DivPivoting(d);
if(t > t_eps && t < (1 - t_eps)) goto intersects;
t = ((se->b).Minus(a)).DivPivoting(d);
if(t > t_eps && t < (1 - t_eps)) goto intersects;
t = a.Minus(se->a).DivPivoting(dse);
if(t > tse_eps && t < (1 - tse_eps)) goto intersects;
t = b.Minus(se->a).DivPivoting(dse);
if(t > tse_eps && t < (1 - tse_eps)) goto intersects;
// So coincident but disjoint, okay.
continue;
}
// Lines are not parallel, so look for an intersection.
pi = Vector::AtIntersectionOfLines(a, b, se->a, se->b,
&skew,
&t, &tse);
if(skew) continue;
inOrEdge0 = (t > -t_eps) && (t < (1 + t_eps));
inOrEdge1 = (tse > -tse_eps) && (tse < (1 + tse_eps));
if(inOrEdge0 && inOrEdge1) {
if((se->a).Equals(a) || (se->b).Equals(a) ||
(se->a).Equals(b) || (se->b).Equals(b))
{
// Not an intersection if we share an endpoint with an edge
continue;
}
// But it's an intersection if a vertex of one edge lies on the
// inside of the other (or if they cross away from either's
// vertex).
if(ppi) *ppi = pi;
if(spl) spl->Add(pi);
goto intersects;
}
continue;
intersects:
if(se->EdgeCrosses(a, b, ppi, spl)) {
cnt++;
// and continue with the loop
}
}
return cnt;
}
@ -273,6 +275,134 @@ void SEdgeList::CullExtraneousEdges(void) {
l.RemoveTagged();
}
//-----------------------------------------------------------------------------
// Make a kd-tree of edges. This is used for O(log(n)) implementations of stuff
// that would naively be O(n).
//-----------------------------------------------------------------------------
SKdNodeEdges *SKdNodeEdges::Alloc(void) {
return (SKdNodeEdges *)AllocTemporary(sizeof(SKdNodeEdges));
}
SEdgeLl *SEdgeLl::Alloc(void) {
return (SEdgeLl *)AllocTemporary(sizeof(SEdgeLl));
}
SKdNodeEdges *SKdNodeEdges::From(SEdgeList *sel) {
SEdgeLl *sell = NULL;
SEdge *se;
for(se = sel->l.First(); se; se = sel->l.NextAfter(se)) {
SEdgeLl *n = SEdgeLl::Alloc();
n->se = se;
n->next = sell;
sell = n;
}
return SKdNodeEdges::From(sell);
}
SKdNodeEdges *SKdNodeEdges::From(SEdgeLl *sell) {
SKdNodeEdges *n = SKdNodeEdges::Alloc();
// Compute the midpoints (just mean, though median would be better) of
// each component.
Vector ptAve = Vector::From(0, 0, 0);
SEdgeLl *flip;
int totaln = 0;
for(flip = sell; flip; flip = flip->next) {
ptAve = ptAve.Plus(flip->se->a);
ptAve = ptAve.Plus(flip->se->b);
totaln++;
}
ptAve = ptAve.ScaledBy(1.0 / (2*totaln));
// For each component, see how it splits.
int ltln[3] = { 0, 0, 0 }, gtln[3] = { 0, 0, 0 };
double badness[3];
for(flip = sell; flip; flip = flip->next) {
for(int i = 0; i < 3; i++) {
if(flip->se->a.Element(i) < ptAve.Element(i) + KDTREE_EPS ||
flip->se->b.Element(i) < ptAve.Element(i) + KDTREE_EPS)
{
ltln[i]++;
}
if(flip->se->a.Element(i) > ptAve.Element(i) - KDTREE_EPS ||
flip->se->b.Element(i) > ptAve.Element(i) - KDTREE_EPS)
{
gtln[i]++;
}
}
}
for(int i = 0; i < 3; i++) {
badness[i] = pow((double)ltln[i], 4) + pow((double)gtln[i], 4);
}
// Choose the least bad coordinate to split along.
if(badness[0] < badness[1] && badness[0] < badness[2]) {
n->which = 0;
} else if(badness[1] < badness[2]) {
n->which = 1;
} else {
n->which = 2;
}
n->c = ptAve.Element(n->which);
if(totaln < 3 || totaln == gtln[n->which] || totaln == ltln[n->which]) {
n->edges = sell;
// and we're a leaf node
return n;
}
// Sort the edges according to which side(s) of the split plane they're on.
SEdgeLl *gtl = NULL, *ltl = NULL;
for(flip = sell; flip; flip = flip->next) {
if(flip->se->a.Element(n->which) < n->c + KDTREE_EPS ||
flip->se->b.Element(n->which) < n->c + KDTREE_EPS)
{
SEdgeLl *selln = SEdgeLl::Alloc();
selln->se = flip->se;
selln->next = ltl;
ltl = selln;
}
if(flip->se->a.Element(n->which) > n->c - KDTREE_EPS ||
flip->se->b.Element(n->which) > n->c - KDTREE_EPS)
{
SEdgeLl *selln = SEdgeLl::Alloc();
selln->se = flip->se;
selln->next = gtl;
gtl = selln;
}
}
n->lt = SKdNodeEdges::From(ltl);
n->gt = SKdNodeEdges::From(gtl);
return n;
}
int SKdNodeEdges::AnyEdgeCrossings(Vector a, Vector b, int cnt,
Vector *pi, SPointList *spl)
{
int inters = 0;
if(gt && lt) {
if(a.Element(which) < c + KDTREE_EPS ||
b.Element(which) < c + KDTREE_EPS)
{
inters += lt->AnyEdgeCrossings(a, b, cnt, pi, spl);
}
if(a.Element(which) > c - KDTREE_EPS ||
b.Element(which) > c - KDTREE_EPS)
{
inters += gt->AnyEdgeCrossings(a, b, cnt, pi, spl);
}
} else {
SEdgeLl *sell;
for(sell = edges; sell; sell = sell->next) {
SEdge *se = sell->se;
if(se->tag == cnt) continue;
if(se->EdgeCrosses(a, b, pi, spl)) {
inters++;
}
se->tag = cnt;
}
}
return inters;
}
//-----------------------------------------------------------------------------
// We have an edge list that contains only collinear edges, maybe with more
// splits than necessary. Merge any collinear segments that join.
@ -541,16 +671,22 @@ bool SPolygon::SelfIntersecting(Vector *intersectsAt) {
SEdgeList el;
ZERO(&el);
MakeEdgesInto(&el);
SKdNodeEdges *kdtree = SKdNodeEdges::From(&el);
int cnt = 1;
el.l.ClearTags();
bool ret = false;
SEdge *se;
for(se = el.l.First(); se; se = el.l.NextAfter(se)) {
int inters = el.AnyEdgeCrossings(se->a, se->b, intersectsAt);
int inters = kdtree->AnyEdgeCrossings(se->a, se->b, cnt, intersectsAt);
if(inters != 1) {
ret = true;
break;
}
cnt++;
}
el.Clear();
return ret;
}

38
polygon.h
View File

@ -15,6 +15,7 @@ public:
Vector a, b;
static SEdge From(Vector a, Vector b);
bool EdgeCrosses(Vector a, Vector b, Vector *pi=NULL, SPointList *spl=NULL);
};
class SEdgeList {
@ -34,6 +35,35 @@ public:
void MergeCollinearSegments(Vector a, Vector b);
};
// A kd-tree element needs to go on a side of a node if it's when KDTREE_EPS
// of the boundary. So increasing this number never breaks anything, but may
// result in more duplicated elements. So it's conservative to be sloppy here.
#define KDTREE_EPS (20*LENGTH_EPS)
class SEdgeLl {
public:
SEdge *se;
SEdgeLl *next;
static SEdgeLl *Alloc(void);
};
class SKdNodeEdges {
public:
int which; // whether c is x, y, or z
double c;
SKdNodeEdges *gt;
SKdNodeEdges *lt;
SEdgeLl *edges;
static SKdNodeEdges *From(SEdgeList *sel);
static SKdNodeEdges *From(SEdgeLl *sell);
static SKdNodeEdges *Alloc(void);
int AnyEdgeCrossings(Vector a, Vector b, int cnt,
Vector *pi=NULL, SPointList *spl=NULL);
};
class SPoint {
public:
int tag;
@ -219,7 +249,6 @@ public:
// A linked list of triangles
class STriangleLl {
public:
int tag;
STriangle *tri;
STriangleLl *next;
@ -229,10 +258,7 @@ public:
class SKdNode {
public:
static const int BY_X = 0;
static const int BY_Y = 1;
static const int BY_Z = 2;
int which;
int which; // whether c is x, y, or z
double c;
SKdNode *gt;
@ -242,7 +268,7 @@ public:
static SKdNode *Alloc(void);
static SKdNode *From(SMesh *m);
static SKdNode *From(STriangleLl *tll, int which);
static SKdNode *From(STriangleLl *tll);
void AddTriangle(STriangle *tr);
void MakeMeshInto(SMesh *m);

4
wishlist.txt
View File

@ -1,9 +1,9 @@
copy and paste
associative entities from solid model, as a special group
n*log(n) intersection finding
spline splitting
-----
associative entities from solid model, as a special group
some kind of import
faster triangulation
loop detection