2f:T53b,md`# Deformation with 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. The demo below allows you to experiment with RBFs by deforming a checkered cloth that represents the Euclidean plane. The red points mark the source points ${tex`X_i`}, whereas the blue points mark the corresponding target points ${tex`Y_i`}. You can drag the blue or the red points to change the RBF. Clicking on a red point removes it and its corresponding blue point from the RBF, while clicking on an empty region defines a new source/target point. The kernel can be selected in the menu below.`30:T726,{ let rbf = computeRbf(); let ctx = DOM.context2d(gridSize, gridSize); function render() { ctx.clearRect(0, 0, gridSize, gridSize); drawGrid(ctx, grid, rbf); drawPoints(ctx, controlPoints.src, 12, "red"); drawPoints(ctx, controlPoints.dst, 8, "blue"); } render(); let sel; let last; let mode; ctx.canvas.onmousedown = function (e) { sel = null; mode = null; let [mousex, mousey] = [e.offsetX, e.offsetY]; let select = (pts) => { for (let p of pts) { let [x, y] = p; let [dx, dy] = [x - mousex, y - mousey]; if (Math.sqrt(dx * dx + dy * dy) < 20) { sel = p; } } return sel; }; if (select(controlPoints.dst)) mode = "dragDst"; else if (select(controlPoints.src)) mode = "dragOrRemoveSrc"; else mode = "create"; rbf = computeRbf(); render(); last = [mousex, mousey]; }; ctx.canvas.onmouseup = function () { if (mode === "dragOrRemoveSrc") { let i = controlPoints.src.findIndex((x) => x == sel); if (i >= 0) { controlPoints.src.splice(i, 1); controlPoints.dst.splice(i, 1); } } else if (mode === "create") { controlPoints.src.push(last.slice()); controlPoints.dst.push(last.slice()); } rbf = computeRbf(); render(); mode = null; sel = null; }; ctx.canvas.onmousemove = function (e) { let [mousex, mousey] = [e.offsetX, e.offsetY]; if (!sel) return; if (mode === "dragOrRemoveSrc") mode = "dragSrc"; if (mode === "dragSrc" || mode === "dragDst") { let [pmousex, pmousey] = last; let [dx, dy] = [mousex - pmousex, mousey - pmousey]; sel[0] += dx; sel[1] += dy; rbf = computeRbf(); render(); } last = [mousex, mousey]; }; return ctx.canvas; }31:Tb38,/** * This is essentially the rbf module by Thibaut Séguy, but modified to use a selectable norm * See https://github.com/thibauts/rbf **/ RBF = { function RBF(points, values, distanceFunction, normFunction, epsilon) { let epsilon_ = epsilon; // `identity` here serves as a shorthand to allocate // the matrix nested array. const M = numeric.identity(points.length); const M_ = Matrix.eye(points.length); // First compute the point to point distance matrix // to allow computing epsilon if it's not provided for (let j = 0; j < points.length; j++) { for (let i = 0; i < points.length; i++) { M[j][i] = normFunction(points[i], points[j]); M_.set(j, i, normFunction(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); epsilon_ = sumM(M_) / (Math.pow(points.length, 2) - points.length); } // console.log(epsilon, epsilon_); // update the matrix to reflect the requested distance function for (let j = 0; j < points.length; j++) { for (let i = 0; i < points.length; i++) { M[j][i] = distanceFunction(M[j][i], epsilon); M_.set(j, i, distanceFunction(M_.get(j, i), epsilon_)); } } // reshape values by dimension/component let valuesByDimension = new Array(2); for (let i = 0; i < 2; i++) { valuesByDimension[i] = values.map((value) => value[i]); } // Compute basis functions weights by solving // the linear system of equations for each target component let w = new Array(2); let w_ = new Array(2); const LU = numeric.LU(M); for (let i = 0; i < 2; i++) { w[i] = numeric.LUsolve(LU, valuesByDimension[i]); w_[i] = mlMatrix.solve(M_, Matrix.columnVector(valuesByDimension[i])); } // console.log(w, w_); // The returned interpolant will compute the value at any point // by summing the weighted contributions of the input points function interpolant(p) { let distances = new Array(points.length); for (let i = 0; i < points.length; i++) { distances[i] = distanceFunction(normFunction(p, points[i]), epsilon); } const distances_ = new Matrix(distances.map((d) => [d])); let sums = new Array(2); let sums_ = new Array(2); for (let i = 0; i < 2; i++) { sums[i] = numeric.sum(numeric.mul(distances, w[i])); // console.log(distances, "x", w[i]); // console.log(distances_, "xx", w_[i]); sums_[i] = sumM(mulM(distances_, w_[i])); // sums_[i] = sumM(w_[i].mmul(distances_)); } // console.log(sums, sums_); return sums_; } interpolant.weights = w_; return interpolant; } return RBF; }32:[] 2d:{"title":"Deformation with RBFs","cells":[{"id":0,"value":"$2f","pinned":false,"hidden":false,"mode":"ojs"},{"id":93,"value":"viewof distanceFunction = Inputs.select(distanceFunctions, {\n format: (fn) => fn.label,\n value: distanceFunctions.find((fn) => fn.label === \"Thin-plate spline\"),\n label: \"RBF Kernel\"\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":40,"value":"$30","pinned":false,"hidden":false,"mode":"ojs"},{"id":166,"value":"viewof gridSize = Inputs.range([100, width], {\n value: width,\n step: 10,\n label: \"Grid size\"\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":157,"value":"viewof cellSize = Inputs.range([2, gridSize], {\n value: 32,\n step: 1,\n label: \"Cell size\"\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":182,"value":"viewof color1 = Inputs.color({ label: \"Color 1\", value: \"#d2c7c8\" })","pinned":false,"hidden":false,"mode":"ojs"},{"id":185,"value":"viewof color2 = Inputs.color({ label: \"Color 2\", value: \"#dcded1\" })","pinned":false,"hidden":false,"mode":"ojs"},{"id":249,"value":"function computeRbf() {\n let anchors = [\n [-gridSize * 10, -gridSize * 10],\n [-gridSize * 10, gridSize * 10],\n [gridSize * 10, gridSize * 10],\n [gridSize * 10, -gridSize * 10]\n ];\n let src = controlPoints.src.concat(anchors);\n let dst = controlPoints.dst.concat(anchors);\n return RBF(src, dst, distanceFunction.fn, euclideanNorm);\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":260,"value":"rbf = computeRbf()","pinned":true,"hidden":false,"mode":"ojs"},{"id":264,"value":"rbf([1, 2])","pinned":true,"hidden":false,"mode":"ojs"},{"id":53,"value":"mutable controlPoints = ({\n src: [\n [gridSize * 0.2, gridSize * 0.2],\n [gridSize * 0.8, gridSize * 0.2],\n [gridSize * 0.8, gridSize * 0.8],\n [gridSize * 0.2, gridSize * 0.8]\n ],\n dst: [\n [gridSize * 0.2, gridSize * 0.2],\n [gridSize * 0.8, gridSize * 0.2],\n [gridSize * 0.8, gridSize * 0.8],\n [gridSize * 0.2, gridSize * 0.8]\n ]\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":37,"value":"function drawGrid(ctx, grid, f) {\n for (let i = 0; i + 1 < grid.length; i++) {\n for (let j = 0; j + 1 < grid[i].length; j++) {\n let p = [grid[i][j], grid[i][j + 1], grid[i + 1][j + 1], grid[i + 1][j]];\n let q = p.map(f);\n\n // Draw first triangle\n ctx.fillStyle = color1;\n ctx.beginPath();\n ctx.moveTo(q[0][0], q[0][1]);\n ctx.lineTo(q[1][0], q[1][1]);\n ctx.lineTo(q[3][0], q[3][1]);\n ctx.fill();\n\n // Draw second triangle\n ctx.fillStyle = color2;\n ctx.beginPath();\n ctx.moveTo(q[1][0], q[1][1]);\n ctx.lineTo(q[2][0], q[2][1]);\n ctx.lineTo(q[3][0], q[3][1]);\n ctx.fill();\n }\n }\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":81,"value":"function drawPoints(ctx, points, radius = 8, color = \"red\") {\n ctx.fillStyle = color;\n for (let [x, y] of points) {\n ctx.beginPath();\n ctx.arc(x, y, radius, 0, 2 * Math.PI, false);\n ctx.fill();\n }\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":30,"value":"grid = {\n let grid = [];\n let d = cellSize;\n for (let x = 0, i = 0; x <= gridSize; i++, x += d) {\n grid[i] = [];\n for (let y = 0, j = 0; y <= gridSize; j++, y += d) {\n grid[i][j] = [x, y];\n }\n }\n return grid;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":340,"value":"## Porting to WebGL","pinned":false,"hidden":false,"mode":"md"},{"id":342,"value":"To run the resulting interpolation in WebGL, we need to port this function to glsl:\n\n```js\nfunction interpolant(p) {\n var distances = new Array(points.length);\n for (var i = 0; i < points.length; i++) {\n distances[i] = distance(norm(p, points[i]), epsilon);\n }\n\n var sums = new Array(D);\n for (var i = 0; i < D; i++) {\n sums[i] = numeric.sum(numeric.mul(distances, w[i]));\n }\n\n return sums;\n}\n```\n","pinned":false,"hidden":false,"mode":"md"},{"id":275,"value":"## Distance functions","pinned":false,"hidden":false,"mode":"md"},{"id":278,"value":"function euclideanNorm(point1, point2) {\n const sub = [point2[0] - point1[0], point2[1] - point1[1]];\n const norm = Math.sqrt(sub[0] ** 2 + sub[1] ** 2);\n return norm;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":282,"value":"\n function distanceLinear(r) {\n return r;\n }\n","pinned":false,"hidden":false,"mode":"ojs"},{"id":284,"value":"function distanceCubic(r) {\n return Math.pow(r, 3);\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":286,"value":"\n function distanceQuintic(r) {\n return Math.pow(r, 5);\n }\n","pinned":false,"hidden":false,"mode":"ojs"},{"id":288,"value":"function distanceThinPlate(r) {\n if (r === 0) return 0;\n return Math.pow(r, 2) * Math.log(r);\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":290,"value":"function distanceGaussian(r, epsilon) {\n return Math.exp(-Math.pow(r / epsilon, 2));\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":281,"value":"function distanceInverseMultiquadric(r, epsilon) {\n return 1.0 / Math.sqrt(Math.pow(r / epsilon, 2) + 1);\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":292,"value":"function distanceMultiquadric(r, epsilon) {\n return Math.sqrt(Math.pow(r / epsilon, 2) + 1);\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":314,"value":"distanceFunctions = [\n { label: \"Linear\", fn: distanceLinear },\n { label: \"Cubic\", fn: distanceCubic },\n { label: \"Quintic\", fn: distanceQuintic },\n { label: \"Thin-plate spline\", fn: distanceThinPlate },\n { label: \"Gaussian\", fn: distanceGaussian },\n { label: \"Inverse multiquadric\", fn: distanceInverseMultiquadric },\n { label: \"Multiquadric\", fn: distanceMultiquadric }\n]","pinned":false,"hidden":false,"mode":"ojs"},{"id":315,"value":"normFunctions = ({\n euclidean: euclideanNorm\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":7,"value":"md`## Libraries`","pinned":false,"hidden":false,"mode":"ojs"},{"id":228,"value":"[numeric](https://github.com/sloisel/numeric) is 11 years old. We can consider moving to [ml-matrix](https://github.com/mljs/matrix) instead. The `RBF` function uses the following functions from numeric:\n\n- `numeric.identity`\n- `numeric.sum`\n- `numeric.mul`\n- `numeric.norm2`\n- `numeric.sub`\n- `numeric.LU`\n- `numeric.LUsolve`\n\nIt seems ml-matrix supports all these functions: \n\n| numeric | ml-matrix |\n|:-------------------|:-------------------|\n| `numeric.identity` | trivial |\n|`numeric.sum` | `Matrix.add(A, B)` |\n| `numeric.mul` | `A.mmul(B)` |\n| `numeric.norm2` | `A.norm()` |\n| `numeric.sub` | `Matrix.sub(A, B)` |\n| `numeric.LU` | `new LuDecomposition(A)` |\n| `numeric.LUsolve` | `solve` |","pinned":false,"hidden":false,"mode":"md"},{"id":420,"value":"// alternative for numeric.sum(M)\nfunction sumM(M) {\n let sum = 0;\n for (let j = 0; j < M.rows; j++) {\n for (let i = 0; i < M.columns; i++) {\n sum += M.get(j, i);\n }\n }\n return sum;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":500,"value":"// alternative for numeric.mul(M)\nfunction mulM(M1, M2) {\n let mul = Matrix.zeros(M1.rows, M1.columns);\n for (let j = 0; j < M1.rows; j++) {\n for (let i = 0; i < M1.columns; i++) {\n mul.set(j, i, M1.get(j, i) * M2.get(j, i));\n }\n }\n return mul;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":13,"value":"$31","pinned":true,"hidden":false,"mode":"ojs"},{"id":550,"value":"gcps2 = [\n {\n world: [4.4091165, 51.9017125],\n image: [1344, 4098]\n },\n {\n world: [4.5029222, 51.9164451],\n image: [4440, 3441]\n },\n {\n world: [4.4764224, 51.897309],\n image: [3549, 4403]\n },\n {\n world: [4.4199066, 51.9391509],\n image: [1794, 2130]\n },\n {\n world: [4.4775683, 51.9324358],\n image: [3656, 2558]\n }\n]","pinned":false,"hidden":false,"mode":"ojs"},{"id":526,"value":"gcps = [\n {\n image: [518, 991],\n world: [4.9516614, 52.4633102]\n },\n {\n image: [4345, 2357],\n world: [5.0480391, 52.5123762]\n },\n {\n image: [2647, 475],\n world: [4.9702906, 52.5035815]\n }\n]","pinned":false,"hidden":false,"mode":"ojs"},{"id":493,"value":"gcpRbf = {\n let src = gcps2.map((gcp) => gcp.image);\n let dst = gcps2.map((gcp) => gcp.world);\n return RBF(src, dst, distanceFunction.fn, euclideanNorm);\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":537,"value":"gcpRbf([518, 991])","pinned":true,"hidden":false,"mode":"ojs"},{"id":15,"value":"numeric = require(\"https://bundle.run/numeric@1.2.6\")","pinned":false,"hidden":false,"mode":"ojs"},{"id":407,"value":"Matrix = mlMatrix.Matrix","pinned":false,"hidden":false,"mode":"ojs"},{"id":410,"value":"mlMatrix = import(\"ml-matrix\")","pinned":false,"hidden":false,"mode":"ojs"}],"files":"$W32","stdlib":"1","titling":"auto"}