I'm into maps, maths and geometry - and many other aspects of life which I'm glad I can commit to without using git.
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;
}
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;
}
computeRbfCustom = {
let src = [
[518, 991],
[4345, 2357],
[2647, 475]
];
let dst = [
[4.9516614, 52.4633102],
[5.0480391, 52.5123762],
[4.9702906, 52.5035815]
];
return RBF(src, dst, distanceThinPlate, euclideanNorm);
}
computeRbfCustom([100, 150])
One platform to build and deploy the best data apps
Experiment and prototype by building visualizations in live JavaScript notebooks. Collaborate with your team and decide which concepts to build out.
Use Observable Framework to build data apps locally. Use data loaders to build in any language or library, including Python, SQL, and R.
Seamlessly deploy to Observable. Test before you ship, use automatic deploy-on-commit, and ensure your projects are always up-to-date.