31:T4c3,# Circle inversion fractals Here we provide code for rendering [Circle Inversion Fractals](https://users.math.yale.edu/public_html/People/frame/Fractals/CircInvFrac/welcome.html). In fact, we focus on fractals generated as [_restricted limit sets_](https://users.math.yale.edu/public_html/People/frame/Fractals/CircInvFrac/RestrLimSet/RestrLimSet.html), since these generate more interesting pictures. We present two implementations, one easier to understand, using canvas, and another, which runs faster but is harder to understand, using GPU programming with the help of the [GPU-IO library](https://github.com/amandaghassaei/gpu-io). Both have controls for selecting two circle datasets which are interpolated linearly using a third control. This allows for some interesting combinations. _**New:**_ A _shuffle_ option has been added, which will shuffle the circle arrays so that the correspondence between the circles in the two sets becomes unpredictable. For interp. values 0 and 1, we still get the same fractals, but for other values we get even more interesting combinations. Notice that the shuffle is determined by a seed, so you can reproduce results. A zero shuffle seed disables shuffling. Have fun!32:T9d4,{ // Init a canvas element. const height = ~~((9 * width) / 16); const canvas = DOM.canvas(width, height, 1); yield visibility(canvas); // Avoid running cell if not visible const scaleFunc = ([x, y, r]) => [ (x * height * generalScale) / 2 + width / 2, (y * height * generalScale) / 2 + height / 2, (r * height * generalScale) / 2 ]; let circlesA = [...circleDatasets[dsA].map(scaleFunc)]; let circlesB = [...circleDatasets[dsB].map(scaleFunc)]; if (shuffleSeed) { let shuff = shuffler(shuffleSeed); if (circlesA.length > circlesB.length) circlesA = shuff(circlesA); else circlesB = shuff(circlesB); } const [smallest, biggest] = circlesA.length > circlesB.length ? [circlesB, circlesA] : [circlesA, circlesB]; const nmin = smallest.length; const n = biggest.length; for (let i = nmin; i < n; i++) { smallest[i] = [...biggest[i].slice(0, 2), missingCircleScale * height]; } const BATCH = 10000; // Inversions const invert = function ([px, py], [cx, cy], r) { let r2 = r * r; let [vx, vy] = [px - cx, py - cy]; let t = r2 / (vx * vx + vy * vy); return [cx + t * vx, cy + t * vy]; }; const ctx = canvas.getContext("2d"); const palette = d3["scheme" + colorPalette].map((c) => { c = d3.color(c); c.opacity = 0.1; return c.formatRgb(); }); let circles = []; let alpha = -1; for (;;) { if (alpha != interp.value) { alpha = interp.value; for (let i = 0; i < n; i++) circles[i] = [0, 1, 2].map( (j) => circlesA[i][j] * (1 - alpha) + circlesB[i][j] * alpha ); ctx.fillStyle = backgroundColor; ctx.fillRect(0, 0, width, height); ctx.strokeStyle = "gray"; if (drawCircles) for (let [x, y, r] of circles) { ctx.beginPath(); ctx.arc(x, y, r, 0, Math.PI * 2); ctx.stroke(); } } for (let j = 0; j < BATCH; j++) { let [px, py] = [Math.random() * width * 20, Math.random() * height * 20]; let icolor; for (let n = 20; n > 0; n--) { let i = Math.floor(Math.random() * circles.length); if (n > 0) icolor = i; let [cx, cy, r] = circles[i]; let [dx, dy] = [px - cx, py - cy]; if (Math.hypot(dx, dy) > r) { [px, py] = invert([px, py], [cx, cy], r); } } ctx.fillStyle = palette[icolor % palette.length]; ctx.fillRect(px, py, 1, 1); ctx.beginPath(); } yield visibility(canvas); } }33:T17a4,fractalRenderGPU = { const { GPUProgram, GPUComposer, GPULayer, setColorProgram, FLOAT, INT } = gpuio; const scaleFunc = (circle) => circle.map((d) => d * generalScale); let circlesA = circleDatasets[circleDatasetNameA].map(scaleFunc); let circlesB = circleDatasets[circleDatasetNameB].map(scaleFunc); if (shuffleSeedGPU) { let shuff = shuffler(shuffleSeedGPU); if (circlesA.length > circlesB.length) circlesA = shuff(circlesA); else circlesB = shuff(circlesB); } const [smallest, biggest] = circlesA.length > circlesB.length ? [circlesB, circlesA] : [circlesA, circlesB]; const nmin = smallest.length; const n = biggest.length; for (let i = nmin; i < n; i++) { smallest[i] = [ ...biggest[i].slice(0, 2), generalScale * missingCircleScale ]; } let circles = circlesA.concat(circlesB); const circlesShaderCode = function (circles) { const spec = circles.map(([x, y, r]) => `vec3(${x},${y},${r})`).join(",\n"); return ` const int n = ${n}; const vec3 circles[${n * 2}] = vec3[${n * 2}](${spec}); uniform float u_t; vec3 getCircle (int i) { return mix (circles[i], circles[i+n], u_t); } `; }; // Init a canvas element. const height = ~~((9 * width) / 16); const canvas = DOM.canvas(width, height, 1); yield visibility(canvas); // Avoid running cell if not visible // Composer const composer = new GPUComposer({ canvas, contextAttributes: { preserveDrawingBuffer: true, premultipliedAlpha: false, antialias: true, alpha: false } }); // Circle renderer const circleRenderProgram = new GPUProgram(composer, { name: "circleRender", fragmentShader: ` in vec2 v_uv; out vec4 out_color; uniform float u_aspect; ${circlesShaderCode(circles)} void main() { float val = 0.0; vec2 coords = v_uv*vec2(2.*u_aspect, 2.)-vec2(u_aspect,1.); for (int i = 0; i < n; i++) { vec3 circle = getCircle (i); vec2 d = coords - circle.xy; float diff = abs(length(d) - circle.z); if (diff < 0.004) val = smoothstep(0.004,0.0,diff); } if (val == 0.) discard; out_color = vec4(val); } `, uniforms: [ { // resolution of the canvas. name: "u_aspect", value: composer._width / composer._height, type: FLOAT }, { // interpolation parameter name: "u_t", value: 0, type: FLOAT } ] }); // Inversion shader circles = circles.map(([x, y, r]) => [ (x * height) / 2 + width / 2, (y * height) / 2 + height / 2, (r * height) / 2 ]); const inversionProgram = new GPUProgram(composer, { name: "inversionProgram", fragmentShader: ` in vec2 v_uv; uniform sampler2D u_positions; out vec4 out_position; ${circlesShaderCode(circles)} vec2 invert (vec2 p, vec3 circle) { float r2 = circle.z*circle.z; vec2 v = p - circle.xy; float t = r2 / dot(v,v); return circle.xy + t * v; } float random (vec2 st) { return fract(sin(dot(st,vec2(12.9898,78.233)))*43758.5453123); } void main() { vec2 coords = texture(u_positions, v_uv).xy; vec2 seed = coords; float clan = 0.; for (int j = 0; j < 40; j++) { int i = int (random(seed)*float(n)); seed = seed+vec2(0.98183,0.2841841); vec3 circle = getCircle(i); if (j == 39) clan = float (i); vec2 d = coords - circle.xy; if (length(d) > circle.z) coords = invert(coords, circle); } clan = clan / float(n); out_position = vec4(coords,clan,clan); } `, uniforms: [ { // The positions array name: "u_positions", value: 0, type: INT }, { // interpolation parameter name: "u_t", value: 0, type: FLOAT } ] }); // // Init particles // const NUM_PARTICLES = 200000; // Init new positions. const positions = new Float32Array(NUM_PARTICLES * 4); // Init particle state const particleState = new GPULayer(composer, { name: "position", dimensions: NUM_PARTICLES, type: FLOAT, numComponents: 4, numBuffers: 2 }); // Random particles init function reinitParticles() { for (let i = 0; i < positions.length; i += 4) { positions[i] = Math.random() * width * 10; positions[i + 1] = Math.random() * height * 10; } particleState.resize(NUM_PARTICLES, positions); } // Render particles program const renderParticles = new GPUProgram(composer, { name: "renderParticles", fragmentShader: ` in vec2 v_uv; in vec2 v_uv_position; uniform sampler2D u_position; out vec4 color; ${glslPalette} void main() { int iColor = int(texture(u_position, v_uv_position).z * float(nPal)); color = vec4(palette[iColor], 0.1); // Premultiply alpha! }` }); // Init a program to draw color. const colorProgram = setColorProgram(composer, { name: "setColor", type: FLOAT }); colorProgram.setUniform("u_color", glslColor(backgroundColor)); let t = -1; let clearFlag = true; for (;;) { if (t != testInterpolationParam.value) { t = testInterpolationParam.value; circleRenderProgram.setUniform("u_t", t); inversionProgram.setUniform("u_t", t); clearFlag = true; } reinitParticles(); composer.step({ program: inversionProgram, input: particleState, output: particleState }); if (clearFlag) { composer.step({ program: colorProgram }); if (showCircles) composer.step({ program: circleRenderProgram }); } composer.drawLayerAsPoints({ layer: particleState, program: renderParticles, input: particleState, blendAlpha: true }); yield visibility(canvas); clearFlag = false; } }34:Tc45,circleDatasets = ({ "3 circles A": [ [0, -Math.sqrt(3) / 4, 0.5], [-0.5, Math.sqrt(3) / 4, 0.5], [0.5, Math.sqrt(3) / 4, 0.5] ], "4 circles A": [ [0.5, 0.5, 0.5], [-0.5, 0.5, 0.5], [0.5, -0.5, 0.5], [-0.5, -0.5, 0.5] ], "4 circles B": [ [0.35, 0.35, 0.4], [-0.35, 0.35, 0.4], [0.35, -0.35, 0.4], [-0.35, -0.35, 0.4] ], "5 circles A": [ [0.5, 0.5, 0.5], [-0.5, 0.5, 0.5], [0.5, -0.5, 0.5], [-0.5, -0.5, 0.5], [0, 0, 0.75 / 2] ], "5 circles B": [ [0.5, 0.5, 0.5], [-0.5, 0.5, 0.5], [0.5, -0.5, 0.5], [-0.5, -0.5, 0.5], [0, 0, 0.2] ], "5 circles C": [ [0, -Math.sqrt(3) / 4, 0.5], [-1, -Math.sqrt(3) / 4, 0.5], [1, -Math.sqrt(3) / 4, 0.5], [-0.5, Math.sqrt(3) / 4, 0.5], [0.5, Math.sqrt(3) / 4, 0.5] ], "5 circles D": [ [0, -Math.sqrt(3) / 4, 0.5], [-0.5, -Math.sqrt(3) / 4, 0.5], [0.5, -Math.sqrt(3) / 4, 0.5], [-0.5, Math.sqrt(3) / 4, 0.5], [0.5, Math.sqrt(3) / 4, 0.5] ], "6 circles A": [ [-0.25, -Math.sqrt(3) / 4, 0.5], [0.25, -Math.sqrt(3) / 4, 0.5], [-0.25, Math.sqrt(3) / 4, 0.5], [0.25, Math.sqrt(3) / 4, 0.5], [-0.75, 0, 0.5], [0.75, 0, 0.5] ], "6 circles B": [ [-0.25, -Math.sqrt(3) / 4, 0.25], [0.25, -Math.sqrt(3) / 4, 0.25], [-0.25, Math.sqrt(3) / 4, 0.25], [0.25, Math.sqrt(3) / 4, 0.25], [-0.5, 0, 0.25], [0.5, 0, 0.25] ], "7 circles A": [ [-0.25, -Math.sqrt(3) / 4, 0.25], [0.25, -Math.sqrt(3) / 4, 0.25], [-0.25, Math.sqrt(3) / 4, 0.25], [0.25, Math.sqrt(3) / 4, 0.25], [-0.5, 0, 0.25], [0.5, 0, 0.25], [0, 0, 0.25] ], "7 circles B": [ [-0.25, -Math.sqrt(3) / 4, 0.5], [0.25, -Math.sqrt(3) / 4, 0.5], [-0.25, Math.sqrt(3) / 4, 0.5], [0.25, Math.sqrt(3) / 4, 0.5], [-0.5, 0, 0.5], [0.5, 0, 0.5], [0, 0, 0.5] ], "7 circles C": [ [-0.25, -Math.sqrt(3) / 4, 0.5], [0.25, -Math.sqrt(3) / 4, 0.5], [-0.25, Math.sqrt(3) / 4, 0.5], [0.25, Math.sqrt(3) / 4, 0.5], [-0.75, 0, 0.5], [0.75, 0, 0.5], [0, 0, 0.5] ], "7 circles D": [ [0, -Math.sqrt(3) / 4, 0.5], [-1, -Math.sqrt(3) / 4, 0.5], [1, -Math.sqrt(3) / 4, 0.5], [-0.5, Math.sqrt(3) / 4, 0.5], [0.5, Math.sqrt(3) / 4, 0.5], [-0.25, 0, 0.3], [0.25, 0, 0.3] ], "8 circles A": [ [0.5, 0.5, 0.5], [-0.5, 0.5, 0.5], [0.5, -0.5, 0.5], [-0.5, -0.5, 0.5], [0.25, 0, 0.2], [0, 0.25, 0.2], [0, -0.25, 0.2], [-0.25, 0, 0.2] ], "8 circles B": [ [0.5, 0.5, 0.5], [-0.5, 0.5, 0.5], [0.5, -0.5, 0.5], [-0.5, -0.5, 0.5], [0.5, 0, 0.5], [0, 0.5, 0.5], [0, -0.5, 0.5], [-0.5, 0, 0.5] ], "8 circles C": [ [0.25, 0, 0.1], [0, 0.25, 0.1], [0, -0.25, 0.1], [-0.25, 0, 0.1], [0.5, 0.5, 0.5], [-0.5, 0.5, 0.5], [0.5, -0.5, 0.5], [-0.5, -0.5, 0.5] ], "9 circles A": [ [0.5, 0.5, 0.5], [-0.5, 0.5, 0.5], [0.5, -0.5, 0.5], [-0.5, -0.5, 0.5], [0.25, 0, 0.2], [0, 0.25, 0.2], [0, -0.25, 0.2], [-0.25, 0, 0.2], [0, 0, 0.2] ] })35:[] 2f:{"title":"Circle inversion fractals","cells":[{"id":0,"value":"$31","pinned":false,"hidden":false,"mode":"md"},{"id":509,"value":"
\n## General settings\n\nThese affect both the CPU-only and the GPU-based implementation.","pinned":false,"hidden":false,"mode":"md"},{"id":513,"value":"//\n// Scale used to zoom the fractal\n//\nviewof generalScale = Inputs.range([0.2, 8], {\n label: \"scale\",\n value: 2,\n step: 0.1\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":511,"value":"//\n// If one dataset used in the interpolation has fewer circles than the other, missing circles\n// are generated with radii multiplied by this factor. The idea is that huge circles have a small effect on the\n// resulting fractal.\n//\nviewof missingCircleScale = Inputs.range([2, 100], {\n label: \"missing circle scale\",\n value: 5,\n step: 1\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":594,"value":"viewof backgroundColor = Inputs.color({\n label: \"background color\",\n value: \"#000000\"\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":555,"value":"viewof colorPalette = Inputs.select(\n [\n \"Accent\",\n \"Tableau10\",\n \"Category10\",\n \"Dark2\",\n \"Paired\",\n \"Pastel1\",\n \"Pastel2\",\n \"Set1\",\n \"Set2\",\n \"Set3\"\n ],\n { label: \"drawing palette\", value: \"Accent\" }\n)","pinned":false,"hidden":false,"mode":"ojs"},{"id":571,"value":"swatches(colorPalette)","pinned":false,"hidden":false,"mode":"ojs"},{"id":518,"value":"
\n## CPU implementation (slower)","pinned":false,"hidden":false,"mode":"md"},{"id":529,"value":"viewof drawCircles = Inputs.toggle({ label: \"display circles\", value: true })","pinned":false,"hidden":false,"mode":"ojs"},{"id":609,"value":"viewof shuffleSeed = Inputs.range([0, 1000], {\n label: html`shuffle seed
(0 = no shuffle)`,\n step: 1,\n value: 0\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":298,"value":"viewof dsA = Inputs.select([...Object.keys(circleDatasets)], {\n label: \"Circle Dataset A\",\n value: \"7 circles A\"\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":478,"value":"viewof dsB = Inputs.select([...Object.keys(circleDatasets)], {\n label: \"Circle Dataset B\",\n value: \"8 circles C\"\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":486,"value":"interp = Inputs.range([0, 1], {\n label: \"Interp. value\",\n step: 0.01,\n value: 0\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":93,"value":"$32","pinned":false,"hidden":false,"mode":"ojs"},{"id":583,"value":"
\n## GPU implementation\n\nUses the [GPU-IO library](https://github.com/amandaghassaei/gpu-io). See also [this notebook](https://observablehq.com/@esperanc/gpu-io).","pinned":false,"hidden":false,"mode":"md"},{"id":640,"value":"viewof shuffleSeedGPU = Inputs.range([0, 1000], {\n label: html`shuffle seed
(0 = no shuffle)`,\n step: 1,\n value: 0\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":287,"value":"viewof showCircles = Inputs.toggle({ label: \"Show Circles\", value: true })","pinned":false,"hidden":false,"mode":"ojs"},{"id":257,"value":"viewof circleDatasetNameA = Inputs.select([...Object.keys(circleDatasets)], {\n label: \"Select Dataset A\",\n value: \"5 circles A\",\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":259,"value":"viewof circleDatasetNameB = Inputs.select([...Object.keys(circleDatasets)], {\n label: \"Select Dataset B\",\n value: \"8 circles A\"\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":230,"value":"testInterpolationParam = Inputs.range([0, 1], {\n label: \"Interp. value\",\n step: 0.01,\n value: 0\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":189,"value":"$33","pinned":false,"hidden":false,"mode":"ojs"},{"id":600,"value":"glslColor = (cssColor) => {\n let c = d3.color(cssColor);\n return [c.r, c.g, c.b].map((x) => x / 255);\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":417,"value":"glslPalette = {\n const scheme = d3[\"scheme\" + colorPalette];\n const n = scheme.length;\n let glslColors = scheme.map((c) => {\n let a = glslColor(c);\n return `vec3(${a[0]},${a[1]},${a[2]})`;\n });\n return `\n const int nPal = ${n};\n const vec3 palette[nPal]= vec3[nPal](${glslColors.join(\",\")});\n `;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":626,"value":"//\n// A shuffle function that gives predictable results depending on the seed\n//\nshuffler = (seed) => d3.shuffler(d3.randomLcg(seed))","pinned":true,"hidden":false,"mode":"ojs"},{"id":550,"value":"## Circle datasets","pinned":false,"hidden":false,"mode":"md"},{"id":125,"value":"$34","pinned":true,"hidden":false,"mode":"ojs"},{"id":576,"value":"## Dependencies","pinned":false,"hidden":false,"mode":"md"},{"id":568,"value":"import { swatches } from \"@d3/color-schemes\"","pinned":true,"hidden":false,"mode":"ojs"},{"id":5,"value":"gpuio = require(\"gpu-io\")","pinned":true,"hidden":false,"mode":"ojs"}],"files":"$W35","stdlib":"1","titling":"auto"}