generalised curve fitting

docs/chapters/curvefitting/comments.txt

@@ -1,3 +0,0 @@
-## Additionals
-
-http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.96.5193&rep=rep1&type=pdf ?

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

@@ -251,13 +251,11 @@ Here, the "to the power negative one" is the notation for the [matrix inverse](h
 
 So before we try that out, how much code is involved in implementing this? Honestly, that answer depends on how much you're going to be writing yourself. If you already have a matrix maths library available, then really not that much code at all. On the other hand, if you are writing this from scratch, you're going to have to write some utility functions for doing your matrix work for you, so it's really anywhere from 50 lines of code to maybe 200 lines of code. Not a bad price to pay for being able to fit curves to prespecified coordinates.
 
-So let's try it out! The following graphic lets you place points, and will start computing exact-fit curves once you've placed at least three. You can click for more points, and the code will simply try to compute an exact fit using a Bézier curve of the appropriate order. Four points? Cubic Bézier. Five points? Quartic. And so on. Of course, this does break down at some point: depending on where you place your points, it might become mighty hard for the fitter to find an exact fit, and things might actually start looking horribly off once you hit 10<sup>th</sup> or higher order curves. But it might not!
+So let's try it out! The following graphic lets you place points, and will start computing exact-fit curves once you've placed at least three. You can click for more points, and the code will simply try to compute an exact fit using a Bézier curve of the appropriate order. Four points? Cubic Bézier. Five points? Quartic. And so on. Of course, this does break down at some point: depending on where you place your points, it might become mighty hard for the fitter to find an exact fit, and things might actually start looking horribly off once there's enough points for compound [floating point rounding errors](https://en.wikipedia.org/wiki/Round-off_error#Floating-point_number_system) to start making a difference (which is around 10~11 points).
 
-<div class="figure">
-  <Graphic title="Fitting a Bézier curve" setup={this.setup} draw={this.draw} onClick={this.onClick}>
-    <button onClick={this.toggle} style="position:absolute; right: 0;">toggle</button>
-    <SliderSet ref={ set => (this.sliders=set) } onChange={this.processTimeUpdate} />
-  </Graphic>
-</div>
+<graphics-element title="Fitting a Bézier curve" width="550" src="./curve-fitting.js" >
+  <button class="toggle">toggle</button>
+  <div class="sliders"><!-- this will contain range inputs, created by the graphic --></div>
+</graphics-element>
 
-You'll note there is a convenient "toggle" buttons that lets you toggle between equidistance `t` values, and distance ratio along the polygon. Arguably more interesting is that once you have points to abstract a curve, you also get <em>direct control</em> over the time values through sliders for each, because if the time values are our degree of freedom, you should be able to freely manipulate them and see what the effect on your curve is.
+You'll note there is a convenient "toggle" buttons that lets you toggle between equidistant `t` values, and distance ratio along the polygon formed by the points. Arguably more interesting is that once you have points to abstract a curve, you also get <em>direct control</em> over the time values through sliders for each, because if the time values are our degree of freedom, you should be able to freely manipulate them and see what the effect on your curve is.

docs/chapters/curvefitting/curve-fitter.js

@@ -0,0 +1,148 @@
+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

@@ -0,0 +1,113 @@
+import fit from "./curve-fitter.js";
+
+let points = [], curve, sliders;
+
+setup() {
+  let btn = find(`.toggle`);
+  if (btn) btn.listen(`click`, evt => this.toggle());
+  sliders = find(`.sliders`);
+  this.mode = 0;
+  this.label = `Using equidistant t values`;
+}
+
+toggle() {
+  this.mode = (this.mode + 1) % 2;
+  if (sliders) this.setSliderValues(this.mode);
+  redraw();
+}
+
+draw() {
+  clear();
+  setColor('lightgrey');
+  drawGrid(10);
+
+  setColor('black');
+  setFontSize(16);
+  setTextStroke(`white`, 4);
+  if (points.length > 2) {
+    curve = this.fitCurve(points);
+    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]}
+      ));
+  return new Bezier(this, bpoints);
+}
+
+updateSliders() {
+  if (sliders && points.length > 2) {
+    sliders.innerHTML = ``;
+    sliders.values = [];
+    this.sliders = points.map((p,i) => {
+      // TODO: this should probably be built into the graphics API as a
+      //       things that you can do, e.g. clearSliders() and addSlider()
+      let s = document.createElement(`input`);
+      s.setAttribute(`type`, `range`);
+      s.setAttribute(`min`, `0`);
+      s.setAttribute(`max`, `1`);
+      s.setAttribute(`step`, `0.01`);
+      s.classList.add(`slide-control`);
+      sliders.values[i] =  i/(points.length-1);
+      s.setAttribute(`value`, sliders.values[i]);
+      s.addEventListener(`input`, evt => {
+        this.label = `Using custom t values`;
+        sliders.values[i] = parseFloat(evt.target.value);
+        redraw();
+      });
+      sliders.append(s);
+    });
+  }
+}
+
+setSliderValues(mode) {
+  let n = points.length;
+
+  // equidistant
+  if (mode === 0) {
+    this.label = `Using equidistant t values`;
+    sliders.values = [...new Array(n)].map((_,i) =>i/(n-1));
+  }
+
+  // polygonal distance
+  if (mode === 1) {
+    this.label = `Using polygonal distance t values`;
+    const D = [0];
+    for(let i = 1; i<n; i++) {
+      D[i] = D[i-1] + dist(
+        points[i-1].x, points[i-1].y,
+        points[i].x, points[i].y
+      );
+    }
+    const S = [], len = D[n-1];
+    D.forEach((v,i) => { S[i] = v/len; });
+    sliders.values = S;
+  }
+
+  findAll(`.sliders input[type=range]`).forEach((s,i) => {
+    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/handler.js

@@ -1,87 +0,0 @@
-var fit = require('../../../lib/curve-fitter.js');
-
-module.exports = {
-  setup: function(api) {
-    this.api = api;
-    api.noDrag = true; // do not allow points to be dragged around
-    this.reset();
-  },
-
-  reset: function() {
-    this.points = [];
-    this.sliders.setOptions([]);
-    this.curveset = false;
-    this.mode = 0;
-    if (this.api) {
-      let api = this.api;
-      api.setCurve(false);
-      api.reset();
-      api.redraw();
-    }
-  },
-
-  toggle: function() {
-    if (this.api) {
-      this.customTimeValues = false;
-      this.mode = (this.mode + 1) % fit.modes.length;
-      this.fitCurve(this.api);
-      this.api.redraw();
-    }
-  },
-
-  draw: function(api, curve) {
-    api.setPanelCount(1);
-    api.reset();
-    api.setColor('lightgrey');
-    api.drawGrid(10,10);
-
-    api.setColor('black');
-
-    if (!this.curveset && this.points.length > 2) {
-      curve = this.fitCurve(api);
-    }
-
-    if (curve) {
-      api.drawCurve(curve);
-      api.drawSkeleton(curve);
-    }
-
-    api.drawPoints(this.points);
-
-    if (!this.customTimeValues) {
-      api.setFill(0);
-      api.text("using "+fit.modes[this.mode]+" t values", {x: 5, y: 10});
-    }
-  },
-
-  processTimeUpdate(sliderid, timeValues) {
-    var api = this.api;
-    this.customTimeValues = true;
-    this.fitCurve(api, timeValues);
-    api.redraw();
-  },
-
-  fitCurve(api, timeValues) {
-    let bestFitData = fit(this.points, timeValues || this.mode),
-        x = bestFitData.C.x,
-        y = bestFitData.C.y,
-        bpoints = [];
-    x.forEach((r,i) => {
-      bpoints.push({
-        x: r[0],
-        y: y[i][0]
-      });
-    });
-    var curve = new api.Bezier(bpoints);
-    api.setCurve(curve);
-    this.curveset = true;
-    this.sliders.setOptions(bestFitData.S);
-    return curve;
-  },
-
-  onClick: function(evt, api) {
-    this.curveset = false;
-    this.points.push({x: api.mx, y: api.my });
-    api.redraw();
-  }
-};

docs/chapters/curvefitting/matrix-invert.js

@@ -0,0 +1,110 @@
+// 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;
+};