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));

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((!onFace && ((m->keepInsideOtherShell && !pos2) ||

(!m->keepInsideOtherShell && pos2))) ||

(onFace && ((m->keepCoplanar && m->flipNormal && !sameNormal) ||

(m->keepCoplanar && !m->flipNormal && sameNormal)))) {

// We have decided that we need to keep a triangle either inside,

// outside or on the other shell. So add it and flip it if requested.

if(!(m->flipNormal)) {

m->AddTriangle(tr->meta, tr->a, tr->b, tr->c);

} else {

m->AddTriangle(tr->meta, tr->c, tr->b, tr->a);

}

} else {

m->atLeastOneDiscarded = true;

}

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(((BspClass::POS == how) && !instead->keepInsideOtherShell) ||

((BspClass::NEG == how) && instead->keepInsideOtherShell)) {

// We have decided that we need to keep a triangle (either inside or

// outside the other shell. So add it and flip it if requested.

if(!(instead->flipNormal)) {

instead->AddTriangle(tr->meta, tr->a, tr->b, tr->c);

} else {

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;

}

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;

}

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);

}

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);

}

// 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);

}

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));

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");

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");

}

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;