In [1]:
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.special as spsc
import scipy.ndimage as spi
import scipy.spatial as spsp
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
from random import uniform
from tqdm import tqdm_notebook as tqdm
from functools import reduce
%matplotlib inline
In [2]:
%%javascript # disable scrolling for output
IPython.OutputArea.prototype._should_scroll = function(lines) {
return false;
}
In [3]:
np.set_printoptions(precision=3)
# make matplotlib figures larger
matplotlib.rcParams['figure.figsize'] = [10, 10]
Load Images¶
In [4]:
im_A = skio.imread('./ims/mid_1.jpg')
im_B = skio.imread('./ims/mid_2.jpg')
Step 1a: Detect Harris Corners¶
In [5]:
from skimage.feature import corner_harris, peak_local_max
def get_harris_corners(im, edge_discard=20):
"""
This function takes a b&w image and an optional amount to discard
on the edge (default is 5 pixels), and finds all harris corners
in the image. Harris corners near the edge are discarded and the
coordinates of the remaining corners are returned. A 2d array (h)
containing the h value of every pixel is also returned.
h is the same shape as the original image, im.
coords is 2 x n (ys, xs).
"""
assert edge_discard >= 20
# find harris corners
h = corner_harris(im, method='eps', sigma=1)
coords = peak_local_max(h, min_distance=1, indices=True, threshold_rel=0.02)
# discard points on edge
edge = edge_discard # pixels
mask = (coords[:, 0] > edge) & \
(coords[:, 0] < im.shape[0] - edge) & \
(coords[:, 1] > edge) & \
(coords[:, 1] < im.shape[1] - edge)
coords = coords[mask].T
return h, coords
def dist2(x, c):
"""
dist2 Calculates squared distance between two sets of points.
Description
D = DIST2(X, C) takes two matrices of vectors and calculates the
squared Euclidean distance between them. Both matrices must be of
the same column dimension. If X has M rows and N columns, and C has
L rows and N columns, then the result has M rows and L columns. The
I, Jth entry is the squared distance from the Ith row of X to the
Jth row of C.
Adapted from code by Christopher M Bishop and Ian T Nabney.
"""
ndata, dimx = x.shape
ncenters, dimc = c.shape
assert(dimx == dimc, 'Data dimension does not match dimension of centers')
return (np.ones((ncenters, 1)) * np.sum((x**2).T, axis=0)).T + \
np.ones(( ndata, 1)) * np.sum((c**2).T, axis=0) - \
2 * np.inner(x, c)
Compute¶
In [6]:
im_A_corners_vals, im_A_corners_coords = get_harris_corners(sk.color.rgb2gray(im_A), edge_discard=80)
im_B_corners_vals, im_B_corners_coords = get_harris_corners(sk.color.rgb2gray(im_B), edge_discard=80)
im_A_corners_coords = im_A_corners_coords.T
im_B_corners_coords = im_B_corners_coords.T
Visualize¶
In [7]:
RADIUS = 2
im_A_corners_overlay = np.copy(im_A)
for p in im_A_corners_coords:
for x_off in range(-RADIUS,RADIUS+1):
for y_off in range(-RADIUS,RADIUS+1):
im_A_corners_overlay[p[0]+y_off, p[1]+x_off] = (255, 0, 0)
im_B_corners_overlay = np.copy(im_B)
for p in im_B_corners_coords:
for x_off in range(-RADIUS,RADIUS+1):
for y_off in range(-RADIUS,RADIUS+1):
im_B_corners_overlay[p[0]+y_off, p[1]+x_off] = (255, 0, 0)
In [8]:
skio.imshow(im_A_corners_overlay)
plt.show()
skio.imsave('./out/im_A_corners_overlay.jpg', im_A_corners_overlay, quality=100)
skio.imshow(im_B_corners_overlay)
plt.show()
skio.imsave('./out/im_B_corners_overlay.jpg', im_B_corners_overlay, quality=100)
Step 1b: Adaptive Non-Maximal Suppression¶
In [9]:
im_A_corners_coords.shape
Out[9]:
In [10]:
def ANMS(pts, H, num_pts=500, c_robust = 0.9):
new_pts = np.array([[0, 0], 0])
for p_i in tqdm(pts):
strength_i = H[p_i[0], p_i[1]]
ops = pts[ c_robust * H[ pts[:,0],pts[:,1] ] > strength_i ]
if ops.shape[0] == 0:
new_pts = np.vstack((new_pts, np.array([p_i, float('inf')])))
else:
dists = dist2(np.array([p_i]), np.array(ops))
new_pts = np.vstack( ( new_pts,np.array([p_i, np.amin(dists)]) ) )
new_pts = new_pts[1:]
sort_i = np.argsort(new_pts[:,1])[-num_pts:]
return np.stack(new_pts[sort_i][:num_pts,0])
Compute¶
In [11]:
im_A_corners_coords_anms = ANMS(im_A_corners_coords, im_A_corners_vals)
im_B_corners_coords_anms = ANMS(im_B_corners_coords, im_B_corners_vals)
Visualize¶
In [12]:
RADIUS = 5
im_A_corners_overlay_anms = np.copy(im_A)
for p in im_A_corners_coords_anms:
for x_off in range(-RADIUS,RADIUS+1):
for y_off in range(-RADIUS,RADIUS+1):
im_A_corners_overlay_anms[p[0]+y_off, p[1]+x_off] = (255, 0, 0)
im_B_corners_overlay_anms = np.copy(im_B)
for p in im_B_corners_coords_anms:
for x_off in range(-RADIUS,RADIUS+1):
for y_off in range(-RADIUS,RADIUS+1):
im_B_corners_overlay_anms[p[0]+y_off, p[1]+x_off] = (255, 0, 0)
In [13]:
skio.imshow(im_A_corners_overlay_anms)
plt.show()
skio.imsave('./out/im_A_corners_overlay_anms.jpg', im_A_corners_overlay_anms, quality=100)
skio.imshow(im_B_corners_overlay_anms)
plt.show()
skio.imsave('./out/im_B_corners_overlay_anms.jpg', im_B_corners_overlay_anms, quality=100)
Step 2: Extract Feature Descriptors¶
In [14]:
def get_feature_descriptors(im, pts, sample_res=40, output_res=8):
feat_descs = []
for p in pts:
sample_raw = im[p[0]-(sample_res//2):p[0]+(sample_res//2), p[1]-(sample_res//2):p[1]+(sample_res//2)]
sample_downscale = cv2.resize(sample_raw, (output_res,output_res))
sample_downscale_norm = (sample_downscale - sample_downscale.mean()) / sample_downscale.std()
sample_downscale_norm_bw = sk.color.rgb2gray(sample_downscale_norm)
feat_descs.append(sample_downscale_norm_bw)
return np.array(feat_descs)
Compute¶
In [15]:
im_A_feature_descriptors = get_feature_descriptors(im_A, im_A_corners_coords_anms, sample_res=160, output_res=16)
im_B_feature_descriptors = get_feature_descriptors(im_B, im_B_corners_coords_anms, sample_res=160, output_res=16)
Visualize¶
In [16]:
# n/a
Step 3: Match Feature Descriptors¶
In [17]:
def match_feature_descriptors(feat_desc_A, feat_desc_B, threshold=0.65):
assert feat_desc_A.shape == feat_desc_B.shape
num_feats, feat_size, _ = feat_desc_A.shape
feat_desc_vectorized_A = np.resize(feat_desc_A, (num_feats,feat_size*feat_size))
feat_desc_vectorized_B = np.resize(feat_desc_B, (num_feats,feat_size*feat_size))
dists = dist2(feat_desc_vectorized_A, feat_desc_vectorized_B)
e_1nn = np.amin(dists, axis=1)
nn_j = np.argmin(dists, axis=1)
e_2nn = np.partition(dists, 1)[:,1]
nn_i = np.nonzero((e_1nn/e_2nn) < threshold)[0]
feature_matchings = (nn_i, nn_j[nn_i])
return feature_matchings
Compute¶
In [18]:
feature_matchings = match_feature_descriptors(im_A_feature_descriptors, im_B_feature_descriptors)
In [19]:
num_feats = feature_matchings[0].shape[0]
print("Number of Matches = {}".format(num_feats))
Visualize¶
In [20]:
RADIUS = 25
im_A_matchings_overlay = np.copy(im_A)
im_B_matchings_overlay = np.copy(im_B)
rand_color = lambda : int(uniform(0,255))
for i in range(num_feats):
pa_i = feature_matchings[0][i]
pb_i = feature_matchings[1][i]
pa = im_A_corners_coords_anms[pa_i]
pb = im_B_corners_coords_anms[pb_i]
color = np.array([rand_color(), rand_color(), rand_color()])
for x_off in range(-RADIUS,RADIUS+1):
for y_off in range(-RADIUS,RADIUS+1):
im_A_matchings_overlay[pa[0]+y_off, pa[1]+x_off] = color
im_B_matchings_overlay[pb[0]+y_off, pb[1]+x_off] = color
In [21]:
skio.imshow(im_A_matchings_overlay)
plt.show()
skio.imsave('./out/im_A_matchings_overlay.jpg', im_A_matchings_overlay, quality=100)
skio.imshow(im_B_matchings_overlay)
plt.show()
skio.imsave('./out/im_B_matchings_overlay.jpg', im_B_matchings_overlay, quality=100)
Step 4: RANSAC for Homography¶
In [22]:
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 x.reshape(3,3)
In [23]:
def ransac(pts_A, pts_B, feature_matchings, check_perc=2, threshold=80):
num_pts = feature_matchings[0].shape[0]
epochs = int(check_perc * spsc.comb(num_pts,4))
pts_A_matched = pts_A[feature_matchings[0]] # (n, 2)
pts_B_matched = pts_B[feature_matchings[1]] # (n, 2)
best_inliers_i = None
best_num_inliers = float('-inf')
best_anchor_pts_A = None
best_anchor_pts_B = None
for _ in tqdm(range(epochs)):
pts_i = []
while len(pts_i) < 4:
i = int(uniform(0, num_pts-1))
if i not in pts_i:
pts_i.append(i)
pts_i = np.array(sorted(pts_i)) # (1, n)
pts_A_subset = pts_A_matched[pts_i] # (k, 2)
pts_B_subset = pts_B_matched[pts_i] # (k, 2)
H = computeH(pts_A_subset, pts_B_subset) # (3, 3) [x (2/3, k)]
b_prep = np.hstack( (pts_B_matched, np.ones( (num_pts,1) )) ) # (n, 3)
Hb = np.matmul(H, b_prep.T)[:2,] # (2, n)
if False:
b_prep = np.hstack( (pts_B_subset, np.ones( (4,1) )) ) # (4, 3)
test = np.matmul(H, b_prep.T) # (2, 4)
print('bprep=')
print(b_prep)
print('test=')
print(test)
print('A=')
print(pts_A_subset.T)
print('diff')
print(test[:2] - pts_A_subset.T)
error = np.sum((test[:2] - pts_A_subset.T)**2, axis=0) # (1, 4)
print('ssdvec={}'.format(error))
error = np.sum(np.sum((test[:2] - pts_A_subset.T)**2, axis=0)) # (1, 4)
print('ssd={}'.format(error))
return None
ssd = np.sum((Hb - pts_A_matched.T)**2, axis=0) # (1, n)
inliers_i = np.nonzero(ssd < threshold)
num_inliers = np.count_nonzero(ssd < threshold)
if num_inliers > best_num_inliers:
best_inliers_i = inliers_i
best_num_inliers = num_inliers
best_anchor_pts_A = pts_A_subset
best_anchor_pts_B = pts_B_subset
print('Found Homography with {} Inliers'.format(num_inliers))
pts_A_matched_inliers = pts_A_matched[best_inliers_i]
pts_B_matched_inliers = pts_B_matched[best_inliers_i]
H = computeH(pts_A_matched_inliers, pts_B_matched_inliers)
return H, pts_A_matched_inliers, pts_B_matched_inliers, best_anchor_pts_A, best_anchor_pts_B
Compute¶
In [24]:
H, pAmi, pBmi, bapA, bapB = ransac(im_A_corners_coords_anms, im_B_corners_coords_anms, feature_matchings)
Visualize¶
In [25]:
RADIUS = 25
im_A_ransac_overlay = np.copy(im_A)
im_B_ransac_overlay = np.copy(im_B)
im_A_ransac_overlay_warped_pts = np.copy(im_A)
rand_color = lambda : int(uniform(0,255))
for i in range(pAmi.shape[0]):
pa = pAmi[i]
pb = pBmi[i]
pc_prep = np.hstack( (pb, 1) ).T
pc = np.rint(np.matmul(H, pc_prep)).astype(int)[:2]
color = np.array([0, rand_color(), rand_color()])
for x_off in range(-RADIUS,RADIUS+1):
for y_off in range(-RADIUS,RADIUS+1):
im_A_ransac_overlay[pa[0]+y_off, pa[1]+x_off] = color
im_B_ransac_overlay[pb[0]+y_off, pb[1]+x_off] = color
im_A_ransac_overlay_warped_pts[pc[0]+y_off, pc[1]+x_off] = color
rand_color = lambda : int(uniform(100,255))
for i in range(bapA.shape[0]):
pa = bapA[i]
pb = bapB[i]
pc_prep = np.hstack( (pb, 1) ).T
pc = np.rint(np.matmul(H, pc_prep)).astype(int)[:2]
color = np.array([rand_color(), 0, 0])
for x_off in range(-RADIUS,RADIUS+1):
for y_off in range(-RADIUS,RADIUS+1):
im_A_ransac_overlay[pa[0]+y_off, pa[1]+x_off] = color
im_B_ransac_overlay[pb[0]+y_off, pb[1]+x_off] = color
im_A_ransac_overlay_warped_pts[pc[0]+y_off, pc[1]+x_off] = color
In [26]:
print('Actual Points im A')
skio.imshow(im_A_ransac_overlay)
plt.show()
skio.imsave('./out/im_A_ransac_overlay.jpg', im_A_ransac_overlay, quality=100)
print('Actual Points im B')
skio.imshow(im_B_ransac_overlay)
plt.show()
skio.imsave('./out/im_B_ransac_overlay.jpg', im_B_ransac_overlay, quality=100)
print('Warped Points im B -> A')
skio.imshow(im_A_ransac_overlay_warped_pts)
plt.show()
skio.imsave('./out/im_A_ransac_overlay_warped_pts.jpg', im_A_ransac_overlay_warped_pts, quality=100)
Step 5a: Warp into Mosaic¶
In [27]:
def warpImage(im, H_mat):
# Warp Corners
im_B_H, im_B_W, im_B_D = im.shape
im_corners = np.array([(0, 0, 1),
(im_B_H, 0, 1),
(im_B_H, im_B_W, 1),
(0, im_B_W, 1)
])
new_corners = np.matmul(H_mat, im_corners.T)
# Fix Bounding Box
HMAX, HMIN = max(new_corners[0]), min(new_corners[0])
WMAX, WMIN = max(new_corners[1]), min(new_corners[1])
warp_H, warp_W = int(HMAX), int(WMAX)
t_vec = np.zeros(2)
if HMIN < 0:
t_vec += np.array([-1 * int(HMIN) + 1, 0])
if WMIN < 0:
t_vec += np.array([0, -1 * int(WMIN) + 1])
t_vec = np.rint(t_vec).astype(int)
print('im_B dims = {}'.format(im.shape[:2]))
print('t_vec = {}'.format(t_vec))
print('warped im size = {}'.format((warp_H + t_vec[0], warp_W + t_vec[1])))
print('corners (TL BL BR TR) =')
print(np.rint(new_corners).astype(int)[:2].T)
# Points from warped canvas to be Calculated
im_points_rr, im_points_cc = skd.polygon(new_corners[0]+t_vec[0], new_corners[1]+t_vec[1], (warp_H + t_vec[0], warp_W + t_vec[1]))
im_points_rr -= t_vec[0]
im_points_cc -= t_vec[1]
im_points_one = np.column_stack((im_points_rr, im_points_cc, np.ones(im_points_cc.shape[0])))
# Find Points from Image B to be Used
raw_pts = np.matmul(np.linalg.inv(H_mat), im_points_one.T)
im_pts = np.rint(raw_pts[:2]).astype(int).T
good_i_1 = np.nonzero(im_pts[:,0] >= 0)
good_i_2 = np.nonzero(im_pts[:,0] < im_B_H)
good_i_3 = np.nonzero(im_pts[:,1] >= 0)
good_i_4 = np.nonzero(im_pts[:,1] < im_B_W)
good_i = reduce(np.intersect1d, (good_i_1, good_i_2, good_i_3, good_i_4))
#print(good_i)
im_pts = im_pts[good_i]
im_points_rr = im_points_rr[good_i]
im_points_cc = im_points_cc[good_i]
#print(np.amin(im_pts[:,0]))
#print(np.amax(im_pts[:,0]))
#print(np.amin(im_pts[:,1]))
#print(np.amax(im_pts[:,1]))
# Warp Points
im_warped = np.zeros((warp_H + t_vec[0], warp_W + t_vec[1], im_B_D), im.dtype)
im_warped[im_points_rr+t_vec[0], im_points_cc+t_vec[1]] = im[im_pts[:,0],[im_pts[:,1]]]
return im_warped, t_vec, np.rint(new_corners).astype(int)[:2].T
Compute¶
In [28]:
im_B_warped, translation_vec, im_B_warped_corners = warpImage(im_B, H)
Visualize¶
In [47]:
skio.imshow(im_B_warped)
plt.show()
skio.imsave('./out/im_B_warped.jpg', im_B_warped, quality=100)
Step 5b: Combine Images and Gaussian+Laplacian Spectrum Blending (NOT USED)¶
In [48]:
def blend(im_1, im_2, mask, N):
assert im_1.shape == im_2.shape == mask.shape
im_curr_gauss_1 = im_1
im_curr_gauss_2 = im_2
mask_curr = mask
im_combined = np.zeros(mask.shape)
# Loop through Levels of Stack
for n in range(0, N + 1):
# Gaussian Constants
k = 100
s = k/3
# Create Gaussian Filter
gaussian_1d = cv2.getGaussianKernel(k, s)
gaussian_2d = np.outer(gaussian_1d, gaussian_1d)
gaussian_2d = gaussian_2d / np.sum(gaussian_2d)
# Apply Gaussian
im_curr_gauss_prev_1 = im_curr_gauss_1
im_curr_gauss_prev_2 = im_curr_gauss_2
im_curr_gauss_1 = spi.convolve(im_curr_gauss_1, gaussian_2d)
im_curr_gauss_2 = spi.convolve(im_curr_gauss_2, gaussian_2d)
mask_curr = spi.convolve(mask_curr, gaussian_2d)
# Find Laplacian
if n == N:
im_curr_lapl_1 = im_curr_gauss_prev_1
im_curr_lapl_2 = im_curr_gauss_prev_2
else:
im_curr_lapl_1 = im_curr_gauss_prev_1 - im_curr_gauss_1
im_curr_lapl_2 = im_curr_gauss_prev_2 - im_curr_gauss_2
# Construct Combined Image
im_combined_curr = np.multiply(mask_curr, im_curr_lapl_1) + np.multiply(1 - mask_curr, im_curr_lapl_2)
im_combined = im_combined + im_combined_curr
# Show Current Stack Layer
print_img('Laplacian Stack N={}'.format(n))
skio.imshow(im_combined_curr)
plt.show()
return im_combined
Create Mask¶
In [ ]:
Compute¶
In [49]:
"""
# Calculate Blending
im_blend = blend(im_1, im_B_warped, border_mask, 3)
"""
Out[49]:
Step 5b: Combine Images and Alpha Blend¶
In [50]:
canvas = np.zeros(im_B_warped.shape)
In [51]:
H, W, D = im_A.shape
rr = np.array([0, H, H, 0])
cc = np.array([0, 0, W, W])
im_points_rr_old, im_points_cc_old = skd.polygon(rr, cc)
im_points_rr_new, im_points_cc_new = im_points_rr_old + translation_vec[0], im_points_cc_old + translation_vec[1]
good_i_1 = np.nonzero(im_points_rr_old >= 0)
good_i_2 = np.nonzero(im_points_rr_old < H)
good_i_3 = np.nonzero(im_points_cc_old >= 0)
good_i_4 = np.nonzero(im_points_cc_old < W)
good_i = reduce(np.intersect1d, (good_i_1, good_i_2, good_i_3, good_i_4))
im_points_cc_old = im_points_cc_old[good_i]
im_points_rr_old = im_points_rr_old[good_i]
im_points_cc_new = im_points_cc_new[good_i]
im_points_rr_new = im_points_rr_new[good_i]
canvas[im_points_rr_new, im_points_cc_new] = im_A[im_points_rr_old,im_points_cc_old]
In [52]:
mask = np.logical_and(canvas, im_B_warped)
canvas += im_B_warped
canvas[mask > 0] = canvas[mask > 0] / 2
canvas = canvas.astype(np.uint8)
Visualize¶
In [53]:
skio.imshow(canvas)
plt.show()
skio.imsave('./out/final_result.jpg', canvas, quality=100)
In [ ]: