30:T106c,chart = { const context = DOM.context2d(width, height); const $ELEMS = [context.canvas]; //------------------------------------- const filteredGradient = gradient.filter((d, i) => !isNaN(P[i])); const m = 7 / Math.sqrt(d3.max(filteredGradient, ([x, y]) => x * x + y * y)); //------------------------------------- if (display.includes("potential")) { context.drawImage(viewof image1, 0, 0); } //------------------------------------- if (display.includes("charges")) { const chargeColor = d3.scaleDiverging( [d3.min(C), 0, d3.max(C)], ["red", "transparent", "steelblue"] ); for (let i = 0; i < points.length; i++) { if (isNaN(C[i])) continue; context.beginPath(); voronoi.renderCell(i, context); context.fillStyle = chargeColor(C[i]); context.fill(); } } //------------------------------------- if (display.includes("streaklines")) { const line = d3.line().context(context).curve(d3.curveCardinal); for (const l of streaklines()) { context.lineWidth = 0.5; context.strokeStyle = "black"; context.beginPath(); line(l); context.stroke(); } } //------------------------------------- if (display.includes("isolines")) { const path = d3.geoPath().context(context); for (const l of isolines) { context.lineWidth = 0.5; context.strokeStyle = "black"; context.beginPath(); path(l); context.stroke(); } } //------------------------------------- if (display.includes("gradient")) { context.beginPath(); for (let i = 0; i < points.length; i++) { if (!gradient[i] || isNaN(P[i])) continue; const p = points[i]; const d = gradient[i].map((d) => d * m); const l = Math.hypot(...d); if (l > 0.15) { context.moveTo(p[0] - d[0], p[1] - d[1]); context.lineTo(p[0] + d[0], p[1] + d[1]); if (l > 2) { context.moveTo( p[0] + d[0] - (3 * d[0] + 1.5 * d[1]) / l, p[1] + d[1] - (3 * d[1] - 1.5 * d[0]) / l ); context.lineTo(p[0] + d[0], p[1] + d[1]); context.lineTo( p[0] + d[0] - (3 * d[0] - 1.5 * d[1]) / l, p[1] + d[1] - (3 * d[1] + 1.5 * d[0]) / l ); } } } context.strokeStyle = "black"; context.lineWidth = 0.5; context.stroke(); } //------------------------------------- if (display.includes("centers")) { const r = 1.3; for (let i = 0; i < points.length; i++) { const p = points[i]; context.fillStyle = isNaN(P[i]) ? "black" : "#777"; context.beginPath(); context.moveTo(p[0] + r, p[1]); context.arc(p[0], p[1], r, 0, 6.25); context.fill(); } } //------------------------------------- if (display.includes("particles")) { const f = 2 / d3.max(gradient, (d) => (d ? Math.hypot(...d) : NaN)); const context = DOM.context2d(width, height); const simulation = d3 .forceSimulation(nodes) .alphaDecay(0) .force("wind", () => { nodes.forEach((d) => { const i = voronoi.delaunay.find(d.x, d.y, d.last); if (i < 0 || Math.random() < 0.01 || !gradient[i] || isNaN(P[i])) { d.x = width * Math.random(); d.y = height * Math.random(); d.vx = 0; d.vy = 0; d.last = undefined; } else { const vx = gradient[i][0] * f; const vy = gradient[i][1] * f; d.vx = vx + 0.0 * (Math.random() - 0.5); d.vy = vy + 0.0 * (Math.random() - 0.5); d.last = i; } }); }) .on("tick", () => { fade(context, 0.1); context.beginPath(); context.strokeStyle = "black"; context.lineWidth = 0.7; for (const { x, y, vx, vy } of nodes) { context.moveTo(x, y); context.lineTo(x + vx, y + vy); } context.stroke(); }); invalidation.then(() => simulation.stop()); context.canvas.style.position = "absolute"; context.canvas.style.top = "0px"; $ELEMS.unshift(context.canvas); } return html`

${$ELEMS}`; }31:T5f2,The Poisson potential is the solution of [Poisson’s equation](https://en.wikipedia.org/wiki/Poisson%27s_equation), a partial differential equation usually written as: > ${tex`\nabla^2 P = -c`} where *c* is the charge. In whole words, this could be read as: > “the charge *c(x)* in location *x* is the difference between the potential *P(x)* and its mean in the neighborhood of *x*” Another way to understand this equation is to say that, if particles appear in *x* at a rate of *c(x)* and diffuse uniformly, the resulting flow follows the gradient of *P*. This equation can be solved analytically in some cases (universal gravitation!), but for the general case, we have to resort to finite element methods, where space is discretized by, say, a rectangular grid, and finite differences between neighbors on the grid replace the derivatives. This hints to a graph version of the Poisson potential, which can be described as: > “the charge *c(x)* in location *x* is the difference between *P(x)* and its mean across the neighbors of *x*” The Poisson equation does not care if we add to P a constant value —or even a “flat” function that reflects an external potential coming from outside sources. But add a few hypotheses, and the solution is unique: in this case, we consider that the graph has no exchanges with its exterior (which is equivalent to [Neumann’s boundary condition](https://en.wikipedia.org/wiki/Neumann_boundary_condition)), and that its value averages to zero over the whole domain.32:T642,function potential( graph, charges, { alpha = 0.3, precision = 1e-8, maxSteps = 20000 } = {} ) { const n = charges.length; const D = new Float64Array(n); // differences const A = new Float64Array(n); const q = d3.mean(charges); // normalize the flow by the count of valid neighbors const nneighbors = new Uint16Array(n); for (let s = 0; s < graph.sources.length; s++) { const i = graph.sources[s]; const j = graph.targets[s]; if (!isNaN(charges[j])) nneighbors[i]++; if (!isNaN(charges[i])) nneighbors[j]++; } // impedance normalization const conductance = graph.costs ? Float32Array.from(graph.costs, (d) => 1 / d) : null; const Q = charges.map( (d, i) => alpha * (d - q) * (!nneighbors[i] ? 0 // special case isolated points : 1 / nneighbors[i]) ); const E = precision * A.length; let steps = 0, error, c; const { sources, targets } = graph; do { // add charges for (let i = 0; i < D.length; i++) D[i] = Q[i]; // diffuse for (let s = 0; s < sources.length; s++) { const i = sources[s]; const j = targets[s]; // let d = alpha * (A[j] - A[i]); // Jacobi let d = alpha * (D[j] + A[j] - D[i] - A[i]); // SOR if (conductance) d *= conductance[s]; if (!isNaN(d)) { D[i] += d; D[j] -= d; } } for (let i = 0; i < D.length; i++) A[i] += D[i]; // measure change error = d3.sum(D, (d) => d * d); steps++; } while (error > E && steps < maxSteps); return Object.assign(A, { steps, error, E }); }33:T4aa,**API**: **potential**(*graph*, *charges*, *options*) - **graph** is an object with an array of sources […i] and an array of targets […j] denoting an connection between node *i* and node *j*; an impedance can be specified as an array of costs —allowing to specify “more difficult” connections on the graph, with a higher impedance (for example, the impedance can be relative to the arc’s length). This is not normalized and might result in a failure of convergence.. - **charges** is an array of floats denoting the charge on each node *i*; nodes for which the charge is NaN will be ignored. - *options*.**alpha** controls the amount of movement at each step; the lower the more precise the result—defaults to 0.3. - The algorithm stops whenever the computation step results in a mean square error less than *options*.**precision** (which default to 1e-6), or when the number of steps exceeds *options*.**maxSteps**, which defaults to 20,000. *Note:* the graph’s API design is identical to the one I use in [dijkstra](/@fil/dijkstra), but here each link is processed symmetrically, so there is no need to duplicate the arcs to get a symmetric graph. I might revisit this later.34:T67e,viewof image0 = { const V = new Float32Array(width * height).fill(NaN); // fill all triangles, interpolate via barycentric coordinates const delaunay = d3.Delaunay.from(points); for (let i = 0; i < delaunay.triangles.length; i += 3) { const [a, b, c] = delaunay.triangles.slice(i, i + 3); const [pa, pb, pc] = Array.from([a, b, c], (i) => P[i]); if (isNaN(pa) || isNaN(pb) || isNaN(pc)) continue; const [Ax, Bx, Cx] = Array.from([a, b, c], (i) => points[i][0]); const [Ay, By, Cy] = Array.from([a, b, c], (i) => points[i][1]); const [x0, x1] = d3.extent([Ax, Bx, Cx]); const [y0, y1] = d3.extent([Ay, By, Cy]); const z = (By - Cy) * (Ax - Cx) + (Ay - Cy) * (Cx - Bx); if (!z) continue; for (let x = Math.floor(x0); x < x1; x++) { for (let y = Math.floor(y0); y < y1; y++) { if (x < 0 || x >= width || y < 0 || y >= height) continue; const ga = ((By - Cy) * (x - Cx) + (y - Cy) * (Cx - Bx)) / z; const gb = ((Cy - Ay) * (x - Cx) + (y - Cy) * (Ax - Cx)) / z; const gc = 1 - ga - gb; if (ga >= 0 && gb >= 0 && gc >= 0) V[x + width * y] = ga * pa + gb * pb + gc * pc; } } } // fill the missing parts with the value of the closest dot const values = points .map((d, i) => [...d, P[i]]) .filter(([, , p]) => !isNaN(p)); const quad = d3.quadtree().addAll(values); for (let x = 0; x < width; x++) { for (let y = 0; y < height; y++) { if (isNaN(V[x + width * y])) V[x + width * y] = quad.find(x, y)[2]; } } return heatmap(V, { width, color: sourcesOnly ? d3.interpolateInferno : d3.interpolateRdYlBu }); }35:T47e,streaklines = () => { const visited = new Uint8Array(points.length); let maxstep = 200000; const streaklines = []; for (const i of d3.shuffle(d3.range(points.length))) { // already visited? drop if (visited[i]) continue; const streakline = [[...points[i], 0]]; visited[i] = 1; // logically useless but clearer let delta; let x; let y; let j = i; let vx, vy, dx, dy; for (const direction of [-1, 1]) { x = points[i][0]; y = points[i][1]; vx = 0; vy = 0; do { dx = vx; dy = vy; j = voronoi.delaunay.find(x, y, j); if (isNaN(P[j])) break; visited[j] = 1; [vx, vy] = grad([x, y], image1).map((d) => d * direction); const delta = Math.hypot(vx, vy); x += vx / delta; y += vy / delta; streakline.push([x, y]); } while ( vx * dx + vy * dy >= 0 && //delta > 1e-12 && x > 0 && y > 0 && x < width && y < height && maxstep-- > 0 ); streakline.reverse(); } streaklines.push(streakline); } return streaklines; } 36:[] 2e:{"title":"Poisson potential","cells":[{"id":530,"value":"# Poisson potential\n\n*For a cartographic application, see [Potential dollar flows](/@visionscarto/potential-dollar-flows).*","pinned":false,"hidden":false,"mode":"md"},{"id":930,"value":"$30","pinned":false,"hidden":false,"mode":"ojs"},{"id":961,"value":"viewof display = Inputs.checkbox(\n [\n \"potential\",\n \"charges\",\n \"gradient\",\n \"isolines\",\n \"streaklines\",\n \"centers\",\n \"particles\"\n ],\n {\n label: \"display\",\n value: [\"potential\"]\n }\n)","pinned":false,"hidden":false,"mode":"ojs"},{"id":596,"value":"viewof sources = Inputs.range([1, 20], {\n label: \"sources and sinks\",\n step: 1,\n value: 8\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":1417,"value":"viewof blurRadius = Inputs.range([0, 50], {\n label: \"blur radius\",\n step: 0.1,\n value: 5\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":602,"value":"viewof sourcesOnly = Inputs.toggle({\n label: \"sources only\"\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":649,"value":"viewof distributionType = Inputs.select(\n [\n \"poisson\",\n \"hex\",\n \"grid\", // note: the grid distribution is not well-served by the delaunay triangulation\n \"random\",\n \"normal\"\n ],\n { label: \"point distribution\" }\n)","pinned":false,"hidden":false,"mode":"ojs"},{"id":611,"value":"viewof npoints = Inputs.range([1000, 20000], {\n step: 10,\n label: \"number of points\",\n value: 5000\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":0,"value":"$31","pinned":false,"hidden":false,"mode":"md"},{"id":716,"value":"The following “relaxation” algorithm converges to the solution, by diffusing the charges through the graph:","pinned":false,"hidden":false,"mode":"md"},{"id":18,"value":"$32","pinned":true,"hidden":false,"mode":"ojs"},{"id":1790,"value":"It is far from being a fast algorithm (if you have a rectangular raster image, doing a Fourier transform is tremendously faster), but it's enough for now. (My suspicion is that a way of making it faster would be to find good spectral cuts in the graph, and diffuse the value first on large aggregates of points, then refine the aggregates. Which is what the Fourier transform does, in a sense, when we work on a regular grid. More on this later…)","pinned":false,"hidden":false,"mode":"md"},{"id":675,"value":"*Result, obtained in ${P.steps} steps:*","pinned":false,"hidden":false,"mode":"md"},{"id":195,"value":"P = potential(graph, C)","pinned":true,"hidden":false,"mode":"ojs"},{"id":806,"value":"$33","pinned":false,"hidden":false,"mode":"md"},{"id":680,"value":"---\n\n*Thanks to all the attendants of the [Transcarto](https://transcarto.sciencesconf.org/) workshop for the discussions on* transformative cartography; *in particular, see Étienne Côme’s [winds](https://www.comeetie.fr/galerie/wind/) project for inspiration (and its associated [observable notebook](https://observablehq.com/@comeetie/ttt-poisson-potential)).*\n\n### Annex\n\nThe following defines a random graph; if you already have a graph, you can pass it directly to the potential function above—and there’s no need to use that part of the code:","pinned":false,"hidden":false,"mode":"md"},{"id":638,"value":"points = [...pick2d(width, height, npoints, distributionType)]","pinned":true,"hidden":false,"mode":"ojs"},{"id":313,"value":"C = {\n const C = new Float32Array(points.length);\n\n for (const i of pointIndex.slice(0, sources)) {\n const v = 1 + (!sourcesOnly && i % 2 ? -3 : 1);\n C[i] += v;\n graph.sources.forEach((k, s) => {\n let j = graph.targets[s];\n if (j === i) (j = k), (k = i);\n if (k === i) C[j] += v * 0.6;\n });\n }\n\n return C;\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":908,"value":"C.filter((d) => d)","pinned":true,"hidden":false,"mode":"ojs"},{"id":659,"value":"pointIndex = d3.shuffle(d3.range(points.length))","pinned":true,"hidden":false,"mode":"ojs"},{"id":731,"value":"// the voronoi partition defines the graph *and* is useful to draw the chart\nvoronoi = d3.Delaunay.from(points).voronoi([0, 0, width, height])","pinned":true,"hidden":false,"mode":"ojs"},{"id":170,"value":"graph = {\n let links = d3\n .range(points.length)\n .flatMap((i) => Array.from(voronoi.neighbors(i), (j) => [i, j]));\n\n // on the square grid, remove diagonals\n if (distributionType === \"grid\") {\n let stride;\n points.find((d, i) => d[0] > 0 && (stride = i));\n links = links.filter(([i, j]) => [1, stride].includes(Math.abs(i - j)));\n }\n\n return {\n sources: Uint16Array.from(links, ([i]) => i),\n targets: Uint16Array.from(links, ([, j]) => j)\n // costs (defaults to ...1)\n };\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":636,"value":"import { pick2d } from \"@fil/2d-point-distributions\"","pinned":true,"hidden":false,"mode":"ojs"},{"id":1100,"value":"To draw the particles system, we need a fade function, and an array of positions","pinned":false,"hidden":false,"mode":"md"},{"id":1042,"value":"function fade(context, alpha) {\n const { globalAlpha, globalCompositeOperation } = context;\n context.globalCompositeOperation = \"destination-out\";\n context.fillStyle = \"white\";\n context.globalAlpha = alpha || 1;\n context.fillRect(0, 0, context.canvas.width, context.canvas.height);\n context.globalAlpha = globalAlpha;\n context.globalCompositeOperation = globalCompositeOperation;\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":1104,"value":"nodes = {\n const random = d3.randomNormal(0, 1000);\n return Array.from({ length: Math.min(width * 10, npoints * 2) }, () => ({\n x: width / 2 + random(),\n y: height / 2 + random()\n }));\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":1363,"value":"An image is computed for the whole map, with barycentric interpolation","pinned":false,"hidden":false,"mode":"md"},{"id":1570,"value":"$34","pinned":true,"hidden":false,"mode":"ojs"},{"id":1532,"value":"import { heatmap } from \"@fil/heatmap\"","pinned":true,"hidden":false,"mode":"ojs"},{"id":1374,"value":"We can blur it to smooth out the jaggedness of the graph","pinned":false,"hidden":false,"mode":"md"},{"id":1379,"value":"B = require(\"array-blur\")","pinned":true,"hidden":false,"mode":"ojs"},{"id":1381,"value":"viewof image1 = heatmap(B.blur().radius(blurRadius).width(width)(image0), {\n width,\n color: sourcesOnly ? d3.interpolateInferno : d3.interpolateRdYlBu\n})","pinned":true,"hidden":false,"mode":"ojs"},{"id":1396,"value":"contours = d3.contours().size([width, height])(image1)","pinned":true,"hidden":false,"mode":"ojs"},{"id":1491,"value":"// filtered contours\nisolines = contours.map(({ value, coordinates }) => ({\n type: \"MultiLineString\",\n value,\n coordinates: coordinates.flatMap((polygon) =>\n polygon.flatMap((ring) => {\n const lines = [];\n let line = [];\n let i;\n for (const p of ring) {\n i = voronoi.delaunay.find(...p, i);\n if (isNaN(P[i])) {\n if (line.length > 1) {\n lines.push(line);\n line = [];\n }\n } else line.push(p);\n }\n if (line.length > 1) lines.push(line);\n return lines;\n })\n )\n}))","pinned":true,"hidden":false,"mode":"ojs"},{"id":1408,"value":"function grad(p, image) {\n let x = Math.floor(p[0]),\n y = Math.floor(p[1]);\n if (x < 0) x = 0;\n if (y < 0) y = 0;\n if (x === width - 1) x--;\n if (y === height - 1) y--;\n return [\n image[x + width * y] - image[x + 1 + width * y],\n image[x + width * y] - image[x + width * (y + 1)]\n ];\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":1384,"value":"gradient = points.map((d, i) => (isNaN(P[i]) ? null : grad(d, image1)))","pinned":true,"hidden":false,"mode":"ojs"},{"id":1305,"value":"The streaklines (or streamlines) method is adapted from a heuristic developed by Whittaker in 1977 ([PDF](https://www.semanticscholar.org/paper/Automated-Streamline-Analysis-Whittaker/768839999d6156f555cf1978a8ca1b4963fdcc86)). *Note to self: there is a better version of this function in the [Pantelleria](https://observablehq.com/@visionscarto/deriving-pantelleria) notebook; however it works only on a grid.*","pinned":false,"hidden":false,"mode":"md"},{"id":1401,"value":"$35","pinned":true,"hidden":false,"mode":"ojs"},{"id":7,"value":"height = 600","pinned":true,"hidden":false,"mode":"ojs"}],"files":"$W36","stdlib":"1","titling":"auto"}