1、实验内容
实验题目：Implement Panorama Stitching with Harris corner detector, RANSAC and HOG descriptor.
实验步骤：
- 1.使用Harris焦点检测器寻找关键点。
- 2.构建描述算子来描述图中的每个关键点，比较两幅图像的两组描述子，并进行匹配。
- 3.根据一组匹配关键点，使用最小二乘法进行仿射变换矩阵的计算。
- 4.使用RANSAC计算一个更加稳定的仿射变换的矩阵，然后将第二幅图变换过来并覆盖在第一幅图上，形成一个全景。
- 5.实现不同的描述子，并得到不同的拼接结果。
2、实验设计（略）
3、实验环境及实验数据集
实验环境:
- Windows10的操作系统、anaconda3、python3.6、jupyter notebook。
数据集:
- 任务给出的原图片以及处理结果后的图片，如image1.jpg、image2.jpg、solution_harris.png等等。
四、实验过程及结果
4.1 Harris角点检测器寻找关键点
结果展示：
实现代码：
def harris_corners(img, window_size=3, k=0.04): """ Compute Harris corner response map. Follow the math equation R=Det(M)-k(Trace(M)^2). Hint: You may use the function scipy.ndimage.filters.convolve, which is already imported above. Args: img: Grayscale image of shape (H, W) window_size: size of the window function k: sensitivity parameter Returns: response: Harris response image of shape (H, W) """ H, W = img.shape window = np.ones((window_size, window_size)) response = np.zeros((H, W)) dx = filters.sobel_v(img) . dy = filters.sobel_h(img) ### YOUR CODE HERE Ix2 = np.multiply(dx,dx) Iy2 = np.multiply(dy,dy) #高斯加权 A = convolve(Ix2,window) B = convolve(Iy2,window) C = convolve(IxIy,window) detM = np.multiply(A, B) - np.multiply(C, C) .traceM = A + B response = detM - k * np.square(traceM) ### END YOUR CODE return response
4.2 构建描述算子来描述图中的每个关键点，比较两幅图像的两组描述子，并进行匹配。
结果展示：
实现代码：
构建描述算子
def simple_descriptor(patch): """ Describe the patch by normalizing the image values into a standard normal distribution (having mean of 0 and standard deviation of 1) and then flattening into a 1D array. The normalization will make the descriptor more robust to change in lighting condition. Hint: If a denominator is zero, divide by 1 instead. Args: patch: grayscale image patch of shape (H, W) Returns: feature: 1D array of shape (H * W) """ feature = [] ### YOUR CODE HERE mean = np.mean(patch) sigma = np.std(patch) . if sigma == 0.0: sigma = 1 #均值0，标准差为1 normalized = (patch - mean) / sigma #生成一维特征向量 feature = normalized.flatten() ### END YOUR CODE return feature
根据描述子进行匹配
#注意，这里引入了heapq堆排列模块算法 import heapq def match_descriptors(desc1, desc2, threshold=0.5): """ Match the feature descriptors by finding distances between them. A match is formed when the distance to the closest vector is much smaller than the distance to the second-closest, that is, the ratio of the distances should be smaller than the threshold. Return the matches as pairs of vector indices. Hint: The Numpy functions np.sort, np.argmin, np.asarray might be useful Args: desc1: an array of shape (M, P) holding descriptors of size P about M keypoints desc2: an array of shape (N, P) holding descriptors of size P about N keypoints Returns: matches: an array of shape (Q, 2) where each row holds the indices of one pair of matching descriptors """ matches = [] N = desc1.shape[0] dists = cdist(desc1, desc2) ### YOUR CODE HERE for i in range(N): #找到最近的两个值 min2 = heapq.nsmallest(2, dists[i,:]) #保证最小的和第二小的有一定的距离 if min2[0] / min2[1] < threshold: matches.append([i, np.argmin(dists[i,:])]) matches = np.asarray(matches) ### END YOUR CODE return matches
4.3 根据一组匹配关键点，使用最小二乘法进行仿射变换矩阵的计算
实验结果：
实验代码：
def fit_affine_matrix(p1, p2): """ Fit affine matrix such that p2 * H = p1 Hint: You can use np.linalg.lstsq function to solve the problem. Args: p1: an array of shape (M, P) p2: an array of shape (M, P) Return: H: a matrix of shape (P, P) that transform p2 to p1. """ assert (p1.shape[0] == p2.shape[0]),\ 'Different number of points in p1 and p2' p1 = pad(p1) p2 = pad(p2) ### YOUR CODE HERE #仿射变换矩阵 H, residuals, rank, s = np.linalg.lstsq(p2, p1, rcond=None) ### END YOUR CODE # Sometimes numerical issues cause least-squares to produce the last # column which is not exactly [0, 0, 1] H[:,2] = np.array([0, 0, 1]) return H
4.4 使用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 """ # Copy matches array, to avoid overwriting it orig_matches = matches.copy() matches = matches.copy() N = matches.shape[0] print(N) n_samples = int(N * 0.2) matched1 = pad(keypoints1[matches[:,0]]) matched2 = pad(keypoints2[matches[:,1]]) max_inliers = np.zeros(N) n_inliers = 0 # RANSAC iteration start ### YOUR CODE HERE for i in range(n_iters): inliersArr = np.zeros(N) idx = np.random.choice(N, n_samples, replace=False) p1 = matched1[idx, :] p2 = matched2[idx, :] H, residuals, rank, s = np.linalg.lstsq(p2, p1, rcond=None) H[:, 2] = np.array([0, 0, 1]) output = np.dot(matched2, H) inliersArr = np.linalg.norm(output-matched1, axis=1)**2 < threshold inliersCount = np.sum(inliersArr) if inliersCount > n_inliers: . max_inliers = inliersArr.copy() #这里需要注意深拷贝和浅拷贝 n_inliers = inliersCount # 迭代完成，拿最大数目的匹配点对进行估计变换矩阵 H, residuals, rank, s = np.linalg.lstsq(matched2[max_inliers], matched1[max_inliers], rcond=None) H[:, 2] = np.array([0, 0, 1]) ### END YOUR CODE print(H) return H, orig_matches[max_inliers]
4.5 实现不同的描述子HOG，得到不同的拼接结果
实验结果：
实验代码：
def hog_descriptor(patch, pixels_per_cell=(8,8)): """ Generating hog descriptor by the following steps: Compute the gradient image in x and y directions (already done for you) Compute gradient histograms for each cell Flatten block of histograms into a 1D feature vector Here, we treat the entire patch of histograms as our block 4Normalize flattened block Normalization makes the descriptor more robust to lighting variations Args: patch: grayscale image patch of shape (H, W) pixels_per_cell: size of a cell with shape (M, N) Returns: block: 1D patch descriptor 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) # Unsigned gradients G = np.sqrt(Gx**2 + Gy**2) theta = (np.arctan2(Gy, Gx) * 180 / np.pi) % 180 # Group entries of G and theta into cells of shape pixels_per_cell, (M, N) # G_cells.shape = theta_cells.shape = (H//M, W//N) # G_cells[0, 0].shape = theta_cells[0, 0].shape = (M, N) 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] # For each cell, keep track of gradient histrogram of size n_bins . cells = np.zeros((rows, cols, n_bins)) # Compute histogram per cell ### YOUR CODE HERE # 首先需要遍历Cell for row in range(rows): for col in range(cols): # 遍历cell中的像素 for y in range(pixels_per_cell[0]): for x in range(pixels_per_cell[1]): # 计算该像素的梯度方向（在n_bins个方向中属于第几个区间） angle = theta_cells[row,col,y,x] order = int(angle) // degrees_per_bin if order == 9: order = 8 # 统计该cell中每个方向区间内的强度 # 累加强度也可以 cells[row,col,order] += G_cells[row,col,y,x] # 最后做归一化处理，直方图归一化 cells = (cells - cells.mean()) / (cells.std()) block = cells.reshape(-1) # 一维 . ### YOUR CODE HERE return block
完整实验见下载连接：
