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.berkeley.qfat.geom.HasBindingGroup;
  10 import edu.wlu.cs.levy.CG.KDTree;
  11 import edu.berkeley.qfat.geom.Point;
  12 import com.infomatiq.jsi.IntProcedure;
  13
  14 public class Mesh implements Iterable<Mesh.T> {
  15
  16     public static final float EPSILON = (float)0.0001;
  17     public static final Random random = new Random();
  18
  19     private RTree<T>         triangles = new RTree<T>();
  20     private PointSet<Vertex> vertices  = new PointSet<Vertex>();
  21
  22     public boolean option_wireframe    = false;
  23     public boolean option_errorNormals = false;
  24     public boolean option_selectable   = true;
  25
  26     public void render(GL gl, Matrix m) {
  27         if (option_wireframe) {
  28             gl.glDisable(GL.GL_LIGHTING);
  29             gl.glBegin(GL.GL_LINES);
  30             gl.glColor3f(1, 1, 1);
  31             for (T t : this) {
  32                 m.times(t.e1().p1.goodp).glVertex(gl);
  33                 m.times(t.e1().p2.goodp).glVertex(gl);
  34                 m.times(t.e2().p1.goodp).glVertex(gl);
  35                 m.times(t.e2().p2.goodp).glVertex(gl);
  36                 m.times(t.e3().p1.goodp).glVertex(gl);
  37                 m.times(t.e3().p2.goodp).glVertex(gl);
  38             }
  39             gl.glEnd();
  40             gl.glEnable(GL.GL_LIGHTING);
  41             return;
  42         }
  43         for(T t : this) {
  44             gl.glColor4f((float)(0.25+(0.05*t.color)),
  45                          (float)(0.25+(0.05*t.color)),
  46                          (float)(0.75+(0.05*t.color)),
  47                          (float)0.3);
  48             t.glTriangle(gl, m);
  49         }
  50         if (option_errorNormals)
  51             for(T t : this)
  52                 for(Mesh.Vertex p : new Mesh.Vertex[] { t.v1(), t.v2(), t.v3() }) {
  53                     if (p.ok) {
  54                         gl.glBegin(GL.GL_LINES);
  55                         gl.glColor3f(1, 1, 1);
  56                         p.p.glVertex(gl);
  57                         p.p.plus(p.norm().times((float)p.error()*10)).glVertex(gl);
  58                         gl.glEnd();
  59                     }
  60                 }
  61     }
  62
  63     public boolean immutableVertices;
  64     public Mesh    error_against      = null;
  65     public double  error              = 0;
  66
  67     public Mesh(boolean immutableVertices) { this.immutableVertices = immutableVertices; }
  68
  69     public void makeVerticesImmutable() { this.immutableVertices = true; }
  70     public float error() { return (float)error; }
  71
  72     public int size() { return vertices.size(); }
  73     public Iterable<Vertex> vertices() { return vertices; }
  74     public Iterator<T> iterator() { return triangles.iterator(); }
  75
  76     public void rebindPoints() {
  77         // unbind all points
  78         for(Mesh.T t : this) {
  79             t.v1().unbind();
  80             t.v2().unbind();
  81             t.v3().unbind();
  82         }
  83         // ask edges to re-implement their bindings
  84         for(Mesh.T t : this) {
  85             t.e1().dobind();
  86             t.e2().dobind();
  87             t.e3().dobind();
  88         }
  89         System.out.println("rebound!");
  90     }
  91
  92     public void transform(Matrix m) {
  93         ArrayList<Vertex> set = new ArrayList<Vertex>();
  94         for(Vertex v : vertices) set.add(v);
  95         for(Vertex v : set) v.transform(m.times(v.p), true, null);
  96         for(Vertex v : set) v.goodp = v.p;
  97     }
  98
  99     public void rebuild() { /*vertices.rebuild();*/ }
 100     public Vec diagonal() { return vertices.diagonal(); }
 101     public Point centroid() { return vertices.centroid(); }
 102     public Vertex nearest(Point p) { return vertices.nearest(p); }
 103
 104     /** compute the volume of the mesh */
 105     public float volume() {
 106         double total = 0;
 107         for(T t : this) {
 108             double area = t.area();
 109             Vec origin_to_centroid = new Vec(new Point(0, 0, 0), t.centroid());
 110             boolean facingAway = t.norm().dot(origin_to_centroid) > 0;
 111             double height = Math.abs(t.norm().dot(origin_to_centroid));
 112             total += ((facingAway ? 1 : -1) * area * height) / 3.0;
 113         }
 114         return (float)total;
 115     }
 116
 117
 118     public void subdivide() {
 119         for (T t : this) t.old = true;
 120         for (Vertex v : vertices()) v.original = true;
 121         OUTER: while(true) {
 122             for (T t : this)
 123                 if (t.old) {
 124                     Point p = t.e1.midpoint();
 125                     System.out.println("shatter " + t.e1);
 126                     t.e1.shatter();
 127                     nearest(p).edge = true;
 128                     continue OUTER;
 129                 }
 130             break;
 131         }
 132     }
 133
 134     // Vertexices //////////////////////////////////////////////////////////////////////////////
 135
 136     /** a vertex in the mesh */
 137     public final class Vertex extends HasQuadric implements Visitor {
 138         public Point p, goodp;
 139         public Point oldp;
 140         E e;                // some edge *leaving* this point
 141
 142         public boolean original = false;
 143         public boolean edge = false;
 144
 145         private boolean illegal = false;
 146
 147         public boolean visible = false;
 148
 149         public Point getPoint() { return p; }
 150         public float error() { return olderror; }
 151
 152         private Vertex(Point p) {
 153             this.p = p;
 154             this.goodp = p;
 155             this.oldp = p;
 156             if (vertices.get(p) != null) throw new Error();
 157             vertices.add(this);
 158         }
 159
 160         public void reinsert() {
 161             vertices.remove(this);
 162             vertices.add(this);
 163             for(E e = this.e; e!=null; e=e.pair.next==this.e?null:e.pair.next) e.t.reinsert();
 164         }
 165
 166         public float olderror = 0;
 167         public void setError(float nerror) {
 168             error -= olderror;
 169             olderror = nerror;
 170             error += olderror;
 171         }
 172
 173         /*
 174         public Vertex hack(GL gl, Point mouse) {
 175             double dist = Double.MAX_VALUE;
 176             Vertex cur = null;
 177             for(E e = this.e; e!=null; e=e.pair.next==this.e?null:e.pair.next) {
 178                 Vertex v = e.getOther(this);
 179                 double dist2 = v.getPoint().glProject(gl).distance(mouse);
 180                 if ((cur==null || dist2 < dist) && v.visible) {
 181                     dist = dist2;
 182                     cur = v;
 183                 }
 184             }
 185             return cur;
 186         }
 187         */
 188
 189         public float averageTriangleArea() {
 190             int count = 0;
 191             float ret = 0;
 192             for(E e = this.e; e!=null; e=e.pair.next==this.e?null:e.pair.next) {
 193                 ret += e.t.area();
 194                 count++;
 195             }
 196             return ret/count;
 197         }
 198         public float averageEdgeLength() {
 199             int count = 0;
 200             float ret = 0;
 201             for(E e = this.e; e!=null; e=e.pair.next==this.e?null:e.pair.next) {
 202                 ret += e.length();
 203                 count++;
 204             }
 205             return ret/count;
 206         }
 207
 208         public Matrix _recomputeFundamentalQuadric() {
 209             Matrix m = Matrix.ZERO;
 210             int count = 0;
 211             for(E e = this.e; e!=null; e=e.pair.next==this.e?null:e.pair.next) {
 212                 m = m.plus(e.t.norm().fundamentalQuadric(e.t.centroid()));
 213                 count++;
 214             }
 215             if (count > 0) {
 216                 m = m.plus(norm().fundamentalQuadric(this.p).times(count));
 217                 count *= 2;
 218             }
 219             return m.times(1/(float)count);
 220         }
 221
 222         public HasQuadric nearest() { return error_against==null ? null : error_against.vertices.nearest(p, this); }
 223         public void computeError() {
 224             if (error_against==null) return;
 225             float nerror =
 226                 nearest_in_other_mesh != null
 227                 ? nearest_in_other_mesh.fundamentalQuadric().preAndPostMultiply(p)
 228                 : nearest().fundamentalQuadric().preAndPostMultiply(p);
 229             if (quadric_count != 0)
 230                 nerror = (nerror + quadric.preAndPostMultiply(p))/(quadric_count+1);
 231
 232             if (!immutableVertices && quadric_count == 0) {
 233                 //nerror = Math.max(nerror, 0.4f);
 234                 //nerror *= 2;
 235             }
 236             //System.out.println(nerror);
 237             for(E e = this.e; e!=null; e=e.pair.next==this.e?null:e.pair.next) {
 238                 double ang = e.dihedralAngle();
 239                 if (ang > Math.PI) throw new Error();
 240                 if (ang < -Math.PI) throw new Error();
 241                 float minangle = (float)(Math.PI * 0.8);
 242                 //nerror += ((ang / Math.PI)*(ang/Math.PI)) * e.length() * 0.05;
 243
 244                 nerror += (1-e.t.quality())*0.0001;
 245                 if (ang > minangle) nerror += (ang - minangle);
 246
 247                 //System.out.println(((ang / Math.PI)*(ang/Math.PI)) * 0.000001);
 248                 /*
 249                 if (e.t.aspect() < 0.2) {
 250                     nerror += (0.2-e.t.aspect()) * 10;
 251                 }
 252                 */
 253             }
 254             if (!immutableVertices) {
 255                 Vertex n = (Vertex)nearest();
 256                 float d = norm().dot(n.norm());
 257                 if (d > 1 || d < -1) throw new Error();
 258                 if (d >= 0) {
 259                     nerror *= (2.0f - d);
 260                 } else {
 261                     nerror += 0.0003 * (2.0f + d);
 262                     nerror *= (2.0f + d);
 263                 }
 264             }
 265
 266             setError(nerror);
 267         }
 268
 269         public boolean move(Vec vv, boolean ignoreProblems) {
 270
 271             boolean good = true;
 272
 273             //     t1' = M * t1
 274             //     t2' = t2.getMatrix(t1) * t1'
 275             //     t2' = t2.getMatrix(t1) * M * t1
 276             //     t1 =     t1.getMatrix(t2) * t2
 277             // M * t1 = M * t1.getMatrix(t2) * t2
 278
 279             if (bindingGroup!=null && this != bindingGroup.getMaster()) {
 280                 Matrix m2 = getBindingMatrix(bindingGroup.getMaster());
 281                 Vec v2 = m2.times(vv.plus(getPoint())).minus(m2.times(getPoint()));
 282                 return ((Vertex)bindingGroup.getMaster()).move(v2, ignoreProblems);
 283             }
 284
 285             Point op = this.p;
 286             Point pp = vv.plus(getPoint());
 287             if (bindingGroup != null) {
 288                 /*
 289                 for(int i=0; i<20 ; i++) {
 290                     Point p2 = getConstraint().times(pp);
 291                     pp = pp.midpoint(p2);
 292                     //System.out.println(m.minus(m2));
 293                 }
 294             */
 295                     pp = getConstraint().times(pp);
 296             }
 297             pp = pp.minus(op).norm().times(vv.mag()).plus(op);
 298             ok = false;
 299             Point pt = pp;
 300             for(Vertex v : (Iterable<Vertex>)getBoundPeers()) {
 301                 Point pt2 = v.getBindingMatrix(this).times(pt);
 302                 /*
 303                 if (Math.abs( v.p.minus(pt2).mag() / pt.minus(op).mag() ) > 5)
 304                     throw new Error(v.p+" "+pt2+"\n"+op+" "+pt+"\n"+v.getBindingMatrix(this));
 305                 if (Math.abs( v.p.minus(pt2).mag() / pt.minus(op).mag() ) < 1/5) throw new Error();
 306                 */
 307                 good &= v.transform(pt2, ignoreProblems, v.getBindingMatrix(this));
 308             }
 309
 310             if (!good && !ignoreProblems) {
 311                 for(Vertex v : (Iterable<Vertex>)getBoundPeers())
 312                     v.transform(v.oldp, true, null);
 313             }
 314
 315             for(Vertex v : (Iterable<Vertex>)getBoundPeers())
 316                 v.recomputeFundamentalQuadricIfNeighborChanged();
 317             for(Vertex v : (Iterable<Vertex>)getBoundPeers())
 318                 v.reComputeErrorAround();
 319             ok = true;
 320             return good;
 321         }
 322         public boolean ok = true;
 323
 324         /** does NOT update bound pairs! */
 325         private boolean transform(Point newp, boolean ignoreProblems, Matrix yes) {
 326             this.oldp = this.p;
 327             if (immutableVertices) throw new Error();
 328
 329             unApplyQuadricToNeighbor();
 330
 331             boolean illegalbefore = illegal;
 332             illegal = false;
 333             /*
 334             if (this.p.minus(newp).mag() > 0.1 && !ignoreProblems) {
 335                 try {
 336                     throw new Exception(""+this.p.minus(newp).mag()+" "+ignoreProblems+" "+yes);
 337                 } catch(Exception e) {
 338                     e.printStackTrace();
 339                 }
 340                 illegal = true;
 341             }
 342             */
 343
 344             this.p = newp;
 345             reinsert();
 346             applyQuadricToNeighbor();
 347
 348             if (!ignoreProblems) {
 349                 checkLegality();
 350             }
 351             for(E e = this.e; e!=null; e=e.pair.next==this.e?null:e.pair.next)
 352                 e.p2.quadricStale = true;
 353             return !illegal || (illegalbefore && illegal);
 354         }
 355
 356         public void checkLegality() {
 357             /*
 358             for(E e = this.e; e!=null; e=e.pair.next==this.e?null:e.pair.next) {
 359                 if (Math.abs(e.dihedralAngle()) > (Math.PI * 0.9) ||
 360                     Math.abs(e.next.dihedralAngle()) > (Math.PI * 0.9)) illegal = true;
 361                 if (e.t.aspect() < 0.2) illegal = true;
 362             }
 363             */
 364             if (!illegal) triangles.range(oldp, this.p, (Visitor<T>)this);
 365         }
 366
 367         public void reComputeErrorAround() {
 368             reComputeError();
 369             if (nearest_in_other_mesh != null)
 370                 nearest_in_other_mesh.reComputeError();
 371             for(E e = this.e; e!=null; e=e.pair.next==this.e?null:e.pair.next)
 372                 e.p2.reComputeError();
 373         }
 374
 375         public boolean visit(Object o) {
 376             if (o instanceof Vertex)
 377                 return ((Vertex)o).e != null && ((Vertex)o).norm().dot(Vertex.this.norm()) >= 0;
 378             T t = (T)o;
 379             if (illegal) return false;
 380             for(E e = Vertex.this.e; e!=null; e=e.pair.next==Vertex.this.e?null:e.pair.next) {
 381                 if (!t.has(e.p1) && !t.has(e.p2) && e.intersects(t)) { illegal = true; }
 382                 if (e.t != null) {
 383                     if (!e.t.has(t.e1().p1) && !e.t.has(t.e1().p2) && t.e1().intersects(e.t)) { illegal = true; }
 384                     if (!e.t.has(t.e2().p1) && !e.t.has(t.e2().p2) && t.e2().intersects(e.t)) { illegal = true; }
 385                     if (!e.t.has(t.e3().p1) && !e.t.has(t.e3().p2) && t.e3().intersects(e.t)) { illegal = true; }
 386                 }
 387             }
 388             return !illegal;
 389         }
 390
 391         public E getEdge() { return e; }
 392         public E getFreeIncident() {
 393             E ret = getFreeIncident(e, e);
 394             if (ret != null) return ret;
 395             for(E e = this.e; e!=null; e=e.pair.next==this.e?null:e.pair.next)
 396                 System.out.println(e + " " + e.t);
 397             throw new Error("unable to find free incident to " + this);
 398         }
 399
 400         public E getFreeIncident(E start, E before) {
 401             for(E e = start; e!=null; e=e.pair.next==before?null:e.pair.next)
 402                 if (e.pair.p2 == this && e.pair.t == null && e.pair.next.t == null)
 403                     return e.pair;
 404             return null;
 405         }
 406
 407         public E getE(Point p2) {
 408             Vertex v = vertices.get(p2);
 409             if (v==null) return null;
 410             return getE(v);
 411         }
 412         public E getE(Vertex p2) {
 413             for(E e = this.e; e!=null; e=e.pair.next==this.e?null:e.pair.next)
 414                 if (e.p1 == this && e.p2 == p2) return e;
 415             return null;
 416         }
 417
 418         private void glNormal(GL gl) {
 419             Vec norm = norm();
 420             gl.glNormal3f(norm.x, norm.y, norm.z);
 421         }
 422         public Vec norm() {
 423             Vec norm = new Vec(0, 0, 0);
 424             for(E e = this.e; e!=null; e=e.pair.next==this.e?null:e.pair.next)
 425                 if (e.t != null)
 426                     norm = norm.plus(e.t.norm().times((float)e.prev.angle()));
 427             return norm.norm();
 428         }
 429
 430         public void bindTo(Vertex p) { bindTo(Matrix.ONE, p); }
 431     }
 432
 433
 434     /** [UNIQUE] an edge */
 435     public final class E extends HasBindingGroup implements Comparable<E> {
 436
 437         public final Vertex p1, p2;
 438         T t;     // triangle to our "left"
 439         E prev;  // previous half-edge
 440         E next;  // next half-edge
 441         E pair;  // partner half-edge
 442         boolean shattered = false;
 443
 444         public boolean intersects(T t) { return t.intersects(p1.p, p2.p); }
 445
 446         public Segment getSegment() { return new Segment(p1.getPoint(), p2.getPoint()); }
 447
 448         public void bindingGroupChanged(edu.berkeley.qfat.geom.BindingGroup newBindingGroup_) {
 449
 450             edu.berkeley.qfat.geom.BindingGroup<E> newBindingGroup =
 451                 (edu.berkeley.qfat.geom.BindingGroup<E>)newBindingGroup_;
 452             if (newBindingGroup==null) return;
 453             if (this==newBindingGroup.getMaster()) return;
 454             HashSet<E> nbg = new HashSet<E>();
 455             for(E eother : (Iterable<E>)newBindingGroup) nbg.add(eother);
 456             for(E eother : nbg) {
 457                 if (next==null || prev==null) continue;
 458                 if (eother.next==null || eother.prev==null) continue;
 459
 460                 if (next.isBoundTo(eother.pair.prev.pair) && !prev.isBoundTo(eother.pair.next.pair))
 461                     prev.bindTo(next.getBindingMatrix(eother.pair.prev.pair), eother.pair.next.pair);
 462                 if (!next.isBoundTo(eother.pair.prev.pair) && prev.isBoundTo(eother.pair.next.pair))
 463                     next.bindTo(prev.getBindingMatrix(eother.pair.next.pair), eother.pair.prev.pair);
 464
 465                 /*
 466                 if (next.isBoundTo(eother.prev) && !prev.isBoundTo(eother.next))
 467                     prev.bindTo(next.getBindingMatrix(eother.prev), eother.next);
 468                 if (!next.isBoundTo(eother.prev) && prev.isBoundTo(eother.next))
 469                     next.bindTo(prev.getBindingMatrix(eother.next), eother.prev);
 470                 */
 471                 if (next.isBoundTo(eother.next) && !prev.isBoundTo(eother.prev))
 472                     prev.bindTo(next.getBindingMatrix(eother.next), eother.prev);
 473                 if (!next.isBoundTo(eother.next) && prev.isBoundTo(eother.prev))
 474                     next.bindTo(prev.getBindingMatrix(eother.prev), eother.next);
 475             }
 476
 477         }
 478
 479         public float stretchRatio() {
 480             Vertex nearest = error_against.nearest(midpoint());
 481             float nearest_distance = midpoint().distance(nearest.p);
 482             float other_distance =
 483                 (p1.p.distance(error_against.nearest(p1.p).p)+
 484                  p2.p.distance(error_against.nearest(p2.p).p))/2;
 485             return nearest_distance/other_distance;
 486         }
 487         public float comparator() {
 488             return length();
 489         }
 490         public int compareTo(E e) {
 491             return e.comparator() > comparator() ? 1 : -1;
 492         }
 493         public void bindEdge(E e, Matrix m) {
 494             _bindEdge(e, m);
 495             pair._bindEdge(e.pair, m);
 496         }
 497         public void _bindEdge(E e, Matrix m) {
 498             e = e.pair;
 499             /*
 500             //assumes edges are identical length at binding time
 501             Vec reflectionPlaneNormal = e.p2.p.minus(e.p1.p).norm();
 502             float a = reflectionPlaneNormal.x;
 503             float b = reflectionPlaneNormal.y;
 504             float c = reflectionPlaneNormal.z;
 505             Matrix reflectionMatrix =
 506                 new Matrix( 1-2*a*a,  -2*a*b,  -2*a*c, 0,
 507                             -2*a*b,  1-2*b*b,  -2*b*c, 0,
 508                             -2*a*c,   -2*b*c, 1-2*c*c, 0,
 509                             0,       0,       0,       1);
 510             m = m.times(Matrix.translate(e.midpoint().minus(Point.ORIGIN))
 511                         .times(reflectionMatrix)
 512                         .times(Matrix.translate(Point.ORIGIN.minus(e.midpoint()))));
 513             System.out.println(reflectionPlaneNormal);
 514             System.out.println("  " + p1.p + " " + m.times(e.p1.p));
 515             System.out.println("  " + p2.p + " " + m.times(e.p2.p));
 516             */
 517             /*
 518             if (m.times(e.p1.p).minus(p1.p).mag() > EPSILON) throw new Error();
 519             if (m.times(e.p2.p).minus(p2.p).mag() > EPSILON) throw new Error();
 520             */
 521             this.bindTo(m, e);
 522         }
 523
 524         public void dobind() {
 525             for(E e : (Iterable<E>)getBoundPeers()) {
 526                 if (e==this) continue;
 527                 p1.bindTo(getBindingMatrix(e), e.p1);
 528                 p2.bindTo(getBindingMatrix(e), e.p2);
 529                 e.p1.setConstraint(getConstraint());
 530                 e.p2.setConstraint(getConstraint());
 531             }
 532         }
 533
 534         public Point shatter() {
 535             if (shattered || destroyed) return null;
 536             shattered = true;
 537             E first = null;
 538             E firste = null;
 539             E firstx = null;
 540             E firstq = null;
 541             for(E e : (Iterable<E>)getBoundPeers()) {
 542                 E enext = e.next;
 543                 E eprev = e.prev;
 544                 E pnext = e.pair.next;
 545                 E pprev = e.pair.prev;
 546                 Point mid = e.midpoint();
 547                 Vertex r = e.next.p2;
 548                 Vertex l = e.pair.next.p2;
 549                 if (!e.destroyed) {
 550                     e.destroy();
 551                     e.pair.destroy();
 552                     newT(r.p, e.p1.p, mid,    null, 0);
 553                     newT(r.p, mid,    e.p2.p, null, 0);
 554                     newT(l.p, mid,    e.p1.p, null, 0);
 555                     newT(l.p, e.p2.p, mid,    null, 0);
 556                 }
 557             }
 558             for(E e : (Iterable<E>)getBoundPeers()) {
 559                 Point mid = e.midpoint();
 560                 if (first==null) {
 561                     first = e.p1.getE(mid);
 562                     firste = e;
 563                     firstx = e.pair;
 564                     firstq = e.p2.getE(mid).pair;
 565                     continue;
 566                 }
 567                 e.p1.getE(mid).          bindTo(e.getBindingMatrix(firste), first);
 568                 e.p1.getE(mid).pair.     bindTo(e.getBindingMatrix(firste), first.pair);
 569                 e.p2.getE(mid).pair.     bindTo(e.getBindingMatrix(firste), firstq);
 570                 e.p2.getE(mid).pair.pair.bindTo(e.getBindingMatrix(firste), firstq.pair);
 571             }
 572             /*
 573             first.setConstraint(firste.getConstraint());
 574             firstq.setConstraint(firste.getConstraint());
 575             */
 576             return null;
 577         }
 578
 579         public boolean destroyed = false;
 580         public void destroy() {
 581             if (destroyed) return;
 582             destroyed = true;
 583             pair.destroyed = true;
 584
 585             if (t != null) t.destroy();
 586             t = null;
 587
 588             if (pair.t != null) pair.t.destroy();
 589             pair.t = null;
 590
 591             if (next.t != null) next.t.destroy();
 592             if (prev.t != null) prev.t.destroy();
 593             next.t = null;
 594             prev.t = null;
 595
 596             if (pair.next.t != null) pair.next.t.destroy();
 597             if (pair.prev.t != null) pair.next.t.destroy();
 598             pair.next.t = null;
 599             pair.prev.t = null;
 600
 601             pair.prev.next = next;
 602             next.prev = pair.prev;
 603             prev.next = pair.next;
 604             pair.next = prev;
 605             if (p1.e == this) p1.e = prev.next;
 606             if (pair.p1.e == pair) pair.p1.e = pair.prev.next;
 607         }
 608
 609         private void sync() {
 610             this.prev.next = this;
 611             this.next.prev = this;
 612             this.pair.pair = this;
 613             if (this.next.p1 != p2) throw new Error();
 614             if (this.prev.p2 != p1) throw new Error();
 615             if (this.p1.e == null) this.p1.e = this;
 616             if (!added) added = true;
 617         }
 618         private boolean added = false;
 619
 620         public T makeT(int colorclass) { return t==null ? (t = new T(this, colorclass)) : t; }
 621
 622         public double dihedralAngle() {
 623             Vec v1 = t.norm().times(-1);
 624             Vec v2 = pair.t.norm().times(-1);
 625             double prod = v1.norm().dot(v2.norm());
 626             prod = Math.min(1,prod);
 627             prod = Math.max(-1,prod);
 628             double ret = Math.acos(prod);
 629             if (Double.isNaN(ret)) throw new Error("nan! " + prod);
 630             return ret;
 631         }
 632
 633         /** angle between this half-edge and the next */
 634         public double angle() {
 635             Vec v1 = next.p2.p.minus(p2.p);
 636             Vec v2 = this.p1.p.minus(p2.p);
 637             return Math.acos(v1.norm().dot(v2.norm()));
 638         }
 639
 640         public Vertex getOther(Vertex v) {
 641             if (this.p1 == v) return p2;
 642             if (this.p2 == v) return p1;
 643             throw new Error();
 644         }
 645
 646         public void makeAdjacent(E e) {
 647             if (this.next == e) return;
 648             if (p2 != e.p1) throw new Error("cannot make adjacent -- no shared vertex");
 649             if (t != null || e.t != null) throw new Error("cannot make adjacent -- edges not both free ");
 650
 651             E freeIncident = p2.getFreeIncident(e, this);
 652
 653             e.prev.next = freeIncident.next;
 654             freeIncident.next.prev = e.prev;
 655
 656             freeIncident.next = this.next;
 657             this.next.prev = freeIncident;
 658
 659             this.next = e;
 660             e.prev = this;
 661
 662             sync();
 663             freeIncident.sync();
 664         }
 665
 666         /** creates an isolated edge out in the middle of space */
 667         public E(Point p1, Point p2) {
 668             if (vertices.get(p1) != null) throw new Error();
 669             if (vertices.get(p2) != null) throw new Error();
 670             this.p1 = new Vertex(p1);
 671             this.p2 = new Vertex(p2);
 672             this.prev = this.next = this.pair = new E(this, this, this);
 673             this.p1.e = this;
 674             this.p2.e = this.pair;
 675             sync();
 676         }
 677
 678         /** adds a new half-edge from prev.p2 to p2 */
 679         public E(E prev, Point p) {
 680             Vertex p2;
 681             p2 = vertices.get(p);
 682             if (p2 == null) p2 = new Vertex(p);
 683             this.p1 = prev.p2;
 684             this.p2 = p2;
 685             this.prev = prev;
 686             if (p2.getE(p1) != null) throw new Error();
 687             if (p2.e==null) {
 688                 this.next = this.pair = new E(this, this, prev.next);
 689             } else {
 690                 E q = p2.getFreeIncident();
 691                 this.next = q.next;
 692                 this.next.prev = this;
 693                 E z = prev.next;
 694                 this.prev.next = this;
 695                 this.pair = new E(q, this, z);
 696             }
 697             if (p2.e==null) p2.e = this.pair;
 698             sync();
 699         }
 700
 701         /** adds a new half-edge to the mesh with a given predecessor, successor, and pair */
 702         public E(E prev, E pair, E next) {
 703             this.p1 = prev.p2;
 704             this.p2 = next.p1;
 705             this.prev = prev;
 706             this.next = next;
 707             this.pair = pair;
 708             sync();
 709         }
 710         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); }
 711         public boolean has(Vertex v) { return v==p1 || v==p2; }
 712         public float length() { return p1.p.minus(p2.p).mag(); }
 713         public String toString() { return p1+"->"+p2; }
 714
 715     }
 716
 717     public E makeE(Point p1, Point p2) {
 718         Vertex v1 = vertices.get(p1);
 719         Vertex v2 = vertices.get(p2);
 720         if (v1 != null && v2 != null) {
 721             E e = v1.getE(v2);
 722             if (e != null) return e;
 723             e = v2.getE(v1);
 724             if (e != null) return e;
 725         }
 726         if (v1 != null) return new E(v1.getFreeIncident(), p2);
 727         if (v2 != null) return new E(v2.getFreeIncident(), p1).pair;
 728         return new E(p1, p2);
 729     }
 730     public boolean coalesce = false;
 731     private static float round(float f) {
 732         return Math.round(f*1000)/1000f;
 733     }
 734     public T newT(Point p1, Point p2, Point p3, Vec norm, int colorclass) {
 735         if (coalesce) {
 736
 737             for(Vertex v : vertices) { if (p1.distance(v.p) < EPSILON) { p1 = v.p; break; } }
 738             for(Vertex v : vertices) { if (p2.distance(v.p) < EPSILON) { p2 = v.p; break; } }
 739             for(Vertex v : vertices) { if (p3.distance(v.p) < EPSILON) { p3 = v.p; break; } }
 740             /*
 741             p1 = new Point(round(p1.x), round(p1.y), round(p1.z));
 742             p2 = new Point(round(p2.x), round(p2.y), round(p2.z));
 743             p3 = new Point(round(p3.x), round(p3.y), round(p3.z));
 744             */
 745         }
 746         if (norm != null) {
 747             Vec norm2 = p3.minus(p1).cross(p2.minus(p1));
 748             float dot = norm.dot(norm2);
 749             //if (Math.abs(dot) < EPointSILON) throw new Error("dot products within evertsilon of each other: "+norm+" "+norm2);
 750             if (dot < 0) { Point p = p1; p1=p2; p2 = p; }
 751         }
 752         E e12 = makeE(p1, p2);
 753         E e23 = makeE(p2, p3);
 754         E e31 = makeE(p3, p1);
 755         while(e12.next != e23 || e23.next != e31 || e31.next != e12) {
 756             e12.makeAdjacent(e23);
 757             e23.makeAdjacent(e31);
 758             e31.makeAdjacent(e12);
 759         }
 760         T ret = e12.makeT(colorclass);
 761         if (e12.t == null) throw new Error();
 762         if (e23.t == null) throw new Error();
 763         if (e31.t == null) throw new Error();
 764         return ret;
 765     }
 766
 767     private int max_serial = 0;
 768     /** [UNIQUE] a triangle (face) */
 769     public final class T extends Triangle {
 770         public final E e1;
 771         public final int color;
 772         public final int colorclass;
 773
 774         public boolean old;
 775
 776         public final int serial = max_serial++;
 777         public boolean occluded;
 778
 779         T(E e1, int colorclass) {
 780             this.e1 = e1;
 781             E e2 = e1.next;
 782             E e3 = e2.next;
 783             if (e1==e2 || e1==e3) throw new Error();
 784             if (e3.next!=e1) throw new Error();
 785             if (e1.t!=null || e2.t!=null || e3.t!=null) throw new Error("non-manifold surface or disagreeing normals");
 786             e1.t = this;
 787             e1.next.t = this;
 788             e1.next.next.t = this;
 789
 790             // FIXME: check for sealed/watertight surface once construction is complete (and infer normal(s)?)
 791
 792             int color = Math.abs(random.nextInt());
 793             while(true) {
 794                 color = color % 4;
 795                 if (e1().pair.t != null && color == e1().pair.t.color) { color++; continue; }
 796                 if (e2().pair.t != null && color == e2().pair.t.color) { color++; continue; }
 797                 if (e3().pair.t != null && color == e3().pair.t.color) { color++; continue; }
 798                 break;
 799             }
 800             this.color = color;
 801             this.colorclass = colorclass;
 802             triangles.add(this);
 803         }
 804         public E e1() { return e1; }
 805         public E e2() { return e1.next; }
 806         public E e3() { return e1.prev; }
 807         public Vertex v1() { return e1.p1; }
 808         public Vertex v2() { return e1.p2; }
 809         public Vertex v3() { return e1.next.p2; }
 810         public Point p1() { return e1.p1.p; }
 811         public Point p2() { return e1.p2.p; }
 812         public Point p3() { return e1.next.p2.p; }
 813         public boolean hasE(E e) { return e1==e || e1.next==e || e1.prev==e; }
 814         public boolean has(Vertex v) { return v1()==v || v2()==v || v3()==v; }
 815
 816         public void removeFromRTree() { triangles.remove(this); }
 817         public void addToRTree() { triangles.insert(this); }
 818         public void destroy() { triangles.remove(this); }
 819         public void reinsert() { triangles.remove(this); triangles.add(this); }
 820
 821         public boolean shouldBeDrawn() {
 822
 823             if (e1().bindingGroupSize() <= 1) return false;
 824             if (e2().bindingGroupSize() <= 1) return false;
 825             if (e3().bindingGroupSize() <= 1) return false;
 826
 827             return true;
 828         }
 829
 830         public void glTriangle(GL gl, Matrix m) {
 831             gl.glPushName(serial);
 832             gl.glBegin(GL.GL_TRIANGLES);
 833             glVertices(gl, m);
 834             gl.glEnd();
 835             gl.glPopName();
 836         }
 837
 838         /** issue gl.glVertex() for each of the triangle's points */
 839         public void glVertices(GL gl, Matrix m) {
 840             if (!shouldBeDrawn()) return;
 841             super.glVertices(gl, m);
 842         }
 843     }
 844 }