portable / efficient JS version of OpenCV getPerspectiveTransform()

getPerspectiveTransform.mts

/**
 * solve a 2x2 linear system, given its augmented matrix, in place.
 * first 2 columns left unchanged.
 */
function solve2x2([ x, y ]: [number[], number[]]) {
	let [ a, b, c, d ] = [ x[0], x[1], y[0], y[1] ]
	const det = a * d - b * c
	if (!det) throw Error('inversion failed')
	const k = 1 / det
	; (a *= k, b *= k, c *= k, d *= k)
	for (let i = 2; i < x.length; i++)
		[ x[i], y[i] ] = [ d * x[i] - b * y[i], a * y[i] - c * x[i] ]
}

type Vec2 = [number, number]

/**
 * given 4 input / output pairs, compute its 3x3 perspective matrix
 * (last slot is always 1)
 */
export default function getPerspectiveTransform(points: [Vec2, Vec2][]): number[][] {
	// prepare coefficients from points
	if (points.length !== 4)
		throw Error('expected 4 pairs of points')
	const M = points.map(([ [xi, yi], [xo, yo] ]) =>
		[ xi, yi, xo, yo, -xi*xo, -xi*yo, -yi*xo, -yi*yo ])
	for (let r = 0; r < 3; r++)
		for (let i = 0; i < 8; i++)
			M[r][i] -= M[3][i]
	; [ M[2], M[3] ] = [ M[3], M[2] ]

	// solve first system (M)
	solve2x2(M)
	for (let r = 2; r < 4; r++)
		for (let i = 2; i < 8; i++)
			for (let v = 0; v < 2; v++)
				M[r][i] -= M[r][v] * M[v][i]

	// solve second system (transpose of M[3].slice(2))
	const L = [0, 1].map(i => [4, 6, 2].map(l => M[3][l + i]))
	solve2x2(L)
	for (let r = 0; r < 3; r++)
		for (let i = 2; i < 4; i++)
			for (let v = 0; v < 2; v++)
				M[r][i] -= M[r][2 + 2 * v + i] * L[v][2]

	return [ [ M[0][2], M[1][2], M[2][2] ] ,
		 [ M[0][3], M[1][3], M[2][3] ] ,
		 [ L[0][2], L[1][2],       1 ] ]
}