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
  16 // FIXME: re-orient goal (how?)
  17
  18 public class Main extends MeshViewer {
  19
  20     public static int verts = 0;
  21
  22     public static final Random random = new Random();
  23
  24     /** magnification factor */
  25     private static final float MAG = 1;
  26
  27     public Main(StlFile stlf, Frame f) {
  28         super(f);
  29
  30         for(int i=0; i<stlf.coordArray.length; i+=3) {
  31             Point p0 = new Point(stlf.coordArray[i+0].x * MAG, stlf.coordArray[i+0].y * MAG, stlf.coordArray[i+0].z * MAG);
  32             Point p1 = new Point(stlf.coordArray[i+1].x * MAG, stlf.coordArray[i+1].y * MAG, stlf.coordArray[i+1].z * MAG);
  33             Point p2 = new Point(stlf.coordArray[i+2].x * MAG, stlf.coordArray[i+2].y * MAG, stlf.coordArray[i+2].z * MAG);
  34             Vec n  = new Vec(stlf.normArray[i/3].x * MAG, stlf.normArray[i/3].y  * MAG, stlf.normArray[i/3].z * MAG);
  35             Mesh.T t  = goal.newT(p0, p1, p2, n, 0);
  36         }
  37
  38         // rotate to align major axis -- this probably needs to be done by a human.
  39         goal.transform(new Matrix(new Vec(0, 0, 1), (float)(Math.PI/2)));
  40
  41         float goal_width  = goal.diagonal().dot(new Vec(1, 0, 0));
  42         float goal_height = goal.diagonal().dot(new Vec(0, 1, 0));
  43         float goal_depth  = goal.diagonal().dot(new Vec(0, 0, 1));
  44
  45         /*
  46         float width  = (float)0.6;
  47         float height = (float)0.08;
  48         float depth  = (float)0.3;
  49         */
  50         float width  = (float)0.6;
  51         float depth  = (float)0.08;
  52         float height = (float)0.3;
  53
  54         float rshift = width/2;
  55         float lshift = -(width/2);
  56
  57         //float halfup = height/2;
  58         float halfup = 0;
  59
  60         translations = new Matrix[] {
  61
  62             new Matrix(new Vec(lshift,  depth,    halfup)),
  63             new Matrix(new Vec(rshift,  depth,    halfup)),
  64             new Matrix(new Vec(lshift, -depth,    halfup)),
  65             new Matrix(new Vec(rshift, -depth,    halfup)),
  66             /*
  67             new Matrix(new Vec(lshift,  depth,   -halfup)),
  68             new Matrix(new Vec(rshift,  depth,   -halfup)),
  69             new Matrix(new Vec(lshift, -depth,   -halfup)),
  70             new Matrix(new Vec(rshift, -depth,   -halfup)),
  71             */
  72
  73             new Matrix(new Vec(lshift,       0,  height)),
  74             new Matrix(new Vec(rshift,       0,  height)),
  75             new Matrix(new Vec(lshift,       0, -height)),
  76             new Matrix(new Vec(rshift,       0, -height)),
  77
  78             new Matrix(new Vec( width,           0,    0)),
  79             new Matrix(new Vec(-width,           0,    0)),
  80         };
  81
  82
  83         Point ltf = new Point(lshift,  (depth/2),  (height/2));
  84         Point mtf = new Point( 0.0,    (depth/2),  (height/2));
  85         Point rtf = new Point(rshift,  (depth/2),  (height/2));
  86         Point lbf = new Point(lshift, -(depth/2),  (height/2));
  87         Point mbf = new Point( 0.0,   -(depth/2),  (height/2));
  88         Point rbf = new Point(rshift, -(depth/2),  (height/2));
  89
  90         Point ltc = new Point(lshift,  (depth/2), 0);
  91         Point mtc = new Point( 0.0,    (depth/2), 0);
  92         Point rtc = new Point(rshift,  (depth/2), 0);
  93         Point lbc = new Point(lshift, -(depth/2), 0);
  94         Point mbc = new Point( 0.0,   -(depth/2), 0);
  95         Point rbc = new Point(rshift, -(depth/2), 0);
  96
  97         Point ltn = new Point(lshift,  (depth/2), -(height/2));
  98         Point mtn = new Point( 0.0,    (depth/2), -(height/2));
  99         Point rtn = new Point(rshift,  (depth/2), -(height/2));
 100         Point lbn = new Point(lshift, -(depth/2), -(height/2));
 101         Point mbn = new Point( 0.0,   -(depth/2), -(height/2));
 102         Point rbn = new Point(rshift, -(depth/2), -(height/2));
 103
 104
 105         Point[] points = new Point[] {
 106             ltf,
 107             mtf,
 108             rtf,
 109             lbf,
 110             mbf,
 111             rbf,
 112
 113             ltc,
 114             mtc,
 115             rtc,
 116             lbc,
 117             mbc,
 118             rbc,
 119
 120             ltn,
 121             mtn,
 122             rtn,
 123             lbn,
 124             mbn,
 125             rbn
 126         };
 127
 128
 129         // top
 130         tile.newT(ltf, mtf, mtc, null, 1);
 131         tile.newT(mtc, ltc, ltf, null, 1);
 132         tile.newT(mtf, rtf, rtc, null, 1);
 133         tile.newT(rtc, mtc, mtf, null, 1);
 134
 135         tile.newT(ltc, mtc, mtn, null, 1);
 136         tile.newT(mtn, ltn, ltc, null, 1);
 137         tile.newT(mtc, rtc, rtn, null, 1);
 138         tile.newT(rtn, mtn, mtc, null, 1);
 139
 140         // bottom (swap normals)
 141         tile.newT(mbf, lbf, mbc, null, 2);
 142         tile.newT(lbc, mbc, lbf, null, 2);
 143         tile.newT(rbf, mbf, rbc, null, 2);
 144         tile.newT(mbc, rbc, mbf, null, 2);
 145
 146         tile.newT(mbc, lbc, mbn, null, 2);
 147         tile.newT(lbn, mbn, lbc, null, 2);
 148
 149         tile.newT(rbc, mbc, rbn, null, 2);
 150         tile.newT(mbn, rbn, mbc, null, 2);
 151
 152
 153         // left
 154         tile.newT(ltf, ltc, lbc, null, 3);
 155         tile.newT(lbc, lbf, ltf, null, 3);
 156         tile.newT(ltc, ltn, lbn, null, 3);
 157         tile.newT(lbn, lbc, ltc, null, 3);
 158
 159         // right (swap normals)
 160         tile.newT(rtc, rtf, rbc, null, 4);
 161         tile.newT(rbf, rbc, rtf, null, 4);
 162         tile.newT(rtn, rtc, rbn, null, 4);
 163         tile.newT(rbc, rbn, rtc, null, 4);
 164
 165         // front
 166         tile.newT(ltn, mtn, mbn, null, 5);
 167         tile.newT(ltn, mbn, lbn, null, 5);
 168         tile.newT(mtn, rtn, rbn, null, 5);
 169         tile.newT(mtn, rbn, mbn, null, 5);
 170
 171         // back
 172         tile.newT(mtf, ltf, mbf, null, 6);
 173         tile.newT(mbf, ltf, lbf, null, 6);
 174         tile.newT(rtf, mtf, rbf, null, 6);
 175         tile.newT(rbf, mtf, mbf, null, 6);
 176
 177         for(Matrix m : translations) {
 178             for(Mesh.T t1 : tile) {
 179                 for(Mesh.T t2 : tile) {
 180                     if (t1==t2) continue;
 181
 182                     if ((t1.v1().p.times(m).minus(t2.v1().p).mag() < Mesh.EPSILON) &&
 183                         (t1.v2().p.times(m).minus(t2.v3().p).mag() < Mesh.EPSILON) &&
 184                         (t1.v3().p.times(m).minus(t2.v2().p).mag() < Mesh.EPSILON)) {
 185                         t1.e1().bindEdge(t2.e3());
 186                         t1.e2().bindEdge(t2.e2());
 187                         t1.e3().bindEdge(t2.e1());
 188                     }
 189                     if ((t1.v2().p.times(m).minus(t2.v1().p).mag() < Mesh.EPSILON) &&
 190                         (t1.v3().p.times(m).minus(t2.v3().p).mag() < Mesh.EPSILON) &&
 191                         (t1.v1().p.times(m).minus(t2.v2().p).mag() < Mesh.EPSILON)) {
 192                         t1.e2().bindEdge(t2.e3());
 193                         t1.e3().bindEdge(t2.e2());
 194                         t1.e1().bindEdge(t2.e1());
 195                     }
 196                     if ((t1.v3().p.times(m).minus(t2.v1().p).mag() < Mesh.EPSILON) &&
 197                         (t1.v1().p.times(m).minus(t2.v3().p).mag() < Mesh.EPSILON) &&
 198                         (t1.v2().p.times(m).minus(t2.v2().p).mag() < Mesh.EPSILON)) {
 199                         t1.e3().bindEdge(t2.e3());
 200                         t1.e1().bindEdge(t2.e2());
 201                         t1.e2().bindEdge(t2.e1());
 202                     }
 203
 204                 }
 205             }
 206         }
 207
 208         //xMesh.Vert mid = lbf.getE(mbn).shatter();
 209
 210         // rescale to match volume
 211         float factor = (float)Math.pow(tile.volume() / goal.volume(), 1.0/3.0);
 212         goal.transform(new Matrix(factor));
 213
 214         // translate to match centroid
 215         goal.transform(new Matrix(tile.centroid().minus(goal.centroid())));
 216
 217         //tx.e2.shatter();
 218         //tx.e3.shatter();
 219
 220
 221         tile.rebindPoints();
 222
 223         //mid.move(new Vec((float)0,0,(float)-0.05));
 224         //ltn.move(new Vec((float)0,0,(float)-0.05));
 225
 226         //mtf.move(new Vec(0, (float)-0.05, (float)0.05));
 227
 228
 229         System.out.println("tile volume: " + tile.volume());
 230         System.out.println("goal volume: " + goal.volume());
 231
 232         tile.score_against = goal;
 233         goal.score_against = tile;
 234     }
 235
 236     public synchronized void breakit() {
 237         if (verts > 800) return;
 238         //while(verts < 800) {
 239         PriorityQueue<Mesh.E> es = new PriorityQueue<Mesh.E>();
 240         for(Mesh.E e : tile.edges()) es.add(e);
 241         for(int i=0; i<8; i++) {
 242             Mesh.E e = es.poll();
 243             verts++;
 244             System.out.println("shatter " + e);
 245             e.shatter();
 246         }
 247         tile.rebindPoints();
 248         //}
 249     }
 250
 251     public synchronized void rand(float temp, Mesh.Vert p) {
 252         double tile_score = tile.score();
 253         double goal_score = goal.score();
 254
 255         //p.reComputeError();
 256         p.reComputeErrorAround();
 257
 258         Vec v = new Vec((random.nextFloat() - (float)0.5) / 400,
 259                         (random.nextFloat() - (float)0.5) / 400,
 260                         (random.nextFloat() - (float)0.5) / 400);
 261         /*
 262         Matrix inv = p.errorQuadric();
 263         Vec v = new Vec(inv.d, inv.h, inv.l).norm().times(1/(float)1000);
 264         */
 265         boolean good = p.move(v);
 266
 267         p.reComputeErrorAround();
 268
 269         double new_tile_score = tile.score();
 270         double new_goal_score = goal.score();
 271         double tile_delta = (new_tile_score - tile_score) / tile_score;
 272         double goal_delta = (new_goal_score - goal_score) / goal_score;
 273         double delta = tile_delta + goal_delta;
 274         double swapProbability = Math.exp((-1 * delta) / temp);
 275         boolean doSwap = good && (Math.random() < swapProbability);
 276         //boolean doSwap = good && (tile_delta <= 0 && goal_delta <= 0);
 277         if (doSwap) {
 278             tile_score = new_tile_score;
 279             goal_score = new_goal_score;
 280             //System.out.println("score: " + tile_score + " / " + goal_score);
 281             hits++;
 282         } else {
 283             p.move(v.times(-1));
 284             misses++;
 285         }
 286     }
 287
 288     float hits = 0;
 289     float misses = 0;
 290     public void anneal() throws Exception {
 291         int verts = 0;
 292         float hightemp = 10;
 293         float temp = hightemp;
 294         float last = 10;
 295         while(true) {
 296             double ratio = (hits+misses==0) ? 1 : (hits / (hits+misses));
 297             System.out.println("temp="+temp + " ratio="+(Math.ceil(ratio*100)));
 298             hits = 0;
 299             misses = 0;
 300             float gamma = 0;
 301             double acceptance = ratio;
 302             if (breaks) {
 303                 breaks = false;
 304                     breakit();
 305                     gamma = 2;
 306                     //temp = last * 0.8f;
 307                     //last = temp;
 308                     //temp = hightemp;
 309             } else
 310             if (acceptance > 0.96) gamma = 0.6f;
 311             else if (acceptance > 0.9) gamma = 0.7f;
 312             else if (acceptance > 0.8) gamma = 0.75f;
 313             else if (acceptance > 0.6) gamma = 0.8f;
 314             else {
 315                 if (acceptance > 0.3) {
 316                     gamma = 0.9f;
 317                 } else if (acceptance > 0.15) {
 318                     gamma = 0.95f;
 319                 } else {
 320                     breakit();
 321                     gamma = 2;
 322                     //temp = last * 0.8f;
 323                     //last = temp;
 324                     //temp = hightemp;
 325                 }
 326             }
 327             temp = temp * gamma;
 328
 329             HashSet<Mesh.Vert> hs = new HashSet<Mesh.Vert>();
 330             for(Mesh.Vert p : tile.vertices()) hs.add(p);
 331             for(int i=0; i<10; i++) {
 332                 repaint();
 333                 for(Mesh.Vert v : hs) rand(temp,v);
 334             }
 335             tile.rebuildPointSet();
 336             repaint();
 337             //breakit();
 338             repaint();
 339             goal.unApplyQuadricToNeighborAll();
 340             repaint();
 341             tile.recomputeAllFundamentalQuadrics();
 342             repaint();
 343             goal.applyQuadricToNeighborAll();
 344        }
 345     }
 346
 347     public static void main(String[] s) throws Exception {
 348         StlFile stlf = new StlFile();
 349         stlf.load("simplefish.stl");
 350         Frame f = new Frame();
 351         Main main = new Main(stlf, f);
 352         f.pack();
 353         f.show();
 354         f.setSize(900, 900);
 355         f.doLayout();
 356         main.anneal();
 357     }
 358
 359 }