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         goal.ignorecollision = true;
  67
  68         // rotate to align major axis -- this probably needs to be done by a human.
  69         goal.transform(Matrix.rotate(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             Matrix.translate(new Vec(lshift,  depth,    halfup)),
  93             Matrix.translate(new Vec(rshift,  depth,    halfup)),
  94             Matrix.translate(new Vec(lshift, -depth,    halfup)),
  95             Matrix.translate(new Vec(rshift, -depth,    halfup)),
  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                         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.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 synchronized 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         }
 280         for(int i=0; i<oldverts; i++) {
 281             Mesh.E e = es.poll();
 282             verts++;
 283             //System.out.println("shatter " + e);
 284             e.shatter();
 285         }
 286         tile.rebindPoints();
 287     }
 288
 289     public synchronized void rand(float temp, Mesh.Vertex p) {
 290
 291         //p.reComputeError();
 292         p.reComputeErrorAround();
 293         double tile_error = tile.error();
 294         double goal_error = goal.error();
 295
 296         Vec v;
 297         /*
 298           Matrix inv = p.errorQuadric();
 299           v = new Vec(inv.d, inv.h, inv.l).norm().times(1/(float)300);
 300           if (p.quadric_count == 0) {
 301           v = goal.nearest(p.p).p.minus(p.p).norm().times(1/(float)300);
 302           }
 303         */
 304         Vec v2 = new Vec((random.nextFloat() - (float)0.5) / 500,
 305                          (random.nextFloat() - (float)0.5)  / 500,
 306                          (random.nextFloat() - (float)0.5)  / 500);
 307         //v = v.plus(v2.norm().times(1/(float)300));
 308         v = v2.norm().times(1/(float)300);
 309
 310         boolean good = p.move(v);
 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(v.times(-1));
 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 }