Julia version of cv.getPerspectiveTransform?

vshesh · September 21, 2021, 10:33pm

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 · September 22, 2021, 4:00am

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

vshesh · September 22, 2021, 10:10pm

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 · September 22, 2021, 11:14pm

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 · September 23, 2021, 3:37am

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 · September 23, 2021, 4:40am

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 · September 23, 2021, 5:15am

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 · September 23, 2021, 5:37am

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 · September 23, 2021, 8:43am

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 · September 23, 2021, 2:56pm

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 · September 23, 2021, 3:35pm

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 · September 23, 2021, 4:34pm

Okay, that definitely solves the problem! Thank you.