+ public T newT(P p12, P p23, P p31, Vec norm) {
+ Vec 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();
+ Vec origin_to_centroid = new Vec(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;
+ }
+
+ public M aspect = new M();
+ public M invaspect = new M();
+
+ public final class V extends P {
+ public V(P p) {
+ super(p.x, p.y, p.z);
+ if (ps.get(p) != null) throw new Error();
+ ps.put(this, this);
+ }
+ }
+
+ public class P {
+ 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 V register() {
+ V v = ps.get(this);
+ if (v==null) v = new V(this);
+ return (V)v;
+ }
+
+ 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((V)this,(V)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(Vec 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 Vec minus(P p) { return new Vec(x-p.x, y-p.y, z-p.z); }
+ public P plus(Vec v) { return newP(x+v.x, y+v.y, z+v.z); }