视觉实验二特征检测与匹配

可以看到这次实验大致上可以分为三部分：
特征检测、特征描述、特征匹配

一、特征检测

角点检测需要我们实现Harris角点检测功能
计算Harris的2*2矩阵,如图所示

使用3*3的Sobel算子对图上每个点计算x方向与y方向导数，也就是

imx = np.zeros(srcImage.shape)
ndimage.sobel(srcImage, 1, imx)
imy = np.zeros(srcImage.shape)
ndimage.sobel(srcImage, 0, imy)

超出边界默认为reflect模式。有了Ixp与Iyp，这之后直接使用高斯滤波，sigma设为0.5.

1
2
3

wxx = ndimage.gaussian_filter(imx * imx, sigma=0.5)
wxy = ndimage.gaussian_filter(imx * imy, sigma=0.5)
wyy = ndimage.gaussian_filter(imy * imy, sigma=0.5)

这样就有了每个点的Harris矩阵值，根据要求计算角点强度与方向.

wdet = wxx * wyy - wxy ** 2
wtr = wxx + wyy
harrisImage = wdet - 0.1 * (wtr ** 2)
orientationImage = np.arctan2(imy, imx)* 180.0 / np.pi

注意这里的角度单位，np.arctan2计算出的是弧度制，转化为角度制。求出每个像素的角点强度与方向后，需要找到在每一个上下左右7*7矩阵中角点强度最大的点，这些点就是我们需要的特征点.（TODO2,3）
用到了最大过滤器maximum_filter

destImage = np.zeros_like(harrisImage, np.bool)
height = harrisImage.shape[0]
width = harrisImage.shape[1]
maxImage = ndimage.maximum_filter(harrisImage,7)
for h in range(height):
    for w in range(width):
        if(harrisImage[h][w]==maxImage[h][w]):
            destImage[h][w]=True
        else:
            destImage[h][w]=False
return destImage

使用cv2.KeyPoint()存储检测到的每一个特征点.

二、特征描述

实现两种特征描述符，一种是只考虑平移的简单描述符，一种是MOPS（方形，考虑平移旋转）
1.简单描述符，存储5*5邻域的强度值

image = image.astype(np.float32)
image /= 255.
grayImage = cv2.cvtColor(image, cv2.COLOR_BGR2GRAY)
height, width = image.shape[:2]
desc = np.zeros((len(keypoints), 5 * 5))
for i, f in enumerate(keypoints):
    x, y = f.pt
    x, y = int(x), int(y)
    num = 0;
    for j in range(-2,3):
        for k in range(-2,3):
            if(y+j<0 or y+j>height-1 or x+k<0 or x+k>width-1):
                desc[i][num] = 0
            else:
                desc[i][num] = grayImage[y+j][x+k]
            num+=1
return desc

2.MOPS描述符，用8*8矩阵存储角度强度朝向右方的40*40邻域

分别求出平移、旋转、缩放、平移的仿射变换矩阵，采用齐次坐标
第一步平移：让位于图像左上角的坐标(0,0)点移动到计算的中心像素上
第二步旋转：将图像绕中心像素旋转，让该点的角点方向朝右
第三步缩放：将图像缩小为1/8
第四步平移：将原点坐标移动到中心像素8*8邻域的左上角像素
这样就可以读出该特征点的MOPS描述符（8*8矩阵）

注意，平移时是对坐标值进行操作，所以让(-4,-4)变为(0,0)是要加4，旋转角度用弧度制。将上述四个仿射矩阵倒序相乘，用cv2.warpAffine实现变换，将结果进行检验，也就是实验说明提到的将 8*8 图像块规范化为零均值和单位方差。

这句话很让人纠结，零均值和单位方差完全是字面意思上的翻译。查了后发现这在数据处理中叫做中心化（又叫零均值化）和标准化（又叫归一化）。前者是指，让每个数据减去平均值；后者则将中心化后的数据除以标准差。这样做的目的是得到均值为0，标准差为1的标准正态分布。

for i, f in enumerate(keypoints):
    x, y = f.pt
    x, y = int(x), int(y)
    transMx = np.zeros((2, 3))

    #平移到目标像素
    t1 = np.array([-x,-y,1])
    get_t1 = np.eye(3)
    get_t1[:3, 2] = t1

    #旋转
    rot = -f.angle/180*np.pi
    get_rot = np.array([[math.cos(rot), -math.sin(rot), 0],
                         [math.sin(rot), math.cos(rot), 0],
                         [0, 0, 1]])
    #缩放(缩小五倍)
    get_scale = np.array([[0.2, 0, 0],
                         [0, 0.2, 0],
                         [0, 0, 1]])

    #平移
    t2 = np.array([4,4,1])
    get_t2 = np.eye(3)
    get_t2[:3, 2] = t2

    a = np.matmul(get_t2,get_scale)
    b = np.matmul(a,get_rot)
    c = np.matmul(b,get_t1)
    transMx = c[: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)
    mean = np.mean(destImage)
    variance = np.std(destImage)
    newdest = np.zeros(destImage.shape)
    if variance<0.00001:
        pass
    else:
        newdest = (destImage-mean)/variance
    desc[i]=newdest.flatten()
return desc

三、特征匹配

下一步是编写匹配它们的代码（即，在一个图像中给定一个特征，在另一个图像中找到最佳匹配特征）。最简单的方法如下：比较两个特征并计算它们之间的标量距离，最佳匹配是距离最小的特征。你将实现两个距离函数：
平方差之和（SSD）：这是两个特征向量之间的欧氏距离平方。
比率测试：通过 SSD距离查找最近和次近的两个特征，比率测试距离是它们的比率（即，最接近的特征匹配的SSD距离除以第二接近的特征匹配的 SSD 距离）。
求欧氏距离的平方，使用scipy.spatial.distance.cdist()，再用np.argmin找到每个特征点的最小距离（也就是最佳匹配）

de = spatial.distance.cdist(desc1, desc2, metric='euclidean')
for n, dee in enumerate(de):
    d = cv2.DMatch()
    d.queryIdx = n
    m = np.argmin(de[n])
    d.trainIdx = int(m)
    d.distance = de[n][m]
    matches.append(d)
return matches

如果是比例测试，则用最小除以第二小作为距离

de = spatial.distance.cdist(desc1, desc2, metric='euclidean')
for n, dee in enumerate(de):
    d = cv2.DMatch()
    d.queryIdx = n
    m = np.argmin(de[n])
    d.trainIdx = int(m)
    f1 = de[n][m]
    de[n][m] = float("inf")
    m = np.argmin(de[n])
    f2 = de[n][m]
    d.distance = f1/f2
    matches.append(d)
return matches