How to calculate a transform matrix from set of tie points for perspective transform

Question

I am trying to implement an equivalent of imagemagick perspective distort using libvips http://www.imagemagick.org/Usage/distorts/#perspective

As per this answer How to perform perspective distort transformation in VIPS? it can be achieved using mapim

However I don't understand how can I convert a set of tie points as used in imagemagick into a trasform matrix of this kind

T = [1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0003, 0.0001]

Answer 1

The transform vector can be found by solving this system of linear equations:

/ x0 y0  1  0  0  0 -x0*u0 -y0*u0 \ /c00\ /u0\
| x1 y1  1  0  0  0 -x1*u1 -y1*u1 | |c01| |u1|
| x2 y2  1  0  0  0 -x2*u2 -y2*u2 | |c02| |u2|
| x3 y3  1  0  0  0 -x3*u3 -y3*u3 |.|c10|=|u3|
|  0  0  0 x0 y0  1 -x0*v0 -y0*v0 | |c11| |v0|
|  0  0  0 x1 y1  1 -x1*v1 -y1*v1 | |c12| |v1|
|  0  0  0 x2 y2  1 -x2*v2 -y2*v2 | |c20| |v2|
\  0  0  0 x3 y3  1 -x3*v3 -y3*v3 / \c21/ \v3/

where x0.., y0.. are the target coordinates, u0.. v0.. are source coordinates.

Here is how to solve it using gonum package in go:

func calculateTransformation(coordinates []float64) mat.VecDense {

    u0, v0, x0, y0 := coordinates[0], coordinates[1], coordinates[2], coordinates[3]
    u1, v1, x1, y1 := coordinates[4], coordinates[5], coordinates[6], coordinates[7]
    u2, v2, x2, y2 := coordinates[8], coordinates[9], coordinates[10], coordinates[11]
    u3, v3, x3, y3 := coordinates[12], coordinates[13], coordinates[14], coordinates[15]

    // The data must be arranged in row-major order, i.e. the (i*c + j)-th
    // element in the data slice is the {i, j}-th element in the matrix.
    Adata := []float64{
        x0, y0, 1, 0, 0, 0, -x0 * u0, -y0 * u0,
        x1, y1, 1, 0, 0, 0, -x1 * u1, -y1 * u1,
        x2, y2, 1, 0, 0, 0, -x2 * u2, -y2 * u2,
        x3, y3, 1, 0, 0, 0, -x3 * u3, -y3 * u3,
        0, 0, 0, x0, y0, 1, -x0 * v0, -y0 * v0,
        0, 0, 0, x1, y1, 1, -x1 * v1, -y1 * v1,
        0, 0, 0, x2, y2, 1, -x2 * v2, -y2 * v2,
        0, 0, 0, x3, y3, 1, -x3 * v3, -y3 * v3,
    }
    A := mat.NewDense(8, 8, Adata)
    b := mat.NewVecDense(8, []float64{u0, u1, u2, u3, v0, v1, v2, v3})
    result := mat.VecDense{}
    var err error
    err = result.SolveVec(A, b)
    if err != nil {
        fmt.Print(err)
    }
    return result
}