import numpy as np
import skimage as sk
import skimage.draw as skd
import skimage.io as skio
import skimage.transform as sktr
import scipy.ndimage as spi
import scipy.spatial as sps
import scipy.interpolate as spin
import matplotlib
matplotlib.use('TkAgg')
import matplotlib.pyplot as plt
import pandas as pd
import cv2
import time
import math

%matplotlib inline

%%javascript # disable scrolling for output
IPython.OutputArea.prototype._should_scroll = function(lines) {
    return false;
}

# make matplotlib figures larger
matplotlib.rcParams['figure.figsize'] = [10, 10]

Load Images

im_main = skio.imread('./ims/backyard_left.jpg')
im_right = skio.imread('./ims/backyard_right.jpg')

Recover Homographies

def computeH(im1_pts, im2_pts):
    pts_1 = np.hstack((im1_pts, np.ones((im1_pts.shape[0],1))))
    pts_2 = np.hstack((im2_pts, np.ones((im2_pts.shape[0],1))))
    
    y = np.zeros(1)
    for p in pts_1:
        p = p.reshape((3,1))
        y = np.vstack((y,p))
    y = y[1:]
    
    A = np.zeros((1,9))
    for p in pts_2:
        r = np.hstack((p, np.zeros(6)))
        A = np.vstack((A, r))
        r = np.hstack((np.zeros(3), p, np.zeros(3)))
        A = np.vstack((A, r))
        r = np.hstack((np.zeros(6), p))
        A = np.vstack((A, r))
    A = A[1:,]
    
    x, r, rank, _ = np.linalg.lstsq(A, y, rcond=None)
    
    print('residuals={}'.format(r))
    print('rank={}'.format(rank))
    
    return np.reshape(x,(3,3))

%matplotlib qt

NUM_POINTS = 7

plt.imshow(im_main)
pts_main = plt.ginput(NUM_POINTS, timeout=-1)
plt.close()
pts_main = np.array(pts_main)
pts_main = pts_main[1:,]

plt.imshow(im_right)
pts_right = plt.ginput(NUM_POINTS, timeout=-1)
plt.close()
pts_right = np.array(pts_right)
pts_right = pts_right[1:,]

%matplotlib inline

H = computeH(pts_main, pts_right)
print(H)

residuals=[8030.60095437]
rank=9
[[ 9.94987087e-01 -4.50759034e-02  1.13514846e+03]
 [ 1.31141194e-01  9.75187735e-01  1.31858212e+02]
 [-1.54374785e-13 -2.98333491e-14  1.00000000e+00]]

d = np.matmul(H, np.hstack((pts_right,np.ones((6,1)))).T).T - np.hstack((pts_main,np.ones((6,1))))
for p in d:
    print('x: {}    y: {}'.format(p[0], p[1]))

x: -26.688596085044082    y: 13.107334516294486
x: 53.44590283863454    y: 28.527873797912434
x: -39.1869095707616    y: 6.035253389310128
x: 17.47017162010934    y: -30.375839067824927
x: 13.505294653174587    y: -8.20982998183672
x: -18.545863456110737    y: -9.08479265410051

Warp the Images

def warpImage(im, H_mat):
    
    # Warp Corners
    H, W, D = im.shape
    im_corners = np.array([(0, 0, 1),
                           (W, 0, 1),
                           (W, H, 1),
                           (0, H, 1)
                          ])
    new_corners = np.matmul(H_mat, im_corners.T)
    
    # Fix Bounding Box
    HMAX, HMIN = max(new_corners[1]), min(new_corners[1])
    WMAX, WMIN = max(new_corners[0]), min(new_corners[0])
    H, W = int(HMAX), int(WMAX)
    
    # Points from Image to be Calculated
    im_points_rr, im_points_cc = skd.polygon(new_corners[1], new_corners[0], (H, W))
    im_points_one = np.column_stack((im_points_cc, im_points_rr, np.ones(im_points_cc.shape[0])))

    # Find Points from Image to be Used
    im_pts = np.matmul(np.linalg.inv(H_mat), im_points_one.T).T[:,:2].astype(int)

    # Warp Points
    im_warped = np.zeros((H, W, D), im.dtype)
    im_warped[im_points_rr, im_points_cc] = im[im_pts[:,1],[im_pts[:,0]]]
    
    return im_warped

warped = warpImage(im_right, H)

skio.imshow(warped)
plt.show()
skio.imsave('./out/warped.jpg', warped, quality=100)

Blend the Images into a Mosaic|

canvas = np.zeros(warped.shape)

warped_mask = np.zeros(canvas.shape)
warped_mask[:warped.shape[0],:warped.shape[1]] = (warped != 0)
main_mask = np.zeros(canvas.shape)
main_mask[:im_main.shape[0],:im_main.shape[1]] = (im_main != 0)

shared_mask = np.logical_and(main_mask, warped_mask)

skio.imshow(shared_mask.astype(int))
plt.show()

canvas += warped
canvas[:im_main.shape[0], :im_main.shape[1]] += im_main

canvas[:,:] = canvas[:,:] / (1 + (shared_mask[:,:] == 1))

canvas = canvas.astype(np.uint8)

skio.imshow(canvas)
plt.show()
skio.imsave('./out/combined.jpg', canvas, quality=100)

C:\Users\Atte\Anaconda3\lib\site-packages\skimage\io\_plugins\matplotlib_plugin.py:77: UserWarning: Low image data range; displaying image with stretched contrast.
  warn("Low image data range; displaying image with "

Image Rectification

NUM_POINTS = 5

im_main = skio.imread('./ims/desktop.jpg')

%matplotlib qt

plt.imshow(im_main)
pts_main_rect = plt.ginput(NUM_POINTS, timeout=-1)
plt.close()
pts_main_rect = np.array(pts_main_rect[1:])

%matplotlib inline

rect_pts = np.array([
    [0, 0],
    [0, 2000], 
    [2000, 2000],
    [2000, 0]
])

H_rect = computeH(rect_pts, pts_main_rect)

rectified = warpImage(im_main, H_rect)

skio.imshow(rectified)
plt.show()
skio.imsave('./out/rectified.jpg', rectified, quality=100)

residuals=[39923.64545301]
rank=9