Obtain interpolated y-values without drawing a line

Community Help
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

2 Likes
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.

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.

1 Like
4

@mbostock 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.

5

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

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:

23%20PM
Screen Shot 2019-04-10 at 11.02.23 PM.png1680×1148 88.9 KB

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
9

Here is another attempt, targeting a d3 free, generic solution, based on .getPointAtLength():

and imported here to demonstrate potential use:

Although SVG based, it does not require drawing anything. It does generate a svg element but that is hidden/detached. It just uses the path string which can be easily generated with d3 and typically would be anyways.

It is a little expensive since it is iterative (although 20 iteration or so are generally enough per interpolation) but should work with any SVG path, eg. any d3 curve.

It does not attempt to somehow deal with multivalued functions which should not be the common case and was not initially considered here. It will just return the first point on the path it finds.
The function to import is d3 free and has no dependencies. It just does simple bisections until it gets close enough to the desired point.

It would be straighforward to add a x(y) interpolator which should be also useful.

Ideas to optimize or any other feedback very welcome.

1 Like
10

Nice work, but I would recommend against this as a general solution; most of the D3 “interpolators” which output (x, y) as a function of t are not designed to guarantee that y is properly a function of x. Using bisection to find a value of t where x matches the target and then reading off y(t) is going to fall over in edge cases (which in practice occur quite commonly in real-world D3 charts), leaving the user mystified when the method breaks.

Folks who have (x, y) pairs intended to lie on a function y(x) should use a method designed for functional interpolation.

When we get right down to it, folks making statistical graphics should very often also be using an approximation for y(x) instead of (x(t), y(t)). Or for example when making a “radar chart”, the outline should be created via an approximation of a function r(θ), instead of (x(t), y(t)). Tools like D3 and Vega should be providing a wider variety of tools for plotting approximations to functions instead of only 2D curve fitting.

11

The background for me was to find a way to easily use d3 curves in a renderer which is not bezier or spline aware, but webgl based, for pleasing, smooth lines such as in this

Out of the curves, monotoneX seems the most suitable.

Since svg already does the spline computation, probably faster than js could, why not rake advantage ?

I hope svg properly calculates path length for interrupted paths. If so, I should add a provision to break out of the search early if there are no improvements even if the target is not sufficiently close.

Another idea may be to accept an initial hint for y, and do more of a grid search.

Actually, it may work pretty well to just use the rasterized line as a lookup grid, on a sufficiently large canvas. And perhaps refine on a zoomed in grid.

Moving median for despiking would be good for noisy data, polynomial fitting for sparse, kriging for geo-maps, but where do you draw the line and use dedicated libraries/systems ?

12

I added gaps and the idea to exit early if there is no improvement:


But it may be preferable to finish iterating to the maximum allowed iteration and return a null or NaN if no sufficiently close point was found. I guess it would depend on the use case. Presumably, if the sampled line was defined discontinuously in the first place, there would be an expectation to not be able to get interpolated values in the gaps, or at least an understanding that this could happen.

13

The idea to detect lack of improvement was flawed since bisection often does not get closer to a solution. So now the path interpolator just finishes iterating and returns a false value if it does not converge in time. This (usually?) means that the target point on the line is in a gap. In the example this is used to hide the cursor and the interpolated dot in interpolation gaps:

14

Hi all, I have put in a proposal pull request to D3 core which may be of interest here: https://github.com/d3/d3-interpolate/pull/67

This PR proposes 4 new functions:

  • d3.interpolateFromCurve()
  • d3.interpolateCardinal()
  • d3.interpolateCatmullRom()
  • d3interpolateMonotoneX()

It would be awesome, if of interest, to see your comments over on this PR. At the moment it is looking like this will have to be a separate module, but I would like to collect a ‘show of hands’ for people who would find this useful if it were in core D3.

1 Like
15

Further to my last comment, if any use to people on this thread, I have released the above as a standalone D3 module:

Create interpolator from D3 curve functions.

2 Likes