checkpoint
[anneal.git] / src / Geom.java
1 import java.awt.*;
2 import java.util.*;
3 import java.awt.event.*;
4 import javax.swing.*;
5 import javax.media.opengl.*;
6 import javax.media.opengl.glu.*;
7 import edu.wlu.cs.levy.CG.KDTree;
8
9 public class Geom implements Iterable<Geom.T> {
10
11     private KDTree kd = new KDTree(3);
12
13     public static float EPSILON = (float)0.000001;
14     public static Random random = new Random();
15
16     private HashMap<P,P> ps = new HashMap<P,P>();
17     public PriorityQueue<E>   es = new PriorityQueue<E>();
18     //private HashSet<T>   ts = new HashSet<T>();
19     public ArrayList<T>   ts = new ArrayList<T>();
20
21     public Iterator<T> iterator() { return ts.iterator(); }
22
23     public P newP(double x, double y, double z) { return newP((float)x, (float)y, (float)z); }
24     public P newP(float x, float y, float z) {
25         P p = new P(x, y, z);
26         P p2 = ps.get(p);
27         if (p2 != null) return p2;
28         ps.put(p,p);
29         p.name = allname++;
30         //try { kd.insert(new double[]{p.x,p.y,p.z},p); } catch (Exception e) { throw new Error(e); }
31         return p;
32     }
33     
34     public T newT(P p12, P p23, P p31, V norm) {
35         V norm2 = p31.minus(p12).cross(p23.minus(p12));
36         float dot = norm.dot(norm2);
37         if (Math.abs(dot) < EPSILON) throw new Error("dot products within epsilon of each other: "+norm+" "+norm2);
38         if (dot < 0) { P p = p12; p12=p23; p23 = p; }
39         return newT(p12, p23, p31);
40     }
41
42     public T newT(P p1, P p2, P p3) {
43         E e12 = p1.makeE(p2);
44         E e23 = p2.makeE(p3);
45         E e31 = p3.makeE(p1);
46         while(e12.next != e23 || e23.next != e31 || e31.next != e12) {
47             e12.makeAdjacent(e23);
48             e23.makeAdjacent(e31);
49             e31.makeAdjacent(e12);
50         }
51         T ret = e12.makeT();
52         if (e12.t == null) throw new Error();
53         if (e23.t == null) throw new Error();
54         if (e31.t == null) throw new Error();
55         return ret;
56     }
57
58     private char allname = 'A';
59
60     /** [UNIQUE] point in 3-space */
61     public final class P {
62         char name;
63
64         float x, y, z;
65
66         private E e;                // some edge *leaving* this point
67         private M binding = new M();
68         private P bound_to = this;
69
70         public E makeE(P p2) {
71             E e = getE(p2);
72             if (e != null) return e;
73             e = p2.getE(this);
74             if (this.e == null && p2.e == null) return this.e = new E(this, p2);
75             if (this.e == null && p2.e != null) return p2.makeE(this).pair;
76             return new E(getFreeIncident(), p2);
77         }
78
79         public E getFreeIncident() {
80             E ret = getFreeIncident(e, e);
81             if (ret != null) return ret;
82             ret = getFreeIncident(e.pair.next, e.pair.next);
83             if (ret == null) throw new Error("unable to find free incident to " + this);
84             return ret;
85         }
86
87         public E getFreeIncident(E start, E before) {
88             E e = start;
89             do {
90                 if (e.pair.p2 == this && e.pair.t == null && e.pair.next.t == null) return e.pair;
91                 e = e.pair.next;
92             } while(e != before);
93             return null;
94         }
95
96         public E getE(P p2) {
97             E e = this.e;
98             do {
99                 if (e==null) return null;
100                 if (e.p1 == this && e.p2 == p2) return e;
101                 e = e.pair.next;
102             } while (e!=this.e);
103             return null;
104         }
105
106         public boolean isBoundTo(P p) {
107             P px = p;
108             do {
109                 if (px==this) return true;
110                 px = px.bound_to;
111             } while(px != p);
112             return false;
113         }
114
115         public void unbind() { bound_to = null; binding = null; }
116         public void bind(P p) { bind(p, new M()); }
117         public void bind(P p, M binding) {
118             if (isBoundTo(p)) return;
119             P temp_bound_to = p.bound_to;
120             M temp_binding = p.binding;
121             p.bound_to = this.bound_to;
122             p.binding = binding.times(this.binding); // FIXME: may have order wrong here
123             this.bound_to = temp_bound_to;
124             this.binding = temp_binding.times(temp_binding); // FIXME: may have order wrong here
125         }
126
127         public void move(V v) {
128             P p = this;
129             do {
130                 p.x = p.x+v.x;
131                 p.y = p.y+v.y;
132                 p.z = p.z+v.z;
133                 v = v.times(binding);
134                 p = p.bound_to;
135             } while (p != this);
136         }
137
138         public P(float x, float y, float z) {
139             this.x = x; this.y = y; this.z = z;
140         }
141         /*
142         public P nearest() {
143             Object[] results;
144             try { results = kd.nearest(new double[]{x,y,z},2); } catch (Exception e) { throw new Error(e); }
145             if (results[0] != this) throw new Error();
146             return (P)results[1];
147         }
148         */
149         public V minus(P p) { return new V(x-p.x, y-p.y, z-p.z); }
150         public P plus(V v) { return newP(x+v.x, y+v.y, z+v.z); }
151         public P times(M m) { return m.apply(this); }
152         public boolean equals(Object o) {
153             if (o==null || !(o instanceof P)) return false;
154             P p = (P)o;
155             return p.x==x && p.y==y && p.z==z;
156         }
157         // FIXME: moving a point alters its hashCode
158         public int hashCode() {
159             return
160                 Float.floatToIntBits(x) ^
161                 Float.floatToIntBits(y) ^
162                 Float.floatToIntBits(z);
163         }
164         public void glVertex(GL gl) { gl.glVertex3f(x, y, z); }
165         public String toString() { return "("+x+","+y+","+z+")"; }
166         public V norm() {
167             V norm = new V(0, 0, 0);
168             E e = this.e;
169             do {
170                 if (e.t != null) norm = norm.plus(e.t.norm().times((float)e.prev.angle()));
171                 e = e.pair.next;
172             } while(e != this.e);
173             return norm.norm();
174         }
175     }
176
177     /** vector in 3-space */
178     public final class V {
179         public final float x, y, z;
180         public V(double x, double y, double z) { this((float)x, (float)y, (float)z); }
181         public V(float x, float y, float z) { this.x = x; this.y = y; this.z = z; }
182         public V cross(V v) { return new V(y*v.z-z*v.y, z*v.x-x*v.z, x*v.y-y*v.x); }
183         public V plus(V v) { return new V(x+v.x, y+v.y, z+v.z); }
184         public V norm() { return div(mag()); }
185         public V times(M m) { return m.apply(this); }
186         public float mag() { return (float)Math.sqrt(x*x+y*y+z*z); }
187         public float dot(V v) { return x*v.x + y*v.y + z*v.z; }
188         public V times(float mag) { return new V(x*mag, y*mag, z*mag); }
189         public V div(float mag) { return new V(x/mag, y/mag, z/mag); }
190         public String toString() { return "<"+x+","+y+","+z+">"; }
191     }
192
193     public class BindingGroup {
194         public HashSet<E> es = new HashSet<E>();
195         public BindingGroup() { }
196         public BindingGroup(E e) { es.add(e); }
197         public void add(E e) {
198             if (e.bg != null) { merge(e.bg); return; }
199             es.add(e);
200             e.bg = this;
201         }
202         public void merge(BindingGroup bg) {
203             for(E e : bg.es) {
204                 e.bg = null;
205                 add(e);
206             }
207         }
208     }
209
210     /** [UNIQUE] an edge */
211     public final class E implements Comparable<E> {
212
213         public int compareTo(E e) {
214             return e.length() > length() ? 1 : -1;
215         }
216
217         public final P p1, p2;
218         T t;     // triangle to our "left"
219         E prev;  // previous half-edge
220         E next;  // next half-edge
221         E pair;  // partner half-edge
222
223
224         public BindingGroup bg = new BindingGroup(this);
225
226         public void bind(E e) { bind(e, new M()); }
227         public void bind(E e, M m) { e.bg.add(this); }
228
229         public void dobind() {
230             if (bg==null) return;
231             for(E ex : bg.es) {
232                 if (ex==this) continue;
233                 p1.bind(ex.p1);
234                 p2.bind(ex.p2);
235             }
236         }
237
238         boolean shattered = false;
239         public P shatter() { return shatter(midpoint(), null, null); }
240         public P shatter(P mid, BindingGroup bg1, BindingGroup bg2) {
241             if (shattered) return mid;
242             shattered = true;
243
244             P r = next.p2;
245             E next = this.next;
246             E prev = this.prev;
247
248             if (bg1==null) bg1 = new BindingGroup();
249             if (bg2==null) bg2 = new BindingGroup();
250             for(E e : bg.es) e.shatter(e.midpoint(), bg1, bg2);
251             pair.shatter();
252             destroy();
253
254             newT(r, p1, mid);
255             newT(r, mid, p2);
256             bg1.add(p1.getE(mid));
257             bg2.add(mid.getE(p2));
258             return mid;
259         }
260
261         public boolean destroyed = false;
262         public void destroy() {
263             if (destroyed) return;
264             destroyed = true;
265             pair.destroyed = true;
266             if (next.t != null) ts.remove(next.t);
267             if (prev.t != null) ts.remove(prev.t);
268             if (pair.next.t != null) ts.remove(pair.next.t);
269             if (pair.prev.t != null) ts.remove(pair.prev.t);
270             next.t = null;
271             prev.t = null;
272             pair.next.t = null;
273             pair.prev.t = null;
274             this.bg = null;
275             pair.bg = null;
276             pair.prev.next = next;
277             next.prev = pair.prev;
278             prev.next = pair.next;
279             pair.next = prev;
280             if (p1.e == this) p1.e = prev.next;
281             if (pair.p1.e == pair) pair.p1.e = pair.prev.next;
282             es.remove(this);
283             es.remove(pair);
284         }
285
286         private void sync() {
287             this.prev.next = this;
288             this.next.prev = this;
289             this.pair.pair = this;
290             if (this.next.p1 != p2) throw new Error();
291             if (this.prev.p2 != p1) throw new Error();
292             if (this.p1.e == null) this.p1.e = this;
293             es.add(this);
294         }
295
296         public T makeT() { return t==null ? (t = new T(this)) : t; }
297
298         /** angle between this half-edge and the next */
299         public double angle() {
300             V v1 = next.p2.minus(p2);
301             V v2 = this.p1.minus(p2);
302             return Math.acos(v1.norm().dot(v2.norm()));
303         }
304
305         public void makeAdjacent(E e) {
306             if (this.next == e) return;
307             if (p2 != e.p1) throw new Error("cannot make adjacent -- no shared vertex");
308             if (t != null || e.t != null) throw new Error("cannot make adjacent -- edges not both free");
309
310             E freeIncident = p2.getFreeIncident(e, this);
311
312             e.prev.next = freeIncident.next;
313             freeIncident.next.prev = e.prev;
314
315             freeIncident.next = this.next;
316             this.next.prev = freeIncident;
317             
318             this.next = e;
319             e.prev = this;
320
321             sync();
322             freeIncident.sync();
323         }
324
325         /** creates an isolated edge out in the middle of space */
326         public E(P p1, P p2) {
327             if (p1==p2) throw new Error("attempt to create edge with single vertex: " + p1);
328             this.p1 = p1;
329             this.p2 = p2;
330             this.prev = this.next = this.pair = new E(this, this, this);
331             sync();
332         }
333
334         /** adds a new half-edge from prev.p2 to p2 */
335         public E(E prev, P p2) {
336             this.p1 = prev.p2;
337             this.p2 = p2;
338             this.prev = prev;
339             if (p2.getE(p1) != null) throw new Error();
340             if (p2.e==null) {
341                 this.next = this.pair = new E(this, this, prev.next);
342             } else {
343                 E q = p2.getFreeIncident();
344                 this.next = q.next;
345                 this.next.prev = this;
346                 E z = prev.next;
347                 this.prev.next = this;
348                 this.pair = new E(q, this, z);
349             }
350             sync();
351         }
352
353         /** adds a new half-edge to the mesh with a given predecessor, successor, and pair */
354         public E(E prev, E pair, E next) {
355             this.p1 = prev.p2;
356             this.p2 = next.p1;
357             this.prev = prev;
358             this.next = next;
359             this.pair = pair;
360             sync();
361         }
362         public P midpoint() { return newP((p1.x+p2.x)/2, (p1.y+p2.y)/2, (p1.z+p2.z)/2); }
363         public boolean has(P p) { return p==p1 || p==p2; }
364         public float length() { return p1.minus(p2).mag(); }
365         public String toString() { return p1+"->"+p2; }
366     }
367
368     /** [UNIQUE] a triangle (face) */
369     public final class T {
370         public final E e1;
371         public final int color;
372
373         T(E e1) {
374             this.e1 = e1;
375             E e2 = e1.next;
376             E e3 = e2.next;
377             if (e1==e2     || e1==e3) throw new Error();
378             if (e3.next!=e1) throw new Error();
379             if (e1.t!=null || e2.t!=null || e3.t!=null)
380                 throw new Error("non-manifold surface or disagreeing normals");
381             e1.t = this;
382             e1.next.t = this;
383             e1.next.next.t = this;
384
385             // FIXME: check for sealed/watertight surface once construction is complete (and infer normal(s)?)
386
387             int color = Math.abs(random.nextInt());
388             while(true) {
389                 color = color % 4;
390                 if (e1().pair.t != null && color == e1().pair.t.color) { color++; continue; }
391                 if (e2().pair.t != null && color == e2().pair.t.color) { color++; continue; }
392                 if (e3().pair.t != null && color == e3().pair.t.color) { color++; continue; }
393                 break;
394             }
395
396             // FIXME unnecssary
397             ts.add(this);
398
399             this.color = color;
400         }
401         public P p1() { return e1.p1; }
402         public P p2() { return e1.p2; }
403         public P p3() { return e1.next.p2; }
404         public E e1() { return e1; }
405         public E e2() { return e1.next; }
406         public E e3() { return e1.prev; }
407         public V norm() { return p2().minus(p1()).cross(p3().minus(p1())).norm(); }
408         public boolean hasE(E e) { return e1==e || e1.next==e || e1.prev==e; }
409         public boolean has(P p) { return p1()==p || p2()==p || p3()==p; }
410         public void glVertices(GL gl) {
411             p1().glVertex(gl);
412             p2().glVertex(gl);
413             p3().glVertex(gl);
414         }
415
416         public P centroid() { return newP((p1().x+p2().x+p3().x)/3,
417                                           (p1().y+p2().y+p3().y)/3, 
418                                           (p1().z+p2().z+p3().z)/3); }
419         public float diameter() {
420             // FIXME: what is this supposed to be?
421             return Math.max(Math.max(e1().length(), e2().length()), e3().length()) / 2;
422         }
423
424
425     }
426
427
428     /** matrix */
429     public class M {
430         public M() { }
431         public P apply(P p) { return p; }
432         public V apply(V v) { return v; }
433         public M invert() { return this; }
434         public M times(M m) { return this; }
435     }
436
437     public void unbind() {
438         /*
439         for(Geom.T t : this) {
440             t.p1().unbind();
441             t.p2().unbind();
442             t.p3().unbind();
443         }
444         */
445     }
446     public void bind() {
447         for(Geom.T t : this) {
448             t.e1().dobind();
449             t.e2().dobind();
450             t.e3().dobind();
451         }
452     }
453 }