# スタンフォード/CS131/宿題3-4 RANSACとHOG

スポンサーリンク

## Part 4 RANSAC¶

Rather than directly feeding all our keypoint matches into fit_affine_matrix function, we can instead use RANSAC (“RANdom SAmple Consensus”) to select only “inliers” to use to compute the transformation matrix.
fit_affine_matrix関数に全てのキーポイントマッチを直接フィードする代わりに、変換行列の算出に使うインライアだけを選択するのにRANSAC(ランダムサンプルコンセンサス)を利用できる。

The steps of RANSAC are(RANSACステップは):

1. Select random set of matches
ランダムなマッチした組を選択する。
2. Compute affine transformation matrix
アフィン変換行列を求める。
3. Find inliers using the given threshold

4. Repeat and keep the largest set of inliers

5. Re-compute least-squares estimate on all of the inliers

Implement ransac in panorama.py, run through the following code to get a panorama. You can see the difference from the result we get without RANSAC.
panorama.pyにransacを実装して、下のコードを走らせてパノラマ画像を得る。RANSACなし画像との違いを見ることができる。

def ransac(keypoints1, keypoints2, matches, n_iters=200, threshold=20):
"""
Use RANSAC to find a robust affine transformation
1. Select random set of matches
2. Compute affine transformation matrix
3. Compute inliers
4. Keep the largest set of inliers
5. Re-compute least-squares estimate on all of the inliers
Args:
keypoints1: M1 x 2 matrix, each row is a point
keypoints2: M2 x 2 matrix, each row is a point
matches: N x 2 matrix, each row represents a match
[index of keypoint1, index of keypoint 2]
n_iters: the number of iterations RANSAC will run
threshold: the number of threshold to find inliers
Returns:
H: a robust estimation of affine transformation from keypoints2 to
keypoints 1
"""
N = matches.shape[0]
n_samples = int(N * 0.2)
max_inliers = np.zeros(N)
n_inliers = 0
# RANSAC iteration start
idx = np.random.choice(N, n_samples, replace=False)
for i in range(n_iters):
H = np.linalg.lstsq(matched2[idx],matched1[idx], rcond=None)[0]
H[:, 2] = np.array([0, 0, 1])
c = np.linalg.norm(np.dot(matched2,H)-matched1,axis=-1)**2<threshold
n_inliers = np.sum(c)
max_inliers = idx

return H, matches[max_inliers]

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
plt.rcParams['figure.figsize'] = (15.0, 12.0) # set default size of plots
plt.rcParams["font.size"] = "17"
plt.rcParams['image.interpolation'] = 'nearest'
plt.rcParams['image.cmap'] = 'gray'

from utils import get_output_space,warp_image,pad,plot_matches
from panorama import harris_corners,describe_keypoints
from skimage.feature import corner_peaks
from panorama import simple_descriptor,match_descriptors,fit_affine_matrix

# Detect keypoints in two images
keypoints1 = corner_peaks(harris_corners(img1, window_size=3),
threshold_rel=0.05,
exclude_border=8)
keypoints2 = corner_peaks(harris_corners(img2, window_size=3),
threshold_rel=0.05,
exclude_border=8)

patch_size = 5
# Extract features from the corners
desc1 = describe_keypoints(img1, keypoints1,
desc_func=simple_descriptor,
patch_size=patch_size)
desc2 = describe_keypoints(img2, keypoints2,
desc_func=simple_descriptor,
patch_size=patch_size)

# Match descriptors in image1 to those in image2
matches = match_descriptors(desc1, desc2, 0.7)

# Extract matched keypoints
p1 = keypoints1[matches[:,0]]
p2 = keypoints2[matches[:,1]]

# Find affine transformation matrix H that maps p2 to p1
H = fit_affine_matrix(p1, p2)

H, robust_matches = ransac(keypoints1, keypoints2, matches, threshold=1)

# Visualize robust matches
fig, ax = plt.subplots(1, 1, figsize=(18, 14))
plot_matches(ax, img1, img2, keypoints1, keypoints2, robust_matches)
plt.axis('off')
plt.show()

plt.axis('off')
plt.title('RANSAC Solution')
plt.show()

output_shape, offset = get_output_space(img1, [img2], [H])

# Warp images into output sapce
img1_warped = warp_image(img1, np.eye(3), output_shape, offset)
img1_mask = (img1_warped != -1) # Mask == 1 inside the image
img1_warped[~img1_mask] = 0     # Return background values to 0

img2_warped = warp_image(img2, H, output_shape, offset)
img2_mask = (img2_warped != -1) # Mask == 1 inside the image
img2_warped[~img2_mask] = 0     # Return background values to 0

# Plot warped images
plt.subplot(1,2,1)
plt.imshow(img1_warped)
plt.title('Image 1 warped')
plt.axis('off')

plt.subplot(1,2,2)
plt.imshow(img2_warped)
plt.title('Image 2 warped')
plt.axis('off')
plt.show()

merged = img1_warped + img2_warped

overlap = (img1_mask * 1.0 +  # Multiply by 1.0 for bool -> float conversion

# Normalize through division by overlap - but ensure the minimum is 1
normalized = merged / np.maximum(overlap, 1)
plt.imshow(normalized)
plt.axis('off')
plt.show()

plt.axis('off')
plt.title('RANSAC Panorama Solution')
plt.show()

スポンサーリンク

## Part 5 Histogram of Oriented Gradients (HOG)¶

In the above code, you are using the simple_descriptor, and in this section, you are going to implement a simplified version of HOG descriptor.

HOG stands for Histogram of Oriented Gradients. In HOG descriptor, the distribution ( histograms ) of directions of gradients ( oriented gradients ) are used as features. Gradients ( x and y derivatives ) of an image are useful because the magnitude of gradients is large around edges and corners ( regions of abrupt intensity changes ) and we know that edges and corners pack in a lot more information about object shape than flat regions.
HOGは方向勾配ヒストグラムを意味する。HOG記述子では、勾配の方向(方向勾配)の分布(ヒストグラム)は特徴として使われる。画像の勾配(x・y微分係数)は、勾配の大きさがエッジとコーナー周辺で大きく(急激に強度が変わる)、エッジとコーナーが平坦領域よりも物体の形状に関する情報を多く詰め込んでいるので有益である。

The steps of HOG are(以下がHOGのステップ):

1. compute the gradient image in x and y directions
Use the sobel filter provided by skimage.filters
xとy方向の勾配画像を求める。
skimage.filtersのsobel filterを利用する。
Divide image into cells, and calculate histogram of gradients in each cell.

3. Flatten block of histograms into feature vector
ヒストグラムのブロックを特徴ベクトルに平坦化する。
4. Normalize flattened block

Implement hog_descriptor in panorama.py, and run through the following code to get a panorama image.
panorama.pyにhog_descriptorを実装し、以下のコードを全て実行してパノラマ画像を獲得する。

def hog_descriptor(patch, pixels_per_cell=(8,8)):
"""
Generating hog descriptor by the following steps:
1. compute the gradient image in x and y (already done for you)
3. normalize across block
4. flattening block into a feature vector
Args:
patch: grayscale image patch of shape (h, w)
pixels_per_cell: size of a cell with shape (m, n)
Returns:
block: 1D array of shape ((h*w*n_bins)/(m*n))
"""
assert (patch.shape[0] % pixels_per_cell[0] == 0),\
'Heights of patch and cell do not match'
assert (patch.shape[1] % pixels_per_cell[1] == 0),\
'Widths of patch and cell do not match'
n_bins = 9
degrees_per_bin = 180 // n_bins
Gx = filters.sobel_v(patch)
Gy = filters.sobel_h(patch)
G = np.sqrt(Gx**2 + Gy**2)
theta = (np.arctan2(Gy, Gx) * 180 / np.pi) % 180
G_cells = view_as_blocks(G, block_shape=pixels_per_cell)
theta_cells = view_as_blocks(theta, block_shape=pixels_per_cell)
rows = G_cells.shape[0]
cols = G_cells.shape[1]
cells = np.zeros((rows, cols, n_bins))
# Compute histogram per cell
for i in range(rows):
for j in range(cols):
hist,_ = np.histogram(theta_cells[i,j], \
bins=n_bins,range=(0,180),weights=G_cells[i,j])
cells[i, j] = hist
cells = (cells - np.mean(cells))/np.std(cells)
block = cells.ravel()
return block

img1 = imread('uttower1.jpg', as_gray=True)

# Detect keypoints in both images
keypoints1 = corner_peaks(harris_corners(img1, window_size=3),
threshold_rel=0.05,
exclude_border=8)
keypoints2 = corner_peaks(harris_corners(img2, window_size=3),
threshold_rel=0.05,
exclude_border=8)

from skimage import filters
from skimage.util.shape import view_as_blocks
# Extract features from the corners

desc1 = describe_keypoints(img1, keypoints1,
desc_func=hog_descriptor,
patch_size=16)
desc2 = describe_keypoints(img2, keypoints2,
desc_func=hog_descriptor,
patch_size=16)
# Match descriptors in image1 to those in image2
matches = match_descriptors(desc1, desc2, 0.7)

# Plot matches
fig, ax = plt.subplots(1, 1, figsize=(15, 12))
ax.axis('off')
plot_matches(ax, img1, img2, keypoints1, keypoints2, matches)
plt.show()
plt.axis('off')
plt.title('HOG descriptor Solution')
plt.show()

H, robust_matches = ransac(keypoints1, keypoints2, matches, threshold=1)

# Plot matches
fig, ax = plt.subplots(1, 1, figsize=(15, 12))
plot_matches(ax, img1, img2, keypoints1, keypoints2, robust_matches)
plt.axis('off')
plt.show()

plt.axis('off')
plt.title('HOG descriptor Solution')
plt.show()

output_shape, offset = get_output_space(img1, [img2], [H])
print(output_shape, offset)

# Warp images into output sapce
img1_warped = warp_image(img1, np.eye(3), output_shape, offset)
img1_mask = (img1_warped != -1) # Mask == 1 inside the image
img1_warped[~img1_mask] = 0     # Return background values to 0

img2_warped = warp_image(img2, H, output_shape, offset)
img2_mask = (img2_warped != -1) # Mask == 1 inside the image
img2_warped[~img2_mask] = 0     # Return background values to 0

# Plot warped images
plt.subplot(1,2,1)
plt.imshow(img1_warped)
plt.title('Image 1 warped')
plt.axis('off')

plt.subplot(1,2,2)
plt.imshow(img2_warped)
plt.title('Image 2 warped')
plt.axis('off')
plt.show()

[438 913] [0. 0.]

merged = img1_warped + img2_warped

overlap = (img1_mask * 1.0 +  # Multiply by 1.0 for bool -> float conversion

# Normalize through division by overlap - but ensure the minimum is 1
normalized = merged / np.maximum(overlap, 1)
plt.imshow(normalized)
plt.axis('off')
plt.show()

plt.axis('off')
plt.title('HOG Descriptor Panorama Solution')
plt.show()

スポンサーリンク

## Extra Credit: Better Image Merging¶

You will notice the blurry region and unpleasant lines in the middle of the final panoramic image. In the cell below, come up with a better merging scheme to make the panorama look more natural. Be creative!

# Modify the code below
# np.set_printoptions(threshold='nan')
merged = img1_warped+img2_warped
output = merged/np.maximum(overlap,1)
plt.imshow(output)
plt.axis('off')
plt.title('method1')
plt.show()

merged = img1_warped + img2_warped
for i in range(img2_warped.shape[1]):
break

RL = i
RR = img1.shape[1]
RW = RR-RL+1

for i in range(RL,RR+1):
for j in range(img2_warped.shape[0]):
a = (i-RL)/(RW)
a = 1-a
merged[j,i]=a*img1_warped[j,i] \
+ (1-a)*img2_warped[j,i]

output = merged