30:T485,FIND = { const m = map(projection), render = m.render; let mouse = null, found; m.render = function() { render(); const context = m.getContext("2d"), path = d3.geoPath(m.projection).context(context); context.beginPath(), path({ type: "MultiLineString", coordinates: edges.map(l => l.map(i => points[i])) }), (context.lineWidth = 1), (context.strokeStyle = "lime"), context.stroke(); if (mouse) { context.beginPath(), path({ type: "Point", coordinates: mouse }), (context.fillStyle = "red"), context.fill(); } if (found !== undefined) { context.beginPath(), path( d3 .geoCircle() .center(points[found]) .radius(radius)() ), (context.lineWidth = 2), (context.strokeStyle = "green"), context.stroke(); } context.lineWidth = 1; }; m.render(); d3.select(m).on("mousemove", function() { mouse = m.projection.invert(d3.mouse(m)); found = find(...mouse, radius); m.render(); }); return m; }31:T96d,polygons = { const polygons = []; const centers = circumcenters.slice(); // supplementary centers for degenerate cases like n = 1,2,3 const supplements = []; if (triangles.length === 0) { if (points.length < 2) return (polygons.centers = centers), polygons; if (points.length === 2) { // two hemispheres const a = cartesian(points[0]), b = cartesian(points[1]), m = normalize(vecsum(a, b)), d = normalize(cross(a, b)), c = cross(m, d); const poly = [ m, cross(m, c), cross(cross(m, c), c), cross(cross(cross(m, c), c), c) ] .map(spherical) .map(supplement); return ( polygons.push(poly), polygons.push(poly.slice().reverse()), (polygons.centers = centers), polygons ); } } triangles.forEach((tri, t) => { for (let j = 0; j < 3; j++) { const a = tri[j], b = tri[(j + 1) % 3], c = tri[(j + 2) % 3]; polygons[a] = polygons[a] || []; polygons[a].push([b, c, t, [a, b, c]]); } }); // reorder each polygon const reordered = polygons.map(poly => { const p = [poly[0][2]]; // t let k = poly[0][1]; // k = c for (let i = 1; i < poly.length; i++) { // look for b = k for (let j = 0; j < poly.length; j++) { if (poly[j][0] == k) { k = poly[j][1]; p.push(poly[j][2]); break; } } } if (p.length > 2) { return p; } else if (p.length == 2) { const R0 = o_midpoint( points[poly[0][3][0]], points[poly[0][3][1]], centers[p[0]] ), R1 = o_midpoint( points[poly[0][3][2]], points[poly[0][3][0]], centers[p[0]] ); const i0 = supplement(R0), i1 = supplement(R1); return [p[0], i1, p[1], i0]; } else { // I don't think we'll ever reach this(?) console && console.warn({ here: "unreachable", poly }); } }); function supplement(point) { let f = -1; centers.slice(triangles.length, Infinity).forEach((p, i) => { if (p[0] === point[0] && p[1] === point[1]) f = i + triangles.length; }); if (f < 0) (f = centers.length), centers.push(point); return f; } // prevents a side effect in degenerate cases reordered.centers = centers; return reordered; }32:[] 2e:{"title":"Geo Delaunay","cells":[{"id":0,"value":"md`# Geo Delaunay\n\n\n## UPDATE\n\n**This code is now available as [d3-geo-voronoi@1](https://github.com/Fil/d3-geo-voronoi). Happy tessellation!**\n\n--------\n\nUsing [d3-delaunay](https://github.com/d3/d3-delaunay) to build a geographic Delaunay/Voronoi diagram. We project the points from the sphere to the plane with a stereographic projection, adding [“fake points”](/@fil/geo-mesh-fake-points) on the infinite horizon.\n\nThis appears to be _much_ more scalable than the method we used in d3-geo-voronoi (See [#1](https://github.com/Fil/d3-geo-voronoi/issues/1)).\n\nThis is a work in progress, please send comments etc. the usual way!`","pinned":false,"hidden":false,"mode":"ojs"},{"id":205,"value":"n = 400","pinned":false,"hidden":false,"mode":"ojs"},{"id":952,"value":"projection = \"EqualEarth\"","pinned":false,"hidden":false,"mode":"ojs"},{"id":225,"value":"function geo_delaunay_from(points) {\n if (points.length < 2) return {};\n\n const r = d3.geoRotation(points[0]),\n projection = d3\n .geoStereographic()\n .translate([0, 0])\n .scale(1)\n .rotate(r.invert([180, 0]));\n points = points.map(projection);\n\n const zeros = [0];\n let max2 = 1;\n for (let i = 1, n = points.length; i < n; i++) {\n let m = points[i][0] * points[i][0] + points[i][1] * points[i][1];\n if (isNaN(m)) zeros.push(i);\n if (m > max2) max2 = m;\n }\n const FAR = 1e6 * Math.sqrt(max2);\n\n zeros.forEach((_, i) => (points[i] = [FAR / 2, 0]));\n\n // infinite horizon points\n for (let i = 0; i < 4; i++) {\n points.push([\n FAR * Math.cos((i / 2) * Math.PI),\n FAR * Math.sin((i / 2) * Math.PI)\n ]);\n }\n\n const delaunay = d3.Delaunay.from(points);\n delaunay.projection = projection;\n return delaunay;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":187,"value":"delaunay = geo_delaunay_from(points)","pinned":false,"hidden":false,"mode":"ojs"},{"id":557,"value":"md`### Edges\na list of Delaunay edges. No overlap. No half-edges (everything is connected).`","pinned":false,"hidden":false,"mode":"ojs"},{"id":109,"value":"edges = {\n const edges = [];\n const { halfedges, triangles } = delaunay;\n if (!halfedges) return edges;\n\n const seen = {};\n for (let i = 0, n = halfedges.length; i < n; ++i) {\n const j = halfedges[i];\n if (j < i) continue;\n let [a, b] = d3.extent([triangles[i], triangles[j]]);\n if (b >= points.length && a < points.length) (b = a), (a = 0);\n if (\n b > 0 &&\n b < points.length &&\n (a > 0 || (!seen[b]++ && (seen[b] = true)))\n )\n edges.push([a, b]);\n }\n\n return edges;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":710,"value":"EDGES = {\n const m = map(projection),\n render = m.render;\n m.render = function() {\n render();\n const context = m.getContext(\"2d\"),\n path = d3.geoPath(m.projection).context(context);\n context.beginPath(),\n path({\n type: \"MultiLineString\",\n coordinates: edges.map(l => l.map(i => points[i]))\n }),\n (context.strokeStyle = \"orange\"),\n context.stroke();\n };\n m.render();\n return m;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1158,"value":"md`### Find\n\n_radius_ in degrees.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1197,"value":"function find(x, y, radius) {\n let cell,\n dist,\n found = 0,\n next = 0;\n\n radius *= radians;\n do {\n cell = next;\n next = null;\n dist = d3.geoDistance([x, y], points[cell]);\n neighbors[cell].forEach(i => {\n var ndist = d3.geoDistance([x, y], points[i]);\n if (ndist < dist) {\n dist = ndist;\n next = i;\n found = i;\n return;\n }\n });\n } while (next !== null);\n\n if (radius == null || dist < radius) return found;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1222,"value":"radius = 10","pinned":false,"hidden":false,"mode":"ojs"},{"id":1160,"value":"$30","pinned":false,"hidden":false,"mode":"ojs"},{"id":705,"value":"md`### Triangles\n\nWe take the non-trivial triangles from the projected delaunay; to draw them, ignore those that cover more than a hemisphere (but we'll need them later for their centers, in order to build the Voronoi).`","pinned":false,"hidden":false,"mode":"ojs"},{"id":703,"value":"triangles = {\n if (!delaunay.triangles) return [];\n const triangles = Array.from(delaunay.triangles).map(\n d => (d >= points.length ? 0 : d)\n );\n const tri = [];\n for (let i = 0, n = triangles.length / 3; i < n; i++) {\n const a = triangles[3 * i],\n b = triangles[3 * i + 1],\n c = triangles[3 * i + 2];\n if (a !== b && b !== c && c !== a) {\n tri.push([a, c, b]);\n }\n }\n return tri;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":8,"value":"TRIANGLES = {\n const m = map(projection),\n render = m.render;\n m.render = function() {\n render();\n const context = m.getContext(\"2d\"),\n path = d3.geoPath(m.projection).context(context);\n\n triangles.forEach((tri, i) => {\n const p = {\n type: \"Polygon\",\n coordinates: [\n [points[tri[0]], points[tri[1]], points[tri[2]], points[tri[0]]]\n ]\n };\n \n if (d3.geoArea(p) > 2*Math.PI) return;\n \n context.beginPath(),\n path(p),\n (context.strokeStyle = context.fillStyle = d3.hcl(\n (360 * i) / triangles.length,\n 60,\n 70,\n 0.6\n )),\n context.fill(),\n context.stroke();\n });\n\n if (points.length < 100) {\n context.fillStyle = \"#000\";\n points.forEach((p, i) => context.fillText(i, ...m.projection(p)));\n }\n };\n m.render();\n return m;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":788,"value":"md`### Urquhart\n\nRemove the largest edge from each triangle.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1106,"value":"urquhart = {\n const _index = {};\n edges.forEach(e => (_index[`${e[0]}-${e[1]}`] = true));\n\n triangles.forEach(tri => {\n var l = 0,\n length = 0,\n remove,\n e;\n for (var j = 0; j < 3; j++) {\n e = d3.extent([tri[j], tri[(j + 1) % 3]]);\n length = d3.geoDistance(points[e[0]], points[e[1]]);\n if (length > l) {\n l = length;\n remove = `${e[0]}-${e[1]}`;\n }\n }\n _index[remove] = false;\n });\n\n const urquhart = [];\n for (let p in _index) {\n if (_index[p]) urquhart.push(p.split(\"-\").map(Number));\n }\n return urquhart;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1104,"value":"URQUHART = {\n const m = map(projection),\n render = m.render;\n m.render = function() {\n render();\n const context = m.getContext(\"2d\"),\n path = d3.geoPath(m.projection).context(context);\n context.beginPath(),\n path({\n type: \"MultiLineString\",\n coordinates: urquhart.map(l => l.map(i => points[i]))\n }),\n (context.strokeStyle = \"red\"),\n context.stroke();\n };\n m.render();\n return m;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":790,"value":"md`### Centers\n\nThe triangles’ circumcenters.\n\nFormula : V = a × b + b × c + c × a\n\n(In some applications, centroids might be preferred. We’re not talking about the Voronoi tessellation anymore, but why not…)\n\nNote that this list _may_ contain more points than there are triangles, so that we can still build the Voronoi tesselation in degenerate cases (example: with n=3 points, there are 2 triangles and 5 centers).\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":854,"value":"circumcenters = {\n if (!use_centroids) {\n return triangles.map(tri => {\n const c = tri.map(i => points[i]).map(cartesian),\n V = vecsum(cross(c[1], c[0]), cross(c[2], c[1]), cross(c[0], c[2]));\n return spherical(normalize(V));\n });\n } else {\n return triangles.map(tri => {\n return d3.geoCentroid({\n type: \"MultiPoint\",\n coordinates: tri.map(i => points[i])\n });\n });\n }\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1655,"value":"// in degenerate cases (n=2, n=3) extra centers are added after the circumcenters\ncenters = polygons.centers","pinned":false,"hidden":false,"mode":"ojs"},{"id":1057,"value":"CIRCUMCENTERS = {\n const m = map(projection),\n render = m.render;\n m.render = function() {\n render();\n const context = m.getContext(\"2d\"),\n path = d3.geoPath(m.projection).context(context);\n\n context.beginPath(),\n path({\n type: \"MultiLineString\",\n coordinates: edges.map(l => l.map(i => points[i]))\n }),\n (context.strokeStyle = \"#aaa\"),\n context.stroke();\n\n context.beginPath(),\n path({\n type: \"MultiPoint\",\n coordinates: circumcenters\n }),\n (context.fillStyle = \"#acc\"),\n context.fill();\n };\n m.render();\n return m;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1047,"value":"viewof use_centroids = md``","pinned":true,"hidden":false,"mode":"ojs"},{"id":1126,"value":"md`### Neighbors(i)\n\nGo through all the triangles [a,b,c], and list _b_ as neighbor of _a_, _c_ as a neighbor of _b_, and _a_ as a neighbor of _c_.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1128,"value":"neighbors = {\n const neighbors = [];\n triangles.forEach((tri, i) => {\n for (let j = 0; j < 3; j++) {\n const a = tri[j],\n b = tri[(j + 1) % 3],\n c = tri[(j + 2) % 3];\n neighbors[a] = neighbors[a] || [];\n neighbors[a].push(b);\n }\n });\n\n // degenerate cases\n if (triangles.length === 0) {\n if (points.length === 2) (neighbors[0] = [1]), (neighbors[1] = [0]);\n else if (points.length === 1) neighbors[0] = [];\n }\n\n return neighbors;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":972,"value":"md`### Polygons\n\nGo through all the triangles, and add their centers to each of their points’ polygons. Then order each polygon so it is walked clockwise.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":807,"value":"$31","pinned":false,"hidden":false,"mode":"ojs"},{"id":1593,"value":"function o_midpoint(a, b, c) {\n a = cartesian(a);\n b = cartesian(b);\n c = cartesian(c);\n const s = Math.sign(dot(cross(a, b), c));\n return spherical(normalize(vecsum(a, b)).map(d => -s * d));\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":830,"value":"POLYGONS = {\n const m = map(projection),\n render = m.render;\n m.render = function() {\n render();\n const context = m.getContext(\"2d\"),\n path = d3.geoPath(m.projection).context(context);\n\n polygons.forEach((poly, i) => {\n if (!poly) return;\n\n const p = {\n type: \"Polygon\",\n coordinates: [[...poly, poly[0]].map(i => centers[i])]\n };\n\n context.beginPath(),\n path(p),\n (context.strokeStyle = context.fillStyle = d3.hcl(\n (360 * i) / (1 + polygons.length),\n 60,\n 70,\n 0.6\n )),\n context.fill(),\n context.stroke();\n });\n\n context.beginPath(),\n path.pointRadius(0.5),\n path({\n type: \"MultiPoint\",\n coordinates: points\n }),\n (context.strokeStyle = \"#000\"),\n context.stroke();\n\n if (points.length < 10) {\n context.fillStyle = \"#000\";\n points.forEach((p, i) => context.fillText(i, ...m.projection(p)));\n }\n };\n m.render();\n return m;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1034,"value":"md`### Voronoi Mesh\n\na little variant on the above…`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1070,"value":"mesh = {\n const mesh = [];\n polygons.forEach(poly => {\n if (!poly) return;\n let p = poly[poly.length - 1];\n for (let q of poly) {\n if (q > p) mesh.push([p, q]);\n p = q;\n }\n });\n return mesh;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1038,"value":"MESH = {\n const m = map(projection),\n render = m.render;\n m.render = function() {\n render();\n const context = m.getContext(\"2d\"),\n path = d3.geoPath(m.projection).context(context);\n\n context.beginPath(),\n path({\n type: \"MultiLineString\",\n coordinates: mesh.map(l => l.map(i => centers[i]))\n }),\n (context.strokeStyle = \"lime\"),\n context.stroke();\n };\n m.render();\n return m;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1681,"value":"md`#### Hull\n\nFollow the edges that we see only once when running all the small triangles. (Click on the checkbox below to have a proper dataset that does not cover the whole globe).`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1816,"value":"viewof see_hull = md``","pinned":true,"hidden":false,"mode":"ojs"},{"id":1818,"value":"hull = {\n const _hull = {},\n hull = [];\n triangles.map(tri => {\n const p = {\n type: \"Polygon\",\n coordinates: [\n [points[tri[0]], points[tri[1]], points[tri[2]], points[tri[0]]]\n ]\n };\n\n if (d3.geoArea(p) > 2 * Math.PI) return;\n for (let i = 0; i < 3; i++) {\n let e = [tri[i], tri[(i + 1) % 3]],\n code = `${e[1]}-${e[0]}`;\n if (_hull[code]) delete _hull[code];\n else _hull[e.join(\"-\")] = true;\n }\n });\n\n const _index = {};\n let start;\n Object.keys(_hull).forEach(e => {\n e = e.split(\"-\").map(Number);\n _index[e[0]] = e[1];\n start = e[0];\n });\n\n if (start === undefined) return hull;\n\n let next = start;\n do {\n hull.push(next);\n next = _index[next];\n } while (next !== start);\n\n return hull;\n}\n\n\n\n","pinned":false,"hidden":false,"mode":"ojs"},{"id":1820,"value":"HULL = {\n const m = map(projection),\n render = m.render;\n m.render = function() {\n render();\n const context = m.getContext(\"2d\"),\n path = d3.geoPath(m.projection).context(context);\n\n context.beginPath(),\n path(\n hull.length === 0\n ? null /*{ type: \"Sphere\" }*/\n : {\n type: \"Polygon\",\n coordinates: [[...hull.map(i => points[i]), points[hull[0]]]]\n }\n ),\n ((context.fillStyle = \"pink\"), (context.strokeStyle = \"red\")),\n context.fill(),\n context.stroke();\n\n context.fillStyle = \"brown\";\n hull.forEach(i => context.fillText(i, ...m.projection(points[i])));\n };\n m.render();\n return m;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":772,"value":"md`### Check standard hull\nOn the sphere, the standard delaunay hull should be empty! Check that, in our computation, it only contains fake points.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":424,"value":"emptyhull = {\n let { hull } = delaunay,\n node = hull;\n const H = [];\n\n if (!node) return [];\n\n do {\n H.push(node.i);\n node = node.next;\n } while (node !== hull);\n\n return H.filter(i => i < points.length);\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1235,"value":"md`## TODO\n\nif we are to replace d3-geo-voronoi:\n- [x] mesh of edges\n- [x] triangles\n- [x] urquhart\n- [x] circumcenters\n- [x] neighbors(i)\n- [x] voronoi polygon(i)\n- [x] voronoi mesh\n- [x] find\n- [x] hull\n- [x] tests for edge cases (FAR should be large enough when points are close together)\n- [x] degenerate cases n=1 (Sphere), n=2 (Hemispheres), n=3 (the 2 centers are antipodal)\n- [x] remove the side effect in polygons= (which affects centers)\n- [ ] API https://beta.observablehq.com/d/8f757e3c06690145 & compatibility\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1004,"value":"md`-----\n\n_boring zone_`","pinned":false,"hidden":false,"mode":"ojs"},{"id":392,"value":"points = d3\n .range(n)\n .map((_, i) => [\n (see_hull ? 55 : 360) * (0.5 - Math.random()),\n (see_hull ? 25 : 90) * (Math.random() - Math.random())\n ])","pinned":false,"hidden":false,"mode":"ojs"},{"id":10,"value":"import {map} from \"@fil/base-map\"","pinned":false,"hidden":false,"mode":"ojs"},{"id":6,"value":"d3 = require(\"d3\", \"d3-delaunay@2\")","pinned":false,"hidden":false,"mode":"ojs"},{"id":862,"value":"// Converts 3D Cartesian to spherical coordinates (degrees).\nfunction spherical(cartesian) {\n return [\n Math.atan2(cartesian[1], cartesian[0]) * degrees,\n Math.asin(Math.max(-1, Math.min(1, cartesian[2]))) * degrees\n ];\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":865,"value":"// Converts spherical coordinates (degrees) to 3D Cartesian.\nfunction cartesian(coordinates) {\n var lambda = coordinates[0] * radians,\n phi = coordinates[1] * radians,\n cosphi = Math.cos(phi);\n return [cosphi * Math.cos(lambda), cosphi * Math.sin(lambda), Math.sin(phi)];\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1398,"value":"function dot(a, b) {\n for (var i = 0, n = a.length, s = 0; i < n; ++i) s += a[i] * b[i];\n return s;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":860,"value":"function cross(a, b) {\n return [\n a[1] * b[2] - a[2] * b[1],\n a[2] * b[0] - a[0] * b[2],\n a[0] * b[1] - a[1] * b[0]\n ];\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":881,"value":"function vecsum(...args) {\n return args[0].map((d, i) => args.reduce((a, b) => a + b[i], 0));\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":1395,"value":"function norm(v) {\n return Math.sqrt(dot(v, v));\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1394,"value":"function normalize(v) {\n var n = norm(v);\n if (n === 0) return [0,0,0];\n return v.map(d => d / n);\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":875,"value":"degrees = 1 / radians","pinned":false,"hidden":false,"mode":"ojs"},{"id":871,"value":"radians = Math.PI / 180","pinned":false,"hidden":false,"mode":"ojs"}],"files":"$W32","stdlib":"1","titling":"auto"}