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