28:null 29:["$","title",null,{"children":"Disk Sampling on a Sphere | Jacob Rus | Observable"}] 2e:T752,{ const background_color = "#d1f2ff"; const circle_color = 'rgba(250, 10, 60, .45)'; const grid_color = "#9ac"; const land_color = 'rgba(40, 70, 0, .25)'; const sphere = { type: "Sphere" }; const w = Math.min(width, 550) const context = DOM.context2d(w, w); const projection = d3.geoOrthographic() .rotate([10, -60, 0]) .scale(250) .center([0, 90]) .fitExtent([[10, 10], [w - 10, w - 10]], sphere) .precision(0.05);; context.canvas.classList.add('draggable'); const path = d3.geoPath(projection, context); let land = land50; let current_circles = [] function render(grid=true) { context.fillStyle = background_color; context.beginPath(), path(sphere), context.fill(); context.strokeStyle = grid_color; context.lineWidth = .3 + .07 * radius; if (grid) context.beginPath(), path(lonlat_hexgrid), context.stroke(); context.fillStyle = land_color; context.beginPath(), path(land), context.fill(); context.lineWidth = 1.0; context.strokeStyle = "#999"; context.beginPath(), path(sphere), context.stroke(); context.fillStyle = circle_color; context.beginPath(), current_circles.forEach(path), context.fill(); } let i = 0; let k = circles.length / 200; function tick() { context.fillStyle = circle_color; context.beginPath(); for (let j = 0; j < k; j++) { const circle = circles[i++]; path(circle); current_circles[i] = circle; } context.fill(); if (i < circles.length) requestAnimationFrame(tick); } render(); tick(); return d3.select(context.canvas) .call(drag(projection) .on("start.render", () => {land = land110; render(radius > 4);}) .on("drag.render", () => render(radius > 4)) .on("end.render", () => {land = land50; render();})) .call(() => render()) .node(); }2f:T5f0,{ const background_color = "#d1f2ff"; const circle_color = 'rgba(250, 10, 60, .45)'; const grid_color = "#9ac"; const land_color = 'rgba(40, 70, 0, .25)'; const w = Math.min(width, 640) const sphere = { type: "Sphere" }; const projection = quincuncial_projection .rotate([0, 0, 0]) .scale(250) .center([0, 90]) .fitExtent([[0, 0], [w, w]], sphere) .precision(0.05); const context = DOM.context2d(w, w); const path = d3.geoPath(projection, context); let land = land50; let current_circles = []; function render(x, y) { context.fillStyle = background_color; context.fillRect(0, 0, w, w); context.strokeStyle = grid_color; context.lineWidth = .3 + .06 * radius; context.beginPath(), path(lonlat_hexgrid), context.stroke() context.fillStyle = land_color; context.beginPath(), path(land), context.fill(); context.lineWidth = 1.0; context.strokeStyle = "#999"; context.beginPath(), path(sphere), context.stroke(); context.fillStyle = circle_color; context.beginPath(), current_circles.forEach(path), context.fill(); } let i = 0; let k = circles.length / 200; function tick() { context.fillStyle = circle_color; context.beginPath(); for (let j = 0; j < k; j++) { const circle = circles[i++]; path(circle); current_circles[i] = circle; } context.fill(); if (i < circles.length) requestAnimationFrame(tick); } render(); tick(); return context.canvas; }30:T135c,{ var w0 = (width > 400) ? Math.min(width, 600)/2 - 20 : width - 40; var images = html`

`; var w1 = Math.min(width - 20, 300); var hexagon_coverage = html`

`; return md` ### Differences from planar disk sampling The main idea of the algorithm is the same between the planar and spherical cases. But some of the details change to fit the context (and the resulting implementation here is about twice as much code). * Instead of ${tex`\langle x, y \rangle`} coordinates in the rectangle ${tex`[0, w] \times [0, h]`} representing planar points, our ${tex`\langle x, y \rangle`} can lie anywhere in the plane ${tex`\mathbb{R} \times \mathbb{R}`}, representing the [stereographic projection of the sphere](/@jrus/planisphere). Starting from geographic coordinates the stereographic projection is ${tex`\quad \langle \lambda, \phi \rangle \mapsto \left\langle \tan \tfrac12 \bigl(\tfrac12\pi - \phi\bigr)\;\! \cos \lambda,\, \tan \tfrac12 \bigl(\tfrac12\pi - \phi\bigr)\;\! \sin \lambda \right\rangle,`} with inverse ${tex`\quad \langle x, y \rangle \mapsto \left\langle \operatorname{atan2}(y, x) ,\, \tfrac12\pi - 2\arctan\!\sqrt{x^2 + y^2} \right\rangle.`} * Instead of using a grid of squares covering the rectangle as a spatial-index data structure, we instead project the sphere into the square ${tex`[-1, 1] \times [-1, 1]`} using [spherical area coordinates](/@jrus/sac-quincuncial). In the northern hemisphere ${tex`x^2 + y^2 < 1,`} this map projection is ${tex`\quad \langle x, y \rangle \mapsto \left\langle \dfrac4\pi \arctan \dfrac{x}{1 + |y|},\, \dfrac4\pi \arctan \dfrac{y}{1 + |x|} \right\rangle.`} Then we use a grid of hexagons. It requires a bit of tricky arithmetic to round a pair of floating point coordinates to the index of its enclosing hexagon: a variant of [hexagon grid rounding](/@jrus/hexround) whose behavior depends on the quadrant. ${images} * Because the sphere has no boundary, grid cells at the edges of the square are adjacent (on the sphere) to distant cells. We handle this by duplicating a few rows and columns near the edge of the square into an enlarged grid (suitably flipping them). Every sample point near an edge is placed into 2 grid cells representing the same hexagon on the sphere (and points near the corners are placed into 4 grid cells). Then checking new candidates is fast because we can just look into the cells beyond the edge and don’t need tricky index-wrapping logic. * Because spherical area coordinates are not completely uniform on the sphere, we need to size our grid so that in the trickiest areas (halfway along octant edges) a circle of radius ${tex`r`} at a grid vertex entirely covers the three adjacent hexagons. We need a ${tex`k \times k`} grid with ${tex`\quad k = 2\left\lceil \dfrac{\pi\bigl(1 + \sqrt2\bigr)}{3r} \right\rceil.`} ${hexagon_coverage} If we used a square grid, we would need ${tex`k`} to be 50% larger, for a grid with more than twice as many total cells. * Because of our grid non-uniformity we need to check every candidate point against hexagons up to 3 hops away in the grid, compared to 2 hops for the planar version of Bridson’s algorithm. (If we used a squares instead of hexagons in our spherical grid, we would need to look 4 hops.) In theory we could limit ourselves to examining the 37 hexagons within 3 hops, but this would need separate logic for each quadrant and some special-case logic for handling the edges between quadrants. Instead we simply examine the 49 cells in a 7 by 7 block. * Instead of measuring planar (squared) distance between points ${tex`\begin{aligned}\quad &\bigl(\langle x_B, y_B \rangle - \langle x_A, y_A \rangle\bigr){}^2 \\[5mu] &\qquad= \langle x_B - x_A, y_B - y_A \rangle^2 \\[5mu] &\qquad= (x_B - x_A)^2 + (y_B - y_A)^2 \\[5mu] &\qquad\leq r^2, \end{aligned}`} we need to measure [spherical distance, using the stereographic projection](/@jrus/planisphere#Distance). If ${tex`\rho`} is the (great-circle arclength) radius, ${tex`\begin{aligned} \quad&\bigl(\langle x_B, y_B \rangle \boxminus \langle x_A, y_A \rangle\bigr){}^2 \\[8mu] &\qquad=\dfrac{\langle x_B - x_A, y_B - y_A \rangle^2}{\bigl| 1 + x_Ax_B + y_Ay_B + I(x_Ay_B - y_Ax_B) \bigr|{}^2} \\[18mu] &\qquad= \dfrac{(x_B - x_A)^2 + (y_B - y_A)^2}{(1 + x_Ax_B + y_Ay_B)^2 + (x_Ay_B - y_Ax_B)^2} \\[18mu] &\qquad\leq s^2 = \tan^2 \tfrac12 \rho. \end{aligned}`} * Instead of moving new candidate points near the origin to the vicinity of a parent point via translation, we move them by [rotating the sphere](/@jrus/planisphere#Rotation_North). `}31:T16a1,function* SphereDiskSample(radius) { const k = 11; // maximum number of samples before rejection const m = 4; // a number mutually prime to k const eps = 1e-6; const rotx = Math.cos((2 * Math.PI * m) / k); const roty = Math.sin((2 * Math.PI * m) / k); const sradius = Math.tan(0.5 * radius); const sradius2 = sradius * sradius; const gridsize = Math.ceil((1 - eps) * Math.PI * (1 + Math.sqrt(2)) / 3 / radius); const grid_padded_rowsize = 2 * gridsize + 8; const grid = new Float64Array(8 * (gridsize + 4) * (gridsize + 4)).fill(NaN); const queue = []; // pick a random point near the north pole for the first sample const [Px, Py] = random_northern_point(); yield sample(Px, Py, far(Px, Py)); // pick a random existing sample from the queue. pick: while (queue.length) { const i = (Math.random() * queue.length) | 0; const parent = queue[i]; const t = tanpi_2(2 * Math.random() - 1); const q = 1 / (1 + t * t); let dw, dx, dy; dx = q ? (1 - t * t) * q : -1; // [dx, dy] = random unit vector dy = q ? 2 * t * q : 0; // make a new candidate for (let j = 0; j < k; ++j) { dw = dx; // temp name for dx dx = dw * rotx - dy * roty; dy = dw * roty + dy * rotx; const r = candidate_radius(); // stereo rotate a point near the north pole to be near the parent const ax = parent[0], ay = parent[1]; const bx = r * dx, by = r * dy; const px = ax + bx, py = ay + by; const qr = 1 - ax*bx - ay*by, qi = ax*by - ay*bx; const q2 = 1 / (qr*qr + qi*qi); if (q2 === Infinity) continue; // skip the point at infinity const x = q2 * (px * qr - py * qi); const y = q2 * (py * qr + px * qi); // accept candidates that are far enough from all existing samples const indices = far(x, y) if (indices) { yield sample(x, y, indices); continue pick; } } // if none of k candidates were accepted, remove the point from the queue const r = queue.pop(); if (i < queue.length) queue[i] = r; } function candidate_radius () { const rand0 = Math.random(); const rand1 = Math.random(); return sradius * (1 + eps + 0.65 * rand0 * rand1); } // fast approximation of tan(πx/2) function tanpi_2(a) { let b = 1 - a * a; return a * (-0.0187108 * b + 0.31583526 + 1.27365776 / b); } // random stereographic point near the north pole function random_northern_point() { const t = tanpi_2(2*Math.random() - 1), u = 0.25*Math.random(), v = Math.sqrt(u/(1 - u)), s = 1 / (1 + t*t); return (s !== 0) ? [(1 - t*t)*(v*s), 2*t*(v*s)] : [-v, 0]; } function stereo_compare (ax, ay, bx, by, rr) { let dx = bx - ax, dy = by - ay; let qr = 1 + ax*bx + ay*by, qi = ax*by - ay*bx; return (dx*dx + dy*dy) <= rr * (qr*qr + qi*qi) } function far(x, y) { const stereoX = x, stereoY = y; // COMPUTE SAC QUINCUNCIAL PROJECTION GRID CELL const sx = (x >= 0) - (x < 0); const sy = (y >= 0) - (y < 0); const s = sx * sy; // reflect into positive quadrant x *= sx, y *= sy; const r_2 = 1 / (x*x + y*y); const rbig = r_2 < 1; if ((r_2 !== 0) | (r_2 !== Infinity)) { // invert large values across unit circle x *= r_2*rbig + !rbig; y *= r_2*rbig + !rbig; } else { x = 0, y = 0; } // project fundamental octant let xt = x; x = (4 / Math.PI) * Math.atan2(xt, 1 + y); y = (4 / Math.PI) * Math.atan2(y, 1 + xt); // TODO: investigate fast approximations for atan2 // invert across x + y = 1 xt = x; x = !rbig * x + rbig * (1 - y); y = !rbig * y + rbig * (1 - xt); // round to nearest grid cell x *= gridsize, y *= gridsize; let xg = x | 0, yg = y | 0; // grid index const xr = x - xg, yr = y - yg; // remainder xg -= s * (s*(2*xr + yr) < 1.5*s - 0.5) * (s * (yr - xr) > 0); yg -= s * (s*(2*yr + xr) < 1.5*s - 0.5) * (s * (yr - xr) <= 0); // reflect back into negative quadrants xg = sx * xg - (sx < 0); yg = sy * yg - (sy < 0); // CHECK NEAREST 7x7 CELLS FROM THE EXISTING GRID const j0 = yg + gridsize + 1; // yg - 3 translated and padded const i0 = xg + gridsize + 1; // xg - 3 translated and padded for (let j = j0; j < j0 + 7; ++j) { const index0 = 2 * (j * grid_padded_rowsize + i0); for (let i = index0; i < index0 + 14; i += 2) { const gridx = grid[i]; if (gridx !== gridx) continue; // skip NaNs if (stereo_compare(stereoX, stereoY, gridx, grid[i+1], sradius2)) { return false; // too close to an existing point } } } return [xg, yg]; // far enough from every existing point } function populate_grid(i, j, x, y) { const jpad = j + gridsize + 4, ipad = i + gridsize + 4 const index = 2 * (grid_padded_rowsize * jpad + ipad); grid[index] = x; grid[index + 1] = y; } // flip across the edge either gridsize or -gridsize function flip(xg) { const sx = (xg > 0) - (xg <= 0); return ~(xg - sx*gridsize) + sx*gridsize; } function sample(x, y, [xg, yg]) { const w = gridsize; const sample_point = [x, y]; populate_grid(xg, yg, x, y); // also add the sample to grid cells in the padded area // representing the same part of the sphere, if applicable const xg_near_side = (xg > w - 4) | (xg <= -w + 4); const yg_near_side = (yg > w - 4) | (yg <= -w + 4); if (xg_near_side) { populate_grid(flip(xg), ~yg, x, y); if (yg_near_side) populate_grid(~flip(xg), ~flip(yg), x, y); } if (yg_near_side) populate_grid(~xg, flip(yg), x, y); queue.push(sample_point); return sample_point; } }32:T65f,stereo_quincuncial = { const stereo_quincuncial = function stereo_quincuncial(x, y) { const sx = Math.sign(x) || Math.sign(1/x); // -0 -> -1; +0 -> +1 const sy = Math.sign(y) || Math.sign(1/y); const r2 = x*x + y*y; x *= sx; y *= sy; // reflect into positive quadrant const rbig = r2 > 1; if (rbig) { // invert across unit circle if (r2 === Infinity) { x = 0; y = 0 } x *= 1/r2, y *= 1/r2; } // project fundamental octant const t = x; x = (4 / Math.PI) * Math.atan2(x, 1 + y); y = (4 / Math.PI) * Math.atan2(y, 1 + t); if (rbig) { const t = x; x = 1 - y; y = 1 - t; } return [x * sx, y * sy] }; stereo_quincuncial.invert = function stereo_quincuncial_inverse(x, y) { // reduce areas outside [-1, 1] square and reflect to positive quadrant const xg = Math.round(0.5*x); const yg = Math.round(0.5*y); x -= 2*xg; y -= 2*yg; let sx = Math.sign(x) || Math.sign(1/x); let sy = Math.sign(y) || Math.sign(1/y); x *= sx; y *= sy; const g = 2 * ((xg + yg) & 1); sx *= 1 - g; sy *= 1 - g; const rbig = x + y > 1; if (rbig) { const t = x; x = 1 - y; y = 1 - t; } // invert across line x + y = 1 // project fundamental triangle const tx = Math.tan(Math.PI/4 * x); const ty = Math.tan(Math.PI/4 * y); x = 1 / (1 - tx*ty) * (tx * (1 + ty)); y = 1 / (1 - tx*ty) * (ty * (1 + tx)); if (rbig) { const r = 1 / (x*x + y*y); if (r === Infinity) {x = y = 1;} x *= r; y *= r; } return [x * sx, y * sy]; }; return stereo_quincuncial; }33:[{"id":"437ea6624ac4b6da2f75e72a33da11706302fb87759aef5f6b22195bcc3871b7c94f0f9e8d74363a08f074d895f8e4c132f9e3a4f3286862c25d5a23f23f9987","href":"https://static.observableusercontent.com/files/437ea6624ac4b6da2f75e72a33da11706302fb87759aef5f6b22195bcc3871b7c94f0f9e8d74363a08f074d895f8e4c132f9e3a4f3286862c25d5a23f23f9987","name":"Quincuncial hexagons.png","type":"image/png","lastModified":1654154881243.536,"size":216348},{"id":"a4035d0f81dca69c14d66a9cb6d62aef2df937a79ef2e107617d593ff84a61f24bef45012fc31174816b80cab9935ae6e54d08f358ae95811feb7a7f38c0e01a","href":"https://static.observableusercontent.com/files/a4035d0f81dca69c14d66a9cb6d62aef2df937a79ef2e107617d593ff84a61f24bef45012fc31174816b80cab9935ae6e54d08f358ae95811feb7a7f38c0e01a","name":"Quincuncial squares@2.png","type":"image/png","lastModified":1654182213741.432,"size":77749},{"id":"698d3739fdbc22710118e88c093705ffea9372d6203d1afe178f0d6fc68ff4eec509e1d45241a778b5a8ececc91c62ab523160eb1c8611eeeb05f3b0b9860bd5","href":"https://static.observableusercontent.com/files/698d3739fdbc22710118e88c093705ffea9372d6203d1afe178f0d6fc68ff4eec509e1d45241a778b5a8ececc91c62ab523160eb1c8611eeeb05f3b0b9860bd5","name":"hexagon coverage.png","type":"image/png","lastModified":1654193632765.864,"size":312150},{"id":"e3f152756b24537a6cb5b6829d229c41a58848cf248a9cca442b9bef05799ad22fd40c7ad41a71c6abc45d275173e6e72537f164944dad74c74994aef0c3824f","href":"https://static.observableusercontent.com/files/e3f152756b24537a6cb5b6829d229c41a58848cf248a9cca442b9bef05799ad22fd40c7ad41a71c6abc45d275173e6e72537f164944dad74c74994aef0c3824f","name":"hexagon coverage@1.png","type":"image/png","lastModified":1654195358517.569,"size":318902}] 2c:{"id":"b0f5eca78588ea7a","slug":"spheredisksample","title":"Disk Sampling on a Sphere","thumbnail":"dca92e29a94aeced0d9cbc1a4c7dd311868b33f251e60e8607589cd2f0a76940","access_level":"public","create_time":"$D2022-05-30T01:55:52.734Z","update_time":"$D2022-06-06T22:30:04.973Z","view_time":null,"first_public_time":"$D2022-06-02T19:04:06.415Z","has_unpublished_edits":null,"fork":{"id":"53160c21ad15efa7","slug":"sac-quincuncial","title":"Spherical area coordinates quincuncial","owner":{"id":"9b8a8e177a49059a","login":"jrus","name":"Jacob Rus"}},"stars":19,"starred":false,"owner":{"type":"team","id":"9b8a8e177a49059a","login":"jrus","name":"Jacob Rus","avatar_url":"https://avatars-next.observableusercontent.com/c33dbeed76527c50e64208b0e005e0b6e54a3bef2a14ea6a823f1c7de4867662","tier":"starter_2024"},"version":3480,"publish_version":3480,"role":null,"authored":false,"authors":[{"type":"individual","id":"ddd7235babc79ce0","login":"jrus","name":"Jacob Rus","avatar_url":"https://avatars-next.observableusercontent.com/c33dbeed76527c50e64208b0e005e0b6e54a3bef2a14ea6a823f1c7de4867662","tier":null}],"members":[{"type":"individual","id":"ddd7235babc79ce0","login":"jrus","name":"Jacob Rus","avatar_url":"https://avatars-next.observableusercontent.com/c33dbeed76527c50e64208b0e005e0b6e54a3bef2a14ea6a823f1c7de4867662","tier":null,"role":"editor"}],"comments":[],"custom_thumbnail":true,"trashed":false,"lapsed":null,"subscribed":null,"forks":6,"collection":{"id":"6ad523ac4dcb4b45","slug":"instagram-shares","title":"Observable's Instagram Shares","owner":{"id":"f35c755083683fe5","login":"observablehq","name":"Observable"}},"type":"notebook","body":{"title":"Disk Sampling on a Sphere","cells":[{"id":0,"value":"Title = md`\n# Disk Sampling on a Sphere\n\n[Bridson’s algorithm](/@mbostock/poisson-disc-distribution) quickly generates a blue-noise distribution of points in the plane similar to Poisson disk sampling (randomly choosing points from a uniform distribution subject to the condition that each new point must be a certain minimum distance from any previous point). See Bridson (2007) [“Fast Poisson disk sampling in arbitrary dimensions”](https://doi.org/10.1145/1278780.1278807) and [Mike Bostock’s 2014 animation](https://bl.ocks.org/mbostock/dbb02448b0f93e4c82c3). The algorithm [can be tweaked slightly to run several times faster](/@fil/poisson-distribution-generators).\n\nIn this notebook we implement more or less the same idea on the sphere.\n\n`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":56,"value":"$2e","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":2899,"value":"viewof randomizeCircles = {\n return Inputs.button(\"Randomize and Re-animate\")\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":2093,"value":"viewof radius_slider = Inputs.range([1.5, 30], {value: 5, step: .1, label: \"Radius (degrees)\", transform: Math.log })","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":2693,"value":"$2f","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":3085,"value":"$30","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":3414,"value":"md`* * *\n\n### Speed\n\nIn comparison to the planar implementation, the spherical implementation is about 2–5 times slower, depending on the browser. Shown here is the time to generate around 5000 points; we run the code 5 times to show variation.`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":3425,"value":"Plot.barX(times, {\n x: \"time\",\n y: \"generator\",\n stroke: \"white\",\n fill: \"generator\",\n sort: \"time\",\n reverse: true\n}).plot({\n x: { nice: true, label: 'time (ms) →'},\n width\n})","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":3448,"value":"viewof replayTimes = {\n await visibility();\n return Inputs.button(\"replay\");\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":3428,"value":"import {C as PlaneDiskSample} from '@fil/poisson-distribution-generators'","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":3417,"value":"times = {\n replayTimes;\n const k = 5; // number of times to repeat each generator\n const timer = (f) => {\n const p = performance.now();\n f();\n return performance.now() - p;\n };\n let n = 5000;\n let S = () => SphereDiskSampleNPoints(n);\n let P = () => PlaneDiskSample(500, 500, n);\n return Promise.all(\n d3.shuffle(d3.range(k).flatMap(i => [S, P])).map(\n (f) =>\n new Promise(async (resolve) => {\n await Promises.delay(500 * Math.random());\n resolve({\n generator: f.name,\n time: timer(() => Array.from(f()))\n });\n })\n )\n );\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":3229,"value":"md`* * *`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":3014,"value":"radius = onchange(viewof radius_slider)","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":2520,"value":"circles = randomizeCircles, ([...SphereDiskSample(radius * Math.PI/180)])\n .map(stereo_inverse)\n .map((pt) =>\n d3.geoCircle().center(pt).radius(radius * .5).precision(10)()\n)","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":3257,"value":"function SphereDiskSampleNPoints(n) {\n return SphereDiskSample(2.949 / Math.sqrt(n));\n}","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":2390,"value":"$31","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":62,"value":"md`---`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":3010,"value":"// from Fabian Iwand, https://talk.observablehq.com/t/_/5112/7\nfunction onchange(input) {\n return Generators.observe(change => {\n const onchange = () => change(input.value);\n input.addEventListener('change', onchange);\n change(input.value);\n return () => input.removeEventListener('change', onchange);\n });\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":2528,"value":"stereo_inverse = ([x, y]) => {\n // invert stereographic projection\n const longitude = Math.atan2(y, x);\n const latitude = - 2 * Math.atan(Math.hypot(x, y)) + 0.5 * Math.PI;\n return [180/Math.PI * longitude, 180/Math.PI * latitude];\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":2395,"value":"hexgrid = {\n var map2 = (arr, cb) => arr.map(pts => pts.map(cb));\n // draw hexagon edges within one prototypical square:\n var square = [[[0, 0], [1/3, 1/3]], [[1, 0], [1/3, 1/3]], [[0, 1], [1/3, 1/3]]];\n return (n) => map2([].concat(\n ...d3.range(n).map(i => [].concat(\n ...d3.range(n).map(j => [].concat(\n map2(square, ([x, y]) => [i+x, j+y]), map2(square, ([x, y]) => [-i-x, -j-y]),\n map2(square, ([x, y]) => [i+1-x, -j-1+y]), map2(square, ([x, y]) => [-i-1+x, j+1-y]) )) ))\n ), ([x, y]) => [x/n, y/n]) }","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":2404,"value":"quincuncial_degrees = {\n var radians = Math.PI/180;\n var degrees = 180/Math.PI;\n var projection = (λ, φ) => quincuncial(λ*radians, φ*radians);\n projection.invert = (x, y) => {x = quincuncial.invert(x, y); return [x[0]*degrees, x[1]*degrees]}\n return projection;\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":2397,"value":"lonlat_hexgrid = inverseResampleJSON(quincuncial_degrees, 1)({\n type: 'MultiLineString',\n coordinates: hexgrid(gridsize_for_r(radius * Math.PI/180)) })","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":2217,"value":"gridsize_for_r = (r) => {\n return Math.ceil(Math.PI * (1 + Math.sqrt(2)) / 2 / r * 2/3)\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":2400,"value":"import {inverseResampleJSON} from '@jrus/sphere-resample'","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":1245,"value":"$32","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":788,"value":"quincuncial = {\n const {sqrt, sign, cos, sin, tan, atan, atan2, hypot, SQRT2} = Math;\n const HALFPI = 0.5 * Math.PI;\n \n let quincuncial = function quincuncial(longitude, latitude) {\n let x, y;\n \n // stereographic projection with north pole at the origin\n { const r = tan(0.5 * (HALFPI - latitude));\n x = cos(longitude) * r;\n y = sin(longitude) * r;\n }\n\n return stereo_quincuncial(x, y)\n };\n \n quincuncial.invert = function quincunxial_inverse(x, y) {\n const z = stereo_quincuncial.invert(x, y);\n x = z[0]; y = z[1];\n \n // invert stereographic projection\n const longitude = atan2(y, x);\n const latitude = - 2 * atan(hypot(x, y)) + HALFPI;\n\n z[0] = longitude, z[1] = latitude;\n return z;\n };\n\n return quincuncial;\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":106,"value":"// 'sphere' stream courtesy of @fil\n\nquincuncial_projection = {\n const p = d3.geoProjection(quincuncial),\n projectionStream = p.stream;\n\n p.stream = function(stream) {\n var rotate = p.rotate(),\n rotateStream = projectionStream(stream),\n sphereStream = (p.rotate([0, 0, 0]), projectionStream(stream));\n p.rotate(rotate);\n rotateStream.sphere = function() {\n var p0 = p.translate(), s = p.scale();\n stream.polygonStart();\n stream.lineStart();\n stream.point(p0[0] + s, p0[1] - s); // top-right\n stream.point(p0[0] - s, p0[1] - s); // top-left\n stream.point(p0[0] - s, p0[1] + s); // bottom-left\n stream.point(p0[0] + s, p0[1] + s); // bottom-right\n stream.lineEnd();\n stream.polygonEnd();\n };\n return rotateStream;\n };\n\n return p\n .rotate([0, 0, 0])\n .scale(250)\n .center([0, 90])\n .fitExtent([[0 , 0], [width - 0, width - 0]], { type: \"Sphere\" })\n .precision(0.05);\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":15,"value":"land50 = {\n const\n url = \"https://cdn.jsdelivr.net/npm/world-atlas@2/land-50m.json\",\n raw = await fetch(url),\n topo = await raw.json();\n return topojson.feature(topo, topo.objects.land);\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":1542,"value":"land110 = {\n const\n url = \"https://cdn.jsdelivr.net/npm/world-atlas@2/land-110m.json\",\n raw = await fetch(url),\n topo = await raw.json();\n return topojson.feature(topo, topo.objects.land);\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":2,"value":"topojson = require(\"topojson-client@3\")","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":7,"value":"d3 = require('d3', \"d3-geo@1\", \"d3-selection@1\", \"d3-drag@1\", \"d3-geo-projection@2\")","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":556,"value":"function drag(projection) {\n let v0, q0, r0;\n \n function dragstarted() {\n v0 = versor.cartesian(projection.invert([d3.event.x, d3.event.y]));\n q0 = versor(r0 = projection.rotate());\n }\n \n function dragged() {\n const v1 = versor.cartesian(projection.rotate(r0).invert([d3.event.x, d3.event.y]));\n const q1 = versor.multiply(q0, versor.delta(v0, v1));\n projection.rotate(versor.rotation(q1));\n }\n \n return d3.drag()\n .on(\"start\", dragstarted)\n .on(\"drag\", dragged);\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":558,"value":"versor = require(\"versor@0.0.3\")","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":1584,"value":"html`\n`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"}],"files":"$W33","stdlib":"1","titling":"auto"}} 2d:"$2c:body"