catmull-rom

docs/chapters/curvefitting/curve-fitter.js

@@ -1,148 +0,0 @@
-import invert from "./matrix-invert.js";
-
-var binomialCoefficients = [[1],[1,1]];
-
-function binomial(n,k) {
-  if (n===0) return 1;
-  var lut = binomialCoefficients;
-  while(n >= lut.length) {
-    var s = lut.length;
-    var nextRow = [1];
-    for(var i=1,prev=s-1; i<s; i++) {
-      nextRow[i] = lut[prev][i-1] + lut[prev][i];
-    }
-    nextRow[s] = 1;
-    lut.push(nextRow);
-  }
-  return lut[n][k];
-}
-
-function transpose(M) {
-  var Mt = [];
-  M.forEach(row => Mt.push([]));
-  M.forEach((row,r) => row.forEach((v,c) => Mt[c][r] = v));
-  return Mt;
-}
-
-function row(M,i) {
-  return M[i];
-}
-
-function col(M,i) {
-  var col = [];
-  for(var r=0, l=M.length; r<l; r++) {
-    col.push(M[r][i]);
-  }
-  return col;
-}
-
-function multiply(M1, M2) {
-  // prep
-  var M = [];
-  var dims = [M1.length, M1[0].length, M2.length, M2[0].length];
-  // work
-  for (var r=0, c; r<dims[0]; r++) {
-    M[r] = [];
-    var _row = row(M1, r);
-    for (c=0; c<dims[3]; c++) {
-      var _col = col(M2,c);
-      var reducer = (a,v,i) => a + _col[i]*v;
-      M[r][c] = _row.reduce(reducer, 0);
-    }
-  }
-  return M;
-}
-
-function getValueColumn(P, prop) {
-  var col = [];
-  P.forEach(v => col.push([v[prop]]));
-  return col;
-}
-
-function computeBasisMatrix(n) {
-  /*
-    We can form any basis matrix using a generative approach:
-
-     - it's an M = (n x n) matrix
-     - it's a lower triangular matrix: all the entries above the main diagonal are zero
-     - the main diagonal consists of the binomial coefficients for n
-     - all entries are symmetric about the antidiagonal.
-
-    What's more, if we number rows and columns starting at 0, then
-    the value at position M[r,c], with row=r and column=c, can be
-    expressed as:
-
-      M[r,c] = (r choose c) * M[r,r] * S,
-
-      where S = 1 if r+c is even, or -1 otherwise
-
-    That is: the values in column c are directly computed off of the
-    binomial coefficients on the main diagonal, through multiplication
-    by a binomial based on matrix position, with the sign of the value
-    also determined by matrix position. This is actually very easy to
-    write out in code:
-  */
-
-  // form the square matrix, and set it to all zeroes
-  var M = [], i = n;
-  while (i--) { M[i] = "0".repeat(n).split('').map(v => parseInt(v)); }
-
-  // populate the main diagonal
-  var k = n - 1;
-  for (i=0; i<n; i++) { M[i][i] = binomial(k, i); }
-
-  // compute the remaining values
-  for (var c=0, r; c<n; c++) {
-    for (r=c+1; r<n; r++) {
-      var sign = (r+c)%2 ? -1 : 1;
-      var value = binomial(r, c) * M[r][r];
-      M[r][c] = sign * value; }}
-
-  return M;
-}
-
-function raiseRowPower(row, i) {
-  return row.map(v => Math.pow(v,i));
-}
-
-function formTMatrix(S, n) {
-  n = n || S.length;
-  var Tp = [];
-  // it's easier to generate the transposed matrix:
-  for(var i=0; i<n; i++) Tp.push( raiseRowPower(S, i));
-  return {
-    Tt: Tp,
-    T: transpose(Tp) // and then transpose "again" to get the real matrix
-  };
-}
-
-function computeBestFit(P, M, S, n) {
-  n = n || P.length;
-  var tm = formTMatrix(S, n),
-      T = tm.T,
-      Tt = tm.Tt,
-      M1 = invert(M),
-      TtT1 = invert(multiply(Tt,T)),
-      step1 = multiply(TtT1, Tt),
-      step2 = multiply(M1, step1),
-      X = getValueColumn(P,'x'),
-      Cx = multiply(step2, X),
-      Y = getValueColumn(P,'y'),
-      Cy = multiply(step2, Y);
-  return { x: Cx, y: Cy };
-}
-
-function fit(points, tvalues) {
-  // mode could be an int index to fit.modes, below,
-  // which are used to abstract time values, OR it
-  // could be a prespecified array of time values to
-  // be used in the final curve fitting step.
-  const n = points.length,
-        P = Array.from(points),
-        M = computeBasisMatrix(n),
-        S = tvalues,
-        C = computeBestFit(P, M, S, n);
-  return { n, P, M, S, C };
-}
-
-export default fit;

docs/chapters/curvefitting/curve-fitting.js

@@ -1,5 +1,3 @@
-import fit from "./curve-fitter.js";
-
 let points = [], curve, sliders;
 
 setup() {
@@ -24,29 +22,90 @@ draw() {
   setColor('black');
   setFontSize(16);
   setTextStroke(`white`, 4);
-  if (points.length > 2) {
-    curve = this.fitCurve(points);
+  const n = points.length;
+  if (n > 2 && sliders && sliders.values) {
+    curve = this.fitCurveToPoints(n);
     curve.drawSkeleton(`blue`);
     curve.drawCurve();
-
     text(this.label, this.width/2, 20, CENTER);
   }
 
   points.forEach(p => circle(p.x, p.y, 3));
 }
 
-fitCurve(points) {
-  let n = points.length;
-  let tvalues = sliders ? sliders.values : [...new Array(n)].map((_,i) =>i/(n-1));
-  let bestFitData = fit(points, tvalues),
-      x = bestFitData.C.x,
-      y = bestFitData.C.y,
-      bpoints = x.map((r,i) => (
-        {x: r[0], y: y[i][0]}
-      ));
+fitCurveToPoints(n) {
+  // alright, let's do this thing:
+  const tm = this.formTMatrix(sliders.values, n),
+        T = tm.T,
+        Tt = tm.Tt,
+        M = this.generateBasisMatrix(n),
+        M1 = M.invert(),
+        TtT1 = Tt.multiply(T).invert(),
+        step1 = TtT1.multiply(Tt),
+        step2 = M1.multiply(step1),
+        // almost there...
+        X = new Matrix(points.map((v) => [v.x])),
+        Cx = step2.multiply(X),
+        x = Cx.data,
+        // almost...
+        Y = new Matrix(points.map((v) => [v.y])),
+        Cy = step2.multiply(Y),
+        y = Cy.data,
+        // last step!
+        bpoints = x.map((r,i) => ({x: r[0], y: y[i][0]}));
+
   return new Bezier(this, bpoints);
 }
 
+formTMatrix(row, n) {
+  // it's actually easier to create the transposed
+  // version, and then (un)transpose that to get T!
+  let data = [];
+  for (var i = 0; i < n; i++) {
+    data.push(row.map((v) => v ** i));
+  }
+  const Tt = new Matrix(n, n, data);
+  const T = Tt.transpose();
+  return { T, Tt };
+}
+
+generateBasisMatrix(n) {
+  const M = new Matrix(n, n);
+
+  // populate the main diagonal
+  var k = n - 1;
+  for (let i = 0; i < n; i++) {
+    M.set(i, i, binomial(k, i));
+  }
+
+  // compute the remaining values
+  for (var c = 0, r; c < n; c++) {
+    for (r = c + 1; r < n; r++) {
+      var sign = (r + c) % 2 === 0 ? 1 : -1;
+      var value = binomial(r, c) * M.get(r, r);
+      M.set(r, c, sign * value);
+    }
+  }
+
+  return M;
+}
+
+onMouseDown() {
+  if (!this.currentPoint) {
+    const {x, y} = this.cursor;
+    points.push({ x, y });
+    resetMovable(points);
+    this.updateSliders();
+    redraw();
+  }
+}
+
+
+// -------------------------------------
+// The rest of this code is slider logic
+// -------------------------------------
+
+
 updateSliders() {
   if (sliders && points.length > 2) {
     sliders.innerHTML = ``;
@@ -100,14 +159,4 @@ setSliderValues(mode) {
     s.setAttribute(`value`, sliders.values[i]);
     s.value = sliders.values[i];
   });
-}
-
-onMouseDown() {
-  if (!this.currentPoint) {
-    const {x, y} = this.cursor;
-    points.push({ x, y });
-    resetMovable(points);
-    this.updateSliders();
-    redraw();
-  }
-}
+}

docs/chapters/curvefitting/matrix-invert.js

@@ -1,110 +0,0 @@
-// Copied from http://blog.acipo.com/matrix-inversion-in-javascript/
-
-// Returns the inverse of matrix `M`.
-export default function matrix_invert(M) {
-  // I use Guassian Elimination to calculate the inverse:
-  // (1) 'augment' the matrix (left) by the identity (on the right)
-  // (2) Turn the matrix on the left into the identity by elemetry row ops
-  // (3) The matrix on the right is the inverse (was the identity matrix)
-  // There are 3 elemtary row ops: (I combine b and c in my code)
-  // (a) Swap 2 rows
-  // (b) Multiply a row by a scalar
-  // (c) Add 2 rows
-
-  //if the matrix isn't square: exit (error)
-  if (M.length !== M[0].length) {
-    console.log('not square');
-    return;
-  }
-
-  //create the identity matrix (I), and a copy (C) of the original
-  var i = 0,
-      ii = 0,
-      j = 0,
-      dim = M.length,
-      e = 0,
-      t = 0;
-  var I = [],
-      C = [];
-  for (i = 0; i < dim; i += 1) {
-    // Create the row
-    I[I.length] = [];
-    C[C.length] = [];
-    for (j = 0; j < dim; j += 1) {
-      //if we're on the diagonal, put a 1 (for identity)
-      if (i == j) {
-        I[i][j] = 1;
-      } else {
-        I[i][j] = 0;
-      }
-
-      // Also, make the copy of the original
-      C[i][j] = M[i][j];
-    }
-  }
-
-  // Perform elementary row operations
-  for (i = 0; i < dim; i += 1) {
-    // get the element e on the diagonal
-    e = C[i][i];
-
-    // if we have a 0 on the diagonal (we'll need to swap with a lower row)
-    if (e == 0) {
-      //look through every row below the i'th row
-      for (ii = i + 1; ii < dim; ii += 1) {
-        //if the ii'th row has a non-0 in the i'th col
-        if (C[ii][i] != 0) {
-          //it would make the diagonal have a non-0 so swap it
-          for (j = 0; j < dim; j++) {
-            e = C[i][j]; //temp store i'th row
-            C[i][j] = C[ii][j]; //replace i'th row by ii'th
-            C[ii][j] = e; //repace ii'th by temp
-            e = I[i][j]; //temp store i'th row
-            I[i][j] = I[ii][j]; //replace i'th row by ii'th
-            I[ii][j] = e; //repace ii'th by temp
-          }
-          //don't bother checking other rows since we've swapped
-          break;
-        }
-      }
-      //get the new diagonal
-      e = C[i][i];
-      //if it's still 0, not invertable (error)
-      if (e == 0) {
-        return;
-      }
-    }
-
-    // Scale this row down by e (so we have a 1 on the diagonal)
-    for (j = 0; j < dim; j++) {
-      C[i][j] = C[i][j] / e; //apply to original matrix
-      I[i][j] = I[i][j] / e; //apply to identity
-    }
-
-    // Subtract this row (scaled appropriately for each row) from ALL of
-    // the other rows so that there will be 0's in this column in the
-    // rows above and below this one
-    for (ii = 0; ii < dim; ii++) {
-      // Only apply to other rows (we want a 1 on the diagonal)
-      if (ii == i) {
-        continue;
-      }
-
-      // We want to change this element to 0
-      e = C[ii][i];
-
-      // Subtract (the row above(or below) scaled by e) from (the
-      // current row) but start at the i'th column and assume all the
-      // stuff left of diagonal is 0 (which it should be if we made this
-      // algorithm correctly)
-      for (j = 0; j < dim; j++) {
-        C[ii][j] -= e * C[i][j]; //apply to original matrix
-        I[ii][j] -= e * I[i][j]; //apply to identity
-      }
-    }
-  }
-
-  //we've done all operations, C should be the identity
-  //matrix I should be the inverse:
-  return I;
-};