31:T65d, md`This notebook walks through simulating a strange attractor on the GPU and then rendering particle tracks as lines. We do this without any of the data ever touching the CPU! The notebook presents the details below, but the high level summary is that the state of all particle tracks is stored in a single WebGL texture where each row of the texture is a ring buffer representing the history of each particle. On each frame we step a differential equation and overwrite the oldest element of each particle's history with a new value. Finally, we render particle tracks directly from the texture data. Each frame of the simulation incurs just four draw calls: 1. Integrate the differential equation, writing a column containing the new state of all particles into a temporary texture. 2. Copy the temporary slice back into full state texture, overwriting the oldest values. 3. Draw all line segments and line joins using instanced triangle strip geometry. Sentinel values separate the lines into multiple contiguous segments. 4. Draw all line caps, using instanced triangle geometry. Doing this requires being fancy about how we move data around, but I think the end result is actually pretty extensible and reusable for both 2D and 3D high performance line rendering. The line rendering itself is implemented in the [regl-gpu-lines](https://github.com/rreusser/regl-gpu-lines) npm module. This notebook walks through all of the above steps in detail. To see more interesting attractors in action, see the companion notebook, [Strange Attractors on the GPU, Part 2: Fun! 🔗](https://observablehq.com/d/df2e1a97f7ab9a46).`32:T16c1,{ const w = Math.min(600, width); const h = 150; const svg = d3.select(DOM.svg(w, h)); const colors = ['#cc2255', '#2255dd']; const xScale = d3 .scaleLinear() .domain([-2.5, 2.5]) .range([0, w]); const yScale = d3 .scaleLinear() .domain([-0.5, 0.9]) .range([h, 0]); const line = d3 .line(d => xScale(d[0]), d => yScale(d[1])) .defined(d => !!d && !isNaN(d[0])); const pointsPerSeg = 5; const points = [ [-2, -0.2], [-1.5, 0.5], [-1, 0.2], [-0.5, 0.2], [0, 0.6], null, [0, -0.2], [0.5, -0.1], [1, 0.2], [1.5, 0.3], [2, 0.1] ]; svg .append('path') .attr('stroke', colors[0]) .attr('fill', 'none') .attr('d', line(points.slice(0, pointsPerSeg))); svg .append('path') .attr('stroke', colors[1]) .attr('fill', 'none') .attr('d', line(points.slice(pointsPerSeg + 1))); svg .append('g') .selectAll('circle') .data(points.filter(d => d)) .enter() .append('circle') .attr('cx', d => xScale(d[0])) .attr('cy', d => yScale(d[1])) .attr('r', 2.5) .attr('fill', (d, i) => colors[Math.floor(i / pointsPerSeg)]); const entry = svg .append('g') .selectAll('text') .data(points.filter(d => d)) .enter() .append('text') .attr('x', d => xScale(d[0])) .attr('y', d => yScale(d[1])) .attr('dy', d => 7) .text('p') .attr('dominant-baseline', d => 'hanging') .attr('font-family', 'KaTeX_Main,Times New Roman,serif') .attr('text-anchor', d => 'start') .attr('font-weight', 700) .attr('fill', (d, i) => colors[Math.floor(i / pointsPerSeg)]); entry .append('tspan') .attr('font-weight', 400) .attr('font-size', '.8em') .attr('dy', '-.4em') .text((d, i) => `(${Math.floor(i / pointsPerSeg)})`); entry .append('tspan') .attr('font-weight', 400) .attr('font-size', '.8em') .attr('dy', '1.4em') .attr('dx', '-1.1em') .text((d, i) => i % pointsPerSeg); const particles = 3; const steps = 5; const texture = html`

← particle i

step ${tex`j`} →
${[...Array(particles).keys()].map( (particle, i) => { return html`${[ ...Array(steps).keys() ].map(step => { if (particle === 2) return html``; let eqn; const linenum = Math.floor(step / (pointsPerSeg + 1)); eqn = tex`\mathbf{p}_${step}^{(${particle})}`; return html``; })}`; } )}

${tex`\vdots`}

${eqn}

`; const endpoints = [0, 1, 2, 4, 3, 2, 0, 1, 2, 4, 3, 2]; const lookupTable = html` ${[ ...Array(points.length).keys() ].map(step => { const p = Math.floor(step / (pointsPerSeg + 1)); const i = step % (pointsPerSeg + 1); if (i === pointsPerSeg) return html``; const css = `color:${colors[p]}`; return html``; })}

vertices

${tex`\textsf{NaN}`}

${tex`uv_{${i}}^{(${p})}`}

${endpoints.map( (endpoint, endpointIndex) => { const seg = Math.floor(endpointIndex / 3); const line = Math.floor(endpointIndex / 6); const css = `color:${colors[line]};background-color:${ Math.floor(endpointIndex / 3) % 2 === 0 ? 'none' : '#eee' }`; return html``; } )} ${[1, 0, 1, 0].map( (isStart, idx) => { const css = `color:${ colors[Math.floor(idx / 2)] };background-color:${isStart ? 'none' : '#eee'}`; return html``; } )}

endpoints	${tex`uv_{${endpoint}}^{(${line})}`}
isstart	${tex`${isStart}`}

`; const css = html``; return html`

${svg.node()}

Line geometry

${texture}

State texture representation

${css}

The layout of the ${tex`j^{th}`} step of particle ${tex`i`}, ${tex`p_j^{(i)}`}, in a WebGL texture.

`; }33:Tcef,{ const w = Math.min(400, width); const h = (w * 3) / 4; const svg = d3.select(DOM.svg(w, h)); const xScale = d3 .scaleLinear() .domain([(-5 * w) / h, (5 * w) / h]) .range([0, w]); const yScale = d3 .scaleLinear() .domain([-5, 5]) .range([h, 0]); const line = d3.line(d => xScale(d[0]), d => yScale(d[1])).defined(x => !!x); const axes = svg.append('g'); axes .append('path') .attr('stroke', 'black') .attr('d', line([[-6, 0], [6, 0]])); axes .append('path') .attr('stroke', 'black') .attr('d', line([[0, -5], [0, 5]])); const quad = [[-1, -1], [1, -1], [1, 1], [-1, 1]]; svg .append('path') .attr('stroke', '#888888') .attr('fill', 'rgba(0,0,0,0.05)') .attr('d', line(quad) + 'Z'); const res = 12; const gridValues = [...Array(res - 1).keys()].map( i => ((i + 1) / res) * 2 - 1 ); const yGrid = gridValues.map(y => [[-1, y], [1, y], null]).flat(); const xGrid = gridValues.map(x => [[x, -1], [x, 1], null]).flat(); svg .append('path') .attr('stroke', 'rgba(0,0,0,0.2)') .attr('fill', 'none') .attr('d', line(xGrid)); svg .append('path') .attr('stroke', 'rgba(0,0,0,0.2)') .attr('fill', 'none') .attr('d', line(yGrid)); const bigTriangle = [[-4, -4], [4, -4], [0, 4]]; svg .append('path') .attr('stroke', 'rgba(128,0,0,0.5)') .attr('fill', 'rgba(255,0,0,0.05)') .attr('d', line(bigTriangle) + 'Z'); svg .append('g') .selectAll('circle') .data(bigTriangle) .enter() .append('circle') .attr('cx', d => xScale(d[0])) .attr('cy', d => yScale(d[1])) .attr('r', 2.5) .attr('fill', 'red'); svg .append('g') .selectAll('circle') .data(quad) .enter() .append('circle') .attr('cx', d => xScale(d[0])) .attr('cy', d => yScale(d[1])) .attr('r', 2.5) .attr('fill', 'black'); const labels = [ { text: '(4, -4)', loc: [4, -4], baseline: 'hanging', anchor: 'start' }, { text: '(-4, -4)', loc: [-4, -4], baseline: 'hanging', anchor: 'end' }, { text: '(0, 4)', loc: [0, 4], baseline: 'auto', anchor: 'start' }, { text: '(-1, -1)', loc: [-1, -1], baseline: 'hanging', anchor: 'end' }, { text: '(1,1)', loc: [1, 1], baseline: 'auto', anchor: 'start' } ]; svg .append('g') .selectAll('text') .data(labels) .enter() .append('text') .attr('x', d => xScale(d.loc[0])) .attr('y', d => yScale(d.loc[1])) .attr('dx', d => (d.anchor === 'start' ? 5 : -5)) .attr('dy', d => d.dy || 0) .text(d => d.text) .attr('dominant-baseline', d => d.baseline || 'auto') .attr('text-anchor', d => d.anchor || 'start'); return html`

${svg.node()}

A single big triangle covers the entire ${tex`[-1, 1] \times [-1, 1]`} viewport in the ${tex`x-y`} plane of Normalized Device Coordinates (NDC), triggering processing of all fragments. This forms our basic map operation for texture-based computation.

`; }34:T5ed,initializeState = regl({ vert: ` precision mediump float; attribute vec2 xy; varying vec2 uv; void main () { uv = 0.5 * xy + 0.5; gl_Position = vec4(xy, 0, 1); }`, frag: ` precision mediump float; varying vec2 uv; uniform vec2 resolution; uniform vec3 origin; uniform float scale; uniform sampler2D src; // http://extremelearning.com.au/unreasonable-effectiveness-of-quasirandom-sequences/ vec3 quasirandom (float n) { const float g = 1.22074408460575947536; return fract(0.5 + n * vec3(1.0 / g, 1.0 / (g * g), 1.0 / (g * g * g))).zyx; } vec3 sphericalRandom (float n) { vec3 rand = quasirandom(n); float u = rand.x * 2.0 - 1.0; float theta = (2.0 * ${Math.PI}) * rand.y; return vec3(sqrt(1.0 - u * u) * vec2(cos(theta), sin(theta)), u) * sqrt(rand.z); } void main () { // gl_FragCoord is offset by a half-texel, but it's actually nice to avoid // identically zero states, since they tend to get stuck along zero // axes and follow unstable paths off to infinity. float particle = gl_FragCoord.y; // Output (x, y, z, 1): gl_FragColor = vec4(origin + scale * sphericalRandom(particle), 1); }`, uniforms: { resolution: ctx => [ctx.framebufferWidth, ctx.framebufferHeight], origin: (ctx, props) => [0, 1, 0], scale: (ctx, props) => 1 }, attributes: { xy: [[-4, -4], [0, 4], [4, -4]] }, framebuffer: regl.prop('dst'), count: 3 })35:T6dd,integrate = regl({ vert: ` precision highp float; attribute vec2 xy; varying vec2 srcTexCoord; uniform float srcTexCoordU; void main () { // Map clip coords ([-1, 1] x [-1, 1]) to texture coords // ([0, 1] x [0, 1]) and store in a varying to be used // in the fragment shader: srcTexCoord = vec2(srcTexCoordU, 0.5 * xy.y + 0.5); gl_Position = vec4(xy, 0, 1); }`, frag: ` precision highp float; uniform sampler2D srcTex; uniform float dt, t; varying vec2 srcTexCoord; ${attractorGLSL} vec3 deriv (vec3 p, float t) { // Unpack, for convenience return derivative(p.x, p.y, p.z, t); } void main () { vec3 p = texture2D(srcTex, srcTexCoord).xyz; // Runge-Kutta 4th order integration vec3 k1 = deriv(p, t); vec3 k2 = deriv(p + (0.5 * dt) * k1, t + 0.5 * dt); vec3 k3 = deriv(p + (0.5 * dt) * k2, t + 0.5 * dt); vec3 k4 = deriv(p + dt * k3, t + dt); // Evaluate the derivative there and use for a whole step: gl_FragColor = vec4(p + (dt / 6.0) * (k1 + k4 + 2.0 * (k2 + k3)), 1); // If the particle diverges off to infinity, place it back near the origin if (dot(gl_FragColor.xyz, gl_FragColor.xyz) > 1e6) { gl_FragColor.xyz *= 0.0001; } }`, attributes: { xy: [[-4, -4], [0, 4], [4, -4]] }, uniforms: { dt: regl.prop('dt'), t: regl.prop('t'), srcTex: regl.prop('src'), // Convert the integer column index into a texture coordinate in [0, 1], // keeping in mind that texture coordinates are centered a the *center* of each texel srcTexCoordU: (ctx, props) => (props.srcColumn + 0.5) / stepCount }, framebuffer: regl.prop('dst'), count: 3 })36:T1849,{ const w = Math.min(600, width); const h = 150; const svg = d3.select(DOM.svg(w, h)); const colors = ['#cc2255', '#2255dd']; const xScale = d3 .scaleLinear() .domain([-2.5, 2.5]) .range([0, w]); const yScale = d3 .scaleLinear() .domain([-0.5, 0.9]) .range([h, 0]); const line = d3 .line(d => xScale(d[0]), d => yScale(d[1])) .defined(d => !!d && !isNaN(d[0])); const pointsPerSeg = 5; const points = [ [-2, -0.2], [-1.5, 0.5], [-1, 0.2], [-0.5, 0.2], [0, 0.6], null, [0, -0.2], [0.5, -0.1], [1, 0.2], [1.5, 0.3], [2, 0.1] ]; svg .append('path') .attr('stroke', colors[0]) .attr('fill', 'none') .attr('d', line(points.slice(0, pointsPerSeg))); svg .append('path') .attr('stroke', colors[1]) .attr('fill', 'none') .attr('d', line(points.slice(pointsPerSeg + 1))); svg .append('g') .selectAll('circle') .data(points.filter(d => d)) .enter() .append('circle') .attr('cx', d => xScale(d[0])) .attr('cy', d => yScale(d[1])) .attr('r', 2.5) .attr('fill', (d, i) => colors[Math.floor(i / pointsPerSeg)]); const entry = svg .append('g') .selectAll('text') .data(points.filter(d => d)) .enter() .append('text') .attr('x', d => xScale(d[0])) .attr('y', d => yScale(d[1])) .attr('dy', d => 7) .attr('dy', 18) .attr('dominant-baseline', d => 'auto') .attr('font-family', 'KaTeX_Main,Times New Roman,serif') .attr('text-anchor', d => 'start') .attr('fill', (d, i) => colors[Math.floor(i / pointsPerSeg)]); entry .append('tspan') .text('p') .attr('font-weight', 700); entry.append('tspan').text('('); entry .append('tspan') .text('u') .attr('font-weight', 700) .attr('dx', '-0.1em'); entry .append('tspan') .attr('font-size', '.6em') .attr('dy', '-.8em') .text((d, i) => `(${Math.floor(i / pointsPerSeg)})`); entry .append('tspan') .attr('font-size', '.6em') .attr('dy', '1.2em') .attr('dx', '-1.1em') .text((d, i) => i % pointsPerSeg); entry .append('tspan') .attr('dx', '0.35em') .attr('dy', '-0.25em') .text(')'); const particles = 2; const steps = 5; const texture = html`

←v

texture coord u→
${[...Array(particles).keys()].map( (particle, i) => { return html`${[ ...Array(steps).keys() ].map(step => { let eqn; const linenum = Math.floor(step / (pointsPerSeg + 1)); eqn = tex`\mathbf{u}_${step}^{(${particle})}`; return html``; })}`; } )}

${eqn}

vertex

${tex`\textsf{NaN}`}

${tex`\mathbf{u}_{${i}}^{(${p})}`}

endpoint	${tex`\mathbf{u}_{${endpoint}}^{(${line})}`}
isstart	${tex`${isStart}`}

`; const css = html``; return html`

${lookupTable}

Texture lookup vertex attribute buffers

${texture}

State texture coordinates

${svg.node()}

Line geometry

${css}

Data layout for line rendering. We start with vertex attributes containing texture coordinates and sample the state texture in the line drawing vertex shader to position line vertices.

`; }37:T5ec,md`Drawing lines proceeds in two steps. First, we draw interior segments. For this step, we pass the *vertex* attribute above. Each contiguous particle track is separated by a sentinel value. (Here I've called that \`NaN\`, though GLSL is not particularly IEEE754 compliant, so what I've *acutally* done is to pass it the out-of-bounds texture coordinate ${tex`(-1, -1)`} which is then converted into a \`vec4\` position with \`w = 0.0\`, which the function interprets as a line break.) Rendering is accomplished by passing the same attribute aliased *four* times with a different offset. At each vertex we draw a segment and one adjacent join, sliding a window along the line and drawing a separate geometry instance for each segment. [Rye Terrell](https://wwwtyro.net/2021/10/01/instanced-lines-part-2.html) has an excellent post which walks through this process very clearly. The second step is to draw end caps. For this step, we have to be a bit more clever and pass blocks of the three points nearest each end cap. Since GLSL ES 1.00 lacks the \`gl_InstanceID\` keyword, we're a bit limited in our options and need to pass a second per-endpoint-instance attribute which defines whether it's a starting or ending cap. The demo below illustrates the resulting geometry, importing its implementation from the [\`regl\` Instanced Lines](https://observablehq.com/d/a5f8f4c894c77aaf) notebook. Line segments from the first pass are illustrated in blue, while line caps from the second pass are illustrated in red.`38:T670,drawLines = reglLines(regl, { vert: ` precision highp float; #pragma lines: attribute vec2 lookupCoord; #pragma lines: attribute float orientation; #pragma lines: position = projectPoint(lookupCoord); #pragma lines: width = pointWidth(lookupCoord); #pragma lines: varying vec3 color = getColor(lookupCoord); #pragma lines: orientation = getCapOrientation(orientation); uniform mat4 uProjectionView; uniform sampler2D src; uniform float scale, texOffset, width, pixelRatio, particleCount; vec4 projectPoint (vec2 point) { if (point.x < 0.0) return vec4(0); vec3 xyz = texture2D(src, fract(point.xy + vec2(texOffset, 0))).xyz; vec4 p = uProjectionView * vec4(xyz, 1); float w = p.w; p /= w; p.z -= 0.0002 * abs(point.y - 0.5); return p * w; } float getCapOrientation(float orientation) { return orientation; } float pointWidth (vec2 point) { return width; } vec3 getColor (vec2 point) { const vec3 color1 = vec3(0.55, 0.89, 0.65); const vec3 color2 = vec3(0.11, 0.32, 0.65); return mix(color1, color2, point.y); }`, frag: ` precision highp float; varying vec3 lineCoord; varying vec3 color; void main () { const vec3 borderColor = vec3(1); gl_FragColor = vec4(mix( color, borderColor, smoothstep(0.5, 0.9, length(lineCoord.xy)) ), 1); }`, uniforms: { width: (ctx, props) => ctx.pixelRatio * 5, src: regl.prop('src'), texOffset: regl.prop('texOffset') }, depth: { enable: true }, cull: { enable: false } })39:T467,mainDrawLoop = { initializeState; camera.taint(); const simulate = ~opts.indexOf('Simulate'); let frame = regl.frame(() => { try { if (simulate) { integrate({ src: state, dst: tmpStateColumn, srcColumn: mutable currentColumn, dt, t: mutable t }); mutable t += dt; mutable currentColumn = (mutable currentColumn + 1) % stepCount; copyStateColumn({ src: tmpStateColumn, dst: state, dstColumn: mutable currentColumn }); } camera(({ dirty }) => { // Only render if we're currently simulating or if the camera is moved if (!simulate && !dirty) return; regl.clear({ color: [1, 1, 1, 1], depth: 1 }); drawLines({ ...lineData, src: state, texOffset: (mutable currentColumn - stepCount + 1) / stepCount }); }); } catch (e) { console.error(e); frame && frame.cancel(); frame = null; } }); invalidation.then(() => { frame && frame.cancel(); frame = null; }); }3a:[] 2f:{"title":"Strange Attractors on the GPU, Part 1: Implementation","cells":[{"id":880,"value":"md`# Strange Attractors on the GPU, Part 1: Implementation`\n","pinned":false,"hidden":false,"mode":"ojs"},{"id":882,"value":"*[Read this notebook ➡️](https://rreusser.github.io/notebooks/strange-attractors-in-webgpu/) on my personal website via [Observable Notebook Kit](https://observablehq.com/notebook-kit/kit).*","pinned":false,"hidden":false,"mode":"md"},{"id":453,"value":"$31","pinned":false,"hidden":false,"mode":"ojs"},{"id":697,"value":"viewof regl = reglCanvas(this, {\n width: width,\n height: Math.max(400, width * 0.6),\n pixelRatio,\n attributes: { antialias: ~contextOpts.indexOf('Antialiasing') },\n extensions: ['ANGLE_instanced_arrays', 'OES_texture_float']\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":215,"value":"viewof restart = Inputs.button('Restart')","pinned":false,"hidden":false,"mode":"ojs"},{"id":199,"value":"viewof opts = Inputs.checkbox(['Simulate'], {\n value: ['Simulate']\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":40,"value":"viewof particleCount = Inputs.range([1, 4096], {\n value: 20,\n label: 'Particle count',\n step: 1\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":41,"value":"viewof stepCount = Inputs.range([1, 1024], {\n label: 'Track length',\n value: 100,\n transform: Math.log,\n step: 1\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":174,"value":"viewof dt = Inputs.range([0.001, 0.1], {\n value: 0.02,\n label: html`Time step, ${tex`\\Delta t`}`\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":590,"value":"viewof contextOpts = Inputs.checkbox(['Antialiasing'], {\n value: ['Antialiasing']\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":349,"value":"section_problemSetup = md`## Simulation method\n\nTo simulate strange attractors, here we use the [regl](https://github.com/regl-project/regl) WebGL library. I like \\`regl\\` because it's mostly just a convenience wrapper around the low-level WebGL API, adding few of its own abstractions while making WebGL fun and easy to use—especially for creative simulations like this.\n\nHowever, there's a big downside. When you rely on low-level graphics APIs, even simple things can become very difficult. Lines are a good example. All browsers these days limit the built-in [OpenGL line primitive](https://www.khronos.org/opengl/wiki/Primitive#Line_primitives) to a single pixel width, so if you want any reasonably well-rendered lines, you really just have to implement it yourself.\n\nThus we have two challenges: to simulate the equations and to draw lines.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":363,"value":"md`## WebGL context creation\n\nTo get started, we import a helper for instantiating the [\\`regl\\`](https://github.com/regl-project/regl) library and create a context in the [\\`regl\\`](#regl) cell above (i.e. the main visualization).`","pinned":false,"hidden":false,"mode":"ojs"},{"id":688,"value":"import { reglCanvas } from '@rreusser/regl-canvas'","pinned":false,"hidden":false,"mode":"ojs"},{"id":4,"value":"pixelRatio = window.devicePixelRatio","pinned":false,"hidden":false,"mode":"ojs"},{"id":353,"value":"md`The reason for this helper is that browsers limit the total number of WebGL contexts before the oldest start getting removed—usually 8 or 16 contexts. If we're not careful, dependence on the notebook's \\`width\\` value will trigger instantiaton of tens of new contexts when we resize the page, causing figures to start to go blank. To avoid this, the \\`reglCanvas\\` helper prevents recreating the context when only the size has changed.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":472,"value":"md`Next, we instantiate a camera. I wish I had a better camera module to offer, but the only module I have on hand is the [inertial-turntable-camera](https://github.com/standardCyborg/inertial-turntable-camera) module I wrote some time ago, and for which I have a lot of regrets. It needs a major overhaul which I haven't yet attacked. Really though, any camera model which exposes a standard, combined projection-view matrix would do just fine.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":18,"value":"import { createReglCamera, createInteractions } from '@rreusser/regl-tools'","pinned":false,"hidden":false,"mode":"ojs"},{"id":22,"value":"camera = {\n restart;\n const camera = createReglCamera(regl, {\n distance: 7,\n near: 0.01,\n far: 100,\n phi: 0,\n theta: -1,\n center: [0, 2, 0]\n });\n createInteractions(camera);\n return camera;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":726,"value":"updateCameraOnResize = {\n // Trigger a dependence of resizing on the page size and regl instance:\n width, regl;\n\n const aspectRatio = viewof regl.width / viewof regl.height;\n camera.resize(aspectRatio);\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":356,"value":"section_dataModel = md`## State layout\nNext, we define how we represent particle track data in a WebGL texture. The state of our ordinary differential equation (ODE) is represented by the three-component vector ${tex`(x, y, z)`}. WebGL generally gives us four color channels to work with, so we choose to represent the state in a four-channel floating point texture. The ${tex`j^{th}`} time step of the ${tex`i^{th}`} particle is thus represented by the four-component vector, ${tex.block`\\mathbf{p}_j^{(i)} = (x_j^{(i)}, y_j^{(i)}, z_j^{(i)}, 1).`}\n\nWe pack these vectors into a texture, each row representing the history of a single particle. `","pinned":false,"hidden":false,"mode":"ojs"},{"id":332,"value":"$32","pinned":false,"hidden":false,"mode":"ojs"},{"id":370,"value":"md`As we step the ODE, we compute one new history point for each particle track. To avoid having to shift the entire history over one texel on every iteration just to append a new step, we instead treat each row as a ring buffer. At each time step ${tex`j`}, we use the previous column of the history, ${tex`p_{j-1}^{(i)}`}, to compute the next time step, ${tex`p_{j}^{(i)}`}. When we reach the end of the row, we loop back to the start, at each step overwriting the oldest time step with the newest.\n\nThe cells below define empty textures which will hold our state data.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":38,"value":"state = {\n const fbo = regl.framebuffer({\n depthStencil: false,\n width: stepCount,\n height: particleCount,\n colorType: 'float',\n colorFormat: 'rgba'\n });\n\n // Clean up resources when Observable reevaluates this cell\n invalidation.then(() => fbo.destroy());\n\n return fbo;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":66,"value":"tmpStateColumn = {\n const fbo = regl.framebuffer({\n depthStencil: false,\n width: 1,\n height: particleCount,\n colorType: 'float',\n colorFormat: 'rgba'\n });\n invalidation.then(() => fbo.destroy());\n return fbo;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":84,"value":"mutable currentColumn = 0","pinned":false,"hidden":false,"mode":"ojs"},{"id":149,"value":"mutable t = 0","pinned":false,"hidden":false,"mode":"ojs"},{"id":286,"value":"md`## Computation\n\nWe now define our WebGL drawing commands. Our fundamental computational primitive will be similar to a *map* operation. WebGL 1 can't modify textures in-place, so this operation will read from one texture (or zero or two textures, if we like), perform some computation in a fragment shader, and write to an output texture.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":329,"value":"$33","pinned":false,"hidden":false,"mode":"ojs"},{"id":378,"value":"md`Modern graphics APIs more directly treat data as data, but in WebGL 1 we must deal with data as if it's colors, triangles, and textures. So in order to step the ODE and update the state, we therefore need to draw some sort of geometric primitive and assign each state texel a color. Many people use a fullscreen quad composed of two triangles for this operation, but I find it simpler (and [potentially better for performance](https://www.npmjs.com/package/a-big-triangle)) to use a single big triangle which covers the entire viewport.\n\nWhichever you prefer, when we write such commands, in a sense we need to think backwards: we focus on the *output* texture we're evaluating and then work backward, sampling our input textures as needed and executing a fragment shader to determine the resulting output color.\n\nThus, the *vertex shader* will trivially position this triangle, and the *fragment shader* is where we will perform the real computation.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":381,"value":"section_initialization = md`### Initialization\n\nWe start with a simple draw command. Initialization uniformly assigns each row of the state texture (all of a particle's history) to some random location within a sphere. Generating good pseudorandom numbers on a GPU is tricky, so to initialize our state, we use a low-discrepancy quasirandom number generator described by Martin Roberts in *[The Unreasonable Effectiveness\nof Quasirandom Sequences](http://extremelearning.com.au/unreasonable-effectiveness-of-quasirandom-sequences/)*.\n\nThe command below implements this step and writes the result to an output framebuffer.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":45,"value":"$34","pinned":false,"hidden":false,"mode":"ojs"},{"id":47,"value":"triggerInitialization = {\n restart;\n regl.poll();\n mutable currentColumn = 0;\n mutable t = 0;\n initializeState({ dst: state });\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":386,"value":"section_integration = md`### Integration\nNext, we define a shader to integrate the ODE. Our output destination is a texture just large enough to hold the single time step we are computing. Its dimension is therefore 1 texel wide and \\`particleCount\\` texels tall. In the fragment shader, we sample the current state from our full state texture, perform an integration step, and output the updated state as the fragment color. Finally, we invoke an additional \\`copyStateColumn\\` draw command to transfer this single history slice back into our main state texture.\n\nFor integration, we use the [fourth-order Runge-Kutta](https://en.wikipedia.org/wiki/Runge%E2%80%93Kutta_methods#The_Runge%E2%80%93Kutta_method) (RK4) method and directly substitute a GLSL derivative function into the shader.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":144,"value":"attractorGLSL = `\nvec3 derivative (float x, float y, float z, float t) {\n const float alpha = 3.0;\n const float beta = 2.20;\n const float gamma = 1.0;\n const float mu = 1.510;\n return vec3(\n alpha * x * (1.0 - y) - beta * z,\n -gamma * y * (1.0 - x * x),\n mu * x\n );\n}`","pinned":false,"hidden":false,"mode":"ojs"},{"id":69,"value":"$35","pinned":false,"hidden":false,"mode":"ojs"},{"id":389,"value":"md`After running this command, we will have computed our updated state, except it will be in the the wrong place—our temporary texture. The following command implements the copy operation so that we can transfer it back to the correct place in our state texture. It uses scissoring to limit this copy to a single column, and we use a mutable counter to track the current time step column.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":80,"value":"copyStateColumn = regl({\n vert: `\n precision mediump float;\n attribute vec2 xy;\n varying vec2 uv;\n void main () {\n uv = 0.5 * xy + 0.5;\n gl_Position = vec4(xy, 0, 1);\n }`,\n frag: `\n precision mediump float;\n varying vec2 uv;\n uniform sampler2D srcTex;\n void main () {\n gl_FragColor = texture2D(srcTex, uv);\n }`,\n attributes: { xy: [[-4, -4], [0, 4], [4, -4]] },\n uniforms: {\n srcTex: regl.prop('src')\n },\n framebuffer: regl.prop('dst'),\n scissor: {\n enable: true,\n box: {\n x: regl.prop('dstColumn'),\n y: 0,\n height: particleCount,\n width: 1\n }\n },\n count: 3\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":340,"value":"md`## Line Rendering\n\nIn the timeless words of Matt DesLauriers, *[Drawing Lines is Hard](https://mattdesl.svbtle.com/drawing-lines-is-hard)*. This notebook didn't really come together until [Rye Terrell](https://twitter.com/wwwtyro)'s post, *[Instanced Line Rendering Part II: Alpha blending](https://wwwtyro.net/2021/10/01/instanced-lines-part-2.html)*, spurred me to resurrect my long-running attempts to come up with a clean solution for line rendering.\n\nThe resulting module is published as [regl-gpu-lines](https://github.com/rreusser/regl-gpu-lines). At the end of the day, it's really nothing more than a basic line rendering function, but it's written to be as general possible so that you can connect \\`attribute\\` and \\`varying\\` inputs to your data in the particular way which is meaningful to you.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":167,"value":"reglLines = require('https://unpkg.com/regl-gpu-lines@2.2.0')","pinned":true,"hidden":false,"mode":"ojs"},{"id":423,"value":"md`The general flow of vertex data required to render lines is illustrated below. Instead of particle positions in the vertex attributes, we instead store *texture coordinates*. Then in the vertex shader, we perform a texture lookup to retrieve the corresponding state vector, which we plot as the vertex position. \n\nThe texture coordinates of the ${tex`i^{th}`} particle at the ${tex`j^{th}`} step are ${tex`(u_j^{(i)}, v_j^{(i)})`}, denoted ${tex`\\mathbf{u}_j^{(i)}`}. They are equal to\n\n${tex.block`\\mathbf{u}_j^{(i)} = \\left(\\frac{j + \\frac{1}{2}}{M}, \\frac{i + \\frac{1}{2}}{N}\\right),`}\nwhere there are ${tex`M`} total time steps and ${tex`N`} total particle tracks. Texture sampler coordinates are located at the center of each texel, accounting for the half-texel offset.\n\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":323,"value":"$36","pinned":false,"hidden":false,"mode":"ojs"},{"id":426,"value":"$37","pinned":false,"hidden":false,"mode":"ojs"},{"id":747,"value":"lineExampleRegl._gl.canvas","pinned":false,"hidden":false,"mode":"ojs"},{"id":437,"value":"viewof lineParams = multiplexInputs(\n {\n join: Inputs.select(['bevel', 'miter', 'round'], {\n value: 'miter',\n label: 'Joins'\n }),\n cap: Inputs.select(['none', 'square', 'round'], {\n value: 'round',\n label: 'Caps'\n }),\n capResolution: Inputs.range([1, 10], {\n value: 4,\n label: 'Cap resolution',\n step: 1\n //disabled: cap !== 'round'\n }),\n joinResolution: Inputs.range([1, 10], {\n value: 3,\n label: 'Join resolution',\n step: 1\n //disabled: join !== 'round'\n }),\n miterLimit: Inputs.range([1, 8], {\n value: 2,\n step: 0.1,\n label: 'Miter limit'\n //disabled: join !== 'miter'\n })\n },\n { label: 'Line parameters', open: true }\n)","pinned":false,"hidden":false,"mode":"ojs"},{"id":759,"value":"renderExample = renderLineExample(lineParams)","pinned":false,"hidden":false,"mode":"ojs"},{"id":469,"value":"md`One small trick is to use the texture coordinate corresponding to the particle ID to offset the [Normalized Device Coordinates](https://learnopengl.com/Getting-started/Coordinate-Systems) z-coordinate of the projected line vertices. At the cost of strictly correct z-ordering, this prevents a bit of the excessive z-fighting which we otherwise encounter.\n\nThe function below defines the line rendering command. The \\`#pragma\\` statements are just something I've made up as a means of defining how it connects to the attribute inputs and converts those to positions and widths.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":170,"value":"$38","pinned":false,"hidden":false,"mode":"ojs"},{"id":465,"value":"md`The lines below construct the vertex attributes containing texture coordinate which the line drawing command uses to sample our state texture.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":90,"value":"vertexAttributes = {\n const positions = [];\n for (let i = 0; i < particleCount; i++) {\n for (let j = 0; j < stepCount; j++) {\n positions.push([(j + 0.5) / stepCount, (i + 0.5) / particleCount]);\n }\n if (i < particleCount - 1) {\n positions.push([-1, -1]);\n }\n }\n\n let lookupCoord = regl.buffer(positions);\n invalidation.then(() => lookupCoord.destroy());\n return { lookupCoord, count: positions.length };\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":92,"value":"endpointAttributes = {\n const endpoints = [];\n const orientation = [];\n let count = 0;\n for (let j = 0; j < particleCount; j++) {\n // Start cap\n for (let i = 0; i < 3; i++) {\n endpoints.push([(i + 0.5) / stepCount, (j + 0.5) / particleCount]);\n }\n\n // End cap\n for (let i = stepCount - 1; i >= stepCount - 3; i--) {\n endpoints.push([(i + 0.5) / stepCount, (j + 0.5) / particleCount]);\n }\n\n // One start, one end\n orientation.push(reglLines.START_CAP, reglLines.END_CAP);\n count += 2;\n }\n\n let endpointBuffer = regl.buffer(endpoints);\n let orientationBuffer = regl.buffer(new Uint8Array(orientation));\n invalidation.then(() => {\n endpointBuffer.destroy();\n orientationBuffer.destroy();\n });\n return { lookupCoord: endpointBuffer, orientation: orientationBuffer, count };\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":100,"value":"lineData = ({\n cap: 'round',\n join: 'bevel',\n capResolution: 4,\n vertexAttributes,\n endpointAttributes,\n vertexCount: vertexAttributes.count,\n endpointCount: endpointAttributes.count\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":834,"value":"md`## Draw loop`","pinned":false,"hidden":false,"mode":"ojs"},{"id":396,"value":"md`Finally we put it all together into a main drawing loop which performs a full step of the ODE and then configures the camera and renders lines.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":189,"value":"$39","pinned":true,"hidden":false,"mode":"ojs"},{"id":475,"value":"md`And that's it! You can see the result at the top of this page. Stay tuned for [Strange Attractors on the GPU, Part 2 🔗](https://observablehq.com/d/df2e1a97f7ab9a46), in which we copy/paste the exact same code, except focusing on making it look as good as possible rather than clear exposition.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":415,"value":"md`## Imports`","pinned":false,"hidden":false,"mode":"ojs"},{"id":439,"value":"import { multiplexInputs } from '@rreusser/multiplex-inputs'","pinned":false,"hidden":false,"mode":"ojs"},{"id":435,"value":"import {\n regl2 as lineExampleRegl,\n renderExample as renderLineExample\n} from 'a5f8f4c894c77aaf'","pinned":false,"hidden":false,"mode":"ojs"}],"files":"$W3a","stdlib":"1","titling":"auto"}