checkpoint
[anneal.git] / src / edu / berkeley / qfat / Main.java
1 package edu.berkeley.qfat;
2 import java.awt.*;
3 import java.awt.event.*;
4 import javax.swing.*;
5 import javax.media.opengl.*;
6 import javax.media.opengl.glu.*;
7 import java.util.*;
8 import edu.berkeley.qfat.geom.*;
9 import edu.berkeley.qfat.geom.Point;
10
11 // TO DO:
12 // - real anneal
13 // - solve self-intersection problem
14 // - get a better test model?
15 // - symmetry constraints withing the tile
16 // - rotation matrices
17 // - overbinding results in forced equational constraints on the leader
18 // - shatter in invertd-triforce pattern brian mentioned
19 // - aspect ratio?  non-uniform deformation?
20 // - rotational alignment
21
22 // - movie-style user interface like
23 //      http://www.coleran.com/markcoleranreell.html ?
24
25 // - consider recasting the Shewchuk predicates in Java?
26 //    http://www.cs.cmu.edu/afs/cs/project/quake/public/code/predicates.c
27
28 /*
29   blender keys
30   - middle mouse = option+click
31   - right mouse = command+click
32
33   3,7,1 = view along axes (control for opp direction)
34   4, 8, 7, 2 = rotate in discrete increments (+control to translate)
35   middle trag: rotate space
36   shift+middle drag: translate space
37   wheel: zoom
38   home: home view: take current angle, zoom to whole scnee
39   5 = ortho vs non-ortho
40
41 */
42
43
44 // FIXME: re-orient goal (how?)
45
46 public class Main extends MeshViewer {
47
48     public static int verts = 1;
49
50     public static final Random random = new Random();
51     
52     /** magnification factor */
53     private static final float MAG = 1;
54
55     public Main(StlFile stlf, Frame f) {
56         super(f);
57
58         for(int i=0; i<stlf.coordArray.length; i+=3) {
59             Point p0 = new Point(stlf.coordArray[i+0].x * MAG, stlf.coordArray[i+0].y * MAG, stlf.coordArray[i+0].z * MAG);
60             Point p1 = new Point(stlf.coordArray[i+1].x * MAG, stlf.coordArray[i+1].y * MAG, stlf.coordArray[i+1].z * MAG);
61             Point p2 = new Point(stlf.coordArray[i+2].x * MAG, stlf.coordArray[i+2].y * MAG, stlf.coordArray[i+2].z * MAG);
62             Vec n  = new Vec(stlf.normArray[i/3].x * MAG, stlf.normArray[i/3].y  * MAG, stlf.normArray[i/3].z * MAG);
63             Mesh.T t  = goal.newT(p0, p1, p2, n, 0);
64         }
65
66         // rotate to align major axis -- this probably needs to be done by a human.
67         goal.transform(Matrix.rotate(new Vec(0, 0, 1), (float)(Math.PI/2)));
68
69         float goal_width  = goal.diagonal().dot(new Vec(1, 0, 0));
70         float goal_height = goal.diagonal().dot(new Vec(0, 1, 0));
71         float goal_depth  = goal.diagonal().dot(new Vec(0, 0, 1));
72
73         /*
74           float width  = (float)0.6;
75           float height = (float)0.08;
76           float depth  = (float)0.3;
77         */
78         float width  = (float)0.7;
79         float depth  = (float)0.08;
80         float height = (float)0.4;
81
82         float rshift =   width/2;
83         float lshift = -(width/2);
84
85         //float halfup = height/2;
86         float halfup = 0;
87
88         translations = new Matrix[] {
89
90             Matrix.translate(new Vec(lshift,  depth,    halfup)),
91             Matrix.translate(new Vec(rshift,  depth,    halfup)),
92             /*
93             Matrix.translate(new Vec(lshift, -depth,    halfup)),
94             Matrix.translate(new Vec(rshift, -depth,    halfup)),
95             */
96
97             /*
98               Matrix.translate(new Vec(0,  depth,    halfup)),
99               Matrix.translate(new Vec(0, -depth,    halfup)),
100             */
101             /*
102             Matrix.translate(new Vec(lshift,       0,  height)),
103             Matrix.translate(new Vec(rshift,       0,  height)),
104             Matrix.translate(new Vec(lshift,       0, -height)),
105             Matrix.translate(new Vec(rshift,       0, -height)),
106
107
108             Matrix.translate(new Vec( width,           0,    0)),
109             Matrix.translate(new Vec(-width,           0,    0)),
110             */
111         };
112
113         //   
114
115
116
117         Point ltf = new Point(lshift,  (depth/2),  (height/2));
118         Point mtf = new Point( 0.0,    (depth/2),  (height/2));
119         Point rtf = new Point(rshift,  (depth/2),  (height/2));
120         Point lbf = new Point(lshift, -(depth/2),  (height/2));
121         Point mbf = new Point( 0.0,   -(depth/2),  (height/2));
122         Point rbf = new Point(rshift, -(depth/2),  (height/2));
123
124         Point ltc = new Point(lshift,  (depth/2), 0);
125         Point mtc = new Point( 0.0,    (depth/2), 0);
126         Point rtc = new Point(rshift,  (depth/2), 0);
127         Point lbc = new Point(lshift, -(depth/2), 0);
128         Point mbc = new Point( 0.0,   -(depth/2), 0);
129         Point rbc = new Point(rshift, -(depth/2), 0);
130
131         Point ltn = new Point(lshift,  (depth/2), -(height/2));
132         Point mtn = new Point( 0.0,    (depth/2), -(height/2));
133         Point rtn = new Point(rshift,  (depth/2), -(height/2));
134         Point lbn = new Point(lshift, -(depth/2), -(height/2));
135         Point mbn = new Point( 0.0,   -(depth/2), -(height/2));
136         Point rbn = new Point(rshift, -(depth/2), -(height/2));
137
138         
139         Point[] points = new Point[] {
140             ltf,
141             mtf,
142             rtf,
143             lbf,
144             mbf,
145             rbf,
146
147             ltc,
148             mtc,
149             rtc,
150             lbc,
151             mbc,
152             rbc,
153
154             ltn,
155             mtn,
156             rtn,
157             lbn,
158             mbn,
159             rbn
160         };
161
162
163         // top
164         tile.newT(ltf, mtf, mtc, null, 1);
165         tile.newT(mtc, ltc, ltf, null, 1);
166         tile.newT(mtf, rtf, rtc, null, 1);
167         tile.newT(rtc, mtc, mtf, null, 1);
168
169         tile.newT(ltc, mtc, mtn, null, 1);
170         tile.newT(mtn, ltn, ltc, null, 1);
171         tile.newT(mtc, rtc, rtn, null, 1);
172         tile.newT(rtn, mtn, mtc, null, 1);
173
174         // bottom (swap normals)
175         tile.newT(mbf, lbf, mbc, null, 2);
176         tile.newT(lbc, mbc, lbf, null, 2);
177         tile.newT(rbf, mbf, rbc, null, 2);
178         tile.newT(mbc, rbc, mbf, null, 2);
179
180         tile.newT(mbc, lbc, mbn, null, 2);
181         tile.newT(lbn, mbn, lbc, null, 2);
182
183         tile.newT(rbc, mbc, rbn, null, 2);
184         tile.newT(mbn, rbn, mbc, null, 2);
185
186
187         // left
188         tile.newT(ltf, ltc, lbc, null, 3);
189         tile.newT(lbc, lbf, ltf, null, 3);
190         tile.newT(ltc, ltn, lbn, null, 3);
191         tile.newT(lbn, lbc, ltc, null, 3);
192
193         // right (swap normals)
194         tile.newT(rtc, rtf, rbc, null, 4);
195         tile.newT(rbf, rbc, rtf, null, 4);
196         tile.newT(rtn, rtc, rbn, null, 4);
197         tile.newT(rbc, rbn, rtc, null, 4);
198
199         // front
200         tile.newT(ltn, mtn, mbn, null, 5);
201         tile.newT(ltn, mbn, lbn, null, 5);
202         tile.newT(mtn, rtn, rbn, null, 5);
203         tile.newT(mtn, rbn, mbn, null, 5);
204
205         // back
206         tile.newT(mtf, ltf, mbf, null, 6);
207         tile.newT(mbf, ltf, lbf, null, 6);
208         tile.newT(rtf, mtf, rbf, null, 6);
209         tile.newT(rbf, mtf, mbf, null, 6);
210
211         for(Matrix m : translations) {
212             for(Mesh.T t1 : tile) {
213                 for(Mesh.T t2 : tile) {
214                     if (t1==t2) continue;
215
216                     if ((t1.v1().p.times(m).minus(t2.v1().p).mag() < Mesh.EPSILON) &&
217                         (t1.v2().p.times(m).minus(t2.v3().p).mag() < Mesh.EPSILON) &&
218                         (t1.v3().p.times(m).minus(t2.v2().p).mag() < Mesh.EPSILON)) {
219                         t2.e3().bindEdge(t1.e1(), m);
220                         t2.e2().bindEdge(t1.e2(), m);
221                         t2.e1().bindEdge(t1.e3(), m);
222                     }
223                     if ((t1.v2().p.times(m).minus(t2.v1().p).mag() < Mesh.EPSILON) &&
224                         (t1.v3().p.times(m).minus(t2.v3().p).mag() < Mesh.EPSILON) &&
225                         (t1.v1().p.times(m).minus(t2.v2().p).mag() < Mesh.EPSILON)) {
226                         t2.e3().bindEdge(t1.e2(), m);
227                         t2.e2().bindEdge(t1.e3(), m);
228                         t2.e1().bindEdge(t1.e1(), m);
229                     }
230                     if ((t1.v3().p.times(m).minus(t2.v1().p).mag() < Mesh.EPSILON) &&
231                         (t1.v1().p.times(m).minus(t2.v3().p).mag() < Mesh.EPSILON) &&
232                         (t1.v2().p.times(m).minus(t2.v2().p).mag() < Mesh.EPSILON)) {
233                         t2.e3().bindEdge(t1.e3(), m);
234                         t2.e2().bindEdge(t1.e1(), m);
235                         t2.e1().bindEdge(t1.e2(), m);
236                     }
237
238                 }
239             }
240         }
241
242         //xMesh.Vertex mid = lbf.getE(mbn).shatter();
243
244         // rescale to match volume
245         float factor = (float)Math.pow(tile.volume() / goal.volume(), 1.0/3.0);
246         goal.transform(Matrix.scale(factor));
247
248         // translate to match centroid
249         goal.transform(Matrix.translate(tile.centroid().minus(goal.centroid())));
250         goal.makeVerticesImmutable();
251
252         //tx.e2.shatter();
253         //tx.e3.shatter();
254
255
256         tile.rebindPoints();
257
258         //mid.move(new Vec((float)0,0,(float)-0.05));
259         //ltn.move(new Vec((float)0,0,(float)-0.05));
260
261         //mtf.move(new Vec(0, (float)-0.05, (float)0.05));
262
263
264         System.out.println("tile volume: " + tile.volume());
265         System.out.println("goal volume: " + goal.volume());
266
267         tile.error_against = goal;
268         goal.error_against = tile;
269     }
270
271     public void breakit() {
272         int oldverts = verts;
273         System.out.println("doubling vertices.");
274         PriorityQueue<Mesh.E> es = new PriorityQueue<Mesh.E>();
275         for(Mesh.T t : tile) {
276             es.add(t.e1());
277             es.add(t.e2());
278             es.add(t.e3());
279             Thread.yield();
280             repaint();
281         }
282         for(int i=0; i<Math.min(oldverts,200); i++) {
283             Mesh.E e = es.poll();
284             verts++;
285             //System.out.println("shatter " + e);
286             //e.shatter(e.midpoint(), null, null, true, true);
287             
288             //e.shatter(e.midpoint(), null, null, true, false);
289             Thread.yield();
290             repaint();
291         }
292         tile.rebindPoints();
293     }
294
295     public synchronized void rand(float temp, Mesh.Vertex p) {
296
297         p.reComputeErrorAround();
298         double tile_error = tile.error();
299         double goal_error = goal.error();
300
301         float max = p.averageEdgeLength()/10;
302         Vec v = new Vec(random.nextFloat(), random.nextFloat(), random.nextFloat());
303         v = v.norm().times((random.nextFloat() - 0.5f) * max);
304         //System.out.println(max + " " + p.averageEdgeLength() + " " + v.mag());
305         Matrix m = Matrix.translate(v);
306
307         boolean good = p.move(m, false);
308         if (!good) { /*misses++;*/ return; }
309
310         double new_tile_error = tile.error();
311         double new_goal_error = goal.error();
312         double tile_delta = (new_tile_error - tile_error) / tile_error;
313         double goal_delta = (new_goal_error - goal_error) / goal_error;
314         double delta = tile_delta + goal_delta;
315         double swapProbability = Math.exp((-1 * delta) / (((double)temp)/1000000));
316         boolean doSwap = good && (Math.random() < swapProbability);
317         //boolean doSwap = good && (tile_delta <= 0 && goal_delta <= 0);
318         //boolean doSwap = good && (tile_delta + goal_delta <= 0);
319         if (doSwap) {
320             tile_error = new_tile_error;
321             goal_error = new_goal_error;
322             //System.out.println("error: " + tile_error + " / " + goal_error);
323             hits++;
324             p.goodp = p.p;
325         } else {
326             p.move(Matrix.translate(v.times(-1)), true);
327             misses++;
328         }
329     }
330
331     float hits = 0;
332     float misses = 0;
333     public void anneal() throws Exception {
334         float hightemp = 1;
335         float temp = hightemp;
336         float last = 10;
337         boolean seek_upward = false;
338         double acceptance = 1;
339         while(true) {
340             synchronized(this) {
341                 double ratio = (hits+misses==0) ? 1 : (hits / (hits+misses));
342                 hits = 0;
343                 misses = 0;
344                 float gamma = 1;
345                 acceptance = (ratio+acceptance)/2;
346                 accepts = (int)(Math.ceil(ratio*100));
347                 temps = (int)(Math.ceil(temp*1000));
348                 vertss = tile.size();
349                 if (breaks > 0) {
350                     while (breaks>0) {
351                         breaks--;
352                         breakit();
353                     }
354                     seek_upward = true;
355                 } else if (acceptance > 0.96) gamma = 0.4f;
356                 else if (acceptance > 0.9)    gamma = 0.5f;
357                 else if (acceptance > 0.8)    gamma = 0.65f;
358                 else if (acceptance > 0.6)    gamma = 0.7f;
359                 else if (acceptance > 0.3)    gamma = 0.8f;
360                 else if (acceptance > 0.15)   gamma = 0.9f;
361                 else if (acceptance > 0.05)   gamma = 0.95f;
362                 else if (acceptance > 0.01)   gamma = 0.98f;
363                 else { /*breaks++;*/ }
364
365                 if (seek_upward) {
366                     if (acceptance > 0.2) seek_upward = false;
367                     else gamma = 2-gamma;
368                 }
369
370                 temp = temp * gamma;
371
372
373                 HashSet<Mesh.Vertex> hs = new HashSet<Mesh.Vertex>();
374                 for(Mesh.Vertex p : tile.vertices()) hs.add(p);
375                 Mesh.Vertex[] pts = (Mesh.Vertex[])hs.toArray(new Mesh.Vertex[0]);
376
377                 int count = 0;
378                 long then = System.currentTimeMillis();
379                 for(int i=0; i<40; i++) {
380                     if (anneal) {
381                         count++;
382                         Mesh.Vertex v = pts[Math.abs(random.nextInt()) % pts.length];
383                         rand(temp,v);
384                         v.recomputeFundamentalQuadricIfStale();
385                         v.recomputeFundamentalQuadricIfNeighborChanged();
386                     }
387                     Thread.yield();
388                     repaint();
389                 }
390                 PriorityQueue<Mesh.E> es = new PriorityQueue<Mesh.E>();
391                 for(Mesh.T t : tile) {
392                     float max = 5;
393                     for(Mesh.E e : new Mesh.E[] { t.e1(), t.e2(), t.e3() }) {
394                         if (e==null) continue;
395                         if (e.stretchRatio() > max) es.add(e);
396                         if (t.aspect() < 0.1 && e.length()>e.next.length() && e.length()>e.prev.length()) es.add(e);
397                     }
398                 }
399                 /*
400                 for(int i=0; i<5; i++) {
401                     Mesh.E e = es.poll();
402                     if (e==null) break;
403                     e.shatter();
404                 }
405                 */
406                 tile.rebindPoints();
407
408                 System.out.println("temp="+temp + " ratio="+(Math.ceil(acceptance*100)) + " " +
409                                    "points_per_second=" +
410                                    (count*1000)/((double)(System.currentTimeMillis()-then)));
411
412                 for(Mesh.Vertex p : goal.vertices()) p.recomputeFundamentalQuadricIfNeighborChanged();
413
414                 synchronized(safeTriangles) {
415                     safeTriangles.clear();
416                     for(Mesh.T t : tile) 
417                         if (t.shouldBeDrawn())
418                             safeTriangles.add(t);
419                 }
420             }
421         }
422     }
423
424
425     public static void main(String[] s) throws Exception {
426         StlFile stlf = new StlFile();
427         stlf.load("fish.stl");
428         //stlf.load("monkey.stl");
429         Frame f = new Frame();
430         Main main = new Main(stlf, f);
431         f.pack();
432         f.show();
433         f.setSize(900, 900);
434         f.doLayout();
435         main.anneal();
436     }
437
438 }