31:T6e5,function drawSpheres(spheres, canvas, options = {}) { const { width, height } = canvas; const { background = "lightgray" } = options; const renderer = new THREE.WebGLRenderer({ canvas, antialias: true }); renderer.setSize(width, height); renderer.shadowMap.enabled = true; renderer.shadowMap.type = THREE.PCFSoftShadowMap; const scene = new THREE.Scene(); scene.background = new THREE.Color(background); const camera = new THREE.PerspectiveCamera(40, width / height, 1, 10000); camera.position.set(0, 2, 3); const controls = new THREE.OrbitControls(camera, renderer.domElement); const light = new THREE.DirectionalLight(0xffffff, 1); light.position.set(10, 10, 10); const ambLight = new THREE.AmbientLight(0xffffff, 0.5); camera.add(light, ambLight); light.castShadow = true; scene.add(camera); const materials = []; const colors = []; const sphereGeom = new THREE.SphereGeometry(1, 32, 16); for (let { center, radius, color } of spheres) { if (radius < 0) continue; color = new THREE.Color(color || "white"); let colorIndex = colors.findIndex((c) => color.equals(c)); let mat; if (colorIndex == -1) { colors.push(color); mat = new THREE.MeshStandardMaterial({ metalness: 0.5, roughness: 0.8, color }); materials.push(mat); } else mat = materials[colorIndex]; const mesh = new THREE.Mesh(sphereGeom, mat); mesh.position.set(...center); mesh.scale.set(radius, radius, radius); mesh.castShadow = true; mesh.receiveShadow = true; scene.add(mesh); } return { draw: function () { renderer.render(scene, camera); }, dispose: function () { renderer.dispose(); controls.dispose(); } }; }32:[{"id":"9e1276d9f0c0e165d4556d641abf0b17162ba23f60c9dc6f2e994d8b9d2b664646d5a1680a1f19207e69f89d965c2ffdc0efbac49fb4ff7443db5d69b00c3a4a","href":"https://static.observableusercontent.com/files/9e1276d9f0c0e165d4556d641abf0b17162ba23f60c9dc6f2e994d8b9d2b664646d5a1680a1f19207e69f89d965c2ffdc0efbac49fb4ff7443db5d69b00c3a4a","name":"Screen Shot 2023-01-23 at 18.08.12.png","type":"image/png","lastModified":1674508125478.792,"size":2357877}] 2f:{"title":"3D Apollonian Sphere Packings","cells":[{"id":0,"value":"# 3D Apollonian Sphere Packings\n\n

\n\nThis is a notebook for working out some of the ideas of the paper [THE FRACTAL DIMENSION OF THE APOLLONIAN SPHERE PACKING](https://doi.org/10.1142/S0218348X94000739) by Martin Borkovec, Wilbur De Paris and Ronald Peikert. The code here derives from this \n[Solution by Kean Walmsley](https://www.keanw.com/2012/02/sphere-packing-in-autocad-creating-an-apollonian-packing-using-f-part-2.html) based on the original code of Ron Peikert.","pinned":false,"hidden":false,"mode":"md"},{"id":8,"value":"## Inversive coordinates\n\nThe first idea is that of [_sphere inversion_](https://en.wikipedia.org/wiki/Inversive_geometry). The authors show how the problem of sphere packing can be elegantly expressed using inverted spheres. Furthermore, they claim that the computation becomes easier using\n_inversive coordinates_ for spheres. In these coordinates, the inversion becomes a linear transformation of the coordinate vector. \nFor a sphere in ${tex`n`} dimensions having its center at ${tex`(x_1 \\cdots x_n)`} and radius ${tex`r`}, this is given by coordinates ${tex`(a_1, a_2 \\cdots a_{n+2})`}, where\n\n>${tex`a_i = \\frac{ x_i}{r}`} for ${tex`i \\leq n`},\n>\n>${tex`a_{n+1} = \\frac{(\\sum x^2_i) - r^2-1 }{2r}`},\n>\n>${tex`a_{n+2} = \\frac{(\\sum x^2_i) - r^2+1 }{2r}`}.\n\nAs an interesting fact, the curvature ${tex`\\kappa=\\frac{1}{r}`} of a sphere can be expressed linearly as ${tex`a_{n+2}-a_{n+1}`}:","pinned":false,"hidden":false,"mode":"md"},{"id":187,"value":"function curv(a) {\n return a[a.length - 1] - a[a.length - 2];\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":185,"value":"Thus, we can write the following two functions that compute the inversive and regular coordinates of a sphere.","pinned":false,"hidden":false,"mode":"md"},{"id":24,"value":"function invCoords(x, r) {\n const kappa = 1 / r;\n const dot2 = x.reduce((acc, d) => acc + d * d, 0);\n return [\n ...x.map((d) => d / r),\n ((dot2 - r * r - 1) * kappa) / 2,\n ((dot2 - r * r + 1) * kappa) / 2\n ];\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":27,"value":"function regCoords(a) {\n const n = a.length - 2;\n const kappa = curv(a);\n const r = 1 / kappa;\n return [a.slice(0, n).map((d) => d * r), r];\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":46,"value":"## A basis for the Apollonian packing\nBy applying a special transformation we can obtain a set of five spheres with integer inversive coordinates that are mutually tangent and form a basis for the Apollonian packing. This means that the Cartesian coordinates of a sphere in the basis is obtained by first transforming it into inversive coordinates and then back to regular coordinates. Note that this basis, as well as the special transformation appear in Ron Peikert's code, and are simpler than the one in the paper.","pinned":false,"hidden":false,"mode":"md"},{"id":49,"value":"basis = [\n [1, 0, 0, 0, 0],\n [0, 1, 0, 0, 0],\n [0, 0, 1, 0, 0],\n [0, 0, 0, 1, 0],\n [0, 0, 0, 0, 1]\n]","pinned":true,"hidden":false,"mode":"ojs"},{"id":53,"value":"specialToCartesian = {\n const rt2o2 = Math.sqrt(2) / 2;\n const rt6o2 = Math.sqrt(6) / 2;\n return (A, B, C, D, E) =>\n regCoords([\n (-B + C + D - E) * rt2o2,\n (B - C + D - E) * rt2o2,\n (B + C - D - E) * rt2o2,\n A - B - C - D - E,\n (B + C + D + E) * rt6o2\n ]);\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":143,"value":"Here's a plot of the 4 inner spheres. Note that the 5th sphere is not plotted by our utility function since it has negative radius. You can use the _mouse_ to move the camera around and zoom in/out (wheel).","pinned":false,"hidden":false,"mode":"md"},{"id":128,"value":"showSpheres(\n basis.map((s) => {\n let [center, radius] = specialToCartesian(...s);\n return { center, radius };\n }),\n visibility,\n invalidation\n)","pinned":false,"hidden":false,"mode":"ojs"},{"id":639,"value":"## Linear mappings for the spherical inversions\nThe special coordinates used in the basis have the added benefit that the next generation of Apollonian spheres can be expressed as linear transformations. Here are the 5 linear transformations that yield spheres in the Apollonian packing when applied to each of the 5 basis spheres.","pinned":false,"hidden":false,"mode":"md"},{"id":879,"value":"I = [\n (A, B, C, D, E) => [-A, A + B, A + C, A + D, A + E],\n (A, B, C, D, E) => [B + A, -B, B + C, B + D, B + E],\n (A, B, C, D, E) => [C + A, C + B, -C, C + D, C + E],\n (A, B, C, D, E) => [D + A, D + B, D + C, -D, D + E],\n (A, B, C, D, E) => [E + A, E + B, E + C, E + D, -E]\n]","pinned":true,"hidden":false,"mode":"ojs"},{"id":1000,"value":"## Duplicate removal\nIn order to produce the remaining spheres of the Apollonian packing, it suffices to apply the transformations `I` to each sphere in the basis, and repeat the process recursively. The process, however, may generate the same sphere multiple times. The code below implements additional tests that filter duplicates.","pinned":false,"hidden":false,"mode":"md"},{"id":844,"value":"function generate(basis) {\n const result = [];\n const rt6o2 = Math.sqrt(6) / 2;\n for (let j = 0; j < basis.length; j++) {\n let [A, B, C, D, E] = basis[j];\n const curvJ = (B + C + D + E) * rt6o2 - (A - B - C - D - E);\n for (let i of [0, 1, 2, 3, 4]) {\n const tmp = I[i](A, B, C, D, E);\n const [A1, B1, C1, D1, E1] = tmp;\n const curvI = (B1 + C1 + D1 + E1) * rt6o2 - (A1 - B1 - C1 - D1 - E1);\n if (\n curvI <= curvJ ||\n (i == 0 && (B1 < A1 || C1 < A1 || D1 < A1 || E1 < A1)) ||\n (i == 1 && (A1 <= B1 || C1 < B1 || D1 < B1 || E1 < B1)) ||\n (i == 2 && (A1 <= C1 || B1 <= C1 || D1 < C1 || E1 < C1)) ||\n (i == 3 && (A1 <= D1 || B1 <= D1 || C1 <= D1 || E1 < D1)) ||\n (i == 4 && (A1 <= E1 || B1 <= E1 || C1 <= E1 || D1 <= E1))\n )\n continue;\n result.push([A1, B1, C1, D1, E1]);\n }\n }\n return result;\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":1002,"value":"Here's the first and the second generation of spheres. Notice that the original basis could be considered the 0th generation.","pinned":false,"hidden":false,"mode":"md"},{"id":645,"value":"firstGen = generate(basis)","pinned":true,"hidden":false,"mode":"ojs"},{"id":1008,"value":"secondGen = generate(firstGen)","pinned":true,"hidden":false,"mode":"ojs"},{"id":667,"value":"Here's a display of these.","pinned":false,"hidden":false,"mode":"md"},{"id":993,"value":"showSpheres(\n [\n ...firstGen.map((s) => {\n let [center, radius] = specialToCartesian(...s);\n return { center, radius, color: \"red\" };\n }),\n ...secondGen.map((s) => {\n let [center, radius] = specialToCartesian(...s);\n return { center, radius, color: \"cyan\" };\n })\n ],\n visibility,\n invalidation\n)","pinned":false,"hidden":false,"mode":"ojs"},{"id":1085,"value":"
\n## Recursive display\n\nWith all pieces in place, we can compute and display the Apollonian packing up to a maximum number of number of recursions. We can also curb computation to discard spheres that are too small.","pinned":false,"hidden":false,"mode":"md"},{"id":429,"value":"viewof minRadiusCompute = Inputs.range([0, 0.2], {\n label: \"min radius to compute\",\n value: 0.01,\n step: 0.005\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":1072,"value":"viewof maxRadiusDisplay = Inputs.range([0, 0.5], {\n label: \"max radius for display\",\n value: 0.5,\n step: 0.01\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":433,"value":"viewof generations = Inputs.range([1, 10], {\n label: \"generations\",\n step: 1,\n value: 7\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":451,"value":"viewof maxCenterY = Inputs.range([-1, 1], {\n label: \"max center y\",\n step: 0.01\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":1095,"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 {\n label: \"drawing palette\",\n value: \"Accent\"\n }\n)","pinned":false,"hidden":false,"mode":"ojs"},{"id":308,"value":"{\n const colors = d3[\"scheme\" + colorPalette];\n const ncolors = colors.length;\n return showSpheres(\n apollonian\n .map((gen, igen) =>\n gen\n .filter((s) => s[1] <= maxRadiusDisplay && s[0][1] <= maxCenterY)\n .map((s) => ({\n center: s[0],\n radius: s[1],\n color: colors[igen % ncolors]\n }))\n )\n .flat(),\n visibility,\n invalidation\n );\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":1032,"value":"apollonian = {\n let gen = basis;\n let result = [basis.map((s) => specialToCartesian(...s))];\n let i = 0;\n for (let i = 1; i <= generations; i++) {\n const spheres = [];\n const nextGen = [];\n for (let s of gen) {\n const sph = specialToCartesian(...s);\n if (Math.abs(sph[1]) >= minRadiusCompute) {\n spheres.push(sph);\n nextGen.push(s);\n }\n }\n if (spheres.length == 0) break;\n gen = generate(nextGen);\n result.push(spheres);\n }\n return result;\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":1108,"value":"//\n// Uncomment to animate the max center y slider \n// { \n// for (let v = -1; v < 1; v += 0.01) {\n// viewof maxCenterY.value = v;\n// viewof maxCenterY.dispatchEvent(new CustomEvent(\"input\"));\n// yield;\n// await Promises.delay(100);\n// }\n// }","pinned":true,"hidden":false,"mode":"ojs"},{"id":349,"value":"## Sphere drawing\n\nHere's an utility function to display a bunch of spheres using [Three.js](https://threejs.org):","pinned":false,"hidden":false,"mode":"md"},{"id":355,"value":"THREE = {\n const THREE = (window.THREE = await require(\"three@0.119.1\"));\n await require(\"three@0.119.1/examples/js/controls/OrbitControls.js\").catch(\n () => {}\n );\n return THREE;\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":352,"value":"$31","pinned":true,"hidden":false,"mode":"ojs"},{"id":359,"value":"async function* showSpheres(spheres, visibility, invalidation) {\n let canvas = DOM.canvas(width, (width * 9) / 16);\n let {draw, dispose} = drawSpheres(spheres, canvas);\n invalidation.then(dispose)\n draw();\n yield canvas;\n while (true) {\n await visibility();\n await Promises.delay(1000 / 60);\n draw();\n yield canvas;\n }\n}","pinned":true,"hidden":false,"mode":"ojs"}],"files":"$W32","stdlib":"1","titling":"auto"}