31:Tb38,md`# H Curve Job van der Zwan introduced me to the H-Curve [in this notebook](https://observablehq.com/@jobleonard/a-simple-algorithm-to-generate-h-curves?collection=@jobleonard/space-filling-curves). As he says, > Whenever strong locality is the reason to resort to space-filling curves (which is almost every case where space-filling curves are used, as far as I know), one might expect the H-curve to be an ideal candidate. However, the only other direct mention I found of this curve was in "A Note on Space-Filling Visualizations and Space-Filling Curves", a paper from 2005 by Martin Wattenberg. So why hasn't it been used more since it was discovered? > > Well, I suspect the reason is that the construction scheme in the original paper is rather cryptic. He then goes on to lay a creative geometric way of building H-Curves in a single pass. But while this algorithm is elegant, it relies on bitwise operations that can't be used in WebGL, and (I suspect?) is better suited to laying out multiple points in tandem than calculating the position of thousands of individual points one-at-a-time, which is what you need to do on WebGL. As he said, the authors' code is a bit inscrutable, and also specifies many of the positions manually. Plus, it uses recursion, which isn't allowed in the shader language. I wanted to try the same strategy I used for the [Hilbert Curve](https://observablehq.com/@bmschmidt/hilbert-curve-layout-for-webgl) of storing the rotation state in WebGL matrices. But the H-Curve presents some interesting problems. One is the question of the origin. Since the H-Curve is closed, you could start anywhere! If you are dividing into quads, it makes sense to consider the origin at one of the corners. But then it doesn't divide neatly up into powers of four, because each square is bisected diagonally. The second is the odd down-down-down-right pattern of the basic curve that makes it distinctive. At the smallest level, the H curve is twiddlier than Hilbert or [the Peano curves I made initially](https://observablehq.com/@bmschmidt/boustrophedon-curves?collection=@bmschmidt/space-filling-curves-in-webgl). Fortunately, these problems can help solve each other. To get up to a power of four, we can essentially consider the direction of the last unit (up or down) and the half of the triangle as two binary flags. Since swapping back and forth across the origin would lose local proximity, I adopt a convention that the origin point [0, 0] in 2d corresponds to the location of the point [0] in 1d, and that the upper right half of the square is in negative space. Here's the basic unit of the H curve--the left side of the H and the bottom. I have the last segment toggling because its direction depends on the orientation of the next bit-- it's only kind of part of the basic block. (This block is 9 points/8 segments long.) `32:T939,H_Curve = glsl` float modI(float a, float b) { float m=a-floor((a+0.5)/b)*b; return floor(m+0.5); } ${rotation_logic.value} ${H_Position.value} // Two necessary functions. vec2 H_Curve(float raw_ix) { // Start near the origin with an identity transform; vec2 xy = vec2(0., 0.); // Numbers greater than zero are bool negative = raw_ix < -.5; float sign = 1.; float abs_ix = abs(raw_ix); if (negative) { sign = -1.; abs_ix -= 1.; // Waste a cycle and do zero from both directions. // Could cause some problems? } mat2 rotation = mat2(sign, 0., 0., sign ); float MAX_ITER = 14.; // Is the point even or odd? Used to position the up or down direction. float jag = modI(abs_ix, 2.); float ix = floor(abs_ix/ 2.0 + 0.1); // Where we are in the L shape (0, 1, 2, 3); float a; for (float i = 16.; i > -.5; i -= 1.) { float rank_in_layer = floor((ix + .01) / (pow(4., i))); a = modI(rank_in_layer, 4.); // This will change a in place. if (rotation[0][0] * rotation[1][1] * sign < 0.) { a = 3. - a; } vec2 p = H_position(a, rotation); rotation = rotation * rotations(a); xy += (p) * pow(2., i - 1.); } /* Position the last element */ // No clue what's going on here. The rotation // functions look like they're a proxy for some kind of logic // gate or something, but I got it nearly by trial and error. float last_direction = rotation[1][1]; float is_rotated = rotation[1][1] * rotation[0][0]; if (is_rotated > 0.) { last_direction *= -1.; } if (jag < .5) {last_direction *= -1.;} xy += vec2(0., .25 * last_direction * sign); // Some magic involving 2 ** 15 and 2**16 to reshift // points somewhere close to the origin. // If not using 16 as the highest point in the loop, this will have to change, too. float y_shift = 0.; float x_shift = 0.; if (negative) { xy += vec2(pow(2., 16.) + pow(2., 15.) + x_shift, -pow(2., 16.) + y_shift); } else { xy += vec2(pow(2., 15.) + x_shift, -pow(2., 16.)); } return xy; } `33:T49a,copy_from_fbo = function(fbo) { regl({ profile: true, blend: { enable: true, func: { srcRGB: 'one', srcAlpha: 'one', dstRGB: 'one minus src alpha', dstAlpha: 'one minus src alpha', } }, frag: ` precision mediump float; varying vec2 uv; uniform sampler2D tex; uniform float wRcp, hRcp; void main() { vec4 here = texture2D(tex, uv); vec4 oneup = texture2D(tex, uv + vec2(wRcp, 0.)); vec4 twoup = texture2D(tex, uv + vec2(0., hRcp)); gl_FragColor = max(here, max(oneup, twoup)); } `, vert: ` precision mediump float; attribute vec2 position; varying vec2 uv; void main() { uv = 0.5 * (position + 1.0); gl_Position = vec4(position, 0, 1); } `, attributes: { position: () => fill_buffer }, depth: { enable: false }, count: 3, uniforms: { tex: fbo, wRcp: ({viewportWidth}) => 1.0 / viewportWidth, hRcp: ({viewportHeight}) => 1.0 / viewportHeight }, })() }34:[] 2f:{"title":"H-Curve in webgl","cells":[{"id":3562,"value":"md`# H-Curve in webgl`","pinned":false,"hidden":false,"mode":"ojs"},{"id":10,"value":"viewof regl = {\n let div = DOM.canvas(width * .75, width * .75 );\n let regl = wrapREGL(div);\n div.value = regl;\n return div;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1738,"value":"{while (true) {\n yield md`${d3.format(\",\")(Math.floor(z.points))} points`\n}}","pinned":false,"hidden":false,"mode":"ojs"},{"id":783,"value":"viewof animate = Checkbox([\"Animate?\"], {value: [\"Animate?\"]})","pinned":false,"hidden":false,"mode":"ojs"},{"id":260,"value":"viewof n_points = Range([0, 22], {label: \"Number of points (log)\", value: 16, step:.05})","pinned":false,"hidden":false,"mode":"ojs"},{"id":753,"value":"viewof primitive = Radio([\"points\", \"line strip\"],{value: \"line strip\", label: \"WebGL Primitive\"})","pinned":false,"hidden":false,"mode":"ojs"},{"id":0,"value":"$31","pinned":false,"hidden":false,"mode":"ojs"},{"id":3573,"value":"{\n while (true) {\n const c = DOM.context2d(100, 100)\n c.translate(10, 10);\n c.scale(20, 20)\n c.font = \"1px Arial\";\n c.fillText(\"↓\", 0, 0)\n c.fillText(\"↓\", 0, 2)\n c.fillText(\"↑\", 1, 2)\n c.fillText(\"↓\", 2, 2)\n \n c.fillText(\"↓\", 0, 1)\n c.fillText(\"→\", 0.2, 2.5)\n c.fillText(\"→\", 1.2, 1.5)\n \n if (Math.floor(now/500) % 2 == 1) {\n c.fillText(\"→\", 2.2, 2.5)\n } else {\n c.fillText(\"↓\", 2, 3)\n }\n return c.canvas\n }\n\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":3578,"value":"md`We can simplify this into two groups of four.\n\nThe red arrows here are notable in that they only ever go up and down;\nthis means the rotations will be much simpler with them.\n\nWe can consider each red arrow as a dyad.\n\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":3582,"value":"{\n while (true) {\n const c = DOM.context2d(120, 120)\n c.translate(20, 20);\n c.scale(30, 30)\n c.globalAlpha = 1\n c.font = \"1px Arial\";\n c.fillStyle = \"tomato\"\n\n c.fillText(\"↓\", 0, 0)\n c.fillText(\"↓\", 0, 2)\n c.fillText(\"↑\", 1, 2)\n c.fillText(\"↓\", 2, 2)\n c.fillStyle = \"black\"\n c.globalAlpha = 0.35\n c.fillText(\"↓\", 0, 1)\n c.fillText(\"→\", 0.2, 2.5)\n c.fillText(\"→\", 1.2, 1.5)\n \n if (Math.floor(now/500) % 2 == 1) {\n c.fillText(\"→\", 2.2, 2.5)\n } else {\n c.fillText(\"↓\", 2, 3)\n }\n return c.canvas\n }\n\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":3585,"value":"md`\nTaking just the red arrows, they make a simple L pattern from the origin.\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":3587,"value":"{\n const c = DOM.context2d(120, 120)\n c.translate(20, 20);\n c.scale(30, 30)\n c.globalAlpha = 1\n c.font = \"1px Arial\";\n c.fillStyle = \"tomato\"\n\n c.fillText(\"↓\", 0, 0)\n c.fillText(\"↓\", 0, 2)\n c.fillText(\"↑\", 1, 2)\n c.fillText(\"↓\", 2, 2)\n return c.canvas\n\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":3591,"value":"md`The H curve is composed of units of these. These groups of four tessellate together into a regular pattern.\n\nThis shape is 4x7; the upper right quadrant is identical and also 4x7. They interlock to form a 4x8 rectangle.\n\nSince each \"point\" here is actually two vertical points, that means that the following shape can fill an 8x8 square.\n\n`\n\n","pinned":false,"hidden":false,"mode":"ojs"},{"id":3594,"value":"{\n const c = DOM.context2d(320, 320)\n c.translate(20, 20);\n c.scale(30, 30)\n c.globalAlpha = 1\n c.font = \"1px Arial\";\n c.fillStyle = \"tomato\"\n\n c.fillText(\"↓\", 0, 0)\n c.fillText(\"↓\", 0, 2)\n c.fillText(\"↑\", 1, 2)\n c.fillText(\"↓\", 2, 2)\n \n c.fillStyle = \"steelblue\"\n c.translate(0, 5.5);\n c.rotate(Math.PI)\n c.scale(-1, 1)\n \n c.fillText(\"↑\", 0, 0)\n c.fillText(\"↑\", 0, 2)\n c.fillText(\"↓\", 1, 2)\n c.fillText(\"↑\", 2, 2)\n\n c.fillStyle = \"purple\"\n\n c.translate(3.5, 1.5);\n c.rotate(Math.PI)\n c.beginPath()\n c.fillText(\"↑\", 0, 0)\n c.fillText(\"↑\", 0, 2)\n c.fillText(\"↓\", 1, 2)\n c.fillText(\"↑\", 2, 2)\n c.stroke()\n\n c.scale(1, -1)\n c.translate(-0.5, 0);\n c.rotate(Math.PI)\n\n c.fillStyle = \"forestgreen\"\n\n c.fillText(\"↓\", 0, 0)\n c.fillText(\"↓\", 0, 2)\n c.fillText(\"↑\", 1, 2)\n c.fillText(\"↓\", 2, 2)\n\n return c.canvas\n\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":3598,"value":"md`this shape is entirely understandable as a fractal.\nPulling back to think of these only as L shapes and considering the \ny orientation after entry, it is composed \nof a down piece, a down piece, an up piece, and a down piece.\n\nFrom here, it's trivial to extend out infinitely.\n\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":3601,"value":"{\n const { path } = d3\n function L (context) {\n context.moveTo(0, 0); \n context.lineTo(0, 2); \n context.lineTo(1, 2); \n context.lineTo(2, 2);\n context.lineTo(3, 2);\n return context\n }\n \n const c = DOM.context2d(100, 100)\n c.translate(10, 10);\n c.scale(10, 10);\n c.lineWidth = 2/10\n \n c.beginPath()\n L(c)\n c.stroke()\n \n c.translate(0, 5);\n c.rotate(Math.PI)\n c.scale(-1, 1)\n \n c.beginPath()\n L(c)\n c.stroke()\n\n c.translate(4, 2);\n c.rotate(Math.PI)\n c.beginPath()\n L(c)\n c.stroke()\n\n c.scale(1, -1)\n c.translate(-1, 0);\n c.rotate(Math.PI)\n\n c.beginPath()\n L(c)\n c.stroke()\n\n return c.canvas\n \n }","pinned":false,"hidden":false,"mode":"ojs"},{"id":1358,"value":"md`## WebGL notes\n\nThe general strategy here is to start with the position in the largest possible square, and then increment the position by each next-smaller square below that. For each of the four positions \\`a\\` in the L shape above, there is a different position and rotation matrix.\n\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1361,"value":"rotation_logic = glsl `\n mat2 rotations(float a) {\n if (a < 0.5) {\n return mat2( 1., 0.,\n 0., 1.);\n }\n if (a <= 1.5) {\n return mat2( 1. , 0.,\n 0., -1.);\n }\n if (a <= 2.5) {\n return mat2( -1., 0.,\n 0., 1.);\n }\n return mat2( 1., 0.,\n 0., 1.);\n}\n`","pinned":true,"hidden":false,"mode":"ojs"},{"id":1370,"value":"md`\n...similarly, the position logic places each of the four points on the curve equal distance from the origin.\nInside this loop, the rotation matrix is actually transformed after calculating the position of the point; that\nrotation will persist for the next smaller square.\n\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1372,"value":"H_Position = glsl`\nvec2 H_position(in float a, mat2 rot) {\n // input is in 0., 1., 2., 3.\n\n vec2 p;\n\n // Check both the rotations at once. This will be 2 or zero.\n \n if (a < 0.5) { p = vec2(-.5, 1.);}\n else if (a < 1.5) {p = vec2(-.5, -1.);}\n else if (a < 2.5) {p = vec2(0.5, -1.);}\n else {p = vec2(1.5, -1.);} \n return p * rot;\n }\n`\n","pinned":true,"hidden":false,"mode":"ojs"},{"id":1394,"value":"md`## implementation\n\nAs with Hilbert, those two functions make it possible to calculate any finite number of points.\n\nThe direction at which you enter the L (long or short) varies depending on rotation, so that is coded into the loop below.\n\n`","pinned":true,"hidden":false,"mode":"ojs"},{"id":136,"value":"$32","pinned":true,"hidden":false,"mode":"ojs"},{"id":20,"value":"canvas = {\n // This actual canvas is shunted to regl to use as the colorscale texture, which is kinda nifty.\n // const width = 512\n const canvas = DOM.canvas(width, 20)\n const ctx = canvas.getContext(\"2d\")\n for (let i = 0; i < width; i++) {\n ctx.fillStyle = d3.interpolateSinebow(i/width)\n\n ctx.fillRect(i, 0, 1, 20) \n }\n return canvas\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":51,"value":"vertex_shader = glsl`\n precision highp float;\n\n attribute float ix;\n\n uniform float u_time;\n uniform float u_max_ix;\n uniform float tick;\n uniform sampler2D u_color_buffer;\n uniform float u_separation;\n varying vec4 fill;\n\n \n ${H_Curve.value}\n\n void main() {\n vec2 p_null = vec2(0., 0.);\n gl_PointSize = 5.;\n float relative_pos = ix - floor(u_max_ix/2.);\n if (relative_pos < 0.) {relative_pos;}\n float scale = 1. / sqrt(u_max_ix) * 2.5;\n vec2 pos = (H_Curve(relative_pos) - p_null) * scale + vec2(-.99, .5);\n gl_Position = vec4(pos.x, pos.y, 0., 1.);\n float prop_through = ix / u_max_ix;\n float offset = -fract(u_time / 5.);\n fill = texture2D(u_color_buffer, vec2(fract(offset + ix/u_max_ix), 0.5));\n// float n_colors = 8.;\n// float m = 8.;\n// fill = texture2D(u_color_buffer, vec2(fract(modI(floor((ix)/m), n_colors) / n_colors), 0.5));\n// fill = texture2D(u_color_buffer, vec2(1000. * floor(ix/m), .5));\n\n }\n `","pinned":true,"hidden":false,"mode":"ojs"},{"id":3540,"value":"","pinned":true,"hidden":false,"mode":"ojs"},{"id":3646,"value":"fbo = regl.framebuffer({width: Math.floor(width * .75), height: Math.floor(width * .75)})\n","pinned":true,"hidden":false,"mode":"ojs"},{"id":3661,"value":"z.points","pinned":true,"hidden":false,"mode":"ojs"},{"id":83,"value":"draw_tick = {\n const tick = regl.frame(({time}) => {\n // clear contents of the drawing buffer\n regl.clear({\n color: [0, 0, 0, 0],\n depth: 1\n })\n drawPoints({points: z.points, rounded_points: Math.round(z.points), primitive: z.primitive}) \n /*\n fbo.use( () => {\n regl.clear({\n color: [0, 0, 0, 0],\n depth: 1\n })\n drawPoints({points: z.points, rounded_points: Math.round(z.points), primitive: z.primitive}) \n }) \n drawPoints({points: z.points, rounded_points: Math.round(z.points), primitive: z.primitive}) \n\n copy_from_fbo(fbo)\n })\n */\n })\n invalidation.then(() => tick.cancel())\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":3651,"value":"$33","pinned":true,"hidden":false,"mode":"ojs"},{"id":3652,"value":"fill_buffer = regl.buffer(\n {data: [ -4, -4, 4, -4, 0, 4 ]})\n ","pinned":true,"hidden":false,"mode":"ojs"},{"id":1309,"value":"p = {\n while (true) {\n const s = d3.scaleLinear().domain([-1, 1]).range([4, 18])\n if (animate.length) {\n z.points = 2 ** (s(Math.sin(now/10000)))\n } else {\n z.points = (Math.pow(2, n_points))\n }\n z.primitive = primitive\n yield z.points\n }\n}\n","pinned":true,"hidden":false,"mode":"ojs"},{"id":58,"value":"fragment_shader = glsl`\n precision mediump float;\n varying vec4 fill;\n \n\n void main() {\n gl_FragColor = fill;\n }\n `","pinned":true,"hidden":false,"mode":"ojs"},{"id":1264,"value":"z = ({\n rounded_points: 16,\n primitive: \"line\"\n})","pinned":true,"hidden":false,"mode":"ojs"},{"id":1450,"value":"z","pinned":true,"hidden":false,"mode":"ojs"},{"id":2,"value":"drawPoints = {\n const regl_params = {\n\n depth: { enable: false },\n stencil: { enable: false },\n blend: {\n enable: true,\n func: {\n srcRGB: 'one',\n srcAlpha: 'one',\n dstRGB: 'one minus src alpha',\n dstAlpha: 'one minus src alpha',\n },\n },\n frag: fragment_shader.value,\n vert: vertex_shader.value,\n attributes: {\n ix: {\n buffer: points,\n stride: 4,\n offset: 0 * 4\n },\n },\n uniforms: {\n tick: ({ tick }) => tick,\n u_max_ix: regl.prop(\"points\"),\n u_boust: 1, // Left older from previous version\n u_time: regl.context(\"time\"),\n u_aspect_ratio: 1,\n u_color_buffer: colorscheme,\n u_separation: 0,\n },\n // specify the number of points to draw from the buffer.\n count: (_, props) => props[\"rounded_points\"],\n primitive: regl.prop(\"primitive\")\n }\n\n return regl(regl_params)\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":41,"value":"points = regl.buffer(vals)","pinned":true,"hidden":false,"mode":"ojs"},{"id":543,"value":"vals = {\n const vals = new Float32Array(2**24)\n for (let i = 0; i < vals.length; i++) {\n vals[i] = i \n }\n return vals\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":16,"value":"colorscheme = regl.texture(canvas)","pinned":true,"hidden":false,"mode":"ojs"},{"id":5,"value":"wrapREGL = require('regl')","pinned":true,"hidden":false,"mode":"ojs"},{"id":7,"value":"import { glsl } from \"@stwind/glsl-chunk-tag\";","pinned":true,"hidden":false,"mode":"ojs"},{"id":17,"value":"d3 = require(\"d3@v6\")","pinned":true,"hidden":false,"mode":"ojs"},{"id":463,"value":"import {Checkbox, Radio, Range} from \"@observablehq/inputs\"\n","pinned":true,"hidden":false,"mode":"ojs"}],"files":"$W34","stdlib":"1","titling":"auto"}