this rename is absolutely stupid

docs/chapters/pointvectors3d/content.en-GB.md

@@ -0,0 +1,116 @@
+# Working with 3D normals
+
+Before we move on to the next section we need to spend a little bit of time on the difference between 2D and 3D. While for many things this difference is irrelevant and the procedures are identical (for instance, getting the 3D tangent is just doing what we do for 2D, but for x, y, and z, instead of just for x and y), when it comes to normals things are a little more complex, and thus more work. Mind you, it's not "super hard", but there are more steps involved and we should have a look at those.
+
+Getting normals in 3D is in principle the same as in 2D: we take the normalised tangent vector, and then rotate it by a quarter turn. However, this is where things get that little more complex: we can turn in quite a few directions, since "the normal" in 3D is a plane, not a single vector, so we basically need to define what "the" normal is in the 3D case.
+
+The "naïve" approach is to construct what is known as the [Frenet normal](https://en.wikipedia.org/wiki/Frenet%E2%80%93Serret_formulas), where we follow a simple recipe that works in many cases (but does super bizarre things in some others). The idea is that even though there are infinitely many vectors that are perpendicular to the tangent (i.e. make a 90 degree angle with it), the tangent itself sort of lies on its own plane already: since each point on the curve (no matter how closely spaced) has its own tangent vector, we can say that each point lies in the same plane as the local tangent, as well as the tangents "right next to it".
+
+Even if that difference in tangent vectors is minute, "any difference" is all we need to find out what that plane is - or rather, what the vector perpendicular to that plane is. Which is what we need: if we can calculate that vector, and we have the tangent vector that we know lies on a plane, then we can rotate the tangent vector over the perpendicular, and presto. We have computed the normal using the same logic we used for the 2D case: "just rotate it 90 degrees".
+
+So let's do that! And in a twist surprise, we can do this in four lines:
+
+- **a** = normalize(B'(t))
+- **b** = normalize(**a** + B''(t))
+- **r** = normalize(**b** × **a**)
+- **normal** = normalize(**r** × **a**)
+
+Let's unpack that a little:
+
+- We start by taking the [normalized vector](https://en.wikipedia.org/wiki/Unit_vector) for the derivative at some point on the curve. We normalize it so the maths is less work. Less work is good.
+- Then, we compute **b** which represents what a next point's tangent would be if the curve stopped changing at our point and just had the same derivative and second derivative from that point on.
+- This lets us find two vectors (the derivative, and the second derivative added to the derivative) that lie on the same plane, which means we can use them to compute a vector perpendicular to that plane, using an elementary vector operation called the [cross product](https://en.wikipedia.org/wiki/Cross_product). (Note that while that operation uses the × operator, it's most definitely not a multiplication!) The result of that gives us a vector that we can use as the "axis of rotation" for turning the tangent a quarter circle to get our normal, just like we did in the 2D case.
+- Since the cross product lets us find a vector that is perpendicular to some plane defined by two other vectors, and since the normal vector should be perpendicular to the plane that the tangent and the axis of rotation lie in, we can use the cross product a second time, and immediately get our normal vector.
+
+And then we're done, we found "the" normal vector for a 3D curve. Let's see what that looks like for a sample curve, shall we? You can move your cursor across the graphic from left to right, to show the normal at a point with a t value that is based on your cursor position: all the way on the left is 0, all the way on the right = 1, midway is t=0.5, etc:
+
+<graphics-element title="Some known and unknown vectors" width="350" height="300" src="./frenet.js"></graphics-element>
+
+However, if you've played with that graphic a bit, you might have noticed something odd. The normal seems to "suddenly twist around" around between t=0.5 and t=0.9... Why is doing that?
+
+As it turns out, it's doing that because that's how the maths works, and that's the problem with Frenet normals: while they are "mathematically correct", they are "practically problematic", and so for any kind of graphics work what we really want is a way to compute normals that just... look good.
+
+Thankfully, Frenet normals are not our only option.
+
+Another option is to take a slightly more algorithmic approach and compute a form of [Rotation Minimising Frame](https://www.microsoft.com/en-us/research/wp-content/uploads/2016/12/Computation-of-rotation-minimizing-frames.pdf) (also known as "parallel transport frame" or "Bishop frame") instead, where a "frame" is a set made up of the tangent, the rotational axis, and the normal vector, centered on an on-curve point.
+
+These type of frames are computed based on "the previous frame", so we cannot simply compute these "on demand" for single points, as we could for Frenet frames; we have to compute them for the entire curve. Thankfully, the procedure is pretty simple, and can be performed at the same time that you're building lookup tables for your curve.
+
+The idea is to take a starting "tangent/rotation axis/normal" frame at t=0, and then compute what the next frame "should" look like by applying some rules that yield a good looking next frame. In the case of the RMF paper linked above, those rules are:
+
+- Take a point on the curve for which we know the RM frame already,
+- take a next point on the curve for which we don't know the RM frame yet, and
+- reflect the known frame onto the next point, by treating the plane through the curve at the point exactly between the next and previous points as a "mirror".
+- This gives the next point a tangent vector that's essentially pointing in the opposite direction of what it should be, and a normal that's slightly off-kilter, so:
+- reflect the vectors of our "mirrored frame" a second time, but this time using the plane through the "next point" itself as "mirror".
+- Done: the tangent and normal have been fixed, and we have a good looking frame to work with.
+
+So, let's write some code for that!
+
+<div class="howtocode">
+
+### Implementing Rotation Minimising Frames
+
+We first assume we have a function for calculating the Frenet frame at a point, which we already discussed above, inn a way that it yields a frame with properties:
+
+```
+{
+  o: origin of all vectors, i.e. the on-curve point,
+  t: tangent vector,
+  r: rotational axis vector,
+  n: normal vector
+}
+```
+
+Then, we can write a function that generates a sequence of RM frames in the following manner:
+
+```
+generateRMFrames(steps) -> frames:
+  step = 1.0/steps
+
+  // Start off with the standard tangent/axis/normal frame
+  // associated with the curve at t=0:
+  frames.add(getFrenetFrame(0))
+
+  // start constructing RM frames:
+  for t0 = 0, t0 < 1.0, t0 += step:
+    // start with the previous, known frame
+    x0 = frames.last
+
+    // get the next frame: we're going to keep its position and tangent,
+    // but we're going to recompute the axis and normal.
+    t1 = t0 + step
+    x1 = { o: getPoint(t1), t: getDerivative(t) }
+
+    // First we reflect x0's tangent and axis of rotation onto x1,
+    // through/ the plane of reflection at the point between x0 x1
+    v1 = x1.o - x0.o
+    c1 = v1 · v1
+    riL = x0.r - v1 * 2/c1 * v1 · x0.r
+    tiL = x0.t - v1 * 2/c1 * v1 · x0.t
+
+    // note that v1 is a vector, but 2/c1 and (v1 · ...) are just
+    // plain numbers, so we're just scaling v1 by some constant.
+
+    // Then we reflect a second time, over a plane at x1, so that
+    // the frame tangent is aligned with the curve tangent again:
+    v2 = x1.t - tiL
+    c2 = v2 · v2
+
+    // and we're done here:
+    x1.r = riL - v2 * 2/c2 * v2 · riL
+    x1.n = x1.r × x1.t
+    frames.add(x1)
+```
+
+Ignoring comments, this is certainly more code than when we were just computing a single Frenet frame, but it's not a crazy amount more code to get much better looking normals.
+
+</div>
+
+Speaking of better looking, what does this actually look like? Let's revisit that earlier curve, but this time use rotation minimising frames rather than Frenet frames:
+
+<graphics-element title="Some known and unknown vectors" width="350" height="300"  src="./rotation-minimizing.js"></graphics-element>
+
+That looks much better!
+
+For those reading along with the code: we don't even strictly speaking need a Frenet frame to start with: we could, for instance, treat the z-axis as our initial axis of rotation, so that our initial normal is **(0,0,1) × tangent**, and then take things from there, but having that initial "mathematically correct" frame so that the initial normal seems to line up based on the curve's orientation in 3D space is just nice.

docs/chapters/pointvectors3d/frenet.js

@@ -0,0 +1,157 @@
+import vec from "./vector-lib.js";
+import { project, projectXY, projectXZ, projectYZ } from "./projection.js";
+
+let d, cube;
+
+setup() {
+  d = this.width/2 + 25;
+  cube = [
+    {x:0, y:0, z:0},
+    {x:d, y:0, z:0},
+    {x:d, y:d, z:0},
+    {x:0, y:d, z:0},
+    {x:0, y:0, z:d},
+    {x:d, y:0, z:d},
+    {x:d, y:d, z:d},
+    {x:0, y:d, z:d}
+  ].map(p => project(p));
+  const points = this.points = [
+    {x:120, y:  0, z:  0},
+    {x:120, y:220, z:  0},
+    {x: 30, y:  0, z: 30},
+    {x:  0, y:  0, z:200}
+  ];
+  this.curve = new Bezier(this, points.map(p => project(p)));
+  this.cxy = new Bezier(this, points.map(p => projectXY(p)));
+  this.cxz = new Bezier(this, points.map(p => projectXZ(p)));
+  this.cyz = new Bezier(this, points.map(p => projectYZ(p)));
+  this.t = 0;
+}
+
+draw() {
+  resetTransform();
+  clear();
+  translate(this.width/2 - 60, this.height/2 + 75);
+  const curve = this.curve;
+
+  this.drawCurveProjections();
+
+  this.drawCubeBack();
+
+  curve.drawCurve(`grey`);
+  setStroke(`grey`)
+  line(curve.points[0].x, curve.points[0].y, curve.points[1].x, curve.points[1].y);
+  line(curve.points[2].x, curve.points[2].y, curve.points[3].x, curve.points[3].y);
+  curve.points.forEach(p => circle(p.x, p.y, 2));
+
+  this.drawPoint(this.t);
+
+  this.drawCubeFront();
+}
+
+drawCubeBack() {
+  const c = cube;
+
+  // x axis
+  setStroke(`blue`);
+  line(c[0].x, c[0].y, c[1].x, c[1].y);
+
+  // y axis
+  setStroke(`red`);
+  line(c[3].x, c[3].y, c[0].x, c[0].y);
+
+  // z axis
+  setStroke(`green`);
+  line(c[0].x, c[0].y, c[4].x, c[4].y);
+}
+
+drawCurveProjections() {
+  this.cxy.drawCurve(`#EEF`);
+  this.cxz.drawCurve(`#EEF`);
+  this.cyz.drawCurve(`#EEF`);
+}
+
+
+drawPoint(t) {
+  const {o, r, n, dt} = this.getFrenetVectors(t, this.points);
+
+  setStroke(`red`);
+  setFill(`red`);
+  const p = project(o);
+  circle(p.x, p.y, 3);
+
+  this.drawVector(p, vec.normalize(r),  40, `blue`, `r`);
+  this.drawVector(p, vec.normalize(n),  40, `red`, `n`);
+  this.drawVector(p, vec.normalize(dt), 40, `green`, `t′`);
+
+  setFill(`black`)
+  text(`t = ${t.toFixed(2)}`, p.x+10, p.y+15);
+}
+
+drawCubeFront() {
+  const c = cube;
+
+  // rest of the cube
+  setStroke("lightgrey");
+  line(c[1].x, c[1].y, c[2].x, c[2].y);
+  line(c[2].x, c[2].y, c[3].x, c[3].y);
+  line(c[1].x, c[1].y, c[5].x, c[5].y);
+  line(c[2].x, c[2].y, c[6].x, c[6].y);
+  line(c[3].x, c[3].y, c[7].x, c[7].y);
+  line(c[4].x, c[4].y, c[5].x, c[5].y);
+  line(c[5].x, c[5].y, c[6].x, c[6].y);
+  line(c[6].x, c[6].y, c[7].x, c[7].y);
+  line(c[7].x, c[7].y, c[4].x, c[4].y);
+}
+
+getFrenetVectors(t, originalPoints) {
+  const curve = new Bezier(this, originalPoints);
+  const d1curve = new Bezier(this, curve.dpoints[0]);
+  const dt = d1curve.get(t);
+  const ddt = d1curve.derivative(t);
+  const o = curve.get(t);
+  const b = vec.normalize(vec.plus(dt, ddt));
+  const r = vec.normalize(vec.cross(b, dt));
+  const n = vec.normalize(vec.cross(r, dt));
+  return { o, dt, r, n };
+}
+
+drawVector(from, vec, length, color, label) {
+  setStroke(color);
+  setFill(`black`);
+
+  let pt = project({
+    x: length * vec.x,
+    y: length * vec.y,
+    z: length * vec.z
+  });
+  line(from.x, from.y, from.x + pt.x, from.y + pt.y);
+
+  let txt = project({
+    x: (length+15) * vec.x,
+    y: (length+15) * vec.y,
+    z: (length+15) * vec.z
+  });
+  text(label, from.x + txt.x, from.y + txt.y);
+}
+
+onMouseMove() {
+  this.t = constrain(
+    map(this.cursor.x,0,this.width,-0.1, 1.1)
+    ,0,1
+  );
+  redraw();
+}
+
+onKeyDown() {
+  let key = this.keyboard.currentKey;
+  if (key === `ArrowUp`) {
+    this.t += 0.01;
+    if (this.t > 1) this.t = 1;
+  }
+  if (key === `ArrowDown`) {
+    this.t -= 0.01;
+    if (this.t < 0) this.t = 0;
+  }
+  redraw();
+}

docs/chapters/pointvectors3d/projection.js

@@ -0,0 +1,35 @@
+/**
+ * A cabinet projection utility
+ */
+
+// Universal projector function
+function project(point3d, offset = { x: 0, y: 0 }, phi = -Math.PI / 6) {
+  // what they rarely tell you: if you want Z to "go up",
+  // X to "come out of the screen", and Y to be the "left/right",
+  // we need to switch some coordinates around:
+  const x = point3d.y,
+    y = -point3d.z,
+    z = -point3d.x;
+
+  return {
+    x: offset.x + x + (z / 2) * Math.cos(phi),
+    y: offset.y + y + (z / 2) * Math.sin(phi),
+  };
+}
+
+// and some rebuilt planar projectors
+function projectXY(p, offset, phi) {
+  return project({ x: p.x, y: p.y, z: 0 }, offset, phi);
+}
+
+function projectXZ(p, offset, phi) {
+  return project({ x: p.x, y: 0, z: p.z }, offset, phi);
+}
+
+function projectYZ(p, offset, phi) {
+  return project({ x: 0, y: p.y, z: p.z }, offset, phi);
+}
+
+
+export { project, projectXY, projectXZ, projectYZ }
+

docs/chapters/pointvectors3d/rotation-minimizing.js

@@ -0,0 +1,213 @@
+import vec from "./vector-lib.js";
+import { project, projectXY, projectXZ, projectYZ } from "./projection.js";
+
+let d, cube;
+
+setup() {
+  d = this.width/2 + 25;
+  cube = [
+    {x:0, y:0, z:0},
+    {x:d, y:0, z:0},
+    {x:d, y:d, z:0},
+    {x:0, y:d, z:0},
+    {x:0, y:0, z:d},
+    {x:d, y:0, z:d},
+    {x:d, y:d, z:d},
+    {x:0, y:d, z:d}
+  ].map(p => project(p));
+  const points = this.points = [
+    {x:120, y:  0, z:  0},
+    {x:120, y:220, z:  0},
+    {x: 30, y:  0, z: 30},
+    {x:  0, y:  0, z:200}
+  ];
+  this.curve = new Bezier(this, points.map(p => project(p)));
+  this.cxy = new Bezier(this, points.map(p => projectXY(p)));
+  this.cxz = new Bezier(this, points.map(p => projectXZ(p)));
+  this.cyz = new Bezier(this, points.map(p => projectYZ(p)));
+  this.t = 0;
+}
+
+draw() {
+  resetTransform();
+  clear();
+  translate(this.width/2 - 60, this.height/2 + 75);
+  const curve = this.curve;
+
+  this.drawCurveProjections();
+
+  this.drawCubeBack();
+
+  curve.drawCurve(`grey`);
+  setStroke(`grey`)
+  line(curve.points[0].x, curve.points[0].y, curve.points[1].x, curve.points[1].y);
+  line(curve.points[2].x, curve.points[2].y, curve.points[3].x, curve.points[3].y);
+  curve.points.forEach(p => circle(p.x, p.y, 2));
+
+  this.drawPoint(this.t);
+
+  this.drawCubeFront();
+}
+
+drawCubeBack() {
+  const c = cube;
+
+  // x axis
+  setStroke(`blue`);
+  line(c[0].x, c[0].y, c[1].x, c[1].y);
+
+  // y axis
+  setStroke(`red`);
+  line(c[3].x, c[3].y, c[0].x, c[0].y);
+
+  // z axis
+  setStroke(`green`);
+  line(c[0].x, c[0].y, c[4].x, c[4].y);
+}
+
+drawCurveProjections() {
+  this.cxy.drawCurve(`#EEF`);
+  this.cxz.drawCurve(`#EEF`);
+  this.cyz.drawCurve(`#EEF`);
+}
+
+drawPoint(t) {
+  const {o, r, n, dt } = this.getRMF(t, this.points);
+
+  setStroke(`red`);
+  setFill(`red`);
+  const p = project(o);
+  circle(p.x, p.y, 3);
+
+  this.drawVector(p, vec.normalize(r),  40, `blue`, `r`);
+  this.drawVector(p, vec.normalize(n),  40, `red`, `n`);
+  this.drawVector(p, vec.normalize(dt), 40, `green`, `t′`);
+
+  setFill(`black`)
+  text(`t = ${t.toFixed(2)}`, p.x+10, p.y+15);
+}
+
+drawCubeFront() {
+  const c = cube;
+
+  // rest of the cube
+  setStroke("lightgrey");
+  line(c[1].x, c[1].y, c[2].x, c[2].y);
+  line(c[2].x, c[2].y, c[3].x, c[3].y);
+  line(c[1].x, c[1].y, c[5].x, c[5].y);
+  line(c[2].x, c[2].y, c[6].x, c[6].y);
+  line(c[3].x, c[3].y, c[7].x, c[7].y);
+  line(c[4].x, c[4].y, c[5].x, c[5].y);
+  line(c[5].x, c[5].y, c[6].x, c[6].y);
+  line(c[6].x, c[6].y, c[7].x, c[7].y);
+  line(c[7].x, c[7].y, c[4].x, c[4].y);
+}
+
+getRMF(t, originalPoints) {
+  const curve = new Bezier(this, originalPoints);
+  const d1curve = new Bezier(this, curve.dpoints[0]);
+
+  if (!this.rmf_LUT) {
+    this.rmf_LUT = this.generateRMF(originalPoints, curve, d1curve);
+  }
+
+  // find the frame for "t".
+  const last = this.rmf_LUT.length - 1;
+  const f =  t * last;
+  const i = Math.floor(f);
+
+  // intenger index, or last index: we're done.
+  if (f === i) return this.rmf_LUT[i];
+
+  // no integer index: interpolate based on the adjacent frames.
+  const j = i + 1, ti = i/last, tj = j/last, ratio = (t - ti) / (tj - ti);
+  return this.lerpFrames(ratio, this.rmf_LUT[i], this.rmf_LUT[j]);
+}
+
+generateRMF(originalPoints, curve, d1curve) {
+  const frames = []
+  frames.push(this.getFrenetVectors(0, originalPoints));
+
+  for(let i=0, steps=24; i<steps; i++) {
+    const x0 = frames[i],
+          // get the next frame
+          t = (i+1)/steps,
+          x1 = {
+            o: curve.get(t),
+            dt: d1curve.get(t)
+          },
+          // then mirror the rotational axis and tangent
+          v1 = vec.minus(x1.o, x0.o),
+          c1 = vec.dot(v1, v1),
+          riL = vec.minus(x0.r, vec.scale(v1, 2/c1 * vec.dot(v1, x0.r))),
+          dtiL = vec.minus(x0.dt, vec.scale(v1, 2/c1 * vec.dot(v1, x0.dt))),
+          // and use those to compute a more stable rotational axis
+          v2 = vec.minus(x1.dt, dtiL),
+          c2 = vec.dot(v2, v2);
+    x1.r = vec.minus(riL, vec.scale(v2, 2/c2 * vec.dot(v2, riL)));
+    // and with that stable axis, a new normal.
+    x1.n = vec.cross(x1.r, x1.dt);
+    frames.push(x1);
+  }
+
+  return frames;
+}
+
+getFrenetVectors(t, originalPoints) {
+  const curve = new Bezier(this, originalPoints);
+  const d1curve = new Bezier(this, curve.dpoints[0]);
+  const dt = d1curve.get(t);
+  const ddt = d1curve.derivative(t);
+  const o = curve.get(t);
+  const b = vec.normalize(vec.plus(dt, ddt));
+  const r = vec.normalize(vec.cross(b, dt));
+  const n = vec.normalize(vec.cross(r, dt));
+  return { o, dt, r, n };
+}
+
+
+lerpFrames(t, f1, f2) {
+  var frame = {};
+  [`o`, `dt`, `r`, `n`].forEach(type => frame[type] = vec.lerp(t, f1[type], f2[type]));
+  return frame;
+}
+
+drawVector(from, vec, length, color, label) {
+  setStroke(color);
+  setFill(`black`);
+
+  let pt = project({
+    x: length * vec.x,
+    y: length * vec.y,
+    z: length * vec.z
+  });
+  line(from.x, from.y, from.x + pt.x, from.y + pt.y);
+
+  let txt = project({
+    x: (length+15) * vec.x,
+    y: (length+15) * vec.y,
+    z: (length+15) * vec.z
+  });
+  text(label, from.x + txt.x, from.y + txt.y);
+}
+
+onMouseMove() {
+  this.t = constrain(
+    map(this.cursor.x,0,this.width,-0.1, 1.1)
+    ,0,1
+  );
+  redraw();
+}
+
+onKeyDown() {
+  let key = this.keyboard.currentKey;
+  if (key === `ArrowUp`) {
+    this.t += 0.01;
+    if (this.t > 1) this.t = 1;
+  }
+  if (key === `ArrowDown`) {
+    this.t -= 0.01;
+    if (this.t < 0) this.t = 0;
+  }
+  redraw();
+}

docs/chapters/pointvectors3d/vector-lib.js

@@ -0,0 +1,71 @@
+function normalize(v) {
+    let z = v.z || 0;
+    var d = Math.sqrt(v.x*v.x + v.y*v.y + z*z);
+    let r = { x:v.x/d, y:v.y/d };
+    if (v.z !== undefined) r.z = z/d;
+    return r;
+}
+
+function dot(v1, v2) {
+    let z1 = v1.z || 0;
+    let z2 = v2.z || 0;
+    return v1.x * v2.x + v1.y * v2.y + z1 * z2;
+}
+
+function scale(v, s) {
+    let r = {
+        x: s * v.x,
+        y: s * v.y
+    }
+    if (v.z !== undefined) {
+        r.z = s * v.z
+    }
+    return r;
+}
+
+function plus(v1, v2) {
+    let r = {
+        x: v1.x + v2.x,
+        y: v1.y + v2.y
+    };
+    if (v1.z !== undefined || v2.z !== undefined) {
+        r.z = (v1.z||0) + (v2.z||0);
+    };
+    return r;
+}
+
+function minus(v1, v2) {
+    let r = {
+        x: v1.x - v2.x,
+        y: v1.y - v2.y
+    };
+    if (v1.z !== undefined || v2.z !== undefined) {
+        r.z = (v1.z||0) - (v2.z||0);
+    };
+    return r;
+}
+
+function cross(v1, v2) {
+    if (v1.z === undefined || v2.z === undefined) {
+        throw new Error(`Cross product is not defined for 2D vectors.`);
+    }
+    return {
+        x: v1.y * v2.z - v1.z * v2.y,
+        y: v1.z * v2.x - v1.x * v2.z,
+        z: v1.x * v2.y - v1.y * v2.x
+    };
+}
+
+function lerp(t, v1, v2) {
+    let r = {
+        x: (1-t)*v1.x + t*v2.x,
+        y: (1-t)*v1.y + t*v2.y
+    };
+    if (v1.z !== undefined || v2.z !== undefined) {
+        r.z = (1-t)*(v1.z||0) + t*(v2.z||0);
+    };
+    return r;
+}
+
+
+export default { normalize, dot, scale, plus, minus, cross, lerp }