In [9]:
import numpy as np
import cv2
from matplotlib import pyplot as plt
import random
import math
import warnings
warnings.filterwarnings("ignore")
def main(I):
img1 = cv2.imread(I[0],0) # queryImage
img2 = cv2.imread(I[1],0) # trainImage
# Initiate SIFT detector
sift = cv2.SIFT_create()
# find the keypoints and descriptors with SIFT
kp1, des1 = sift.detectAndCompute(img1,None)
kp2, des2 = sift.detectAndCompute(img2,None)
# BFMatcher with default params
bf = cv2.BFMatcher()
matches = bf.knnMatch(des1,des2, k=2)
# Apply ratio test
good = []
for m,n in matches:
if m.distance < 0.7*n.distance:
good.append([m])
# cv2.drawMatchesKnn expects list of lists as matches.
img3 = cv2.drawMatchesKnn(img1,kp1,img2,kp2,good,None,flags=2) #NOTE: 'None' parameter has to be added (not in documentation)
plt.figure(figsize=(12,8))
plt.imshow(img3),plt.show()
def MatrixMul(m1,m2):
multiplication = [[sum(a * b for a, b in zip(X_row, Y_col)) for Y_col in zip(*m2)] for X_row in m1]
return multiplication
#RANSAC Loop
#Variables to be updated to find the bestH
lowestsum = 10000000000
bestH = []
for iteration in range(5000):
#Picking 4 random indices in the list of good matches
randompts = random.sample(range(0,len(good)), 4)
a = []
#Picking 4 pairs to create matrix a
#this loop will used the 4 calculated random indices above and create the 8x9 matrix a
for index in range(0,4):
i = randompts[index]
qIdx = good[i][0].queryIdx
tIdx = good[i][0].trainIdx
x1 = kp1[qIdx].pt[0]
y1 = kp1[qIdx].pt[1]
x2 = kp2[tIdx].pt[0]
y2 = kp2[tIdx].pt[1]
#Each of the matching pair = 2 rows of the 8 total row of the matrix
a.append([0, 0, 0, -x1, -y1, -1, y2*x1, y2*y1, y2])
a.append([x1, y2, 1, 0, 0, 0, -x2*x1, -x2*y1, -x2])
#RANSAC
U, s, V = np.linalg.svd(a, full_matrices=True)
hCol = np.zeros((9,1), np.float64)
hCol = V[8,:]
# Creating h
hMatrix = [[hCol[0], hCol[1], hCol[2]], [hCol[3], hCol[4], hCol[5]], [hCol[6], hCol[7], hCol[8]]]
totalsum = 0
# Calculating sum of distance
for pairs in range(len(good)):
qIdx = good[pairs][0].queryIdx
tIdx = good[pairs][0].trainIdx
x1 = kp1[qIdx].pt[0]
y1 = kp1[qIdx].pt[1]
x2 = kp2[tIdx].pt[0]
y2 = kp2[tIdx].pt[1]
pt1Matrix = [[x1], [y1], [1]]
result = MatrixMul(hMatrix,pt1Matrix)
#Dividing all 3 rows by homogenous number
for item in range(len(result)):
result[item][0] = result[item][0] / result[2][0]
totalsum += abs(math.sqrt(((result[0][0] - x2) ** 2) + ((result[1][0] - y2) ** 2)))
if totalsum < lowestsum:
lowestsum = totalsum
bestH = hMatrix
print(bestH)
#Panorama Size
mc1 = [[0], [0], [1]]
mc1 = MatrixMul(bestH,mc1)
for item in range(len(mc1)):
mc1[item][0] = mc1[item][0] / mc1[2][0]
mc2 = [[img1.shape[1]], [0], [1]]
mc2 = MatrixMul(bestH,mc2)
for item in range(len(mc2)):
mc2[item][0] = mc2[item][0] / mc2[2][0]
mc3 = [[0], [img1.shape[0]], [1]]
mc3 = MatrixMul(bestH,mc3)
for item in range(len(mc3)):
mc3[item][0] = mc3[item][0] / mc3[2][0]
mc4 = [[img1.shape[1]], [img1.shape[0]], [1]]
mc4 = MatrixMul(bestH,mc4)
for item in range(len(mc4)):
mc4[item][0] = mc4[item][0] / mc4[2][0]
minX = round(min(0, mc1[0][0], mc3[0][0]))
maxX = round(max(img2.shape[1], mc2[0][0], mc4[0][0]))
minY = round(min(0, mc1[1][0], mc2[1][0]))
maxY = round(max(img2.shape[0], mc3[1][0], mc4[1][0]))
#Creating blank Panorama and re-reading colored version of images
panorama = np.zeros(((maxY-minY),(maxX-minX), 3), np.float32)
repeatedpixel = np.zeros(((maxY-minY),(maxX-minX)), np.float32)
img1pix = np.zeros(((maxY-minY),(maxX-minX)), np.float32)
totalred = 0
totalgreen = 0
totalblue = 0
totalpix = 0
img1 = cv2.imread(I[0])
img2 = cv2.imread(I[1])
shiftx = min(0, minX)
shifty = min(0, minY)
#placing img2 into panorama
for row in range(img2.shape[0]-1):
for col in range(img2.shape[1]-1):
panorama[row+abs(shifty)][col+abs(shiftx)][0] = img2[row][col][0]
panorama[row+abs(shifty)][col+abs(shiftx)][1] = img2[row][col][1]
panorama[row+abs(shifty)][col+abs(shiftx)][2] = img2[row][col][2]
repeatedpixel[row+abs(shifty)][col+abs(shiftx)] = 100
#stitching phase
for row in range(panorama.shape[0]):
for col in range(panorama.shape[1]):
# multiplying each pixel by inverse H and dividing by homogenous value
Hinv = np.linalg.inv(bestH)
temp = [[col + shiftx],[row + shifty], [1]]
mappedpt1 = MatrixMul(Hinv, temp)
for item in range(len(mappedpt1)):
mappedpt1[item][0] = mappedpt1[item][0] / mappedpt1[2][0]
roundx = round(mappedpt1[0][0])
roundy = round(mappedpt1[1][0])
# If calculated value of x and y in boundaries of img1, copy to panorama
if 0 < roundy < img1.shape[0] and 0 < roundx < img1.shape[1]:
if repeatedpixel[row][col] > 50:
totalblue += panorama[row][col][2] - img1[roundy][roundx][2]
totalgreen += panorama[row][col][1] - img1[roundy][roundx][1]
totalred += panorama[row][col][0] - img1[roundy][roundx][0]
repeatedpixel[row][col] = 200
panorama[row][col][0] = img1[roundy][roundx][0]
panorama[row][col][1] = img1[roundy][roundx][1]
panorama[row][col][2] = img1[roundy][roundx][2]
img1pix[row][col] = 100
totalpix += 1
totalred = totalred / totalpix
totalgreen = totalgreen / totalpix
totalblue = totalblue / totalpix
#image blending
for row in range(repeatedpixel.shape[0]):
for col in range(repeatedpixel.shape[1]):
if img1pix[row][col] > 50:
panorama[row][col][0] += totalblue
panorama[row][col][1] += totalgreen
panorama[row][col][2] += totalred
panorama = cv2.cvtColor(panorama, cv2.COLOR_RGB2BGR)
plt.figure(figsize=(12,8))
plt.imshow(panorama / 255.0),plt.show()
if __name__ == "__main__":
waterfall = ('./inputs/waterfall1.jpg', './inputs/waterfall2.jpg', 'waterfall')
mountain = ('./inputs/mountainLeft.jpg', './inputs/mountainRight.jpg', 'mountain')
florence = ('./inputs/florenceLeft.jpg', './inputs/florenceRight.jpg', 'florence')
queenstown = ('./inputs/queenstownLeft.jpg', './inputs/queenstownRight.jpg', 'queenstown')
main(waterfall)
main(mountain)
main(florence)
main(queenstown)
[[0.0018994347138088417, 3.4953667426202495e-06, -0.7302152545041634], [-7.093636555920987e-05, 0.0018498829577031782, -0.6832116619027682], [8.798819151990983e-09, 2.8689743369393095e-08, 0.0006871434363526835]]
Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).
[[-0.007557390067003188, 0.00019841994839076834, 0.9992503981730043], [-0.0007068878836423777, -0.006960985605317211, 0.03680230881080681], [-4.657502704080183e-06, 9.61055667475412e-07, -0.0061744043025477955]]
Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).
[[-0.011580961089875288, 0.0017236247943429584, 0.9244139934152562], [-7.638272189096399e-05, -0.012870734865743869, 0.38086819841251396], [-1.2342911133354924e-06, -6.715906131430786e-06, -0.009768980342758034]]
Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).
[[-0.004150859567082871, -0.0006210051005411741, 0.999496190322013], [0.00024229396018350563, -0.003994395984906326, -0.031039559206188748], [-2.289834529982891e-06, -4.2190939774314937e-07, -0.003206571770345762]]
Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).
