Python机器视觉（2）——特征检测

角点检测

       角点检测可以提取点的邻域特征，一般是根据梯度大小和方向判断。Harris自相关检测提供了一种角点检测思路。

       它希望用一个窗函数（可以是中值滤波算子，也可以是高斯算子）来扫描整幅图像，并加一定的偏移量（u,v）计算邻域像素点的差值（平方是避免差值为负数），最后求窗函数所包含的像素点的加权值，然后移动窗函数，继续计算。只要找到了一个极小值点，那么它应该为角点。

       但是，这样扫描会十分耗时，Harris用数学的方法，将上方公式进行泰勒二阶展开：

       整理可得：

       M包含了像素点在x和y方向的导数，并与窗函数求加权值。由于M实对称矩阵的特性，可以对其求特征值，来表达平坦区域、边缘区域和角点：

       为了更好表达，采用角点响应函数

       R为大数值负数时，是边缘区；R为小数值时，是平坦区；R为大数值正数时，是角点。

       由于R的计算会出现很多不合理的点，需要进行非极大值抑制。

       由于R的值会十分大，需要进行拉伸变换到[0,255]区间。

       效果如下：

       代码如下，后期会补一下注释，滤波器已经在Python机器视觉（1）中解释过。

import cv2
import numpy as np
import math
import matplotlib.pyplot as plt

def Filter(img,direction):
    core_size = 3
    # 零填充
    pad = core_size//2
    if direction == 'gaussian':
        filter = np.zeros((img.shape[0]+2*pad,img.shape[1]+2*pad,3),dtype=np.uint8)
        filter[pad:pad+img.shape[0],pad:pad+img.shape[1],:] = img.astype(np.uint8)
    else:
        filter = np.zeros((img.shape[0]+2*pad,img.shape[1]+2*pad),dtype=np.uint8)
        filter[pad:pad+img.shape[0],pad:pad+img.shape[1]] = img.astype(np.uint8)
    #滤波计算
    tmp = filter.copy()
    if direction == 'x':
        core = np.array([[1,0,-1],[2,0,-2],[1,0,-1]])/9.0
    elif direction == 'y':
        core = np.array([[1,2,1],[0,0,0],[-1,-2,-1]])/9.0
    elif direction == 'gaussian':
        core = np.zeros((3,3))
        kesai = 1
        for dx in range(-1,2):
            for dy in range(-1,2):
                core[dx+1,dy+1] = (1/(2*kesai*math.pi))*math.exp(-((dx**2+dy**2)/(2*kesai**2)))
        core = core/np.sum(core)
    else:
        print('error: please input the right string of direction')
    for y in range(img.shape[0]):
        for x in range(img.shape[1]):
            if direction == 'gaussian':
                for channel in range(3):
                    filter[pad+y,pad+x,channel] = np.sum(tmp[y:y+core_size,x:x+core_size,channel]*core)
            else:
                filter[pad+y,pad+x] = np.sum(tmp[y:y+core_size,x:x+core_size]*core)
    filter = filter[pad:pad+img.shape[0],pad:pad+img.shape[1]].astype(np.uint8)
    return filter

def NMS(img):
    for row in range(2,img.shape[0]-2):
        for col in range(2,img.shape[1]-2):
            for i in range(row-2,row+3):
                for j in range(col-2,col+3):
                    if i!=row and j!=col and img[i,j]<=img[row,col] and img[i,j]>0:
                        img[row,col]=0
    return img

def plt_cornor(img,array):
    plt.imshow(img)
    for row in range(array.shape[0]):
        for col in range(array.shape[1]):
            if array[row,col]>=200:
                plt.plot(col,row,color='r',marker='o',markersize=0.1)
            else:
                continue
    plt.show()


image = cv2.imread("image.png")
img_gray = cv2.cvtColor(image,cv2.COLOR_BGR2GRAY)
filterx = Filter(img_gray,'x')
filtery = Filter(img_gray,'y')

M = np.zeros((img_gray.shape[0],img_gray.shape[1],3),dtype=np.float32)
for row in range(img_gray.shape[0]):
    for col in range(img_gray.shape[1]):  
        M[row,col,0] = filterx[row,col]*filtery[row,col]    #Ix*Ix
        M[row,col,1] = filterx[row,col]*filtery[row,col]    #Ix*Iy 
        M[row,col,2] = filtery[row,col]*filtery[row,col]    #Iy*Iy 
gaussian = Filter(M,'gaussian')
k = 0.04
R = np.zeros((gaussian.shape[0],gaussian.shape[1]),dtype=np.float32)
for row in range(gaussian.shape[0]):
    for col in range(gaussian.shape[1]):       
        A = gaussian[row,col,0];  #Ix*Ix            
        C = gaussian[row,col,1];  #Ix*Iy            
        B = gaussian[row,col,2];  #Iy*Iy            
        #响应值计算公式            
        R[row,col] =(A*B-C*C)-(k*(A+B)*(A+B)); 
R = NMS(R)
R = (R/np.max(R))*255
R = R.astype(np.uint8)
print(R)
plt_cornor(image,R)
'''cv2.imshow("filter",R)
cv2.waitKey(0)
cv2.destroyAllWindows()'''

边缘检测

       可以直接用sobel、scharr求图片一阶导数，即为边缘，但准确度低（scharr优于sobel）；laplace求二阶导可以检测边缘，且具有旋转不变性的优点，但对噪声更加敏感，降噪方法不好的话，很影响边缘检测效果；canny检测效果很好，但计算量大，实现复杂，非最大值抑制需要检测邻域中沿梯度方向的最大值。 下面主要介绍canny边缘检测。

       首先在x和y方向上进行高斯差分滤波

       计算梯度方向和幅值

       进行非极大值抑制，

       滞后阈值法进行边缘连接，即定义高低阈值，使用高阈值检测边缘，使用低阈值继续检测边缘，分出强弱边缘，再从强边缘像素开始的“跟随”边缘，如果弱边缘和强边缘连通，则认为弱边缘是实际边缘的一部分。

直线检测

       hough变换即可。本质是笛卡尔系转参数系，将极坐标下表达直线的（角度、距离）作为新坐标系的（x，y）轴，可绘出正弦波，求这些正弦波最多的交点，该交点就是最佳拟合的直线。不过，异常点对检测结果影响很大，可以通过ransac随机抽样检测，但对更高级的圆、椭圆检测，ransac也达不到预期效果。

       hough变换起初对直线的表达是y=mx+b，但将其从笛卡尔系（图像空间）转到参数系（霍夫空间）后，可能会出现无限多的交点（代价太大了），而且无法检测垂直直线的检测。

       所以，用极坐标系表达直线，伪代码如下：

       除此之外，还有一些广义霍夫变换，不过参数更多，难度也更大。

import numpy as np 
import cv2 
from matplotlib import pyplot as plt 

def abline (slope, intercept):
    axes = plt.gca()
    x_vals = np.array (axes.get_xlim())
    y_vals = intercept + slope * x_vals
    plt.plot (x_vals, y_vals, '--')

# 读取图像 
img = cv2.imread('2.png')
gray = cv2.cvtColor(img,cv2.COLOR_BGR2GRAY)
edges = cv2.Canny(gray,50,150,apertureSize=3)  

lens = np.ceil(np.sqrt(edges.shape[0]**2 + edges.shape[1]**2))

acc_count = np.zeros((int(2*lens),180),dtype = np.uint64)

for row in range(edges.shape[0]):
    for col in range(edges.shape[1]):
        if edges[row,col] == 255:
            for theta in range(0,180,10):
                d = round(col*np.cos(np.deg2rad(theta))+row*np.sin(np.deg2rad(theta)))+lens
                if isinstance(d,int):
                    acc_count[d,theta] += 1
                else:
                    acc_count[d.astype(int),theta] += 1

img = cv2.cvtColor(img, cv2.COLOR_BGR2RGB)
for row in range(acc_count.shape[0]):
    for col in range(acc_count.shape[1]):
        if acc_count[row,col] > 100:
            k = -np.cos(np.deg2rad(col))/np.sin(np.deg2rad(col))
            b = (row-lens)/np.sin(np.deg2rad(col))
            plt.imshow(img)
            abline(k,b)
plt.show()