@allmaps/transform - npm Package Compare versions

Comparing version 1.0.0-beta.14 to 1.0.0-beta.15

dist/index.d.ts

@@ -1,4 +0,9 @@
-export * from './transformer.js';
-export * from './geojson.js';
+import GCPTransformer from './transformer.js';
+import HelmertGCPTransformer from './transformers/helmert-transformer.js';
+import PolynomialGCPTransformer from './transformers/polynomial-transformer.js';
+import ProjectiveGCPTransformer from './transformers/projective-transformer.js';
+import RadialBasisFunctionGCPTransformer from './transformers/radial-basis-function-transformer.js';
+import { thinPlateKernel } from './shared/kernel-functions.js';
+export { GCPTransformer, HelmertGCPTransformer, PolynomialGCPTransformer, ProjectiveGCPTransformer, RadialBasisFunctionGCPTransformer };
+export { thinPlateKernel };
 export * from './shared/types.js';
-export type { GCPTransformInfo } from './gdaltransform.js';

dist/index.js

@@ -1,3 +0,11 @@
-export * from './transformer.js';
-export * from './geojson.js';
+import GCPTransformer from './transformer.js';
+import HelmertGCPTransformer from './transformers/helmert-transformer.js';
+import PolynomialGCPTransformer from './transformers/polynomial-transformer.js';
+import ProjectiveGCPTransformer from './transformers/projective-transformer.js';
+import RadialBasisFunctionGCPTransformer from './transformers/radial-basis-function-transformer.js';
+import { thinPlateKernel } from './shared/kernel-functions.js';
+// TODO: consider only exporting GCPTransformer
+export { GCPTransformer, HelmertGCPTransformer, PolynomialGCPTransformer, ProjectiveGCPTransformer, RadialBasisFunctionGCPTransformer };
+// TODO: can this be removed?
+export { thinPlateKernel };
 export * from './shared/types.js';

dist/shared/gdaltransform.js

@@ -96,3 +96,3 @@ const HUGE_VAL = Number.POSITIVE_INFINITY;
     // let nCRSresult
-    let sPoints = new ControlPoints();
+    const sPoints = new ControlPoints();
     let x1Sum = 0;
@@ -147,7 +147,7 @@ let y1Sum = 0;
         /* -------------------------------------------------------------------- */
-        let padfGeoX = [];
-        let padfGeoY = [];
-        let padfRasterX = [];
-        let padfRasterY = [];
-        let panStatus = [];
+        const padfGeoX = [];
+        const padfGeoY = [];
+        const padfRasterX = [];
+        const padfRasterY = [];
+        const panStatus = [];
         for (let iGCP = 0; iGCP < nGCPCount; iGCP++) {
@@ -212,4 +212,4 @@ panStatus[iGCP] = 1;
     const m = new MATRIX();
-    let a = [];
-    let b = [];
+    const a = [];
+    const b = [];
     let numactive = 0; /* NUMBER OF ACTIVE CONTROL POINTS */
@@ -320,3 +320,3 @@ let i = 0;
     for (let i = 1; i <= m.n; i++) {
-        let j = i;
+        const j = i;
         /* find row with largest magnitude value for pivot value */
@@ -323,0 +323,0 @@ let pivot = m.getM(i, j); /* ACTUAL VALUE OF THE LARGEST PIVOT CANDIDATE */

dist/shared/geojson.d.ts

@@ -8,1 +8,3 @@ import type { Position, GeoJSONPoint, GeoJSONLineString, GeoJSONPolygon, GCPTransformerInterface, OptionalTransformOptions } from './types.js';
 export declare function fromGeoJSONPolygon(transformer: GCPTransformerInterface, polygon: GeoJSONPolygon, options?: OptionalTransformOptions): Position[];
+export declare function toGeoPolygon(transformer: GCPTransformerInterface, points: Position[], options?: OptionalTransformOptions): Position[];
+export declare function toResourcePolygon(transformer: GCPTransformerInterface, points: Position[], options?: OptionalTransformOptions): Position[];

dist/shared/geojson.js

@@ -6,3 +6,3 @@ // TODO: consider implementing these functions in this module instead of using dependencies
 function mergeDefaultOptions(options) {
-    let mergedOptions = {
+    const mergedOptions = {
         close: false,
@@ -43,3 +43,3 @@ maxOffsetRatio: 0,
 export function toGeoJSONPoint(transformer, point) {
-    return makeGeoJSONPoint(transformer.toWorld(point));
+    return makeGeoJSONPoint(transformer.toGeo(point));
 }
@@ -51,27 +51,27 @@ export function toGeoJSONLineString(transformer, points, options) {
     }
-    const imageWorldPoints = points.map((point) => ({
-        image: point,
-        world: transformer.toWorld(point)
+    const resourceGeoPoints = points.map((point) => ({
+        resource: point,
+        geo: transformer.toGeo(point)
     }));
-    const segments = getSegments(imageWorldPoints);
+    const segments = getSegments(resourceGeoPoints);
     const segmentsWithMidpoints = recursivelyAddWorldMidpoints(transformer, segments, mergedOptions);
     return makeGeoJSONLineString([
-        segmentsWithMidpoints[0].from.world,
-        ...segmentsWithMidpoints.map((segment) => segment.to.world)
+        segmentsWithMidpoints[0].from.geo,
+        ...segmentsWithMidpoints.map((segment) => segment.to.geo)
     ]);
 }
 export function toGeoJSONPolygon(transformer, points, options) {
-    const mergedOtions = mergeDefaultOptions(options);
+    const mergedOptions = mergeDefaultOptions(options);
     if (!Array.isArray(points) || points.length < 3) {
         throw new Error('Polygon should contain at least 3 points');
     }
-    const imageWorldPoints = points.map((point) => ({
-        image: point,
-        world: transformer.toWorld(point)
+    const resourceGeoPoints = points.map((point) => ({
+        resource: point,
+        geo: transformer.toGeo(point)
     }));
-    const segments = getSegments(imageWorldPoints, true);
-    const segmentsWithMidpoints = recursivelyAddWorldMidpoints(transformer, segments, mergedOtions);
+    const segments = getSegments(resourceGeoPoints, true);
+    const segmentsWithMidpoints = recursivelyAddWorldMidpoints(transformer, segments, mergedOptions);
     return makeGeoJSONPolygon([
-        segmentsWithMidpoints[0].from.world,
-        ...segmentsWithMidpoints.map((segment) => segment.to.world)
+        segmentsWithMidpoints[0].from.geo,
+        ...segmentsWithMidpoints.map((segment) => segment.to.geo)
     ]);
@@ -83,3 +83,3 @@ }
 export function fromGeoJSONLineString(transformer, lineString, options) {
-    const mergedOtions = mergeDefaultOptions(options);
+    const mergedOptions = mergeDefaultOptions(options);
     const coordinates = lineString.coordinates;
@@ -89,15 +89,15 @@ if (!Array.isArray(coordinates) || coordinates.length < 2) {
     }
-    const imageWorldPoints = coordinates.map((coordinate) => ({
-        image: transformer.toResource(coordinate),
-        world: coordinate
+    const resourceGeoPoints = coordinates.map((coordinate) => ({
+        resource: transformer.toResource(coordinate),
+        geo: coordinate
     }));
-    const segments = getSegments(imageWorldPoints);
-    const segmentsWithMidpoints = recursivelyAddImageMidpoints(transformer, segments, mergedOtions);
+    const segments = getSegments(resourceGeoPoints);
+    const segmentsWithMidpoints = recursivelyAddResourceMidpoints(transformer, segments, mergedOptions);
     return [
-        segmentsWithMidpoints[0].from.image,
-        ...segmentsWithMidpoints.map((segment) => segment.to.image)
+        segmentsWithMidpoints[0].from.resource,
+        ...segmentsWithMidpoints.map((segment) => segment.to.resource)
     ];
 }
 export function fromGeoJSONPolygon(transformer, polygon, options) {
-    const mergedOtions = mergeDefaultOptions(options);
+    const mergedOptions = mergeDefaultOptions(options);
     const coordinates = polygon.coordinates[0];
@@ -107,16 +107,48 @@ if (!Array.isArray(coordinates) || coordinates.length < 4) {
     }
-    const imageWorldPoints = coordinates.map((coordinates) => ({
-        image: transformer.toResource(coordinates),
-        world: coordinates
+    const resourceGeoPoints = coordinates.map((coordinates) => ({
+        resource: transformer.toResource(coordinates),
+        geo: coordinates
     }));
-    const segments = getSegments(imageWorldPoints);
-    const segmentsWithMidpoints = recursivelyAddImageMidpoints(transformer, segments, mergedOtions);
+    const segments = getSegments(resourceGeoPoints);
+    const segmentsWithMidpoints = recursivelyAddResourceMidpoints(transformer, segments, mergedOptions);
     return [
-        segmentsWithMidpoints[0].from.image,
-        ...segmentsWithMidpoints.map((segment) => segment.to.image)
+        segmentsWithMidpoints[0].from.resource,
+        ...segmentsWithMidpoints.map((segment) => segment.to.resource)
     ];
 }
+export function toGeoPolygon(transformer, points, options) {
+    const mergedOptions = mergeDefaultOptions(options);
+    if (!Array.isArray(points) || points.length < 3) {
+        throw new Error('Polygon should contain at least 3 points');
+    }
+    const resourceGeoPoints = points.map((point) => ({
+        resource: point,
+        geo: transformer.toGeo(point)
+    }));
+    const segments = getSegments(resourceGeoPoints, false);
+    const segmentsWithMidpoints = recursivelyAddWorldMidpoints2(transformer, segments, mergedOptions);
+    return [
+        segmentsWithMidpoints[0].from.geo,
+        ...segmentsWithMidpoints.map((segment) => segment.to.geo)
+    ];
+}
+export function toResourcePolygon(transformer, points, options) {
+    const mergedOptions = mergeDefaultOptions(options);
+    if (!Array.isArray(points) || points.length < 3) {
+        throw new Error('Polygon should contain at least 3 points');
+    }
+    const resourceGeoPoints = points.map((point) => ({
+        resource: transformer.toResource(point),
+        geo: point
+    }));
+    const segments = getSegments(resourceGeoPoints, false);
+    const segmentsWithMidpoints = recursivelyAddResourceMidpoints2(transformer, segments, mergedOptions);
+    return [
+        segmentsWithMidpoints[0].from.resource,
+        ...segmentsWithMidpoints.map((segment) => segment.to.resource)
+    ];
+}
 function getSegments(points, close = false) {
     const segmentCount = points.length - (close ? 0 : 1);
-    let segments = [];
+    const segments = [];
     for (let index = 0; index < segmentCount; index++) {
@@ -138,3 +170,3 @@ segments.push({
 }
-function recursivelyAddImageMidpoints(transformer, segments, options) {
+function recursivelyAddResourceMidpoints(transformer, segments, options) {
     if (options.maxDepth <= 0 || options.maxOffsetRatio <= 0) {
@@ -144,10 +176,26 @@ return segments;
     return segments
-        .map((segment) => addImageMidpoints(transformer, segment, options, 0))
+        .map((segment) => addResourceMidpoints(transformer, segment, options, 0))
         .flat(1);
 }
+function recursivelyAddWorldMidpoints2(transformer, segments, options) {
+    if (options.maxDepth <= 0 || options.maxOffsetRatio <= 0) {
+        return segments;
+    }
+    return segments
+        .map((segment) => addWorldMidpoints2(transformer, segment, options, 0))
+        .flat(1);
+}
+function recursivelyAddResourceMidpoints2(transformer, segments, options) {
+    if (options.maxDepth <= 0 || options.maxOffsetRatio <= 0) {
+        return segments;
+    }
+    return segments
+        .map((segment) => addResourceMidpoints2(transformer, segment, options, 0))
+        .flat(1);
+}
 function addWorldMidpoints(transformer, segment, options, depth) {
-    const imageMidpoint = getImageMidpoint(segment.from.image, segment.to.image);
-    const segmentWorldMidpoint = getWorldMidpoint(makeGeoJSONPoint(segment.from.world), makeGeoJSONPoint(segment.to.world));
-    const actualWorldMidpoint = makeGeoJSONPoint(transformer.toWorld(imageMidpoint));
-    const distanceSegment = getWorldDistance(makeGeoJSONPoint(segment.from.world), makeGeoJSONPoint(segment.to.world));
+    const imageMidpoint = getMidpoint(segment.from.resource, segment.to.resource);
+    const segmentWorldMidpoint = getWorldMidpoint(makeGeoJSONPoint(segment.from.geo), makeGeoJSONPoint(segment.to.geo));
+    const actualWorldMidpoint = makeGeoJSONPoint(transformer.toGeo(imageMidpoint));
+    const distanceSegment = getWorldDistance(makeGeoJSONPoint(segment.from.geo), makeGeoJSONPoint(segment.to.geo));
     const distanceMidpoints = getWorldDistance(segmentWorldMidpoint, actualWorldMidpoint);
@@ -160,4 +208,4 @@ if (depth < options.maxDepth &&
         const newSegmentMidpoint = {
-            image: imageMidpoint,
-            world: actualWorldMidpoint.coordinates
+            resource: imageMidpoint,
+            geo: actualWorldMidpoint.coordinates
         };
@@ -173,8 +221,8 @@ return [
 }
-function addImageMidpoints(transformer, segment, options, depth) {
-    const worldMidpoint = getWorldMidpoint(segment.from.world, segment.to.world);
-    const segmentImageMidpoint = getImageMidpoint(segment.from.image, segment.to.image);
+function addResourceMidpoints(transformer, segment, options, depth) {
+    const worldMidpoint = getWorldMidpoint(segment.from.geo, segment.to.geo);
+    const segmentImageMidpoint = getMidpoint(segment.from.resource, segment.to.resource);
     const actualImageMidpoint = transformer.toResource(worldMidpoint.geometry.coordinates);
-    const distanceSegment = getImageDistance(segment.from.image, segment.to.image);
-    const distanceMidpoints = getImageDistance(segmentImageMidpoint, actualImageMidpoint);
+    const distanceSegment = getDistance(segment.from.resource, segment.to.resource);
+    const distanceMidpoints = getDistance(segmentImageMidpoint, actualImageMidpoint);
     if (depth < options.maxDepth &&
@@ -186,8 +234,8 @@ (options.maxOffsetRatio
         const newSegmentMidpoint = {
-            image: actualImageMidpoint,
-            world: worldMidpoint.geometry.coordinates
+            resource: actualImageMidpoint,
+            geo: worldMidpoint.geometry.coordinates
         };
         return [
-            addImageMidpoints(transformer, { from: segment.from, to: newSegmentMidpoint }, options, depth + 1),
-            addWorldMidpoints(transformer, { from: newSegmentMidpoint, to: segment.to }, options, depth + 1)
+            addResourceMidpoints(transformer, { from: segment.from, to: newSegmentMidpoint }, options, depth + 1),
+            addResourceMidpoints(transformer, { from: newSegmentMidpoint, to: segment.to }, options, depth + 1)
         ].flat(1);
@@ -199,3 +247,51 @@ }
 }
-function getImageMidpoint(point1, point2) {
+function addWorldMidpoints2(transformer, segment, options, depth) {
+    const imageMidpoint = getMidpoint(segment.from.resource, segment.to.resource);
+    const segmentWorldMidpoint = getMidpoint(segment.from.geo, segment.to.geo);
+    const actualWorldMidpoint = transformer.toGeo(imageMidpoint);
+    const distanceSegment = getDistance(segment.from.geo, segment.to.geo);
+    const distanceMidpoints = getDistance(segmentWorldMidpoint, actualWorldMidpoint);
+    if (depth < options.maxDepth &&
+        (options.maxOffsetRatio
+            ? distanceMidpoints / distanceSegment > options.maxOffsetRatio
+            : false) &&
+        distanceSegment > 0) {
+        const newSegmentMidpoint = {
+            resource: imageMidpoint,
+            geo: actualWorldMidpoint
+        };
+        return [
+            addWorldMidpoints2(transformer, { from: segment.from, to: newSegmentMidpoint }, options, depth + 1),
+            addWorldMidpoints2(transformer, { from: newSegmentMidpoint, to: segment.to }, options, depth + 1)
+        ].flat(1);
+    }
+    else {
+        return segment;
+    }
+}
+function addResourceMidpoints2(transformer, segment, options, depth) {
+    const worldMidpoint = getMidpoint(segment.from.geo, segment.to.geo);
+    const segmentImageMidpoint = getMidpoint(segment.from.resource, segment.to.resource);
+    const actualImageMidpoint = transformer.toResource(worldMidpoint);
+    const distanceSegment = getDistance(segment.from.resource, segment.to.resource);
+    const distanceMidpoints = getDistance(segmentImageMidpoint, actualImageMidpoint);
+    if (depth < options.maxDepth &&
+        (options.maxOffsetRatio
+            ? distanceMidpoints / distanceSegment > options.maxOffsetRatio
+            : false) &&
+        distanceSegment > 0) {
+        const newSegmentMidpoint = {
+            resource: actualImageMidpoint,
+            geo: worldMidpoint
+        };
+        return [
+            addResourceMidpoints2(transformer, { from: segment.from, to: newSegmentMidpoint }, options, depth + 1),
+            addResourceMidpoints2(transformer, { from: newSegmentMidpoint, to: segment.to }, options, depth + 1)
+        ].flat(1);
+    }
+    else {
+        return segment;
+    }
+}
+function getMidpoint(point1, point2) {
     return [
@@ -206,4 +302,4 @@ (point2[0] - point1[0]) / 2 + point1[0],
 }
-function getImageDistance(from, to) {
+function getDistance(from, to) {
     return Math.sqrt((to[0] - from[0]) ** 2 + (to[1] - from[1]) ** 2);
 }

dist/shared/radial-basis-function.d.ts

 import { Matrix } from 'ml-matrix';
-import type { DistanceFunction, NormFunction, Position } from './types.js';
-/**
- * This is essentially the rbf module by Thibaut Séguy, but modified to use a selectable norm
- * See https://github.com/thibauts/rbf
- **/
+import type { KernelFunction, NormFunction, Position } from './types.js';
 export default class RBF {
-    points: Position[];
-    values: Position[];
-    distanceFunction: DistanceFunction;
+    sourcePoints: Position[];
+    destinationPoints: Position[];
+    kernelFunction: KernelFunction;
     normFunction: NormFunction;
-    epsilon: number;
-    weights?: [Matrix, Matrix];
-    constructor(points: Position[], values: Position[], distanceFunction: DistanceFunction, normFunction: NormFunction, epsilon?: number);
-    interpolant(point: Position): Position;
+    weightsMatrices: [Matrix, Matrix];
+    nPoints: number;
+    epsilon?: number;
+    constructor(sourcePoints: Position[], destinationPoints: Position[], kernelFunction: KernelFunction, normFunction: NormFunction, epsilon?: number);
+    interpolant(newSourcePoint: Position): Position;
 }

dist/shared/radial-basis-function.js

@@ -1,82 +0,102 @@
-import { Matrix, solve } from 'ml-matrix';
-function sumMatrix(matrix) {
-    let sum = 0;
-    for (let j = 0; j < matrix.rows; j++) {
-        for (let i = 0; i < matrix.columns; i++) {
-            sum += matrix.get(j, i);
-        }
-    }
-    return sum;
-}
-function mulitplyMatricesValues(matrix1, matrix2) {
-    let mul = Matrix.zeros(matrix1.rows, matrix1.columns);
-    for (let j = 0; j < matrix1.rows; j++) {
-        for (let i = 0; i < matrix1.columns; i++) {
-            mul.set(j, i, matrix1.get(j, i) * matrix2.get(j, i));
-        }
-    }
-    return mul;
-}
-/**
- * This is essentially the rbf module by Thibaut Séguy, but modified to use a selectable norm
- * See https://github.com/thibauts/rbf
- **/
+import { Matrix, inverse } from 'ml-matrix';
+import { multiplyMatricesElementwise } from './matrix.js';
 export default class RBF {
-    constructor(points, values, distanceFunction, normFunction, epsilon) {
-        this.points = points;
-        this.values = values;
-        this.distanceFunction = distanceFunction;
+    constructor(sourcePoints, destinationPoints, kernelFunction, normFunction, epsilon) {
+        this.sourcePoints = sourcePoints;
+        this.destinationPoints = destinationPoints;
+        this.kernelFunction = kernelFunction;
         this.normFunction = normFunction;
-        // `identity` here serves as a shorthand to allocate
-        // the matrix nested array.
-        const M = Matrix.eye(this.points.length);
-        // First compute the point to point distance matrix
-        // to allow computing epsilon if it's not provided
-        for (let j = 0; j < this.points.length; j++) {
-            for (let i = 0; i < this.points.length; i++) {
-                M.set(j, i, normFunction(this.points[i], this.points[j]));
+        this.nPoints = this.sourcePoints.length;
+        if (this.nPoints < 3) {
+            throw new Error(`Not enough controle points. A thin plate spline transformation (with affine component) requires a minimum of 3 points, but ${this.nPoints} are given.`);
+        }
+        // 2D Radial Basis Function interpolation
+        // See notebook https://observablehq.com/d/0b57d3b587542794 for code source and explanation
+        // The system of equations is solved for x and y separately (because they are independant)
+        // Hence destinationPointsMatrices and weightsMatrices are an Array of two column vector matrices
+        // Since they both use the same coefficients, there is only one kernelsAndAffineCoefsMatrix
+        const destinationPointsMatrices = [
+            Matrix.columnVector([...destinationPoints, [0, 0], [0, 0], [0, 0]].map((value) => value[0])),
+            Matrix.columnVector([...destinationPoints, [0, 0], [0, 0], [0, 0]].map((value) => value[1]))
+        ];
+        // Pre-compute kernelsMatrix: fill it with the point to point distances between all controle points
+        const kernelsMatrix = Matrix.zeros(this.nPoints, this.nPoints);
+        for (let i = 0; i < this.nPoints; i++) {
+            for (let j = 0; j < this.nPoints; j++) {
+                kernelsMatrix.set(i, j, normFunction(sourcePoints[i], sourcePoints[j]));
             }
         }
-        // if needed, compute espilon as the average distance between points
+        // If it's not provided, and if it's an input to the kernelFunction, compute epsilon as the average distance between the controle points
         if (epsilon === undefined) {
-            epsilon =
-                sumMatrix(M) / (Math.pow(this.points.length, 2) - this.points.length);
+            epsilon = kernelsMatrix.sum() / (Math.pow(this.nPoints, 2) - this.nPoints);
         }
         this.epsilon = epsilon;
-        // update the matrix to reflect the requested distance function
-        for (let j = 0; j < this.points.length; j++) {
-            for (let i = 0; i < this.points.length; i++) {
-                M.set(j, i, distanceFunction(M.get(j, i), this.epsilon));
+        // Finish the computation of kernelsMatrix by applying the requested kernel function
+        for (let i = 0; i < this.nPoints; i++) {
+            for (let j = 0; j < this.nPoints; j++) {
+                kernelsMatrix.set(i, j, kernelFunction(kernelsMatrix.get(i, j), epsilon));
             }
         }
-        // reshape values by dimension/component
-        let valuesByDimension = new Array(2);
-        for (let i = 0; i < 2; i++) {
-            valuesByDimension[i] = this.values.map((value) => value[i]);
+        // Extend kernelsMatrix to include the affine transformation
+        const affineCoefsMatrix = Matrix.zeros(this.nPoints, 3);
+        const kernelsAndAffineCoefsMatrix = Matrix.zeros(this.nPoints + 3, this.nPoints + 3);
+        // Construct Nx3 Matrix affineCoefsMatrix
+        // 1 x0 y0
+        // 1 x1 y1
+        // 1 x2 y2
+        // ...
+        for (let i = 0; i < this.nPoints; i++) {
+            affineCoefsMatrix.set(i, 0, 1);
+            affineCoefsMatrix.set(i, 1, sourcePoints[i][0]);
+            affineCoefsMatrix.set(i, 2, sourcePoints[i][1]);
         }
-        // Compute basis functions weights by solving
-        // the linear system of equations for each target component
-        // let this.weights: [Matrix, Matrix]
-        this.weights = [
-            solve(M, Matrix.columnVector(valuesByDimension[0])),
-            solve(M, Matrix.columnVector(valuesByDimension[1]))
+        // Combine kernelsMatrix and affineCoefsMatrix into new matrix kernelsAndAffineCoefsMatrix
+        // Note: mlMatrix has no knowledge of block matrices, but this approach is good enough
+        // To speed this up, we could maybe use kernelsMatrix.addRow() and kernelsMatrix.addVector()
+        for (let i = 0; i < this.nPoints + 3; i++) {
+            for (let j = 0; j < this.nPoints + 3; j++) {
+                if (i < this.nPoints && j < this.nPoints) {
+                    kernelsAndAffineCoefsMatrix.set(i, j, kernelsMatrix.get(i, j));
+                }
+                else if (i >= this.nPoints && j < this.nPoints) {
+                    kernelsAndAffineCoefsMatrix.set(i, j, affineCoefsMatrix.transpose().get(i - this.nPoints, j));
+                }
+                else if (i < this.nPoints && j >= this.nPoints) {
+                    kernelsAndAffineCoefsMatrix.set(i, j, affineCoefsMatrix.get(i, j - this.nPoints));
+                }
+            }
+        }
+        // Compute basis functions weights and the affine parameters by solving the linear system of equations for each target component
+        // Note: the same kernelsAndAffineCoefsMatrix is used for both solutions
+        const inverseKernelsAndAffineCoefsMatrix = inverse(kernelsAndAffineCoefsMatrix);
+        this.weightsMatrices = [
+            inverseKernelsAndAffineCoefsMatrix.mmul(destinationPointsMatrices[0]),
+            inverseKernelsAndAffineCoefsMatrix.mmul(destinationPointsMatrices[1])
         ];
     }
-    // The returned interpolant will compute the value at any point
-    // by summing the weighted contributions of the input points
-    interpolant(point) {
-        if (!this.weights) {
+    // The interpolant function will compute the value at any point.
+    interpolant(newSourcePoint) {
+        if (!this.weightsMatrices) {
             throw new Error('Weights not computed');
         }
-        let distances = new Array(this.points.length);
-        for (let i = 0; i < this.points.length; i++) {
-            distances[i] = this.distanceFunction(this.normFunction(point, this.points[i]), this.epsilon);
+        // Make a column matrix with the distances of that point to all controle points
+        const newDistancesMatrix = Matrix.zeros(this.nPoints, 1);
+        for (let i = 0; i < this.nPoints; i++) {
+            newDistancesMatrix.set(i, 0, this.kernelFunction(this.normFunction(newSourcePoint, this.sourcePoints[i]), this.epsilon));
         }
-        const distancesM = new Matrix(distances.map((d) => [d]));
-        let sums = [0, 0];
+        // Compute the interpolated value by summing the weighted contributions of the input point
+        const newDestinationPoint = [0, 0];
         for (let i = 0; i < 2; i++) {
-            sums[i] = sumMatrix(mulitplyMatricesValues(distancesM, this.weights[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] = multiplyMatricesElementwise(newDistancesMatrix, this.weightsMatrices[i].selection([...Array(this.nPoints).keys()], [0])).sum();
+            // Add the affine part
+            const a0 = this.weightsMatrices[i].get(this.nPoints, 0);
+            const ax = this.weightsMatrices[i].get(this.nPoints + 1, 0);
+            const ay = this.weightsMatrices[i].get(this.nPoints + 2, 0);
+            newDestinationPoint[i] +=
+                a0 + ax * newSourcePoint[0] + ay * newSourcePoint[1];
         }
-        return sums;
+        return newDestinationPoint;
     }
 }

dist/shared/types.d.ts

@@ -0,1 +1,2 @@
+export type TransformationType = 'helmert' | 'polynomial' | 'polynomial1' | 'polynomial2' | 'polynomial3' | 'projective' | 'thinPlateSpline';
 export type Position = [number, number];
@@ -16,11 +17,12 @@ export type BBox = [number, number, number, number];
 };
-export type ImageWorldPosition = {
-    image: Position;
-    world: Position;
+export type GCP = {
+    resource: Position;
+    geo: Position;
 };
 export type Segment = {
-    from: ImageWorldPosition;
-    to: ImageWorldPosition;
+    from: GCP;
+    to: GCP;
 };
 export type TransformOptions = {
+    close: boolean;
     maxOffsetRatio: number;
@@ -30,1 +32,8 @@ maxDepth: number;
 export type OptionalTransformOptions = Partial<TransformOptions>;
+export type KernelFunction = (r: number, epsilon?: number) => number;
+export type NormFunction = (point1: Position, point2: Position) => number;
+export type GCPTransformerInterface = {
+    gcps: GCP[];
+    toGeo(point: Position): Position;
+    toResource(point: Position): Position;
+};

dist/transformer.d.ts

@@ -1,5 +0,20 @@
-import type { Position, ImageWorldPosition } from './shared/types.js';
-import type { GCPTransformInfo } from './gdaltransform.js';
-export declare function toWorld(transformArgs: GCPTransformInfo, point: Position): Position;
-export declare function toImage(transformArgs: GCPTransformInfo, point: Position): Position;
-export declare function createTransformer(gcps: ImageWorldPosition[]): GCPTransformInfo;
+import type { GCPTransformerInterface, TransformationType, OptionalTransformOptions, Position, GCP, GeoJSONGeometry, GeoJSONPoint, GeoJSONLineString, GeoJSONPolygon } from './shared/types.js';
+export default class GCPTransformer implements GCPTransformerInterface {
+    gcps: GCP[];
+    type: TransformationType;
+    transformer: GCPTransformerInterface;
+    constructor(gcps: GCP[], type?: TransformationType);
+    toGeo(point: Position): Position;
+    toResource(point: Position): Position;
+    toGeoJSON(point: Position): GeoJSONGeometry;
+    toGeoJSON(points: Position[], options?: OptionalTransformOptions): GeoJSONGeometry;
+    toGeoJSONPoint(point: Position): GeoJSONPoint;
+    toGeoJSONLineString(points: Position[], options?: OptionalTransformOptions): GeoJSONLineString;
+    toGeoJSONPolygon(points: Position[], options?: OptionalTransformOptions): GeoJSONPolygon;
+    fromGeoJSON(geometry: GeoJSONGeometry, options?: OptionalTransformOptions): Position | Position[];
+    fromGeoJSONPoint(geometry: GeoJSONPoint): Position;
+    fromGeoJSONLineString(geometry: GeoJSONLineString, options?: OptionalTransformOptions): Position[];
+    fromGeoJSONPolygon(geometry: GeoJSONPolygon, options?: OptionalTransformOptions): Position[];
+    toGeoPolygon(polygon: Position[], options?: OptionalTransformOptions): Position[];
+    toResourcePolygon(polygon: Position[], options?: OptionalTransformOptions): Position[];
+}

dist/transformer.js

@@ -1,18 +0,99 @@
-import { GCP, GDALCreateGCPTransformer, GDALGCPTransform } from './gdaltransform.js';
-export function toWorld(transformArgs, point) {
-    const bInverse = false;
-    const input = [{ x: point[0], y: point[1] }];
-    const output = GDALGCPTransform(transformArgs, bInverse, input);
-    return [output[0].y, output[0].x];
+import HelmertGCPTransformer from './transformers/helmert-transformer.js';
+import PolynomialGCPTransformer from './transformers/polynomial-transformer.js';
+import ProjectiveGCPTransformer from './transformers/projective-transformer.js';
+import RadialBasisFunctionGCPTransformer from './transformers/radial-basis-function-transformer.js';
+import { thinPlateKernel } from './shared/kernel-functions.js';
+import { toGeoJSONPoint, toGeoJSONLineString, toGeoJSONPolygon, fromGeoJSONPoint, fromGeoJSONLineString, fromGeoJSONPolygon, toGeoPolygon, toResourcePolygon } from './shared/geojson.js';
+export default class GCPTransformer {
+    constructor(gcps, type = 'polynomial'
+    // options: TransformationOptions
+    ) {
+        this.gcps = gcps;
+        this.type = type;
+        if (type === 'helmert') {
+            this.transformer = new HelmertGCPTransformer(gcps);
+        }
+        else if (type === 'polynomial1' || type === 'polynomial') {
+            this.transformer = new PolynomialGCPTransformer(gcps);
+        }
+        else if (type === 'polynomial2') {
+            this.transformer = new PolynomialGCPTransformer(gcps, 2);
+        }
+        else if (type === 'polynomial3') {
+            this.transformer = new PolynomialGCPTransformer(gcps, 3);
+        }
+        else if (type === 'projective') {
+            this.transformer = new ProjectiveGCPTransformer(gcps);
+        }
+        else if (type === 'thinPlateSpline') {
+            this.transformer = new RadialBasisFunctionGCPTransformer(gcps, thinPlateKernel);
+        }
+        else {
+            throw new Error(`Unsupported transformation type: ${type}`);
+        }
+    }
+    toGeo(point) {
+        return this.transformer.toGeo(point);
+    }
+    toResource(point) {
+        return this.transformer.toResource(point);
+    }
+    toGeoJSON(pointOrPoints, options) {
+        // TODO: also support empty point and points arrays by throwing error
+        const isPoints = Array.isArray(pointOrPoints[0]);
+        let point;
+        let points;
+        if (isPoints) {
+            points = pointOrPoints;
+            const close = options && options.close;
+            if (close) {
+                return toGeoJSONPolygon(this, points, options);
+            }
+            else {
+                return toGeoJSONLineString(this, points, options);
+            }
+        }
+        else {
+            point = pointOrPoints;
+            return toGeoJSONPoint(this, point);
+        }
+    }
+    toGeoJSONPoint(point) {
+        return toGeoJSONPoint(this, point);
+    }
+    toGeoJSONLineString(points, options) {
+        return toGeoJSONLineString(this, points, options);
+    }
+    toGeoJSONPolygon(points, options) {
+        return toGeoJSONPolygon(this, points, options);
+    }
+    fromGeoJSON(geometry, options) {
+        if (geometry.type === 'Point') {
+            return fromGeoJSONPoint(this, geometry);
+        }
+        else if (geometry.type === 'LineString') {
+            return fromGeoJSONLineString(this, geometry, options);
+        }
+        else if (geometry.type === 'Polygon') {
+            return fromGeoJSONPolygon(this, geometry, options);
+        }
+        else {
+            throw new Error(`Unsupported geometry`);
+        }
+    }
+    fromGeoJSONPoint(geometry) {
+        return fromGeoJSONPoint(this, geometry);
+    }
+    fromGeoJSONLineString(geometry, options) {
+        return fromGeoJSONLineString(this, geometry, options);
+    }
+    fromGeoJSONPolygon(geometry, options) {
+        return fromGeoJSONPolygon(this, geometry, options);
+    }
+    toGeoPolygon(polygon, options) {
+        return toGeoPolygon(this, polygon, options);
+    }
+    toResourcePolygon(polygon, options) {
+        return toResourcePolygon(this, polygon, options);
+    }
 }
-export function toImage(transformArgs, point) {
-    const bInverse = true;
-    const input = [{ x: point[1], y: point[0] }];
-    const output = GDALGCPTransform(transformArgs, bInverse, input);
-    return [output[0].x, output[0].y];
-}
-export function createTransformer(gcps) {
-    const pasGCPs = gcps.map((gcp) => new GCP(gcp.image[0], gcp.image[1], gcp.world[1], gcp.world[0]));
-    const nOrder = 0;
-    return GDALCreateGCPTransformer(pasGCPs, nOrder, false);
-}

dist/transformers/polynomial-transformer.d.ts

@@ -1,9 +0,15 @@
-import type { GCPTransformInfo } from '../shared/gdaltransform.js';
-import type { GCPTransformerInterface, Position, ImageWorldPosition } from '../shared/types.js';
+import Polynomial from '../shared/polynomial.js';
+import type { GCPTransformerInterface, Position, GCP } from '../shared/types.js';
 export default class PolynomialGCPTransformer implements GCPTransformerInterface {
-    gcps: ImageWorldPosition[];
-    transformArgs: GCPTransformInfo;
-    constructor(gcps: ImageWorldPosition[]);
-    toWorld(point: Position): Position;
+    gcps: GCP[];
+    gcpGeoCoords: Position[];
+    gcpResourceCoords: Position[];
+    toGeoPolynomial?: Polynomial;
+    toResourcePolynomial?: Polynomial;
+    order?: number;
+    constructor(gcps: GCP[], order?: number);
+    createToGeoPolynomial(): Polynomial;
+    createToResourcePolynomial(): Polynomial;
+    toGeo(point: Position): Position;
     toResource(point: Position): Position;
 }

dist/transformers/polynomial-transformer.js

@@ -1,21 +0,29 @@
-import { GCP, GDALCreateGCPTransformer, GDALGCPTransform } from '../shared/gdaltransform.js';
+import Polynomial from '../shared/polynomial.js';
 export default class PolynomialGCPTransformer {
-    constructor(gcps) {
+    constructor(gcps, order) {
         this.gcps = gcps;
-        const pasGCPs = gcps.map((gcp) => new GCP(gcp.image[0], gcp.image[1], gcp.world[1], gcp.world[0]));
-        const nOrder = 0;
-        this.transformArgs = GDALCreateGCPTransformer(pasGCPs, nOrder, false);
+        this.gcpGeoCoords = gcps.map((gcp) => gcp.geo);
+        this.gcpResourceCoords = gcps.map((gcp) => gcp.resource);
+        if (order) {
+            this.order = order;
+        }
     }
-    toWorld(point) {
-        const bInverse = false;
-        const input = [{ x: point[0], y: point[1] }];
-        const output = GDALGCPTransform(this.transformArgs, bInverse, input);
-        return [output[0].y, output[0].x];
+    createToGeoPolynomial() {
+        return new Polynomial(this.gcpResourceCoords, this.gcpGeoCoords, this.order);
     }
+    createToResourcePolynomial() {
+        return new Polynomial(this.gcpGeoCoords, this.gcpResourceCoords, this.order);
+    }
+    toGeo(point) {
+        if (!this.toGeoPolynomial) {
+            this.toGeoPolynomial = this.createToGeoPolynomial();
+        }
+        return this.toGeoPolynomial.interpolant(point);
+    }
     toResource(point) {
-        const bInverse = true;
-        const input = [{ x: point[1], y: point[0] }];
-        const output = GDALGCPTransform(this.transformArgs, bInverse, input);
-        return [output[0].x, output[0].y];
+        if (!this.toResourcePolynomial) {
+            this.toResourcePolynomial = this.createToResourcePolynomial();
+        }
+        return this.toResourcePolynomial.interpolant(point);
     }
 }

dist/transformers/radial-basis-function-transformer.d.ts

 import RBF from '../shared/radial-basis-function.js';
-import type { GCPTransformerInterface, DistanceFunction, NormFunction, Position, ImageWorldPosition } from '../shared/types.js';
+import type { GCPTransformerInterface, KernelFunction, NormFunction, Position, GCP } from '../shared/types.js';
 export default class RadialBasisFunctionGCPTransformer implements GCPTransformerInterface {
-    gcps: ImageWorldPosition[];
-    worldGcps: Position[];
-    resourceGcps: Position[];
-    distanceFunction: DistanceFunction;
+    gcps: GCP[];
+    gcpGeoCoords: Position[];
+    gcpResourceCoords: Position[];
+    kernelFunction: KernelFunction;
     normFunction: NormFunction;
-    toWorldRbf?: RBF;
+    toGeoRbf?: RBF;
     toResourceRbf?: RBF;
-    constructor(gcps: ImageWorldPosition[], distanceFunction: DistanceFunction);
-    createToWorldRbf(): RBF;
+    constructor(gcps: GCP[], kernelFunction: KernelFunction);
+    createToGeoRbf(): RBF;
     createToResourceRbf(): RBF;
-    toWorld(point: Position): Position;
+    toGeo(point: Position): Position;
     toResource(point: Position): Position;
 }

dist/transformers/radial-basis-function-transformer.js

 import RBF from '../shared/radial-basis-function.js';
 import { euclideanNorm } from '../shared/norm-functions.js';
 export default class RadialBasisFunctionGCPTransformer {
-    constructor(gcps, distanceFunction) {
+    constructor(gcps, kernelFunction) {
         this.gcps = gcps;
-        this.worldGcps = gcps.map((gcp) => gcp.world);
-        this.resourceGcps = gcps.map((gcp) => gcp.image);
-        this.distanceFunction = distanceFunction;
+        this.gcpGeoCoords = gcps.map((gcp) => gcp.geo);
+        this.gcpResourceCoords = gcps.map((gcp) => gcp.resource);
+        this.kernelFunction = kernelFunction;
         this.normFunction = euclideanNorm;
     }
-    createToWorldRbf() {
-        return new RBF(this.resourceGcps, this.worldGcps, this.distanceFunction, this.normFunction);
+    createToGeoRbf() {
+        return new RBF(this.gcpResourceCoords, this.gcpGeoCoords, this.kernelFunction, this.normFunction);
     }
     createToResourceRbf() {
-        return new RBF(this.worldGcps, this.resourceGcps, this.distanceFunction, this.normFunction);
+        return new RBF(this.gcpGeoCoords, this.gcpResourceCoords, this.kernelFunction, this.normFunction);
     }
-    toWorld(point) {
-        if (!this.toWorldRbf) {
-            this.toWorldRbf = this.createToWorldRbf();
+    toGeo(point) {
+        if (!this.toGeoRbf) {
+            this.toGeoRbf = this.createToGeoRbf();
         }
-        return this.toWorldRbf.interpolant(point);
+        return this.toGeoRbf.interpolant(point);
     }
@@ -23,0 +23,0 @@ toResource(point) {

package.json

 {
   "name": "@allmaps/transform",
-  "version": "1.0.0-beta.14",
-  "author": {
-    "name": "Bert Spaan",
-    "email": "[e-mail]",
-    "url": "https://bertspaan.nl"
-  },
-  "description": "Coordinate transformation functions, based on gdaltransform",
+  "version": "1.0.0-beta.15",
+  "contributors": [
+    {
+      "name": "Bert Spaan",
+      "email": "[e-mail]",
+      "url": "https://bertspaan.nl"
+    },
+    {
+      "name": "Manuel Claeys Bouuaert",
+      "email": "[e-mail]",
+      "url": "https://manuelclaeysbouuaert.be"
+    }
+  ],
+  "description": "Coordinate transformation functions",
   "type": "module",
@@ -16,3 +23,4 @@ "types": "./dist/index.d.ts",
     ".": {
-      "import": "./dist/index.js"
+      "types": "./dist/index.d.ts",
+      "default": "./dist/index.js"
     }
@@ -39,6 +47,8 @@ },
     "@turf/distance": "^6.3.0",
-    "@turf/midpoint": "^6.3.0"
+    "@turf/midpoint": "^6.3.0",
+    "ml-matrix": "^6.10.4"
   },
   "scripts": {
     "watch": "tsc --watch",
+    "bench": "node bench/index.js",
     "build": "tsc",
@@ -64,3 +74,3 @@ "test": "NODE_ENV=test mocha",
   },
-  "gitHead": "e1277b2bc286652d62751976b39b9f81e7281225"
+  "gitHead": "45c0da2e8fa40e3d97a6a5d2f3531d6c03bd34aa"
 }

README.md

 # @allmaps/transform
 
-JavaScript port of [gdaltransform](https://gdal.org/programs/gdaltransform.html).
+This module serves to **transform points**, polygons and other spatial features from one cartesian `(x, y)` plane to another **using a set of control points** identified in both planes.
 
-See https://observablehq.com/@bertspaan/gdaltransform?collection=@bertspaan/iiif-maps.
+It is used in [@allmaps/render](../../packages/render/) and [@allmaps/tileserver](../../apps/tileserver/), two packages where we produce a georeferenced image by triangulating a IIIF image and drawing these triangles on a map in a specific new location, with the triangle's new vertex location computed by the transformer of this package. The transformer is constructed from controle points in the annotation and transforms coordinates **from the pixel coordinates space** of a IIIF Resource **to the map coordinate space** of an interactive map.
 
-## See also
+## How it works
 
-- https://khanhha.github.io/posts/Thin-Plate-Splines-Warping/
-- http://eeweb.poly.edu/~yao/EL5123/lecture12_ImageWarping.pdf
-- https://www.comp.nus.edu.sg/~cs4340/lecture/imorph.pdf
-- http://www.corrmap.com/features/rubber-sheeting_transformation.php
-- https://stackoverflow.com/questions/32266408/warping-an-image-using-control-points#32498207
-- https://stackoverflow.com/questions/34884779/whats-a-simple-way-of-warping-an-image-with-a-given-set-of-points
-- https://www.fil.ion.ucl.ac.uk/spm/doc/books/hbf2/pdfs/Ch4.pdf
-- https://en.wikipedia.org/wiki/Multivariate_interpolation#Irregular_grid_.28scattered_data.29
-- https://observablehq.com/@esperanc/deformation-with-rbfs
-- https://en.wikipedia.org/wiki/Radial_basis_function
-- https://stackoverflow.com/questions/34884779/whats-a-simple-way-of-warping-an-image-with-a-given-set-of-points
+This library started out as a JavaScript port of [gdaltransform](https://gdal.org/programs/gdaltransform.html) (as described in [this notebook](https://observablehq.com/@bertspaan/gdaltransform?collection=@bertspaan/iiif-maps)) and initially only implemented polynomial transformations of order 1. Later Thin Plate Spline transformations were added (see [this notebook](https://observablehq.com/d/0b57d3b587542794)) amongst other transformations, which lead to a refactoring using the [`ml-matrix`](https://github.com/mljs/matrix) library applied for creating and solving the linear systems of equations which are the essence of each of these transformations.
+
+The algorithms correspond to those of **GDAL** and the results are (nearly) identical as can be checked in the [tests](./test/test-transform.js).
+
+These are the **supported transformations**:
+
+| Type              | Options    | Description                                                              | Properties                                                                       | Minimum number of GCPs |
+| ----------------- | ---------- | ------------------------------------------------------------------------ | -------------------------------------------------------------------------------- | ---------------------- |
+| `helmert`         |            | Helmert transformation or 'similarity transformation'                    | Preserves shape and angle                                                        | 2                      |
+| `polynomial`      | `order: 1` | First order polynomial transformation                                    | Preserves lines and parallelism                                                  | 3                      |
+| `polynomial`      | `order: 2` | Second order polynomial transformation                                   | Some bending flexibility                                                         | 6                      |
+| `polynomial`      | `order: 3` | Third order polynomial transformation                                    | More bending flexibility                                                         | 10                     |
+| `thinPlateSpline` |            | Thin Plate Spline transformation or 'rubber sheeting' (with affine part) | Exact, smooth (see [this notebook](https://observablehq.com/d/0b57d3b587542794)) | 3                      |
+| `projective`      |            | Projective or 'perspective' transformation, used for aerial images       | Preserves lines and cross-ratios                                                 | 4                      |
+
+## Usage
+
+```js
+import { GCPTransformer } from '@allmaps/transform'
+
+// ground control points
+const gcps = [
+  {
+    image: [518, 991],
+    world: [4.9516614, 52.4633102]
+  },
+  {
+    image: [4345, 2357],
+    world: [5.0480391, 52.5123762]
+  },
+  {
+    image: [2647, 475],
+    world: [4.9702906, 52.5035815]
+  }
+]
+
+// create a transformation of a specific type using the ground control points
+const transformer = new GCPTransformer(gcps, 'helmert')
+
+// forward transform
+const result = transformer.toGeo([100, 100])
+// result = [4.9385700843392435, 52.46580484503631]
+
+// backward transform
+const result = transformer.toResource([4.9385700843392435, 52.46580484503631])
+// result = [100, 100]
+```
+
+The transformer Class also exports the function `transformer.toGeoJSON()` for forward transforms (Arrays of) \[x, y] position(s) to the geometry object of the appropriate GeoJSON object (Point, LineString, Polygon), and the function `transformer.fromGeoJSON()` to take in such a GeoJSON object and backwards transform it to an (Array of) \[x, y] position(s).
+
+```js
+const geoJSONPointGeometry = transformer.toGeoJSON([100, 100])
+// geojsonPointGeometry = {
+//     type: "Point",
+//     coordinates: [4.9385700843392435, 52.46580484503631]
+//   }
+const geoJSONPoint = {
+  type: 'Feature',
+  geometry: geojsonPointGeometry
+}
+// geojsonPoint = {
+//   "type": "Feature",
+//   "geometry": {
+//     type: "Point",
+//     coordinates: [4.9385700843392435, 52.46580484503631]
+//   }
+// }
+
+const point = transformer.fromGeoJSON(geoJSONPointGeometry)
+// point = [100, 100]
+```
+
+Both of these functions can take a second argument with options for this transformation
+
+| Option           | Description                                                     | Default |
+| :--------------- | :-------------------------------------------------------------- | :------ |
+| `maxOffsetRatio` | Maximum offset ratio between original and transformed midpoints | 0       |
+| `maxDepth`       | Maximum recursion depth                                         | 6       |
+
+```js
+const options = {
+  maxOffsetRatio: 0.01,
+  maxDepth: 6
+}
+const geoJSONLineGeometry = transformer.toGeoJSON(
+  [
+    [100, 100],
+    [110, 115],
+    [120, 170]
+  ],
+  options
+)
+/// geoJSONLineGeometry = ...
+```
+
+### CLI
+
+The [@allmaps/cli](../../apps/cli/) package exports an interface to transform **SVG** objects from the pixel coordinates space of a IIIF Resource to **GeoJSON** objects in the map coordinate space of an interactive map or vise versa **given (the ground control points and transformation type from) a Georeference Annotation**, and to export the SVG resource mask included in a Georeference Annotation as a GeoJSON object.
+
+In the future, we hope to add an interface to manually enter coordinates of points to transform forward or backwards and easily specify GCPs and transformation type.
+
+### Notes
+
+- Only **linearly independent control points** should be considered when checking if the criterion for the minimum number of control points is met. For example, three control points that are collinear (one the same line) only count as two linearly independent points. The current implementation doesn't check such linear (in)dependance, but building a transformer with insufficient linearly independent control points will result in a badly conditioned matrix (no error but diverging results) or non-invertible matrix (**error when inverting matrix**).
+- The examples here use `(longitude, latitude)` coordinates in the destination space, but it is important to stress that the functions `toGeo()` and `toResource()` are currently map-projection agnostic forward and backward transformations built for the general case of a transformation for one cartesian `(x, y)` plane to another (even though their name suggests otherwise). **When working in a geographic context** one can use control points with `(longitude, latitude)` coordinates in the destination space, which will build a transformation to the cartesian plane of an equidistant projection. One could also transform such coordinates to a specific projection (for example Mercator) first and use control points with such projection coordinates in the destination space. This will build a transformation to the cartesian space of a Mercator projection. The output of the `toResource()` function is then also in these coordinates. Finally, since `toGeoJSON()` and `fromGeoJSON()` were specifically built with `(longitude, latitude)` coordinates in mind (as they are used in GeoJSON objects), these functions should only be used for such coordinates.
+
+### Benchmark
+
+Here are some benchmarks on building and using a transformer.
+
+Create transformer with 10 points (and transform 1 point)
+
+| Type              | Options    | Ops/s |
+| ----------------- | ---------- | ----- |
+| `helmert`         |            | 20813 |
+| `polynomial`      | `order: 1` | 32631 |
+| `polynomial`      | `order: 2` | 20927 |
+| `polynomial`      | `order: 3` | 6682  |
+| `thinPlateSpline` |            | 4532  |
+| `projective`      |            | 7044  |
+
+Use transformer made with with 10 points
+
+| Type              | Options    | Ops/s   |
+| ----------------- | ---------- | ------- |
+| `helmert`         |            | 7561796 |
+| `polynomial`      | `order: 1` | 5546792 |
+| `polynomial`      | `order: 2` | 4811662 |
+| `polynomial`      | `order: 3` | 887390  |
+| `thinPlateSpline` |            | 154310  |
+| `projective`      |            | 5265895 |
+
+See `./bench/index.js`

New alerts

Non-permissive License

License

(Experimental) A license not known to be considered permissive was found.

Found 1 instance in 1 package

Fixed alerts

Mixed license

License

(Experimental) Package contains multiple licenses.

Found 1 instance in 1 package

Non-permissive License

License

(Experimental) A license not known to be considered permissive was found.

Found 1 instance in 1 package

Improved metrics

Total package byte prevSize

110490

increased by27.1%

Number of package files

46

increased by39.39%

Number of low license alerts

2

decreased by-33.33%

Lines of code

2479

increased by37.34%

Number of lines in readme file

142

increased by647.37%

Worsened metrics

Dependency count

4

increased by33.33%

Dependency changes

• + Addedml-matrix@^6.10.4

• + Addedis-any-array@2.0.1(transitive)

• + Addedml-array-max@1.2.4(transitive)

• + Addedml-array-min@1.2.3(transitive)

• + Addedml-array-rescale@1.3.7(transitive)

• + Addedml-matrix@6.12.0(transitive)