checkpoint
[anneal.git] / src / edu / berkeley / qfat / Mesh.java
index 64724a6..57eda26 100644 (file)
@@ -51,7 +51,8 @@ public class Mesh implements Iterable<Mesh.T> {
     public void transform(Matrix m) {
         ArrayList<Vertex> set = new ArrayList<Vertex>();
         for(Vertex v : vertices) set.add(v);
-        for(Vertex v : set) v.transform(m.times(v.p), true);
+        for(Vertex v : set) v.transform(m.times(v.p), true, null);
+        for(Vertex v : set) v.goodp = v.p;
     }
 
     public void rebuild() { /*vertices.rebuild();*/ }
@@ -77,7 +78,8 @@ public class Mesh implements Iterable<Mesh.T> {
 
     /** a vertex in the mesh */
     public final class Vertex extends HasQuadric implements Visitor {
-        public Point p, oldp, goodp;
+        public Point p, goodp;
+        public Point oldp;
         E e;                // some edge *leaving* this point
 
         private boolean illegal = false;
@@ -88,6 +90,7 @@ public class Mesh implements Iterable<Mesh.T> {
         private Vertex(Point p) {
             this.p = p;
             this.goodp = p;
+            this.oldp = p;
             if (vertices.get(p) != null) throw new Error();
             vertices.add(this);
         }
@@ -144,11 +147,17 @@ public class Mesh implements Iterable<Mesh.T> {
             if (quadric_count != 0)
                 nerror = (nerror + quadric.preAndPostMultiply(p))/(quadric_count+1);
 
+            if (!immutableVertices && quadric_count == 0)
+                nerror *= 2;
+            //System.out.println(nerror);
             for(E e = this.e; e!=null; e=e.pair.next==this.e?null:e.pair.next) {
-                double ang = Math.abs(e.dihedralAngle());
+                double ang = e.dihedralAngle();
                 if (ang > Math.PI) throw new Error();
-                float minangle = (float)(Math.PI * 0.8);
-                if (ang > minangle) nerror += (ang - minangle);
+                if (ang < -Math.PI) throw new Error();
+                //float minangle = (float)(Math.PI * 0.8);
+                nerror += ((ang / Math.PI)*(ang/Math.PI)) * e.length() * 0.05;
+                //System.out.println(((ang / Math.PI)*(ang/Math.PI)) * 0.000001);
+                //if (ang > minangle) nerror += (ang - minangle);
                 /*
                 if (e.t.aspect() < 0.2) {
                     nerror += (0.2-e.t.aspect()) * 10;
@@ -167,33 +176,76 @@ public class Mesh implements Iterable<Mesh.T> {
             //     t2' = t2.getMatrix(t1) * M * t1
             //     t1 =     t1.getMatrix(t2) * t2
             // M * t1 = M * t1.getMatrix(t2) * t2
-            for(Vertex v : (Iterable<Vertex>)getBoundPeers())
-                good &= v.transform(v.getBindingMatrix(this).times(m).times(this.p),
-                                    ignoreProblems);
 
-            for(Vertex v : (Iterable<Vertex>)getBoundPeers())
-                if (good || ignoreProblems)  v.reComputeErrorAround();
-                else                         v.transform(v.oldp, true);
+            if (bindingGroup!=null && this != bindingGroup.getMaster()) {
+                Matrix v = getBindingMatrix(bindingGroup.getMaster());
+                return ((Vertex)bindingGroup.getMaster()).move(v.inverse().times(m).times(v), ignoreProblems);
+            }
+
+            if (bindingGroup != null) {
+                Matrix m2 = null;
+                for(int i=0; i<20 && !m.equals(m2); i++) {
+                    m2 = m.times(getConstraint());
+                    //System.out.println(m.minus(m2));
+                }
+                if (!m.equals(m2)) return true;
+            }
+
+            Point op = this.p;
+            Point pt = m.times(this.p);
+            for(Vertex v : (Iterable<Vertex>)getBoundPeers()) {
+                Point pt2 = v.getBindingMatrix(this).times(pt);
+                /*
+                if (Math.abs( v.p.minus(pt2).mag() / pt.minus(op).mag() ) > 5)
+                    throw new Error(v.p+" "+pt2+"\n"+op+" "+pt+"\n"+v.getBindingMatrix(this));
+                if (Math.abs( v.p.minus(pt2).mag() / pt.minus(op).mag() ) < 1/5) throw new Error();
+                */
+                good &= v.transform(pt2, ignoreProblems, v.getBindingMatrix(this));
+            }
+
+            if (!good && !ignoreProblems) {
+                for(Vertex v : (Iterable<Vertex>)getBoundPeers()) 
+                    v.transform(v.oldp, true, null);
+            }
 
+            for(Vertex v : (Iterable<Vertex>)getBoundPeers())
+                v.recomputeFundamentalQuadricIfNeighborChanged();
+            for(Vertex v : (Iterable<Vertex>)getBoundPeers())
+                v.reComputeErrorAround();
             return good;
         }
 
         /** does NOT update bound pairs! */
-        private boolean transform(Point newp, boolean ignoreProblems) {
+        private boolean transform(Point newp, boolean ignoreProblems, Matrix yes) {
             this.oldp = this.p;
             if (immutableVertices) throw new Error();
 
             unApplyQuadricToNeighbor();
+
+
+            boolean illegalbefore = illegal;
+            illegal = false;
+            /*
+            if (this.p.minus(newp).mag() > 0.1 && !ignoreProblems) {
+                try {
+                    throw new Exception(""+this.p.minus(newp).mag()+" "+ignoreProblems+" "+yes);
+                } catch(Exception e) {
+                    e.printStackTrace();
+                }
+                illegal = true;
+            }
+            */
+
             this.p = newp;
             reinsert();
             applyQuadricToNeighbor();
 
             if (!ignoreProblems) {
-                illegal = false;
                 checkLegality();
             }
-            for(E e = this.e; e!=null; e=e.pair.next==this.e?null:e.pair.next) e.p2.quadricStale = true;
-            return !illegal;
+            for(E e = this.e; e!=null; e=e.pair.next==this.e?null:e.pair.next)
+                e.p2.quadricStale = true;
+            return !illegal || (illegalbefore && illegal);
         } 
 
         public void checkLegality() {
@@ -201,7 +253,7 @@ public class Mesh implements Iterable<Mesh.T> {
             for(E e = this.e; e!=null; e=e.pair.next==this.e?null:e.pair.next) {
                 if (Math.abs(e.dihedralAngle()) > (Math.PI * 0.9) ||
                     Math.abs(e.next.dihedralAngle()) > (Math.PI * 0.9)) illegal = true;
-                if (e.t.aspect() < 0.1) illegal = true;
+                if (e.t.aspect() < 0.2) illegal = true;
             }
             */
             if (!illegal) triangles.range(oldp, this.p, (Visitor<T>)this);
@@ -272,40 +324,6 @@ public class Mesh implements Iterable<Mesh.T> {
         public void bindTo(Vertex p) { bindTo(Matrix.ONE, p); }
     }
 
-    public class BindingGroup {
-        private HashSet<E> set = new HashSet<E>();
-        public BindingGroup bind_others;
-        public BindingGroup other() { return bind_others; }
-        public BindingGroup(BindingGroup bind_others) { this.bind_others = bind_others; }
-        public BindingGroup() { this.bind_others = new BindingGroup(this); }
-        public BindingGroup(E e) { this(); set.add(e); }
-        public void add(E e) {
-            if (set.contains(e)) return;
-            set.add(e);
-            BindingGroup e_bind_peers = e.bind_peers;
-            BindingGroup e_bind_to    = e.bind_to;
-            e.bind_peers = this;
-            e.bind_to    = bind_others;
-            for (E epeer  : e_bind_peers.set) add(epeer);
-            for (E eother : e_bind_to.set)    bind_others.add(eother);
-            /*
-            for(E eother : bind_others.set) {
-                if (e.next.bind_to.set.contains(eother.prev)) {
-                    e.next.next.bindEdge(eother.prev.prev);
-                }
-                if (e.prev.bind_to.set.contains(eother.next)) {
-                    e.prev.prev.bindEdge(eother.next.next);
-                }
-            }
-            */
-        }
-        public void dobind(E e) {
-            for(E ebound : set) {
-                e.p1.bindTo(Matrix.ONE, ebound.p2);
-                e.p2.bindTo(Matrix.ONE, ebound.p1);
-            }
-        }
-    }
 
     /** [UNIQUE] an edge */
     public final class E extends HasBindingGroup implements Comparable<E> {
@@ -315,8 +333,6 @@ public class Mesh implements Iterable<Mesh.T> {
         E prev;  // previous half-edge
         E next;  // next half-edge
         E pair;  // partner half-edge
-        public BindingGroup bind_peers  = new BindingGroup(this);
-        public BindingGroup bind_to     = bind_peers.other();
         boolean shattered = false;
 
         public boolean intersects(T t) { return t.intersects(p1.p, p2.p); }
@@ -327,13 +343,21 @@ public class Mesh implements Iterable<Mesh.T> {
                 (edu.berkeley.qfat.geom.BindingGroup<E>)newBindingGroup_;
             if (newBindingGroup==null) return;
             if (this==newBindingGroup.getMaster()) return;
-            for(E eother : (Iterable<E>)newBindingGroup) {
+            HashSet<E> nbg = new HashSet<E>();
+            for(E eother : (Iterable<E>)newBindingGroup) nbg.add(eother);
+            for(E eother : nbg) {
                 if (next==null || prev==null) continue;
                 if (eother.next==null || eother.prev==null) continue;
+
                 if (next.isBoundTo(eother.pair.prev.pair) && !prev.isBoundTo(eother.pair.next.pair))
                     prev.bindTo(next.getBindingMatrix(eother.pair.prev.pair), eother.pair.next.pair);
                 if (!next.isBoundTo(eother.pair.prev.pair) && prev.isBoundTo(eother.pair.next.pair))
                     next.bindTo(prev.getBindingMatrix(eother.pair.next.pair), eother.pair.prev.pair);
+
+                if (next.isBoundTo(eother.next) && !prev.isBoundTo(eother.prev))
+                    prev.bindTo(next.getBindingMatrix(eother.next), eother.prev);
+                if (!next.isBoundTo(eother.next) && prev.isBoundTo(eother.prev))
+                    next.bindTo(prev.getBindingMatrix(eother.prev), eother.next);
             }
 
         }
@@ -358,8 +382,6 @@ public class Mesh implements Iterable<Mesh.T> {
             pair._bindEdge(e.pair, m);
         }
         public void _bindEdge(E e, Matrix m) {
-            //bind_to.add(e);
-
             e = e.pair;
             /*
             //assumes edges are identical length at binding time
@@ -385,11 +407,12 @@ public class Mesh implements Iterable<Mesh.T> {
         }
         
         public void dobind() {
-            //bind_to.dobind(this);
             for(E e : (Iterable<E>)getBoundPeers()) {
                 if (e==this) continue;
                 p1.bindTo(getBindingMatrix(e), e.p1);
                 p2.bindTo(getBindingMatrix(e), e.p2);
+                e.p1.setConstraint(getConstraint());
+                e.p2.setConstraint(getConstraint());
             }
         }
 
@@ -408,14 +431,14 @@ public class Mesh implements Iterable<Mesh.T> {
                 Point mid = e.midpoint();
                 Vertex r = e.next.p2;
                 Vertex l = e.pair.next.p2;
-                e.destroy();
-                e.pair.destroy();
-                newT(r.p, e.p1.p, mid,    null, 0);
-                newT(r.p, mid,    e.p2.p, null, 0);
-                newT(l.p, mid,    e.p1.p, null, 0);
-                newT(l.p, e.p2.p, mid,    null, 0);
-                /*
-                */
+                if (!e.destroyed) {
+                    e.destroy();
+                    e.pair.destroy();
+                    newT(r.p, e.p1.p, mid,    null, 0);
+                    newT(r.p, mid,    e.p2.p, null, 0);
+                    newT(l.p, mid,    e.p1.p, null, 0);
+                    newT(l.p, e.p2.p, mid,    null, 0);
+                }
             }
             for(E e : (Iterable<E>)getBoundPeers()) {
                 Point mid = e.midpoint();
@@ -430,8 +453,11 @@ public class Mesh implements Iterable<Mesh.T> {
                 e.p1.getE(mid).pair.     bindTo(e.getBindingMatrix(firste), first.pair);
                 e.p2.getE(mid).pair.     bindTo(e.getBindingMatrix(firste), firstq);
                 e.p2.getE(mid).pair.pair.bindTo(e.getBindingMatrix(firste), firstq.pair);
-                //e.p1.getE(mid).pair.prev.bindTo(e.pair.getBindingMatrix(firstx), first.pair.prev);
             }
+            /*
+            first.setConstraint(firste.getConstraint());
+            firstq.setConstraint(firste.getConstraint());
+            */
             return null;
         }
 
@@ -457,10 +483,6 @@ public class Mesh implements Iterable<Mesh.T> {
             pair.next.t = null;
             pair.prev.t = null;
 
-            this.bind_to = null;
-            pair.bind_to = null;
-            this.bind_peers = null;
-            pair.bind_peers = null;
             pair.prev.next = next;
             next.prev = pair.prev;
             prev.next = pair.next;
@@ -473,7 +495,6 @@ public class Mesh implements Iterable<Mesh.T> {
             this.prev.next = this;
             this.next.prev = this;
             this.pair.pair = this;
-            bind_peers.add(this);
             if (this.next.p1 != p2) throw new Error();
             if (this.prev.p2 != p1) throw new Error();
             if (this.p1.e == null) this.p1.e = this;
@@ -486,7 +507,12 @@ public class Mesh implements Iterable<Mesh.T> {
         public double dihedralAngle() {
             Vec v1 = t.norm().times(-1);
             Vec v2 = pair.t.norm().times(-1);
-            return Math.acos(v1.norm().dot(v2.norm()));
+            double prod = v1.norm().dot(v2.norm());
+            prod = Math.min(1,prod);
+            prod = Math.max(-1,prod);
+            double ret = Math.acos(prod);
+            if (Double.isNaN(ret)) throw new Error("nan! " + prod);
+            return ret;
         }
 
         /** angle between this half-edge and the next */
@@ -499,7 +525,7 @@ public class Mesh implements Iterable<Mesh.T> {
         public void makeAdjacent(E e) {
             if (this.next == e) return;
             if (p2 != e.p1) throw new Error("cannot make adjacent -- no shared vertex");
-            if (t != null || e.t != null) throw new Error("cannot make adjacent -- edges not both free");
+            if (t != null || e.t != null) throw new Error("cannot make adjacent -- edges not both free ");
 
             E freeIncident = p2.getFreeIncident(e, this);
 
@@ -580,7 +606,22 @@ public class Mesh implements Iterable<Mesh.T> {
         if (v2 != null) return new E(v2.getFreeIncident(), p1).pair;
         return new E(p1, p2);
     }
+    public boolean coalesce = false;
+    private static float round(float f) {
+        return Math.round(f*1000)/1000f;
+    }
     public T newT(Point p1, Point p2, Point p3, Vec norm, int colorclass) {
+        if (coalesce) {
+
+            for(Vertex v : vertices) { if (p1.distance(v.p) < EPSILON) { p1 = v.p; break; } }
+            for(Vertex v : vertices) { if (p2.distance(v.p) < EPSILON) { p2 = v.p; break; } }
+            for(Vertex v : vertices) { if (p3.distance(v.p) < EPSILON) { p3 = v.p; break; } }
+            /*
+            p1 = new Point(round(p1.x), round(p1.y), round(p1.z));
+            p2 = new Point(round(p2.x), round(p2.y), round(p2.z));
+            p3 = new Point(round(p3.x), round(p3.y), round(p3.z));
+            */
+        }
         if (norm != null) {
             Vec norm2 = p3.minus(p1).cross(p2.minus(p1));
             float dot = norm.dot(norm2);
@@ -651,9 +692,11 @@ public class Mesh implements Iterable<Mesh.T> {
         public void reinsert() { triangles.remove(this); triangles.add(this); }
 
         public boolean shouldBeDrawn() {
+
             if (e1().bindingGroupSize() <= 1) return false;
             if (e2().bindingGroupSize() <= 1) return false;
             if (e3().bindingGroupSize() <= 1) return false;
+
             return true;
         }