In [1]:

    

import cv2
import matplotlib.pylab as plt
import numpy as np
import PIL
from cStringIO import StringIO


    

In [2]:

    

filename = 'G0113773.JPG'
#define the GCP points in pixels, srcpoints, and the GCP points in RD coordinates, dstpoints
src_points  = np.array([[1312,1048],
                        [2794,796],
                        [1441,1877],
                        [2836,1826],
                        [2056,1440]]).astype(np.float32)
dst_points  = np.array([[85670.316,445806.969],
                        [85700.903,445821.268],
                        [85678.655,445788.205],
                        [85709.590,445799.925],
                        [85689.973,445803.060]]).astype(np.float32)
# slice for performance
s = np.s_[::10,::10]


    

In [3]:

    

testimg = cv2.imread(filename)
m, n, _ = testimg.shape
testimg = cv2.cvtColor(testimg, cv2.COLOR_BGR2RGBA)


    

In [4]:

    

fig, axes = subplots(1,2, figsize=(12,5))
axes[0].imshow(testimg)
_ = axes[0].set_title('Undistorted GoPro image of Architecture square with Photoshop')
axes[0].plot(src_points[:,0], src_points[:,1], 'r--o')
axes[1].plot(dst_points[:,0], dst_points[:,1], 'g--o')
axes[1].axis('equal')
axes[1].set_title('Coordinates in EPSG:28992')


    

    
Out[4]:





<matplotlib.text.Text at 0x2e3fdd0>






    

In [5]:

    

#find the transformation matrix
M,  dc  = cv2.findHomography(src_points,dst_points)
M       = M.astype(float32)
#set up the new coordinate systems from pixels to RD coordinates
# m,n => row, columns -> y, x
Y, X = np.mgrid[:m, :n]


XY = np.concatenate([X[:,:,np.newaxis],Y[:,:,np.newaxis]], axis=2).astype('float32')
RD  = cv2.perspectiveTransform(XY,M)
RDx, RDy = RD[...,0], RD[...,1]


    

In [6]:

    

RDx[s].shape, RDy[s].shape, testimg[s].shape


    

    
Out[6]:





((300, 400), (300, 400), (300, 400, 4))

In [7]:

    

fig, axes = subplots(1,2, figsize=(12,5))
axes[0].imshow(testimg)
_ = axes[0].set_title('Undistorted GoPro image of Architecture square with Photoshop')
axes[0].plot(src_points[:,0], src_points[:,1], 'r--o')
axes[1].plot(dst_points[:,0], dst_points[:,1], 'g--o')
axes[1].axis('equal')
axes[1].set_title('Coordinates in EPSG:28992')
# Pcolormesh does not accept images, so just pass greyscale for now
axes[1].pcolormesh(RDx[s], RDy[s], testimg[s].mean(-1), cmap='Greys_r')


    

    
Out[7]:





<matplotlib.collections.QuadMesh at 0x3264950>






    

In [10]:

    

img = testimg.copy()
# this can be done faster using resize/reshpae
C = np.c_[img[...,0].flat,img[...,1].flat, img[...,2].flat, img[...,3].flat].astype('float32')/255.0


    

In [11]:

    

fig, ax = plt.subplots(figsize=(15,10))

# quadmesh expects one extra cell, so create X,Y in one dimension longer than image size 
Y, X = np.mgrid[:(m+1), :(n+1)]
XY = np.concatenate([X[:,:,np.newaxis],Y[:,:,np.newaxis]], axis=2).astype('float32')
RD  = cv2.perspectiveTransform(XY,M)
RDx, RDy = RD[...,0], RD[...,1]

im = ax.pcolormesh(RDx, RDy, np.zeros(RDx.shape)) 
im.set_array(None) # remove the array
im.set_edgecolor('none')
im.set_facecolor(C)


    

In [9]: