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