2d:null 2e:["$","title",null,{"children":"Cubic basis vs. Hermite interpolation | Ricky Reusser | Observable"}] 33:T1cde,{ let selected = mutable points[0]; const xScale = d3.scaleLinear().domain([-0.5, 5.5]).range([0, width]); const yScale = d3.scaleLinear().domain([-2, 5]).range([height, 0]); const svg = d3 .create("svg") .attr("viewBox", [-14, 0, width + 28, height]) .attr("tabindex", 1) .attr("pointer-events", "all") .call( d3 .drag() .subject(dragsubject) .on("start", dragstarted) .on("drag", dragged) ); svg .append("rect") .attr("fill", "none") .attr("width", width) .attr("height", height); svg .append("path") .attr("class", "nearest1") .attr("fill", "none") .attr("stroke", "grey") .attr("stroke-dasharray", "4 4") .attr("stroke-width", 1); svg .append("path") .attr("class", "nearest2") .attr("fill", "none") .attr("stroke", "black") .attr("stroke-width", 2); svg .append("path") .attr("class", "cubic") .attr("fill", "none") .attr("stroke", "#38f") .attr("stroke-width", 3); svg .append("path") .attr("class", "monotone") .attr("fill", "none") .attr("stroke", "#3d8") .attr("stroke-width", 3); svg .append("path") .attr("class", "cubicspline") .attr("fill", "none") .attr("stroke", "#f38") .attr("stroke-width", 3); if (~show.indexOf("derivatives")) { svg .append("path") .attr("class", "axis") .attr("fill", "none") .attr("stroke", "gray") .attr("stroke-dasharray", "2 2") .attr("stroke-width", 2) .datum([ [-10, 0], [10, 0] ]) .attr( "d", d3 .line() .x((d) => xScale(d[0])) .y((d) => yScale(d[1])) ); svg .append("path") .attr("class", "cubic-slope") .attr("fill", "none") .attr("stroke", "#38f") .attr("stroke-dasharray", "2 2") .attr("stroke-width", 2); svg .append("path") .attr("class", "monotone-slope") .attr("fill", "none") .attr("stroke", "#3d8") .attr("stroke-dasharray", "2 2") .attr("stroke-width", 2); svg .append("path") .attr("class", "cubicspline-slope") .attr("fill", "none") .attr("stroke", "#f38") .attr("stroke-dasharray", "2 2") .attr("stroke-width", 2); } update(); function update() { updateNearest(); updateCubic(); const circle = svg.selectAll("g").data(mutable points); circle .enter() .append("g") .call((g) => g.append("circle").attr("r", 30).attr("fill", "none")) .call((g) => g .append("circle") .attr("r", 0) .attr("stroke", "black") .attr("stroke-width", 1.5) .transition() .duration(750) .ease(d3.easeElastic) .attr("r", 5) ) .merge(circle) .attr("transform", (d) => `translate(${xScale(d[0])},${yScale(d[1])})`) .select("circle:last-child") .attr("fill", (d) => (d === selected ? "lightblue" : "black")); circle.exit().remove(); } function updateNearest() { if (~show.indexOf("Nearest")) { svg .select("path.nearest1") .datum( mutable points .map((p) => [ [p[0] - 0.5, p[1]], [p[0] + 0.5, p[1]] ]) .flat() ) .attr( "d", d3 .line() .x((d) => xScale(d[0])) .y((d) => yScale(d[1])) .defined((d) => !!d) ); svg .select("path.nearest2") .datum( mutable points .map((p) => [ [p[0] - 0.5, p[1]], [p[0] + 0.5, p[1]] ]) .map((seg) => [seg, null]) .flat() .flat() ) .attr( "d", d3 .line() .x((d) => xScale(d[0])) .y((d) => yScale(d[1])) .defined((d) => !!d) ); } } function evalCubic(x, basisfunc) { const p = mutable points; const i = [-1, 0, 1, 2].map((i) => Math.min(Math.max(0, Math.floor(x + i)), mutable points.length - 1) ); const basis = basisfunc((x + 1) % 1); return ( basis[0] * p[i[0]][1] + basis[1] * p[i[1]][1] + basis[2] * p[i[2]][1] + basis[3] * p[i[3]][1] ); } function evalMonotone(x, basisfunc) { const p = mutable points; const i = [-1, 0, 1, 2].map((i) => Math.min(Math.max(0, Math.floor(x + i)), mutable points.length - 1) ); return monotoneCubic( p[i[0]][1], p[i[1]][1], p[i[2]][1], p[i[3]][1], (x + 1) % 1 ); } function deriv(f, x, h = 0.01) { return (f(x + h) - f(x - h)) / (2 * h); } function updateCubic() { if (~show.indexOf("B-spline")) { svg .select("path.cubicspline") .datum(linspace(...xScale.domain())) .attr( "d", d3 .line() .x((d) => xScale(d)) .y((d) => yScale(evalCubic(d, (u) => bsplineBasis({ u })))) ); svg .select("path.cubicspline-slope") .datum(linspace(...xScale.domain())) .attr( "d", d3 .line() .x((d) => xScale(d)) .y((d) => yScale(deriv((x) => evalCubic(x, (u) => bsplineBasis({ u })), d)) ) ); } if (~show.indexOf("Monotonic")) { svg .select("path.monotone") .datum(linspace(...xScale.domain())) .attr( "d", d3 .line() .x((d) => xScale(d)) .y((d) => yScale(evalMonotone(d))) ); svg .select("path.monotone-slope") .datum(linspace(...xScale.domain())) .attr( "d", d3 .line() .x((d) => xScale(d)) .y((d) => yScale(deriv((x) => evalMonotone(x), d))) ); } if (~show.indexOf("Hermite")) { svg .select("path.cubic") .datum(linspace(...xScale.domain())) .attr( "d", d3 .line() .x((d) => xScale(d)) .y((d) => yScale(evalCubic(d, (u) => hermiteBasis({ u })))) ); svg .select("path.cubic-slope") .datum(linspace(...xScale.domain())) .attr( "d", d3 .line() .x((d) => xScale(d)) .y((d) => yScale(deriv((x) => evalCubic(x, (u) => hermiteBasis({ u })), d)) ) ); } } function dragsubject(event) { let subject = event.sourceEvent.target.__data__; return subject; } function dragstarted({ subject }) { selected = subject; update(); } function dragged(event) { event.subject[1] = yScale.invert(event.y); mutable points = mutable points; update(); } return html`

${svg.node()}

Three styles of cubic interpolation: cubic basis which undershoots the data, Hermite which overshoots the data, and monotone cubic which clips the slopes of Hermite interpolation to enforce monotonicity and offer a compromise.

`; }34:[] 31:{"id":"fab87daf989e8adb","slug":"cubic-basis-vs-hermite-interpolation","title":"Cubic basis vs. Hermite interpolation","thumbnail":"63bb8630b777bb2df8cee9c64ce2c5219041c1a5196ba514ccb2cfa6a2f0a2f3","access_level":"public","create_time":"$D2023-05-10T15:04:14.080Z","update_time":"$D2024-12-18T17:37:01.614Z","view_time":null,"first_public_time":"$D2023-06-25T15:38:44.301Z","has_unpublished_edits":null,"fork":false,"stars":12,"starred":false,"owner":{"type":"team","id":"f1b8352b963a77b4","login":"rreusser","name":"Ricky Reusser","avatar_url":"https://avatars-next.observableusercontent.com/49c18a2ac0092176d4ea8348a1dd538e0b460b086b89d602b8d21ffe0821d92e","tier":"starter_2024"},"version":874,"publish_version":null,"role":null,"authored":false,"authors":[{"type":"individual","id":"908c4be2279a714e","login":"rreusser","name":"Ricky Reusser","avatar_url":"https://avatars-next.observableusercontent.com/fa41989f47bb986723af5b9c22176cb503748426f56e565000ea93acd9eab86c","tier":null}],"members":[{"type":"individual","id":"908c4be2279a714e","login":"rreusser","name":"Ricky Reusser","avatar_url":"https://avatars-next.observableusercontent.com/fa41989f47bb986723af5b9c22176cb503748426f56e565000ea93acd9eab86c","tier":null,"role":"editor"}],"comments":[],"custom_thumbnail":false,"trashed":false,"lapsed":null,"subscribed":null,"forks":1,"collection":{"id":"dcd095ad05c2b722","slug":"writeups","title":"Fun Notebooks","owner":{"id":"f1b8352b963a77b4","login":"rreusser","name":"Ricky Reusser"}},"type":"notebook","body":{"title":"Cubic basis vs. Hermite interpolation","cells":[{"id":237,"value":"# Cubic basis vs. Hermite interpolation","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":361,"value":"The diagram below illustrates three different methods of cubic spline interpolation which in the style of D3 I'll call cubic basis and Hermite splines. (In addition, I've tested out _monotone_ Hermite interpolation.)\n\nClick to plot the derivatives as well, which provides some interesting insight to just how smooth each method is and how monotone interpolation works.","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":15,"value":"$33","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":640,"value":"viewof show = Inputs.checkbox(\n [\"Nearest\", \"B-spline\", \"Hermite\", \"Monotonic\", \"derivatives\"],\n {\n value: [\"Nearest\", \"B-spline\", \"Hermite\", \"Monotonic\"]\n }\n)","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":352,"value":"Cubic basis spline interpolation uses data as spline control points. A simple method for computing it is repeated linear interpolation, referred to as [De Casteljau's algorithm](https://en.wikipedia.org/wiki/De_Casteljau%27s_algorithm). Since the weight of points in each linear interpolation is in the range ${tex`[0, 1]`}, the weight of points in the repeated sequence of linear interpolations must also be in the range ${tex`[0, 1]`}. This form lends itself well to efficient interpolation using built-in linear filtering capability of GPU hardware. This method is described in *[GPU Gems Chapter 20: Fast Third-Order Texture Filtering](https://developer.nvidia.com/gpugems/gpugems2/part-iii-high-quality-rendering/chapter-20-fast-third-order-texture-filtering)* and tested out in the notebook [Bicubic Texture Interpolation using Linear Filtering](https://observablehq.com/@rreusser/bicubic-texture-interpolation-using-linear-filtering).","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":747,"value":"The four basis function weights in the interval ${tex`[0, 1]`} are","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":708,"value":"viewof bsplineBasis = eqn`f = [\n 1/6 * (1 - u)**3,\n 1/6 * (3 * u**3 - 6 * u**2 + 4),\n 1/6 * (-3 * u**3 + 3 * u**2 + 3 * u + 1),\n 1/6 * u**3\n]`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":263,"value":"{\n const x = linspace(0, 1);\n const basis = bsplineBasis.tabulate({ u: [0, 1] }, 101);\n return htl.html`

\n ${Plot.plot({\n width: 384,\n height: 256,\n grid: true,\n x: { domain: [0, 1] },\n y: { domain: [-0.1, 1] },\n marks: [\n Plot.ruleX([0, 1]),\n Plot.ruleY([0]),\n Plot.line(basis, { x: \"u\", y: \"0\", stroke: \"#d38\", strokeWidth: 2 }),\n Plot.line(basis, { x: \"u\", y: \"1\", stroke: \"#8d3\", strokeWidth: 2 }),\n Plot.line(basis, { x: \"u\", y: \"2\", stroke: \"#38d\", strokeWidth: 2 }),\n Plot.line(basis, { x: \"u\", y: \"3\", stroke: \"#d83\", strokeWidth: 2 })\n ]\n })}\n

Basis functions for B-spline interpolation. Note the maximum value of 2/3, which results in undershoot.

\n `;\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":795,"value":"html`

` /*tex`f(u) = \\frac{1}{6}\n\\begin{bmatrix}\n(1-u)^3 \\\\\n3u^3 - 6u^2 + 4 \\\\\n-3u^3 + 3u^2 + 3u + 1 \\\\\nu^3\n\\end{bmatrix}^T\n\\begin{bmatrix} f_0 \\\\ f_1 \\\\ f_2 \\\\ f_3 \\end{bmatrix}\n`*/","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":364,"value":"Hermite interpolation represents a cubic spline with the derivative specified at two endpoints through which it passes. For one-dimensional interpolation, this derivative can be computed as a simple finite difference of adjacent points. This method results in negative contribution of points, which results in overshoot.","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":793,"value":"The four basis function weights in the interval ${tex`[0, 1]`} are","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":765,"value":"viewof hermiteBasis = eqn`f = [\n 1/2 * (-(u^3) + 2*u^2 - u),\n 1/2 * (3 * u^3 - 5 * u^2 + 2),\n 1/2 * (-3 * u^3 + 4 * u^2 + u),\n 1/2 * (u^3 - u^2)\n]`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":286,"value":"{\n const x = linspace(0, 1);\n const basis = hermiteBasis.tabulate({ u: [0, 1] }, 101);\n return htl.html`

${Plot.plot({\n width: 384,\n height: 256,\n grid: true,\n x: { domain: [0, 1] },\n y: { domain: [-0.1, 1] },\n marks: [\n Plot.ruleX([0, 1]),\n Plot.ruleY([0]),\n Plot.line(basis, { x: \"u\", y: \"0\", stroke: \"#d38\", strokeWidth: 2 }),\n Plot.line(basis, { x: \"u\", y: \"1\", stroke: \"#8d3\", strokeWidth: 2 }),\n Plot.line(basis, { x: \"u\", y: \"2\", stroke: \"#38d\", strokeWidth: 2 }),\n Plot.line(basis, { x: \"u\", y: \"3\", stroke: \"#d83\", strokeWidth: 2 })\n ]\n })}\n

Basis functions for Hermite interpolation. Note the slightly negative weights which result in overshoot.

\n `;\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":814,"value":"html`

` /*tex`f(u) = \\frac{1}{2}\n\\begin{bmatrix}\n-u^3 + 2u^2 - u \\\\\n3u^3 - 5u^2 + 2 \\\\\n-3u^3 + 4u^2 + u \\\\\nu^3 - u^2\n\\end{bmatrix}^T\n\\begin{bmatrix} f_0 \\\\ f_1 \\\\ f_2 \\\\ f_3 \\end{bmatrix}\n`*/","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":632,"value":"Finally, there's one more option using the concepts from [`curveMonotoneX`](https://github.com/d3/d3-shape#curveMonotoneX) in d3. It's a Hermite spline (passes through endpoints with a derivative) but clips the slope to ensure that local extrema may only occur at the points. I'm holding out this will be a good balance between cubic basis splines which are too fuzzy and Hermite splines which overshoot the data.","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":822,"value":"tex`\n\\begin{aligned}\ns_0 &= f_1 - f_0 \\\\\ns_1 &= f_2 - f_1 \\\\\nf'_1 &= (\\operatorname{sign}(s_0) + \\operatorname{sign}(s_1))\n\\operatorname{min}\\left(\\left|s_0\\right|, \\left|s_1\\right|, \\left|\\frac{s_0 + s_1}{4}\\right|\\right)\n\\end{aligned}`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":853,"value":"tex`\n\\begin{aligned}\n\\Delta &= f_2 - f_1 \\\\\nf(u) &= f_1 + f'_1 u + (3\\Delta - 2f'_1 - f'_2) u^2 + (f'_1 + f'_2 - 2\\Delta) u^3\n\\end{aligned}`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":389,"value":"monotoneCubic = {\n function sign(x) {\n return x < 0 ? -1 : 1;\n }\n // Calculate the slopes of the tangents (Hermite-type interpolation) based on\n // the following paper: Steffen, M. 1990. A Simple Method for Monotonic\n // Interpolation in One Dimension. Astronomy and Astrophysics, Vol. 239, NO.\n // NOV(II), P. 443, 1990.\n function slope3(y0, y1, y2) {\n var s0 = y1 - y0;\n var s1 = y2 - y1;\n return (\n (sign(s0) + sign(s1)) *\n Math.min(Math.abs(s0), Math.abs(s1), 0.25 * Math.abs(s0 + s1))\n );\n }\n return function monotoneCubic(f0, f1, f2, f3, u) {\n let d1 = slope3(f0, f1, f2);\n let d2 = slope3(f1, f2, f3);\n let D = f2 - f1;\n return (\n f1 + d1 * u + (3 * D - 2 * d1 - d2) * u ** 2 + (d1 + d2 - 2 * D) * u ** 3\n );\n };\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":0,"value":"mutable points = [\n [0, 0],\n [1, 4],\n [2, 0],\n [3, 0],\n [4, 2],\n [5, 4]\n]","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":48,"value":"function linspace(a, b, n = 200) {\n const arr = new Array(n);\n for (let i = 0; i < n; i++) {\n arr[i] = a + ((b - a) * i) / (n - 1);\n }\n return arr;\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":13,"value":"height = 400","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":10,"value":"d3 = require(\"d3@7\")","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":221,"value":"import { html } from \"@observable/htl\"","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":684,"value":"import { eqn } from \"@rreusser/eqn\"","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"}],"files":"$W34","stdlib":"1","titling":"auto"}} 32:"$31:body"