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                         t2.e3().bindEdge(t1.e1(), m);
 218                         t2.e2().bindEdge(t1.e2(), m);
 219                         t2.e1().bindEdge(t1.e3(), m);
 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                         t2.e3().bindEdge(t1.e2(), m);
 225                         t2.e2().bindEdge(t1.e3(), m);
 226                         t2.e1().bindEdge(t1.e1(), m);
 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                         t2.e3().bindEdge(t1.e3(), m);
 232                         t2.e2().bindEdge(t1.e1(), m);
 233                         t2.e1().bindEdge(t1.e2(), m);
 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(e.midpoint(), null, null, true, true);
 285
 286             e.shatter();
 287             Thread.yield();
 288             repaint();
 289         }
 290         tile.rebindPoints();
 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             p.goodp = p.p;
 323         } else {
 324             p.move(Matrix.translate(v.times(-1)), true);
 325             misses++;
 326         }
 327     }
 328
 329     float hits = 0;
 330     float misses = 0;
 331     public void anneal() throws Exception {
 332         float hightemp = 1;
 333         float temp = hightemp;
 334         float last = 10;
 335         boolean seek_upward = false;
 336         double acceptance = 1;
 337         while(true) {
 338             synchronized(this) {
 339                 double ratio = (hits+misses==0) ? 1 : (hits / (hits+misses));
 340                 hits = 0;
 341                 misses = 0;
 342                 float gamma = 1;
 343                 acceptance = (ratio+acceptance)/2;
 344                 accepts = (int)(Math.ceil(ratio*100));
 345                 temps = (int)(Math.ceil(temp*1000));
 346                 vertss = tile.size();
 347                 if (breaks > 0) {
 348                     while (breaks>0) {
 349                         breaks--;
 350                         breakit();
 351                     }
 352                     seek_upward = true;
 353                 } else if (acceptance > 0.96) gamma = 0.4f;
 354                 else if (acceptance > 0.9)    gamma = 0.5f;
 355                 else if (acceptance > 0.8)    gamma = 0.65f;
 356                 else if (acceptance > 0.6)    gamma = 0.7f;
 357                 else if (acceptance > 0.3)    gamma = 0.8f;
 358                 else if (acceptance > 0.15)   gamma = 0.9f;
 359                 else if (acceptance > 0.05)   gamma = 0.95f;
 360                 else if (acceptance > 0.01)   gamma = 0.98f;
 361                 else { /*breaks++;*/ }
 362
 363                 if (seek_upward) {
 364                     if (acceptance > 0.2) seek_upward = false;
 365                     else gamma = 2-gamma;
 366                 }
 367
 368                 temp = temp * gamma;
 369
 370
 371                 HashSet<Mesh.Vertex> hs = new HashSet<Mesh.Vertex>();
 372                 for(Mesh.Vertex p : tile.vertices()) hs.add(p);
 373                 Mesh.Vertex[] pts = (Mesh.Vertex[])hs.toArray(new Mesh.Vertex[0]);
 374
 375                 int count = 0;
 376                 long then = System.currentTimeMillis();
 377                 for(int i=0; i<40; i++) {
 378                     if (anneal) {
 379                         count++;
 380                         Mesh.Vertex v = pts[Math.abs(random.nextInt()) % pts.length];
 381                         rand(temp,v);
 382                         v.recomputeFundamentalQuadricIfStale();
 383                         v.recomputeFundamentalQuadricIfNeighborChanged();
 384                     }
 385                     Thread.yield();
 386                     repaint();
 387                 }
 388                 PriorityQueue<Mesh.E> es = new PriorityQueue<Mesh.E>();
 389                 for(Mesh.T t : tile) {
 390                     float max = 5;
 391                     for(Mesh.E e : new Mesh.E[] { t.e1(), t.e2(), t.e3() }) {
 392                         if (e==null) continue;
 393                         if (e.stretchRatio() > max) es.add(e);
 394                         if (t.aspect() < 0.1 && e.length()>e.next.length() && e.length()>e.prev.length()) es.add(e);
 395                     }
 396                 }
 397                 /*
 398                 for(int i=0; i<5; i++) {
 399                     Mesh.E e = es.poll();
 400                     if (e==null) break;
 401                     e.shatter();
 402                 }
 403                 */
 404                 tile.rebindPoints();
 405
 406                 System.out.println("temp="+temp + " ratio="+(Math.ceil(acceptance*100)) + " " +
 407                                    "points_per_second=" +
 408                                    (count*1000)/((double)(System.currentTimeMillis()-then)));
 409
 410                 for(Mesh.Vertex p : goal.vertices()) p.recomputeFundamentalQuadricIfNeighborChanged();
 411
 412                 synchronized(safeTriangles) {
 413                     safeTriangles.clear();
 414                     for(Mesh.T t : tile)
 415                         if (t.shouldBeDrawn())
 416                             safeTriangles.add(t);
 417                 }
 418             }
 419         }
 420     }
 421
 422
 423     public static void main(String[] s) throws Exception {
 424         StlFile stlf = new StlFile();
 425         stlf.load("fish.stl");
 426         //stlf.load("monkey.stl");
 427         Frame f = new Frame();
 428         Main main = new Main(stlf, f);
 429         f.pack();
 430         f.show();
 431         f.setSize(900, 900);
 432         f.doLayout();
 433         main.anneal();
 434     }
 435
 436 }