30:T75d,md`# 3D RBFs A [Radial basis function (RBF)](https://en.wikipedia.org/wiki/Radial_basis_function) provides a neat way to solve the so-called *scattered data interpolation* problem. In a nutshell, suppose you want a function ${tex`f(X)`} about which you nothing other than its value at a few points of the the domain, i.e., you know that ${tex`f(X_i) = Y_i`} for sets of points ${tex`\{X_0, X_1 \cdots X_n \}`} and ${tex`\{Y_0, Y_1 \cdots Y_n\}`}. The "radial" in RBF is because ${tex`f(X)`} depends on the distance from ${tex`X`} to all other domain points ${tex`X_i`}. The *kernel* of a RBF controls how distances are measured. For instance, a \`linear\` kernel means that the Euclidean norm ${tex`d_i (X) = \| X - X_i \|`} is used, whereas a \`quintic\` kernel employs the Euclidean norm to the fifth power, and so forth. A surface can be modeled using RBFs by adding points ${tex`X_i`} that we wish to interpolate and associating them with ${tex`Y_i = 0`}. In the implementation below, these are generated at random within a cube ${tex`-1 \leq x,y,z \leq 1`}. We must also add at least one point not associated with zero so that the equation system is fully determined. Here we add 8 points at the corners of a bigger cube and map them to the distance to the closest point that is mapped to zero. In this notebook, the isosurfaces ${tex`f(X) = 0`} are polygonized using an approach called [Surface nets](https://0fps.net/2012/07/12/smooth-voxel-terrain-part-2/) as implemented by [Mikola Lysenko](https://github.com/mikolalysenko/isosurface). This approach requires a rather fine sampling of the field, and since RBF evaluation is costly, we try to speed up things a bit by using a coarse 3D grid sampling in order to discard regions of space that are unlikely to contain zeroes of the function. You can turn off this optimization by unchecking *Use adaptive sampling* below. `31:T75d,// // Returns a scalar field sampling function f(p) that speeds up the computation of func(p). // The idea is to presample the field within the rectangular region between min and max // with a n x n x n grid. If p is in a cell that likely does not contain a zero, an // approximation is returned. // function adaptiveSampler(min, max, n, func) { const dn = d3.zip(min, max).map(([xmin, xmax]) => (xmax - xmin) / n); // Sample the coarse grid const grid = []; for (let i = 0; i <= n; i++) { const x = i * dn[0] + min[0]; const a = []; grid.push(a); for (let j = 0; j <= n; j++) { const y = j * dn[1] + min[1]; const b = []; a.push(b); for (let k = 0; k <= n; k++) { const z = k * dn[2] + min[2]; b.push(func([x, y, z])[0]); } } } // Collect cell addresses straddling the zero values of func const hash = (i, j, k) => (i * (n + 1) + j) * (n + 1) + k; const zeroCells = new Set(); for (let i = 0; i < n; i++) { for (let j = 0; j < n; j++) { for (let k = 0; k < n; k++) { let vmin, vmax; vmin = vmax = grid[i][j][k]; for (let [di, dj, dk] of [ [0, 0, 1], [0, 1, 0], [0, 1, 1], [1, 0, 0], [1, 0, 1], [1, 1, 0], [1, 1, 1] ]) { let v = grid[i + di][j + dj][k + dk]; vmin = Math.min(vmin, v); vmax = Math.max(vmax, v); } if (vmin < 0 && vmax > 0) { zeroCells.add(hash(i, j, k)); } } } } return (p) => { const index = [], rem = []; for (let i of [0, 1, 2]) { index[i] = ~~((p[i] - min[i]) / dn[i]); rem[i] = (p[i] - min[i]) % dn[i]; } const [i, j, k] = index; if (Math.max(...rem) < 0.01 || !zeroCells.has(hash(i, j, k))) return [grid[i][j][k]]; return func(p); }; }32:T456,testAdaptiveSampler = { return; // Comment this line to enable const [min, max] = [ [-1, -1, -1], [1, 1, 1] ]; const n = gridSize; const adaptive = adaptiveSampler(min, max, Math.max(8, n / 8), interpolant); const dn = d3.zip(min, max).map(([xmin, xmax]) => (xmax - xmin) / n); let timea = performance.now(); for (let i = 0; i <= n; i++) { const x = i * dn[0] + min[0]; for (let j = 0; j <= n; j++) { const y = j * dn[1] + min[1]; for (let k = 0; k <= n; k++) { const z = k * dn[2] + min[2]; const a = adaptive([x, y, z])[0]; //const b = interpolant([x, y, z])[0]; } } } timea = performance.now() - timea; let timeb = performance.now(); for (let i = 0; i <= n; i++) { const x = i * dn[0] + min[0]; for (let j = 0; j <= n; j++) { const y = j * dn[1] + min[1]; for (let k = 0; k <= n; k++) { const z = k * dn[2] + min[2]; //const a = adaptive([x, y, z])[0]; const b = interpolant([x, y, z])[0]; } } } timeb = performance.now() - timeb; return { timea, timeb }; }33:T6f7,isoSurfaceApp = { const frag = `#version 300 es precision highp float; in vec3 vertex_normal; out vec4 fragColor; uniform vec4 color; uniform vec3 light_dir; void main () { fragColor = vec4(color.xyz *(0.2+0.8*(max(0.,dot(vertex_normal,light_dir)))),1.0); }`; const vert = `#version 300 es layout(location=0) in vec4 position; layout(location=1) in vec3 normal; out vec3 vertex_normal; uniform mat4 modelview; uniform mat4 projection; uniform float scaleFactor; void main () { vec4 worldpos = modelview*(scaleFactor*position); gl_Position = projection*worldpos; vertex_normal = normalize((modelview*vec4(normal,0.0)).xyz); }`; const canvas = isoSurfaceCanvas; const { width: w, height: h } = canvas; const app = picogl .createApp(canvas) .clearColor(0.5, 0.5, 0.5, 1) .enable(picogl.DEPTH_TEST); function vtxArray(mesh) { const positions = new Float32Array(mesh.positions.flat()); const normals = new Float32Array(mesh.vertexNormals.flat()); const indices = new Uint16Array(mesh.faces.flat()); const positionsBuffer = app.createVertexBuffer(picogl.FLOAT, 3, positions); const normalsBuffer = app.createVertexBuffer(picogl.FLOAT, 3, normals); const indexBuffer = app.createIndexBuffer(picogl.UNSIGNED_SHORT, indices); return app .createVertexArray() .vertexAttributeBuffer(0, positionsBuffer) .vertexAttributeBuffer(1, normalsBuffer) .indexBuffer(indexBuffer); } const program = app.createProgram(vert, frag); const drawCallSurface = app.createDrawCall(program, vtxArray(mesh)); const drawCallBall = app.createDrawCall(program, vtxArray(ballMesh)); return { app, drawCallSurface, drawCallBall }; }34:T63d,renderIsoSurface = { const canvas = isoSurfaceCanvas; const { app, drawCallSurface, drawCallBall } = isoSurfaceApp; const arcball = new Arcball(canvas); const refresh = () => { app.clear(); const modelview = mat4.mul( [], mat4.fromTranslation([], [0, 0, -5]), arcball.transform ); const projection = mat4.perspective( [], (45 * Math.PI) / 180, canvas.width / canvas.height, 0.01, 100 ); drawCallSurface .uniform("modelview", modelview) .uniform("light_dir", vec3.normalize([], [1, 1, 1])) .uniform("color", [1.0, 0.0, 0.0, 1.0]) .uniform("scaleFactor", 1.0) .uniform("projection", projection) .draw(); data.points.forEach((p, i) => { drawCallBall .uniform( "modelview", mat4.multiply( [], modelview, mat4.multiply( [], mat4.fromTranslation([], p), mat4.fromScaling([], [0.05, 0.05, 0.05]) ) ) ) .uniform("light_dir", vec3.normalize([], [1, 1, 1])) .uniform("color", i < 8 ? [0, 0, 1, 1] : [1.0, 1.0, 0.0, 1.0]) .uniform("scaleFactor", 1.0) .uniform("projection", projection) .draw(); }); }; canvas.onmousedown = (e) => arcball.mousedown(e.offsetX, canvas.height - e.offsetY); canvas.onmousemove = (e) => { if (e.buttons) { arcball.mousemove(e.offsetX, canvas.height - e.offsetY); refresh(); } }; canvas.onmouseup = (e) => arcball.mouseup(); refresh(); }35:Te9f,/** * This is essentially the rbf module by Thibaut Séguy, but modified to return the RBF weights, rather than a * function **/ RBF = { var distanceFunctions = { linear: distanceLinear, cubic: distanceCubic, quintic: distanceQuintic, "thin-plate": distanceThinPlate, gaussian: distanceGaussian, "inverse-multiquadric": distanceInverseMultiquadric, multiquadric: distanceMultiquadric }; var normFunctions = { euclidean: euclideanNorm }; var distance = distanceFunctions.linear; var norm = normFunctions.euclidean; function RBF(points, values, distanceFunction, normFunction, epsilon) { if (distanceFunction) { distance = typeof distanceFunction !== "string" ? distanceFunction : distanceFunctions[distanceFunction]; } if (normFunction) { norm = typeof normFunction !== "string" ? normFunction : normFunctions[normFunction]; } // `identity` here serves as a shorthand to allocate // the matrix nested array. var M = numeric.identity(points.length); // First compute the point to point distance matrix // to allow computing epsilon if it's not provided for (var j = 0; j < points.length; j++) { for (var i = 0; i < points.length; i++) { M[j][i] = norm(points[i], points[j]); } } // if needed, compute espilon as the average distance between points if (epsilon === undefined) { epsilon = numeric.sum(M) / (Math.pow(points.length, 2) - points.length); } // update the matrix to reflect the requested distance function for (var j = 0; j < points.length; j++) { for (var i = 0; i < points.length; i++) { M[j][i] = distance(M[j][i], epsilon); } } // determine dimensionality of target values var sample = values[0]; var D = typeof sample === "number" ? 1 : sample.length; // generalize to vector values if (D === 1) { values = values.map(function (value) { return [value]; }); } // reshape values by component var tmp = new Array(D); for (var i = 0; i < D; i++) { tmp[i] = values.map(function (value) { return value[i]; }); } values = tmp; // Compute basis functions weights by solving // the linear system of equations for each target component var w = new Array(D); var LU = numeric.LU(M); for (var i = 0; i < D; i++) { w[i] = numeric.LUsolve(LU, values[i]); } // The returned interpolant will compute the value at any point // by summing the weighted contributions of the input points function interpolant(p) { var distances = new Array(points.length); for (var i = 0; i < points.length; i++) { distances[i] = distance(norm(p, points[i]), epsilon); } var sums = new Array(D); for (var i = 0; i < D; i++) { sums[i] = numeric.sum(numeric.mul(distances, w[i])); } return sums; } return { w, interpolant }; } function euclideanNorm(pa, pb) { return numeric.norm2(numeric.sub(pb, pa)); } function distanceLinear(r) { return r; } function distanceCubic(r) { return Math.pow(r, 3); } function distanceQuintic(r) { return Math.pow(r, 5); } function distanceThinPlate(r) { if (r === 0) return 0; return Math.pow(r, 2) * Math.log(r); } function distanceGaussian(r, epsilon) { return Math.exp(-Math.pow(r / epsilon, 2)); } function distanceInverseMultiquadric(r, epsilon) { return 1.0 / Math.sqrt(Math.pow(r / epsilon, 2) + 1); } function distanceMultiquadric(r, epsilon) { return Math.sqrt(Math.pow(r / epsilon, 2) + 1); } return RBF; }36:T475,class Arcball { constructor(canvas) { this.canvas = canvas; this.center = [canvas.width / 2, canvas.height / 2]; this.radius = (Math.min(canvas.width, canvas.height) / 2) * 0.9; this.prev = []; this.transform = mat4.create(); this.mouseIsDown = false; } ballPos(p) { let u = vec2.sub([], p, this.center); if (vec2.len(u) >= this.radius) return vec3.normalize([], [...u, 0]); let z = Math.sqrt(this.radius ** 2 - u[0] ** 2 - u[1] ** 2); return vec3.normalize([], [...u, z]); } mousedown(x, y) { this.mouseIsDown = true; this.prev = [x, y]; } mousemove(x, y) { if (!this.mouseIsDown) return this.transform; let mouse = [x, y]; let u = this.ballPos(this.prev); let v = this.ballPos(mouse); this.prev = mouse; let axis = vec3.cross([], u, v); let ang = Math.asin(vec3.len(axis)); let rot = mat4.fromRotation([], ang, vec3.normalize([], axis)); [this.u, this.v, this.ang, this.axis] = [u, v, ang, axis]; if (rot) mat4.mul(this.transform, rot, this.transform); return this.transform; } mouseup(x, y) { this.mouseIsDown = false; } }37:[] 2e:{"title":"3D RBFs","cells":[{"id":0,"value":"$30","pinned":false,"hidden":false,"mode":"ojs"},{"id":307,"value":"viewof kernel = Inputs.select(\n [\n \"linear\",\n \"cubic\",\n \"quintic\",\n \"thin-plate\",\n \"gaussian\",\n \"multiquadric\",\n \"inverse-multiquadric\"\n ],\n { label: \"RBF Kernel\", value: \"linear\" }\n)","pinned":false,"hidden":false,"mode":"ojs"},{"id":585,"value":"viewof gridSize = Inputs.select([16, 32, 64, 128], {\n label: \"grid size\",\n value: 32\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":589,"value":"viewof useAdaptive = Inputs.toggle({\n label: \"Use adaptive sampling\",\n value: true\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":643,"value":"viewof npoints = Inputs.select([4, 8, 16, 32], { label: \"points\", value: 8 })","pinned":false,"hidden":false,"mode":"ojs"},{"id":597,"value":"viewof randomize = html`