diff --git a/src/marks/raster.js b/src/marks/raster.js
index fd601e3ced..61abbc2eea 100644
--- a/src/marks/raster.js
+++ b/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(
diff --git a/test/output/rasterCa55Barycentric.svg b/test/output/rasterCa55Barycentric.svg
index 463c98d850..a252b93c47 100644
--- a/test/output/rasterCa55Barycentric.svg
+++ b/test/output/rasterCa55Barycentric.svg
@@ -14,6 +14,6 @@
}
-
+
\ No newline at end of file
diff --git a/test/output/rasterFacet.svg b/test/output/rasterFacet.svg
new file mode 100644
index 0000000000..a339e632b3
--- /dev/null
+++ b/test/output/rasterFacet.svg
@@ -0,0 +1,88 @@
+
\ No newline at end of file
diff --git a/test/output/rasterPenguinsBarycentric.svg b/test/output/rasterPenguinsBarycentric.svg
index fb5f3197e5..1d3a0de8ee 100644
--- a/test/output/rasterPenguinsBarycentric.svg
+++ b/test/output/rasterPenguinsBarycentric.svg
@@ -66,7 +66,7 @@
body_mass_g →
-
+
diff --git a/test/output/rasterPrecision.svg b/test/output/rasterPrecision.svg
new file mode 100644
index 0000000000..58317959bb
--- /dev/null
+++ b/test/output/rasterPrecision.svg
@@ -0,0 +1,59 @@
+
\ No newline at end of file
diff --git a/test/output/rasterVaporEqualEarthBarycentric.svg b/test/output/rasterVaporEqualEarthBarycentric.svg
index 97351c8168..68d95c0ef5 100644
--- a/test/output/rasterVaporEqualEarthBarycentric.svg
+++ b/test/output/rasterVaporEqualEarthBarycentric.svg
@@ -17,7 +17,7 @@
-
+
diff --git a/test/output/rasterWalmartBarycentric.svg b/test/output/rasterWalmartBarycentric.svg
index e5d72b6d2a..7cf25a853b 100644
--- a/test/output/rasterWalmartBarycentric.svg
+++ b/test/output/rasterWalmartBarycentric.svg
@@ -14,7 +14,7 @@
}
-
+
diff --git a/test/output/rasterWalmartBarycentricOpacity.svg b/test/output/rasterWalmartBarycentricOpacity.svg
index 02629ee3da..a4517c0413 100644
--- a/test/output/rasterWalmartBarycentricOpacity.svg
+++ b/test/output/rasterWalmartBarycentricOpacity.svg
@@ -14,7 +14,7 @@
}
-
+
diff --git a/test/plots/index.ts b/test/plots/index.ts
index b0544b26aa..fcfc9f075b 100644
--- a/test/plots/index.ts
+++ b/test/plots/index.ts
@@ -232,6 +232,7 @@ export * from "./random-quantile.js";
export * from "./random-walk.js";
export * from "./raster-ca55.js";
export * from "./raster-penguins.js";
+export * from "./raster-precision.js";
export * from "./raster-vapor.js";
export * from "./raster-walmart.js";
export * from "./rect-band.js";
diff --git a/test/plots/raster-precision.ts b/test/plots/raster-precision.ts
new file mode 100644
index 0000000000..b0ae134dd2
--- /dev/null
+++ b/test/plots/raster-precision.ts
@@ -0,0 +1,40 @@
+import * as Plot from "@observablehq/plot";
+import * as d3 from "d3";
+
+// Test for floating point precision issue in interpolateBarycentric.
+export async function rasterPrecision() {
+ const data = d3.range(4).map((i) => {
+ const x = i % 2;
+ const y = Math.floor(i / 2);
+ return [49.4 + 100 * (x + y), 150.4 + 100 * (x - y)];
+ });
+ return Plot.plot({
+ x: {type: "identity"},
+ y: {type: "identity"},
+ color: {scheme: "Sinebow"},
+ marks: [
+ Plot.raster(data, {
+ fill: (d, i) => i,
+ interpolate: "barycentric"
+ }),
+ Plot.dot(data, {fill: (d, i) => i, stroke: "white"})
+ ]
+ });
+}
+
+export async function rasterFacet() {
+ const points = d3.range(0, 2 * Math.PI, Math.PI / 10).map((d) => [Math.cos(d), Math.sin(d)]);
+ return Plot.plot({
+ aspectRatio: 1,
+ inset: 100,
+ color: {scheme: "Sinebow"},
+ marks: [
+ Plot.raster(points, {
+ fill: "0",
+ fx: (d, i) => i % 2,
+ interpolate: "barycentric"
+ }),
+ Plot.dot(points, {fx: (d, i) => i % 2, fill: "0", stroke: "white"})
+ ]
+ });
+}