Vector ab = b.Minus(a), bc = c.Minus(b);

Vector np = ab.Cross(bc);

if(np.Magnitude() < 1e-10) {

// ugh; gl sometimes tesselates to collinear triangles

return;

}

if(np.Dot(n) > 0) {

AddTriangle(meta, a, b, c);

} else {

AddTriangle(meta, c, b, a);

}

STriangle t = {};

t.meta = meta;

t.a = a;

t.b = b;

t.c = c;

AddTriangle(&t);

vmax->x = max(vmax->x, v.x);

vmax->y = max(vmax->y, v.y);

vmax->z = max(vmax->z, v.z);

vmin->x = min(vmin->x, v.x);

vmin->y = min(vmin->y, v.y);

vmin->z = min(vmin->z, v.z);

int i;

*vmin = Vector::From( 1e12, 1e12, 1e12);

*vmax = Vector::From(-1e12, -1e12, -1e12);

for(i = 0; i < l.n; i++) {

STriangle *st = &(l.elem[i]);

DoBounding(st->a, vmax, vmin);

DoBounding(st->b, vmax, vmin);

DoBounding(st->c, vmax, vmin);

}

SMesh m = {};

m.MakeFromCopyOf(this);

// Delete all triangles in the mesh that do not lie in our export plane.

m.l.ClearTags();

int i;

for(i = 0; i < m.l.n; i++) {

STriangle *tr = &(m.l.elem[i]);

if((fabs(n.Dot(tr->a) - d) >= LENGTH_EPS) ||

(fabs(n.Dot(tr->b) - d) >= LENGTH_EPS) ||

(fabs(n.Dot(tr->c) - d) >= LENGTH_EPS))

{

tr->tag = 1;

}

}

m.l.RemoveTagged();

// Select the naked edges in our resulting open mesh.

SKdNode *root = SKdNode::From(&m);

root->SnapToMesh(&m);

root->MakeCertainEdgesInto(sel, EdgeKind::NAKED_OR_SELF_INTER,

/*coplanarIsInter=*/false, NULL, NULL);

m.Clear();

int maxTriangles = (l.n - start) + 10;

STriMeta meta = l.elem[start].meta;

STriangle *tout = (STriangle *)AllocTemporary(maxTriangles*sizeof(*tout));

int toutc = 0;

Vector n = Vector::From(0, 0, 0);

Vector *conv = (Vector *)AllocTemporary(maxTriangles*3*sizeof(*conv));

int convc = 0;

int start0 = start;

int i, j;

for(i = start; i < l.n; i++) {

STriangle *tr = &(l.elem[i]);

if(tr->MinAltitude() < LENGTH_EPS) {

tr->tag = 1;

} else {

tr->tag = 0;

}

}

for(;;) {

bool didAdd;

convc = 0;

for(i = start; i < l.n; i++) {

STriangle *tr = &(l.elem[i]);

if(tr->tag) continue;

tr->tag = 1;

n = (tr->Normal()).WithMagnitude(1);

conv[convc++] = tr->a;

conv[convc++] = tr->b;

conv[convc++] = tr->c;

start = i+1;

break;

}

if(i >= l.n) break;

do {

didAdd = false;

for(j = 0; j < convc; j++) {

Vector a = conv[WRAP((j-1), convc)],

b = conv[j],

d = conv[WRAP((j+1), convc)],

e = conv[WRAP((j+2), convc)];

Vector c;

for(i = start; i < l.n; i++) {

STriangle *tr = &(l.elem[i]);

if(tr->tag) continue;

if((tr->a).Equals(d) && (tr->b).Equals(b)) {

c = tr->c;

} else if((tr->b).Equals(d) && (tr->c).Equals(b)) {

c = tr->a;

} else if((tr->c).Equals(d) && (tr->a).Equals(b)) {

c = tr->b;

} else {

continue;

}

// The vertex at C must be convex; but the others must

// be tested

Vector ab = b.Minus(a);

Vector bc = c.Minus(b);

Vector cd = d.Minus(c);

Vector de = e.Minus(d);

double bDot = (ab.Cross(bc)).Dot(n);

double dDot = (cd.Cross(de)).Dot(n);

bDot /= min(ab.Magnitude(), bc.Magnitude());

dDot /= min(cd.Magnitude(), de.Magnitude());

if(fabs(bDot) < LENGTH_EPS && fabs(dDot) < LENGTH_EPS) {

conv[WRAP((j+1), convc)] = c;

// and remove the vertex at j, which is a dup

memmove(conv+j, conv+j+1,

(convc - j - 1)*sizeof(conv[0]));

convc--;

} else if(fabs(bDot) < LENGTH_EPS && dDot > 0) {

conv[j] = c;

} else if(fabs(dDot) < LENGTH_EPS && bDot > 0) {

conv[WRAP((j+1), convc)] = c;

} else if(bDot > 0 && dDot > 0) {

// conv[j] is unchanged, conv[j+1] goes to [j+2]

memmove(conv+j+2, conv+j+1,

(convc - j - 1)*sizeof(conv[0]));

conv[j+1] = c;

convc++;

} else {

continue;

}

didAdd = true;

tr->tag = 1;

break;

}

}

} while(didAdd);

// I need to debug why this is required; sometimes the above code

// still generates a convex polygon

for(i = 0; i < convc; i++) {

Vector a = conv[WRAP((i-1), convc)],

b = conv[i],

c = conv[WRAP((i+1), convc)];

Vector ab = b.Minus(a);

Vector bc = c.Minus(b);

double bDot = (ab.Cross(bc)).Dot(n);

bDot /= min(ab.Magnitude(), bc.Magnitude());

if(bDot < 0) return; // XXX, shouldn't happen

}

for(i = 0; i < convc - 2; i++) {

STriangle tr = STriangle::From(meta, conv[0], conv[i+1], conv[i+2]);

if(tr.MinAltitude() > LENGTH_EPS) {

tout[toutc++] = tr;

}

}

}

l.n = start0;

for(i = 0; i < toutc; i++) {

AddTriangle(&(tout[i]));

}

FreeTemporary(tout);

FreeTemporary(conv);

int i;

for(i = 0; i < srcm->l.n; i++) {

STriangle *st = &(srcm->l.elem[i]);

int pn = l.n;

atLeastOneDiscarded = false;

SBsp3::InsertOrCreate(bsp3, st, this);

if(!atLeastOneDiscarded && (l.n != (pn+1))) {

l.n = pn;

if(flipNormal) {

AddTriangle(st->meta, st->c, st->b, st->a);

} else {

AddTriangle(st->meta, st->a, st->b, st->c);

}

}

if(l.n - pn > 1) {

Simplify(pn);

}

}

SBsp3 *bspa = SBsp3::FromMesh(a);

SBsp3 *bspb = SBsp3::FromMesh(b);

flipNormal = false;

keepCoplanar = false;

AddAgainstBsp(b, bspa);

flipNormal = false;

keepCoplanar = true;

AddAgainstBsp(a, bspb);

SBsp3 *bspa = SBsp3::FromMesh(a);

SBsp3 *bspb = SBsp3::FromMesh(b);

flipNormal = true;

keepCoplanar = true;

AddAgainstBsp(b, bspa);

flipNormal = false;

keepCoplanar = false;

AddAgainstBsp(a, bspb);

int i;

for(i = 0; i < a->l.n; i++) {

AddTriangle(&(a->l.elem[i]));

}

Quaternion q, double scale)

{

STriangle *tr;

for(tr = a->l.First(); tr; tr = a->l.NextAfter(tr)) {

STriangle tt = *tr;

tt.a = (tt.a).ScaledBy(scale);

tt.b = (tt.b).ScaledBy(scale);

tt.c = (tt.c).ScaledBy(scale);

if(scale < 0) {

// The mirroring would otherwise turn a closed mesh inside out.

swap(tt.a, tt.b);

}

tt.a = (q.Rotate(tt.a)).Plus(trans);

tt.b = (q.Rotate(tt.b)).Plus(trans);

tt.c = (q.Rotate(tt.c)).Plus(trans);

AddTriangle(&tt);

}

Vector rayPoint = SS.GW.UnProjectPoint3(Vector::From(mp.x, mp.y, 0.0));

Vector rayDir = SS.GW.UnProjectPoint3(Vector::From(mp.x, mp.y, 1.0)).Minus(rayPoint);

uint32_t face = 0;

double faceT = VERY_NEGATIVE;

for(int i = 0; i < l.n; i++) {

const STriangle &tr = l.elem[i];

if(tr.meta.face == 0) continue;

double t;

if(!tr.Raytrace(rayPoint, rayDir, &t, NULL)) continue;

if(t > faceT) {

face = tr.meta.face;

faceT = t;

}

}

return face;

Vector center = {};

double vol = 0.0;

for(int i = 0; i < l.n; i++) {

STriangle &tr = l.elem[i];

double tvol = tr.SignedVolume();

center = center.Plus(tr.a.Plus(tr.b.Plus(tr.c)).ScaledBy(tvol / 4.0));

vol += tvol;

}

return center.ScaledBy(1.0 / vol);

if(gt && lt) {

gt->ClearTags();

lt->ClearTags();

} else {

STriangleLl *ll;

for(ll = tris; ll; ll = ll->next) {

ll->tri->tag = 0;

}

}

if(gt && lt) {

double ta = (tr->a).Element(which),

tb = (tr->b).Element(which),

tc = (tr->c).Element(which);

if(ta < c + KDTREE_EPS ||

tb < c + KDTREE_EPS ||

tc < c + KDTREE_EPS)

{

lt->AddTriangle(tr);

}

if(ta > c - KDTREE_EPS ||

tb > c - KDTREE_EPS ||

tc > c - KDTREE_EPS)

{

gt->AddTriangle(tr);

}

} else {

STriangleLl *tn = STriangleLl::Alloc();

tn->tri = tr;

tn->next = tris;

tris = tn;

}

if(gt) gt->MakeMeshInto(m);

if(lt) lt->MakeMeshInto(m);

STriangleLl *ll;

for(ll = tris; ll; ll = ll->next) {

if(ll->tri->tag) continue;

m->AddTriangle(ll->tri);

ll->tri->tag = 1;

}

if(gt) gt->ListTrianglesInto(tl);

if(lt) lt->ListTrianglesInto(tl);

STriangleLl *ll;

for(ll = tris; ll; ll = ll->next) {

if(ll->tri->tag) continue;

tl->push_back(ll->tri);

ll->tri->tag = 1;

}

if(gt && lt) {

double vc = v.Element(which);

if(vc < c + KDTREE_EPS) {

lt->SnapToVertex(v, extras);

}

if(vc > c - KDTREE_EPS) {

gt->SnapToVertex(v, extras);

}

// Nothing bad happens if the triangle to be split appears in both

// branches; the first call will split the triangle, so that the

// second call will do nothing, because the modified triangle will

// already contain v

} else {

STriangleLl *ll;

for(ll = tris; ll; ll = ll->next) {

STriangle *tr = ll->tri;

// Do a cheap bbox test first

int k;

bool mightHit = true;

for(k = 0; k < 3; k++) {

if((tr->a).Element(k) < v.Element(k) - KDTREE_EPS &&

(tr->b).Element(k) < v.Element(k) - KDTREE_EPS &&

(tr->c).Element(k) < v.Element(k) - KDTREE_EPS)

{

mightHit = false;

break;

}

if((tr->a).Element(k) > v.Element(k) + KDTREE_EPS &&

(tr->b).Element(k) > v.Element(k) + KDTREE_EPS &&

(tr->c).Element(k) > v.Element(k) + KDTREE_EPS)

{

mightHit = false;

break;

}

}

if(!mightHit) continue;

if(tr->a.Equals(v)) { tr->a = v; continue; }

if(tr->b.Equals(v)) { tr->b = v; continue; }

if(tr->c.Equals(v)) { tr->c = v; continue; }

if(v.OnLineSegment(tr->a, tr->b)) {

STriangle nt = STriangle::From(tr->meta, tr->a, v, tr->c);

extras->AddTriangle(&nt);

tr->a = v;

continue;

}

if(v.OnLineSegment(tr->b, tr->c)) {

STriangle nt = STriangle::From(tr->meta, tr->b, v, tr->a);

extras->AddTriangle(&nt);

tr->b = v;

continue;

}

if(v.OnLineSegment(tr->c, tr->a)) {

STriangle nt = STriangle::From(tr->meta, tr->c, v, tr->b);

extras->AddTriangle(&nt);

tr->c = v;

continue;

}

}

}

int i, j, k;

for(i = 0; i < m->l.n; i++) {

STriangle *tr = &(m->l.elem[i]);

for(j = 0; j < 3; j++) {

Vector v = tr->vertices[j];

SMesh extra = {};

SnapToVertex(v, &extra);

for(k = 0; k < extra.l.n; k++) {

STriangle *tra = (STriangle *)AllocTemporary(sizeof(*tra));

*tra = extra.l.elem[k];

AddTriangle(tra);

}

extra.Clear();

}

}

SEdgeList seln = {};

Vector tn = tr->Normal().WithMagnitude(1);

double td = tn.Dot(tr->a);

// Consider front-facing triangles only.

if(tn.z > LENGTH_EPS) {

// If the edge crosses our triangle's plane, then split into above

// and below parts. Note that we must preserve auxA, which contains

// the style associated with this line, as well as the tag, which

// contains the occlusion status.

SEdge *se;

for(se = sel->l.First(); se; se = sel->l.NextAfter(se)) {

double da = (se->a).Dot(tn) - td,

db = (se->b).Dot(tn) - td;

if((da < -LENGTH_EPS && db > LENGTH_EPS) ||

(db < -LENGTH_EPS && da > LENGTH_EPS))

{

Vector m = Vector::AtIntersectionOfPlaneAndLine(

tn, td,

se->a, se->b, NULL);

seln.AddEdge(m, se->b, se->auxA, 0, se->tag);

se->b = m;

}

}

for(se = seln.l.First(); se; se = seln.l.NextAfter(se)) {

sel->AddEdge(se->a, se->b, se->auxA, 0, se->tag);

}

seln.Clear();

for(se = sel->l.First(); se; se = sel->l.NextAfter(se)) {

Vector pt = ((se->a).Plus(se->b)).ScaledBy(0.5);

if(pt.Dot(tn) - td > -LENGTH_EPS) {

// Edge is in front of or on our plane (remember, tn.z > 0)

// so it is exempt from further splitting

se->auxB = 1;

} else {

// Edge is behind our plane, needs further splitting

se->auxB = 0;

}

}

// Considering only the (x, y) coordinates, split the edge against our

// triangle.

Point2d a = (tr->a).ProjectXy(),

b = (tr->b).ProjectXy(),

c = (tr->c).ProjectXy();

Point2d n[3] = { (b.Minus(a)).Normal().WithMagnitude(1),

(c.Minus(b)).Normal().WithMagnitude(1),

(a.Minus(c)).Normal().WithMagnitude(1) };

double d[3] = { n[0].Dot(b),

n[1].Dot(c),

n[2].Dot(a) };

// Split all of the edges where they intersect the triangle edges

int i;

for(i = 0; i < 3; i++) {

for(se = sel->l.First(); se; se = sel->l.NextAfter(se)) {

if(se->auxB) continue;

Point2d ap = (se->a).ProjectXy(),

bp = (se->b).ProjectXy();

double da = n[i].Dot(ap) - d[i],

db = n[i].Dot(bp) - d[i];

if((da < -LENGTH_EPS && db > LENGTH_EPS) ||

(db < -LENGTH_EPS && da > LENGTH_EPS))

{

double dab = (db - da);

Vector spl = ((se->a).ScaledBy( db/dab)).Plus(

(se->b).ScaledBy(-da/dab));

seln.AddEdge(spl, se->b, se->auxA, 0, se->tag);

se->b = spl;

}

}

for(se = seln.l.First(); se; se = seln.l.NextAfter(se)) {

// The split pieces are all behind the triangle, since only

// edges behind the triangle got split. So their auxB is 0.

sel->AddEdge(se->a, se->b, se->auxA, 0, se->tag);

}

seln.Clear();

}

for(se = sel->l.First(); se; se = sel->l.NextAfter(se)) {

bool occluded;

if(se->auxB) {

// Lies above or on the triangle plane, so triangle doesn't

// occlude it.

occluded = false;

} else {

// Test the segment to see if it lies outside the triangle

// (i.e., outside wrt at least one edge), and keep it only

// then.

Point2d pt = ((se->a).Plus(se->b).ScaledBy(0.5)).ProjectXy();

occluded = true;

for(i = 0; i < 3; i++) {

// If the test point lies on the boundary of our triangle,

// then we still discard the edge.

if(n[i].Dot(pt) - d[i] > LENGTH_EPS) occluded = false;

}

}

if(occluded) {

se->tag = 1;

}

}

}

if(gt && lt) {

double ac = (orig.a).Element(which),

bc = (orig.b).Element(which);

// We can ignore triangles that are separated in x or y, but triangles

// that are separated in z may still contribute

if(ac < c + KDTREE_EPS ||

bc < c + KDTREE_EPS ||

which == 2)

{

lt->OcclusionTestLine(orig, sel, cnt);

}

if(ac > c - KDTREE_EPS ||

bc > c - KDTREE_EPS ||

which == 2)

{

gt->OcclusionTestLine(orig, sel, cnt);

}

} else {

STriangleLl *ll;

for(ll = tris; ll; ll = ll->next) {

STriangle *tr = ll->tri;

if(tr->tag == cnt) continue;

SplitLinesAgainstTriangle(sel, tr);

tr->tag = cnt;

}

}

EdgeOnInfo *info) const

{

if(gt && lt) {

double ac = a.Element(which),

bc = b.Element(which);

if(ac < c + KDTREE_EPS ||

bc < c + KDTREE_EPS)

{

lt->FindEdgeOn(a, b, cnt, coplanarIsInter, info);

}

if(ac > c - KDTREE_EPS ||

bc > c - KDTREE_EPS)

{

gt->FindEdgeOn(a, b, cnt, coplanarIsInter, info);

}

return;

}

// We are a leaf node; so we iterate over all the triangles in our

// linked list.

STriangleLl *ll;

for(ll = tris; ll; ll = ll->next) {

STriangle *tr = ll->tri;

if(tr->tag == cnt) continue;

// Test if this triangle matches up with the given edge

if((a.Equals(tr->b) && b.Equals(tr->a)) ||

(a.Equals(tr->c) && b.Equals(tr->b)) ||

(a.Equals(tr->a) && b.Equals(tr->c)))

{

info->count++;

// Record whether this triangle is front- or back-facing.

if(tr->Normal().z > LENGTH_EPS) {

info->frontFacing = true;

} else {

info->frontFacing = false;

}

// Record the triangle

info->tr = tr;

// And record which vertexes a and b correspond to

info->ai = a.Equals(tr->a) ? 0 : (a.Equals(tr->b) ? 1 : 2);

info->bi = b.Equals(tr->a) ? 0 : (b.Equals(tr->b) ? 1 : 2);

} else if(((a.Equals(tr->a) && b.Equals(tr->b)) ||

(a.Equals(tr->b) && b.Equals(tr->c)) ||

(a.Equals(tr->c) && b.Equals(tr->a))))

{

// It's an edge of this triangle, okay.

} else {

// Check for self-intersection

Vector n = (tr->Normal()).WithMagnitude(1);

double d = (tr->a).Dot(n);

double pa = a.Dot(n) - d, pb = b.Dot(n) - d;

// It's an intersection if neither point lies in-plane,

// and the edge crosses the plane (should handle in-plane

// intersections separately but don't yet).

if((pa < -LENGTH_EPS || pa > LENGTH_EPS) &&

(pb < -LENGTH_EPS || pb > LENGTH_EPS) &&

(pa*pb < 0))

{

// The edge crosses the plane of the triangle; now see if

// it crosses inside the triangle.

if(tr->ContainsPointProjd(b.Minus(a), a)) {

if(coplanarIsInter) {

info->intersectsMesh = true;

} else {

Vector p = Vector::AtIntersectionOfPlaneAndLine(

n, d, a, b, NULL);

Vector ta = tr->a,

tb = tr->b,

tc = tr->c;

if((p.DistanceToLine(ta, tb.Minus(ta)) < LENGTH_EPS) ||

(p.DistanceToLine(tb, tc.Minus(tb)) < LENGTH_EPS) ||

(p.DistanceToLine(tc, ta.Minus(tc)) < LENGTH_EPS))

{

// Intersection lies on edge. This happens when

// our edge is from a triangle coplanar with

// another triangle in the mesh. We don't test

// the edge against triangles whose plane contains

// that edge, but we do end up testing against

// the coplanar triangle's neighbours, which we

// will intersect on their edges.

} else {

info->intersectsMesh = true;

}

}

}

}

}

// Ensure that we don't count this triangle twice if it appears

// in two buckets of the kd tree.

tr->tag = cnt;

}

STriangle *a, STriangle *b)

{

if(pairs->find(std::make_pair(a, b)) != pairs->end() ||

pairs->find(std::make_pair(b, a)) != pairs->end())

return true;

pairs->emplace(a, b);

return false;

bool *inter, bool *leaky, int auxA) const

{

if(inter) *inter = false;

if(leaky) *leaky = false;

std::vector<STriangle *> tris;

ClearTags();

ListTrianglesInto(&tris);

std::set<std::pair<STriangle *, STriangle *>> edgeTris;

int cnt = 1234;

for(STriangle *tr : tris) {

for(int j = 0; j < 3; j++) {

Vector a = tr->vertices[j];

Vector b = tr->vertices[(j + 1) % 3];

SKdNode::EdgeOnInfo info = {};

FindEdgeOn(a, b, cnt, coplanarIsInter, &info);

switch(how) {

case EdgeKind::NAKED_OR_SELF_INTER:

if(info.count != 1) {

sel->AddEdge(a, b, auxA);

if(leaky) *leaky = true;

}

if(info.intersectsMesh) {

sel->AddEdge(a, b, auxA);

if(inter) *inter = true;

}

break;

case EdgeKind::SELF_INTER:

if(info.intersectsMesh) {

sel->AddEdge(a, b, auxA);

if(inter) *inter = true;

}

break;

case EdgeKind::TURNING:

if((tr->Normal().z < LENGTH_EPS) &&

(info.count == 1) &&

info.frontFacing)

{

if(CheckAndAddTrianglePair(&edgeTris, tr, info.tr))

break;

// This triangle is back-facing (or on edge), and

// this edge has exactly one mate, and that mate is

// front-facing. So this is a turning edge.

sel->AddEdge(a, b, auxA);

}

break;

case EdgeKind::EMPHASIZED:

if(info.count == 1 && tr->meta.face != info.tr->meta.face) {

if(CheckAndAddTrianglePair(&edgeTris, tr, info.tr))

break;

// The two triangles that join at this edge come from

// different faces; either really different faces,

// or one is from a face and the other is zero (i.e.,

// not from a face).

sel->AddEdge(a, b, auxA);

}

break;

case EdgeKind::SHARP:

if(info.count == 1) {

Vector na0 = tr->normals[j].WithMagnitude(1.0);

Vector nb0 = tr->normals[(j + 1) % 3].WithMagnitude(1.0);

Vector na1 = info.tr->normals[info.ai].WithMagnitude(1.0);

Vector nb1 = info.tr->normals[info.bi].WithMagnitude(1.0);

if(!((na0.Equals(na1) && nb0.Equals(nb1)) ||

(na0.Equals(nb1) && nb0.Equals(na1)))) {