From c450fadcff609f5a555337b1e45d6f41e95eb04f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Philippe=20Rivi=C3=A8re?= Date: Sun, 18 Jun 2023 18:55:18 +0200 Subject: [PATCH 1/7] fix barycentric extrapolation closes #1700 --- src/marks/raster.js | 87 +++++++++++++------ test/output/rasterCa55Barycentric.svg | 2 +- test/output/rasterPenguinsBarycentric.svg | 2 +- .../rasterVaporEqualEarthBarycentric.svg | 2 +- test/output/rasterWalmartBarycentric.svg | 2 +- .../rasterWalmartBarycentricOpacity.svg | 2 +- 6 files changed, 64 insertions(+), 33 deletions(-) diff --git a/src/marks/raster.js b/src/marks/raster.js index fd601e3ced..90beb21a53 100644 --- a/src/marks/raster.js +++ b/src/marks/raster.js @@ -1,7 +1,7 @@ -import {blurImage, Delaunay, randomLcg, rgb} from "d3"; +import {blurImage, Delaunay, least, minIndex, 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,30 +282,14 @@ 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 {points, triangles, hull} = Delaunay.from( + index, + (i) => X[i], + (i) => Y[i] + ); const W = new V.constructor(width * height); + const I = new Uint8Array(width * height); const mix = mixer(V, random); for (let i = 0; i < triangles.length; i += 3) { const ta = triangles[i]; @@ -323,9 +307,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 +321,61 @@ 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); + I[i] = 1; } } } + + // Extrapolate by projection on the hull + const projections = Array.from(hull, (a, i) => { + const b = index[hull.at(i - 1)]; + a = index[a]; + return {a, b, project: segmentProject(X[b], Y[b], X[a], Y[a])}; + }); + const closest = (x, y) => + least( + projections.map(({a, b, project}) => ({a, b, p: project(x, y)})), + ({p: {distance2}}) => distance2 + ); + for (let i = 0; i < I.length; ++i) { + if (!I[i]) { + const x = i % width; + const y = Math.floor(i / width); + const { + a, + b, + p: {value} + } = closest(x + 0.5, y + 0.5); + W[i] = mix(V[a], 1 - value, V[b], value, V[a], 0, x, y); + } + } + return W; }; } +function segmentProject(x1, y1, x2, y2) { + const xm = (x1 + x2) / 2; + const ym = (y1 + y2) / 2; + const dx = x2 - xm; + const dy = y2 - ym; + return dx === 0 && dy === 0 + ? (x, y) => ({value: 0, distance2: (x - xm) ** 2 + (y - ym) ** 2}) + : (x, y) => { + const tx = x - xm; + const ty = y - ym; + const D = [(tx - 2 * dx) ** 2 + (ty - 2 * dy) ** 2, tx ** 2 + ty ** 2, (tx + 2 * dx) ** 2 + (ty + 2 * dy) ** 2]; + const c = minIndex(D); + const value = c === 0 ? 0 : c === 2 ? 1 : (D[0] - D[1]) / (D[0] + D[2] - 2 * D[1]); + return { + value, + distance2: (tx + (2 * value - 1) * dx) ** 2 + (ty + (2 * value - 1) * dy) ** 2 + }; + }; +} + 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/rasterPenguinsBarycentric.svg b/test/output/rasterPenguinsBarycentric.svg index fb5f3197e5..c316911345 100644 --- a/test/output/rasterPenguinsBarycentric.svg +++ b/test/output/rasterPenguinsBarycentric.svg @@ -66,7 +66,7 @@ body_mass_g → - + 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 @@ } - + From 5436702e01549ce5402d7130518054573d10b7c7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Philippe=20Rivi=C3=A8re?= Date: Sun, 18 Jun 2023 20:57:18 +0200 Subject: [PATCH 2/7] arithmetic --- src/marks/raster.js | 25 +++++++++---------------- 1 file changed, 9 insertions(+), 16 deletions(-) diff --git a/src/marks/raster.js b/src/marks/raster.js index 90beb21a53..3b38a00d89 100644 --- a/src/marks/raster.js +++ b/src/marks/raster.js @@ -1,4 +1,4 @@ -import {blurImage, Delaunay, least, minIndex, randomLcg, rgb} from "d3"; +import {blurImage, Delaunay, least, randomLcg, rgb} from "d3"; import {valueObject} from "../channel.js"; import {create} from "../context.js"; import {map, first, second, third, isTuples, isNumeric, isTemporal, identity} from "../options.js"; @@ -337,18 +337,14 @@ export function interpolatorBarycentric({random = randomLcg(42)} = {}) { const closest = (x, y) => least( projections.map(({a, b, project}) => ({a, b, p: project(x, y)})), - ({p: {distance2}}) => distance2 + (d) => d.p.dist2 ); for (let i = 0; i < I.length; ++i) { if (!I[i]) { const x = i % width; const y = Math.floor(i / width); - const { - a, - b, - p: {value} - } = closest(x + 0.5, y + 0.5); - W[i] = mix(V[a], 1 - value, V[b], value, V[a], 0, x, y); + const {a, b, p} = closest(x + 0.5, y + 0.5); + W[i] = mix(V[a], 1 - p.t, V[b], p.t, V[a], 0, x, y); } } @@ -362,17 +358,14 @@ function segmentProject(x1, y1, x2, y2) { const dx = x2 - xm; const dy = y2 - ym; return dx === 0 && dy === 0 - ? (x, y) => ({value: 0, distance2: (x - xm) ** 2 + (y - ym) ** 2}) + ? (x, y) => ({t: 0, dist2: (x - xm) ** 2 + (y - ym) ** 2}) : (x, y) => { const tx = x - xm; const ty = y - ym; - const D = [(tx - 2 * dx) ** 2 + (ty - 2 * dy) ** 2, tx ** 2 + ty ** 2, (tx + 2 * dx) ** 2 + (ty + 2 * dy) ** 2]; - const c = minIndex(D); - const value = c === 0 ? 0 : c === 2 ? 1 : (D[0] - D[1]) / (D[0] + D[2] - 2 * D[1]); - return { - value, - distance2: (tx + (2 * value - 1) * dx) ** 2 + (ty + (2 * value - 1) * dy) ** 2 - }; + const a = dx * (dx - tx) + dy * (dy - ty); + const b = dx * (dx + tx) + dy * (dy + ty); + const t = a > 0 && b > 0 ? a / (a + b) : +(a > b); + return {t, dist2: (tx + (2 * t - 1) * dx) ** 2 + (ty + (2 * t - 1) * dy) ** 2}; }; } From d013420bc960d2e675b643712bc6161eb35fa4ae Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Philippe=20Rivi=C3=A8re?= Date: Sun, 18 Jun 2023 21:41:07 +0200 Subject: [PATCH 3/7] more arithmetic --- src/marks/raster.js | 27 +++++++++++++-------------- 1 file changed, 13 insertions(+), 14 deletions(-) diff --git a/src/marks/raster.js b/src/marks/raster.js index 3b38a00d89..d8e22ab35e 100644 --- a/src/marks/raster.js +++ b/src/marks/raster.js @@ -353,20 +353,19 @@ export function interpolatorBarycentric({random = randomLcg(42)} = {}) { } function segmentProject(x1, y1, x2, y2) { - const xm = (x1 + x2) / 2; - const ym = (y1 + y2) / 2; - const dx = x2 - xm; - const dy = y2 - ym; - return dx === 0 && dy === 0 - ? (x, y) => ({t: 0, dist2: (x - xm) ** 2 + (y - ym) ** 2}) - : (x, y) => { - const tx = x - xm; - const ty = y - ym; - const a = dx * (dx - tx) + dy * (dy - ty); - const b = dx * (dx + tx) + dy * (dy + ty); - const t = a > 0 && b > 0 ? a / (a + b) : +(a > b); - return {t, dist2: (tx + (2 * t - 1) * dx) ** 2 + (ty + (2 * t - 1) * dy) ** 2}; - }; + const dx = x2 - x1; + const dy = y2 - y1; + if (dx === 0 && dy === 0) { + const xm = x1 + dx / 2; + const ym = y1 + dy / 2; + return (x, y) => ({t: 0, dist2: (x - xm) ** 2 + (y - ym) ** 2}); + } + return (x, y) => { + const a = dx * (x2 - x) + dy * (y2 - y); + const b = dx * (x - x1) + dy * (y - y1); + const t = a > 0 && b > 0 ? a / (a + b) : +(a > b); + return {t, dist2: (x - x2 + t * dx) ** 2 + (y - y2 + t * dy) ** 2}; + }; } export function interpolateNearest(index, width, height, X, Y, V) { From 9c33c894dfb388cf1bb4a9d6c11d3c6fd4a83703 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Philippe=20Rivi=C3=A8re?= Date: Sun, 18 Jun 2023 22:04:49 +0200 Subject: [PATCH 4/7] inline and simplify --- src/marks/raster.js | 45 ++++++++++++++++++--------------------------- 1 file changed, 18 insertions(+), 27 deletions(-) diff --git a/src/marks/raster.js b/src/marks/raster.js index d8e22ab35e..2fa4d03fd4 100644 --- a/src/marks/raster.js +++ b/src/marks/raster.js @@ -1,4 +1,4 @@ -import {blurImage, Delaunay, least, randomLcg, rgb} from "d3"; +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, identity} from "../options.js"; @@ -329,22 +329,20 @@ export function interpolatorBarycentric({random = randomLcg(42)} = {}) { } // Extrapolate by projection on the hull - const projections = Array.from(hull, (a, i) => { - const b = index[hull.at(i - 1)]; - a = index[a]; - return {a, b, project: segmentProject(X[b], Y[b], X[a], Y[a])}; - }); - const closest = (x, y) => - least( - projections.map(({a, b, project}) => ({a, b, p: project(x, y)})), - (d) => d.p.dist2 - ); for (let i = 0; i < I.length; ++i) { if (!I[i]) { - const x = i % width; - const y = Math.floor(i / width); - const {a, b, p} = closest(x + 0.5, y + 0.5); - W[i] = mix(V[a], 1 - p.t, V[b], p.t, V[a], 0, x, y); + const x0 = i % width; + const y0 = Math.floor(i / width); + const x = x0 + 0.5; + const y = y0 + 0.5; + let r = {dist2: Infinity}; + for (let j = 0; j < hull.length; ++j) { + const a = index[hull[j]]; + const b = index[hull.at(j - 1)]; + const {dist2, t} = segmentProject(X[b], Y[b], X[a], Y[a], x, y); + if (dist2 < r.dist2) r = {dist2, a, b, t}; + } + W[i] = mix(V[r.a], 1 - r.t, V[r.b], r.t, V[r.a], 0, x0, y0); } } @@ -352,20 +350,13 @@ export function interpolatorBarycentric({random = randomLcg(42)} = {}) { }; } -function segmentProject(x1, y1, x2, y2) { +function segmentProject(x1, y1, x2, y2, x, y) { const dx = x2 - x1; const dy = y2 - y1; - if (dx === 0 && dy === 0) { - const xm = x1 + dx / 2; - const ym = y1 + dy / 2; - return (x, y) => ({t: 0, dist2: (x - xm) ** 2 + (y - ym) ** 2}); - } - return (x, y) => { - const a = dx * (x2 - x) + dy * (y2 - y); - const b = dx * (x - x1) + dy * (y - y1); - const t = a > 0 && b > 0 ? a / (a + b) : +(a > b); - return {t, dist2: (x - x2 + t * dx) ** 2 + (y - y2 + t * dy) ** 2}; - }; + const a = dx * (x2 - x) + dy * (y2 - y); + const b = dx * (x - x1) + dy * (y - y1); + const t = a > 0 && b > 0 ? a / (a + b) : +(a > b); + return {t, dist2: (x - x2 + t * dx) ** 2 + (y - y2 + t * dy) ** 2}; } export function interpolateNearest(index, width, height, X, Y, V) { From 4434fbda9f5f8b3e9679e28f91a074f7615e0dea Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Philippe=20Rivi=C3=A8re?= Date: Mon, 19 Jun 2023 08:59:51 +0200 Subject: [PATCH 5/7] optimize --- src/marks/raster.js | 62 +++++++++++++------ test/output/rasterCa55Barycentric.svg | 2 +- test/output/rasterPenguinsBarycentric.svg | 2 +- .../rasterVaporEqualEarthBarycentric.svg | 2 +- test/output/rasterWalmartBarycentric.svg | 2 +- .../rasterWalmartBarycentricOpacity.svg | 2 +- 6 files changed, 48 insertions(+), 24 deletions(-) diff --git a/src/marks/raster.js b/src/marks/raster.js index 2fa4d03fd4..21e9598e81 100644 --- a/src/marks/raster.js +++ b/src/marks/raster.js @@ -328,35 +328,59 @@ export function interpolatorBarycentric({random = randomLcg(42)} = {}) { } } - // Extrapolate by projection on the hull - for (let i = 0; i < I.length; ++i) { - if (!I[i]) { - const x0 = i % width; - const y0 = Math.floor(i / width); - const x = x0 + 0.5; - const y = y0 + 0.5; - let r = {dist2: Infinity}; - for (let j = 0; j < hull.length; ++j) { - const a = index[hull[j]]; - const b = index[hull.at(j - 1)]; - const {dist2, t} = segmentProject(X[b], Y[b], X[a], Y[a], x, y); - if (dist2 < r.dist2) r = {dist2, a, b, t}; - } - W[i] = mix(V[r.a], 1 - r.t, V[r.b], r.t, V[r.a], 0, x0, y0); - } - } + extrapolateBarycentric(W, I, X, Y, V, width, height, hull, index, mix); return W; }; } +// Extrapolate by finding the closest point on the hull. Optimized with a +// Delaunay search on the interpolated hull. +function extrapolateBarycentric(W, I, 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 points = []; + const refs = []; + for (let j = 0; j < hull.length; ++j) { + const xa = X.at(j - 1); + const ya = Y.at(j - 1); + const dx = X.at(j) - xa; + const dy = Y.at(j) - ya; + const s = 1 / (1 + (Math.hypot(dx, dy) << 1)); + for (let d = 0; d < 1; d += s) { + points.push(xa + d * dx, ya + d * dy); + refs.push(j); + } + } + const delaunay = new Delaunay(points); + let iy, ix; + for (let y = 0; y < height; ++y) { + const yp = y + 0.5; + ix = iy; + for (let x = 0; x < width; ++x) { + const i = x + width * y; + const xp = x + 0.5; + if (!I[i]) { + ix = delaunay.find(xp, yp, ix); + if (x === 0) iy = ix; + const j = refs[ix]; + 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); + } + } + } +} + +// 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); - const t = a > 0 && b > 0 ? a / (a + b) : +(a > b); - return {t, dist2: (x - x2 + t * dx) ** 2 + (y - y2 + t * dy) ** 2}; + return a > 0 && b > 0 ? a / (a + b) : +(a > b); } export function interpolateNearest(index, width, height, X, Y, V) { diff --git a/test/output/rasterCa55Barycentric.svg b/test/output/rasterCa55Barycentric.svg index a252b93c47..ae8c4ed3f2 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/rasterPenguinsBarycentric.svg b/test/output/rasterPenguinsBarycentric.svg index c316911345..1d3a0de8ee 100644 --- a/test/output/rasterPenguinsBarycentric.svg +++ b/test/output/rasterPenguinsBarycentric.svg @@ -66,7 +66,7 @@ body_mass_g → - + diff --git a/test/output/rasterVaporEqualEarthBarycentric.svg b/test/output/rasterVaporEqualEarthBarycentric.svg index 68d95c0ef5..8fb9001ab9 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 7cf25a853b..44c22fb2ff 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 a4517c0413..f542d1b370 100644 --- a/test/output/rasterWalmartBarycentricOpacity.svg +++ b/test/output/rasterWalmartBarycentricOpacity.svg @@ -14,7 +14,7 @@ } - + From d72e63203873b6befe72ea1c4cafe60e46fee2c2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Philippe=20Rivi=C3=A8re?= Date: Fri, 23 Jun 2023 15:29:53 +0200 Subject: [PATCH 6/7] fix barycentric extrapolation variant (#1705) * An exact algorithm, which doesn't make use of a delaunay search * repeat const * fix a floating point precision issue 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. * test barycentric data binding --------- Co-authored-by: Mike Bostock --- src/marks/raster.js | 86 +++++++++++------- test/output/rasterCa55Barycentric.svg | 2 +- test/output/rasterFacet.svg | 88 +++++++++++++++++++ test/output/rasterPrecision.svg | 59 +++++++++++++ .../rasterVaporEqualEarthBarycentric.svg | 2 +- test/output/rasterWalmartBarycentric.svg | 2 +- .../rasterWalmartBarycentricOpacity.svg | 2 +- test/plots/index.ts | 1 + test/plots/raster-precision.ts | 40 +++++++++ 9 files changed, 247 insertions(+), 35 deletions(-) create mode 100644 test/output/rasterFacet.svg create mode 100644 test/output/rasterPrecision.svg create mode 100644 test/plots/raster-precision.ts diff --git a/src/marks/raster.js b/src/marks/raster.js index 21e9598e81..504b4a8972 100644 --- a/src/marks/raster.js +++ b/src/marks/raster.js @@ -288,9 +288,10 @@ export function interpolatorBarycentric({random = randomLcg(42)} = {}) { (i) => X[i], (i) => Y[i] ); - const W = new V.constructor(width * height); - const I = new Uint8Array(width * height); + 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,51 +324,38 @@ export function interpolatorBarycentric({random = randomLcg(42)} = {}) { if (gc < 0) continue; const i = x + width * y; W[i] = mix(va, ga, vb, gb, vc, gc, x, y); - I[i] = 1; + S[i] = 1; } } } - - extrapolateBarycentric(W, I, X, Y, V, width, height, hull, index, mix); - + extrapolateBarycentric(W, S, X, Y, V, width, height, hull, index, mix); return W; }; } -// Extrapolate by finding the closest point on the hull. Optimized with a -// Delaunay search on the interpolated hull. +// Extrapolate by finding the closest point on the hull. function extrapolateBarycentric(W, I, 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 points = []; - const refs = []; - for (let j = 0; j < hull.length; ++j) { - const xa = X.at(j - 1); - const ya = Y.at(j - 1); - const dx = X.at(j) - xa; - const dy = Y.at(j) - ya; - const s = 1 / (1 + (Math.hypot(dx, dy) << 1)); - for (let d = 0; d < 1; d += s) { - points.push(xa + d * dx, ya + d * dy); - refs.push(j); - } - } - const delaunay = new Delaunay(points); - let iy, ix; + 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; - ix = iy; for (let x = 0; x < width; ++x) { const i = x + width * y; - const xp = x + 0.5; if (!I[i]) { - ix = delaunay.find(xp, yp, ix); - if (x === 0) iy = ix; - const j = refs[ix]; - 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); + 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; + } + } } } } @@ -383,6 +371,42 @@ function segmentProject(x1, y1, x2, y2, x, y) { 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 ae8c4ed3f2..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 @@ + + + + + 0 + + + 1 + + + + + + + + + + + + −1 + 0 + 1 + + + + + + + + + + + + + 0 + + + 0 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + \ No newline at end of file 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 @@ + + + + + + + + + + + + + 350 + 300 + 250 + 200 + 150 + 100 + 50 + + + + + + + + + + + 100 + 200 + 300 + 400 + 500 + 600 + + + + + + + + + + + \ No newline at end of file diff --git a/test/output/rasterVaporEqualEarthBarycentric.svg b/test/output/rasterVaporEqualEarthBarycentric.svg index 8fb9001ab9..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 44c22fb2ff..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 f542d1b370..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 93bc831398..b15aa96a04 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"}) + ] + }); +} From 0f5a9b4e6e723547caa4af076eb5f18fd6d37b2c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Philippe=20Rivi=C3=A8re?= Date: Fri, 23 Jun 2023 23:24:14 +0200 Subject: [PATCH 7/7] S --- src/marks/raster.js | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/marks/raster.js b/src/marks/raster.js index 504b4a8972..61abbc2eea 100644 --- a/src/marks/raster.js +++ b/src/marks/raster.js @@ -334,7 +334,7 @@ export function interpolatorBarycentric({random = randomLcg(42)} = {}) { } // Extrapolate by finding the closest point on the hull. -function extrapolateBarycentric(W, I, X, Y, V, width, height, hull, index, mix) { +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]]); @@ -345,7 +345,7 @@ function extrapolateBarycentric(W, I, X, Y, V, width, height, hull, index, mix) const yp = y + 0.5; for (let x = 0; x < width; ++x) { const i = x + width * y; - if (!I[i]) { + 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;