🏠 back to Observable

Obtain interpolated y-values without drawing a line

#1

I’m looking for a function that would use one of d3’s interpolate functions (ex. cardinal) to return a y value given a set of x, y data and a specific x value within the domain but not already in the data set.

We can do this using path.getPointAtLength but it requires an SVG path element. I don’t want to actually draw the line or rely on SVG elements.

Reference on getPointAtLength: https://stackoverflow.com/questions/11503151/in-d3-how-to-get-the-interpolated-line-data-from-a-svg-line/39442651#39442651

1 Like

#2

You can find the source here, but note that it’s all dependent on the CanvasContextRendering2D API (context.bezierCurveTo):

And so you’ll need something that can evaluate cubic Bézier curves if you want to extract points without using SVG or Canvas. (And getting the value of y given x will be difficult in general with piecewise cubic Bézier curves, as you’ll need to compute the intersection of a line at the given x with each curve segment.) Bezier.js is a library that will do this (curve.intersects), but you’ll need to do some plumbing to connect the output of d3-shape’s curves to Bezier.js curve objects.

0 Likes

#3

This is a bit limited (it only interpolates the d3.curveNatural curve) but it might be good enough for your purposes:

Based on the `cubic-spline` package.

Edit: I thought that the spline agreed with d3.curveNatural, but that was wrong! I’ve updated the notebook so that you can drag the data points and compare the two spline curves.

0 Likes

#4

@mike Thanks for this confirmation that there doesn’t appear to be an existing implementation to this problem (apart from using SVG/Canvas). I’ve come across it several times in the last few years, but just opted for a linear approach for simplicity. I had hoped it would be easier to take advantage of the existing interpolaters without without the dependencies on other libraries.

0 Likes

#5

Thanks @bgchen I’ll take a look at this approach.

0 Likes

#6

This one is using an O(n^3) algorithm to solve a tridiagonal matrix, which can be done in O(n) time. For a small number of points on fast modern computers it will be okay, but it still makes me sad.

I’m working on a notebook explaining some of the math behind cubic splines. Code isn’t there yet, and the explanation isn’t finished either, but someone might find it useful already. https://observablehq.com/d/fb4d519e95cdbd83

2 Likes

#7

As for what D3’s “natural” spline is doing: it is approximating a 2-dimensional function where we know the values at t = [0, 1, 2, 3, …, n], which are [(x(0), y(0)), (x(1), y(1)), … (x(n), y(n))]. It is solving a tridiagonal matrix to find each of the two functions x(t) and y(t) which interpolate the given points.

From what I understand this can work somewhat better if the values of t are based on the chord length between points instead of just being uniform, at the expense of needing to compute an extra inverse square root per interval to scale the coefficients to match the chord length. I will try implementing that one at some point.

Edit: if the points are not equally spaced there is quite a dramatic difference between uniform vs. chord-length parametrization:

If what you are trying to approximate is a 1-dimensional function y(x), it’s generally not the best idea to use a method based on finding (x(t), y(t)), though obviously in many cases it works out more or less fine.

@john-clarke if you can be a bit more specific about your context / use case, there are many different methods out there of interpolating data. This has been a subject of active research for centuries now, and there is a significant amount of theory and many practical implementations.

3 Likes

#8

I still have a ways to go with more explanatory diagrams but I added an amusing picture to the top of

… so it’s not just formulas anymore. Warning: moving the points around and watching the little rainbow lines move in response is kind of addicting, and might distract you from whatever work you were otherwise trying to do.


@john-clarke

Try this self-contained version:

// Given an array of knots and an array of values,
// return a function which evaluates a natural cubic spline
// extrapolated by linear functions on either end.
NCS = function NCS(knots, values) {
  const n = knots.length - 1;
  const t = new Float64Array(n+3); t.set(knots, 1);
  const y = new Float64Array(n+3); y.set(values, 1);
  const s = new Float64Array(n+2); // horizontal scales
  const m = new Float64Array(n+3); // slope at each knot
  const d = new Float64Array(n+3); // diagonal matrix
  
  // Natural cubic spline algorithm
  for (let i = 1; i < n+1; i++) {
    s[i] = 1 / (t[i+1] - t[i]);
    m[i] += (m[i+1] = 3*s[i]*s[i]*(y[i+1] - y[i])); }
  d[1] = 0.5 / s[1];
  m[1] = d[1] * m[1];
  for (let i = 2; i <= n+1; i++) {
    d[i] = 1 / (2*(s[i] + s[i-1]) - s[i-1]*s[i-1]*d[i-1]);
    m[i] = d[i]*(m[i] - s[i-1]*m[i-1]); }
  for (let i = n; i >= 1; i--) {
    m[i] -= d[i]*s[i]*m[i+1]; }

  // Linear extrapolation
  t[0] = t[1] - 1; t[n+2] = t[n+1] + 1;
  y[0] = y[1] - m[1]; y[n+2] = y[n+1] + m[n+1];
  s[0] = s[n+1] = 1;
  m[0] = m[1]; m[n+2] = m[n+1];
  
  // Set up Chebyshev coefficients
  const coeffs = new Float64Array(4*(n+2));
  for (let i = 0; i < n+2; i++) {
    const h = t[i+1] - t[i];
    const y0 = y[i], y1 = y[i+1], m0 = h*m[i], m1 = h*m[i+1], j = 4*i;
    coeffs[j] = 0.5*(y0 + y1 + 0.125*(m0 - m1));
    coeffs[j+1] = 0.03125*(18*(y1 - y0) - m0 - m1);
    coeffs[j+2] = 0.0625*(m1 - m0);
    coeffs[j+3] = 0.03125*(2*(y0 - y1) + m0 + m1); }

  return Object.assign(
    function eval_spline(u) {
      // Binary search
      let nn = t.length - 1, low = 0, half;
      while ((half = (nn/2)|0) > 0) {
        low += half * (t[low + half] <= u);
        nn -= half; }
      const i = low, j = 4*i;
      u = (2*u - t[i] - t[i+1]) * s[i]; // scale to [–1, 1]
      // Clenshaw's method.
      let b1 = coeffs[j+3], b2 = coeffs[j+2] + 2 * u * b1;
      b1 = coeffs[j+1] + 2 * u * b2 - b1;
      return coeffs[j] + u * b1 - b2; },
    { knots: t, values: y, dvalues: m, scales: s, coeffs: coeffs });
}

Example usage:

{
  const
    t = [1, 2.5, 3, 4, 5.5, 6, 7, 7.2, 9],
    y = [1, 1.5, 1, 2, 1, 1.5, 1.3, 1.1, 1],
    spline = NCS(t, y),
    x = [...Array(500).keys()].map(i=> i/50),
    div = DOM.element('div');

  Plotly.newPlot(div, [
    { x: x, y: x.map(spline), mode: 'lines', type: 'scatter'},
    { x: t, y: y, mode: 'markers', type: 'scatter'}
  ], {width: width});
  return div;
}
2 Likes