Perspective correction demo

perspective.py

#
#  Copyright (C) 2015 Victor Arribas
#  Copyright (C) 2015 JDE Developers Team
#
#  This program is free software: you can redistribute it and/or modify
#  it under the terms of the GNU General Public License as published by
#  the Free Software Foundation, either version 3 of the License, or
#  (at your option) any later version.
#
#  This program is distributed in the hope that it will be useful,
#  but WITHOUT ANY WARRANTY; without even the implied warranty of
#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#  GNU General Public License for more details.
#
#  You should have received a copy of the GNU General Public License
#  along with this program.  If not, see http://www.gnu.org/licenses/.
#  Authors :
#       Victor Arribas <v.arribas.urjc@gmail.com>

""" This is a ugly test of perspective correction.
 Ground true is a 4 pairs of points taken from OpenCV getPerspectiveTransform (first attempt)
 Here you will found 3 resolution approaches for method based into A x = B
    * Direct resolution of A x = B
    * SVD decomposition
    * Pseudo-inverse

 Epic pitfall:
    `x1 = y1 = [None]*4` implies shared array. So each implemented method was failing once again.
"""

''' OpenCV output:
points= [[[126  51]] [[170  65]] [[167 105]] [[126 101]]]
ref_points= [[[126  51]] [[171  51]] [[171 106]] [[126 106]]]
tf= [[  2.12783461e-01  -1.33259284e-01   6.00299175e+01]
 [ -3.29632744e-01   5.86693388e-01   4.40112909e+01]
 [ -2.46657637e-03  -1.05761337e-03   1.00000000e+00]]
'''

import numpy as np

## Define points
if __name__ == '__main__':
    pO = [] # original points
    pO.append((126,51))
    pO.append((170,65))
    pO.append((167,105))
    pO.append((126,101))

    pD = [] # desired points
    pD.append( (126,51) )
    pD.append( (171,51) )
    pD.append( (171,106) )
    pD.append( (126,106) )


def checkit(H):
    for i in range(0,4):
        p1 = np.append(pO[i],1)
        p2 = np.append(pD[i],1)
        P2 = np.dot(H,p1)
        P2 = P2/P2[2]
        print 'Point',p1,'>>',p2,'?',P2,'-',np.allclose(p2,P2)


if __name__=='__main__':

    A = []
    B = []
    for i in range(0,4):
        (x1,y1) = (pO[i][0], pO[i][1])
        (x2,y2) = (pD[i][0], pD[i][1])

        A.append([x1,y1,1, 0, 0, 0, -x2*x1, -x2*y1])
        A.append([0, 0, 0, x1,y1,1, -y2*x1, -y2*y1])
        B.append(x2)
        B.append(y2)

    print 'A=',A
    print 'B=',B

    print
    print 'PSEUDO-INVERSE METHOD>>'
    Ap = np.linalg.pinv(A)
    H = np.dot(Ap,B)
    H = np.append(H,1).reshape((3,3))
    print 'H=',H
    checkit(H)

    print
    print 'DIRECT METHOD>>'
    #https://docs.scipy.org/doc/numpy/reference/generated/numpy.linalg.solve.html#numpy.linalg.solve
    H = np.linalg.solve(A,B)
    H = np.append(H,1).reshape((3,3))
    print 'H=',H
    checkit(H)


    print
    print 'SVD METHOD>>>'
    A = np.zeros((8,9))
    for i in range(0,4):
        (x1,y1) = (pO[i][0], pO[i][1])
        (x2,y2) = (pD[i][0], pD[i][1])

        A[2*i  ] = [x1,y1,1, 0, 0, 0, -x2*x1, -x2*y1, -x2]
        A[2*i+1] = [0, 0, 0, x1,y1,1, -y2*x1, -y2*y1, -y2]

    print 'A=',A
    U,S,V = np.linalg.svd(A)
    H = V[-1].reshape((3,3))
    H = H/H[2,2]
    print 'H=',H
    checkit(H)


def calculePerspectiveTransform(pO, pD):
    n = pO.shape[0]

    A = np.zeros((2*n,9))
    for i in range(0,n):
        (x1,y1) = (pO[i][0], pO[i][1])
        (x2,y2) = (pD[i][0], pD[i][1])

        A[2*i  ] = [x1,y1,1, 0, 0, 0, -x2*x1, -x2*y1, -x2]
        A[2*i+1] = [0, 0, 0, x1,y1,1, -y2*x1, -y2*y1, -y2]

    U,S,V = np.linalg.svd(A)
    H = V[-1].reshape((3,3))
    H = H/H[2,2]

    return H