Add example matching script

doc/examples/plot_matching.py

@@ -0,0 +1,143 @@
+"""
+============================
+Robust matching using RANSAC
+============================
+
+In this simplified example we first generate two synthetic images as if they
+were  taken from different view points.
+
+In the next step we find interest points in both images and find
+correspondencies based on a weighted sum of squared differences measure. Note,
+that this measure is only robust towards linear radiometric and not geometric
+distortions and is thus only usable with slight view point changes.
+
+After finding the correspondencies we end up having a set of source and
+destination coordinates which can be used to estimate the geometric
+transformation between both images. However, many of the correspondencies are
+faulty and simply estimating the parameter set with all coordinates is not
+sufficient. Therefore, the RANSAC algorithm is used on top of the normal model
+to robustly estimate the parameter set by detecting outliers.
+
+"""
+import numpy as np
+from matplotlib import pyplot as plt
+
+from skimage import data
+from skimage.feature import corner_harris, corner_subpix, corner_peaks
+from skimage.transform import warp, AffineTransform
+from skimage.exposure import rescale_intensity
+from skimage.color import rgb2gray
+from skimage.measure import ransac
+
+
+# generate synthetic checkerboard image and add gradient for the later matching
+checkerboard = data.checkerboard()
+img_orig = np.zeros(list(checkerboard.shape) + [3])
+img_orig[..., 0] = checkerboard
+gradient_r, gradient_c = np.mgrid[0:img_orig.shape[0], 0:img_orig.shape[1]]
+img_orig[..., 1] = gradient_r
+img_orig[..., 2] = gradient_c
+img_orig = rescale_intensity(img_orig)
+img_orig_gray = rgb2gray(img_orig)
+
+# warp synthetic image
+tform = AffineTransform(scale=(0.9, 0.9), rotation=0.2, shear=0,
+                        translation=(20, -10))
+img_warped = warp(img_orig, tform.inverse, output_shape=(200, 200))
+img_warped_gray = rgb2gray(img_warped)
+
+# extract corners using Harris' corner measure
+coords_orig = corner_peaks(corner_harris(img_orig_gray), threshold_rel=0.001,
+                           min_distance=5)
+coords_warped = corner_peaks(corner_harris(img_warped_gray),
+                             threshold_rel=0.001, min_distance=5)
+
+# determine subpixel corner position
+coords_orig_subpix = corner_subpix(img_orig_gray, coords_orig, window_size=10)
+coords_warped_subpix = corner_subpix(img_warped_gray, coords_warped,
+                                     window_size=10)
+
+
+def gaussian_weights(window_ext, sigma=1):
+    y, x = np.mgrid[-window_ext:window_ext+1, -window_ext:window_ext+1]
+    g = np.zeros(y.shape, dtype=np.double)
+    g[:] = np.exp(-0.5 * (x**2 / sigma**2 + y**2 / sigma**2))
+    g /= 2 * np.pi * sigma * sigma
+    return g
+
+
+def match_corner(coord, window_ext=5):
+    r, c =  np.round(coord)
+    window_orig = img_orig[r-window_ext:r+window_ext+1,
+                           c-window_ext:c+window_ext+1, :]
+
+    # weight pixels depending on distance to center pixel
+    weights = gaussian_weights(window_ext, 3)
+    weights = np.dstack((weights, weights, weights))
+
+    # compute sum of squared differences to all corners in warped image
+    SSDs = []
+    for cr, cc in coords_warped:
+        window_warped = img_warped[cr-window_ext:cr+window_ext+1,
+                                   cc-window_ext:cc+window_ext+1, :]
+        SSD = np.sum(weights * (window_orig - window_warped)**2)
+        SSDs.append(SSD)
+
+    # use corner with minimum SSD as correspondency
+    min_idx = np.argmin(SSDs)
+    return coords_warped_subpix[min_idx]
+
+
+# find correspondencies using simple weighted sum of squared differences
+src = []
+dst = []
+for coord in coords_orig_subpix:
+    src.append(coord)
+    dst.append(match_corner(coord))
+src = np.array(src)
+dst = np.array(dst)
+
+
+# estimate affine transform model using all coordinates
+model = AffineTransform()
+model.estimate(src, dst)
+
+# robustly estimate affine transform model with RANSAC
+model_robust, inliers = ransac((src, dst), AffineTransform, min_samples=3,
+                               residual_threshold=2, max_trials=100)
+
+
+# compare "true" and estimated transform parameters
+print tform.scale, tform.translation, tform.rotation
+print model.scale, model.translation, model.rotation
+print model_robust.scale, model_robust.translation, model_robust.rotation
+
+
+# visualize correspondencies
+img_combined = np.concatenate((img_orig_gray, img_warped_gray), axis=1)
+
+fig, ax = plt.subplots(nrows=2, ncols=1)
+plt.gray()
+
+ax[0].imshow(img_combined, interpolation='nearest')
+ax[0].axis('off')
+ax[0].axis((0, 400, 200, 0))
+ax[0].set_title('Correct correspondencies')
+ax[1].imshow(img_combined, interpolation='nearest')
+ax[1].axis('off')
+ax[1].axis((0, 400, 200, 0))
+ax[1].set_title('Faulty correspondencies')
+
+for i in range(len(src)):
+    if i in inliers:
+        ax_idx = 0
+        color = 'g'
+    else:
+        ax_idx = 1
+        color = 'r'
+    ax[ax_idx].plot((src[i, 1], dst[i, 1] + 200), (src[i, 0], dst[i, 0]), '-',
+                    color=color)
+    ax[ax_idx].plot(src[i, 1], src[i, 0], '.', markersize=10, color=color)
+    ax[ax_idx].plot(dst[i, 1] + 200, dst[i, 0], '.', markersize=10, color=color)
+
+plt.show()