convex polyhedron

sketch.js

@@ -1,3 +1,6 @@
+let camera;
+let red;
+
 let charges = [];
 
 let sphere_radius = 200;
@@ -7,8 +10,6 @@ let earth = false;
 let skeleton = false;
 let planes = false;
 
-let red;
-let yellow;
 
 function preload() {
   earth_image = loadImage("atlas1.jpg");
@@ -16,8 +17,8 @@ function preload() {
 
 function setup() {
   createCanvas(600, 600, WEBGL);
+  camera = createCamera();
   red = color(0xbf, 0x00, 0x00);
-  yellow = color(0xff, 0xff, 0x00);
 }
 
 function draw() {
@@ -25,76 +26,92 @@ function draw() {
   background(50);
 
   make_lights();
-  if (physics) {
+  if (physics) move_charges(charges);
-    move_charges(charges);
-  }
-  if (skeleton) {
-    draw_skeleton(sphere_radius);
-  }
 
-  if (planes) {
+  draw_charges(sphere_radius);
-    let p_charge = charges[0];
+  draw_sphere(sphere_radius, 25);
-    let a_charge = charges[1];
+  if (skeleton) draw_skeleton(sphere_radius);
-    let b_charge = charges[2];
+  if (planes) draw_planes(sphere_radius);
-    p_charge.color = yellow;
+}
-    a_charge.color = yellow;
+
-    b_charge.color = yellow;
+function face_dist_sq([v1, v2, v3]) {
-    let p = p_charge.position;
+  const center = p5.Vector.add(v1, v2).add(v3).mult(1 / 3);
-    let a = a_charge.position;
+  return createVector(camera.eyeX, camera.eyeY, camera.eyeZ).sub(center).magSq();
-    let b = b_charge.position;
+};
-    let n = p5.Vector.sub(p, a).cross(p5.Vector.sub(p, b));
+
-    n.normalize(); // unnecessary
+function draw_planes(radius) {
-    push();
+  let faces = [];
-    strokeWeight(5);
+  for (let i = 2; i < charges.length; i += 1) {
-    stroke(0x7f);
+    for (let j = 1; j < i; j += 1) {
-    line(
+      for (let k = 0; k < j; k += 1) {
-      0, 0, 0,
+        // Check if p1 p2 p3 form a face of the convex polyhedron
-      n.x * sphere_radius, n.y * sphere_radius, n.z * sphere_radius,
+        // enclosing all vertices ...
-    )
+        const p1 = charges[i].position;
-    fill(0xff);
+        const p2 = charges[j].position;
-    strokeWeight(3);
+        const p3 = charges[k].position;
-    stroke(0x00);
+        const normal = p5.Vector.sub(p2, p1).cross(p5.Vector.sub(p3, p1));
+        // ... by checking if the other vertices are on the same
+        // side of the plane generated by p1 p2 p3
+        let plane_separates_vertices = false;
+        for (let r = 1; r < charges.length; r += 1) {
+          for (let s = 0; s < r; s += 1) {
+            if (
+              r === i || r === j || r === k ||
+              s === i || s === j || s === k
+            ) continue;
+            const q1 = charges[r].position;
+            const q2 = charges[s].position;
+            // Let l(t) := q1 + (q2 - q1) * t.
+            // L := { l(t) : 0 <= t <= 1 } is the line segment
+            // between q1 and q2.  L intersects the plane
+            // generated by p1 p2 p3 iff l(t) intersects with
+            // the plane for some 0 <= t <= 1.
+            //
+            // If L intersects the plane, q1 and q2 are on
+            // opposite sides of the plane generated by p1 p2 p3,
+            // so p1 p2 p3 can't be a face.  If we want the
+            // polyhedron to be convex.  Which we do.
+            //
+            // A point k is on the plane generated by p1 p2 p3 iff
+            // dot(k - p1, normal) = 0.  Let n := normal.
+            //
+            // dot(l(t) - p1, n) = 0
+            // iff dot(q1 + (q2 - q1) * t - p1, n) = 0
+            // iff dot(q1 - p1, n) + dot(q2 - q2, n) * t = 0
+            // iff t = dot(p1 - q1, n) / dot(q2 - q2, n)
+            const t = (
+              p5.Vector.dot(p5.Vector.sub(p1, q1), normal) /
+              p5.Vector.dot(p5.Vector.sub(q2, q1), normal)
+            );
+            plane_separates_vertices ||= t >= 0 && t <= 1;
+            if (plane_separates_vertices) break;
+          }
+          if (plane_separates_vertices) break;
+        }
+        if (!plane_separates_vertices) {
+          faces.push([
+            p5.Vector.mult(p1, radius),
+            p5.Vector.mult(p2, radius),
+            p5.Vector.mult(p3, radius),
+          ]);
+        }
+      }
+    }
+  }
+  // fix OpenGL stacking alpha behaviour
+  faces.sort((a, b) => face_dist_sq(b) - face_dist_sq(a));
+  push();
+  strokeWeight(2);
+  stroke(0x00);
+  for ([v1, v2, v3] of faces) {
+    fill(0xbf, 0x7f);
     beginShape(TRIANGLES);
-    const v1 = p5.Vector.mult(p, sphere_radius);
-    const v2 = p5.Vector.mult(a, sphere_radius);
-    const v3 = p5.Vector.mult(b, sphere_radius);
     vertex(v1.x, v1.y, v1.z);
     vertex(v2.x, v2.y, v2.z);
     vertex(v3.x, v3.y, v3.z);
     endShape();
-    pop();
-    for (let i = 4; i < charges.length; i += 1) {
-      for (let j = 3; j < i; j += 1) {
-        push();
-        const u = charges[i].position;
-        const v = charges[j].position;
-        const t = p5.Vector.dot(p5.Vector.sub(p, u), n) / p5.Vector.dot(p5.Vector.sub(v, u), n);
-        const intersects_plane = t >= 0 && t <= 1;
-        if (intersects_plane) {
-          stroke(0xff, 0x1f, 0x00);
-        } else {
-          stroke(0x00, 0xff, 0x00);
-        }
-        strokeWeight(3);
-        line(
-          charges[i].position.x * sphere_radius,
-          charges[i].position.y * sphere_radius,
-          charges[i].position.z * sphere_radius,
-          charges[j].position.x * sphere_radius,
-          charges[j].position.y * sphere_radius,
-          charges[j].position.z * sphere_radius,
-        );
-        pop();
-      }
-    }
-  } else {
-    for (let charge of charges) {
-      charge.color = red;
-    }
   }
-
+  pop();
-  draw_charges(sphere_radius);
-  draw_sphere(sphere_radius, 25);
 }
 
 function draw_skeleton(radius) {