Python圖像處理之Hough圓形檢測

更新時間：2023年07月27日 16:34:55 作者：菜菜的小粉豬

霍夫變換是一種特征檢測(feature extraction)，被廣泛應(yīng)用在圖像分析，本文將利用Hough變換實(shí)現(xiàn)圓形檢測，感興趣的小伙伴可以跟隨小編一起了解一下

hough檢測原理

點(diǎn)擊圖像處理之Hough變換檢測直線查看

下面直接描述檢測圓形的方法

基于Hough變換的圓形檢測方法

對于一個半徑為r，圓心為（a,b）的圓，我們將其表示為：

(x−a)²+(y−b)²=r²

此時 x=[x,y]^T，a=[a,b,r]^T，其參數(shù)空間為三維。顯然，圖像空間上的一點(diǎn) （x,y），在參數(shù)空間中對應(yīng)著一個圓錐，如下圖所示。

而圖像空間的一個圓就對應(yīng)著這一簇圓錐相交的一個點(diǎn)，這個特定點(diǎn)在參數(shù)空間的三維參數(shù)一定，就表示一定半徑一定圓心坐標(biāo)的圖像空間的那個圓。

上述方法是經(jīng)典的Hough圓檢測方法的原理，它具有精度高，抗干擾能力強(qiáng)等優(yōu)點(diǎn)，但由于該方法的參數(shù)空間為三維，要在三維空間上進(jìn)行證據(jù)累計(jì)的話，需要的時間和空間都是龐大的，在實(shí)際應(yīng)用中不適用。為加快Hough變換檢測圓的速度，學(xué)者們進(jìn)行了大量研究，也出現(xiàn)了很多改進(jìn)的Hough變換檢測圓的方法。如利用圖像梯度信息的Hough變換，對圓的標(biāo)準(zhǔn)方程對x求導(dǎo)得到下式：

從上式看出，此時的參數(shù)空間從半徑r，圓心（a,b）三維，變成了只有圓心（a,b）的二維空間，利用這種方法檢測圓其計(jì)算量明顯減少了。但這種改進(jìn)的Hough變換檢測圓的方法其檢測精度并不高，原因在于，此種方法利用了邊界斜率。從本質(zhì)上講，邊界斜率其實(shí)是用曲線在某一點(diǎn)的弦的斜率來代替的，這種情況下，要保證不存在誤差，只有在弦長為零的情況。但在數(shù)字圖像中，曲線的表現(xiàn)形式是離散的，其在某一點(diǎn)處的斜率指的是此點(diǎn)右向n步斜率或是左向n步斜率。如果弦長過小了，斜率的量化誤差就會增大。這種方法比較適用于干擾較少的完整圓形目標(biāo)。

主要代碼

def AHTforCircles(edge,center_threhold_factor = None,score_threhold = None,min_center_dist = None,minRad = None,maxRad = None,center_axis_scale = None,radius_scale = None,halfWindow = None,max_circle_num = None):
    if center_threhold_factor == None:
        center_threhold_factor = 10.0
    if score_threhold == None:
        score_threhold = 15.0
    if min_center_dist == None:
        min_center_dist = 80.0
    if minRad == None:
        minRad = 0.0
    if maxRad == None:
        maxRad = 1e7*1.0
    if center_axis_scale == None:
        center_axis_scale = 1.0
    if radius_scale == None:
        radius_scale = 1.0
    if halfWindow == None:
        halfWindow = 2
    if max_circle_num == None:
        max_circle_num = 6
    min_center_dist_square = min_center_dist**2
    sobel_kernel_y = np.array([[-1.0, -2.0, -1.0], [0.0, 0.0, 0.0], [1.0, 2.0, 1.0]])
    sobel_kernel_x = np.array([[-1.0, 0.0, 1.0], [-2.0, 0.0, 2.0], [-1.0, 0.0, 1.0]])
    edge_x = convolve(sobel_kernel_x,edge,[1,1,1,1],[1,1])
    edge_y = convolve(sobel_kernel_y,edge,[1,1,1,1],[1,1])
    center_accumulator = np.zeros((int(np.ceil(center_axis_scale*edge.shape[0])),int(np.ceil(center_axis_scale*edge.shape[1]))))
    k = np.array([[r for c in range(center_accumulator.shape[1])] for r in range(center_accumulator.shape[0])])
    l = np.array([[c for c in range(center_accumulator.shape[1])] for r in range(center_accumulator.shape[0])])
    minRad_square = minRad**2
    maxRad_square = maxRad**2
    points = [[],[]]
    edge_x_pad = np.pad(edge_x,((1,1),(1,1)),'constant')
    edge_y_pad = np.pad(edge_y,((1,1),(1,1)),'constant')
    Gaussian_filter_3 = 1.0 / 16 * np.array([(1.0, 2.0, 1.0), (2.0, 4.0, 2.0), (1.0, 2.0, 1.0)])
    for i in range(edge.shape[0]):
        for j in range(edge.shape[1]):
            if not edge[i,j] == 0:
                dx_neibor = edge_x_pad[i:i+3,j:j+3]
                dy_neibor = edge_y_pad[i:i+3,j:j+3]
                dx = (dx_neibor*Gaussian_filter_3).sum()
                dy = (dy_neibor*Gaussian_filter_3).sum()
                if not (dx == 0 and dy == 0):
                    t1 = (k/center_axis_scale-i)
                    t2 = (l/center_axis_scale-j)
                    t3 = t1**2 + t2**2
                    temp = (t3 > minRad_square)&(t3 < maxRad_square)&(np.abs(dx*t1-dy*t2) < 1e-4)
                    center_accumulator[temp] += 1
                    points[0].append(i)
                    points[1].append(j)
    M = center_accumulator.mean()
    for i in range(center_accumulator.shape[0]):
        for j in range(center_accumulator.shape[1]):
            neibor = \
                center_accumulator[max(0, i - halfWindow + 1):min(i + halfWindow, center_accumulator.shape[0]),
                max(0, j - halfWindow + 1):min(j + halfWindow, center_accumulator.shape[1])]
            if not (center_accumulator[i,j] >= neibor).all():
                center_accumulator[i,j] = 0
                                                                        # 非極大值抑制
    plt.imshow(center_accumulator,cmap='gray')
    plt.axis('off')
    plt.show()
    center_threshold = M * center_threhold_factor
    possible_centers = np.array(np.where(center_accumulator > center_threshold))  # 閾值化
    sort_centers = []
    for i in range(possible_centers.shape[1]):
        sort_centers.append([])
        sort_centers[-1].append(possible_centers[0,i])
        sort_centers[-1].append(possible_centers[1,i])
        sort_centers[-1].append(center_accumulator[sort_centers[-1][0],sort_centers[-1][1]])
    sort_centers.sort(key=lambda x:x[2],reverse=True)
    centers = [[],[],[]]
    points = np.array(points)
    for i in range(len(sort_centers)):
        radius_accumulator = np.zeros(
            (int(np.ceil(radius_scale * min(maxRad, np.sqrt(edge.shape[0] ** 2 + edge.shape[1] ** 2)) + 1))),dtype=np.float32)
        if not len(centers[0]) < max_circle_num:
            break
        iscenter = True
        for j in range(len(centers[0])):
            d1 = sort_centers[i][0]/center_axis_scale - centers[0][j]
            d2 = sort_centers[i][1]/center_axis_scale - centers[1][j]
            if d1**2 + d2**2 < min_center_dist_square:
                iscenter = False
                break
        if not iscenter:
            continue
        temp = np.sqrt((points[0,:] - sort_centers[i][0] / center_axis_scale) ** 2 + (points[1,:] - sort_centers[i][1] / center_axis_scale) ** 2)
        temp2 = (temp > minRad) & (temp < maxRad)
        temp = (np.round(radius_scale * temp)).astype(np.int32)
        for j in range(temp.shape[0]):
            if temp2[j]:
                radius_accumulator[temp[j]] += 1
        for j in range(radius_accumulator.shape[0]):
            if j == 0 or j == 1:
                continue
            if not radius_accumulator[j] == 0:
                radius_accumulator[j] = radius_accumulator[j]*radius_scale/np.log(j) #radius_accumulator[j]*radius_scale/j
        score_i = radius_accumulator.argmax(axis=-1)
        if radius_accumulator[score_i] < score_threhold:
            iscenter = False
        if iscenter:
            centers[0].append(sort_centers[i][0]/center_axis_scale)
            centers[1].append(sort_centers[i][1]/center_axis_scale)
            centers[2].append(score_i/radius_scale)
    centers = np.array(centers)
    centers = centers.astype(np.float64)
    return centers