Python手动实现Hough圆变换的示例代码

Hough圆变换的原理很多博客都已经说得非常清楚了，但是手动实现的比较少，所以本文直接贴上手动实现的代码。

这里使用的图片是一堆硬币：

 首先利用通过计算梯度来寻找边缘，代码如下：

def detect_edges(image):
   h = image.shape[0]
   w = image.shape[1]
   sobeling = np.zeros((h, w), np.float64)
   sobelx = [[-3, 0, 3],
             [-10, 0, 10],
             [-3, 0, 3]]
   sobelx = np.array(sobelx)

sobely = [[-3, -10, -3],
             [0, 0, 0],
             [3, 10, 3]]
   sobely = np.array(sobely)
   gx = 0
   gy = 0
   testi = 0
   for i in range(1, h - 1):
       for j in range(1, w - 1):
           edgex = 0
           edgey = 0
           for k in range(-1, 2):
               for l in range(-1, 2):
                   edgex += image[k + i, l + j] * sobelx[1 + k, 1 + l]
                   edgey += image[k + i, l + j] * sobely[1 + k, 1 + l]
           gx = abs(edgex)
           gy = abs(edgey)
           sobeling[i, j] = gx + gy
           # if you want to imshow ,run codes below first
           # if sobeling[i,j]>255:
           #  sobeling[i, j]=255
           # sobeling[i, j] = sobeling[i,j]/255
   return sobeling

需要注意的是，这里使用的kernel内的数值比较大，所以得到了结果图中的某些位置的数值超过255，但并不影响显示，但如果想通过cv2.imshow来显示，就需要将超过255的地方设为255即可（已经在代码中用注释标出），结果如下：

接下来就是要进行Hough圆变换，先看代码：

def hough_circles(edge_image, edge_thresh, radius_values):
   h = edge_image.shape[0]
   w = edge_image.shape[1]
   # print(h,w)
   edgimg = np.zeros((h, w), np.int64)
   for i in range(h):
       for j in range(w):
           if edge_image[i][j] > edge_thresh:
               edgimg[i][j] = 255
           else:
               edgimg[i][j] = 0

accum_array = np.zeros((len(radius_values), h, w))
   # return edgimg , []
   for i in range(h):
       print('Hough Transform进度：', i, '/', h)
       for j in range(w):
           if edgimg[i][j] != 0:
               for r in range(len(radius_values)):
                   rr = radius_values[r]
                   hdown = max(0, i - rr)
                   for a in range(hdown, i):
                       b = round(j+math.sqrt(rr*rr - (a - i) * (a - i)))
                       if b>=0 and b<=w-1:
                           accum_array[r][a][b] += 1
                           if 2 * i - a >= 0 and 2 * i - a <= h - 1:
                               accum_array[r][2 * i - a][b] += 1
                       if 2 * j - b >= 0 and 2 * j - b <= w - 1:
                           accum_array[r][a][2 * j - b] += 1
                       if 2 * i - a >= 0 and 2 * i - a <= h - 1 and 2 * j - b >= 0 and 2 * j - b <= w - 1:
                           accum_array[r][2 * i - a][2 * j - b] += 1

return edgimg, accum_array

其中输入是我们之前得到的边缘图，以及确定强边缘的阈值，以及一个包含着我们估计的半径的数组；返回值是强边缘图以及参数域矩阵。代码中首先遍历边缘图，通过阈值留下那些较强的位置，这里的阈值需要自己根据自己的输入图进行调节。接着就是进行Hough变换，这里的候选半径集合需要根据自己的输入图进行调节。在绘制参数域的过程中，只遍历了所需正方形区域（大小为 r*r）的 1/4，这是因为在坐出参数域上的一个点之后，由于圆的对称性，就可以找到与之对称的另外三个点，无需额外进行遍历。

最后一步就是从参数域矩阵中提取出结果圆，代码如下，其中筛选阈值需要根据你的输入图像自己调节：

def find_circles(image, accum_array, radius_values, hough_thresh):
   returnlist = []
   hlist = []
   wlist = []
   rlist = []
   returnimg = deepcopy(image)
   for r in range(accum_array.shape[0]):
       print('Find Circles 进度：', r, '/', accum_array.shape[0])
       for h in range(accum_array.shape[1]):
           for w in range(accum_array.shape[2]):
               if accum_array[r][h][w] > hough_thresh:

tmp = 0
                   for i in range(len(hlist)):
                       if abs(w-wlist[i])<10 and abs(h-hlist[i])<10:
                           tmp = 1
                           break

if tmp == 0:
                       #print(accum_array[r][h][w])
                       rr = radius_values[r]
                       flag = '(h,w,r)is:(' + str(h) + ',' + str(w) + ',' + str(rr) + ')'
                       returnlist.append(flag)
                       hlist.append(h)
                       wlist.append(w)
                       rlist.append(rr)

print('圆的数量:', len(hlist))

for i in range(len(hlist)):
       center = (wlist[i], hlist[i])
       rr = rlist[i]

color = (0, 255, 0)
       thickness = 2
       cv2.circle(returnimg, center, rr, color, thickness)

return returnlist, returnimg

注意一下在这一步中需要将那些圆心相近的圆剔除掉，只保留一个结果。

接着是main函数，这没啥好说的：

def main(argv):
   img_name = argv[0]

img = cv2.imread('data/' + img_name + '.png', cv2.IMREAD_COLOR)
   # print(img.shape[0], img.shape[1])
   gray_image = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)

# print(gray_image.shape[0], gray_image.shape[1])
   img1 = detect_edges(gray_image)
   cv2.imwrite('output/' + img_name + "_after_find_detect.png", img1)

thresh = 1500
   # 需要注意的是，在img1中有些地方的像素值是高于255的，这是由于之前的kernel内的数更大
   # 但这并不影响图像的显示
   # 因此这里的thresh要大于255
   radius_values = []
   for i in range(10):
       radius_values.append(20 + i)

edgeimg, accum_array = hough_circles(img1, thresh, radius_values)
   cv2.imwrite('output/' + img_name + "_after_binary.png", edgeimg)
   # Findcircle
   hough_thresh = 70
   resultlist, resultimg = find_circles(img, accum_array, radius_values, hough_thresh)

print(resultlist)
   cv2.imwrite('output/' + img_name + "_circles.png", resultimg)

if __name__ == '__main__':
   sys.argv.append("coins")
   main(sys.argv[1:])
   # TODO

下面是我的运行结果：

来源：https://blog.csdn.net/aRossoneri/article/details/122527480