29:null 2a:["$","title",null,{"children":"Thin Plate Spline interpolation in Allmaps/transform - Test for Helmert transform | Manuel Claeys Bouuaert | Observable"}] 2f:T681,## Demo The demo below allows you to experiment with RBF interpolation by deforming a triangulated cloth that represents the Euclidean plane. - The **red points** 🔴 mark the **source points** ${tex`(x_i, y_i)`}, whereas the **blue points** 🔵 mark the corresponding **target points** ${tex`(x'_i, y'_i)`}. - You can **drag** 👇 the blue or the red points to change the interpolation. 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 **RBF interpolation algorithm** ✨ creates an interpolation algorithm based on the ground controle points and selected kernel, then transforms all triangle coordinates and draws the triangles in the right place. - You can **switch between various kernels** in the menu below. - You can turn the **affine component** off and on. Strictly speaking this component only makes sense for Thin Plate Splines, but it smoothens the results for the other kernels as well. - You can turn on and of the usage of **anchor points** far outside the view. These were added in the original code to make the cloth stay in place even for a low number of ground controle points. You will see that adding the affine component reduce the need for such points. In Allmaps, these triangles will be the result of the triangulation of our map image, the red points will be the controle points on the image and the blue points the corresponding controle points on the world map. These will be pixel coordinates and (lon, lat) coordinates respectively (so they'll be in a very different range!). ## Give it a try! ## 🏁 🤸‍♀️ ✨30:T736,{ let rbf = computeHelmert(); 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 = computeHelmert(); 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 = computeHelmert(); 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 = computeHelmert(); render(); } last = [mousex, mousey]; }; return ctx.canvas; }31:T5eb,## Why Thin Plate Spline interpolation A Thin Plate Spline transformation, also known as rubber sheeting, is a technique to draw a map image over the world map and make it fit in the right place. Other techniques exist to do this, like the affine transformation currently implemented in Allmaps. Unlike the affine transformation, however, TPS allows for an **exact** transformation. The transformation is exact in the sense that the control points of the map image (pixel coordinates) end up exactly where their corresponding controle points were identified in on the world map ((lon, lat) coordinates). All other pixels are then interpolated in between. There are many such exact transformations, each with their own characteristics. The specificity about TPSs is that they are **smooth**. They are the result of minimising the so called 'bending energy' of the transformation, and therefore behave like you would expect when manipulating, stretching or compressing a physical object. TPS are the preferred technique for many researchers in historical cartography. Using a TPS transformation allows us transform maps with large distortions (caused by imperfect cartographic techniques and/or originating from bad physical conservation of the map) or maps with unknown reference systems and superimpose then on the current map, **enabling comparative historical research** in each area of the map: the more controle points you add in a certain area, the more the map will correspond to the underlying world map.32:T1090,## The Maths behind the code Let's find out what is happening here. _Don't be shy!_ We start of with a function that describes a transformation from 2D to 2D: ${tex`F(x,y)`}. We know nothing about this function other than some controle points for which we have the values in the domain and codomain, i.e. we know that ${tex`F(x_i,y_i) = (x'_i, y'_i)`} for sets of controle points ${tex`\{(x_0,y_0), (x_1,y_1) \cdots (x_N,y_N) \}`} and their corresponding values ${tex`\{(x'_0,y'_0), (x'_1,y'_1) \cdots (x'_N,y'_N) \}`}. The idea then is to approximate this unknown function by a sum of smaller functions ${tex`\varphi_i(x,y)`}, one such function around each of the ${tex`N`} controle points. The "**radial**" in 'Radial Basis Functions' is because these smaller functions each depend only on the distance from ${tex`(x,y)`} to their corresponding controle point ${tex`(x_i,y_i)`}. This smaller function is also called a 'kernel'. There are many such kernels we could use. For instance, a 'linear kernel' means that the Euclidean norm, which is the standard distance function, ${tex`d_i(x, y) = \| (x,y) - (x_i,y_i)\|`} is used directly as the kernel. More formally, we have an unknown function ${tex`F(x,y): \mathbb{R^2} \to \mathbb{R^2}`}, with it's two parts ${tex`f_x(x,y)`} and ${tex`f_y(x,y): \mathbb{R} \to \mathbb{R^2}`}, which we approximate as ${tex`\hat{f}_x(x,y)`} and ${tex`\hat{f}_y(x,y)`}. We start of knowing ${tex`N`} controle points ${tex`(x_i, y_i) \to (x'_i, y'_i)`} and picking radial kernel ${tex`\varphi_i(r)`} with ${tex`r = \|(x, y)\|`}. Then, we set up the approximating functions in the following way: ${tex` \left\{ \begin{array}{ll} \hat{f}_x(x, y) = \sum_{i=1}^N w_{i, x} \, \varphi_i(\|(x, y) - (x_i, y_i)\|) \\ \hat{f}_y(x, y) = \sum_{i=1}^N w_{i, y} \, \varphi_i(\|(x, y) - (x_i, y_i)\|) \end{array} \right. `} Here there is a basis function ${tex`\varphi_i(x,y)`} for each controle point, and each basis function is multiplied by weight parameter ${tex`w_{i,x}`} or ${tex`w_{i,y}`}. By setting ${tex`\hat{f}_x(x, y) = x'_i`} and ${tex`\hat{f}_y(x, y) = y'_i`} for each of the ${tex`N`} controle points **this results in two systems of ${tex`N`} equations that can be solved**. Solving these systems will result in the weights, which in turn fully define the two approximating functions and hence allow us to compute the location ${tex`(x', y')`} of any random point ${tex`(x, y)`}. Nice! In the case of **Thin Plate Spline** interpolation the RBF that is picked is ${tex`\varphi_i(r) = r^2 \ln(r)`} with ${tex`r = \|(x, y)\|`}. This function is what you get if you look for the function that that [minimises the 'bending energy'](https://en.wikipedia.org/wiki/Thin_plate_spline) of the transformation. For this kernel function it makes sense to add an **affine component** (or first order polynomial transformation) to the approximating functions. The approximation thus becomes: ${tex` \left\{ \begin{array}{ll} \hat{f}_x(x, y) = a_{0, x} + a_{x, x}x + a_{y, x}y + \sum_{i=1}^N w_{i, x} \, \varphi_i(\|(x, y) - (x_i, y_i)\|) \\ \hat{f}_y(x, y) = a_{0, y} + a_{x, y}x + a_{y, y}y + \sum_{i=1}^N w_{i, y} \, \varphi_i(\|(x, y) - (x_i, y_i)\|) \end{array} \right. `} Where ${tex`a_{0, x}, a_{x, x}, ...`} are the six parameters of a 2D [affine transformation](https://en.wikipedia.org/wiki/Affine_transformation). These represent the translation, rotation and sheering which are all the transformations that can be done without adding bending energy. (Since we have three more parameters we also need three more equations to complete the system (for both ${tex`\hat{f}_x`} and ${tex`\hat{f}_y`}): we demand that ${tex`\sum_{i=1}^N w_{i, x} = 0`}, ${tex`\sum_{i=1}^N w_{i, x}x_i = 0`} and ${tex`\sum_{i=1}^N w_{i, x}y_i = 0`}, and analogously for ${tex`\hat{f}_y`}.) **So, implementing this in code** we'll have to get the ground controle point coordinates, evaluate the distances between all pairs of ground control points and then evaluate the kernel function for all these distances, build the system of equations with this information, solve it and build the interpolation functions using the resulting weights (and affine parameters). To see how this is done, go browse through the code below!33:T1c42,RBF = { function RBF( sourcePoints, destinationPoints, distanceFunction, normFunction, epsilon ) { const dimension = sourcePoints.length; // Notes on types: // // 'sourcePoints' and 'destinationPoints' are Arrays // sourcePoints = [[x0, y0], [x1, y1], ...] // destinationPoints = [[x'0, y0], [x'1, y'1], ...] // // 'destinationPointsMatrices' and 'weightsMatrices' is an Array of Matrices // destinationPointsMatrices = [Matrix([[x'0], [x'1], ...]), Matrix([[y'0], [y'1], ...])] // destinationPointsMatrices[i] is a Matrix of (dimension, 1) // Make destinationPointsMatrices let destinationPointsMatrices = new Array(2); if (useAffine) { destinationPoints = [...destinationPoints, [0, 0], [0, 0], [0, 0]]; } for (let i = 0; i < 2; i++) { destinationPointsMatrices[i] = mlMatrix.Matrix.columnVector( destinationPoints.map((value) => value[i]) ); } // Compute distancesMatrix, containing the point to point distances between all controle points const distancesMatrix = mlMatrix.Matrix.zeros(dimension, dimension); for (let i = 0; i < dimension; i++) { for (let j = 0; j < dimension; j++) { distancesMatrix.set( i, j, normFunction(sourcePoints[i], sourcePoints[j]) ); } } // If it's not provided, compute espilon as the average distance between the controle points if (epsilon === undefined) { epsilon = distancesMatrix.sum() / (Math.pow(dimension, 2) - dimension); } // Update the matrix to reflect the requested distance function // Note: for distance functions which do not take an epsilon value this step does not alter the matrix. In production this can be skipped for such distance functions. for (let i = 0; i < dimension; i++) { for (let j = 0; j < dimension; j++) { distancesMatrix.set( i, j, distanceFunction(distancesMatrix.get(i, j), epsilon) ); } } // Extend the point distance matrix to include the affine transformation let affineCoefsMatrix = mlMatrix.Matrix.zeros(dimension, 3); let distancesAndAffineCoefsMatrix = mlMatrix.Matrix.zeros( dimension + 3, dimension + 3 ); if (useAffine) { // Construct Nx3 Matrix affineCoefsMatrix // 1 x0 y0 // 1 x1 y1 // 1 x2 y2 // ... for (let i = 0; i < dimension; i++) { affineCoefsMatrix.set(i, 0, 1); affineCoefsMatrix.set(i, 1, sourcePoints[i][0]); affineCoefsMatrix.set(i, 2, sourcePoints[i][1]); } // Combine distancesMatrix and affineCoefsMatrix into new matrix distancesAndAffineCoefsMatrix // Note: mlMatrix has no knowledge of block matrices, but this approach is good enough // To speed this up, we could maybe use distancesMatrix.addRow() and distancesMatrix.addVector() for (let i = 0; i < dimension + 3; i++) { for (let j = 0; j < dimension + 3; j++) { if (i < dimension && j < dimension) { distancesAndAffineCoefsMatrix.set(i, j, distancesMatrix.get(i, j)); } else if (i >= dimension && j < dimension) { distancesAndAffineCoefsMatrix.set( i, j, affineCoefsMatrix.transpose().get(i - dimension, j) ); } else if (i < dimension && j >= dimension) { distancesAndAffineCoefsMatrix.set( i, j, affineCoefsMatrix.get(i, j - dimension) ); } } } } // Compute basis functions weights by solving the linear system of equations for each target component // Note: for this example the system is solved with and without the affine part to allow comparing between both situation. This is, of course, not efficient in production, where it's best to make a choice. let weightsMatrices = new Array(2); for (let i = 0; i < 2; i++) { if (useAffine) { weightsMatrices[i] = mlMatrix.solve( distancesAndAffineCoefsMatrix, destinationPointsMatrices[i] ); } else { weightsMatrices[i] = mlMatrix.solve( distancesMatrix, destinationPointsMatrices[i] ); } } // The interpolant function will compute the value at any point. function interpolant(newSourcePoint) { // First make a column matrix with the distances of that point to all controle points let newDistancesMatrix = mlMatrix.Matrix.zeros(dimension, 1); for (let i = 0; i < dimension; i++) { newDistancesMatrix.set( i, 0, distanceFunction( normFunction(newSourcePoint, sourcePoints[i]), epsilon ) ); } // Then it computes the interpolated value by summing the weighted contributions of the input points let newDestinationPoint = new Array(2); for (let i = 0; i < 2; i++) { // Apply the weights to the new distances // Note: don't consider the last three weights who are there for the affine part newDestinationPoint[i] = mulElements( newDistancesMatrix, weightsMatrices[i].selection([...Array(dimension).keys()], [0]) ).sum(); // Add the affine part if (useAffine) { let a0 = weightsMatrices[i].get(dimension, 0); let ax = weightsMatrices[i].get(dimension + 1, 0); let ay = weightsMatrices[i].get(dimension + 2, 0); newDestinationPoint[i] += a0 + ax * newSourcePoint[0] + ay * newSourcePoint[1]; } } return newDestinationPoint; } let affineCoefs = [null, null, null, null, null, null]; if (useAffine) { affineCoefs = [ weightsMatrices[0].get(dimension, 0), // a0_x weightsMatrices[0].get(dimension + 1, 0), // ax_x weightsMatrices[0].get(dimension + 2, 0), // ay_x weightsMatrices[1].get(dimension, 0), // a0_y weightsMatrices[1].get(dimension + 1, 0), // ax_y weightsMatrices[1].get(dimension + 2, 0) // ay_y ]; } // This value has to be one for an affine transformation let affineCheck = affineCoefs[1] ** 2 * affineCoefs[5] ** 2 - affineCoefs[3] ** 2 * affineCoefs[4] ** 2; // This value has to be zero for a Helmert transformation, which affine transformations are not in general let helmertCheck = affineCoefs[1] - affineCoefs[5]; // This value has to be zero for a Helmert transformation, which affine transformations are not in general let helmertCheck2 = affineCoefs[2] + affineCoefs[4]; let translation = [null, null]; if (useAffine) { translation = [affineCoefs[0], affineCoefs[3]]; } // Only valid for Helmert transformations // let scale = null; // if (useAffine) { // scale = Math.sqrt(affineCoefs[1] ** 2 + affineCoefs[5] ** 2); // } // Only valid for Helmert transformations // let rotation = null; // if (useAffine) { // rotation = (Math.atan2(affineCoefs[4] / affineCoefs[2]) * 180) / Math.PI; // } return interpolant; } return RBF; }34:T80c,Helmert = { function Helmert(sourcePoints, destinationPoints) { const pointCount = sourcePoints.length; const destinationPointsMatrix = mlMatrix.Matrix.columnVector( destinationPoints.flat() ); // Construct 2Nx4 Matrix helmertCoefsMatrix // 1 0 x0 -y0 // 0 1 y0 x0 // 1 0 x1 -y1 // 0 1 y1 x1 // ... const helmertCoefsMatrix = mlMatrix.Matrix.zeros(2 * pointCount, 4); for (let i = 0; i < pointCount; i++) { helmertCoefsMatrix.set(2 * i, 0, 1); helmertCoefsMatrix.set(2 * i, 1, 0); helmertCoefsMatrix.set(2 * i, 2, sourcePoints[i][0]); helmertCoefsMatrix.set(2 * i, 3, -sourcePoints[i][1]); helmertCoefsMatrix.set(2 * i + 1, 0, 0); helmertCoefsMatrix.set(2 * i + 1, 1, 1); helmertCoefsMatrix.set(2 * i + 1, 2, sourcePoints[i][1]); helmertCoefsMatrix.set(2 * i + 1, 3, sourcePoints[i][0]); } // Compute helmert parameters by solving the linear system of equations for each target component // Will result in a Matrix([[t_x], [t_y], [m], [n]]) const pseudoInverseHelmertCoefsMatrix = mlMatrix.pseudoInverse(helmertCoefsMatrix); const helmertParametersMatrix = pseudoInverseHelmertCoefsMatrix.mmul( destinationPointsMatrix ); // The interpolant function will compute the value at any point. function interpolate(newSourcePoint) { if (!helmertParametersMatrix) { throw new Error("Helmert parameters not computed"); } // Compute the interpolated value by applying the helmert coefficients to the input point const newDestinationPoint = [ helmertParametersMatrix.get(0, 0) + helmertParametersMatrix.get(2, 0) * newSourcePoint[0] - helmertParametersMatrix.get(3, 0) * newSourcePoint[1], helmertParametersMatrix.get(1, 0) + helmertParametersMatrix.get(2, 0) * newSourcePoint[1] + helmertParametersMatrix.get(3, 0) * newSourcePoint[0] ]; return newDestinationPoint; } return interpolate; } return Helmert; }35:[] 2d:{"id":"effebbef2fb9c505","slug":null,"title":"Thin Plate Spline interpolation in Allmaps/transform - Test for Helmert transform","thumbnail":"dfd6f7d93bc79d2a1b88ec997b5b336344616a01a94a95e3285ff83a456b3897","access_level":"public","create_time":"$D2024-01-10T18:09:28.085Z","update_time":"$D2026-02-10T14:03:27.493Z","view_time":null,"first_public_time":"$D2024-01-10T18:09:28.085Z","fork":{"id":"0b57d3b587542794","slug":null,"title":"Thin Plate Spline interpolation in Allmaps/transform","owner":{"id":"8917b8f0b5bcd50a","login":"allmaps","name":"Allmaps"}},"stars":1,"starred":false,"owner":{"type":"team","id":"d151fa908c850adc","login":"mclaeysb","name":"Manuel Claeys Bouuaert","avatar_url":"https://avatars-next.observableusercontent.com/889d708b9f3a10153ae462f64a458a10cdfaea8ed74bc1b6ea12b2926f0238a8","tier":"starter_2024"},"version":1174,"publish_version":null,"role":null,"authored":false,"authors":[{"type":"individual","id":"6edbd541c72a1f56","login":"mclaeysb","name":"Manuel Claeys Bouuaert","avatar_url":"https://avatars-next.observableusercontent.com/889d708b9f3a10153ae462f64a458a10cdfaea8ed74bc1b6ea12b2926f0238a8","tier":null}],"members":[{"type":"individual","id":"6edbd541c72a1f56","login":"mclaeysb","name":"Manuel Claeys Bouuaert","avatar_url":"https://avatars-next.observableusercontent.com/889d708b9f3a10153ae462f64a458a10cdfaea8ed74bc1b6ea12b2926f0238a8","tier":null,"role":"editor"}],"comments":[],"custom_thumbnail":false,"trashed":false,"subscribed":null,"forks":0,"collection":null,"type":"notebook","body":{"title":"Thin Plate Spline interpolation in Allmaps/transform - Test for Helmert transform","cells":[{"id":0,"value":"# Thin Plate Spline interpolation in `Allmaps/transform` - Test for Helmert transform\n\nWe are in the process of implementing flexible georeferencing algorithms in Allmaps. \n\nWe have earlier built a triangulation algorithm in our [previous notebook](https://observablehq.com/d/199e169d58f0bf0d) (and are currently implementing it as the new `Allmaps/triangulate` package). Now it is time for the next step: extending the [`Allmaps/transform`](https://github.com/allmaps/allmaps/tree/main/packages/transform) package with Thin Plate Spline (TPS) interpolation, and it's more general case: Radial Basis Functions (RBF)!\n","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":1170,"value":"## *Note: This is a fork to test the implementation of a Helmert transformation*","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":992,"value":"$2f","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":40,"value":"$30","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":1141,"value":"Notice how when moving around one blue point, the cloth stays in place near the other blue points (the '**exactness**' of TPS's). And how how the cloth is smoothly warped (the '**smoothness**' of TPS's).","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"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","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":620,"value":"viewof useAffine = Inputs.toggle({\n label: \"Use Affine component (only advised for TPS)\",\n value: true\n})","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":625,"value":"viewof useAnchorPoints = Inputs.toggle({label: \"Use Anchor Points (at 10*gridSize)\"})","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":995,"value":"$31","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":987,"value":"$32","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":1113,"value":"## Ressources\n\n- Wikipedia has some great pages on this topic: \n - [Radial Basis Functions](https://en.wikipedia.org/wiki/Radial_basis_function)\n - [Radial Basis Function interpolation](https://en.wikipedia.org/wiki/Radial_basis_function_interpolation)\n - [Thin Plate Splines](https://en.wikipedia.org/wiki/Thin_plate_spline)\n- This [explanation and implementation of Radial Basis Functions in OpenCV](https://khanhha.github.io/posts/Thin-Plate-Splines-Warping/) was really helpful to us.","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":990,"value":"## Credits\n\nThis notebook started as a fork from [this one](https://observablehq.com/@esperanc/deformation-with-rbfs) by Claudio Esperança, who build the very pleasant interactive elements as well as the fundamental RBF implementation using the [RBF module by Thibaut Séguy](https://github.com/thibauts/rbf), but modified to use a selectable norm. We made some modifications to this notebook:\n\n- We swapped out the is 11 years old 'numeric' library for the **more recent 'ml-matrix' library** instead instead.\n- **We added the affine component** of the transformation, crucial to come to a smooth transformation even for a low number of ground controle points. This also allowed us to get rid of the 'anchor points' which were added far outside the plane to increase the number of points and keep the cloth in place even for a low number of points.\n- We cleaned up the naming in the code etc.\n- We added the explanatory text.\n","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":1006,"value":"# Code\n\nThis sample code is now being implemented as part of [Allmaps/transform](https://github.com/allmaps/allmaps/tree/main/packages/transform).","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":249,"value":"function computeRbf() {\n let src;\n let dst;\n if (useAnchorPoints) {\n src = controlPoints.src.concat(anchors);\n dst = controlPoints.dst.concat(anchors);\n } else {\n src = controlPoints.src;\n dst = controlPoints.dst;\n }\n return RBF(src, dst, distanceFunction.fn, euclideanNorm);\n}\n// To use this function in Observable and compute the dst location of a new src point:\n// rbf = computeRbf()\n// rbf([1, 2])","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":1162,"value":"function computeHelmert() {\n let src = controlPoints.src;\n let dst = controlPoints.dst;\n return Helmert(src, dst);\n}\n// To use this function in Observable and compute the dst location of a new src point:\n// rbf = computeRbf()\n// rbf([1, 2])","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":775,"value":"## Controle points and grid","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":166,"value":"gridSize = 400","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":157,"value":"cellSize = 20","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"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","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":750,"value":"anchors = [\n [-gridSize * 10, -gridSize * 10],\n [-gridSize * 10, gridSize * 10],\n [gridSize * 10, gridSize * 10],\n [gridSize * 10, -gridSize * 10]\n]","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"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","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"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","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"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 = \"#d2c7c8\";\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 = \"#dcded1\";\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","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":275,"value":"## Distance and norm functions","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":282,"value":"function distanceLinear(r) {\n return r;\n}\n","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":284,"value":"function distanceCubic(r) {\n return Math.pow(r, 3);\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":286,"value":"function distanceQuintic(r) {\n return Math.pow(r, 5);\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"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","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":290,"value":"function distanceGaussian(r, epsilon) {\n return Math.exp(-Math.pow(r / epsilon, 2));\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"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","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":292,"value":"function distanceMultiquadric(r, epsilon) {\n return Math.sqrt(Math.pow(r / epsilon, 2) + 1);\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"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","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"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","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":315,"value":"normFunctions = ({\n euclidean: euclideanNorm\n})","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":7,"value":"## Libraries and helper functions","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":410,"value":"mlMatrix = require(\"ml-matrix\")","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":500,"value":"// Elementwise product of two matrices, also called Hadamard product or Shur product\n// This function is not provided as part of mlMatrix\n// This is different from a classical matrix product availble in numeric as M1.mul(M2)\nfunction mulElements(M1, M2) {\n M1 = mlMatrix.Matrix.checkMatrix(M1);\n M2 = mlMatrix.Matrix.checkMatrix(M2);\n if (M1.rows != M2.rows || M1.columns != M2.columns)\n throw new Error(\"Matrices have different dimensions\");\n let M = mlMatrix.Matrix.zeros(M1.rows, M1.columns);\n for (let i = 0; i < M1.rows; i++) {\n for (let j = 0; j < M1.columns; j++) {\n M.set(i, j, M1.get(i, j) * M2.get(i, j));\n }\n }\n return M;\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":833,"value":"## Radial Basis Functions","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":13,"value":"$33","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":1149,"value":"# Helmert","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":1151,"value":"$34","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":795,"value":"## RBF of custom points\n\nThe above code works with the points displayed in the interactive window. If you want to quickly test the outcome for some custom points and distance function, you can use the example below.","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":493,"value":"computeRbfCustom = {\n let src = [\n [518, 991],\n [4345, 2357],\n [2647, 475]\n ];\n let dst = [\n [4.9516614, 52.4633102],\n [5.0480391, 52.5123762],\n [4.9702906, 52.5035815]\n ];\n return RBF(src, dst, distanceThinPlate, euclideanNorm);\n}","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":818,"value":"computeRbfCustom([100, 150])","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"}],"files":"$W35","stdlib":"1","titling":"auto"}} 2e:"$2d:body"