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