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             Matrix.translate(new Vec(lshift, -depth,    halfup)),
  93             Matrix.translate(new Vec(rshift, -depth,    halfup)),
  94
  95             /*
  96               Matrix.translate(new Vec(0,  depth,    halfup)),
  97               Matrix.translate(new Vec(0, -depth,    halfup)),
  98             */
  99
 100             Matrix.translate(new Vec(lshift,       0,  height)),
 101             Matrix.translate(new Vec(rshift,       0,  height)),
 102             Matrix.translate(new Vec(lshift,       0, -height)),
 103             Matrix.translate(new Vec(rshift,       0, -height)),
 104
 105
 106             Matrix.translate(new Vec( width,           0,    0)),
 107             Matrix.translate(new Vec(-width,           0,    0)),
 108
 109         };
 110
 111         //
 112
 113
 114
 115         Point ltf = new Point(lshift,  (depth/2),  (height/2));
 116         Point mtf = new Point( 0.0,    (depth/2),  (height/2));
 117         Point rtf = new Point(rshift,  (depth/2),  (height/2));
 118         Point lbf = new Point(lshift, -(depth/2),  (height/2));
 119         Point mbf = new Point( 0.0,   -(depth/2),  (height/2));
 120         Point rbf = new Point(rshift, -(depth/2),  (height/2));
 121
 122         Point ltc = new Point(lshift,  (depth/2), 0);
 123         Point mtc = new Point( 0.0,    (depth/2), 0);
 124         Point rtc = new Point(rshift,  (depth/2), 0);
 125         Point lbc = new Point(lshift, -(depth/2), 0);
 126         Point mbc = new Point( 0.0,   -(depth/2), 0);
 127         Point rbc = new Point(rshift, -(depth/2), 0);
 128
 129         Point ltn = new Point(lshift,  (depth/2), -(height/2));
 130         Point mtn = new Point( 0.0,    (depth/2), -(height/2));
 131         Point rtn = new Point(rshift,  (depth/2), -(height/2));
 132         Point lbn = new Point(lshift, -(depth/2), -(height/2));
 133         Point mbn = new Point( 0.0,   -(depth/2), -(height/2));
 134         Point rbn = new Point(rshift, -(depth/2), -(height/2));
 135
 136
 137         Point[] points = new Point[] {
 138             ltf,
 139             mtf,
 140             rtf,
 141             lbf,
 142             mbf,
 143             rbf,
 144
 145             ltc,
 146             mtc,
 147             rtc,
 148             lbc,
 149             mbc,
 150             rbc,
 151
 152             ltn,
 153             mtn,
 154             rtn,
 155             lbn,
 156             mbn,
 157             rbn
 158         };
 159
 160
 161         // top
 162         tile.newT(ltf, mtf, mtc, null, 1);
 163         tile.newT(mtc, ltc, ltf, null, 1);
 164         tile.newT(mtf, rtf, rtc, null, 1);
 165         tile.newT(rtc, mtc, mtf, null, 1);
 166
 167         tile.newT(ltc, mtc, mtn, null, 1);
 168         tile.newT(mtn, ltn, ltc, null, 1);
 169         tile.newT(mtc, rtc, rtn, null, 1);
 170         tile.newT(rtn, mtn, mtc, null, 1);
 171
 172         // bottom (swap normals)
 173         tile.newT(mbf, lbf, mbc, null, 2);
 174         tile.newT(lbc, mbc, lbf, null, 2);
 175         tile.newT(rbf, mbf, rbc, null, 2);
 176         tile.newT(mbc, rbc, mbf, null, 2);
 177
 178         tile.newT(mbc, lbc, mbn, null, 2);
 179         tile.newT(lbn, mbn, lbc, null, 2);
 180
 181         tile.newT(rbc, mbc, rbn, null, 2);
 182         tile.newT(mbn, rbn, mbc, null, 2);
 183
 184
 185         // left
 186         tile.newT(ltf, ltc, lbc, null, 3);
 187         tile.newT(lbc, lbf, ltf, null, 3);
 188         tile.newT(ltc, ltn, lbn, null, 3);
 189         tile.newT(lbn, lbc, ltc, null, 3);
 190
 191         // right (swap normals)
 192         tile.newT(rtc, rtf, rbc, null, 4);
 193         tile.newT(rbf, rbc, rtf, null, 4);
 194         tile.newT(rtn, rtc, rbn, null, 4);
 195         tile.newT(rbc, rbn, rtc, null, 4);
 196
 197         // front
 198         tile.newT(ltn, mtn, mbn, null, 5);
 199         tile.newT(ltn, mbn, lbn, null, 5);
 200         tile.newT(mtn, rtn, rbn, null, 5);
 201         tile.newT(mtn, rbn, mbn, null, 5);
 202
 203         // back
 204         tile.newT(mtf, ltf, mbf, null, 6);
 205         tile.newT(mbf, ltf, lbf, null, 6);
 206         tile.newT(rtf, mtf, rbf, null, 6);
 207         tile.newT(rbf, mtf, mbf, null, 6);
 208
 209         for(Matrix m : translations) {
 210             for(Mesh.T t1 : tile) {
 211                 for(Mesh.T t2 : tile) {
 212                     if (t1==t2) continue;
 213
 214                     if ((t1.v1().p.times(m).minus(t2.v1().p).mag() < Mesh.EPSILON) &&
 215                         (t1.v2().p.times(m).minus(t2.v3().p).mag() < Mesh.EPSILON) &&
 216                         (t1.v3().p.times(m).minus(t2.v2().p).mag() < Mesh.EPSILON)) {
 217                         t1.e1().bindEdge(t2.e3());
 218                         t1.e2().bindEdge(t2.e2());
 219                         t1.e3().bindEdge(t2.e1());
 220                     }
 221                     if ((t1.v2().p.times(m).minus(t2.v1().p).mag() < Mesh.EPSILON) &&
 222                         (t1.v3().p.times(m).minus(t2.v3().p).mag() < Mesh.EPSILON) &&
 223                         (t1.v1().p.times(m).minus(t2.v2().p).mag() < Mesh.EPSILON)) {
 224                         t1.e2().bindEdge(t2.e3());
 225                         t1.e3().bindEdge(t2.e2());
 226                         t1.e1().bindEdge(t2.e1());
 227                     }
 228                     if ((t1.v3().p.times(m).minus(t2.v1().p).mag() < Mesh.EPSILON) &&
 229                         (t1.v1().p.times(m).minus(t2.v3().p).mag() < Mesh.EPSILON) &&
 230                         (t1.v2().p.times(m).minus(t2.v2().p).mag() < Mesh.EPSILON)) {
 231                         t1.e3().bindEdge(t2.e3());
 232                         t1.e1().bindEdge(t2.e2());
 233                         t1.e2().bindEdge(t2.e1());
 234                     }
 235
 236                 }
 237             }
 238         }
 239
 240         //xMesh.Vertex mid = lbf.getE(mbn).shatter();
 241
 242         // rescale to match volume
 243         float factor = (float)Math.pow(tile.volume() / goal.volume(), 1.0/3.0);
 244         goal.transform(Matrix.scale(factor));
 245
 246         // translate to match centroid
 247         goal.transform(Matrix.translate(tile.centroid().minus(goal.centroid())));
 248         goal.makeVerticesImmutable();
 249
 250         //tx.e2.shatter();
 251         //tx.e3.shatter();
 252
 253
 254         tile.rebindPoints();
 255
 256         //mid.move(new Vec((float)0,0,(float)-0.05));
 257         //ltn.move(new Vec((float)0,0,(float)-0.05));
 258
 259         //mtf.move(new Vec(0, (float)-0.05, (float)0.05));
 260
 261
 262         System.out.println("tile volume: " + tile.volume());
 263         System.out.println("goal volume: " + goal.volume());
 264
 265         tile.error_against = goal;
 266         goal.error_against = tile;
 267     }
 268
 269     public synchronized void breakit() {
 270         int oldverts = verts;
 271         System.out.println("doubling vertices.");
 272         PriorityQueue<Mesh.E> es = new PriorityQueue<Mesh.E>();
 273         for(Mesh.T t : tile) {
 274             es.add(t.e1());
 275             es.add(t.e2());
 276             es.add(t.e3());
 277         }
 278         for(int i=0; i<oldverts; i++) {
 279             Mesh.E e = es.poll();
 280             verts++;
 281             //System.out.println("shatter " + e);
 282             e.shatter();
 283         }
 284         tile.rebindPoints();
 285     }
 286
 287     public synchronized void rand(float temp, Mesh.Vertex p) {
 288
 289         //p.reComputeError();
 290         p.reComputeErrorAround();
 291         double tile_error = tile.error();
 292         double goal_error = goal.error();
 293
 294         Vec v;
 295         /*
 296           Matrix inv = p.errorQuadric();
 297           v = new Vec(inv.d, inv.h, inv.l).norm().times(1/(float)300);
 298           if (p.quadric_count == 0) {
 299           v = goal.nearest(p.p).p.minus(p.p).norm().times(1/(float)300);
 300           }
 301         */
 302         Vec v2 = new Vec((random.nextFloat() - (float)0.5) / 500,
 303                          (random.nextFloat() - (float)0.5)  / 500,
 304                          (random.nextFloat() - (float)0.5)  / 500);
 305         //v = v.plus(v2.norm().times(1/(float)300));
 306         v = v2.norm().times(1/(float)300);
 307
 308         Matrix m = Matrix.translate(v);
 309
 310         boolean good = p.move(m, false);
 311         if (!good) { misses++; return; }
 312
 313         //p.reComputeErrorAround();
 314
 315         double new_tile_error = tile.error();
 316         double new_goal_error = goal.error();
 317         double tile_delta = (new_tile_error - tile_error) / tile_error;
 318         double goal_delta = (new_goal_error - goal_error) / goal_error;
 319         double delta = tile_delta + goal_delta;
 320         double swapProbability = Math.exp((-1 * delta) / temp);
 321         boolean doSwap = good && (Math.random() < swapProbability);
 322         //boolean doSwap = good && (tile_delta <= 0 && goal_delta <= 0);
 323         //boolean doSwap = good && (tile_delta + goal_delta <= 0);
 324         if (doSwap) {
 325             tile_error = new_tile_error;
 326             goal_error = new_goal_error;
 327             //System.out.println("error: " + tile_error + " / " + goal_error);
 328             hits++;
 329         } else {
 330             p.move(Matrix.translate(v.times(-1)), true);
 331             misses++;
 332         }
 333     }
 334
 335     float hits = 0;
 336     float misses = 0;
 337     public void anneal() throws Exception {
 338         float hightemp = 10;
 339         float temp = hightemp;
 340         float last = 10;
 341         float lastbreak = 10;
 342         while(true) {
 343             synchronized(this) {
 344                 double ratio = (hits+misses==0) ? 1 : (hits / (hits+misses));
 345                 hits = 0;
 346                 misses = 0;
 347                 float gamma = 0;
 348                 double acceptance = ratio;
 349                 accepts = (int)(Math.ceil(ratio*100));
 350                 temps = (int)(Math.ceil(temp*1000));
 351                 vertss = tile.size();
 352                 if (breaks > 0) { while (breaks>0) {
 353                         breaks--;
 354                         breakit();
 355                         //gamma = 1;
 356                         float t = temp;
 357                         temp = lastbreak;
 358                         gamma = 1;
 359                         lastbreak = t;
 360                         //temp = last * 0.8f;
 361                         //last = temp;
 362                         //temp = hightemp;
 363                     } } else
 364                     if (acceptance > 0.96) gamma = 0.4f;
 365                     else if (acceptance > 0.9) gamma = 0.5f;
 366                     else if (acceptance > 0.8) gamma = 0.65f;
 367                     else if (acceptance > 0.6) gamma = 0.7f;
 368                     else {
 369                         if (acceptance > 0.3) {
 370                             gamma = 0.9f;
 371                         } else if (acceptance > 0.15) {
 372                             gamma = 0.95f;
 373                         } else if (acceptance > 0.10) {
 374                             gamma = 0.98f;
 375                         } else {
 376                             breakit();
 377                         float t = temp;
 378                         temp = lastbreak;
 379                         gamma = 1;
 380                         lastbreak = t;
 381                             //gamma = 1;
 382                         //gamma = 0.99f;
 383                             //gamma = 1;
 384                             //temp = last * 0.8f;
 385                             //last = temp;
 386                             //temp = hightemp;
 387                         }
 388                     }
 389                 temp = temp * gamma;
 390
 391
 392                 HashSet<Mesh.Vertex> hs = new HashSet<Mesh.Vertex>();
 393                 for(Mesh.Vertex p : tile.vertices()) hs.add(p);
 394                 Mesh.Vertex[] pts = (Mesh.Vertex[])hs.toArray(new Mesh.Vertex[0]);
 395
 396                 int count = 0;
 397                 long then = System.currentTimeMillis();
 398                 for(int i=0; i<40; i++) {
 399                     if (anneal) {
 400                         count++;
 401                         Mesh.Vertex v = pts[Math.abs(random.nextInt()) % pts.length];
 402                         rand(temp,v);
 403                         v.recomputeFundamentalQuadricIfStale();
 404                         v.recomputeFundamentalQuadricIfNeighborChanged();
 405                     }
 406                     Thread.yield();
 407                     repaint();
 408                 }
 409                 System.out.println("temp="+temp + " ratio="+(Math.ceil(ratio*100)) + " " +
 410                                    "points_per_second=" +
 411                                    (count*1000)/((double)(System.currentTimeMillis()-then)));
 412
 413                 for(Mesh.Vertex p : goal.vertices()) p.recomputeFundamentalQuadricIfNeighborChanged();
 414
 415                 synchronized(safeTriangles) {
 416                     safeTriangles.clear();
 417                     for(Mesh.T t : tile) if (t.shouldBeDrawn()) safeTriangles.add(t);
 418                 }
 419             }
 420         }
 421     }
 422
 423
 424     public static void main(String[] s) throws Exception {
 425         StlFile stlf = new StlFile();
 426         stlf.load("fish.stl");
 427         //stlf.load("monkey.stl");
 428         Frame f = new Frame();
 429         Main main = new Main(stlf, f);
 430         f.pack();
 431         f.show();
 432         f.setSize(900, 900);
 433         f.doLayout();
 434         main.anneal();
 435     }
 436
 437 }