1 //-----------------------------------------------------------------------------
2 // Operations on triangle meshes, like our mesh Booleans using the BSP, and
3 // the stuff to check for watertightness.
4 //
5 // Copyright 2008-2013 Jonathan Westhues.
6 //-----------------------------------------------------------------------------
7 #include "solvespace.h"
8
9 #include <set>
10
Clear(void)11 void SMesh::Clear(void) {
12 l.Clear();
13 }
14
AddTriangle(STriMeta meta,Vector n,Vector a,Vector b,Vector c)15 void SMesh::AddTriangle(STriMeta meta, Vector n, Vector a, Vector b, Vector c) {
16 Vector ab = b.Minus(a), bc = c.Minus(b);
17 Vector np = ab.Cross(bc);
18 if(np.Magnitude() < 1e-10) {
19 // ugh; gl sometimes tesselates to collinear triangles
20 return;
21 }
22 if(np.Dot(n) > 0) {
23 AddTriangle(meta, a, b, c);
24 } else {
25 AddTriangle(meta, c, b, a);
26 }
27 }
AddTriangle(STriMeta meta,Vector a,Vector b,Vector c)28 void SMesh::AddTriangle(STriMeta meta, Vector a, Vector b, Vector c) {
29 STriangle t = {};
30 t.meta = meta;
31 t.a = a;
32 t.b = b;
33 t.c = c;
34 AddTriangle(&t);
35 }
AddTriangle(STriangle * st)36 void SMesh::AddTriangle(STriangle *st) {
37 RgbaColor color = st->meta.color;
38 if(color.ToARGB32() != 0 && color.alpha != 255) isTransparent = true;
39 l.Add(st);
40 }
41
DoBounding(Vector v,Vector * vmax,Vector * vmin)42 void SMesh::DoBounding(Vector v, Vector *vmax, Vector *vmin) {
43 vmax->x = max(vmax->x, v.x);
44 vmax->y = max(vmax->y, v.y);
45 vmax->z = max(vmax->z, v.z);
46
47 vmin->x = min(vmin->x, v.x);
48 vmin->y = min(vmin->y, v.y);
49 vmin->z = min(vmin->z, v.z);
50 }
GetBounding(Vector * vmax,Vector * vmin)51 void SMesh::GetBounding(Vector *vmax, Vector *vmin) {
52 int i;
53 *vmin = Vector::From( 1e12, 1e12, 1e12);
54 *vmax = Vector::From(-1e12, -1e12, -1e12);
55 for(i = 0; i < l.n; i++) {
56 STriangle *st = &(l.elem[i]);
57 DoBounding(st->a, vmax, vmin);
58 DoBounding(st->b, vmax, vmin);
59 DoBounding(st->c, vmax, vmin);
60 }
61 }
62
63 //----------------------------------------------------------------------------
64 // Report the edges of the boundary of the region(s) of our mesh that lie
65 // within the plane n dot p = d.
66 //----------------------------------------------------------------------------
MakeEdgesInPlaneInto(SEdgeList * sel,Vector n,double d)67 void SMesh::MakeEdgesInPlaneInto(SEdgeList *sel, Vector n, double d) {
68 SMesh m = {};
69 m.MakeFromCopyOf(this);
70
71 // Delete all triangles in the mesh that do not lie in our export plane.
72 m.l.ClearTags();
73 int i;
74 for(i = 0; i < m.l.n; i++) {
75 STriangle *tr = &(m.l.elem[i]);
76
77 if((fabs(n.Dot(tr->a) - d) >= LENGTH_EPS) ||
78 (fabs(n.Dot(tr->b) - d) >= LENGTH_EPS) ||
79 (fabs(n.Dot(tr->c) - d) >= LENGTH_EPS))
80 {
81 tr->tag = 1;
82 }
83 }
84 m.l.RemoveTagged();
85
86 // Select the naked edges in our resulting open mesh.
87 SKdNode *root = SKdNode::From(&m);
88 root->SnapToMesh(&m);
89 root->MakeCertainEdgesInto(sel, SKdNode::NAKED_OR_SELF_INTER_EDGES,
90 false, NULL, NULL);
91
92 m.Clear();
93 }
94
MakeCertainEdgesAndOutlinesInto(SEdgeList * sel,SOutlineList * sol,int type)95 void SMesh::MakeCertainEdgesAndOutlinesInto(SEdgeList *sel, SOutlineList *sol, int type) {
96 SKdNode *root = SKdNode::From(this);
97 root->MakeCertainEdgesInto(sel, type, false, NULL, NULL);
98 root->MakeOutlinesInto(sol);
99 }
100
101 //-----------------------------------------------------------------------------
102 // When we are called, all of the triangles from l.elem[start] to the end must
103 // be coplanar. So we try to find a set of fewer triangles that covers the
104 // exact same area, in order to reduce the number of triangles in the mesh.
105 // We use this after a triangle has been split against the BSP.
106 //
107 // This is really ugly code; basically it just pastes things together to
108 // form convex polygons, merging collinear edges when possible, then
109 // triangulates the convex poly.
110 //-----------------------------------------------------------------------------
Simplify(int start)111 void SMesh::Simplify(int start) {
112 int maxTriangles = (l.n - start) + 10;
113
114 STriMeta meta = l.elem[start].meta;
115
116 STriangle *tout = (STriangle *)AllocTemporary(maxTriangles*sizeof(*tout));
117 int toutc = 0;
118
119 Vector n = Vector::From(0, 0, 0);
120 Vector *conv = (Vector *)AllocTemporary(maxTriangles*3*sizeof(*conv));
121 int convc = 0;
122
123 int start0 = start;
124
125 int i, j;
126 for(i = start; i < l.n; i++) {
127 STriangle *tr = &(l.elem[i]);
128 if(tr->MinAltitude() < LENGTH_EPS) {
129 tr->tag = 1;
130 } else {
131 tr->tag = 0;
132 }
133 }
134
135 for(;;) {
136 bool didAdd;
137 convc = 0;
138 for(i = start; i < l.n; i++) {
139 STriangle *tr = &(l.elem[i]);
140 if(tr->tag) continue;
141
142 tr->tag = 1;
143 n = (tr->Normal()).WithMagnitude(1);
144 conv[convc++] = tr->a;
145 conv[convc++] = tr->b;
146 conv[convc++] = tr->c;
147
148 start = i+1;
149 break;
150 }
151 if(i >= l.n) break;
152
153 do {
154 didAdd = false;
155
156 for(j = 0; j < convc; j++) {
157 Vector a = conv[WRAP((j-1), convc)],
158 b = conv[j],
159 d = conv[WRAP((j+1), convc)],
160 e = conv[WRAP((j+2), convc)];
161
162 Vector c;
163 for(i = start; i < l.n; i++) {
164 STriangle *tr = &(l.elem[i]);
165 if(tr->tag) continue;
166
167 if((tr->a).Equals(d) && (tr->b).Equals(b)) {
168 c = tr->c;
169 } else if((tr->b).Equals(d) && (tr->c).Equals(b)) {
170 c = tr->a;
171 } else if((tr->c).Equals(d) && (tr->a).Equals(b)) {
172 c = tr->b;
173 } else {
174 continue;
175 }
176 // The vertex at C must be convex; but the others must
177 // be tested
178 Vector ab = b.Minus(a);
179 Vector bc = c.Minus(b);
180 Vector cd = d.Minus(c);
181 Vector de = e.Minus(d);
182
183 double bDot = (ab.Cross(bc)).Dot(n);
184 double dDot = (cd.Cross(de)).Dot(n);
185
186 bDot /= min(ab.Magnitude(), bc.Magnitude());
187 dDot /= min(cd.Magnitude(), de.Magnitude());
188
189 if(fabs(bDot) < LENGTH_EPS && fabs(dDot) < LENGTH_EPS) {
190 conv[WRAP((j+1), convc)] = c;
191 // and remove the vertex at j, which is a dup
192 memmove(conv+j, conv+j+1,
193 (convc - j - 1)*sizeof(conv[0]));
194 convc--;
195 } else if(fabs(bDot) < LENGTH_EPS && dDot > 0) {
196 conv[j] = c;
197 } else if(fabs(dDot) < LENGTH_EPS && bDot > 0) {
198 conv[WRAP((j+1), convc)] = c;
199 } else if(bDot > 0 && dDot > 0) {
200 // conv[j] is unchanged, conv[j+1] goes to [j+2]
201 memmove(conv+j+2, conv+j+1,
202 (convc - j - 1)*sizeof(conv[0]));
203 conv[j+1] = c;
204 convc++;
205 } else {
206 continue;
207 }
208
209 didAdd = true;
210 tr->tag = 1;
211 break;
212 }
213 }
214 } while(didAdd);
215
216 // I need to debug why this is required; sometimes the above code
217 // still generates a convex polygon
218 for(i = 0; i < convc; i++) {
219 Vector a = conv[WRAP((i-1), convc)],
220 b = conv[i],
221 c = conv[WRAP((i+1), convc)];
222 Vector ab = b.Minus(a);
223 Vector bc = c.Minus(b);
224 double bDot = (ab.Cross(bc)).Dot(n);
225 bDot /= min(ab.Magnitude(), bc.Magnitude());
226
227 if(bDot < 0) return; // XXX, shouldn't happen
228 }
229
230 for(i = 0; i < convc - 2; i++) {
231 STriangle tr = STriangle::From(meta, conv[0], conv[i+1], conv[i+2]);
232 if(tr.MinAltitude() > LENGTH_EPS) {
233 tout[toutc++] = tr;
234 }
235 }
236 }
237
238 l.n = start0;
239 for(i = 0; i < toutc; i++) {
240 AddTriangle(&(tout[i]));
241 }
242 FreeTemporary(tout);
243 FreeTemporary(conv);
244 }
245
AddAgainstBsp(SMesh * srcm,SBsp3 * bsp3)246 void SMesh::AddAgainstBsp(SMesh *srcm, SBsp3 *bsp3) {
247 int i;
248
249 for(i = 0; i < srcm->l.n; i++) {
250 STriangle *st = &(srcm->l.elem[i]);
251 int pn = l.n;
252 atLeastOneDiscarded = false;
253 SBsp3::InsertOrCreate(bsp3, st, this);
254 if(!atLeastOneDiscarded && (l.n != (pn+1))) {
255 l.n = pn;
256 if(flipNormal) {
257 AddTriangle(st->meta, st->c, st->b, st->a);
258 } else {
259 AddTriangle(st->meta, st->a, st->b, st->c);
260 }
261 }
262 if(l.n - pn > 1) {
263 Simplify(pn);
264 }
265 }
266 }
267
MakeFromUnionOf(SMesh * a,SMesh * b)268 void SMesh::MakeFromUnionOf(SMesh *a, SMesh *b) {
269 SBsp3 *bspa = SBsp3::FromMesh(a);
270 SBsp3 *bspb = SBsp3::FromMesh(b);
271
272 flipNormal = false;
273 keepCoplanar = false;
274 AddAgainstBsp(b, bspa);
275
276 flipNormal = false;
277 keepCoplanar = true;
278 AddAgainstBsp(a, bspb);
279 }
280
MakeFromDifferenceOf(SMesh * a,SMesh * b)281 void SMesh::MakeFromDifferenceOf(SMesh *a, SMesh *b) {
282 SBsp3 *bspa = SBsp3::FromMesh(a);
283 SBsp3 *bspb = SBsp3::FromMesh(b);
284
285 flipNormal = true;
286 keepCoplanar = true;
287 AddAgainstBsp(b, bspa);
288
289 flipNormal = false;
290 keepCoplanar = false;
291 AddAgainstBsp(a, bspb);
292 }
293
MakeFromCopyOf(SMesh * a)294 void SMesh::MakeFromCopyOf(SMesh *a) {
295 int i;
296 for(i = 0; i < a->l.n; i++) {
297 AddTriangle(&(a->l.elem[i]));
298 }
299 }
300
MakeFromAssemblyOf(SMesh * a,SMesh * b)301 void SMesh::MakeFromAssemblyOf(SMesh *a, SMesh *b) {
302 MakeFromCopyOf(a);
303 MakeFromCopyOf(b);
304 }
305
MakeFromTransformationOf(SMesh * a,Vector trans,Quaternion q,double scale)306 void SMesh::MakeFromTransformationOf(SMesh *a,
307 Vector trans, Quaternion q, double scale)
308 {
309 STriangle *tr;
310 for(tr = a->l.First(); tr; tr = a->l.NextAfter(tr)) {
311 STriangle tt = *tr;
312 tt.a = (tt.a).ScaledBy(scale);
313 tt.b = (tt.b).ScaledBy(scale);
314 tt.c = (tt.c).ScaledBy(scale);
315 if(scale < 0) {
316 // The mirroring would otherwise turn a closed mesh inside out.
317 swap(tt.a, tt.b);
318 }
319 tt.a = (q.Rotate(tt.a)).Plus(trans);
320 tt.b = (q.Rotate(tt.b)).Plus(trans);
321 tt.c = (q.Rotate(tt.c)).Plus(trans);
322 AddTriangle(&tt);
323 }
324 }
325
IsEmpty(void)326 bool SMesh::IsEmpty(void) {
327 return (l.n == 0);
328 }
329
FirstIntersectionWith(Point2d mp)330 uint32_t SMesh::FirstIntersectionWith(Point2d mp) {
331 Vector p0 = Vector::From(mp.x, mp.y, 0);
332 Vector gn = Vector::From(0, 0, 1);
333
334 double maxT = -1e12;
335 uint32_t face = 0;
336
337 int i;
338 for(i = 0; i < l.n; i++) {
339 STriangle tr = l.elem[i];
340 tr.a = SS.GW.ProjectPoint3(tr.a);
341 tr.b = SS.GW.ProjectPoint3(tr.b);
342 tr.c = SS.GW.ProjectPoint3(tr.c);
343
344 Vector n = tr.Normal();
345
346 if(n.Dot(gn) < LENGTH_EPS) continue; // back-facing or on edge
347
348 if(tr.ContainsPointProjd(gn, p0)) {
349 // Let our line have the form r(t) = p0 + gn*t
350 double t = -(n.Dot((tr.a).Minus(p0)))/(n.Dot(gn));
351 if(t > maxT) {
352 maxT = t;
353 face = tr.meta.face;
354 }
355 }
356 }
357 return face;
358 }
359
Alloc(void)360 STriangleLl *STriangleLl::Alloc(void)
361 { return (STriangleLl *)AllocTemporary(sizeof(STriangleLl)); }
Alloc(void)362 SKdNode *SKdNode::Alloc(void)
363 { return (SKdNode *)AllocTemporary(sizeof(SKdNode)); }
364
From(SMesh * m)365 SKdNode *SKdNode::From(SMesh *m) {
366 int i;
367 STriangle *tra = (STriangle *)AllocTemporary((m->l.n) * sizeof(*tra));
368
369 for(i = 0; i < m->l.n; i++) {
370 tra[i] = m->l.elem[i];
371 }
372
373 srand(0);
374 int n = m->l.n;
375 while(n > 1) {
376 int k = rand() % n;
377 n--;
378 swap(tra[k], tra[n]);
379 }
380
381 STriangleLl *tll = NULL;
382 for(i = 0; i < m->l.n; i++) {
383 STriangleLl *tn = STriangleLl::Alloc();
384 tn->tri = &(tra[i]);
385 tn->next = tll;
386 tll = tn;
387 }
388
389 return SKdNode::From(tll);
390 }
391
From(STriangleLl * tll)392 SKdNode *SKdNode::From(STriangleLl *tll) {
393 int which = 0;
394 SKdNode *ret = Alloc();
395
396 int i;
397 int gtc[3] = { 0, 0, 0 }, ltc[3] = { 0, 0, 0 }, allc = 0;
398 double badness[3] = { 0, 0, 0 };
399 double split[3] = { 0, 0, 0 };
400
401 if(!tll) {
402 goto leaf;
403 }
404
405 for(i = 0; i < 3; i++) {
406 int tcnt = 0;
407 STriangleLl *ll;
408 for(ll = tll; ll; ll = ll->next) {
409 split[i] += (ll->tri->a).Element(i);
410 split[i] += (ll->tri->b).Element(i);
411 split[i] += (ll->tri->c).Element(i);
412 tcnt++;
413 }
414 split[i] /= (tcnt*3);
415
416 for(ll = tll; ll; ll = ll->next) {
417 STriangle *tr = ll->tri;
418
419 double a = (tr->a).Element(i),
420 b = (tr->b).Element(i),
421 c = (tr->c).Element(i);
422
423 if(a < split[i] + KDTREE_EPS ||
424 b < split[i] + KDTREE_EPS ||
425 c < split[i] + KDTREE_EPS)
426 {
427 ltc[i]++;
428 }
429 if(a > split[i] - KDTREE_EPS ||
430 b > split[i] - KDTREE_EPS ||
431 c > split[i] - KDTREE_EPS)
432 {
433 gtc[i]++;
434 }
435 if(i == 0) allc++;
436 }
437 badness[i] = pow((double)ltc[i], 4) + pow((double)gtc[i], 4);
438 }
439 if(badness[0] < badness[1] && badness[0] < badness[2]) {
440 which = 0;
441 } else if(badness[1] < badness[2]) {
442 which = 1;
443 } else {
444 which = 2;
445 }
446
447 if(allc < 3 || allc == gtc[which] || allc == ltc[which]) {
448 goto leaf;
449 }
450
451 STriangleLl *ll;
452 STriangleLl *lgt, *llt; lgt = llt = NULL;
453 for(ll = tll; ll; ll = ll->next) {
454 STriangle *tr = ll->tri;
455
456 double a = (tr->a).Element(which),
457 b = (tr->b).Element(which),
458 c = (tr->c).Element(which);
459
460 if(a < split[which] + KDTREE_EPS ||
461 b < split[which] + KDTREE_EPS ||
462 c < split[which] + KDTREE_EPS)
463 {
464 STriangleLl *n = STriangleLl::Alloc();
465 *n = *ll;
466 n->next = llt;
467 llt = n;
468 }
469 if(a > split[which] - KDTREE_EPS ||
470 b > split[which] - KDTREE_EPS ||
471 c > split[which] - KDTREE_EPS)
472 {
473 STriangleLl *n = STriangleLl::Alloc();
474 *n = *ll;
475 n->next = lgt;
476 lgt = n;
477 }
478 }
479
480 ret->which = which;
481 ret->c = split[which];
482 ret->gt = SKdNode::From(lgt);
483 ret->lt = SKdNode::From(llt);
484 return ret;
485
486 leaf:
487 ret->tris = tll;
488 return ret;
489 }
490
ClearTags(void)491 void SKdNode::ClearTags(void) {
492 if(gt && lt) {
493 gt->ClearTags();
494 lt->ClearTags();
495 } else {
496 STriangleLl *ll;
497 for(ll = tris; ll; ll = ll->next) {
498 ll->tri->tag = 0;
499 }
500 }
501 }
502
AddTriangle(STriangle * tr)503 void SKdNode::AddTriangle(STriangle *tr) {
504 if(gt && lt) {
505 double ta = (tr->a).Element(which),
506 tb = (tr->b).Element(which),
507 tc = (tr->c).Element(which);
508 if(ta < c + KDTREE_EPS ||
509 tb < c + KDTREE_EPS ||
510 tc < c + KDTREE_EPS)
511 {
512 lt->AddTriangle(tr);
513 }
514 if(ta > c - KDTREE_EPS ||
515 tb > c - KDTREE_EPS ||
516 tc > c - KDTREE_EPS)
517 {
518 gt->AddTriangle(tr);
519 }
520 } else {
521 STriangleLl *tn = STriangleLl::Alloc();
522 tn->tri = tr;
523 tn->next = tris;
524 tris = tn;
525 }
526 }
527
MakeMeshInto(SMesh * m)528 void SKdNode::MakeMeshInto(SMesh *m) {
529 if(gt) gt->MakeMeshInto(m);
530 if(lt) lt->MakeMeshInto(m);
531
532 STriangleLl *ll;
533 for(ll = tris; ll; ll = ll->next) {
534 if(ll->tri->tag) continue;
535
536 m->AddTriangle(ll->tri);
537 ll->tri->tag = 1;
538 }
539 }
540
ListTrianglesInto(std::vector<STriangle * > * tl)541 void SKdNode::ListTrianglesInto(std::vector<STriangle *> *tl) {
542 if(gt) gt->ListTrianglesInto(tl);
543 if(lt) lt->ListTrianglesInto(tl);
544
545 STriangleLl *ll;
546 for(ll = tris; ll; ll = ll->next) {
547 if(ll->tri->tag) continue;
548
549 tl->push_back(ll->tri);
550 ll->tri->tag = 1;
551 }
552 }
553
554 //-----------------------------------------------------------------------------
555 // If any triangles in the mesh have an edge that goes through v (but not
556 // a vertex at v), then split those triangles so that they now have a vertex
557 // there. The existing triangle is modified, and the new triangle appears
558 // in extras.
559 //-----------------------------------------------------------------------------
SnapToVertex(Vector v,SMesh * extras)560 void SKdNode::SnapToVertex(Vector v, SMesh *extras) {
561 if(gt && lt) {
562 double vc = v.Element(which);
563 if(vc < c + KDTREE_EPS) {
564 lt->SnapToVertex(v, extras);
565 }
566 if(vc > c - KDTREE_EPS) {
567 gt->SnapToVertex(v, extras);
568 }
569 // Nothing bad happens if the triangle to be split appears in both
570 // branches; the first call will split the triangle, so that the
571 // second call will do nothing, because the modified triangle will
572 // already contain v
573 } else {
574 STriangleLl *ll;
575 for(ll = tris; ll; ll = ll->next) {
576 STriangle *tr = ll->tri;
577
578 // Do a cheap bbox test first
579 int k;
580 bool mightHit = true;
581
582 for(k = 0; k < 3; k++) {
583 if((tr->a).Element(k) < v.Element(k) - KDTREE_EPS &&
584 (tr->b).Element(k) < v.Element(k) - KDTREE_EPS &&
585 (tr->c).Element(k) < v.Element(k) - KDTREE_EPS)
586 {
587 mightHit = false;
588 break;
589 }
590 if((tr->a).Element(k) > v.Element(k) + KDTREE_EPS &&
591 (tr->b).Element(k) > v.Element(k) + KDTREE_EPS &&
592 (tr->c).Element(k) > v.Element(k) + KDTREE_EPS)
593 {
594 mightHit = false;
595 break;
596 }
597 }
598 if(!mightHit) continue;
599
600 if(tr->a.Equals(v)) { tr->a = v; continue; }
601 if(tr->b.Equals(v)) { tr->b = v; continue; }
602 if(tr->c.Equals(v)) { tr->c = v; continue; }
603
604 if(v.OnLineSegment(tr->a, tr->b)) {
605 STriangle nt = STriangle::From(tr->meta, tr->a, v, tr->c);
606 extras->AddTriangle(&nt);
607 tr->a = v;
608 continue;
609 }
610 if(v.OnLineSegment(tr->b, tr->c)) {
611 STriangle nt = STriangle::From(tr->meta, tr->b, v, tr->a);
612 extras->AddTriangle(&nt);
613 tr->b = v;
614 continue;
615 }
616 if(v.OnLineSegment(tr->c, tr->a)) {
617 STriangle nt = STriangle::From(tr->meta, tr->c, v, tr->b);
618 extras->AddTriangle(&nt);
619 tr->c = v;
620 continue;
621 }
622 }
623 }
624 }
625
626 //-----------------------------------------------------------------------------
627 // Snap to each vertex of each triangle of the given mesh. If the given mesh
628 // is identical to the mesh used to make this kd tree, then the result should
629 // be a vertex-to-vertex mesh.
630 //-----------------------------------------------------------------------------
SnapToMesh(SMesh * m)631 void SKdNode::SnapToMesh(SMesh *m) {
632 int i, j, k;
633 for(i = 0; i < m->l.n; i++) {
634 STriangle *tr = &(m->l.elem[i]);
635 for(j = 0; j < 3; j++) {
636 Vector v = ((j == 0) ? tr->a :
637 ((j == 1) ? tr->b :
638 tr->c));
639
640 SMesh extra = {};
641 SnapToVertex(v, &extra);
642
643 for(k = 0; k < extra.l.n; k++) {
644 STriangle *tra = (STriangle *)AllocTemporary(sizeof(*tra));
645 *tra = extra.l.elem[k];
646 AddTriangle(tra);
647 }
648 extra.Clear();
649 }
650 }
651 }
652
653 //-----------------------------------------------------------------------------
654 // For all the edges in sel, split them against the given triangle, and test
655 // them for occlusion. Keep only the visible segments. sel is both our input
656 // and our output.
657 //-----------------------------------------------------------------------------
SplitLinesAgainstTriangle(SEdgeList * sel,STriangle * tr,bool removeHidden)658 void SKdNode::SplitLinesAgainstTriangle(SEdgeList *sel, STriangle *tr, bool removeHidden) {
659 SEdgeList seln = {};
660
661 Vector tn = tr->Normal().WithMagnitude(1);
662 double td = tn.Dot(tr->a);
663
664 // Consider front-facing triangles only.
665 if(tn.z > LENGTH_EPS) {
666 // If the edge crosses our triangle's plane, then split into above
667 // and below parts. Note that we must preserve auxA, which contains
668 // the style associated with this line.
669 SEdge *se;
670 for(se = sel->l.First(); se; se = sel->l.NextAfter(se)) {
671 double da = (se->a).Dot(tn) - td,
672 db = (se->b).Dot(tn) - td;
673 if((da < -LENGTH_EPS && db > LENGTH_EPS) ||
674 (db < -LENGTH_EPS && da > LENGTH_EPS))
675 {
676 Vector m = Vector::AtIntersectionOfPlaneAndLine(
677 tn, td,
678 se->a, se->b, NULL);
679 seln.AddEdge(m, se->b, se->auxA);
680 se->b = m;
681 }
682 }
683 for(se = seln.l.First(); se; se = seln.l.NextAfter(se)) {
684 sel->AddEdge(se->a, se->b, se->auxA);
685 }
686 seln.Clear();
687
688 for(se = sel->l.First(); se; se = sel->l.NextAfter(se)) {
689 Vector pt = ((se->a).Plus(se->b)).ScaledBy(0.5);
690 if(pt.Dot(tn) - td > -LENGTH_EPS) {
691 // Edge is in front of or on our plane (remember, tn.z > 0)
692 // so it is exempt from further splitting
693 se->auxB = 1;
694 } else {
695 // Edge is behind our plane, needs further splitting
696 se->auxB = 0;
697 }
698 }
699
700 // Considering only the (x, y) coordinates, split the edge against our
701 // triangle.
702 Point2d a = (tr->a).ProjectXy(),
703 b = (tr->b).ProjectXy(),
704 c = (tr->c).ProjectXy();
705
706 Point2d n[3] = { (b.Minus(a)).Normal().WithMagnitude(1),
707 (c.Minus(b)).Normal().WithMagnitude(1),
708 (a.Minus(c)).Normal().WithMagnitude(1) };
709
710 double d[3] = { n[0].Dot(b),
711 n[1].Dot(c),
712 n[2].Dot(a) };
713
714 // Split all of the edges where they intersect the triangle edges
715 int i;
716 for(i = 0; i < 3; i++) {
717 for(se = sel->l.First(); se; se = sel->l.NextAfter(se)) {
718 if(se->auxB) continue;
719
720 Point2d ap = (se->a).ProjectXy(),
721 bp = (se->b).ProjectXy();
722 double da = n[i].Dot(ap) - d[i],
723 db = n[i].Dot(bp) - d[i];
724 if((da < -LENGTH_EPS && db > LENGTH_EPS) ||
725 (db < -LENGTH_EPS && da > LENGTH_EPS))
726 {
727 double dab = (db - da);
728 Vector spl = ((se->a).ScaledBy( db/dab)).Plus(
729 (se->b).ScaledBy(-da/dab));
730 seln.AddEdge(spl, se->b, se->auxA);
731 se->b = spl;
732 }
733 }
734 for(se = seln.l.First(); se; se = seln.l.NextAfter(se)) {
735 // The split pieces are all behind the triangle, since only
736 // edges behind the triangle got split. So their auxB is 0.
737 sel->AddEdge(se->a, se->b, se->auxA, 0);
738 }
739 seln.Clear();
740 }
741
742 for(se = sel->l.First(); se; se = sel->l.NextAfter(se)) {
743 if(se->auxB) {
744 // Lies above or on the triangle plane, so triangle doesn't
745 // occlude it.
746 se->tag = 0;
747 } else {
748 // Test the segment to see if it lies outside the triangle
749 // (i.e., outside wrt at least one edge), and keep it only
750 // then.
751 Point2d pt = ((se->a).Plus(se->b).ScaledBy(0.5)).ProjectXy();
752 se->tag = 1;
753 for(i = 0; i < 3; i++) {
754 // If the test point lies on the boundary of our triangle,
755 // then we still discard the edge.
756 if(n[i].Dot(pt) - d[i] > LENGTH_EPS) se->tag = 0;
757 }
758 }
759 if(!removeHidden && se->tag == 1)
760 se->auxA = Style::HIDDEN_EDGE;
761 }
762 if(removeHidden)
763 sel->l.RemoveTagged();
764 }
765 }
766
767 //-----------------------------------------------------------------------------
768 // Given an edge orig, occlusion test it against our mesh. We output an edge
769 // list in sel, containing the visible portions of that edge.
770 //-----------------------------------------------------------------------------
OcclusionTestLine(SEdge orig,SEdgeList * sel,int cnt,bool removeHidden)771 void SKdNode::OcclusionTestLine(SEdge orig, SEdgeList *sel, int cnt, bool removeHidden) {
772 if(gt && lt) {
773 double ac = (orig.a).Element(which),
774 bc = (orig.b).Element(which);
775 // We can ignore triangles that are separated in x or y, but triangles
776 // that are separated in z may still contribute
777 if(ac < c + KDTREE_EPS ||
778 bc < c + KDTREE_EPS ||
779 which == 2)
780 {
781 lt->OcclusionTestLine(orig, sel, cnt, removeHidden);
782 }
783 if(ac > c - KDTREE_EPS ||
784 bc > c - KDTREE_EPS ||
785 which == 2)
786 {
787 gt->OcclusionTestLine(orig, sel, cnt, removeHidden);
788 }
789 } else {
790 STriangleLl *ll;
791 for(ll = tris; ll; ll = ll->next) {
792 STriangle *tr = ll->tri;
793
794 if(tr->tag == cnt) continue;
795
796 SplitLinesAgainstTriangle(sel, tr, removeHidden);
797 tr->tag = cnt;
798 }
799 }
800 }
801
802 //-----------------------------------------------------------------------------
803 // Search the mesh for a triangle with an edge from b to a (i.e., the mate
804 // for the edge from a to b), and increment info->count each time that we
805 // find one. If a triangle is found, then report whether it is front- or
806 // back-facing using info->frontFacing. And regardless of whether a mate is
807 // found, report whether the edge intersects the mesh with info->intersectsMesh;
808 // if coplanarIsInter then we count the edge as intersecting if it's coplanar
809 // with a triangle in the mesh, otherwise not.
810 //-----------------------------------------------------------------------------
FindEdgeOn(Vector a,Vector b,int cnt,bool coplanarIsInter,EdgeOnInfo * info)811 void SKdNode::FindEdgeOn(Vector a, Vector b, int cnt, bool coplanarIsInter,
812 EdgeOnInfo *info)
813 {
814 if(gt && lt) {
815 double ac = a.Element(which),
816 bc = b.Element(which);
817 if(ac < c + KDTREE_EPS ||
818 bc < c + KDTREE_EPS)
819 {
820 lt->FindEdgeOn(a, b, cnt, coplanarIsInter, info);
821 }
822 if(ac > c - KDTREE_EPS ||
823 bc > c - KDTREE_EPS)
824 {
825 gt->FindEdgeOn(a, b, cnt, coplanarIsInter, info);
826 }
827 return;
828 }
829
830 // We are a leaf node; so we iterate over all the triangles in our
831 // linked list.
832 STriangleLl *ll;
833 for(ll = tris; ll; ll = ll->next) {
834 STriangle *tr = ll->tri;
835
836 if(tr->tag == cnt) continue;
837
838 // Test if this triangle matches up with the given edge
839 if((a.Equals(tr->b) && b.Equals(tr->a)) ||
840 (a.Equals(tr->c) && b.Equals(tr->b)) ||
841 (a.Equals(tr->a) && b.Equals(tr->c)))
842 {
843 info->count++;
844 // Record whether this triangle is front- or back-facing.
845 if(tr->Normal().z > LENGTH_EPS) {
846 info->frontFacing = true;
847 } else {
848 info->frontFacing = false;
849 }
850 // Record the triangle
851 info->tr = tr;
852 // And record which vertexes a and b correspond to
853 info->ai = a.Equals(tr->a) ? 0 : (a.Equals(tr->b) ? 1 : 2);
854 info->bi = b.Equals(tr->a) ? 0 : (b.Equals(tr->b) ? 1 : 2);
855 } else if(((a.Equals(tr->a) && b.Equals(tr->b)) ||
856 (a.Equals(tr->b) && b.Equals(tr->c)) ||
857 (a.Equals(tr->c) && b.Equals(tr->a))))
858 {
859 // It's an edge of this triangle, okay.
860 } else {
861 // Check for self-intersection
862 Vector n = (tr->Normal()).WithMagnitude(1);
863 double d = (tr->a).Dot(n);
864 double pa = a.Dot(n) - d, pb = b.Dot(n) - d;
865 // It's an intersection if neither point lies in-plane,
866 // and the edge crosses the plane (should handle in-plane
867 // intersections separately but don't yet).
868 if((pa < -LENGTH_EPS || pa > LENGTH_EPS) &&
869 (pb < -LENGTH_EPS || pb > LENGTH_EPS) &&
870 (pa*pb < 0))
871 {
872 // The edge crosses the plane of the triangle; now see if
873 // it crosses inside the triangle.
874 if(tr->ContainsPointProjd(b.Minus(a), a)) {
875 if(coplanarIsInter) {
876 info->intersectsMesh = true;
877 } else {
878 Vector p = Vector::AtIntersectionOfPlaneAndLine(
879 n, d, a, b, NULL);
880 Vector ta = tr->a,
881 tb = tr->b,
882 tc = tr->c;
883 if((p.DistanceToLine(ta, tb.Minus(ta)) < LENGTH_EPS) ||
884 (p.DistanceToLine(tb, tc.Minus(tb)) < LENGTH_EPS) ||
885 (p.DistanceToLine(tc, ta.Minus(tc)) < LENGTH_EPS))
886 {
887 // Intersection lies on edge. This happens when
888 // our edge is from a triangle coplanar with
889 // another triangle in the mesh. We don't test
890 // the edge against triangles whose plane contains
891 // that edge, but we do end up testing against
892 // the coplanar triangle's neighbours, which we
893 // will intersect on their edges.
894 } else {
895 info->intersectsMesh = true;
896 }
897 }
898 }
899 }
900 }
901
902 // Ensure that we don't count this triangle twice if it appears
903 // in two buckets of the kd tree.
904 tr->tag = cnt;
905 }
906 }
907
CheckAndAddTrianglePair(std::set<std::pair<STriangle *,STriangle * >> * pairs,STriangle * a,STriangle * b)908 static bool CheckAndAddTrianglePair(std::set<std::pair<STriangle *, STriangle *>> *pairs,
909 STriangle *a, STriangle *b)
910 {
911 if(pairs->find(std::make_pair(a, b)) != pairs->end() ||
912 pairs->find(std::make_pair(b, a)) != pairs->end())
913 return true;
914
915 pairs->emplace(a, b);
916 return false;
917 }
918
919 //-----------------------------------------------------------------------------
920 // Pick certain classes of edges out from our mesh. These might be:
921 // * naked edges (i.e., edges with no anti-parallel neighbor) and self-
922 // intersecting edges (i.e., edges that cross another triangle)
923 // * turning edges (i.e., edges where a front-facing triangle joins
924 // a back-facing triangle)
925 // * emphasized edges (i.e., edges where a triangle from one face joins
926 // a triangle from a different face)
927 //-----------------------------------------------------------------------------
MakeCertainEdgesInto(SEdgeList * sel,int how,bool coplanarIsInter,bool * inter,bool * leaky,int auxA)928 void SKdNode::MakeCertainEdgesInto(SEdgeList *sel, int how, bool coplanarIsInter,
929 bool *inter, bool *leaky, int auxA)
930 {
931 if(inter) *inter = false;
932 if(leaky) *leaky = false;
933
934 std::vector<STriangle *> tris;
935 ClearTags();
936 ListTrianglesInto(&tris);
937
938 std::set<std::pair<STriangle *, STriangle *>> edgeTris;
939 int cnt = 1234;
940 for(STriangle *tr : tris) {
941 for(int j = 0; j < 3; j++) {
942 Vector a = (j == 0) ? tr->a : ((j == 1) ? tr->b : tr->c);
943 Vector b = (j == 0) ? tr->b : ((j == 1) ? tr->c : tr->a);
944
945 SKdNode::EdgeOnInfo info = {};
946 FindEdgeOn(a, b, cnt, coplanarIsInter, &info);
947
948 switch(how) {
949 case NAKED_OR_SELF_INTER_EDGES:
950 if(info.count != 1) {
951 sel->AddEdge(a, b, auxA);
952 if(leaky) *leaky = true;
953 }
954 if(info.intersectsMesh) {
955 sel->AddEdge(a, b, auxA);
956 if(inter) *inter = true;
957 }
958 break;
959
960 case SELF_INTER_EDGES:
961 if(info.intersectsMesh) {
962 sel->AddEdge(a, b, auxA);
963 if(inter) *inter = true;
964 }
965 break;
966
967 case TURNING_EDGES:
968 if((tr->Normal().z < LENGTH_EPS) &&
969 (info.count == 1) &&
970 info.frontFacing)
971 {
972 if(CheckAndAddTrianglePair(&edgeTris, tr, info.tr))
973 break;
974 // This triangle is back-facing (or on edge), and
975 // this edge has exactly one mate, and that mate is
976 // front-facing. So this is a turning edge.
977 sel->AddEdge(a, b, auxA);
978 }
979 break;
980
981 case EMPHASIZED_EDGES:
982 if(info.count == 1 && tr->meta.face != info.tr->meta.face) {
983 if(CheckAndAddTrianglePair(&edgeTris, tr, info.tr))
984 break;
985 // The two triangles that join at this edge come from
986 // different faces; either really different faces,
987 // or one is from a face and the other is zero (i.e.,
988 // not from a face).
989 sel->AddEdge(a, b, auxA);
990 }
991 break;
992
993 case SHARP_EDGES:
994 if(info.count == 1) {
995 Vector na0 = (j == 0) ? tr->an :
996 ((j == 1) ? tr->bn : tr->cn);
997 Vector nb0 = (j == 0) ? tr->bn :
998 ((j == 1) ? tr->cn : tr->an);
999 Vector na1 = (info.ai == 0) ? info.tr->an :
1000 ((info.ai == 1) ? info.tr->bn : info.tr->cn);
1001 Vector nb1 = (info.bi == 0) ? info.tr->an :
1002 ((info.bi == 1) ? info.tr->bn : info.tr->cn);
1003 na0 = na0.WithMagnitude(1.0);
1004 nb0 = nb0.WithMagnitude(1.0);
1005 na1 = na1.WithMagnitude(1.0);
1006 nb1 = nb1.WithMagnitude(1.0);
1007 if(!((na0.Equals(na1) && nb0.Equals(nb1)) ||
1008 (na0.Equals(nb1) && nb0.Equals(na1)))) {
1009 if(CheckAndAddTrianglePair(&edgeTris, tr, info.tr))
1010 break;
1011 // The two triangles that join at this edge meet at a sharp
1012 // angle. This implies they come from different faces.
1013 sel->AddEdge(a, b, auxA);
1014 }
1015 }
1016 break;
1017
1018 default: oops();
1019 }
1020
1021 cnt++;
1022 }
1023 }
1024 }
1025
MakeOutlinesInto(SOutlineList * sol)1026 void SKdNode::MakeOutlinesInto(SOutlineList *sol)
1027 {
1028 std::vector<STriangle *> tris;
1029 ClearTags();
1030 ListTrianglesInto(&tris);
1031
1032 std::set<std::pair<STriangle *, STriangle *>> edgeTris;
1033 int cnt = 1234;
1034 for(STriangle *tr : tris) {
1035 for(int j = 0; j < 3; j++) {
1036 Vector a = (j == 0) ? tr->a : ((j == 1) ? tr->b : tr->c);
1037 Vector b = (j == 0) ? tr->b : ((j == 1) ? tr->c : tr->a);
1038
1039 SKdNode::EdgeOnInfo info = {};
1040 FindEdgeOn(a, b, cnt, /*coplanarIsInter=*/false, &info);
1041 cnt++;
1042 if(info.count != 1) continue;
1043 if(CheckAndAddTrianglePair(&edgeTris, tr, info.tr))
1044 continue;
1045
1046 sol->AddEdge(a, b, tr->Normal(), info.tr->Normal());
1047 }
1048 }
1049 }
1050
Clear()1051 void SOutlineList::Clear() {
1052 l.Clear();
1053 }
1054
AddEdge(Vector a,Vector b,Vector nl,Vector nr)1055 void SOutlineList::AddEdge(Vector a, Vector b, Vector nl, Vector nr) {
1056 SOutline so = {};
1057 so.a = a;
1058 so.b = b;
1059 so.nl = nl;
1060 so.nr = nr;
1061 l.Add(&so);
1062 }
1063
FillOutlineTags(Vector projDir)1064 void SOutlineList::FillOutlineTags(Vector projDir) {
1065 for(SOutline *so = l.First(); so; so = l.NextAfter(so)) {
1066 so->tag = ((so->nl.Dot(projDir) > 0.0) != (so->nr.Dot(projDir) > 0.0));
1067 }
1068 }
1069
MakeFromCopyOf(SOutlineList * sol)1070 void SOutlineList::MakeFromCopyOf(SOutlineList *sol) {
1071 for(SOutline *so = sol->l.First(); so; so = sol->l.NextAfter(so)) {
1072 l.Add(so);
1073 }
1074 }
1075