caiyu/solvespace
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
27ac3c7b5f
solvespace/src/bsp.cpp
EvilSpirit 7681f6df74 Refactor insertion into BSP tree. 
Before this commit, inserting into BSP tree could easily overflow
the stack because we allocate very large stack frames and, on
convex geometries (e.g. a sphere), the BSP tree degenerates into
a "BSP list", thus requiring one large stack frame per triangle.
This can be reproduced by exporting a 2d shaded view of sphere.

After this commit, the stack frames only contan a pointer to
a supplementary data structure, and moreover it only allocates
its fields on demand, conserving heap memory as well.

As a side effect, an arbitrary classifier limit of 50 vertices
is removed.
2016-11-27 21:04:31 +00:00

726 lines
22 KiB
C++
Raw Blame History

//-----------------------------------------------------------------------------
// Binary space partitioning tree, used to represent a volume in 3-space
// bounded by a triangle mesh. These are used to compute Boolean operations
// on meshes. These aren't used for anything relating to an SShell of
// ratpoly surfaces.
//
// Copyright 2008-2013 Jonathan Westhues.
//-----------------------------------------------------------------------------
#include "solvespace.h"
SBsp2 *SBsp2::Alloc() { return (SBsp2 *)AllocTemporary(sizeof(SBsp2)); }
SBsp3 *SBsp3::Alloc() { return (SBsp3 *)AllocTemporary(sizeof(SBsp3)); }
SBsp3 *SBsp3::FromMesh(const SMesh *m) {
SBsp3 *bsp3 = NULL;
int i;
SMesh mc = {};
for(i = 0; i < m->l.n; i++) {
mc.AddTriangle(&(m->l.elem[i]));
}
srand(0); // Let's be deterministic, at least!
int n = mc.l.n;
while(n > 1) {
int k = rand() % n;
n--;
swap(mc.l.elem[k], mc.l.elem[n]);
}
for(i = 0; i < mc.l.n; i++) {
bsp3 = InsertOrCreate(bsp3, &(mc.l.elem[i]), NULL);
}
mc.Clear();
return bsp3;
}
Vector SBsp3::IntersectionWith(Vector a, Vector b) const {
double da = a.Dot(n) - d;
double db = b.Dot(n) - d;
ssassert(da*db < 0, "Expected segment to intersect BSP node");
double dab = (db - da);
return (a.ScaledBy(db/dab)).Plus(b.ScaledBy(-da/dab));
}
void SBsp3::InsertInPlane(bool pos2, STriangle *tr, SMesh *m) {
Vector tc = ((tr->a).Plus(tr->b).Plus(tr->c)).ScaledBy(1.0/3);
bool onFace = false;
bool sameNormal = false;
double maxNormalMag = -1;
Vector lln, trn = tr->Normal();
SBsp3 *ll = this;
while(ll) {
if((ll->tri).ContainsPoint(tc)) {
onFace = true;
// If the mesh contains almost-zero-area triangles, and we're
// just on the edge of one of those, then don't trust its normal.
lln = (ll->tri).Normal();
if(lln.Magnitude() > maxNormalMag) {
sameNormal = trn.Dot(lln) > 0;
maxNormalMag = lln.Magnitude();
}
}
ll = ll->more;
}
if(m->flipNormal && ((!pos2 && !onFace) ||
(onFace && !sameNormal && m->keepCoplanar)))
{
m->AddTriangle(tr->meta, tr->c, tr->b, tr->a);
} else if(!(m->flipNormal) && ((pos2 && !onFace) ||
(onFace && sameNormal && m->keepCoplanar)))
{
m->AddTriangle(tr->meta, tr->a, tr->b, tr->c);
} else {
m->atLeastOneDiscarded = true;
}
}
void SBsp3::InsertHow(BspClass how, STriangle *tr, SMesh *instead) {
switch(how) {
case BspClass::POS:
if(instead && !pos) goto alt;
pos = InsertOrCreate(pos, tr, instead);
break;
case BspClass::NEG:
if(instead && !neg) goto alt;
neg = InsertOrCreate(neg, tr, instead);
break;
case BspClass::COPLANAR: {
if(instead) goto alt;
SBsp3 *m = Alloc();
m->n = n;
m->d = d;
m->tri = *tr;
m->more = more;
more = m;
break;
}
}
return;
alt:
if(how == BspClass::POS && !(instead->flipNormal)) {
instead->AddTriangle(tr->meta, tr->a, tr->b, tr->c);
} else if(how == BspClass::NEG && instead->flipNormal) {
instead->AddTriangle(tr->meta, tr->c, tr->b, tr->a);
} else if(how == BspClass::COPLANAR) {
if(edges) {
edges->InsertTriangle(tr, instead, this);
} else {
// I suppose this actually is allowed to happen, if the coplanar
// face is the leaf, and all of its neighbors are earlier in tree?
InsertInPlane(/*pos2=*/false, tr, instead);
}
} else {
instead->atLeastOneDiscarded = true;
}
}
class BspUtil {
public:
SBsp3 *bsp;
size_t onc;
size_t posc;
size_t negc;
bool *isPos;
bool *isNeg;
bool *isOn;
// triangle operations
STriangle *tr;
STriangle *btri; // also as alone
STriangle *ctri;
// convex operations
Vector *on;
size_t npos;
size_t nneg;
Vector *vpos; // also as quad
Vector *vneg;
static BspUtil *Alloc() {
return (BspUtil *)AllocTemporary(sizeof(BspUtil));
}
void AllocOn() {
on = (Vector *)AllocTemporary(sizeof(Vector) * 2);
}
void AllocTriangle() {
btri = (STriangle *)AllocTemporary(sizeof(STriangle));
}
void AllocTriangles() {
btri = (STriangle *)AllocTemporary(sizeof(STriangle) * 2);
ctri = &btri[1];
}
void AllocQuad() {
vpos = (Vector *)AllocTemporary(sizeof(Vector) * 4);
}
void AllocClassify(size_t size) {
// Allocate a one big piece is faster than a small ones.
isPos = (bool *)AllocTemporary(sizeof(bool) * size * 3);
isNeg = &isPos[size];
isOn = &isNeg[size];
}
void AllocVertices(size_t size) {
vpos = (Vector *)AllocTemporary(sizeof(Vector) * size * 2);
vneg = &vpos[size];
}
void ClassifyTriangle(STriangle *tri, SBsp3 *node) {
tr = tri;
bsp = node;
onc = 0;
posc = 0;
negc = 0;
AllocClassify(3);
double dt[3] = { (tr->a).Dot(bsp->n), (tr->b).Dot(bsp->n), (tr->c).Dot(bsp->n) };
double d = bsp->d;
// Count vertices in the plane
for(int i = 0; i < 3; i++) {
if(dt[i] > d + LENGTH_EPS) {
posc++;
isPos[i] = true;
} else if(dt[i] < d - LENGTH_EPS) {
negc++;
isNeg[i] = true;
} else {
onc++;
isOn[i] = true;
}
}
}
bool ClassifyConvex(Vector *vertex, size_t cnt, SBsp3 *node, bool insertEdge) {
bsp = node;
onc = 0;
posc = 0;
negc = 0;
AllocClassify(cnt);
AllocOn();
for(size_t i = 0; i < cnt; i++) {
double dt = bsp->n.Dot(vertex[i]);
isPos[i] = isNeg[i] = isOn[i] = false;
if(fabs(dt - bsp->d) < LENGTH_EPS) {
isOn[i] = true;
if(onc < 2) {
on[onc] = vertex[i];
}
onc++;
} else if(dt > bsp->d) {
isPos[i] = true;
posc++;
} else {
isNeg[i] = true;
negc++;
}
}
if(onc != 2 && onc != 1 && onc != 0) return false;
if(onc == 2) {
if(insertEdge) {
Vector e01 = (vertex[1]).Minus(vertex[0]);
Vector e12 = (vertex[2]).Minus(vertex[1]);
Vector out = e01.Cross(e12);
SEdge se = SEdge::From(on[0], on[1]);
bsp->edges = SBsp2::InsertOrCreateEdge(bsp->edges, &se, bsp->n, out);
}
}
return true;
}
bool ClassifyConvexVertices(Vector *vertex, size_t cnt, bool insertEdges) {
Vector inter[2];
int inters = 0;
npos = 0;
nneg = 0;
// Enlarge vertices list to consider two intersections
AllocVertices(cnt + 4);
for(size_t i = 0; i < cnt; i++) {
size_t ip = WRAP((i + 1), cnt);
if(isPos[i]) {
vpos[npos++] = vertex[i];
}
if(isNeg[i]) {
vneg[nneg++] = vertex[i];
}
if(isOn[i]) {
vneg[nneg++] = vertex[i];
vpos[npos++] = vertex[i];
}
if((isPos[i] && isNeg[ip]) || (isNeg[i] && isPos[ip])) {
Vector vi = bsp->IntersectionWith(vertex[i], vertex[ip]);
vpos[npos++] = vi;
vneg[nneg++] = vi;
if(inters >= 2) return false; // triangulate: XXX shouldn't happen but does
inter[inters++] = vi;
}
}
ssassert(npos <= cnt + 1 && nneg <= cnt + 1, "Impossible");
if(insertEdges) {
Vector e01 = (vertex[1]).Minus(vertex[0]);
Vector e12 = (vertex[2]).Minus(vertex[1]);
Vector out = e01.Cross(e12);
if(inters == 2) {
SEdge se = SEdge::From(inter[0], inter[1]);
bsp->edges = SBsp2::InsertOrCreateEdge(bsp->edges, &se, bsp->n, out);
} else if(inters == 1 && onc == 1) {
SEdge se = SEdge::From(inter[0], on[0]);
bsp->edges = SBsp2::InsertOrCreateEdge(bsp->edges, &se, bsp->n, out);
} else if(inters == 0 && onc == 2) {
// We already handled this on-plane existing edge
} else {
return false; //triangulate;
}
}
if(nneg < 3 || npos < 3) return false; // triangulate; // XXX
return true;
}
void ProcessEdgeInsert() {
ssassert(onc == 2, "Impossible");
Vector a, b;
if (!isOn[0]) { a = tr->b; b = tr->c; }
else if(!isOn[1]) { a = tr->c; b = tr->a; }
else if(!isOn[2]) { a = tr->a; b = tr->b; }
else ssassert(false, "Impossible");
SEdge se = SEdge::From(a, b);
bsp->edges = SBsp2::InsertOrCreateEdge(bsp->edges, &se, bsp->n, tr->Normal());
}
bool SplitIntoTwoTriangles(bool insertEdge) {
ssassert(posc == 1 && negc == 1 && onc == 1, "Impossible");
bool bpos;
Vector a, b, c;
// Standardize so that a is on the plane
if (isOn[0]) { a = tr->a; b = tr->b; c = tr->c; bpos = isPos[1];
} else if(isOn[1]) { a = tr->b; b = tr->c; c = tr->a; bpos = isPos[2];
} else if(isOn[2]) { a = tr->c; b = tr->a; c = tr->b; bpos = isPos[0];
} else ssassert(false, "Impossible");
AllocTriangles();
Vector bPc = bsp->IntersectionWith(b, c);
*btri = STriangle::From(tr->meta, a, b, bPc);
*ctri = STriangle::From(tr->meta, c, a, bPc);
if(insertEdge) {
SEdge se = SEdge::From(a, bPc);
bsp->edges = SBsp2::InsertOrCreateEdge(bsp->edges, &se, bsp->n, tr->Normal());
}
return bpos;
}
bool SplitIntoTwoPieces(bool insertEdge) {
Vector a, b, c;
if(posc == 2 && negc == 1) {
// Standardize so that a is on one side, and b and c are on the other.
if (isNeg[0]) { a = tr->a; b = tr->b; c = tr->c;
} else if(isNeg[1]) { a = tr->b; b = tr->c; c = tr->a;
} else if(isNeg[2]) { a = tr->c; b = tr->a; c = tr->b;
} else ssassert(false, "Impossible");
} else if(posc == 1 && negc == 2) {
if (isPos[0]) { a = tr->a; b = tr->b; c = tr->c;
} else if(isPos[1]) { a = tr->b; b = tr->c; c = tr->a;
} else if(isPos[2]) { a = tr->c; b = tr->a; c = tr->b;
} else ssassert(false, "Impossible");
} else ssassert(false, "Impossible");
Vector aPb = bsp->IntersectionWith(a, b);
Vector cPa = bsp->IntersectionWith(c, a);
AllocTriangle();
AllocQuad();
*btri = STriangle::From(tr->meta, a, aPb, cPa);
vpos[0] = aPb;
vpos[1] = b;
vpos[2] = c;
vpos[3] = cPa;
if(insertEdge) {
SEdge se = SEdge::From(aPb, cPa);
bsp->edges = SBsp2::InsertOrCreateEdge(bsp->edges, &se, bsp->n, btri->Normal());
}
return posc == 2 && negc == 1;
}
static SBsp3 *Triangulate(SBsp3 *bsp, const STriMeta &meta, Vector *vertex,
size_t cnt, SMesh *instead) {
for(size_t i = 0; i < cnt - 2; i++) {
STriangle tr = STriangle::From(meta, vertex[0], vertex[i + 1], vertex[i + 2]);
bsp = SBsp3::InsertOrCreate(bsp, &tr, instead);
}
return bsp;
}
};
void SBsp3::InsertConvexHow(BspClass how, STriMeta meta, Vector *vertex, size_t n,
SMesh *instead) {
switch(how) {
case BspClass::POS:
if(pos) {
pos = pos->InsertConvex(meta, vertex, n, instead);
return;
}
break;
case BspClass::NEG:
if(neg) {
neg = neg->InsertConvex(meta, vertex, n, instead);
return;
}
break;
default: ssassert(false, "Unexpected BSP insert type");
}
for(size_t i = 0; i < n - 2; i++) {
STriangle tr = STriangle::From(meta,
vertex[0], vertex[i+1], vertex[i+2]);
InsertHow(how, &tr, instead);
}
}
SBsp3 *SBsp3::InsertConvex(STriMeta meta, Vector *vertex, size_t cnt, SMesh *instead) {
BspUtil *u = BspUtil::Alloc();
if(u->ClassifyConvex(vertex, cnt, this, !instead)) {
if(u->posc == 0) {
InsertConvexHow(BspClass::NEG, meta, vertex, cnt, instead);
return this;
}
if(u->negc == 0) {
InsertConvexHow(BspClass::POS, meta, vertex, cnt, instead);
return this;
}
if(u->ClassifyConvexVertices(vertex, cnt, !instead)) {
InsertConvexHow(BspClass::NEG, meta, u->vneg, u->nneg, instead);
InsertConvexHow(BspClass::POS, meta, u->vpos, u->npos, instead);
return this;
}
}
// We don't handle the special case for this; do it as triangles
return BspUtil::Triangulate(this, meta, vertex, cnt, instead);
}
SBsp3 *SBsp3::InsertOrCreate(SBsp3 *where, STriangle *tr, SMesh *instead) {
if(where == NULL) {
if(instead) {
if(instead->flipNormal) {
instead->atLeastOneDiscarded = true;
} else {
instead->AddTriangle(tr->meta, tr->a, tr->b, tr->c);
}
return NULL;
}
// Brand new node; so allocate for it, and fill us in.
SBsp3 *r = Alloc();
r->n = (tr->Normal()).WithMagnitude(1);
r->d = (tr->a).Dot(r->n);
r->tri = *tr;
return r;
}
where->Insert(tr, instead);
return where;
}
void SBsp3::Insert(STriangle *tr, SMesh *instead) {
BspUtil *u = BspUtil::Alloc();
u->ClassifyTriangle(tr, this);
// All vertices in-plane
if(u->onc == 3) {
InsertHow(BspClass::COPLANAR, tr, instead);
return;
}
// No split required
if(u->posc == 0 || u->negc == 0) {
if(!instead && u->onc == 2) {
u->ProcessEdgeInsert();
}
if(u->posc > 0) {
InsertHow(BspClass::POS, tr, instead);
} else {
InsertHow(BspClass::NEG, tr, instead);
}
return;
}
// The polygon must be split into two triangles, one above, one below.
if(u->posc == 1 && u->negc == 1 && u->onc == 1) {
if(u->SplitIntoTwoTriangles(!instead)) {
InsertHow(BspClass::POS, u->btri, instead);
InsertHow(BspClass::NEG, u->ctri, instead);
} else {
InsertHow(BspClass::POS, u->ctri, instead);
InsertHow(BspClass::NEG, u->btri, instead);
}
return;
}
// The polygon must be split into two pieces: a triangle and a quad.
if(u->SplitIntoTwoPieces(!instead)) {
InsertConvexHow(BspClass::POS, tr->meta, u->vpos, 4, instead);
InsertHow(BspClass::NEG, u->btri, instead);
} else {
InsertConvexHow(BspClass::NEG, tr->meta, u->vpos, 4, instead);
InsertHow(BspClass::POS, u->btri, instead);
}
}
void SBsp3::GenerateInPaintOrder(SMesh *m) const {
// Doesn't matter which branch we take if the normal has zero z
// component, so don't need a separate case for that.
if(n.z < 0) {
if(pos) pos->GenerateInPaintOrder(m);
} else {
if(neg) neg->GenerateInPaintOrder(m);
}
const SBsp3 *flip = this;
while(flip) {
m->AddTriangle(&(flip->tri));
flip = flip->more;
}
if(n.z < 0) {
if(neg) neg->GenerateInPaintOrder(m);
} else {
if(pos) pos->GenerateInPaintOrder(m);
}
}
/////////////////////////////////
Vector SBsp2::IntersectionWith(Vector a, Vector b) const {
double da = a.Dot(no) - d;
double db = b.Dot(no) - d;
ssassert(da*db < 0, "Expected segment to intersect BSP node");
double dab = (db - da);
return (a.ScaledBy(db/dab)).Plus(b.ScaledBy(-da/dab));
}
SBsp2 *SBsp2::InsertOrCreateEdge(SBsp2 *where, SEdge *nedge, Vector nnp, Vector out) {
if(where == NULL) {
// Brand new node; so allocate for it, and fill us in.
SBsp2 *r = Alloc();
r->np = nnp;
r->no = ((r->np).Cross((nedge->b).Minus(nedge->a))).WithMagnitude(1);
if(out.Dot(r->no) < 0) {
r->no = (r->no).ScaledBy(-1);
}
r->d = (nedge->a).Dot(r->no);
r->edge = *nedge;
return r;
}
where->InsertEdge(nedge, nnp, out);
return where;
}
void SBsp2::InsertEdge(SEdge *nedge, Vector nnp, Vector out) {
double dt[2] = { (nedge->a).Dot(no), (nedge->b).Dot(no) };
bool isPos[2] = {}, isNeg[2] = {}, isOn[2] = {};
for(int i = 0; i < 2; i++) {
if(fabs(dt[i] - d) < LENGTH_EPS) {
isOn[i] = true;
} else if(dt[i] > d) {
isPos[i] = true;
} else {
isNeg[i] = true;
}
}
if((isPos[0] && isPos[1])||(isPos[0] && isOn[1])||(isOn[0] && isPos[1])) {
pos = InsertOrCreateEdge(pos, nedge, nnp, out);
return;
}
if((isNeg[0] && isNeg[1])||(isNeg[0] && isOn[1])||(isOn[0] && isNeg[1])) {
neg = InsertOrCreateEdge(neg, nedge, nnp, out);
return;
}
if(isOn[0] && isOn[1]) {
SBsp2 *m = Alloc();
m->np = nnp;
m->no = ((m->np).Cross((nedge->b).Minus(nedge->a))).WithMagnitude(1);
if(out.Dot(m->no) < 0) {
m->no = (m->no).ScaledBy(-1);
}
m->d = (nedge->a).Dot(m->no);
m->edge = *nedge;
m->more = more;
more = m;
return;
}
if((isPos[0] && isNeg[1]) || (isNeg[0] && isPos[1])) {
Vector aPb = IntersectionWith(nedge->a, nedge->b);
SEdge ea = SEdge::From(nedge->a, aPb);
SEdge eb = SEdge::From(aPb, nedge->b);
if(isPos[0]) {
pos = InsertOrCreateEdge(pos, &ea, nnp, out);
neg = InsertOrCreateEdge(neg, &eb, nnp, out);
} else {
neg = InsertOrCreateEdge(neg, &ea, nnp, out);
pos = InsertOrCreateEdge(pos, &eb, nnp, out);
}
return;
}
ssassert(false, "Impossible");
}
void SBsp2::InsertTriangleHow(BspClass how, STriangle *tr, SMesh *m, SBsp3 *bsp3) {
switch(how) {
case BspClass::POS:
if(pos) {
pos->InsertTriangle(tr, m, bsp3);
} else {
bsp3->InsertInPlane(/*pos2=*/true, tr, m);
}
break;
case BspClass::NEG:
if(neg) {
neg->InsertTriangle(tr, m, bsp3);
} else {
bsp3->InsertInPlane(/*pos2=*/false, tr, m);
}
break;
default: ssassert(false, "Unexpected BSP insert type");
}
}
void SBsp2::InsertTriangle(STriangle *tr, SMesh *m, SBsp3 *bsp3) {
double dt[3] = { (tr->a).Dot(no), (tr->b).Dot(no), (tr->c).Dot(no) };
bool isPos[3] = {}, isNeg[3] = {}, isOn[3] = {};
int inc = 0, posc = 0, negc = 0;
for(int i = 0; i < 3; i++) {
if(fabs(dt[i] - d) < LENGTH_EPS) {
isOn[i] = true;
inc++;
} else if(dt[i] > d) {
isPos[i] = true;
posc++;
} else {
isNeg[i] = true;
negc++;
}
}
if(inc == 3) {
// All vertices on-line; so it's a degenerate triangle, to ignore.
return;
}
// No split required
if(posc == 0 || negc == 0) {
if(posc > 0) {
InsertTriangleHow(BspClass::POS, tr, m, bsp3);
} else {
InsertTriangleHow(BspClass::NEG, tr, m, bsp3);
}
return;
}
// The polygon must be split into two pieces, one above, one below.
Vector a, b, c;
if(posc == 1 && negc == 1 && inc == 1) {
bool bpos;
// Standardize so that a is on the plane
if (isOn[0]) { a = tr->a; b = tr->b; c = tr->c; bpos = isPos[1];
} else if(isOn[1]) { a = tr->b; b = tr->c; c = tr->a; bpos = isPos[2];
} else if(isOn[2]) { a = tr->c; b = tr->a; c = tr->b; bpos = isPos[0];
} else ssassert(false, "Impossible");
Vector bPc = IntersectionWith(b, c);
STriangle btri = STriangle::From(tr->meta, a, b, bPc);
STriangle ctri = STriangle::From(tr->meta, c, a, bPc);
if(bpos) {
InsertTriangleHow(BspClass::POS, &btri, m, bsp3);
InsertTriangleHow(BspClass::NEG, &ctri, m, bsp3);
} else {
InsertTriangleHow(BspClass::POS, &ctri, m, bsp3);
InsertTriangleHow(BspClass::NEG, &btri, m, bsp3);
}
return;
}
if(posc == 2 && negc == 1) {
// Standardize so that a is on one side, and b and c are on the other.
if (isNeg[0]) { a = tr->a; b = tr->b; c = tr->c;
} else if(isNeg[1]) { a = tr->b; b = tr->c; c = tr->a;
} else if(isNeg[2]) { a = tr->c; b = tr->a; c = tr->b;
} else ssassert(false, "Impossible");
} else if(posc == 1 && negc == 2) {
if (isPos[0]) { a = tr->a; b = tr->b; c = tr->c;
} else if(isPos[1]) { a = tr->b; b = tr->c; c = tr->a;
} else if(isPos[2]) { a = tr->c; b = tr->a; c = tr->b;
} else ssassert(false, "Impossible");
} else ssassert(false, "Impossible");
Vector aPb = IntersectionWith(a, b);
Vector cPa = IntersectionWith(c, a);
STriangle alone = STriangle::From(tr->meta, a, aPb, cPa);
STriangle quad1 = STriangle::From(tr->meta, aPb, b, c );
STriangle quad2 = STriangle::From(tr->meta, aPb, c, cPa);
if(posc == 2 && negc == 1) {
InsertTriangleHow(BspClass::POS, &quad1, m, bsp3);
InsertTriangleHow(BspClass::POS, &quad2, m, bsp3);
InsertTriangleHow(BspClass::NEG, &alone, m, bsp3);
} else {
InsertTriangleHow(BspClass::NEG, &quad1, m, bsp3);
InsertTriangleHow(BspClass::NEG, &quad2, m, bsp3);
InsertTriangleHow(BspClass::POS, &alone, m, bsp3);
}
return;
}
Reference in New Issue View Git Blame Copy Permalink