Let's prepare our environment and load image pair
import cv2
import scipy.ndimage as nd
figsize(12, 10)
img1 = cv2.imread('data/temp.jpg')
img2 = cv2.imread('data/DSC_0059.JPG')
img1 = cv2.cvtColor(img1, cv2.COLOR_BGR2GRAY)
img2 = cv2.cvtColor(img2, cv2.COLOR_BGR2GRAY)
print img1.shape, img2.shape
(400L, 600L) (4000L, 6016L)
img2
is hige! Downscale it:
for _ in xrange(3):
img2 = cv2.pyrDown(img2)
print img2.shape
(500L, 752L)
gray()
subplot(121)
imshow(img1)
subplot(122)
imshow(img2)
<matplotlib.image.AxesImage at 0x7f6a898>
This image pair is hard to match with conventinal keypoint-based alignment. Let's try to take a look on their DFTs. Don't forget to apply windowing to suppres boundary artifacts.
def win_fft(a):
h, w = a.shape
wnd = hamming
a = a - a.mean()
aw = a * wnd(w)*wnd(h)[:,newaxis]
return fftshift(fft2(aw))
F1 = abs(win_fft(img1))
F2 = abs(win_fft(img2))
imshow(log(F1), cmap='jet')
<matplotlib.image.AxesImage at 0x7fafa90>
imshow(log(F2), cmap='jet')
<matplotlib.image.AxesImage at 0x1056d240>
As we see, spectrums contain aligned series of peaks from harmonics, produced by repetitive field texture. Also they have lines produced by straight roads. Aligning theses spectrums (for example by matching peaks or histograms) can give us rotation and scale. This makes searching for translation vector much easier.
Let's make a few plots to illustrate the idea.
def make_xy(F):
h, w = F.shape
y, x = indices((h, w))
return x-w//2, y-h//2
def ang_hist(F):
x, y = make_xy(F)
ang = arctan2(y, x) % pi
hist(ang.ravel(), 200, weights=F.ravel())
xlabel('Angle')
ylabel('Power')
def r_hist(F):
x, y = make_xy(F)
r = sqrt(x*x + y*y)
hist(r.ravel(), 200, weights=F.ravel())
xlabel('Frequency')
ylabel('Power')
subplot(211)
title('img1')
ang_hist(F1**2)
subplot(212)
title('img2')
ang_hist(F2**2)
The angular energy distribution can be used to match rotaion.
ax = subplot(211)
title('img1')
r_hist(F1**2)
subplot(212, sharex=ax)
title('img2')
r_hist(F2**2)
h, w = F2.shape
F2p = cv2.logPolar(F2, (w//2, h//2), 100.0, cv2.INTER_AREA)
imshow(log(F2p), cmap='jet')
<matplotlib.image.AxesImage at 0x36db4cc0>
h, w = F1.shape
F1p = cv2.logPolar(F1, (w//2+1, h//2+1), 100.0, cv2.INTER_AREA)
imshow(log(F1p), cmap='jet')
<matplotlib.image.AxesImage at 0x3b532dd8>
a1 = cv2.resize(F1p, (512, 512))
a2 = cv2.resize(F2p, (512, 512))
cv2.phaseCorrelate(a1, a2)
((-90.2063829500409, -6.2185609284478005), -0.05441922318895057)
imshow(a1-a2)
<matplotlib.image.AxesImage at 0x371ddb70>
v, W = cv2.findTransformECC(img2, img1, np.float32([[1, 0, 0],[0, 1, 0]]), cv2.MOTION_AFFINE)
print v
print W
0.14068992121 [[ 0.93549651 -0.03079971 12.88512707] [ 0.08168394 1.03515792 -11.24993896]]
h, w = img2.shape
img1t = cv2.warpAffine(img1, W, (w, h), flags=cv2.WARP_INVERSE_MAP)
imshow(img1*0.5+0.5*img2t)
<matplotlib.image.AxesImage at 0xb2145f8>