Minimal code to calculate an affine transformation to match three source points to three destination points. Has no dependencies.

affineTransform.ts

type Point = {
    x: number;
    y: number;
}

type Vector = number[];

type Matrix = Vector[];

type Vector2D = [number, number];

type Matrix2D = [Vector2D, Vector2D];

function getAffineTransformation(inGCPs: [Vector2D, Vector2D, Vector2D], outGCPs: [Vector2D, Vector2D, Vector2D]) {
    // magic from https://stackoverflow.com/a/57696101
    const det = (m: Matrix): number =>
        m.length == 1 ?
            m[0][0] :
            m.length == 2 ?
                m[0][0] * m[1][1] - m[0][1] * m[1][0] :
                m[0].reduce((r, e, i) =>
                    r + (-1) ** (i + 2) * e * det(m.slice(1).map(c =>
                        c.filter((_, j) => i != j))), 0);

    // matrix functions
    const copy = <T>(m: T): T => JSON.parse(JSON.stringify(m));
    const transpose = (m: Matrix) => m[0].map((_, i) => m.map(x => x[i]));
    const ones = (n: number) => new Array(n).fill(1) as 1[];
    const del = (m: Matrix, index: number, axis: number) => {
        if (axis !== 0) throw new Error("Unsupported axis.");
        m = copy(m); // deep copy
        m.splice(index, 1);
        return m;
    };
    const inv = (m: Matrix2D): Matrix2D => {
        if (m.length !== 2 || m[0].length !== 2) throw new Error("Unsupported dimensions.");
        const x = 1 / det(m);
        return [
            [m[1][1] * x, -m[0][1] * x],
            [-m[1][0] * x, m[0][0] * x],
        ] as Matrix2D;
    };

    // magic from https://stackoverflow.com/a/56259087
    const n = inGCPs.length;
    const B = [...transpose(inGCPs), ones(n)];
    const D = 1.0 / det(B);
    const entry = (r: Vector, d: number) => det(del([r, ...B], d + 1, 0));
    const range = [...Array(n).keys()];
    const M = transpose(outGCPs).map(R => range.map(i => (-1) ** i * D * entry(R, i)));
    const A = M.map(x => x.slice(0, n - 1));
    const t = transpose(M.map(x => x.slice(n - 1)))[0];

    // process results
    const matrix = A as Matrix2D;
    const vector = t as Vector2D;
    const matrixInv = inv(matrix);

    const apply = (gcp: Vector2D): Vector2D => [
        matrix[0][0] * gcp[0] + matrix[0][1] * gcp[1] + vector[0],
        matrix[1][0] * gcp[0] + matrix[1][1] * gcp[1] + vector[1],
    ];
    const applyInv = (gcp: Vector2D): Vector2D => {
        gcp = [gcp[0] - vector[0], gcp[1] - vector[1]];
        return [
            matrixInv[0][0] * gcp[0] + matrixInv[0][1] * gcp[1],
            matrixInv[1][0] * gcp[0] + matrixInv[1][1] * gcp[1],
        ]
    };

    return {
        matrix,
        vector,
        apply,
        matrixInv,
        applyInv,
    }
}

function example() {
    const inGCPs_: [Point, Point, Point] = [{ x: 1, y: 1 }, { x: 2, y: 3 }, { x: 3, y: 2 }]  // <- points
    const outGCPs_: [Point, Point, Point] = [{ x: 0, y: 2 }, { x: 1, y: 2 }, { x: -2, y: -1 }]  // <- mapped to
    const inGCPs = inGCPs_.map(gcp => [gcp.x, gcp.y]) as [Vector2D, Vector2D, Vector2D];
    const outGCPs = outGCPs_.map(gcp => [gcp.x, gcp.y]) as [Vector2D, Vector2D, Vector2D];

    const result = getAffineTransformation(inGCPs, outGCPs);
    console.log("Result:\n", result);

    // round test results to 14dp
    const round = (n: number) => Math.round((n + Number.EPSILON) * 1e14) / 1e14;
    const roundVector = (v: Vector2D) => [round(v[0]), round(v[1])];

    inGCPs.forEach(gcp => {
        console.log("Test: Map input point", gcp, "to output point", roundVector(result.apply(gcp)));
    });
    outGCPs.forEach(gcp => {
        console.log("Test: Map output point", gcp, "to input point", roundVector(result.applyInv(gcp)));
    });
}

example();