31:T4a2,## Quadkeys "Quadkeys" seem to me kinda like numbers or strings or whatever you want them to be, just as long as they uniquely identify a quadtree node. The most common form I've seen is to assign strings `"0"`, `"1"`, `"2"`, `"3"` to the four children of each node. The "quadkey" is then a repeated string of numbers from 0-3 which identify which of the four child nodes to follow. For example, `"210123"` would be a node at the sixth level of the tree. Recursing into the tree means appending the appropriate number, e.g. `"2101232"`. Stepping up to a parent simply means chopping off the last character. The parent of `"210123"` is `"21012"`. A downside of common quadkeys is that it's not a perfectly effortless fit for map tiles. Web map tiles operate in terms of (x, y, z) tuples where x and y are raster coordinates and z is the zoom level. It's a pretty trivial matter to traverse the string and perform some powers of two to translate one into the other. It's a bit tedious though, especially in a weird place like postgres database trigger. I had an easier time thinking of quadtree operations as arithmetic operations on numbers rather than loop-wise operations on strings.32:T52b,## Quadtrees At last, the code below implements sort of a fake quadtree. That is, a quadtree with a size and shape, but without any actual data or persisted nodes. The use case was to perform quadtree logic in order to fetch map tiles, but without any locally-available data necessitating actual construction of a quadtree. This logic is a key piece of the [Barnes-Hut approximation](https://observablehq.com/@jheer/the-barnes-hut-approximation) in which we seek to traverse a tree, stopping whenever nodes are sufficiently small and distant. Whenever we stop and don't traverse into the children of a particular node, we use an aggregated approximation of its data rather than the child nodes' more accurate representation. This is common for gravitational simulations, to name one use, in which we replace a collection of masses with a single aggregated body at the center of mass—repeating this process up the tree until we're left with a single mass at the root node which approximates the entire collection of masses. Jeffrey Heer has [a really good visualization of the Barnes-Hut approximation](https://observablehq.com/@jheer/the-barnes-hut-approximation) as well as [a very nice accompanying explanation](https://jheer.github.io/barnes-hut/). My goal was to use it for another computational geometry algorithm.33:T1355,/* * Implements a quadtree-like structure without any actual nodes or data, * so that we can perform basic queries like asking "what are the nearby * tiles this point would need to query?" without actually requiring an * explicit tree. This allows us to precompute the tiles to request from * the server. */ class SyntheticQuadtree { constructor([xmin, ymin, xmax, ymax], zmax) { this.bounds = [xmin, ymin, xmax, ymax]; this.zmax = zmax; } /** * Returns lists of full and partial tiles required to evaluate a particular point. * * @param {number[]} [xq, yq] - x and y coordinates of queried point * @param {number[]} [zstart, zend] - start and end range zoom range of desired tiles * @param {number} theta - Barnes-Hut parameter. Smaller means more accuracy and more recursion * * @returns {object} - Returns an object containing fullTiles and partialTiles lists of quadkeys * fulfilling the desired query. Full tiles are those which should be queried in full since * the full geometry is necessary, while partialTiles are those tiles for which only * aggregated data is required. */ queryBarnesHutTilesForPoint([xq, yq], [zstart, zend] = [0, 0], theta = 1.0) { return this.queryBarnesHutTilesForAabb( [xq, yq, xq, yq], [zstart, zend], theta ); } /** * Returns lists of full and partial tiles required to evaluate a particular axis-aligned * bounding box. * * @param {number[]} aabb - axis-aligned bounding box [xmin, ymin, xmax, ymax] of desired region * @param {number[]} [zstart, zend] - See queryBarnesHutTilesForPoint * @param {number} theta - See queryBarnesHutTilesForPoint * * @returns {object} - See queryBarnesHutTilesForPoint */ queryBarnesHutTilesForAabb(aabb, [zstart, zend] = [0, 0], theta = 1.0) { const fullTiles = []; const partialTiles = []; this.traversePreorder(zend, (xyz) => { const env = this.xyzToAabb(xyz); const size = Math.hypot(env[2] - env[0], env[3] - env[1]); // Distance to node const dist = distToAabb( 0.5 * (env[0] + env[2]), 0.5 * (env[1] + env[3]), aabb ); // Be careful to behave well when dividing by a potential zero const noRecurse = !(dist <= 0 || size / dist >= theta) && xyz[2] >= zstart; if (noRecurse) { partialTiles.push(xyz); return false; } // Should recurse, so check if we can. if (xyz[2] < zend && xyz[2] < this.zmax) return true; // Otherwise want recursion but hit max z, so try requesting full tile if (xyz[2] < zend && xyz[2] === this.zmax) { // We *can* recurse but we're at zmax, so we should request the full tile data fullTiles.push(xyz); } else { // Otherwise we don't need recursion, so push the partial (aggregated) tile partialTiles.push(xyz); } return false; }); return { fullTiles, partialTiles }; } /** * Convert z/y/z tuple to an axis-aligned bounding box * * @param {number[]} [x, y, z] - x/y/z tuple specifying the WebMercator tile * * @returns {number[]} - axis-aligned bounding box */ xyzToAabb([x, y, z]) { const dim = 2 ** z; const dx = (this.bounds[2] - this.bounds[0]) / dim; const dy = (this.bounds[3] - this.bounds[1]) / dim; return [ this.bounds[0] + dx * x, this.bounds[1] + dy * y, this.bounds[0] + dx * (x + 1), this.bounds[1] + dy * (y + 1) ]; } /** * Convert a point to the x/y/z tuple of the tile in which it lies * * @param {number[]} [qx, qy] - query point coordinates * * @returns {number[]} - [x, y, z] tile tuple */ pointToXYZ([qx, qy], z = 0) { const dim = 2 ** z; const x = Math.floor( (dim * (qx - this.bounds[0])) / (this.bounds[2] - this.bounds[0]) ); const y = Math.floor( (dim * (qy - this.bounds[1])) / (this.bounds[3] - this.bounds[1]) ); return [x, y, z]; } /** * Perform preorder traversal of the tree, meaning that it always visits parent * tiles before children. * * @param {number} zMax - maximum traversal depth * @param {function} visitor - function invoked for each tile, receiving the [z, x, y] tuple. * If function returns false, recursionis halted. Otherwise recursion continues until zMax. */ traversePreorder(zMax, visitor) { const stack = []; let cur = [0, 0, 0]; while (stack.length || cur) { while (cur) { const descend = visitor(cur) !== false && cur[2] < zMax; const nx = cur[0] * 2; const ny = cur[1] * 2; const nz = cur[2] + 1; if (descend) { stack.push([nx + 1, ny + 1, nz]); stack.push([nx, ny + 1, nz]); stack.push([nx + 1, ny, nz]); } cur = descend && [nx, ny, nz]; } if (stack.length) cur = stack.pop(); } } }34:T50b,handleInteraction = { viewof regl.style.position = "relative"; function onMove({ offsetX, offsetY }) { const x = (offsetX / viewof regl.offsetWidth) * 2 - 1; const y = (1 - offsetY / viewof regl.offsetHeight) * 2 - 1; mutable queryPoint = [x, y]; } function onTouchmove(e) { const rect = e.target.getBoundingClientRect(); const offsetX = e.targetTouches[0].pageX - rect.left; const offsetY = e.targetTouches[0].pageY - rect.top; onMove({ offsetX, offsetY }); } function onTouchstart(event) { event.preventDefault(); event.stopPropagation(); onTouchmove(event); viewof regl.addEventListener("touchmove", onTouchmove); viewof regl.addEventListener("touchend", onTouchend); } function onTouchend(event) { viewof regl.removeEventListener("touchend", onTouchend); viewof regl.removeEventListener("touchmove", onTouchmove); } viewof regl.addEventListener("mousemove", onMove); viewof regl.addEventListener("touchstart", onTouchstart); invalidation.then(() => { viewof regl.removeEventListener("mousemove", onMove); viewof regl.removeEventListener("touchstart", onTouchstart); viewof regl.removeEventListener("touchmove", onTouchmove); viewof regl.removeEventListener("touchend", onTouchend); }); }35:T79b,drawBoxes = regl({ vert: ` precision highp float; attribute vec2 box; attribute vec4 aabb; varying float radius; varying vec4 color; varying vec2 uv; uniform vec2 range; uniform vec4 color_cast; uniform sampler2D colorscale; void main () { float sizeparam = (log2(abs(aabb.z - aabb.x)) - range.x) / (range.y - range.x); color = texture2D(colorscale, vec2(sizeparam, 0.5)); color = mix(color, color_cast, color_cast.a); uv = box; gl_Position = vec4(vec2( mix(aabb.x, aabb.z, box.x), mix(aabb.y, aabb.w, box.y) ), 0, 1); }`, frag: ` #extension GL_OES_standard_derivatives : enable precision highp float; uniform float stroke; varying float radius; varying vec2 uv; varying vec4 color; float grid (vec2 parameter, float width, float feather) { float w1 = width - feather * 0.5; vec2 d = fwidth(parameter); vec2 looped = 0.5 - abs(mod(parameter, 1.0) - 0.5); vec2 a2 = smoothstep(d * (w1 + feather), d * w1, looped); return max(a2.x, a2.y); } void main () { gl_FragColor = mix(vec4(color.rgb, 1) * color.a, vec4(0, 0, 0, 1), stroke * grid(uv, 1.0, 1.0)); }`, attributes: { box: [ [1, 0], [1, 1], [0, 0], [0, 1] ], aabb: { buffer: regl.prop("aabbs"), offset: 0, stride: 16, divisor: 1 } }, uniforms: { colorscale: regl.prop("colorscale"), range: (ctx, props) => [Math.log2(2 ** -10), Math.log2(2 ** 0)], stroke: (ctx, props) => (!!props.stroke ? 1 : 0), color_cast: (ctx, props) => props.cast || [0, 0, 0, 0] }, depth: { enable: false }, blend: { enable: true, func: { srcRGB: 1, srcAlpha: 1, dstRGB: "one minus src alpha", dstAlpha: "one minus src alpha" } }, primitive: "triangle strip", count: 4, instances: regl.prop("count") })36:[{"id":"f1048eb9b693a993c41b302c35c740b8e27d7fe6e1b93e629bdc0c91cba2543f98147b13a4c09535841510adc52affd0fb019062f4c4d28a1bc5dac601d8a604","href":"https://static.observableusercontent.com/files/f1048eb9b693a993c41b302c35c740b8e27d7fe6e1b93e629bdc0c91cba2543f98147b13a4c09535841510adc52affd0fb019062f4c4d28a1bc5dac601d8a604","name":"diagram.png","type":"image/png","lastModified":1743871402292.207,"size":129841}] 2f:{"title":"Quadkeys and Quadtrees","cells":[{"id":1982,"value":"# Quadkeys and Quadtrees","pinned":false,"hidden":false,"mode":"md"},{"id":1975,"value":"A quadtree is a geometric data structure in which the four quadrants of a rectangular node are recursively split into four child nodes. They have all sorts of uses, but the most applicable to me is web map tiles.","pinned":false,"hidden":false,"mode":"md"},{"id":1020,"value":"Web maps typically use the [Web Mercator projection](https://en.wikipedia.org/wiki/Web_Mercator_projection). When tiled, the top level map tile covers the entire world (at least up to ${tex`\\pm 85.06^\\circ`} latitude). Successive quadtree depths, or *zoom* levels, cover increasingly small square sections of the world. The top level tile is has zoom ${tex`z = 0`}. The four child tiles have ${tex`z = 1`} and are assigned raster coordinates ${tex`(x, y)`} = ${tex`(0, 0)`}, ${tex`(1, 0)`}, ${tex`(0, 1)`}, and ${tex`(1, 1)`}.","pinned":false,"hidden":false,"mode":"md"},{"id":1977,"value":"An important consideration then is how to uniquely address tiles and move around in this space. Let me be clear, this is *not* a challenge in the geospatial industry. You just call them ${tex`(x, y, z)`} tuples and that's the end of it. But that was marginally inconvenient for my needs and presented a tiny little challenge for me when implementing more out-there computation geometry algorithms on web map tiles.","pinned":false,"hidden":false,"mode":"md"},{"id":1979,"value":"This notebook puts down on paper an obvious and kinda fun numbering scheme I used.","pinned":false,"hidden":false,"mode":"md"},{"id":952,"value":"$31","pinned":false,"hidden":false,"mode":"md"},{"id":967,"value":"{\n const img = await FileAttachment(\"diagram.png\").image();\n img.style.maxWidth = \"100%\";\n img.style.width = \"500px\";\n img.style.display = \"block\";\n img.style.margin = \"1em auto\";\n return html`

${img}`;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":955,"value":"The natural answer, then, was to assign serial numbers to quadtree nodes. The scheme starts at ${tex`z=0`} and assigns serial integers to all tiles of each zoom level—in raster order—before proceeding to the next zoom level. The single tile at ${tex`z = 0`} has quadkey ${tex`0`}. The four tiles at ${tex`z=1`} have quadkeys ${tex`1, 2, 3, 4`}, and so on, proceeding next to quadkey ${tex`5`} at ${tex`z = 2`}. To compute the quadkey for a given tile, we therefore need to know where to start numbers at a particular ${tex`z`}. That is, how many _total_ tiles there are for _all_ ${tex`z`} less than a given zoom.","pinned":false,"hidden":false,"mode":"md"},{"id":957,"value":"The number of tiles at a given zoom ${tex`z`} is ${tex`2^{2z}`}, or ${tex`4^z`}. Thus zoom levels ${tex`0, 1, 2, 3, 4, \\ldots`} contain ${tex`1, 4, 16, 64, 256, \\ldots`} tiles, respectively. The cumulative sum of tiles up to zoom Z is ${tex`1, 5, 21, 85, 341, \\ldots`}. Although there's a cute geometric proof, a search of the Online Encyclopedia of Integer Sequences quickly gets us to a general equation for sequence [A002450](https://oeis.org/A002450):","pinned":false,"hidden":false,"mode":"md"},{"id":960,"value":"$${tex.block`\n\\sum_{z=0}^{Z-1} 4^z = \\frac{4^Z - 1}{3}\n`}","pinned":false,"hidden":false,"mode":"md"},{"id":962,"value":"This equation represents _the total number of tiles up until a given zoom ${tex`Z`}_.","pinned":false,"hidden":false,"mode":"md"},{"id":1375,"value":"To convert back and forth between this quadkey and an (x, y, z) tuple, we compute how many tiles come before it at lower zoom levels, subtract that from the quadkey, and then work out the raster coordinates from the remainder. For example, quadkey 18 is at zoom level 2. So we subtract the total number of tiles at zoom level 0 (1) and zoom level 1 (4). Computing 18 - 4 - 1 leaves a raster index of 13. Working out the raster coordinates with some division and a modulo operation, we have x=1, y=3, and z=2.","pinned":false,"hidden":false,"mode":"md"},{"id":2102,"value":"Mathematically, we write for quadkey ${tex`q,`}\n\n${tex.block`q = \\frac{4^Z - 1}{3} + 2^Z Y + X,`}","pinned":false,"hidden":false,"mode":"md"},{"id":2108,"value":"and we invert the logic to write\n\n${tex.block`\n\\begin{aligned}\nZ &= \\lfloor \\log_2 \\frac{3q + 1}{2} \\rfloor \\\\\nY &= \\lfloor \\frac{q - \\frac{4^Z-1}{3}}{2^Z} \\rfloor \\\\\nX &= q - 2^Z Y. \\\\\n\\end{aligned}\n`}","pinned":false,"hidden":false,"mode":"md"},{"id":2091,"value":"What happened next is where this all got a little bit regrettable. This is pretty simple to implement in JavaScript. Since there are so many powers of two involved, it's quite easy to speed up a lot of the operations with bitwise operations. However, 32-bit precision limits us to roughly zoom z=14. The slow version uses double precision floating point, which gets us up through z=25 at the cost of being about 10x slower.","pinned":false,"hidden":false,"mode":"md"},{"id":964,"value":"This is all way over-optimized for such a silly function, but the result is that it's robust, pretty fast, and perfectly adequate for reasonable zoom ranges.","pinned":false,"hidden":false,"mode":"md"},{"id":297,"value":"function xyzToQuadkey(xyz) {\n if (xyz[2] > 14) {\n // Double precision variant (except this first line is safe as bitwise)\n let a = 1 << xyz[2];\n return (a * a - 1) / 3 + xyz[1] * a + xyz[0];\n }\n // 32-bit integer variant\n return ((1 << (xyz[2] << 1)) - 1) / 3 + xyz[1] * (1 << xyz[2]) + xyz[0];\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":435,"value":"quadkeyToXYZ = {\n const MAX_32BIT_MODE_KEY = 2505397588;\n\n return function quadkeyToXYZ(out, key) {\n let a, z;\n if (key > MAX_32BIT_MODE_KEY) {\n // Double-precision variant.\n // log2 is the particularly slow double-precision operation.\n z = Math.floor(Math.log2(3 * key + 1) / 2);\n a = 1 << z;\n key -= (a * a - 1) / 3;\n out[1] = Math.floor(key / a);\n out[0] = key - a * out[1];\n } else {\n // 32-bit variant\n z = (31 - Math.clz32(key * 3 + 1)) >>> 1;\n key -= ((1 << (z << 1)) - 1) / 3;\n out[0] = key & ((1 << z) - 1);\n out[1] = key >>> z;\n }\n out[2] = z;\n return out;\n };\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1995,"value":"A simple back-and-forth looks like this:","pinned":false,"hidden":false,"mode":"md"},{"id":2047,"value":"key = 123456","pinned":true,"hidden":false,"mode":"ojs"},{"id":1807,"value":"xyz = quadkeyToXYZ([], key)","pinned":true,"hidden":false,"mode":"ojs"},{"id":1990,"value":"xyzToQuadkey(xyz)","pinned":true,"hidden":false,"mode":"ojs"},{"id":1999,"value":"Pushing the limits, we can use this up to around z=25:","pinned":false,"hidden":false,"mode":"md"},{"id":1771,"value":"xyzToQuadkey(quadkeyToXYZ([], 987654321098767))","pinned":true,"hidden":false,"mode":"ojs"},{"id":987,"value":"Computing the parent tile requires going back and forth between two representations with a bit of math. It's a bit more work than chopping off a single character. Is this a good or a bad thing? I don't know. It's a thing. And for what I was doing, I found it more useful and natural and easy to work with.","pinned":false,"hidden":false,"mode":"md"},{"id":1850,"value":"Below are a few more variants for different numeric formats. The benchmarking is a bit finicky and takes significantly longer after re-evaluating these cells as opposed to execution when the page loads fresh. I don't know why. And `BigInt` is just way too slow to use.","pinned":false,"hidden":false,"mode":"md"},{"id":1307,"value":"function xyzToQuadkey_DOUBLE(xyz) {\n let a = Math.pow(2, xyz[2]);\n return (a * a - 1) / 3 + xyz[1] * a + xyz[0];\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1298,"value":"function xyzToQuadkey_32BIT(xyz) {\n return ((1 << (xyz[2] << 1)) - 1) / 3 + xyz[1] * (1 << xyz[2]) + xyz[0];\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1658,"value":"function xyzToQuadkey_BIGINT(xyz) {\n const x = BigInt(xyz[0]);\n const y = BigInt(xyz[1]);\n const z = BigInt(xyz[2]);\n const a = 1n << z;\n return Number((a * a - 1n) / 3n + y * a + x);\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1304,"value":"function quadkeyToXYZ_DOUBLE(out, key) {\n const z = Math.floor(Math.log2(3 * key + 1) / 2);\n const w = 2 ** z;\n const i = key - (w * w - 1) / 3;\n out[0] = i % w;\n out[1] = Math.floor(i / w);\n out[2] = z;\n return out;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1301,"value":"function quadkeyToXYZ_32BIT(out, key) {\n const z = (31 - Math.clz32(key * 3 + 1)) >>> 1;\n key -= ((1 << (z << 1)) - 1) / 3;\n out[0] = key & ((1 << z) - 1);\n out[1] = key >>> z;\n out[2] = z;\n return out;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1583,"value":"quadkeyToXYZ_BIGINT = {\n function clz64(x) {\n if (x <= 0n) return 64;\n let n = 0n;\n let bit = 9223372036854775808n;\n while (bit > 0n && (x & bit) === 0n) {\n n++;\n bit >>= 1n;\n }\n return n;\n }\n\n return function quadkeyToXYZ_BIGINT(out, key) {\n const bkey = BigInt(key);\n const z = (63n - clz64(3n * bkey + 1n)) >> 1n;\n const i = bkey - ((1n << (z << 1n)) - 1n) / 3n;\n out[0] = Number(i % (1n << z));\n out[1] = Number(i >> z);\n out[2] = Number(z);\n return out;\n };\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":2097,"value":"You can toggle benchmarks below. For reasons I don't understand, it's faster to run benchmarks when the notebook first loads than to toggle the selector below and then run them. 🤷\n","pinned":false,"hidden":false,"mode":"md"},{"id":1262,"value":"viewof numericFormat = Inputs.radio(\n [\"Double-precision\", \"32-bit integer\", \"BigInt\"],\n {\n value: \"32-bit integer\",\n label: \"Benchmark\"\n }\n)","pinned":false,"hidden":false,"mode":"ojs"},{"id":1061,"value":"{\n rerun;\n\n const toXYZ = {\n \"Double-precision\": quadkeyToXYZ_DOUBLE,\n \"32-bit integer\": quadkeyToXYZ_32BIT,\n BigInt: quadkeyToXYZ_BIGINT\n }[numericFormat];\n const toQuadkey = {\n \"Double-precision\": xyzToQuadkey_DOUBLE,\n \"32-bit integer\": xyzToQuadkey_32BIT,\n BigInt: xyzToQuadkey_BIGINT\n }[numericFormat];\n\n const N = 1e6;\n const xyz = new Uint32Array(3);\n let computedKey;\n\n const t0 = performance.now();\n for (let inputKey = 0; inputKey < N; inputKey++) {\n toXYZ(xyz, inputKey);\n computedKey = toQuadkey(xyz);\n if (inputKey !== computedKey)\n throw new Error(`Incorrect conversion for key=${inputKey}`);\n }\n const t1 = performance.now();\n\n return `${N.toLocaleString()} round-trip quadkey <-> xyz conversions computed in ${(\n t1 - t0\n ).toFixed(1)}ms (${Math.floor(\n (1e3 * N) / (t1 - t0)\n ).toLocaleString()} ops/sec)`;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1293,"value":"viewof rerun = Inputs.button(\"Re-run\")","pinned":false,"hidden":false,"mode":"ojs"},{"id":1715,"value":"testKey = 1234567891234567 // fails with 32-bit\n//testKey = 2505397588 // 32-bit safe","pinned":true,"hidden":false,"mode":"ojs"},{"id":1383,"value":"xyzToQuadkey_DOUBLE(quadkeyToXYZ_DOUBLE([], testKey))","pinned":true,"hidden":false,"mode":"ojs"},{"id":1709,"value":"xyzToQuadkey_32BIT(quadkeyToXYZ_32BIT([], testKey))","pinned":true,"hidden":false,"mode":"ojs"},{"id":1461,"value":"xyzToQuadkey_BIGINT(quadkeyToXYZ_BIGINT([], testKey))","pinned":true,"hidden":false,"mode":"ojs"},{"id":0,"value":"$32","pinned":false,"hidden":false,"mode":"md"},{"id":1048,"value":"A class with some basic operations as well as an visual example is implemented below. The primary purpose of this visualization is to confirm that my computation is correct and results in non-overlapping nodes. To confirm this, you may disable the colorscale and observe that it has a uniform color. If nodes were double-counted, transparency would result in darker regions wherever there is overlap. (To introduce an error, try adding a trailing `+ 1` to [`xyzToQuadkey`](#xyzToQuadkey).)\n","pinned":false,"hidden":false,"mode":"md"},{"id":1852,"value":"The four zoom parameters below are a bit cryptic but relate to the particular operation I was trying to perform. I was fetching full data for adjacent nodes and aggregated data for all non-adjacent beyond a certain distance. And I wasn't just querying a point, I was computing data for an entire tile. To be honest, it all worked pretty well, though the geometric computation I worked on never quite made it into production. Sorry to be a bit cryptic about it all. *C'est la vie*.","pinned":false,"hidden":false,"mode":"md"},{"id":6,"value":"$33","pinned":false,"hidden":false,"mode":"ojs"},{"id":851,"value":"viewof zmax = Inputs.range([0, 10], {\n label: html`z_max`,\n value: 7,\n step: 1\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":309,"value":"viewof zstart = Inputs.range([0, 6], {\n label: html`z_start`,\n step: 1,\n value: 0\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":237,"value":"viewof zend = Inputs.range([0, 10], {\n label: html`z_end`,\n step: 1,\n value: 8\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":195,"value":"viewof theta = Inputs.range([0.1, 2], {\n label: \"Barnes-Hut parameter, θ\",\n step: 0.01,\n value: 0.8\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":688,"value":"viewof queryType = Inputs.radio([\"Point\", \"Tile\"], {\n label: \"Query type\",\n value: \"Tile\"\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":662,"value":"viewof zQuery = Inputs.range([0, 10], {\n label: html`z_query`,\n step: 1,\n value: 5,\n disabled: queryType === \"Point\"\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":15,"value":"viewof regl = reglCanvas(this, {\n pixelRatio: devicePixelRatio,\n width: Math.min(width, 512),\n height: Math.min(width, 512),\n attributes: { antialias: false, depthStencil: false },\n extensions: [\n \"ANGLE_instanced_arrays\",\n \"OES_element_index_uint\",\n \"OES_standard_derivatives\"\n ]\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":413,"value":"html`

Move the mouse (or touch/drag) to modify the query point.

`","pinned":false,"hidden":false,"mode":"ojs"},{"id":355,"value":"viewof config = Inputs.checkbox([\"Stroke\", \"Colorscale\"], {\n label: \"Drawing config\",\n value: [\"Stroke\", \"Colorscale\"]\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":10,"value":"tree = new SyntheticQuadtree([-1, -1, 1, 1], zmax)","pinned":false,"hidden":false,"mode":"ojs"},{"id":511,"value":"mutable queryPoint = [-0.45, 0.33]","pinned":false,"hidden":false,"mode":"ojs"},{"id":772,"value":"queryTile = queryType === \"Point\" ? null : tree.pointToXYZ(queryPoint, zQuery)","pinned":false,"hidden":false,"mode":"ojs"},{"id":914,"value":"quadkeys = {\n let aabb;\n if (queryType === \"Point\") {\n aabb = [...queryPoint, ...queryPoint];\n } else {\n aabb = tree.xyzToAabb(queryTile);\n }\n queryAabbBuffer.subdata(aabb);\n\n const { fullTiles, partialTiles } = tree.queryBarnesHutTilesForAabb(\n aabb,\n [zstart, zend],\n theta\n );\n\n return {\n full: fullTiles.map(xyzToQuadkey).sort((a, b) => a - b),\n partial: partialTiles.map(xyzToQuadkey).sort((a, b) => a - b)\n };\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":128,"value":"$34","pinned":false,"hidden":false,"mode":"ojs"},{"id":463,"value":"## Rendering\nThe remaining code uses the [regl](https://github.com/regl-project/regl) library to render the quadtree. It uses instanced rendering to permit direct input of packed `float32` axis-aligned bounding boxes (*AABB's*) with format ${tex`[x_{min}, y_{min}, x_{max}, y_{max}]`}. For each box instance, it renders lines connecting five vertices to create a square outline, or a triangle strip connecting four vertices to create a filled quad.","pinned":false,"hidden":false,"mode":"md"},{"id":13,"value":"import { reglCanvas } from \"@rreusser/regl-canvas\"","pinned":false,"hidden":false,"mode":"ojs"},{"id":106,"value":"fullAabbsBuffer = regl.buffer(new Float32Array(2 ** 21))","pinned":false,"hidden":false,"mode":"ojs"},{"id":829,"value":"partialAabbsBuffer = regl.buffer(new Float32Array(2 ** 21))","pinned":false,"hidden":false,"mode":"ojs"},{"id":698,"value":"queryAabbBuffer = regl.buffer(new Float32Array([0, 0, 0.5, 0.5]))","pinned":false,"hidden":false,"mode":"ojs"},{"id":20,"value":"repaint = {\n const fullAabbsData = quadkeys.full\n .map((key) => quadkeyToXYZ([], key))\n .map((node) => tree.xyzToAabb(node));\n fullAabbsBuffer.subdata(new Float32Array(fullAabbsData.flat()));\n\n const partialAabbsData = quadkeys.partial\n .map((key) => quadkeyToXYZ([], key))\n .map((node) => tree.xyzToAabb(node));\n partialAabbsBuffer.subdata(new Float32Array(partialAabbsData.flat()));\n\n regl.poll();\n\n drawBoxes([\n {\n aabbs: fullAabbsBuffer,\n count: fullAabbsData.length,\n colorscale,\n zstart,\n zend,\n cast: [0, 0, 255, 0.5],\n stroke: ~config.indexOf(\"Stroke\")\n },\n {\n aabbs: partialAabbsBuffer,\n count: partialAabbsData.length,\n colorscale,\n zstart,\n zend,\n stroke: ~config.indexOf(\"Stroke\")\n }\n ]);\n\n if (queryType === \"Tile\") {\n drawBoxes({\n aabbs: queryAabbBuffer,\n count: 1,\n colorscale: queryColorscale,\n stroke: false\n });\n }\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":22,"value":"$35","pinned":false,"hidden":false,"mode":"ojs"},{"id":248,"value":"d3 = require(\"d3@7\")","pinned":false,"hidden":false,"mode":"ojs"},{"id":704,"value":"queryColorscale = {\n const texture = regl.texture({\n width: 1,\n height: 1,\n data: new Uint8Array([255, 50, 100, 200])\n });\n invalidation.then(() => texture.destroy());\n return texture;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":250,"value":"colorscale = {\n const texture = regl.texture({\n width: 256,\n height: 1,\n min: \"linear\",\n mag: \"linear\",\n data: d3\n .quantize(d3.interpolateRainbow, 256)\n .map((c) => {\n if (~config.indexOf(\"Colorscale\")) {\n return (c = d3.rgb(c)), [c.r, c.g, c.b, 128];\n } else {\n return [255, 128, 128, 100];\n }\n })\n .flat()\n });\n invalidation.then(() => texture.destroy());\n return texture;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":651,"value":"function distToAabb(px, py, [xmin, ymin, xmax, ymax]) {\n const dx = Math.max(xmin - px, 0, px - xmax);\n const dy = Math.max(ymin - py, 0, py - ymax);\n return Math.sqrt(dx * dx + dy * dy);\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":144,"value":"import { inspect, Inspector } from \"@observablehq/inspector\"","pinned":false,"hidden":false,"mode":"ojs"}],"files":"$W36","stdlib":"1","titling":"auto"}