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         /*
 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();
 286             Thread.yield();
 287             repaint();
 288         }
 289         tile.rebindPoints();
 290         */
 291     }
 292
 293     public synchronized void rand(float temp, Mesh.Vertex p) {
 294
 295         p.reComputeErrorAround();
 296         double tile_error = tile.error();
 297         double goal_error = goal.error();
 298
 299         float max = p.averageEdgeLength()/10;
 300         Vec v = new Vec(random.nextFloat(), random.nextFloat(), random.nextFloat());
 301         v = v.norm().times((random.nextFloat() - 0.5f) * max);
 302         //System.out.println(max + " " + p.averageEdgeLength() + " " + v.mag());
 303         Matrix m = Matrix.translate(v);
 304
 305         boolean good = p.move(m, false);
 306         if (!good) { misses++; return; }
 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)/1000000));
 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.8f;
 357                 else if (acceptance > 0.15)   gamma = 0.9f;
 358                 else if (acceptance > 0.05)   gamma = 0.95f;
 359                 else if (acceptance > 0.01)   gamma = 0.98f;
 360                 else breaks++;
 361
 362                 if (seek_upward) {
 363                     if (acceptance > 0.2) seek_upward = false;
 364                     else gamma = 2-gamma;
 365                 }
 366
 367                 temp = temp * gamma;
 368
 369
 370                 HashSet<Mesh.Vertex> hs = new HashSet<Mesh.Vertex>();
 371                 for(Mesh.Vertex p : tile.vertices()) hs.add(p);
 372                 Mesh.Vertex[] pts = (Mesh.Vertex[])hs.toArray(new Mesh.Vertex[0]);
 373
 374                 int count = 0;
 375                 long then = System.currentTimeMillis();
 376                 for(int i=0; i<40; i++) {
 377                     if (anneal) {
 378                         count++;
 379                         Mesh.Vertex v = pts[Math.abs(random.nextInt()) % pts.length];
 380                         rand(temp,v);
 381                         v.recomputeFundamentalQuadricIfStale();
 382                         v.recomputeFundamentalQuadricIfNeighborChanged();
 383                     }
 384                     Thread.yield();
 385                     repaint();
 386                 }
 387                 PriorityQueue<Mesh.E> es = new PriorityQueue<Mesh.E>();
 388                 for(Mesh.T t : tile) {
 389                     float max = 5;
 390                     for(Mesh.E e : new Mesh.E[] { t.e1(), t.e2(), t.e3() }) {
 391                         if (e==null) continue;
 392                         if (e.stretchRatio() > max) es.add(e);
 393                         if (t.aspect() < 0.1) es.add(e);
 394                     }
 395                 }
 396
 397                 for(int i=0; i<1; i++) {
 398                     Mesh.E e = es.poll();
 399                     if (e==null) break;
 400                     e.shatter();
 401                 }
 402                 tile.rebindPoints();
 403
 404                 System.out.println("temp="+temp + " ratio="+(Math.ceil(acceptance*100)) + " " +
 405                                    "points_per_second=" +
 406                                    (count*1000)/((double)(System.currentTimeMillis()-then)));
 407
 408                 for(Mesh.Vertex p : goal.vertices()) p.recomputeFundamentalQuadricIfNeighborChanged();
 409
 410                 synchronized(safeTriangles) {
 411                     safeTriangles.clear();
 412                     for(Mesh.T t : tile) if (t.shouldBeDrawn()) safeTriangles.add(t);
 413                 }
 414             }
 415         }
 416     }
 417
 418
 419     public static void main(String[] s) throws Exception {
 420         StlFile stlf = new StlFile();
 421         stlf.load("fish.stl");
 422         //stlf.load("monkey.stl");
 423         Frame f = new Frame();
 424         Main main = new Main(stlf, f);
 425         f.pack();
 426         f.show();
 427         f.setSize(900, 900);
 428         f.doLayout();
 429         main.anneal();
 430     }
 431
 432 }