fix barycentric extrapolation (#1701)

Fix barycentric extrapolation

Note: When points on the hull are collinear, the extrapolation fails in some regions. I've found a way to fix this by reprojecting the points (to make the hull slightly bulge out, which resolves the ties), but in the end it was too much code for a use case that is pretty bad anyway (the delaunay triangulation itself being unstable in that case). This seems to be quite enough.


closes #1700

---------
Co-authored-by: Mike Bostock <[email protected]>

src/marks/raster.js

@@ -1,7 +1,7 @@
 import {blurImage, Delaunay, randomLcg, rgb} from "d3";
 import {valueObject} from "../channel.js";
 import {create} from "../context.js";
-import {map, first, second, third, isTuples, isNumeric, isTemporal, take, identity} from "../options.js";
+import {map, first, second, third, isTuples, isNumeric, isTemporal, identity} from "../options.js";
 import {maybeColorChannel, maybeNumberChannel} from "../options.js";
 import {Mark} from "../mark.js";
 import {applyAttr, applyDirectStyles, applyIndirectStyles, applyTransform, impliedString} from "../style.js";
@@ -282,31 +282,16 @@ export function interpolateNone(index, width, height, X, Y, V) {
 
 export function interpolatorBarycentric({random = randomLcg(42)} = {}) {
   return (index, width, height, X, Y, V) => {
-    // Flatten the input coordinates to prepare to insert extrapolated points
-    // along the perimeter of the grid (so there’s no blank output).
-    const n = index.length;
-    const nw = width >> 2;
-    const nh = (height >> 2) - 1;
-    const m = n + nw * 2 + nh * 2;
-    const XY = new Float64Array(m * 2);
-    for (let i = 0; i < n; ++i) (XY[i * 2] = X[index[i]]), (XY[i * 2 + 1] = Y[index[i]]);
-
-    // Add points along each edge, making sure to include the four corners for
-    // complete coverage (with no chamfered edges).
-    let i = n;
-    const addPoint = (x, y) => ((XY[i * 2] = x), (XY[i * 2 + 1] = y), i++);
-    for (let j = 0; j <= nw; ++j) addPoint((j / nw) * width, 0), addPoint((j / nw) * width, height);
-    for (let j = 0; j < nh; ++j) addPoint(width, (j / nh) * height), addPoint(0, (j / nh) * height);
-
-    // To each edge point, assign the closest (non-extrapolated) value.
-    V = take(V, index);
-    const delaunay = new Delaunay(XY.subarray(0, n * 2));
-    for (let j = n, ij; j < m; ++j) V[j] = V[(ij = delaunay.find(XY[j * 2], XY[j * 2 + 1], ij))];
-
     // Interpolate the interior of all triangles with barycentric coordinates
-    const {points, triangles} = new Delaunay(XY);
-    const W = new V.constructor(width * height);
+    const {points, triangles, hull} = Delaunay.from(
+      index,
+      (i) => X[i],
+      (i) => Y[i]
+    );
+    const W = new V.constructor(width * height).fill(NaN);
+    const S = new Uint8Array(width * height); // 1 if pixel has been seen.
     const mix = mixer(V, random);
+
     for (let i = 0; i < triangles.length; i += 3) {
       const ta = triangles[i];
       const tb = triangles[i + 1];
@@ -323,9 +308,9 @@ export function interpolatorBarycentric({random = randomLcg(42)} = {}) {
       const y2 = Math.max(Ay, By, Cy);
       const z = (By - Cy) * (Ax - Cx) + (Ay - Cy) * (Cx - Bx);
       if (!z) continue;
-      const va = V[ta];
-      const vb = V[tb];
-      const vc = V[tc];
+      const va = V[index[ta]];
+      const vb = V[index[tb]];
+      const vc = V[index[tc]];
       for (let x = Math.floor(x1); x < x2; ++x) {
         for (let y = Math.floor(y1); y < y2; ++y) {
           if (x < 0 || x >= width || y < 0 || y >= height) continue;
@@ -337,14 +322,91 @@ export function interpolatorBarycentric({random = randomLcg(42)} = {}) {
           if (gb < 0) continue;
           const gc = 1 - ga - gb;
           if (gc < 0) continue;
-          W[x + width * y] = mix(va, ga, vb, gb, vc, gc, x, y);
+          const i = x + width * y;
+          W[i] = mix(va, ga, vb, gb, vc, gc, x, y);
+          S[i] = 1;
         }
       }
     }
+    extrapolateBarycentric(W, S, X, Y, V, width, height, hull, index, mix);
     return W;
   };
 }
 
+// Extrapolate by finding the closest point on the hull.
+function extrapolateBarycentric(W, S, X, Y, V, width, height, hull, index, mix) {
+  X = Float64Array.from(hull, (i) => X[index[i]]);
+  Y = Float64Array.from(hull, (i) => Y[index[i]]);
+  V = Array.from(hull, (i) => V[index[i]]);
+  const n = X.length;
+  const rays = Array.from({length: n}, (_, j) => ray(j, X, Y));
+  let k = 0;
+  for (let y = 0; y < height; ++y) {
+    const yp = y + 0.5;
+    for (let x = 0; x < width; ++x) {
+      const i = x + width * y;
+      if (!S[i]) {
+        const xp = x + 0.5;
+        for (let l = 0; l < n; ++l) {
+          const j = (n + k + (l % 2 ? (l + 1) / 2 : -l / 2)) % n;
+          if (rays[j](xp, yp)) {
+            const t = segmentProject(X.at(j - 1), Y.at(j - 1), X[j], Y[j], xp, yp);
+            W[i] = mix(V.at(j - 1), t, V[j], 1 - t, V[j], 0, x, y);
+            k = j;
+            break;
+          }
+        }
+      }
+    }
+  }
+}
+
+// Projects a point p = [x, y] onto the line segment [p1, p2], returning the
+// projected coordinates p’ as t in [0, 1] with p’ = t p1 + (1 - t) p2.
+function segmentProject(x1, y1, x2, y2, x, y) {
+  const dx = x2 - x1;
+  const dy = y2 - y1;
+  const a = dx * (x2 - x) + dy * (y2 - y);
+  const b = dx * (x - x1) + dy * (y - y1);
+  return a > 0 && b > 0 ? a / (a + b) : +(a > b);
+}
+
+function cross(xa, ya, xb, yb) {
+  return xa * yb - xb * ya;
+}
+
+function ray(j, X, Y) {
+  const n = X.length;
+  const xc = X.at(j - 2);
+  const yc = Y.at(j - 2);
+  const xa = X.at(j - 1);
+  const ya = Y.at(j - 1);
+  const xb = X[j];
+  const yb = Y[j];
+  const xd = X.at(j + 1 - n);
+  const yd = Y.at(j + 1 - n);
+  const dxab = xa - xb;
+  const dyab = ya - yb;
+  const dxca = xc - xa;
+  const dyca = yc - ya;
+  const dxbd = xb - xd;
+  const dybd = yb - yd;
+  const hab = Math.hypot(dxab, dyab);
+  const hca = Math.hypot(dxca, dyca);
+  const hbd = Math.hypot(dxbd, dybd);
+  return (x, y) => {
+    const dxa = x - xa;
+    const dya = y - ya;
+    const dxb = x - xb;
+    const dyb = y - yb;
+    return (
+      cross(dxa, dya, dxb, dyb) > -1e-6 &&
+      cross(dxa, dya, dxab, dyab) * hca - cross(dxa, dya, dxca, dyca) * hab > -1e-6 &&
+      cross(dxb, dyb, dxbd, dybd) * hab - cross(dxb, dyb, dxab, dyab) * hbd <= 0
+    );
+  };
+}
+
 export function interpolateNearest(index, width, height, X, Y, V) {
   const W = new V.constructor(width * height);
   const delaunay = Delaunay.from(