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 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             Thread.yield();
 278             repaint();
 279         }
 280         for(int i=0; i<Math.min(oldverts,200); i++) {
 281             Mesh.E e = es.poll();
 282             verts++;
 283             //System.out.println("shatter " + e);
 284             e.shatter();
 285             Thread.yield();
 286             repaint();
 287         }
 288         tile.rebindPoints();
 289     }
 290
 291     public synchronized void rand(float temp, Mesh.Vertex p) {
 292
 293         p.reComputeErrorAround();
 294         double tile_error = tile.error();
 295         double goal_error = goal.error();
 296
 297         float max = p.averageEdgeLength()/10;
 298         Vec v = new Vec(random.nextFloat(), random.nextFloat(), random.nextFloat());
 299         v = v.norm().times((random.nextFloat() - 0.5f) * max);
 300         //System.out.println(max + " " + p.averageEdgeLength() + " " + v.mag());
 301         Matrix m = Matrix.translate(v);
 302
 303         boolean good = p.move(m, false);
 304         if (!good) { misses++; return; }
 305
 306         p.reComputeErrorAround();
 307
 308         double new_tile_error = tile.error();
 309         double new_goal_error = goal.error();
 310         double tile_delta = (new_tile_error - tile_error) / tile_error;
 311         double goal_delta = (new_goal_error - goal_error) / goal_error;
 312         double delta = tile_delta + goal_delta;
 313         double swapProbability = Math.exp((-1 * delta) / (((double)temp)/10000));
 314         boolean doSwap = good && (Math.random() < swapProbability);
 315         //boolean doSwap = good && (tile_delta <= 0 && goal_delta <= 0);
 316         //boolean doSwap = good && (tile_delta + goal_delta <= 0);
 317         if (doSwap) {
 318             tile_error = new_tile_error;
 319             goal_error = new_goal_error;
 320             //System.out.println("error: " + tile_error + " / " + goal_error);
 321             hits++;
 322         } else {
 323             p.move(Matrix.translate(v.times(-1)), true);
 324             misses++;
 325         }
 326     }
 327
 328     float hits = 0;
 329     float misses = 0;
 330     public void anneal() throws Exception {
 331         float hightemp = 1;
 332         float temp = hightemp;
 333         float last = 10;
 334         boolean seek_upward = false;
 335         double acceptance = 1;
 336         while(true) {
 337             synchronized(this) {
 338                 double ratio = (hits+misses==0) ? 1 : (hits / (hits+misses));
 339                 hits = 0;
 340                 misses = 0;
 341                 float gamma = 1;
 342                 acceptance = (ratio+acceptance)/2;
 343                 accepts = (int)(Math.ceil(ratio*100));
 344                 temps = (int)(Math.ceil(temp*1000));
 345                 vertss = tile.size();
 346                 if (breaks > 0) {
 347                     while (breaks>0) {
 348                         breaks--;
 349                         breakit();
 350                     }
 351                     seek_upward = true;
 352                 } else if (acceptance > 0.96) gamma = 0.4f;
 353                 else if (acceptance > 0.9)    gamma = 0.5f;
 354                 else if (acceptance > 0.8)    gamma = 0.65f;
 355                 else if (acceptance > 0.6)    gamma = 0.7f;
 356                 else if (acceptance > 0.3)    gamma = 0.9f;
 357                 else if (acceptance > 0.15)   gamma = 0.95f;
 358                 else if (acceptance > 0.10)   gamma = 0.98f;
 359                 else breaks++;
 360
 361                 if (seek_upward) {
 362                     if (acceptance > 0.3) seek_upward = false;
 363                     else gamma = 2-gamma;
 364                 }
 365
 366                 temp = temp * gamma;
 367
 368
 369                 HashSet<Mesh.Vertex> hs = new HashSet<Mesh.Vertex>();
 370                 for(Mesh.Vertex p : tile.vertices()) hs.add(p);
 371                 Mesh.Vertex[] pts = (Mesh.Vertex[])hs.toArray(new Mesh.Vertex[0]);
 372
 373                 int count = 0;
 374                 long then = System.currentTimeMillis();
 375                 for(int i=0; i<40; i++) {
 376                     if (anneal) {
 377                         count++;
 378                         Mesh.Vertex v = pts[Math.abs(random.nextInt()) % pts.length];
 379                         rand(temp,v);
 380                         v.recomputeFundamentalQuadricIfStale();
 381                         v.recomputeFundamentalQuadricIfNeighborChanged();
 382                     }
 383                     Thread.yield();
 384                     repaint();
 385                 }
 386                 System.out.println("temp="+temp + " ratio="+(Math.ceil(acceptance*100)) + " " +
 387                                    "points_per_second=" +
 388                                    (count*1000)/((double)(System.currentTimeMillis()-then)));
 389
 390                 for(Mesh.Vertex p : goal.vertices()) p.recomputeFundamentalQuadricIfNeighborChanged();
 391
 392                 synchronized(safeTriangles) {
 393                     safeTriangles.clear();
 394                     for(Mesh.T t : tile) if (t.shouldBeDrawn()) safeTriangles.add(t);
 395                 }
 396             }
 397         }
 398     }
 399
 400
 401     public static void main(String[] s) throws Exception {
 402         StlFile stlf = new StlFile();
 403         stlf.load("fish.stl");
 404         //stlf.load("monkey.stl");
 405         Frame f = new Frame();
 406         Main main = new Main(stlf, f);
 407         f.pack();
 408         f.show();
 409         f.setSize(900, 900);
 410         f.doLayout();
 411         main.anneal();
 412     }
 413
 414 }