Harris角点检测
特征检测一般使用角点检测或边缘检测。对于特征比较明显的图像，角点检测基本足够。
特征一般是图像变化明显的地方，我们可以使用自相关来描述这一特点。但是角点检测更关注角而不是边缘，因为边缘在匹配时仍有一个自由度。
我们先对某点的自相关E(u,v)的计算泰勒展开取一阶近似，以便进一步讨论：
$$ E(u,v)=\sum_{\Delta x, \Delta y}{w(\Delta x,\Delta y)[I(u+\Delta x,y+\Delta v)-I(u,v)]^2}\\ \approx \sum_{\Delta x,\Delta y}{w(\Delta x,\Delta y)(I_x \Delta x+I_y \Delta y)^2} \\ = \begin{pmatrix} u \\ v \end{pmatrix} W \begin{pmatrix} I_x^2 & I_x I_y \\ I_x I_y & I_y^2 \end{pmatrix} \begin{pmatrix} u & v \end{pmatrix} $$
上式中w(x,y)为自相关计算的权重，一般取高斯分布。\(I_x, I_y\)是导数，乘上高斯分布后即Sobel算子。
考虑这个二次型的意义：它表示这个点的梯度情况。对它做特征分解，则特征值代表了两个主要方向的梯度值。所以只要两个特征值均比较大，就可以认为是角点。
这里只用特征值总体大小和它们的相对大小，使用角点响应公式简化计算：$$ R=det(M)-k * (trace(M))^2 = \lambda_1 \lambda_2 -k(\lambda_1 +\lambda_2)^2 $$
计算角点强度及其角点，返回响应强度和角点矩阵：
def computeHarrisValues(self, srcImage): ''' Input: srcImage -- Grayscale input image in a numpy array with values in [0, 1]. The dimensions are (rows, cols). Output: harrisImage -- numpy array containing the Harris score at each pixel. orientationImage -- numpy array containing the orientation of the gradient at each pixel in degrees. ''' height, width = srcImage.shape[:2] harrisImage = np.zeros(srcImage.shape[:2]) orientationImage = np.zeros(srcImage.shape[:2]) sobelx=ndimage.sobel(srcImage, 0) sobely=ndimage.sobel(srcImage, 1) x2=sobelx**2 y2=sobely**2 xy=sobelx*sobely x2=ndimage.gaussian_filter(x2, 0.5) y2=ndimage.gaussian_filter(y2, 0.5) xy=ndimage.gaussian_filter(xy, 0.5) harrisImage=(x2*y2-xy**2)-0.1*(x2+y2)**2 orientationImage=np.arctan2(sobelx,sobely)/np.pi*180 # Save the harris image as harris.png for the website assignment self.saveHarrisImage(harrisImage, srcImage) return harrisImage, orientationImage
该代码直接使用sobel算子计算梯度，若图像较模糊效果会不好。可配合高斯金字塔实现尺度不变性。
非极大值抑制 NMS
上述方式算出来的角点可能太多，而且容易在一些区域聚集，不利于之后的匹配。需要使用非极大值抑制抑制掉相邻的一些点。
return scipy.ndimage.filters.maximum_filter(harrisImage, 7)==harrisImage
上面是使用简单粗暴的7*7窗口排除非极大值点。
更优秀的方法是自适应非极大值抑制(ANMS)，通过不断缩小半径来找到更合适的点：
if USE_ANMS: radius=800 features=sorted(features,key=lambda f:f.response, reverse=True) fans=[] while True: for f in features: for f2 in fans: if math.hypot(f.pt[0]-f2.pt[0],f.pt[1]-f2.pt[1])<radius and f.response<f2.response*0.9: break else: fans.append(f) if (len(fans)>=400): return fans radius*=0.6 return features
特征描述符
找到特征点后需要描述这个点的特征，以便和另一个图中相应的点匹配。
最简单的方式是使用这个点周围的窗口，适合仅平移的情况：
class SimpleFeatureDescriptor(FeatureDescriptor): def describeFeatures(self, image, keypoints): ''' Input: image -- BGR image with values between [0, 255] keypoints -- the detected features, we have to compute the feature descriptors at the specified coordinates Output: desc -- K x 25 numpy array, where K is the number of keypoints ''' image = image.astype(np.float32) image /= 255. grayImage = cv2.cvtColor(image, cv2.COLOR_BGR2GRAY) desc = np.zeros((len(keypoints), 5 * 5)) padImage=np.pad(grayImage, ((2,2),(2,2)), 'constant') for i, f in enumerate(keypoints): x, y = f.pt x, y = int(x), int(y) desc[i]=padImage[y:y+5, x:x+5].reshape(25) return desc
MOPS描述符考虑了图的旋转和光度变化情况。同时可以考虑尺度变换，根据特征点的情况取适合大小的窗口，再缩放为8*8。处理时根据角度旋转，最后将图像归一化：
class MOPSFeatureDescriptor(FeatureDescriptor): # TODO: Implement parts of this function def describeFeatures(self, image, keypoints): ''' Input: image -- BGR image with values between [0, 255] keypoints -- the detected features, we have to compute the feature descriptors at the specified coordinates Output: desc -- K x W^2 numpy array, where K is the number of keypoints and W is the window size ''' image = image.astype(np.float32) image /= 255. # This image represents the window around the feature you need to # compute to store as the feature descriptor (row-major) windowSize = 8 desc = np.zeros((len(keypoints), windowSize * windowSize)) grayImage = cv2.cvtColor(image, cv2.COLOR_BGR2GRAY) grayImage = ndimage.gaussian_filter(grayImage, 0.5) for i, f in enumerate(keypoints): # compute the transform from each pixel in the 40x40 rotated window surrounding transMx = np.zeros((2, 3)) T2=np.array([[1,0,-f.pt[0]+4],[0,1,-f.pt[1]+4],[0,0,1]]) R=np.eye(3) R[:2,:]=cv2.getRotationMatrix2D(f.pt,f.angle,0.2) transMx=(T2 @ R)[:2,:] # Call the warp affine function to do the mapping # It expects a 2x3 matrix destImage = cv2.warpAffine(grayImage, transMx, (windowSize, windowSize), flags=cv2.INTER_LINEAR) if (math.fabs(np.std(destImage))>=1e-5): desc[i]=((destImage-np.mean(destImage))/np.std(destImage)).reshape(64) return desc
除MOPS之外还有一些效果较好的特征描述符，如梯度直方图等。梯度直方图根据之前计算梯度方向，在特征周围的窗口内进行统计。
特征匹配
有了两幅图中的特征，就可以进行匹配了。最简单的方法是计算两个特征的欧氏距离的平方(SSD)， 每个点直接取与它SSD最小的点。匹配的强度取它们的距离，只要距离小于阈值就认为是合适的匹配。
# get SSD match dis=scipy.spatial.distance.cdist(desc1, desc2)**2 for i in range(desc1.shape[0]): to=np.argmin(dis[i]) # DMathh: (from, to, _, threshold) matches.append(cv2.DMatch(i,to,0,dis[i,to]))
一种更优的做法是比率测试(Ratio Test)。如果说一个匹配可信度高，那么它应该与次小值相比近很多。
dis=scipy.spatial.distance.cdist(desc1, desc2)**2 for i in range(desc1.shape[0]): to=np.argmin(dis[i]) near1=dis[i,to] dis[i,to]=dis[i,0] near2=np.min(dis[i]) matches.append(cv2.DMatch(i,to,0,near1/near2))