Julia version of cv.getPerspectiveTransform?

[vshesh]

I want to create and use a projective transform matrix between four points in one coordinate system to four points in another coordinate system. Is there a way to do this in Julia? In OpenCV I can call cv.getPerspectiveTransform and cv.perspectiveTransform to make and use the transformation.

For example:
(0, 0) → (333, 333)
(0, 100) → (345, 600)
(100, 0) → (555, 346)
(100, 100) → (610, 600)

So you can see that the shape may not be perfectly square from the four points.

It’s probably not affine, so the equation here will not work. Couldn’t find a similar one for projective transformations.

[baggepinnen]

Perhaps you’ll find what you’re looking for here https://github.com/JuliaGeometry/CoordinateTransformations.jl#perspective-transformations

[vshesh]

Can you show me how to setup the PerspectiveMap? I don’t see how to generate the matrix from the 4 points in that example …

[vshesh]

Ok one solution is:

using OpenCV 

m = reshape(
  OpenCV.getPerspectiveTransform(
    reshape(Array{Float32}([0 0 0 100 100 0 100 100]), 1, 2 ,4), 
    reshape(Array{Float32}([333 333 345 600 345 550 550 625]), 1, 2, 4) , 
    Int32(0)), 
  3, 3)'

But I would prefer to not depend on OpenCV if possible. Also the matrix that OpenCV returns makes no sense to me - when I use it by taking m * [1, 1, 1] i get something outside the range of the coordinates in the new space. If I put in m * [100, 100, 1] I get something with the wrong signs altogether (the y coordinate is negative somehow).

Another option is to solve the system of equations myself and use the resulting formula … this is the 12 equations I’d need to reduce:

Does Julia have some way of solving these large sets of symbolic equations?

[stevengj]

@vshesh

Does Julia have some way of solving these large sets of symbolic equations?

The resulting symbolic formula could be nasty — it might be more efficient simply to construct the matrix and solve it numerically each time you change the parameters.

[vshesh]

Yeah I tried solving the full thing in a wolfram sandbox and the solver timed out.

Solving for converting from an arbitrary coordinate system to a square with side 100 (to simplify some of the equations) gave me this:

Yeah … no hope for implementing that in code.

This is for converting from a square with side = 100:

Somehow this “side” isn’t so bad. Messy to implement directly though.

@stevengl is right - I will need to do this numerically when necessary.

In my application I may need to make and use ~10-20 perspective transforms at 20fps or so. Is there any approach that will be efficient enough to support this? I imagine trying to run rref on that many matrices will be quite slow…

[vshesh]

So this is what making it myself looks like:

using RowEchelon

function perspective_transform(
  ((x1, y1), (x2, y2), (x3, y3), (x4, y4)),  
  ((X1, Y1), (X2, Y2), (X3, Y3), (X4, Y4)))  

  m = [
    x1 y1 1 0 0 0 0 0   -X1 0 0 0 0
    x2 y2 1 0 0 0 0 0   0 -X2 0 0 0
    x3 y3 1 0 0 0 0 0   0 0 -X3 0 0
    x4 y4 1 0 0 0 0 0   0 0 0 -X4 0
    0 0 0 x1 y1 1 0 0   -Y1 0 0 0 0
    0 0 0 x2 y2 1 0 0   0 -Y2 0 0 0
    0 0 0 x3 y3 1 0 0   0 0 -Y3 0 0
    0 0 0 x4 y4 1 0 0   0 0 0 -Y4 0
    0 0 0 0 0 0 x1 y1    -1 0 0 0 -1
    0 0 0 0 0 0 x2 y2    0 -1 0 0 -1
    0 0 0 0 0 0 x3 y3    0 0 -1 0 -1
    0 0 0 0 0 0 x4 y4    0 0 0 -1 -1
  ]
  values = rref(m)[1:8, end]
  [
    values[1] values[2] values[3]
    values[4] values[5] values[6]
    values[7] values[8] 1
  ]
end

Is this the fastest way to set this up? that matrix is 75% zeros so maybe there’s a faster/cleaner way to do this? I’ve also heard that rref isn’t the most numerically stable, so if there’s a better option please tell me.

As for using it I get

function transform_perspective(m::Matrix, (x1, x2)) 
  t = m * [x1 x2 1]'
  [t[1] / t[3], t[2] / t[3]]
end

Of course both of these functions are specialized to the 2d case, but that’s all I need for now.

[baggepinnen]

@vshesh

In my application I may need to make and use ~10-20 perspective transforms at 20fps or so. Is there any approach that will be efficient enough to support this? I imagine trying to run rref on that many matrices will be quite slow

If you preallocate all the matrices such that you do not allocate new ones for each transform application there is a good chance you’ll make 20fps. For example, this

m = [
    x1 y1 ... 

allocates new memory, consider allocating that once outside the function (remember to use const) and just update the already allocated memory. Then make use of in-place multiplications using mul!.

[yakir12]

Do you know about https://github.com/JuliaImages/ImageTransformations.jl? It might have everything you need. And checkout this answer from Tim Holy: Perspective warp an image in Julia - Stack Overflow

Hope this helps!

[vshesh]

Seems like this resource allows someone to apply a transformation matrix (I don’t see a way to solve for one). So I can write the application function like:

using ImageTransformations

transform_perspective(m::Matrix, point::Vector) = PerspectiveMap() ∘ inv(LinearMap(m)) ∘ \x -> push(x, 1)

I’m not clear what I’m gaining compared to the simpler code above. Isn’t this a lot more operations to create all the transformation objects?

[stillyslalom]

Yes, ImageTransformations.jl currently just allows you to apply transformations, not fit them - this would be a nice addition to the package. Your function can be sped up with static arrays (and by rewriting the rref to use \ to solve the linear equation):

using StaticArrays

function perspective_transform_SA(
         ((x1, y1), (x2, y2), (x3, y3), (x4, y4)),
         ((X1, Y1), (X2, Y2), (X3, Y3), (X4, Y4)))

    m = SA_F64[
        x1 y1 1 0 0 0 0 0   -X1 0 0 0
        x2 y2 1 0 0 0 0 0   0 -X2 0 0
        x3 y3 1 0 0 0 0 0   0 0 -X3 0
        x4 y4 1 0 0 0 0 0   0 0 0 -X4
        0 0 0 x1 y1 1 0 0   -Y1 0 0 0
        0 0 0 x2 y2 1 0 0   0 -Y2 0 0
        0 0 0 x3 y3 1 0 0   0 0 -Y3 0
        0 0 0 x4 y4 1 0 0   0 0 0 -Y4
        0 0 0 0 0 0 x1 y1    -1 0 0 0 
        0 0 0 0 0 0 x2 y2    0 -1 0 0 
        0 0 0 0 0 0 x3 y3    0 0 -1 0 
        0 0 0 0 0 0 x4 y4    0 0 0 -1 
    ]
    v = SA_F64[0, 0, 0, 0, 0, 0, 0, 0, -1, -1, -1, -1]
    values = m \ v
    SA[
    values[1] values[2] values[3]
    values[4] values[5] values[6]
    values[7] values[8] 1
    ]
end

julia> @btime perspective_transform(((0, 0), (0, 100), (100, 0), (100, 100)),
                             ((333, 333), (345,600), (555, 346), (610, 600)))
  3.288 μs (59 allocations: 7.95 KiB)
3×3 Matrix{Float64}:
 2.50406      -0.476459    333.0
 0.307087      1.63268     333.0
 0.000511811  -0.00172887    1.0

julia> @btime perspective_transform_SA(((0, 0), (0, 100), (100, 0), (100, 100)),
                             ((333, 333), (345,600), (555, 346), (610, 600)))
  718.519 ns (0 allocations: 0 bytes)
3×3 SMatrix{3, 3, Float64, 9} with indices SOneTo(3)×SOneTo(3):
 2.50406      -0.476459    333.0
 0.307087      1.63268     333.0
 0.000511811  -0.00172887    1.0

At 718 ns per transform, that’s 1.4 million transforms per second. If you only need 400 transforms per second, this won’t be your limiting factor.

[vshesh]

Okay, that definitely solves the problem! Thank you.