+ public float rescore() {
+ int num = 0;
+ double dist = 0;
+ HashSet<P> done = new HashSet<P>();
+ for(T t : ts)
+ for(P p : new P[] { t.p1(), t.p2(), t.p3() }) {
+ if (done.contains(p)) continue;
+ done.add(p);
+ p.rescore();
+ }
+ for(T t : ts)
+ for(P p : new P[] { t.p1(), t.p2(), t.p3() })
+ p.kdremove();
+ kd = new KDTree(3);
+ for(T t : ts)
+ for(P p : new P[] { t.p1(), t.p2(), t.p3() })
+ p.kdinsert();
+ return (float)(dist/num);
+ }
+
+ public void transform(M m) {
+ ArrayList<P> set = new ArrayList<P>();
+ set.addAll(ps.keySet());
+ for(P p : set) p.transform(m);
+ }
+
+ public V diagonal() {
+ float min_x = Float.MAX_VALUE;
+ float min_y = Float.MAX_VALUE;
+ float min_z = Float.MAX_VALUE;
+ float max_x = Float.MIN_VALUE;
+ float max_y = Float.MIN_VALUE;
+ float max_z = Float.MIN_VALUE;
+ for(P p : ps.keySet()) {
+ if (p.x < min_x) min_x = p.x;
+ if (p.y < min_y) min_y = p.y;
+ if (p.z < min_z) min_z = p.z;
+ if (p.x > max_x) max_x = p.x;
+ if (p.y > max_y) max_y = p.y;
+ if (p.z > max_z) max_z = p.z;
+ }
+ return new V(max_x - min_x, max_y - min_y, max_z - min_z);
+ }
+
+ public P centroid() {
+ float min_x = Float.MAX_VALUE;
+ float min_y = Float.MAX_VALUE;
+ float min_z = Float.MAX_VALUE;
+ float max_x = Float.MIN_VALUE;
+ float max_y = Float.MIN_VALUE;
+ float max_z = Float.MIN_VALUE;
+ for(P p : ps.keySet()) {
+ if (p.x < min_x) min_x = p.x;
+ if (p.y < min_y) min_y = p.y;
+ if (p.z < min_z) min_z = p.z;
+ if (p.x > max_x) max_x = p.x;
+ if (p.y > max_y) max_y = p.y;
+ if (p.z > max_z) max_z = p.z;
+ }
+ return new P((float)(max_x + min_x)/2,
+ (float)(max_y + min_y)/2,
+ (float)(max_z + min_z)/2);
+ }
+
+ public P newP(double x, double y, double z) { return newP((float)x, (float)y, (float)z); }
+ public P newP(float x, float y, float z) { return new P(x, y, z); }
+
+ public T newT(P p12, P p23, P p31, V norm) {
+ V norm2 = p31.minus(p12).cross(p23.minus(p12));
+ float dot = norm.dot(norm2);
+ //if (Math.abs(dot) < EPSILON) throw new Error("dot products within epsilon of each other: "+norm+" "+norm2);
+ if (dot < 0) { P p = p12; p12=p23; p23 = p; }
+ return newT(p12, p23, p31);
+ }
+
+
+ public float volume() {
+ double total = 0;
+ for(T t : ts) {
+ double area = t.area();
+ V origin_to_centroid = new V(newP(0, 0, 0), t.centroid());
+ boolean facingAway = t.norm().dot(origin_to_centroid) > 0;
+ double height = Math.abs(t.norm().dot(origin_to_centroid));
+ total += ((facingAway ? 1 : -1) * area * height) / 3.0;
+ }
+ return (float)total;
+ }
+
+ public P nearest(P p) {
+ Object[] results;
+ try { results = kd.nearest(new double[]{p.x,p.y,p.z},1); } catch (Exception e) { throw new Error(e); }
+ return (P)results[0];
+ }
+
+ public T newT(P p1, P p2, P p3) {
+ p1 = p1.register();
+ p2 = p2.register();
+ p3 = p3.register();
+ E e12 = p1.makeE(p2);
+ E e23 = p2.makeE(p3);
+ E e31 = p3.makeE(p1);
+ while(e12.next != e23 || e23.next != e31 || e31.next != e12) {
+ e12.makeAdjacent(e23);
+ e23.makeAdjacent(e31);
+ e31.makeAdjacent(e12);
+ }
+ T ret = e12.makeT();
+ if (e12.t == null) throw new Error();
+ if (e23.t == null) throw new Error();
+ if (e31.t == null) throw new Error();
+ return ret;
+ }
+
+ private char allname = 'A';
+
+ public M aspect = new M();
+ public M invaspect = new M();
+
+ /** [UNIQUE] point in 3-space */
+ public final class P {
+ char name;
+
+ float x, y, z;
+
+ int watch_count;
+ float watch_x;
+ float watch_y;
+ float watch_z;
+ P watch;
+
+ private E e; // some edge *leaving* this point
+ private M binding = new M();
+ private P bound_to = this;
+
+ private float oldscore = 0;
+
+ private boolean inserted = false;
+
+ public P register() {
+ P p2 = ps.get(this);
+ if (p2==null) { p2 = this; ps.put(this,this); name = allname++; }
+ return p2;
+ }
+
+ public void kdremove() {
+ if (!inserted) return;
+ inserted = false;
+ P p = this;
+ try { kd.delete(new double[]{p.x,p.y,p.z}); } catch (Exception e) { }
+ }
+ public void kdinsert() {
+ if (inserted) return;
+ inserted = true;
+ P p = this;
+ try { kd.insert(new double[]{p.x,p.y,p.z},this); } catch (Exception e) { throw new Error(e); }
+ }
+
+ public float score() { return oldscore; }
+ public void unscore() {
+ if (watch == null) return;
+ watch.watch_x -= x;
+ watch.watch_y -= y;
+ watch.watch_z -= z;
+ watch.watch_count--;
+ if (watch.watch_count==0) {
+ watch.watch_x = 0;
+ watch.watch_y = 0;
+ watch.watch_z = 0;
+ }
+ watch = null;
+ }
+ public P times(M m) { return m.times(this); }
+ public P partner() {
+ if (watch==null) return this;
+ return watch.times(score_against.aspect).times(invaspect);
+ }
+ public P watchback() {
+ if (watch_count==0) return partner();
+ return newP(watch_x/watch_count, watch_y/watch_count, watch_z/watch_count);
+ }
+ public void rescore() {
+ if (score_against == null) return;
+
+ score -= oldscore;
+ oldscore = 0;
+
+ if (watch != null) unscore();
+ P po = this.times(aspect).times(score_against.invaspect);
+ if (watch == null) {
+ watch = score_against.nearest(po);
+
+ // don't attract to vertices that face the other way
+ if (watch.norm().times(score_against.aspect).times(invaspect).dot(norm()) < 0) {
+ watch = null;
+ } else {
+ watch.watch_x += po.x;
+ watch.watch_y += po.y;
+ watch.watch_z += po.z;
+ watch.watch_count++;
+ }
+ }
+
+ double s1, s2;
+ if (watch_count==0) s1 = 0;
+ else s1 = this.distance(watch_x/watch_count, watch_y/watch_count, watch_z/watch_count);
+ s2 = watch==null ? 0 : po.distance(watch);
+ oldscore = (float)(s1 + s2);
+ score += oldscore;
+ }
+
+ public float distance(P p) { return distance(p.x, p.y, p.z); }
+ public float distance(float ox, float oy, float oz) {
+ return
+ (float)Math.sqrt((x-ox)*(x-ox)+
+ (y-oy)*(y-oy)+
+ (z-oz)*(z-oz));
+ }
+
+ /** does NOT update bound pairs! */
+ public boolean transform(M m) {
+ // FIXME: screws up kdtree
+ // FIXME: screws up hashmap
+ unscore();
+ try {
+ if (ps.get(this)==null) throw new Error();
+ ps.remove(this);
+ float newx = m.a*x + m.b*y + m.c*z + m.d;
+ float newy = m.e*x + m.f*y + m.g*z + m.h;
+ float newz = m.i*x + m.j*y + m.k*z + m.l;
+ this.x = newx;
+ this.y = newy;
+ this.z = newz;
+ // FIXME: what if we move onto exactly where another point is?
+ ps.put(this,this);
+ } catch (Exception e) {
+ throw new RuntimeException(e);
+ }
+ rescore();
+ boolean good = true;
+ /*
+ for(T t : ts) {
+ for(E e = this.e; ;) {
+ if (e.intersects(t)) { good = false; break; }
+ e = e.pair.next;
+ if (e == this.e) break;
+ }
+ }
+ */
+ /*
+ if (t==this.t) continue;
+ if (this.intersects(t)) good = false;
+ }
+ */
+ return good;
+ }
+ public boolean move(V v) {
+ M m = new M(v);
+ P p = this;
+ boolean good = true;
+ do {
+ good &= p.transform(m);
+ v = v.times(binding); // bleh wrong
+ p = p.bound_to;
+ } while (p != this);
+ return good;
+ }
+
+ public E makeE(P p2) {
+ p2 = p2.register();
+ E e = getE(p2);
+ if (e != null) return e;
+ e = p2.getE(this);
+ if (this.e == null && p2.e == null) return this.e = new E(this, p2);
+ if (this.e == null && p2.e != null) return p2.makeE(this).pair;
+ return new E(getFreeIncident(), p2);
+ }
+
+ public E getFreeIncident() {
+ E ret = getFreeIncident(e, e);
+ if (ret != null) return ret;
+ ret = getFreeIncident(e.pair.next, e.pair.next);
+ if (ret == null) throw new Error("unable to find free incident to " + this);
+ return ret;
+ }
+
+ public E getFreeIncident(E start, E before) {
+ E e = start;
+ do {
+ if (e.pair.p2 == this && e.pair.t == null && e.pair.next.t == null) return e.pair;
+ e = e.pair.next;
+ } while(e != before);
+ return null;
+ }
+
+ public E getE(P p2) {
+ p2 = p2.register();
+ E e = this.e;
+ do {
+ if (e==null) return null;
+ if (e.p1 == this && e.p2 == p2) return e;
+ e = e.pair.next;
+ } while (e!=this.e);
+ return null;
+ }
+
+ public boolean isBoundTo(P p) {
+ p = p.register();
+ P px = p;
+ do {
+ if (px==this) return true;
+ px = px.bound_to;
+ } while(px != p);
+ return false;
+ }
+
+ public void unbind() { bound_to = this; binding = new M(); }
+ public void bind(P p) { bind(p, new M()); }
+ public void bind(P p, M binding) {
+ p = p.register();
+ if (isBoundTo(p)) return;
+ P temp_bound_to = p.bound_to;
+ M temp_binding = p.binding;
+ p.bound_to = this.bound_to;
+ p.binding = binding.times(this.binding); // FIXME: may have order wrong here
+ this.bound_to = temp_bound_to;
+ this.binding = temp_binding.times(temp_binding); // FIXME: may have order wrong here
+ }
+
+ public P(float x, float y, float z) {
+ this.x = x; this.y = y; this.z = z;
+ }
+
+ public V minus(P p) { return new V(x-p.x, y-p.y, z-p.z); }
+ public P plus(V v) { return newP(x+v.x, y+v.y, z+v.z); }