import numpy as np
import matplotlib.pyplot as plt
import exif
from IPython.core.display import HTML
HTML("""
<style>
div.cell { /* Tunes the space between cells */
margin-top:1em;
margin-bottom:1em;
}
div.text_cell_render h1 { /* Main titles bigger, centered */
font-size: 2.2em;
line-height:0.9em;
}
div.text_cell_render h2 { /* Parts names nearer from text */
margin-bottom: -0.4em;
}
div.text_cell_render { /* Customize text cells */
font-family: 'Georgia';
font-size:1.2em;
line-height:1.4em;
padding-left:3em;
padding-right:3em;
}
.output_png {
display: table-cell;
text-align: center;
vertical-align: middle;
}
</style>
<script>
code_show=true;
function code_toggle() {
if (code_show){
$('div.input').hide();
} else {
$('div.input').show();
}
code_show = !code_show
}
$( document ).ready(code_toggle);
</script>
""")
Recovering Homographies¶
Before we can warp our images into alignment, we need to recover the parameters of the transformation between each pair of images. In our case, the transformation is a homography: $p’=Hp$, where $H \in R^{3\times3}$ matrix with 8 degrees of freedom (lower right corner is a scaling factor and can be set to 1). One way to recover the homography is via a set of $(p’,p)$ pairs of corresponding points taken from the two images.
We collected 20 points for our mapping, and then constructed an overdetermined system of equations in which we solve for our 8 homographic parameters using least squares.
For one pair of points, we can write the transform as $$p'=Hp$$ $$\begin{bmatrix}wx' \\ wy' \\ w\end{bmatrix} = \begin{bmatrix}a & b & c\\ c & d & e \\ f & g & h\end{bmatrix}\begin{bmatrix}x \\ y \\ 1\end{bmatrix}$$
Expanding the linear equations we get $$wx' = ax + by + c$$ $$wy' =dx + ey + f$$ $$w = gx + hy + 1$$
Substituiting our equation for $w$ into our first and second equations to get: $$x' = ax+by+c-gxx'-hyx'$$ $$y' = ax+by+c-gxy'-hyy'$$
We can rewrite this into a matrix vector form: $$\begin{bmatrix}x & y & 1 & 0 & 0 & 0 & -xx' & -yx'\\ 0 & 0 & 0 & x & y & 1 & -xy' & -yy'\end{bmatrix} \begin{bmatrix}a \\ b \\ c \\ d \\ e \\ f \\ g \\ h\end{bmatrix} = \begin{bmatrix} x' \\ y' \end{bmatrix}$$
With more points, we see that this form turns into the general $$Ah=b$$ form that can be solved with least squares to recover $h$.
Image Rectification¶
Now that we can recover Homographies, we need to warp images. Before we start warping for stitching, we can "rectify" an image, as in make an image planar. This is done by defining 4 correspondences on the image, and fixing the corresponding points to be a rectangle.
I chose a picture I took in the Palace of Versailles that was too big to take from directly in front. The image below is the original image with the hand picked correspondences overlayed. The results of the rectification are shown below that.
versailles = plt.imread("data/versailles.png")[...,:-1]
pts_versailles = np.load("out/versailles.npy")
plt.figure(figsize=(10,10))
plt.imshow(versailles)
plt.axis("off")
plt.scatter(*zip(*pts_versailles), color="orange")
plt.show()
def computeH(im1_pts,im2_pts):
assert im1_pts.shape == im2_pts.shape
A = np.zeros((im1_pts.shape[0]*2, 8))
b = np.zeros((im1_pts.shape[0]*2, 1))
for i in range(im1_pts.shape[0]):
x, y = im1_pts[i]
x_p, y_p = im2_pts[i]
b[2*i] = x_p
b[2*i+1] = y_p
A[2*i] = np.array([x, y, 1, 0, 0, 0, -x*x_p, -y*x_p])
A[2*i+1] = np.array([0, 0, 0, x, y, 1, -x*y_p, -y*y_p])
h_vec = np.linalg.lstsq(A, b, rcond=None)[0]
h_vec = np.append(h_vec, 1)
return h_vec.reshape(3,3)
# # # Forward map - bad has holes
# from scipy.interpolate import interp2d
# def warpImage(im, H):
# y_max, x_max = im.shape[:2]
# x_range = range(x_max)
# y_range = range(y_max)
# xx, yy = np.meshgrid(x_range, y_range)
# base_coord = np.vstack((xx.flatten(), yy.flatten(), np.ones(xx.flatten().shape)))
# trans_coord_homo = H.dot(base_coord)
# trans_coord = (trans_coord_homo[:2]/trans_coord_homo[-1])
# bbox_input = np.array([[0, 0, im.shape[1], im.shape[1]],
# [0, im.shape[0], 0, im.shape[0]],
# [1, 1, 1, 1]])
# bbox_output_homo = H.dot(bbox_input)
# bbox_out = bbox_output_homo[:2]/bbox_output_homo[-1]
# shift = -np.minimum(np.min(bbox_out, axis=1), np.zeros((2)))[:,None]
# bbox_out = (bbox_out+shift).astype(int)
# trans_coord = (trans_coord+shift).astype(int)
# base_coord = base_coord[:2].astype(int)
# new_im = np.zeros(tuple(np.max(bbox_out+1, axis=1))[::-1]+(3,))
# values = im[base_coord[1,:], base_coord[0,:]]
# new_im[trans_coord[1,:], trans_coord[0,:]] = values
# # max_x, max_y = np.max(trans_coord, axis=1)
# # bw_im = np.mean(new_im, axis=-1)
# # print(base_coord.shape)
# # print(values.shape)
# # points = base_coord.T
# # values = np.mean(values, axis=-1)
# # grid_x, grid_y = np.mgrid[0:new_im.shape[0], 0:new_im.shape[1]]
# # grid_z0 = griddata(points, values, (grid_x, grid_y), method='linear')
# # interp_splines = [RectBivariateSpline(range(max_y), range(max_x), new_im[...,i]) for i in range(3)]
# # f = interp2d(range(max_y), range(max_x), bw_im.flatten())
# # xnew, ynew = np.meshgrid(range(new_im.shape[0]), range(new_im.shape[1]))
# # tcks = [interpolate.bisplrep(x, y, new_im[...,i], s=0) for i in range(3)]
# # znew = interpolate.bisplev(xnew[:,0], ynew[0,:], tck)
# # print(new_im.shape)
# # print(im.shape)
# return new_im
# Inverse warp - good
import cv2
def warpImage(im, H):
# y_max, x_max = im.shape[:2]
# x_range = range(x_max)
# y_range = range(y_max)
# xx, yy = np.meshgrid(x_range, y_range)
# base_coord = np.vstack((xx.flatten(), yy.flatten(), np.ones(xx.flatten().shape)))
# bbox_input = np.array([[0, 0, im.shape[1]-1, im.shape[1]-1],
# [0, im.shape[0]-1, 0, im.shape[0]-1],
# [1, 1, 1, 1]])
# bbox_output_homo = H.dot(bbox_input)
# bbox_out = bbox_output_homo[:2]/bbox_output_homo[-1]
# shift = (-np.minimum(np.min(bbox_out, axis=1), np.zeros((2)))[:,None]).astype(int)
# shifted_bbox_out = (bbox_out+shift).astype(int)
# max_x, max_y = np.max(shifted_bbox_out, axis=1)
# new_im = np.zeros((max_y+1, max_x+1, 3))
new_im = np.zeros_like(im)
# shift_val = shift[:,0]
# x_range = range(-shift_val[0], new_im.shape[1]-shift_val[0])
# y_range = range(-shift_val[1], new_im.shape[0]-shift_val[1])
# print(x_range, y_range)
# xx, yy = np.meshgrid(x_range, y_range)
xx, yy = np.meshgrid(range(new_im.shape[1]), range(new_im.shape[0]))
# print(xx.flatten().shape)
# print(bbox_output_homo[-1].shape)
base_coord = np.vstack((xx.flatten(), yy.flatten(), np.ones(xx.flatten().shape)))
trans_coord_homo = np.linalg.inv(H).dot(base_coord)
trans_coord = trans_coord_homo[:2]/trans_coord_homo[-1]
shift = -np.minimum(np.min(trans_coord, axis=1), np.zeros((2)))[:,None]
trans_coord = (trans_coord+shift).astype(int)
base_coord = (base_coord[:2]).astype(int)
bounds = np.max(trans_coord, axis=1)[::-1]
# print(bounds)
y_pad = max(0, bounds[0]+1-im.shape[0])
x_pad = max(0, bounds[1]+1-im.shape[1])
padding = ((y_pad//2, y_pad//2+y_pad%2), (x_pad//2, x_pad//2+x_pad%2), (0, 0))
im = np.pad(im, padding)
# im = cv2.resize(im, tuple(bounds[::-1]+1), interpolation = cv2.INTER_AREA)
# print(new_im.shape)
# print(im.shape)
values = im[trans_coord[1,:], trans_coord[0,:]]
new_im[base_coord[1,:], base_coord[0,:]] = values
return new_im
y_lim, x_lim = versailles.shape[:2]
square = np.array([[0, y_lim],
[0, 0],
[x_lim, 0],
[x_lim, y_lim]])
H = computeH(pts_versailles, square)
# H_inv = computeH(square, pts_versailles)
res = warpImage(versailles, H)
plt.figure(figsize=(10, 10))
plt.imshow(res)
plt.axis("off")
plt.show()
plt.imsave("versailles.png", res)
Mosaic Stitching¶
alley_left = plt.imread("data/alley_1.jpg")[::-1, ::-1]
alley_mid = plt.imread("data/alley_2.jpg")[::-1, ::-1]
alley_right = plt.imread("data/alley_3.jpg")[::-1, ::-1]
plt.figure(figsize=(20,20))
ims = [alley_left, alley_mid, alley_right]
for i, im in enumerate(ims, 1):
plt.subplot(1, 3, i)
plt.imshow(im)
plt.axis("off")
plt.show()
pts_1_2 = np.load("out/alley_1_2.npy")
pts_2_1 = np.load("out/alley_2_1.npy")
pts_2_3 = np.load("out/alley_2_3.npy")
pts_3_2 = np.load("out/alley_3_2.npy")
H1_2 = computeH(pts_1_2, pts_2_1)
H3_2 = computeH(pts_3_2, pts_2_3)
alley_1_warp = warpImage(alley_left, H1_2)
alley_3_warp = warpImage(alley_right, H3_2)
plt.figure(figsize=(20,10))
plt.imshow(np.hstack((alley_1_warp, alley_mid, alley_3_warp)))
plt.show()
