ff48f424e3c324493150646b4ebcfafb1349adef
[anneal.git] / src / edu / berkeley / qfat / Mesh.java
1 package edu.berkeley.qfat;
2 import java.awt.*;
3 import java.util.*;
4 import java.awt.event.*;
5 import javax.swing.*;
6 import javax.media.opengl.*;
7 import javax.media.opengl.glu.*;
8 import edu.berkeley.qfat.geom.*;
9 import edu.wlu.cs.levy.CG.KDTree;
10 import edu.berkeley.qfat.geom.Point;
11
12 public class Mesh implements Iterable<Mesh.T> {
13
14     public static final float EPSILON = (float)0.0001;
15     public static final Random random = new Random();
16
17     private RTree<T> tris = new RTree<T>();
18     private PointSet<Vert> vertices = new PointSet<Vert>();
19
20     public boolean tilemesh = false;
21     public boolean ignorecollision = false;
22     public Mesh score_against = null;
23     public double score = 0;
24     public float score() { return (float)score; }
25
26     public int size() { return vertices.size(); }
27     public Iterable<Vert> vertices() { return vertices; }
28
29     public Iterator<T> iterator() {
30         return tris.iterator();
31     }
32
33     public void rebindPoints() {
34         // unbind all points
35         for(Mesh.T t : this) {
36             t.v1().unbind();
37             t.v2().unbind();
38             t.v3().unbind();
39         }
40         // ask edges to re-implement their bindings
41         for(Mesh.T t : this) {
42             t.e1().dobind();
43             t.e2().dobind();
44             t.e3().dobind();
45         }
46     }
47
48     public void unApplyQuadricToNeighborAll() {
49         HashSet<Vert> done = new HashSet<Vert>();
50         for(T t : this)
51             for(Vert p : new Vert[] { t.v1(), t.v2(), t.v3() }) {
52                 if (done.contains(p)) continue;
53                 done.add(p);
54                 p.unApplyQuadricToNeighbor();
55             }
56     }
57     public void recomputeAllFundamentalQuadrics() {
58         HashSet<Vert> done = new HashSet<Vert>();
59         for(T t : this)
60             for(Vert p : new Vert[] { t.v1(), t.v2(), t.v3() }) {
61                 if (done.contains(p)) continue;
62                 done.add(p);
63                 p.recomputeFundamentalQuadric();
64             }
65     }
66     public float applyQuadricToNeighborAll() {
67         int num = 0;
68         double dist = 0;
69         HashSet<Vert> done = new HashSet<Vert>();
70         for(T t : this)
71             for(Vert p : new Vert[] { t.v1(), t.v2(), t.v3() }) {
72                 if (done.contains(p)) continue;
73                 done.add(p);
74                 p.applyQuadricToNeighbor();
75                 
76             }
77         return (float)(dist/num);
78     }
79
80     public void transform(Matrix m) {
81         ArrayList<Vert> set = new ArrayList<Vert>();
82         for (Vert v : vertices)
83             set.add(v);
84         for(Vert v : set) v.transform(m);
85     }
86
87     public float volume() {
88         double total = 0;
89         for(T t : this) {
90             double area = t.area();
91             Vec origin_to_centroid = new Vec(new Point(0, 0, 0), t.centroid());
92             boolean facingAway = t.norm().dot(origin_to_centroid) > 0;
93             double height = Math.abs(t.norm().dot(origin_to_centroid));
94             total += ((facingAway ? 1 : -1) * area * height) / 3.0;
95         }
96         return (float)total;
97     }
98
99     public void rebuild() { /*vertices.rebuild();*/ }
100     public Vec diagonal() { return vertices.diagonal(); }
101     public Point centroid() { return vertices.centroid(); }
102     public Vert nearest(Point p) { return vertices.nearest(p); }
103
104     public final class Vert extends HasPoint {
105         public String toString() { return p.toString(); }
106         public Point p;
107         E e;                // some edge *leaving* this point
108
109         /** the nearest vertex in the "score_against" mesh */
110         Vert   nearest_in_other_mesh;
111         /** the number of vertices in the other mesh for which this is the nearest_in_other_mesh */
112         int    quadric_count;
113         /** the total error quadric (contributions from all vertices in other mesh for which this is nearest) */
114         Matrix quadric = Matrix.ZERO;
115
116         Vert bound_to = this;
117         Matrix binding = Matrix.ONE;
118         float oldscore = 0;
119         boolean quadricStale = false;
120
121         public Matrix errorQuadric() { return quadric; }
122         public Point getPoint() { return p; }
123         public float score() { return oldscore; }
124
125         private Matrix fundamentalQuadric = null;
126         public Matrix fundamentalQuadric() {
127             if (fundamentalQuadric == null) recomputeFundamentalQuadric();
128             return fundamentalQuadric;
129         }
130
131         private Vert(Point p) {
132             this.p = p;
133             if (vertices.get(p) != null) throw new Error();
134             vertices.add(this);
135         }
136
137         private void glNormal(GL gl) {
138             Vec norm = norm();
139             gl.glNormal3f(norm.x, norm.y, norm.z);
140         }
141
142         public void recomputeFundamentalQuadric() {
143             //if (!quadricStale && fundamentalQuadric != null) return;
144             quadricStale = false;
145             unApplyQuadricToNeighbor();
146             Matrix m = Matrix.ZERO;
147             E e = this.e;
148             int count = 0;
149             do {
150                 T t = e.t;
151                 m = m.plus(t.norm().fundamentalQuadric(t.centroid()));
152                 count++;
153                 e = e.pair.next;
154             } while(e != this.e);
155             fundamentalQuadric = m.times(1/(float)count);
156             applyQuadricToNeighbor();
157         }
158
159         public void unApplyQuadricToNeighbor() {
160             if (nearest_in_other_mesh == null) return;
161             if (fundamentalQuadric == null) return;
162             nearest_in_other_mesh.unComputeError();
163             nearest_in_other_mesh.quadric = nearest_in_other_mesh.quadric.minus(fundamentalQuadric);
164             nearest_in_other_mesh.quadric_count--;
165             if (nearest_in_other_mesh.quadric_count==0)
166                 nearest_in_other_mesh.quadric = Matrix.ZERO;
167             nearest_in_other_mesh.computeError();
168             nearest_in_other_mesh = null;
169         }
170
171         public void applyQuadricToNeighbor() {
172             if (score_against == null) return;
173
174             Vert new_nearest = score_against.nearest(p);
175             if (nearest_in_other_mesh != null && new_nearest == nearest_in_other_mesh) return;
176
177             if (nearest_in_other_mesh != null) unApplyQuadricToNeighbor();
178             if (nearest_in_other_mesh != null) throw new Error();
179
180             nearest_in_other_mesh = new_nearest;
181                 
182             // don't attract to vertices that face the other way
183             if (nearest_in_other_mesh.e == null || nearest_in_other_mesh.norm().dot(norm()) < 0) {
184                 nearest_in_other_mesh = null;
185             } else {
186                 nearest_in_other_mesh.unComputeError();
187                 nearest_in_other_mesh.quadric = nearest_in_other_mesh.quadric.plus(fundamentalQuadric());
188                 nearest_in_other_mesh.quadric_count++;
189                 nearest_in_other_mesh.computeError();
190             }
191             reComputeError();
192         }
193
194         public void reComputeErrorAround() {
195             reComputeError();
196             if (nearest_in_other_mesh != null) nearest_in_other_mesh.reComputeError();
197             E e = this.e;
198             do {
199                 e.p2.reComputeError();
200                 e = e.pair.next;
201             } while (e != this.e);
202         }
203         public void reComputeError() {
204             unComputeError();
205             computeError();
206         }
207         public void unComputeError() {
208             score -= oldscore;
209             oldscore = 0;
210         }
211         public void computeError() {
212             if (quadric_count == 0) {
213                 if (!tilemesh) {
214                 }
215                 else if (nearest_in_other_mesh == null) {
216                     if (score_against != null) {
217                         Vert ne = score_against.nearest(p);
218                         oldscore = ne.fundamentalQuadric().preAndPostMultiply(p) * 100 * 10;
219                     } else {
220                         oldscore = 0;
221                     }
222                 } else {
223                     oldscore = nearest_in_other_mesh.fundamentalQuadric().preAndPostMultiply(p) * 100 * 10;
224                 }
225             } else {
226                 oldscore = (quadric.preAndPostMultiply(p) * 100) / quadric_count;
227             }
228
229             oldscore = oldscore;
230
231             int numaspects = 0;
232             float aspects = 0;
233             E e = this.e;
234             do {
235                 //double ang = Math.abs(e.crossAngle());
236                 double ang = Math.abs(e.crossAngle());
237                 if (ang > Math.PI) throw new Error();
238                 /*
239                 if (e.t != null) {
240                     numaspects++;
241                     aspects += e.t.aspect()*e.t.aspect();
242                 }
243                 */
244
245                 float minangle = (float)(Math.PI * 0.8);
246                 if (ang > minangle)
247                     oldscore += (ang - minangle);
248
249                 e = e.pair.next;
250             } while (e != this.e);
251             if (numaspects > 0) oldscore += (aspects / numaspects);
252
253             //System.out.println(oldscore);
254             //oldscore = oldscore*oldscore;
255             score += oldscore;
256         }
257
258         private void removeTrianglesFromRTree() {
259             E e = this.e;
260             do {
261                 if (e.t != null) e.t.removeFromRTree();
262                 e = e.pair.next;
263             } while(e != this.e);
264         }
265         private void addTrianglesToRTree() {
266             E e = this.e;
267             do {
268                 if (e.t != null) e.t.addToRTree();
269                 e = e.pair.next;
270             } while(e != this.e);
271         }
272
273         /** does NOT update bound pairs! */
274         public boolean transform(Matrix m) {
275             unApplyQuadricToNeighbor();
276             Point oldp = this.p;
277             try {
278                 if (vertices.get(this.p)==null) throw new Error();
279                 vertices.remove(this);
280                 removeTrianglesFromRTree();
281                 float newx = m.a*p.x + m.b*p.y + m.c*p.z + m.d;
282                 float newy = m.e*p.x + m.f*p.y + m.g*p.z + m.h;
283                 float newz = m.i*p.x + m.j*p.y + m.k*p.z + m.l;
284                 this.p = new Point(newx, newy, newz);
285                 addTrianglesToRTree();
286                 vertices.add(this);
287             } catch (Exception e) {
288                 throw new RuntimeException(e);
289             }
290             applyQuadricToNeighbor();
291
292             // FIXME: intersection test needed?
293             good = true;
294
295             // should recompute fundamental quadrics of all vertices sharing a face, but we defer...
296             E e = this.e;
297             do {
298                 /*
299                 if (Math.abs(e.crossAngle()) > (Math.PI * 0.9) ||
300                     Math.abs(e.next.crossAngle()) > (Math.PI * 0.9)) {
301                     good = false;
302                 }
303                 if (e.t.aspect() < 0.1) {
304                     good = false;
305                 }
306                 */
307                 e.p2.quadricStale = true;
308                 e = e.pair.next;
309             } while(e != this.e);
310
311
312             if (!ignorecollision && good) {
313
314                 tris.range(new Segment(oldp, this.p),
315                             new Visitor<T>() {
316                                 public void visit(T t) {
317                                     if (!good) return;
318                                     E e = Vert.this.e;
319                                     do {
320                                         if (!t.has(e.p1) && !t.has(e.p2) && e.intersects(t)) { good = false; }
321                                         if (e.t != null) {
322                                             if (!e.t.has(t.e1().p1) && !e.t.has(t.e1().p2) && t.e1().intersects(e.t)) { good = false; }
323                                             if (!e.t.has(t.e2().p1) && !e.t.has(t.e2().p2) && t.e2().intersects(e.t)) { good = false; }
324                                             if (!e.t.has(t.e3().p1) && !e.t.has(t.e3().p2) && t.e3().intersects(e.t)) { good = false; }
325                                         }
326                                         e = e.pair.next;
327                                     } while(e != Vert.this.e);
328                                 }
329                             });
330
331                 /*
332                 for(T t : Mesh.this) {
333                     if (!good) break;
334                     e = this.e;
335                     do {
336                         if (!t.has(e.p1) && !t.has(e.p2) && e.intersects(t)) { good = false; break; }
337                         if (e.t != null) {
338                             if (!e.t.has(t.e1().p1) && !e.t.has(t.e1().p2) && t.e1().intersects(e.t)) { good = false; break; }
339                             if (!e.t.has(t.e2().p1) && !e.t.has(t.e2().p2) && t.e2().intersects(e.t)) { good = false; break; }
340                             if (!e.t.has(t.e3().p1) && !e.t.has(t.e3().p2) && t.e3().intersects(e.t)) { good = false; break; }
341                         }
342                         e = e.pair.next;
343                     } while(e != this.e);
344                 }
345                 */
346             }
347
348
349             reComputeErrorAround();
350             return good;
351         }
352         private boolean good;
353
354         public boolean move(Vec v) {
355             Matrix m = Matrix.translate(v);
356             Vert p = this;
357             boolean good = true;
358             do {
359                 good &= p.transform(m);
360                 p = p.bound_to;
361             } while (p != this);
362             return good;
363         }
364
365         public E getFreeIncident() {
366             E ret = getFreeIncident(e, e);
367             if (ret != null) return ret;
368             ret = getFreeIncident(e.pair.next, e.pair.next);
369             if (ret == null) {
370                 E ex = e;
371                 do {
372                     System.out.println(ex + " " + ex.t);
373                     ex = ex.pair.next;
374                 } while (ex != e);
375                 throw new Error("unable to find free incident to " + this);
376             }
377             return ret;
378         }
379
380         public E getFreeIncident(E start, E before) {
381             E e = start;
382             do {
383                 if (e.pair.p2 == this && e.pair.t == null && e.pair.next.t == null) return e.pair;
384                 e = e.pair.next;
385             } while(e != before);
386             return null;
387         }
388
389         public E getE(Point p2) {
390             Vert v = vertices.get(p2);
391             if (v==null) return null;
392             return getE(v);
393         }
394         public E getE(Vert p2) {
395             E e = this.e;
396             do {
397                 if (e==null) return null;
398                 if (e.p1 == this && e.p2 == p2) return e;
399                 e = e.pair.next;
400             } while (e!=this.e);
401             return null;
402         }
403
404         public Vec norm() {
405             Vec norm = new Vec(0, 0, 0);
406             E e = this.e;
407             do {
408                 if (e.t != null) norm = norm.plus(e.t.norm().times((float)e.prev.angle()));
409                 e = e.pair.next;
410             } while(e != this.e);
411             return norm.norm();
412         }
413
414         public boolean isBoundTo(Vert p) {
415             Vert px = p;
416             do {
417                 if (px==this) return true;
418                 px = px.bound_to;
419             } while(px != p);
420             return false;
421         }
422         public void unbind() { bound_to = this; binding = Matrix.ONE; }
423         public void bind(Vert p) { bind(p, Matrix.ONE); }
424         public void bind(Vert p, Matrix binding) {
425             if (isBoundTo(p)) return;
426             Vert temp_bound_to = p.bound_to;
427             Matrix temp_binding = p.binding;
428             p.bound_to = this.bound_to;
429             p.binding = binding.times(this.binding); // FIXME: may have order wrong here
430             this.bound_to = temp_bound_to;
431             this.binding = temp_binding.times(temp_binding); // FIXME: may have order wrong here
432         }
433     }
434
435     public class BindingGroup {
436         private HashSet<E> set = new HashSet<E>();
437         public BindingGroup bind_others;
438         public BindingGroup other() { return bind_others; }
439         public BindingGroup(BindingGroup bind_others) { this.bind_others = bind_others; }
440         public BindingGroup() { this.bind_others = new BindingGroup(this); }
441         public BindingGroup(E e) { this(); set.add(e); }
442         public void add(E e) {
443             if (set.contains(e)) return;
444             set.add(e);
445             BindingGroup e_bind_peers = e.bind_peers;
446             BindingGroup e_bind_to    = e.bind_to;
447             e.bind_peers = this;
448             e.bind_to    = bind_others;
449             for (E epeer  : e_bind_peers.set) add(epeer);
450             for (E eother : e_bind_to.set)    bind_others.add(eother);
451
452             for(E eother : bind_others.set) {
453                 if (e.next.bind_to.set.contains(eother.prev)) {
454                     e.next.next.bindEdge(eother.prev.prev);
455                 }
456                 if (e.prev.bind_to.set.contains(eother.next)) {
457                     e.prev.prev.bindEdge(eother.next.next);
458                 }
459             }
460
461         }
462         public void dobind(E e) {
463             for(E ebound : set) {
464                 e.p1.bind(ebound.p2);
465                 e.p2.bind(ebound.p1);
466             }
467         }
468         public void shatter(BindingGroup bg1, BindingGroup bg2) {
469             for(E e : set) {
470                 e.shatter(e.midpoint(), bg1, bg2);
471             }
472         }
473     }
474
475     /** [UNIQUE] an edge */
476     public final class E implements Comparable<E> {
477
478         public final Vert p1, p2;
479         T t;     // triangle to our "left"
480         E prev;  // previous half-edge
481         E next;  // next half-edge
482         E pair;  // partner half-edge
483         public BindingGroup bind_peers  = new BindingGroup(this);
484         public BindingGroup bind_to     = bind_peers.other();
485         boolean shattered = false;
486
487         public boolean intersects(T t) { return t.intersects(p1.p, p2.p); }
488         public float comparator() {
489             Vert nearest = score_against.nearest(midpoint());
490             //if (t==null) return length();
491             /*
492             double ang = Math.abs(crossAngle());
493             float minangle = (float)(Math.PI * 0.9);
494             if (ang > minangle)
495                 return 300;
496             */
497             /*
498             if ((length() * length()) / t.area() > 10)
499                 return (float)(length()*Math.sqrt(t.area()));
500             return length()*t.area();
501             */
502             return (float)Math.max(length(), midpoint().distance(nearest.p));
503             //return length();
504         }
505         public int compareTo(E e) {
506             return e.comparator() > comparator() ? 1 : -1;
507         }
508         public void bindEdge(E e) { bind_to.add(e); }
509         public void dobind() { bind_to.dobind(this); }
510
511         public Point shatter() { return shatter(midpoint(), null, null); }
512         public Point shatter(Point mid, BindingGroup bg1, BindingGroup bg2) {
513             if (shattered || destroyed) return mid;
514             shattered = true;
515
516             Vert r = next.p2;
517             E next = this.next;
518             E prev = this.prev;
519
520             int old_colorclass = t==null ? 0 : t.colorclass;
521             if (bg1==null) bg1 = new BindingGroup();
522             if (bg2==null) bg2 = new BindingGroup();
523             BindingGroup old_bind_to = bind_to;
524             bind_peers.shatter(bg1, bg2);
525             old_bind_to.shatter(bg2.other(), bg1.other());
526             pair.shatter();
527             destroy();
528
529             newT(r.p, p1.p, mid, null, old_colorclass);
530             newT(r.p, mid, p2.p, null, old_colorclass);
531             bg1.add(p1.getE(mid));
532             bg2.add(p2.getE(mid).pair);
533             return mid;
534         }
535
536         public boolean destroyed = false;
537         public void destroy() {
538             if (destroyed) return;
539             destroyed = true;
540             pair.destroyed = true;
541
542             if (t != null) t.destroy();
543             t = null;
544
545             if (pair.t != null) pair.t.destroy();
546             pair.t = null;
547
548             if (next.t != null) next.t.destroy();
549             if (prev.t != null) prev.t.destroy();
550             next.t = null;
551             prev.t = null;
552
553             if (pair.next.t != null) pair.next.t.destroy();
554             if (pair.prev.t != null) pair.next.t.destroy();
555             pair.next.t = null;
556             pair.prev.t = null;
557
558             this.bind_to = null;
559             pair.bind_to = null;
560             this.bind_peers = null;
561             pair.bind_peers = null;
562             pair.prev.next = next;
563             next.prev = pair.prev;
564             prev.next = pair.next;
565             pair.next = prev;
566             if (p1.e == this) p1.e = prev.next;
567             if (pair.p1.e == pair) pair.p1.e = pair.prev.next;
568         }
569
570         private void sync() {
571             this.prev.next = this;
572             this.next.prev = this;
573             this.pair.pair = this;
574             bind_peers.add(this);
575             if (this.next.p1 != p2) throw new Error();
576             if (this.prev.p2 != p1) throw new Error();
577             if (this.p1.e == null) this.p1.e = this;
578             if (!added) added = true;
579         }
580         private boolean added = false;
581
582         public T makeT(int colorclass) { return t==null ? (t = new T(this, colorclass)) : t; }
583
584         public double crossAngle() {
585             Vec v1 = t.norm().times(-1);
586             Vec v2 = pair.t.norm().times(-1);
587             return Math.acos(v1.norm().dot(v2.norm()));
588         }
589
590         /** angle between this half-edge and the next */
591         public double angle() {
592             Vec v1 = next.p2.p.minus(p2.p);
593             Vec v2 = this.p1.p.minus(p2.p);
594             return Math.acos(v1.norm().dot(v2.norm()));
595         }
596
597         public void makeAdjacent(E e) {
598             if (this.next == e) return;
599             if (p2 != e.p1) throw new Error("cannot make adjacent -- no shared vertex");
600             if (t != null || e.t != null) throw new Error("cannot make adjacent -- edges not both free");
601
602             E freeIncident = p2.getFreeIncident(e, this);
603
604             e.prev.next = freeIncident.next;
605             freeIncident.next.prev = e.prev;
606
607             freeIncident.next = this.next;
608             this.next.prev = freeIncident;
609             
610             this.next = e;
611             e.prev = this;
612
613             sync();
614             freeIncident.sync();
615         }
616
617         /** creates an isolated edge out in the middle of space */
618         public E(Point p1, Point p2) {
619             if (vertices.get(p1) != null) throw new Error();
620             if (vertices.get(p2) != null) throw new Error();
621             this.p1 = new Vert(p1);
622             this.p2 = new Vert(p2);
623             this.prev = this.next = this.pair = new E(this, this, this);
624             this.p1.e = this;
625             this.p2.e = this.pair;
626             sync();
627         }
628
629         /** adds a new half-edge from prev.p2 to p2 */
630         public E(E prev, Point p) {
631             Vert p2;
632             p2 = vertices.get(p);
633             if (p2 == null) p2 = new Vert(p);
634             this.p1 = prev.p2;
635             this.p2 = p2;
636             this.prev = prev;
637             if (p2.getE(p1) != null) throw new Error();
638             if (p2.e==null) {
639                 this.next = this.pair = new E(this, this, prev.next);
640             } else {
641                 E q = p2.getFreeIncident();
642                 this.next = q.next;
643                 this.next.prev = this;
644                 E z = prev.next;
645                 this.prev.next = this;
646                 this.pair = new E(q, this, z);
647             }
648             if (p2.e==null) p2.e = this.pair;
649             sync();
650         }
651
652         /** adds a new half-edge to the mesh with a given predecessor, successor, and pair */
653         public E(E prev, E pair, E next) {
654             this.p1 = prev.p2;
655             this.p2 = next.p1;
656             this.prev = prev;
657             this.next = next;
658             this.pair = pair;
659             sync();
660         }
661         public Point midpoint() { return new Point((p1.p.x+p2.p.x)/2, (p1.p.y+p2.p.y)/2, (p1.p.z+p2.p.z)/2); }
662         public boolean has(Vert v) { return v==p1 || v==p2; }
663         public float length() { return p1.p.minus(p2.p).mag(); }
664         public String toString() { return p1+"->"+p2; }
665
666     }
667
668     public E makeE(Point p1, Point p2) {
669         Vert v1 = vertices.get(p1);
670         Vert v2 = vertices.get(p2);
671         if (v1 != null && v2 != null) {
672             E e = v1.getE(v2);
673             if (e != null) return e;
674             e = v2.getE(v1);
675             if (e != null) return e;
676         }
677         if (v1 != null) return new E(v1.getFreeIncident(), p2);
678         if (v2 != null) return new E(v2.getFreeIncident(), p1).pair;
679         return new E(p1, p2);
680     }
681     public T newT(Point p1, Point p2, Point p3, Vec norm, int colorclass) {
682         if (norm != null) {
683             Vec norm2 = p3.minus(p1).cross(p2.minus(p1));
684             float dot = norm.dot(norm2);
685             //if (Math.abs(dot) < EPointSILON) throw new Error("dot products within evertsilon of each other: "+norm+" "+norm2);
686             if (dot < 0) { Point p = p1; p1=p2; p2 = p; }
687         }
688         E e12 = makeE(p1, p2);
689         E e23 = makeE(p2, p3);
690         E e31 = makeE(p3, p1);
691         while(e12.next != e23 || e23.next != e31 || e31.next != e12) {
692             e12.makeAdjacent(e23);
693             e23.makeAdjacent(e31);
694             e31.makeAdjacent(e12);
695         }
696         T ret = e12.makeT(colorclass);
697         if (e12.t == null) throw new Error();
698         if (e23.t == null) throw new Error();
699         if (e31.t == null) throw new Error();
700         return ret;
701     }
702
703
704     /** [UNIQUE] a triangle (face) */
705     public final class T extends Triangle {
706         public final E e1;
707         public final int color;
708         public final int colorclass;
709
710         public void removeFromRTree() { tris.remove(this); }
711         public void addToRTree() { tris.insert(this); }
712
713         public void destroy() { tris.remove(this); }
714
715         T(E e1, int colorclass) {
716             this.e1 = e1;
717             E e2 = e1.next;
718             E e3 = e2.next;
719             if (e1==e2 || e1==e3) throw new Error();
720             if (e3.next!=e1) throw new Error();
721             if (e1.t!=null || e2.t!=null || e3.t!=null) throw new Error("non-manifold surface or disagreeing normals");
722             e1.t = this;
723             e1.next.t = this;
724             e1.next.next.t = this;
725
726             // FIXME: check for sealed/watertight surface once construction is complete (and infer normal(s)?)
727
728             int color = Math.abs(random.nextInt());
729             while(true) {
730                 color = color % 4;
731                 if (e1().pair.t != null && color == e1().pair.t.color) { color++; continue; }
732                 if (e2().pair.t != null && color == e2().pair.t.color) { color++; continue; }
733                 if (e3().pair.t != null && color == e3().pair.t.color) { color++; continue; }
734                 break;
735             }
736             this.color = color;
737             this.colorclass = colorclass;
738             tris.add(this);
739         }
740         public E e1() { return e1; }
741         public E e2() { return e1.next; }
742         public E e3() { return e1.prev; }
743         public Vert v1() { return e1.p1; }
744         public Vert v2() { return e1.p2; }
745         public Vert v3() { return e1.next.p2; }
746         public Point p1() { return e1.p1.p; }
747         public Point p2() { return e1.p2.p; }
748         public Point p3() { return e1.next.p2.p; }
749         public boolean hasE(E e) { return e1==e || e1.next==e || e1.prev==e; }
750         public boolean has(Vert v) { return v1()==v || v2()==v || v3()==v; }
751
752         public boolean shouldBeDrawn() {
753             if (e1().bind_to.set.size() == 0) return false;
754             if (e2().bind_to.set.size() == 0) return false;
755             if (e3().bind_to.set.size() == 0) return false;
756             return true;
757         }
758
759     }
760 }