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 = 0;
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         goal.ignorecollision = true;
67
68         // rotate to align major axis -- this probably needs to be done by a human.
69         goal.transform(new Matrix(new Vec(0, 0, 1), (float)(Math.PI/2)));
70
71         float goal_width  = goal.diagonal().dot(new Vec(1, 0, 0));
72         float goal_height = goal.diagonal().dot(new Vec(0, 1, 0));
73         float goal_depth  = goal.diagonal().dot(new Vec(0, 0, 1));
74
75         /*
76         float width  = (float)0.6;
77         float height = (float)0.08;
78         float depth  = (float)0.3;
79         */
80         float width  = (float)0.7;
81         float depth  = (float)0.08;
82         float height = (float)0.4;
83
84         float rshift =   width/2;
85         float lshift = -(width/2);
86
87         //float halfup = height/2;
88         float halfup = 0;
89
90         translations = new Matrix[] {
91             /*
92             new Matrix(new Vec(lshift,  depth,    halfup)),
93             new Matrix(new Vec(rshift,  depth,    halfup)),
94             new Matrix(new Vec(lshift, -depth,    halfup)),
95             new Matrix(new Vec(rshift, -depth,    halfup)),
96             */
97
98             new Matrix(new Vec(0,  depth,    halfup)),
99             new Matrix(new Vec(0, -depth,    halfup)),
100
101
102             new Matrix(new Vec(lshift,       0,  height)),
103             new Matrix(new Vec(rshift,       0,  height)),
104             new Matrix(new Vec(lshift,       0, -height)),
105             new Matrix(new Vec(rshift,       0, -height)),
106
107
108             new Matrix(new Vec( width,           0,    0)),
109             new Matrix(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                         t1.e1().bindEdge(t2.e3());
220                         t1.e2().bindEdge(t2.e2());
221                         t1.e3().bindEdge(t2.e1());
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                         t1.e2().bindEdge(t2.e3());
227                         t1.e3().bindEdge(t2.e2());
228                         t1.e1().bindEdge(t2.e1());
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                         t1.e3().bindEdge(t2.e3());
234                         t1.e1().bindEdge(t2.e2());
235                         t1.e2().bindEdge(t2.e1());
236                     }
237
238                 }
239             }
240         }
241
242         //xMesh.Vert 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(new Matrix(factor));
247
248         // translate to match centroid
249         goal.transform(new Matrix(tile.centroid().minus(goal.centroid())));
250
251         //tx.e2.shatter();
252         //tx.e3.shatter();
253
254
255         tile.rebindPoints();
256
257         //mid.move(new Vec((float)0,0,(float)-0.05));
258         //ltn.move(new Vec((float)0,0,(float)-0.05));
259
260         //mtf.move(new Vec(0, (float)-0.05, (float)0.05));
261
262
263         System.out.println("tile volume: " + tile.volume());
264         System.out.println("goal volume: " + goal.volume());
265
266         tile.score_against = goal;
267         goal.score_against = tile;
268         tile.tilemesh = true;
269     }
270
271     public synchronized void breakit() {
272         if (verts > 800) return;
273         //while(verts < 800) {
274         PriorityQueue<Mesh.E> es = new PriorityQueue<Mesh.E>();
275         for(Mesh.E e : tile.edges()) es.add(e);
276         for(int i=0; i<10; i++) {
277             Mesh.E e = es.poll();
278             verts++;
279             System.out.println("shatter " + e);
280             e.shatter();
281         }
282         tile.rebindPoints();
283         //}
284     }
285
286     public synchronized void rand(float temp, Mesh.Vert p) {
287
288         //p.reComputeError();
289         p.reComputeErrorAround();
290         double tile_score = tile.score();
291         double goal_score = goal.score();
292
293         Vec v;
294         /*
295         Matrix inv = p.errorQuadric();
296         v = new Vec(inv.d, inv.h, inv.l).norm().times(1/(float)300);
297         if (p.quadric_count == 0) {
298             v = goal.nearest(p.p).p.minus(p.p).norm().times(1/(float)300);
299         }
300         */
301         Vec v2 = new Vec((random.nextFloat() - (float)0.5) / 500,
302                         (random.nextFloat() - (float)0.5)  / 500,
303                         (random.nextFloat() - (float)0.5)  / 500);
304         //v = v.plus(v2.norm().times(1/(float)300));
305         v = v2.norm().times(1/(float)300);
306
307         boolean good = p.move(v);
308
309         p.reComputeErrorAround();
310
311         double new_tile_score = tile.score();
312         double new_goal_score = goal.score();
313         double tile_delta = (new_tile_score - tile_score) / tile_score;
314         double goal_delta = (new_goal_score - goal_score) / goal_score;
315         double delta = tile_delta + goal_delta;
316         double swapProbability = Math.exp((-1 * delta) / temp);
317         boolean doSwap = good && (Math.random() < swapProbability);
318         //boolean doSwap = good && (tile_delta <= 0 && goal_delta <= 0);
319         //boolean doSwap = good && (tile_delta + goal_delta <= 0);
320         if (doSwap) {
321             tile_score = new_tile_score;
322             goal_score = new_goal_score;
323             //System.out.println("score: " + tile_score + " / " + goal_score);
324             hits++;
325         } else {
326             p.move(v.times(-1));
327             misses++;
328         }
329     }
330
331     float hits = 0;
332     float misses = 0;
333     public void anneal() throws Exception {
334         int verts = 0;
335         float hightemp = 10;
336         float temp = hightemp;
337         float last = 10;
338         while(true) {
339             synchronized(this) {
340             double ratio = (hits+misses==0) ? 1 : (hits / (hits+misses));
341             System.out.println("temp="+temp + " ratio="+(Math.ceil(ratio*100)));
342             hits = 0;
343             misses = 0;
344             float gamma = 0;
345             double acceptance = ratio;
346             accepts = (int)(Math.ceil(ratio*100));
347             temps = (int)(Math.ceil(temp*1000));
348             vertss = tile.size();
349             if (breaks) {
350                 breaks = false;
351                     breakit();
352                     //gamma = 1;
353                     gamma = 1;
354                     //temp = last * 0.8f;
355                     //last = temp;
356                     //temp = hightemp;
357             } else
358             if (acceptance > 0.96) gamma = 0.4f;
359             else if (acceptance > 0.9) gamma = 0.5f;
360             else if (acceptance > 0.8) gamma = 0.65f;
361             else if (acceptance > 0.6) gamma = 0.7f;
362             else {
363                 if (acceptance > 0.3) {
364                     gamma = 0.9f;
365                 } else if (acceptance > 0.15) {
366                     gamma = 0.95f;
367                 } else {
368                     breakit();
369                     //gamma = 1;
370                     gamma = 0.99f;
371                     //gamma = 1;
372                     //temp = last * 0.8f;
373                     //last = temp;
374                     //temp = hightemp;
375                 }
376             }
377             temp = temp * gamma;
378
379
380             HashSet<Mesh.Vert> hs = new HashSet<Mesh.Vert>();
381             for(Mesh.Vert p : tile.vertices()) hs.add(p);
382             for(int i=0; i<10; i++) {
383                 repaint();
384                 for(Mesh.Vert v : hs) {
385                     if (anneal) rand(temp,v);
386                     Thread.yield();
387                     repaint();
388                 }
389             }
390             tile.rebuildPointSet();
391             repaint();
392             //breakit();
393             repaint();
394             goal.unApplyQuadricToNeighborAll();
395             repaint();
396             tile.recomputeAllFundamentalQuadrics();
397             repaint();
398             goal.applyQuadricToNeighborAll();
399             }
400        }
401     }
402
403     public static void main(String[] s) throws Exception {
404         StlFile stlf = new StlFile();
405         stlf.load("fish.stl");
406         //stlf.load("monkey.stl");
407         Frame f = new Frame();
408         Main main = new Main(stlf, f);
409         f.pack();
410         f.show();
411         f.setSize(900, 900);
412         f.doLayout();
413         main.anneal();
414     }
415
416 }