mutable controlePoints1 = ({
src: [
[gridSize * 0.2, gridSize * 0.1],
[gridSize * 0.8, gridSize * 0.1],
[gridSize * 0.8, gridSize * 0.7],
[gridSize * 0.2, gridSize * 0.7],
// [gridSize * 0.4, gridSize * 0.8],
[gridSize * 0.5, gridSize * 0.9]
// [gridSize * 0.6, gridSize * 1]
],
dst: [
[gridSize * 0.1, gridSize * 0.1],
[gridSize * 0.4, gridSize * 0.1],
[gridSize * 0.4, gridSize * 0.4],
[gridSize * 0.1, gridSize * 0.4]
]
})
mutable controlePoints2 = ({
src: [
[gridSize * 0.2, gridSize * 0.3],
[gridSize * 0.8, gridSize * 0.3],
[gridSize * 0.8, gridSize * 0.9],
[gridSize * 0.2, gridSize * 0.9],
// [gridSize * 0.4, gridSize * 0.0],
[gridSize * 0.5, gridSize * 0.1]
// [gridSize * 0.6, gridSize * 0.2]
],
dst: [
[gridSize * 0.1, gridSize * 0.6],
[gridSize * 0.4, gridSize * 0.6],
[gridSize * 0.4, gridSize * 0.9],
[gridSize * 0.1, gridSize * 0.9]
]
})
function Transformation(
sourcePoints1,
destinationPoints1,
sourcePoints2,
destinationPoints2,
sourceStaplePoints1,
sourceStaplePoints2,
kernelFunction,
normFunction,
epsilon1,
epsilon2
) {
const pointCount1 = sourcePoints1.length;
const pointCount2 = sourcePoints2.length;
const staplePointCount =
enforceStaplePoints && sourceStaplePoints1 ? sourceStaplePoints1.length : 0;
let destinationPoints;
if (transformationType == "thinPlateSpline") {
destinationPoints = [
...destinationPoints1,
[0, 0],
[0, 0],
[0, 0],
...destinationPoints2,
[0, 0],
[0, 0],
[0, 0]
];
} else {
destinationPoints = [...destinationPoints1, ...destinationPoints2];
}
if (staplePointCount) {
destinationPoints = [
...destinationPoints,
...sourceStaplePoints1.map((point) => [0, 0])
];
}
const destinationPointsMatrices = [
mlMatrix.Matrix.columnVector(destinationPoints.map((value) => value[0])),
mlMatrix.Matrix.columnVector(destinationPoints.map((value) => value[1]))
];
let kernelsMatrix1;
let affineCoefsMatrix1;
let staplesKernelsMatrix1;
let staplesAffineCoefsMatrix1;
let kernelsMatrix2;
let affineCoefsMatrix2;
let staplesKernelsMatrix2;
let staplesAffineCoefsMatrix2;
// Pre-compute kernelsMatrix: fill it with the point to point distances between all control points
kernelsMatrix1 = mlMatrix.Matrix.zeros(pointCount1, pointCount1);
for (let i = 0; i < pointCount1; i++) {
for (let j = 0; j < pointCount1; j++) {
kernelsMatrix1.set(
i,
j,
normFunction(sourcePoints1[i], sourcePoints1[j])
);
}
}
// If it's not provided, and if it's an input to the kernelFunction,
// compute epsilon as the average distance between the control points
if (epsilon1 === undefined) {
epsilon1 = kernelsMatrix1.sum() / (Math.pow(pointCount1, 2) - pointCount1);
}
// Finish the computation of kernelsMatrix by applying the requested kernel function
for (let i = 0; i < pointCount1; i++) {
for (let j = 0; j < pointCount1; j++) {
kernelsMatrix1.set(
i,
j,
kernelFunction(kernelsMatrix1.get(i, j), { epsilon: epsilon1 })
);
}
}
////
// Pre-compute kernelsMatrix: fill it with the point to point distances between all control points
kernelsMatrix2 = mlMatrix.Matrix.zeros(pointCount2, pointCount2);
for (let i = 0; i < pointCount2; i++) {
for (let j = 0; j < pointCount2; j++) {
kernelsMatrix2.set(
i,
j,
normFunction(sourcePoints2[i], sourcePoints2[j])
);
}
}
// If it's not provided, and if it's an input to the kernelFunction,
// compute epsilon as the average distance between the control points
if (epsilon2 === undefined) {
epsilon2 = kernelsMatrix2.sum() / (Math.pow(pointCount2, 2) - pointCount2);
}
// Finish the computation of kernelsMatrix by applying the requested kernel function
for (let i = 0; i < pointCount2; i++) {
for (let j = 0; j < pointCount2; j++) {
kernelsMatrix2.set(
i,
j,
kernelFunction(kernelsMatrix2.get(i, j), { epsilon: epsilon2 })
);
}
}
////
// Define the affine transformation
affineCoefsMatrix1 = mlMatrix.Matrix.zeros(pointCount1, 3);
// Construct Nx3 Matrix affineCoefsMatrix
// 1 x0 y0
// 1 x1 y1
// 1 x2 y2
// ...
for (let i = 0; i < pointCount1; i++) {
affineCoefsMatrix1.set(i, 0, 1);
affineCoefsMatrix1.set(i, 1, sourcePoints1[i][0]);
affineCoefsMatrix1.set(i, 2, sourcePoints1[i][1]);
}
////
// Define the affine transformation
affineCoefsMatrix2 = mlMatrix.Matrix.zeros(pointCount2, 3);
// Construct Nx3 Matrix affineCoefsMatrix
// 1 x0 y0
// 1 x1 y1
// 1 x2 y2
// ...
for (let i = 0; i < pointCount2; i++) {
affineCoefsMatrix2.set(i, 0, 1);
affineCoefsMatrix2.set(i, 1, sourcePoints2[i][0]);
affineCoefsMatrix2.set(i, 2, sourcePoints2[i][1]);
}
////
if (staplePointCount) {
// Pre-compute staplesKernelsMatrix: fill it with the point to point distances between all control points
staplesKernelsMatrix1 = mlMatrix.Matrix.zeros(
staplePointCount,
pointCount1
);
for (let i = 0; i < staplePointCount; i++) {
for (let j = 0; j < pointCount1; j++) {
staplesKernelsMatrix1.set(
i,
j,
normFunction(sourceStaplePoints1[i], sourcePoints1[j])
);
}
}
// Finish the computation of staplesKernelsMatrix1 by applying the requested kernel function
for (let i = 0; i < staplePointCount; i++) {
for (let j = 0; j < pointCount1; j++) {
staplesKernelsMatrix1.set(
i,
j,
kernelFunction(staplesKernelsMatrix1.get(i, j), { epsilon: epsilon1 })
);
}
}
// Define the affine transformation
staplesAffineCoefsMatrix1 = mlMatrix.Matrix.zeros(staplePointCount, 3);
// Construct Nx3 Matrix affineCoefsMatrix
// 1 x0 y0
// 1 x1 y1
// 1 x2 y2
// ...
for (let i = 0; i < staplePointCount; i++) {
staplesAffineCoefsMatrix1.set(i, 0, 1);
staplesAffineCoefsMatrix1.set(i, 1, sourceStaplePoints1[i][0]);
staplesAffineCoefsMatrix1.set(i, 2, sourceStaplePoints1[i][1]);
}
////
// Pre-compute staplesKernelsMatrix: fill it with the point to point distances between all control points
staplesKernelsMatrix2 = mlMatrix.Matrix.zeros(
staplePointCount,
pointCount2
);
for (let i = 0; i < staplePointCount; i++) {
for (let j = 0; j < pointCount2; j++) {
staplesKernelsMatrix2.set(
i,
j,
normFunction(sourceStaplePoints2[i], sourcePoints2[j])
);
}
}
// Finish the computation of staplesKernelsMatrix1 by applying the requested kernel function
for (let i = 0; i < staplePointCount; i++) {
for (let j = 0; j < pointCount2; j++) {
staplesKernelsMatrix2.set(
i,
j,
kernelFunction(staplesKernelsMatrix2.get(i, j), { epsilon: epsilon2 })
);
}
}
// Define the affine transformation
staplesAffineCoefsMatrix2 = mlMatrix.Matrix.zeros(staplePointCount, 3);
// Construct Nx3 Matrix affineCoefsMatrix
// 1 x0 y0
// 1 x1 y1
// 1 x2 y2
// ...
for (let i = 0; i < staplePointCount; i++) {
staplesAffineCoefsMatrix2.set(i, 0, 1);
staplesAffineCoefsMatrix2.set(i, 1, sourceStaplePoints2[i][0]);
staplesAffineCoefsMatrix2.set(i, 2, sourceStaplePoints2[i][1]);
}
}
////
function addBlockMatrix(matrix, blockMatrix, atRow, atColumn) {
for (let i = 0; i < blockMatrix.rows; i++) {
for (let j = 0; j < blockMatrix.columns; j++) {
matrix.set(atRow + i, atColumn + j, blockMatrix.get(i, j));
}
}
return matrix;
}
// Combine kernelsMatrix and affineCoefsMatrix into new matrix transformationMatrix
let transformationMatrix;
if (transformationType == "thinPlateSpline") {
transformationMatrix = mlMatrix.Matrix.zeros(
pointCount1 + 3 + pointCount2 + 3 + staplePointCount,
pointCount1 + 3 + pointCount2 + 3
);
transformationMatrix = addBlockMatrix(
transformationMatrix,
kernelsMatrix1,
0,
0
);
transformationMatrix = addBlockMatrix(
transformationMatrix,
affineCoefsMatrix1.transpose(),
pointCount1,
0
);
transformationMatrix = addBlockMatrix(
transformationMatrix,
affineCoefsMatrix1,
0,
pointCount1
);
transformationMatrix = addBlockMatrix(
transformationMatrix,
kernelsMatrix2,
pointCount1 + 3,
pointCount1 + 3
);
transformationMatrix = addBlockMatrix(
transformationMatrix,
affineCoefsMatrix2.transpose(),
pointCount1 + 3 + pointCount2,
pointCount1 + 3
);
transformationMatrix = addBlockMatrix(
transformationMatrix,
affineCoefsMatrix2,
pointCount1 + 3,
pointCount1 + 3 + pointCount2
);
if (staplePointCount) {
transformationMatrix = addBlockMatrix(
transformationMatrix,
staplesKernelsMatrix1,
pointCount1 + 3 + pointCount2 + 3,
0
);
transformationMatrix = addBlockMatrix(
transformationMatrix,
staplesAffineCoefsMatrix1,
pointCount1 + 3 + pointCount2 + 3,
pointCount1
);
transformationMatrix = addBlockMatrix(
transformationMatrix,
mlMatrix.Matrix.multiply(staplesKernelsMatrix2, -1),
pointCount1 + 3 + pointCount2 + 3,
pointCount1 + 3
);
transformationMatrix = addBlockMatrix(
transformationMatrix,
mlMatrix.Matrix.multiply(staplesAffineCoefsMatrix2, -1),
pointCount1 + 3 + pointCount2 + 3,
pointCount1 + 3 + pointCount2
);
}
} else {
transformationMatrix = mlMatrix.Matrix.zeros(
pointCount1 + pointCount2 + staplePointCount,
3 + 3
);
transformationMatrix = addBlockMatrix(
transformationMatrix,
affineCoefsMatrix1,
0,
0
);
transformationMatrix = addBlockMatrix(
transformationMatrix,
affineCoefsMatrix2,
pointCount1,
3
);
if (staplePointCount) {
transformationMatrix = addBlockMatrix(
transformationMatrix,
staplesAffineCoefsMatrix1,
pointCount1 + pointCount2,
0
);
transformationMatrix = addBlockMatrix(
transformationMatrix,
mlMatrix.Matrix.multiply(staplesAffineCoefsMatrix2, -1),
pointCount1 + pointCount2,
3
);
}
}
// Compute basis functions weights and the affine parameters by solving the linear system of equations for each component
// Note: the same transformationMatrix is used for both solutions
const inverseKernelsAndAffineCoefsMatrix =
mlMatrix.inverse(transformationMatrix);
const weightsMatrices = [
inverseKernelsAndAffineCoefsMatrix.mmul(destinationPointsMatrices[0]),
inverseKernelsAndAffineCoefsMatrix.mmul(destinationPointsMatrices[1])
];
// Store rbf and affine parts of the weights more as arrays for more efficient access on evaluation
let rbfWeights1;
let affineWeights1;
let rbfWeights2;
let affineWeights2;
if (transformationType == "thinPlateSpline") {
rbfWeights1 = weightsMatrices.map((matrix) =>
matrix.selection([...Array(pointCount1).keys()], [0]).to1DArray()
);
affineWeights1 = weightsMatrices.map((matrix) =>
matrix
.selection(
[0, 1, 2].map((n) => n + pointCount1),
[0]
)
.to1DArray()
);
rbfWeights2 = weightsMatrices.map((matrix) =>
matrix
.selection(
[
...Array(pointCount2)
.keys()
.map((n) => n + pointCount1 + 3)
],
[0]
)
.to1DArray()
);
affineWeights2 = weightsMatrices.map((matrix) =>
matrix
.selection(
[0, 1, 2].map((n) => n + pointCount1 + 3 + pointCount2),
[0]
)
.to1DArray()
);
} else {
affineWeights1 = weightsMatrices.map((matrix) =>
matrix.selection([0, 1, 2], [0]).to1DArray()
);
affineWeights2 = weightsMatrices.map((matrix) =>
matrix.selection([3, 4, 5], [0]).to1DArray()
);
}
return {
destinationPointsMatrices,
transformationMatrix,
weightsMatrices,
rbfWeights1,
affineWeights1,
rbfWeights2,
affineWeights2,
epsilon1,
epsilon2
};
}
solveTransformation = () =>
Transformation(
sourcePoints1NonStaple,
destinationPoints1NonStaple,
sourcePoints2NonStaple,
destinationPoints2NonStaple,
sourcePoints1Staple,
sourcePoints2Staple,
kernelThinPlate,
euclideanNorm
)
// pt(solveTransformation().destinationPointsMatrices[0], { maxrows: 20 })
pt(solveTransformation().transformationMatrix, { maxrows: 20 })
function evaluateTransformation(
sourcePoints,
destinationPoints,
kernelFunction,
normFunction,
epsilon,
rbfWeights,
affineWeights,
newSourcePoint
) {
// Evaluate the transformation function at a new point
const pointCount = sourcePoints.length;
// Compute the distances of that point to all control points
const newDistances = sourcePoints.map((sourcePoint) =>
normFunction(newSourcePoint, sourcePoint)
);
// Sum the weighted contributions of the input point
const newDestinationPoint = [0, 0];
for (let i = 0; i < 2; i++) {
// Apply the weights to the new distances
if (transformationType == "thinPlateSpline") {
newDestinationPoint[i] += newDistances.reduce(
(sum, dist, index) =>
sum +
kernelFunction(dist, { epsilon: epsilon }) * rbfWeights[i][index],
0
);
}
// Add the affine part
newDestinationPoint[i] +=
affineWeights[i][0] +
affineWeights[i][1] * newSourcePoint[0] +
affineWeights[i][2] * newSourcePoint[1];
}
return newDestinationPoint;
}
function evaluateTransformation1(newSourcePoint) {
const solvedTransformation = solveTransformation();
return evaluateTransformation(
sourcePoints1NonStaple,
destinationPoints1NonStaple,
kernelThinPlate,
euclideanNorm,
solvedTransformation.epsilon1,
solvedTransformation.rbfWeights1,
solvedTransformation.affineWeights1,
newSourcePoint
);
}
function evaluateTransformation1Grid(newSourcePoint) {
let solvedTransformation;
if (intermediateOrFinal == "intermediate") {
solvedTransformation = solveTransformation();
return evaluateTransformation(
sourcePoints1NonStaple,
destinationPoints1NonStaple,
kernelThinPlate,
euclideanNorm,
solvedTransformation.epsilon1,
solvedTransformation.rbfWeights1,
solvedTransformation.affineWeights1,
newSourcePoint
);
} else {
solvedTransformation = Transformation(
sourcePoints1,
destinationPoints1,
sourcePoints2,
destinationPoints2,
undefined,
undefined,
kernelThinPlate,
euclideanNorm
);
return evaluateTransformation(
sourcePoints1,
destinationPoints1,
kernelThinPlate,
euclideanNorm,
solvedTransformation.epsilon1,
solvedTransformation.rbfWeights1,
solvedTransformation.affineWeights1,
newSourcePoint
);
}
}
function evaluateTransformation2(newSourcePoint) {
const solvedTransformation = solveTransformation();
return evaluateTransformation(
sourcePoints2NonStaple,
destinationPoints2NonStaple,
kernelThinPlate,
euclideanNorm,
solvedTransformation.epsilon2,
solvedTransformation.rbfWeights2,
solvedTransformation.affineWeights2,
newSourcePoint
);
}
evaluateTransformation1([100, 180])
evaluateTransformation2([100, 20])
function evaluateTransformation2Grid(newSourcePoint) {
let solvedTransformation;
if (intermediateOrFinal == "intermediate") {
solvedTransformation = solveTransformation();
return evaluateTransformation(
sourcePoints2NonStaple,
destinationPoints2NonStaple,
kernelThinPlate,
euclideanNorm,
solvedTransformation.epsilon2,
solvedTransformation.rbfWeights2,
solvedTransformation.affineWeights2,
newSourcePoint
);
} else {
solvedTransformation = Transformation(
sourcePoints1,
destinationPoints1,
sourcePoints2,
destinationPoints2,
undefined,
undefined,
kernelThinPlate,
euclideanNorm
);
return evaluateTransformation(
sourcePoints2,
destinationPoints2,
kernelThinPlate,
euclideanNorm,
solvedTransformation.epsilon2,
solvedTransformation.rbfWeights2,
solvedTransformation.affineWeights2,
newSourcePoint
);
}
}
Purpose-built for displays of data
Observable is your go-to platform for exploring data and creating expressive data visualizations. Use reactive JavaScript notebooks for prototyping and a collaborative canvas for visual data exploration and dashboard creation.
