questionable patch: merge of a lot of stuff from the svg branch
[org.ibex.core.git] / src / org / ibex / graphics / Mesh.java
index 0680a4c..274f668 100644 (file)
@@ -6,17 +6,19 @@ package org.ibex.graphics;
 import java.util.*;
 import java.util.collections.*;
 import org.ibex.util.*;
+import java.awt.*;
 
 // TODO:
+
+// FIXME:   Not all triangles who get their dirty bit set wind up having fixup() called on them -- bad!!
+
+// FEATURE: Delauanay refinement
+
 //  - allow edge-constraint removal
-//
-//  ~30% of our time spent finding vertices => use a balanced quadtree
-//
 //  - store "which curve is inside me" pointer in Triangle
 //      - split if two curves enter
 //  - go to treating Vertex as a value class (epsilon==0)
 //  - union()
-//  - subtract()
 //  - [??] preserve in/out-ness every time we delete() a triangle
 
 /**
@@ -25,9 +27,11 @@ import org.ibex.util.*;
  */
 public final class Mesh {
 
-    private static final float epsilon = (float)0.0001;
-    private static final float epsilon2 = (float)0.001;
+    private static final float epsilon = (float)0.0;
+    private static final float epsilon2 = (float)0.0;
     private static final boolean debug = false;
+    //private static final boolean check = true;
+    private static final boolean check = false;
 
     private  Vector    triangles   = new Vector();     /* we no longer need this */
     private  Hash      edges       = new Hash();       /* we no longer need this either */
@@ -99,6 +103,7 @@ public final class Mesh {
         iter[numiter++] = triangle0;
         while(numiter > 0) {
             Triangle t = iter[--numiter];
+            if (t==null) return;  // FIXME: is this right?!?
             if (t.tick >= this.tick) continue;
             switch(mode) {
                 case ITERATE_STROKE:          t.stroke(buf, a, color); break;
@@ -107,7 +112,8 @@ public final class Mesh {
                 case ITERATE_INTERSECT:
                 case ITERATE_SUBTRACT: {
                     if (!t.in) break;
-                    boolean oin = m.queryPoint(t.c().multiply(a));
+                    Point p = Imprecise.center(Mesh.this, t.v(1), t.v(2), t.v(3));
+                    boolean oin = m.queryPoint(p.multiply(a));
                     t.in = (mode==ITERATE_SUBTRACT) ? (t.in && !oin) : (t.in && oin);
                     break;
                     }
@@ -121,11 +127,12 @@ public final class Mesh {
                         Point p1 = e.v(1).multiply(a);
                         Point p2 = e.v(2).multiply(a);
                         Vertex v1=null, v2=null;
-                            v1 = m.vertex(p1);
-                            v2 = m.vertex(p2);
+                        v1 = m.vertex(p1);
+                        v2 = m.vertex(p2);
                         if (v1==v2) continue;
                         m.getEdge(v1, v2).lock(v1, 0);
                     }
+                    if (check) checkAllDelaunay();
                     break;
                 }
             }
@@ -138,6 +145,10 @@ public final class Mesh {
             }
         }
     }
+    public void checkAllDelaunay() {
+        for(int i=0; i<triangles.size(); i++)
+            ((Triangle)triangles.get(i)).checkDelaunay();
+    }
 
     public void clipOp(Mesh m, Affine a, boolean subtract) {
         seekTime = 0;
@@ -149,64 +160,111 @@ public final class Mesh {
         if (total > 80) System.out.println("clip in " + (100 * (seek/total)) + "%");
     }
 
-    // Geometry //////////////////////////////////////////////////////////////////////////////
-
-    public static double ddistance(double x1, double y1, double x2, double y2) {
-        return Math.sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));}
-    public static float area(Point p1, Point p2, Point p3) {
-        float x1 = p1.x;
-        float x2 = p2.x;
-        float x3 = p3.x;
-        float y1 = p1.y;
-        float y2 = p2.y;
-        float y3 = p3.y;
-        double a = ddistance(x1,y1,x2,y2);
-        double b = ddistance(x2,y2,x3,y3);
-        double c = ddistance(x3,y3,x1,y1);
-        double s = (a+b+c)/2;
-        double t = s*(s-a)*(s-b)*(s-c);
-        if (t < 0) return 0;
-        return (float)Math.sqrt(t);
-    }
+    // Imprecise Geometry //////////////////////////////////////////////////////////////////////////////
+    // (all of these can be off by a good margin and degrade only performance, not correctness) ////////
+    ////////////////////////////////////////////////////////////////////////////////////////////////////
 
-    public Point intersect(Point v1, Point v2, Point v3, Point v4) {
-        double a1 = v2.y-v1.y;
-        double a2 = v4.y-v3.y;
-        double b1 = v1.x-v2.x;
-        double b2 = v3.x-v4.x;
-        double c1 = -1 * (a1*v1.x+b1*v1.y);
-        double c2 = -1 * (a2*v3.x+b2*v3.y);
-        double x = (b2*c1-c2*b1)/(b1*a2-b2*a1);
-        double y = (a2*c1-c2*a1)/(a1*b2-a2*b1);
-        if (Double.isNaN(x) || Double.isNaN(y)) throw new Error("cannot intersect:\n ");
-        return point((float)x,(float)y);
+    private static class Imprecise {
+        public static Point center(Mesh m, Point v1, Point v2, Point v3) {
+            return m.point((float)((double)v1.x+(double)v2.x+(double)v3.x)/3,
+                           (float)(((double)v1.y+(double)v2.y+(double)v3.y)/3));
+        }
+        public static boolean near(Point a, Point b) {
+            return ddistance(a.x, a.y, b.x, b.y) <= epsilon;
+        }
+        public static boolean incircle(Point v1, Point v2, Point v3, Point p) {
+            /*
+            float a = 0;
+            Q: for(int q=0; q<2; q++) {
+                for(int i=0; i<3; i++) {
+                    if ((a=(v2.y-v3.y)*(v2.x-v1.x)-(v2.y-v1.y)*(v2.x-v3.x))!=0) break Q;
+                    Point t = v2; v2=v3; v3=v1; v1 = t; 
+                }
+                Point t = v2; v2=v3; v3=t;
+            }
+            if (a==0) throw new Error("a==0 for " + v1 + " " + v2 + " " + v3);
+            double a1  = (v1.x+v2.x)*(v2.x-v1.x)+(v2.y-v1.y)*(v1.y+v2.y);
+            double a2  = (v2.x+v3.x)*(v2.x-v3.x)+(v2.y-v3.y)*(v2.y+v3.y);
+            double ccx = (a1*(v2.y-v3.y)-a2*(v2.y-v1.y))/a/2;
+            double ccy = (a2*(v2.x-v1.x)-a1*(v2.x-v3.x))/a/2;
+            double r2  = (v1.x-ccx)*(v1.x-ccx)+(v1.y-ccy)*(v1.y-ccy);
+            double pd  = (p.x-ccx)*(p.x-ccx)+(p.y-ccy)*(p.y-ccy);
+            return r2 > pd;
+            */
+            return Predicates.incircle(v1.x, v1.y, v2.x, v2.y, v3.x, v3.y, p.x, p.y)>0;
+        }
+        public static Point midpoint(Mesh m, Point a, Point b) { return m.point((a.x+b.x)/2,(a.y+b.y)/2); }
+        private static double ddistance(double x1, double y1, double x2, double y2) {
+            return Math.sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));}
+        public static float area(Point p1, Point p2, Point p3) {
+            float x1 = p1.x;
+            float x2 = p2.x;
+            float x3 = p3.x;
+            float y1 = p1.y;
+            float y2 = p2.y;
+            float y3 = p3.y;
+            double a = ddistance(x1,y1,x2,y2);
+            double b = ddistance(x2,y2,x3,y3);
+            double c = ddistance(x3,y3,x1,y1);
+            double s = (a+b+c)/2;
+            double t = s*(s-a)*(s-b)*(s-c);
+            if (t < 0) return 0;
+            return (float)Math.sqrt(t);
+        }
+        
+        public static Point intersect(Mesh m, Point v1, Point v2, Point v3, Point v4) {
+            double a1 = v2.y-v1.y;
+            double a2 = v4.y-v3.y;
+            double b1 = v1.x-v2.x;
+            double b2 = v3.x-v4.x;
+            double c1 = -1 * (a1*v1.x+b1*v1.y);
+            double c2 = -1 * (a2*v3.x+b2*v3.y);
+            double x = (b2*c1-c2*b1)/(b1*a2-b2*a1);
+            double y = (a2*c1-c2*a1)/(a1*b2-a2*b1);
+            if (Double.isNaN(x) || Double.isNaN(y)) throw new Error("cannot intersect:\n ");
+            return m.point((float)x,(float)y);
+        }
+        public static int side(Point p1, Point p2, Point p3) {
+            /*
+            int ret = 0;
+            float x0 = p1.x;
+            float x = p2.x;
+            float x2 = p3.x;
+            float y0 = p1.y;
+            float y = p2.y;
+            float y2 = p3.y;
+            
+            // this MUST be done to double precision
+            double a = y-y0, b = x0-x, c = a*(x0 - x2) + b*(y0 - y2);
+            if      (c > 0) ret = b>=0 ? -1 :  1;
+            else if (c < 0) ret = b>=0 ?  1 : -1;
+            else ret = 0;
+            return ret;
+            */
+            return Predicates.side(p1.x, p1.y, p2.x, p2.y, p3.x, p3.y);
+        }
+        
     }
 
-    public static int side(Point p1, Point p2, Point p3) {
-        float x0 = p1.x;
-        float x = p2.x;
-        float x2 = p3.x;
-        float y0 = p1.y;
-        float y = p2.y;
-        float y2 = p3.y;
-        // this MUST be done to double precision
-        double a = y-y0, b = x0-x, c = a*(x0 - x2) + b*(y0 - y2);
-        if (c > 0) return b>=0 ? -1 :  1;
-        if (c < 0) return b>=0 ?  1 : -1;
-        return 0;
+    // Precise Geometry /////////////////////////////////////////////////////////////////////////////
+
+    static {
+        Runtime.getRuntime().load(new java.io.File("Predicates.jnilib").getAbsolutePath());
     }
 
+
     // Vertex //////////////////////////////////////////////////////////////////////////////
 
+    public static int state = 0;
     public Vertex vertex(Point p, Affine a) { return vertex(p.multiply(a)); }
     public Vertex vertex(Point p) {
         Vertex ret = null;
         switch(numvertices) {
             case 0: return (vertex0 = new Vertex(p));
-            case 1: return vertex0.distance(p)<=epsilon ? vertex0 : (vertex1 = new Vertex(p));
+            case 1: return Imprecise.near(vertex0,p) ? vertex0 : (vertex1 = new Vertex(p));
             case 2: {
-                if (vertex0.distance(p)<=epsilon) return vertex0;
-                if (vertex1.distance(p)<=epsilon) return vertex1;
+                if (Imprecise.near(vertex0,p)) return vertex0;
+                if (Imprecise.near(vertex1,p)) return vertex1;
                 Vertex v2 = new Vertex(p);
                 triangle(newEdge(vertex0,vertex1), newEdge(vertex1,v2), newEdge(v2,vertex0));
                 return v2;
@@ -216,7 +274,23 @@ public final class Mesh {
                     t = triangle0.seek(p);
                 for(int i=1; i<=3; i++)
                     for(int j=1; j<=2; j++)
-                        if (t != null && t.e(i).v(j).distance(p)<=epsilon) return t.e(i).v(j);
+                        if (t != null && Imprecise.near(t.e(i).v(j),p))
+                            return t.e(i).v(j);
+                // this will probably always need to be here since a vertex for which side()==0 could still
+                // be slightly off to one side (and hence part of a neighboring "sliver" triangle.
+                for(int k=1; k<=3; k++)
+                    for(int i=1; i<=3; i++)
+                        for(int j=1; j<=2; j++)
+                            if (t != null && t.t(k)!=null && Imprecise.near(t.t(k).e(i).v(j),p))
+                                return t.t(k).e(i).v(j);
+                for(int i=0; i<vertices.size(); i++)
+                    if (Imprecise.near(((Vertex)vertices.get(i)),p)) 
+                        throw new Error("bah! " + p + " " + vertices.get(i) + " state: " + state +
+                                        "\n  "   + t +
+                                        "\n    " + t.t(1) +
+                                        "\n    " + t.t(2) +
+                                        "\n    " + t.t(3)
+                                        );
                 Vertex v = new Vertex(p);
                 if      (t.e(3).intersects(p))          t.e(3).bisect(v);
                 else if (t.e(1).intersects(p))          t.e(1).bisect(v);
@@ -228,9 +302,10 @@ public final class Mesh {
         }
     }
 
+    Vector vertices = new Vector();
+
     private Point point(float x, float y) { return new Point(x,y); }
-    private Point point(Point a, Point b) { return new Point((a.x+b.x)/2,(a.y+b.y)/2); }
-    private Point point(Point p, Affine a) { return point(p.x(a), p.y(a)); }
+    //private Point point(Point p, Affine a) { return point(p.x(a), p.y(a)); }
     private class Point {
         public float x;
         public float y;
@@ -238,19 +313,13 @@ public final class Mesh {
         public Point multiply(Affine a) { return point(x(a),y(a)); }
         public Point(float x, float y) { this.x = x; this.y = y; }
         public Point(Point p) { this(p.x,p.y); }
-        public boolean equals(float x, float y) { return distance(x,y) <= epsilon; }
-        public boolean equals(Object o) { return (!(o instanceof Point)) ? false : ((Point)o).distance(this) <= epsilon; }
-        public float distance(Point v) { return distance(v.x,v.y); }
-        private float distance(float x, float y) { return (float)Math.sqrt(distance2(x, y)); }
-        public float distance2(Point v) { return distance2(v.x, v.y); }
-        private float distance2(float x, float y) { return (this.x-x)*(this.x-x)+(this.y-y)*(this.y-y); }
         public float x(Affine a) { return a.multiply_px(x,y); }
         public float y(Affine a) { return a.multiply_py(x,y); }
         public int xi(Affine a) { return (int)x(a); }
         public int yi(Affine a) { return (int)y(a); }
         public boolean intersects(Point p1, Point p2) {
             return
-                Mesh.side(p1,p2,this)==0 &&
+                Imprecise.side(p1,p2,this)==0 &&
                 x <= Math.max(p1.x,p2.x) &&
                 x >= Math.min(p1.x,p2.x) &&
                 y <= Math.max(p1.y,p2.y) &&
@@ -259,7 +328,11 @@ public final class Mesh {
     }
 
     private final class Vertex extends Point implements org.ibex.classgen.opt.Arena.Gladiator {
-        public Vertex(Point p) { super(p); numvertices++; }
+        public Vertex(Point p) {
+            super(p);
+            numvertices++;
+            vertices.add(this);
+        }
     }
 
     // Edge //////////////////////////////////////////////////////////////////////////////
@@ -277,11 +350,12 @@ public final class Mesh {
 
     public Edge getEdge(Vertex v1, Vertex v2) {
         if (v1==v2) throw new Error();
-        //Edge ret = (Edge)edges.get(v1,v2);
-        Edge ret = null;
+        Edge ret = (Edge)edges.get(v1,v2);
+        if (ret != null) return ret;
+        //Edge ret = null;
         Triangle t = null;
         if (triangle0 != null) {
-            t = triangle0.seek(point(v1,v2));
+            t = triangle0.seek(Imprecise.midpoint(this, v1,v2));
             if (t != null)
                 for(int i=1; i<=3; i++)
                     if (t.e(i).hasVertex(v1) && t.e(i).hasVertex(v2)) ret = t.e(i);
@@ -309,7 +383,13 @@ public final class Mesh {
             edges.put(v2,v1,null);
         }
 
-        public Vertex v(int i) { return i==1?v1:i==2?v2:null; }
+        public Vertex v(int i) {
+            switch(i) {
+                case 1: return v1;
+                case 2: return v2;
+                default: return null;
+            }
+        }
 
         public Edge rotate(Vertex v, boolean clockwise) {
             Triangle t = v==v1 ? (clockwise?t1:t2) : v==v2 ? (clockwise?t2:t1) : null;
@@ -326,7 +406,7 @@ public final class Mesh {
             return e;
         }
 
-        public boolean isNear(Point p) { return area(v1,v2,p) < epsilon2; }
+        public boolean isNear(Point p) { return Imprecise.area(v1,v2,p) <= epsilon2; }
         public Vertex commonVertex(Edge e) { return v1==e.v(1) || v1==e.v(2) ? v1 : v2==e.v(1) || v2==e.v(2) ? v2 : null; }
         public Vertex unCommonVertex(Edge e) { return v1!=e.v(1) && v1!=e.v(2) ? v1 : v2!=e.v(1) && v2!=e.v(2) ? v2 : null; }
         public Vertex opposingVertex(Vertex v) { return v1==v ? v2 : v1; }
@@ -335,7 +415,7 @@ public final class Mesh {
         public boolean locked() { return locks > 0; }
         public boolean partitions(Point va, Point vb) { return side(va,vb)==-1; }
         public int     side(Point a, Point b) { return side(a) * side(b); }
-        public int     side(Point a)          { return Mesh.side(v(1), v(2), a); }
+        public int     side(Point a)          { return Imprecise.side(v(1), v(2), a); }
         public boolean hasVertex(Vertex v) { return v1==v || v2==v; }
         public boolean hasTriangle(Triangle t) { return t==t1 || t==t2; }
         public String toString() { return v(1) + "--" + v(2); }
@@ -347,26 +427,18 @@ public final class Mesh {
         }
         public boolean convex() { return this.intersects(t1.opposingVertex(t2), t2.opposingVertex(t1)); }
 
-        public boolean colinear(Point v) { return area(v,v1,v2)<=epsilon; }
-
-        public boolean intersects(Point p) {
-            return
-                side(p)==0 &&
-                p.x <= Math.max(v1.x,v2.x) &&
-                p.x >= Math.min(v1.x,v2.x) &&
-                p.y <= Math.max(v1.y,v2.y) &&
-                p.y >= Math.min(v1.y,v2.y);
-        }
+        public boolean intersects(Point p) { return p.intersects(v1, v2); }
         public boolean intersects(Edge e) { return intersects(e.v(1), e.v(2)); }
         public boolean intersects(Point va, Point vb) {
             return
-                !v1.equals(va) &&
-                !v1.equals(vb) &&
-                !v2.equals(va) &&
-                !v2.equals(vb) &&
+                !Imprecise.near(v1,va) &&
+                !Imprecise.near(v1,vb) &&
+                !Imprecise.near(v2,va) &&
+                !Imprecise.near(v2,vb) &&
                 partitions(va, vb) &&
-                Mesh.side(va, vb, v1) * Mesh.side(va, vb, v2) == -1;
+                Imprecise.side(va, vb, v1) * Imprecise.side(va, vb, v2) == -1;
         }
+
         public Triangle opposingTriangle(Triangle t) {
             if (t1 == t) return t2;
             if (t2 == t) return t1;
@@ -400,6 +472,8 @@ public final class Mesh {
             }
         }
         public void bisect(Vertex v) {
+            if (v==this.v(1)) throw new Error("this should never happen");
+            if (v==this.v(2)) throw new Error("this should never happen");
             Edge e = this;
             Triangle t1 = this.t1==null?this.t2:this.t1;
             Triangle t  = t1==this.t1?this.t2:this.t1;
@@ -429,26 +503,36 @@ public final class Mesh {
                 tov = t.opposingVertex(t1);
                 in2 = t.in;
                 tov_v = newEdge(tov, v);
-                if (top == t1) top = null;
-                if (ton == t1) ton = null;
-                if (opposingTriangleLeft == t) opposingTriangleLeft = null;
-                if (opposingTriangleRight == t) opposingTriangleRight = null;
+                if (top == t1) top = null;                                    // is this possible?
+                if (ton == t1) ton = null;                                    // is this possible?
+                if (opposingTriangleLeft == t) opposingTriangleLeft = null;   // is this possible?
+                if (opposingTriangleRight == t) opposingTriangleRight = null; // is this possible?
                 t.delete();
             }
             t1.delete();
-            Triangle ta, tb, tc, td;
+            Triangle ta, tb, tc=null, td=null;
             ta = triangle(right_opposing, opposing_v, right_v);
             tb = triangle(left_opposing, opposing_v, left_v);
             ta.in = in1;
             tb.in = in1;
             if (t != null) {
+                if (tov_v==left_v) throw new Error("barf");
+                if (tov_v==right_v) throw new Error("barf");
+                if (tov_v==left_tov) throw new Error("barf");
+                if (tov_v==right_tov) throw new Error("barf");
+                if (right_v==right_tov) throw new Error("barf");
+                if (left_v==left_tov) throw new Error("barf " + tov + " " + left);
                 tc = triangle(left_tov, tov_v, left_v);
                 td = triangle(right_tov, tov_v, right_v);
                 tc.in = in2;
                 td.in = in2;
             }
             if (locked()) fracture(v);
-            else ta.fixup();
+            else          ta.fixup();
+            if (ta!=null) ta.check();
+            if (tb!=null) tb.check();
+            if (tc!=null) tc.check();
+            if (td!=null) td.check();
         }
         public Edge flip() {
             if (locked()) throw new Error("attempted to remove a locked edge: " + this);
@@ -489,7 +573,7 @@ public final class Mesh {
         }
         public void fracture(Edge e) {
             triangle0=e.t1==null?e.t2:e.t1;
-            Vertex v0 = vertex(Mesh.this.intersect(v1,v2,e.v(1),e.v(2)));
+            Vertex v0 = vertex(Imprecise.intersect(Mesh.this, v1,v2,e.v(1),e.v(2)));
             if (v0 != e.v(1) && v0 != e.v(2) && e.locked()) e.fracture(v0);
             if (v0 != v1 && v0 != v2) fracture(v0);
         }
@@ -533,6 +617,11 @@ public final class Mesh {
                         while(t.intersects(this)) {
                             if (t==told) break;
                             told = t;
+                            /*
+                            System.out.println("I think that " + this + " intersects:\n  "+t);
+                            for(int i=1; i<=3; i++)
+                                System.out.println("    " + t.e(i) + ": " + t.e(i).intersects(this));
+                            */
                             t = t.followVector(v2,v1);
                         }
                         t = told;
@@ -545,8 +634,12 @@ public final class Mesh {
             if (t2!=null) t2.fixup();
         }
 
+        public boolean violated = false;
         public void stroke(PixelBuffer buf, Affine a, int color) {
-            int c = debug
+            int c = 
+                violated
+                ? 0xffff0000
+                : debug
                 ? (weight() == 0 ? color : 0xffff0000)
                 : (weight() != 0 ? color : 0);
             if (c != 0) buf.drawLine(v1.xi(a), v1.yi(a), v2.xi(a), v2.yi(a), c);
@@ -563,29 +656,27 @@ public final class Mesh {
     // Triangle //////////////////////////////////////////////////////////////////////////////
 
     public Triangle triangle(Edge e1, Edge e2, Edge e3) {
-        float x = (e1.v(1).x+e1.v(2).x+e2.v(1).x+e2.v(2).x+e3.v(1).x+e3.v(2).x)/6;
-        float y = (e1.v(1).y+e1.v(2).y+e2.v(1).y+e2.v(2).y+e3.v(1).y+e3.v(2).y)/6;
-        Point p = point(x,y);
-        Triangle t = triangle0==null ? null : triangle0.seek(p);
-        if (t != null &&
-            (t.contains(p) || t.intersects(p)) &&
-            t.hasEdge(e1) &&
-            t.hasEdge(e2) &&
-            t.hasEdge(e3))
-            return triangle0 = t;
-        t = new Triangle(e1, e2, e3);
+        if (e3.t1!=null && e3.t1.hasEdge(e2) && e3.t1.hasEdge(e1)) return e3.t1;
+        if (e3.t2!=null && e3.t2.hasEdge(e2) && e3.t2.hasEdge(e1)) return e3.t2;
+        if (e2.t1!=null && e2.t1.hasEdge(e2) && e2.t1.hasEdge(e1)) return e2.t1;
+        if (e2.t2!=null && e2.t2.hasEdge(e2) && e2.t2.hasEdge(e1)) return e2.t2;
+        if (e1.t1!=null && e1.t1.hasEdge(e2) && e1.t1.hasEdge(e1)) return e1.t1;
+        if (e1.t2!=null && e1.t2.hasEdge(e2) && e1.t2.hasEdge(e1)) return e1.t2;
+        Triangle t = new Triangle(e1, e2, e3);
         if (debug) t.check();
         if (triangle0 == null) triangle0 = t;
         return t;
     }
 
-        public static boolean fixing = false;
+    public static boolean fixing = false;
+    public static int count = 0;
     private final class Triangle implements org.ibex.classgen.opt.Arena.Gladiator {
 
-        final float r2;
-        final Point cc;
+        //final double r2;
+        //final Point cc;
 
         private Edge e1, e2, e3;    // should be final =(
+        private Vertex v1, v2, v3;
         public int tick;
 
         boolean in = false;
@@ -593,14 +684,35 @@ public final class Mesh {
         boolean painted = false;
         boolean dirty = true;
 
-        public Edge e(int i) { return i==1?e1:i==2?e2:i==3?e3:null; }
-        public Vertex v(int i) { return e(i==1?2:i==2?3:i==3?1:0).unCommonVertex(e(i)); }
+        //public Edge e(int i) { return i==1?e1:i==2?e2:i==3?e3:null; }
+        public Edge e(int i) {
+            switch(i) {
+                case 1: return e1;
+                case 2: return e2;
+                case 3: return e3;
+                default: return null;
+            } }
+
+        //public Vertex v(int i) { return e(i==1?2:i==2?3:i==3?1:0).unCommonVertex(e(i)); }
+        public Vertex v(int i) {
+            switch(i) {
+                case 1: return v1;
+                case 2: return v2;
+                case 3: return v3;
+                default: return null;
+            } }
         public Triangle t(int i) { return e(i).t1==this ? e(i).t2 : e(i).t1; }
-
+        public Vertex closestVertex(Point p) {
+            double d1 = Math.sqrt((v(1).x-p.x)*(v(1).x-p.x)+(v(1).y-p.y)*(v(1).y-p.y));
+            double d2 = Math.sqrt((v(2).x-p.x)*(v(2).x-p.x)+(v(2).y-p.y)*(v(2).y-p.y));
+            double d3 = Math.sqrt((v(3).x-p.x)*(v(3).x-p.x)+(v(3).y-p.y)*(v(3).y-p.y));
+            return (d1 < d2 && d1 < d3) ? v(1) : (d2 < d3) ? v(2) : v(3);
+        }
+        
         public boolean encounters(Point p1, Point p2) {
             for(int i=1; i<=3; i++) {
-                if (v(i).equals(p1)) return true;
-                if (v(i).equals(p2)) return true;
+                if (Imprecise.near(v(i),p1)) return true;
+                if (Imprecise.near(v(i),p2)) return true;
                 if (v(i).intersects(p1,p2)) return true;
                 if (e(i).intersects(p1,p2)) return true;
             }
@@ -608,9 +720,7 @@ public final class Mesh {
         }
         public Edge getSharedEdge(Triangle t) { return e(1).t1==t||e(1).t2==t?e(1):e(2).t1==t||e(2).t2==t?e(2):e(3).t1==t||e(3).t2==t?e(3):null; }
         public boolean contains(Point p) { return e(1).side(v(1),p)==1 && e(2).side(v(2),p)==1 && e(3).side(v(3),p)==1; }
-        public Point c()  { return point(cx(),cy()); }
-        public float cx() { return (float)(((double)v(1).x+(double)v(2).x+(double)v(3).x)/3); }
-        public float cy() { return (float)(((double)v(1).y+(double)v(2).y+(double)v(3).y)/3); }
+
         public boolean intersects(Vertex va, Vertex vb){return e(1).intersects(va,vb)||e(2).intersects(va,vb)||e(3).intersects(va,vb);}
         public boolean intersects(Edge e){ return intersects(e.v(1),e.v(2)); }
         public boolean intersects(Point p){ return e(1).intersects(p) || e(2).intersects(p) || e(3).intersects(p); }
@@ -632,7 +742,7 @@ public final class Mesh {
             for(int i=1; i<=3; i++) {
                 Triangle t = t(i);
                 if (t==null) continue;
-                if (t.r2 <= v(i).distance2(t.cc)) continue;
+                if (!t.cc(v(i))) continue;
                 Edge e = e(i);
                 if (e.locked()) { t.fixup(); continue; }
                 return e.flip().t1.fixup();
@@ -666,14 +776,29 @@ public final class Mesh {
             Triangle t = this;
             try {
             while (true) {
-                if      (t.contains(p) || t.intersects(p))  return t;
-                else if (t.e(3).intersects(p))                return (t.t(3)!=null && t.t(3).contains(p)) ? t.t(3) : t;
-                else if (t.e(1).intersects(p))                return (t.t(1)!=null && t.t(1).contains(p)) ? t.t(1) : t;
-                else if (t.e(2).intersects(p))                return (t.t(2)!=null && t.t(2).contains(p)) ? t.t(2) : t;
+                count++;
+                //System.out.println("seek " + t + " -> " + p + " / " + count);
+                if      (t.contains(p)) { state = -1; return t; }
+                else if (t.intersects(p))    { state = 1; return t; }
+                else if (t.e(3).intersects(p))                { state = 2; return (t.t(3)!=null && t.t(3).contains(p)) ? t.t(3) : t; }
+                else if (t.e(1).intersects(p))                { state = 3; return (t.t(1)!=null && t.t(1).contains(p)) ? t.t(1) : t; }
+                else if (t.e(2).intersects(p))                {state = 4;  return (t.t(2)!=null && t.t(2).contains(p)) ? t.t(2) : t; }
                 else {
-                    Triangle t2 = t.followVector(t.c(), p);
-                    if (t2==null || t2==t) return t;
-                    t = t2;
+                    // we "slingshot" back from the centroid in case we're inside of a "sliver" triangle 
+                    //Point p0 = t.c();
+                    //Triangle t2 = t.followVector(p0, p);
+                    Triangle t2 = t.followVector(t.closestVertex(p), p);
+                    if (t2==null || t2==t) {
+                        if      (t.e(1).partitions(p, t.v(1)) && t.t(1)!=null) t = t.t(1);
+                        else if (t.e(2).partitions(p, t.v(2)) && t.t(2)!=null) t = t.t(2);
+                        else if (t.e(3).partitions(p, t.v(3)) && t.t(3)!=null) t = t.t(3);
+                        else {
+                            state = 5;
+                            return t;
+                        }
+                    } else {
+                        t = t2;
+                    }
                 }
             }
             } finally { if (t!=null) triangle0 = t; }
@@ -688,8 +813,9 @@ public final class Mesh {
             return ret;
         }
         public Triangle followVector2(Point p1, Point p2) {
-            if (contains(p2) || intersects(p2) || v(1).equals(p2) || v(2).equals(p2) || v(3).equals(p2)) return this;
-            for(int i=1; i<=3; i++) if (!v(i).equals(p1) && v(i).intersects(p1,p2)) return followVector(v(i),p2);
+            if (contains(p2) || intersects(p2) || Imprecise.near(v(1),p2) || Imprecise.near(v(2),p2) || Imprecise.near(v(3),p2))
+                return this;
+            for(int i=1; i<=3; i++) if (!Imprecise.near(v(i),p1) && v(i).intersects(p1,p2)) return followVector(v(i),p2);
             Triangle t1 = t(1);
             Triangle t2 = t(2);
             Triangle t3 = t(3);
@@ -697,24 +823,67 @@ public final class Mesh {
                 int k1 = i==1?3:i==2?1:i==3?2:0;
                 int k2 = i==1?2:i==2?3:i==3?1:0;
                 int k3 = i==1?1:i==2?2:i==3?3:0;
-                if (v(i).equals(p1)) {
+                if (Imprecise.near(v(i),p1)) {
                     if (e(k1).partitions(v(k1),p2)) return t(k1);
                     if (e(k2).partitions(v(k2),p2)) return t(k2);
                     if (e(k3).partitions(v(k3),p2)) return t(k3);
                     throw new Error("bad!");
                 }
             }
-            if (!e(1).intersects(p1,p2) && !e(2).intersects(p1,p2) && !e(3).intersects(p1,p2))
-                throw new Error("invoked followVector() on a Triangle which it does not encounter:\n" +
-                                "  p1=" + p1 + "\n" +
-                                "  p2=" + p2 + "\n" +
-                                "  t =" + this + " (area "+area(v(1),v(2),v(3))+")\n");
-            for(int i=1; i<=3; i++) if (e(i).intersects(p1,p2) && e(i).side(v(i)) * e(i).side(p2) == -1) return t(i);
-            throw new Error("giving up: \n  "+p1+" -> "+p2+"\n  " + this);
+            //if (!e(1).intersects(p1,p2) && !e(2).intersects(p1,p2) && !e(3).intersects(p1,p2))
+            for(int i=1; i<=3; i++)
+                if (e(i).intersects(p1,p2))
+                    if (e(i).side(v(i)) * e(i).side(p2) == -1)
+                        return t(i);
+            for(int i=1; i<=3; i++)
+                if (e(i).partitions(p1,p2))
+                    return t(i);
+            for(int i=1; i<=3; i++)
+                if (e(i).partitions(v(i),p2))
+                    return t(i);
+            for(int i=1; i<=3; i++)
+                if (v(i).intersects(p1,p2))
+                    throw new Error("bad news: \n  "+p1+" -> "+p2+"\n  " + this);
+
+            System.out.println("slingshot from: " + p1 + " to " + p2 + " on " + this + "\n" +
+                               (e(1).side(v(1)) * e(1).side(p2))+" "+
+                               (e(2).side(v(2)) * e(2).side(p2))+" "+
+                               (e(3).side(v(3)) * e(3).side(p2))
+                               );
+            /*
+            return followVector(new Point(2*p1.x-p2.x, 2*p1.y-p2.y), p2);
+            */
+            //throw new Error("giving up: \n  "+p1+" -> "+p2+"\n  " + this);
+
+            final Point pp1 = p1;
+            final Point pp2 = p2;
+            new Frame() {
+                public void paint(Graphics g) {
+                    g.setColor(java.awt.Color.white);
+                    g.fillRect(0, 0, getWidth(), getHeight());
+                    g.setColor(java.awt.Color.black);
+                    g.drawLine((int)v(1).x+100, (int)v(1).y+100, (int)v(2).x+100, (int)v(2).y+100);
+                    g.drawLine((int)v(3).x+100, (int)v(3).y+100, (int)v(2).x+100, (int)v(2).y+100);
+                    g.drawLine((int)v(1).x+100, (int)v(1).y+100, (int)v(3).x+100, (int)v(3).y+100);
+                    
+                    g.setColor(java.awt.Color.red);
+                    g.drawLine((int)pp1.x+100, (int)pp1.y+100, (int)pp2.x+100, (int)pp2.y+100);
+                }
+            }.show();
+            try { Thread.sleep(100000); } catch (Exception e) { }
+            return null;
+
+            /*
+            throw new Error("invoked followVector() on a Triangle which it does not encounter:\n" +
+                            "  p1=" + p1 + "\n" +
+                            "  p2=" + p2 + "\n" +
+                            "  t =" + this + " (area "+area(v(1)+100,v(2),v(3))+")\n");
+            */
         }
 
         public void check() {
-            if (debug) {
+            if (e1==null && e2==null && e3==null) return;
+            if (check) {
                 for(int i=1; i<=3; i++) {
                     if (e(i).v(1) != v(1) && e(i).v(1) != v(2) && e(i).v(1) != v(3)) throw new Error("inconsistent");
                     if (e(i).t1 != this && e(i).t2 != this) throw new Error("inconsistent");
@@ -722,21 +891,40 @@ public final class Mesh {
                 }
                 if (e(1)==e(2) || e(2)==e(3) || e(3)==e(1)) throw new Error("identical edges");
                 for(int i=1; i<=3; i++) {
-                    if (t(i) != null) if (!t(i).hasEdge(e(i))) throw new Error("t1 doesn't have e(1)");
-                    if (t(i) != null) {
-                        if (t(i).getSharedEdge(this) != e(i)) throw new Error("blark");
-                        if (!e(i).hasTriangle(t(i))) throw new Error("blark2");
-                        if (!e(i).hasTriangle(this)) throw new Error("blark3");
-                    }
+                    if (t(i) == null) continue;
+                    if (!t(i).hasEdge(e(i))) throw new Error("t1 doesn't have e(1)");
+                    if (t(i).getSharedEdge(this) != e(i)) throw new Error("blark");
+                    if (!e(i).hasTriangle(t(i))) throw new Error("blark2");
+                    if (!e(i).hasTriangle(this)) throw new Error("blark3");
+                }
+                for(int i=1; i<=3; i++)
+                    if (e(i).commonVertex(e(i==3?1:(i+1)))==null)
+                        throw new Error("edges have no common vertex");
+                // check that delauanay property is preserved
+            }
+        }
+
+        public void checkDelaunay() {
+            for(int i=1; i<=3; i++) {
+                if (t(i) == null) continue;
+                Vertex v = t(i).opposingVertex(e(i));
+                if (!e(i).locked() && /*Imprecise.incircle(v(1), v(2), v(3), v)*/cc(v) /*&& !dirty && !t(i).dirty*/) {
+                    //throw new Error("Delaunay violation: vertex " + v + "\n  triangle: " + this);
+                    //System.out.println("violation: " + e(i));
+                    e(i).violated = true;
+                } else {
+                    e(i).violated = false;
                 }
-                // check that edges all join up
             }
         }
 
         public void trisect(Vertex v) {
             if (!contains(v)) throw new Error("trisect(v) but I don't contain v = " + v);
             if (hasVertex(v)) throw new Error("attempt to trisect a triangle at one of its own vertices");
-            for(int i=3; i>0; i--) if (e(i).isNear(v)) { e(i).bisect(v); return; }
+            for(int i=3; i>0; i--) if (e(i).intersects(v)/* || e(i).isNear(v)*/) {
+                e(i).bisect(v);
+                return;
+            }
             Triangle a=null,b=null,c=null;
 
             boolean oldIn = in;
@@ -755,6 +943,9 @@ public final class Mesh {
             b.in = oldIn;
             c.in = oldIn;
             a.fixup();
+            a.check();
+            b.check();
+            c.check();
         }
 
         public void setIn(boolean evenOdd, int weight) {
@@ -771,7 +962,7 @@ public final class Mesh {
             for(int i=1; i<=3; i++)
                 if (t(i) != null) {
                     boolean prepaint = t(i).painted;
-                    if (debug) e(i).stroke(buf, a, color);
+                    //if (debug) e(i).stroke(buf, a, color);
                     t(i).fill(buf, a, clip, color, strokeOnly);
                 }
         }
@@ -780,31 +971,25 @@ public final class Mesh {
             this.e1 = e1;
             this.e2 = e2;
             this.e3 = e3;
-            Vertex v1 = e(2).unCommonVertex(e(1));
-            Vertex v2 = e(3).unCommonVertex(e(2));
-            Vertex v3 = e(1).unCommonVertex(e(3));
+            /*Vertex*/ this.v1 = e(2).unCommonVertex(e(1));
+            /*Vertex*/ this.v2 = e(3).unCommonVertex(e(2));
+            /*Vertex*/ this.v3 = e(1).unCommonVertex(e(3));
             if (e(1).intersects(v1)) throw new Error("triangle points are colinear");
             if (e(2).intersects(v2)) throw new Error("triangle points are colinear");
             if (e(3).intersects(v3)) throw new Error("triangle points are colinear");
             e(1).addTriangle(this);
             e(2).addTriangle(this);
             e(3).addTriangle(this);
-
-            float a = 0;
-            Q: for(int q=0; q<2; q++) {
-                for(int i=0; i<3; i++) {
-                    if ((a=(v2.y-v3.y)*(v2.x-v1.x)-(v2.y-v1.y)*(v2.x-v3.x))!=0) break Q;
-                    Vertex t = v2; v2=v3; v3=v1; v1 = t; 
-                }
-                Vertex t = v2; v2=v3; v3=t;
-            }
-            if (a==0) throw new Error("a==0 for " + v1 + " " + v2 + " " + v3);
-            float a1=(v1.x+v2.x)*(v2.x-v1.x)+(v2.y-v1.y)*(v1.y+v2.y);
-            float a2=(v2.x+v3.x)*(v2.x-v3.x)+(v2.y-v3.y)*(v2.y+v3.y);
-            cc=point((a1*(v2.y-v3.y)-a2*(v2.y-v1.y))/a/2, (a2*(v2.x-v1.x)-a1*(v2.x-v3.x))/a/2);
-            r2 = v1.distance2(cc);
             triangles.add(this);
         }
+        public boolean cc(Point p) {
+            /*
+            Vertex v1 = e(2).unCommonVertex(e(1));
+            Vertex v2 = e(3).unCommonVertex(e(2));
+            Vertex v3 = e(1).unCommonVertex(e(3));
+            */
+            return Imprecise.incircle(v1, v2, v3, p);
+        }
         public void clear() {
             if (!painted) return;
             painted = false;
@@ -888,6 +1073,7 @@ public final class Mesh {
         } finally {
             last = vx;
         }
+        if (check) checkAllDelaunay();
     }
 
 }