homography.lua

-- Requires https://github.com/davidm/lua-matrix/blob/master/lua/matrix.lua
local matrix = require("matrix")
local h = {}

-- Adapted from https://github.com/nenoNaninu/HomographySharp/blob/master/HomographySharp/Single/SingleHomography.cs

function SetXCoefficients(mx, srcX, srcY, dstX, row)
    mx[row][1] = srcX
    mx[row][2] = srcY
    mx[row][3] = 1

    mx[row][4] = 0
    mx[row][5] = 0
    mx[row][6] = 0
    
    mx[row][7] = -dstX * srcX
    mx[row][8] = -dstX * srcY
end
function SetYCoefficients(mx, srcX, srcY, dstY, row)
    mx[row][1] = 0
    mx[row][2] = 0
    mx[row][3] = 0

    mx[row][4] = srcX
    mx[row][5] = srcY
    mx[row][6] = 1
    
    mx[row][7] = -dstY * srcX
    mx[row][8] = -dstY * srcY
end

--[[
    w is width, h is height, v is default value
]]
function h.CreateMatrix(w, h, v)
    local m = {}
    for x = 1, w do
        m[x] = {}
        for y = 1, h do
            m[x][y] = v
        end
    end
    return m
end
--[[
Example:
f = h.FindHomography(
    {{176.69, -4.07}, {0, 0}}, -- top left
    {{176.65, 19.30}, {0, h}}, -- bottom left
    {{200.03, 18.25}, {w, h}}, -- bottom right
    {{200.05, -3.92}, {w, 0}} -- top right
)
f(176.69, -4.07) approximately equals {0,0}
]]--

function h.FindHomography(p)
    local coeff = h.CreateMatrix(#p*2,8,0)
    local destination = {}

    for i, point in ipairs(p) do
        SetXCoefficients(coeff, point[1][1], point[1][2], point[2][1], (2 * (i-1))+1 )
        SetYCoefficients(coeff, point[1][1], point[1][2], point[2][2], (2 * (i-1))+2 )

        destination[(2 * (i-1))+1] = {point[2][1]}
        destination[(2 * (i-1))+2] = {point[2][2]}
    end

    local invM = matrix.invert(coeff)
    local param = matrix.mul(invM, destination)
    
    local v = {
        param[1][1], param[2][1], param[3][1], param[4][1],
        param[5][1], param[6][1], param[7][1], param[8][1],
        1
    }

    return function(x, y)
        local dst1 = v[1] * x + v[2] * y + v[3]
        local dst2 = v[4] * x + v[5] * y + v[6]
        local dst3 = v[7] * x + v[8] * y + v[9]

        return { dst1 / dst3, dst2 / dst3}
    end
end

return h