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         double new_tile_error = tile.error();
 307         double new_goal_error = goal.error();
 308         double tile_delta = (new_tile_error - tile_error) / tile_error;
 309         double goal_delta = (new_goal_error - goal_error) / goal_error;
 310         double delta = tile_delta + goal_delta;
 311         double swapProbability = Math.exp((-1 * delta) / (((double)temp)/1000000));
 312         boolean doSwap = good && (Math.random() < swapProbability);
 313         //boolean doSwap = good && (tile_delta <= 0 && goal_delta <= 0);
 314         //boolean doSwap = good && (tile_delta + goal_delta <= 0);
 315         if (doSwap) {
 316             tile_error = new_tile_error;
 317             goal_error = new_goal_error;
 318             //System.out.println("error: " + tile_error + " / " + goal_error);
 319             hits++;
 320             p.goodp = p.p;
 321         } else {
 322             p.move(Matrix.translate(v.times(-1)), true);
 323             misses++;
 324         }
 325     }
 326
 327     float hits = 0;
 328     float misses = 0;
 329     public void anneal() throws Exception {
 330         float hightemp = 1;
 331         float temp = hightemp;
 332         float last = 10;
 333         boolean seek_upward = false;
 334         double acceptance = 1;
 335         while(true) {
 336             synchronized(this) {
 337                 double ratio = (hits+misses==0) ? 1 : (hits / (hits+misses));
 338                 hits = 0;
 339                 misses = 0;
 340                 float gamma = 1;
 341                 acceptance = (ratio+acceptance)/2;
 342                 accepts = (int)(Math.ceil(ratio*100));
 343                 temps = (int)(Math.ceil(temp*1000));
 344                 vertss = tile.size();
 345                 if (breaks > 0) {
 346                     while (breaks>0) {
 347                         breaks--;
 348                         breakit();
 349                     }
 350                     seek_upward = true;
 351                 } else if (acceptance > 0.96) gamma = 0.4f;
 352                 else if (acceptance > 0.9)    gamma = 0.5f;
 353                 else if (acceptance > 0.8)    gamma = 0.65f;
 354                 else if (acceptance > 0.6)    gamma = 0.7f;
 355                 else if (acceptance > 0.3)    gamma = 0.8f;
 356                 else if (acceptance > 0.15)   gamma = 0.9f;
 357                 else if (acceptance > 0.05)   gamma = 0.95f;
 358                 else if (acceptance > 0.01)   gamma = 0.98f;
 359                 else { /*breaks++;*/ }
 360
 361                 if (seek_upward) {
 362                     if (acceptance > 0.2) 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                 PriorityQueue<Mesh.E> es = new PriorityQueue<Mesh.E>();
 387                 for(Mesh.T t : tile) {
 388                     float max = 5;
 389                     for(Mesh.E e : new Mesh.E[] { t.e1(), t.e2(), t.e3() }) {
 390                         if (e==null) continue;
 391                         if (e.stretchRatio() > max) es.add(e);
 392                         if (t.aspect() < 0.1 && e.length()>e.next.length() && e.length()>e.prev.length()) es.add(e);
 393                     }
 394                 }
 395
 396                 for(int i=0; i<5; i++) {
 397                     Mesh.E e = es.poll();
 398                     if (e==null) break;
 399                     e.shatter();
 400                 }
 401                 tile.rebindPoints();
 402
 403                 System.out.println("temp="+temp + " ratio="+(Math.ceil(acceptance*100)) + " " +
 404                                    "points_per_second=" +
 405                                    (count*1000)/((double)(System.currentTimeMillis()-then)));
 406
 407                 for(Mesh.Vertex p : goal.vertices()) p.recomputeFundamentalQuadricIfNeighborChanged();
 408
 409                 synchronized(safeTriangles) {
 410                     safeTriangles.clear();
 411                     for(Mesh.T t : tile)
 412                         if (t.shouldBeDrawn())
 413                             safeTriangles.add(t);
 414                 }
 415             }
 416         }
 417     }
 418
 419
 420     public static void main(String[] s) throws Exception {
 421         StlFile stlf = new StlFile();
 422         stlf.load("fish.stl");
 423         //stlf.load("monkey.stl");
 424         Frame f = new Frame();
 425         Main main = new Main(stlf, f);
 426         f.pack();
 427         f.show();
 428         f.setSize(900, 900);
 429         f.doLayout();
 430         main.anneal();
 431     }
 432
 433 }