Python代码实现霍夫变换检测圆

霍夫变换检测圆

参考资料：计算机视觉 北京邮电大学 鲁鹏

1.简单的介绍一下霍夫变换

我们在使用canny算子检测到图像的边缘之后，我们要如何去描述这些图像的边缘？

我们是否可以通过一些简单的特征来描述这些图像边缘？比如利用直线、圆等简单的集合形状来描述边缘的特征？这就是霍夫变换的目的。

霍夫变换准确来说是一种特征检测，不仅能识别直线，也能够识别任何形状，常见的有圆形、椭圆形。

2.霍夫变换检测圆的原理

一个圆所具有的的基本属性：圆心坐标$(x,y)$，半径r

已知一个圆形边缘，如何得出其圆心坐标和半径？

2.1 简单但不实用的方法

该方法我们仅用来参考其算法思想，实际中并不应用。

如图，我们用像素格来表示图片。

我们想要获得圆心坐标和半径，我们可以选取圆上任意一点A，然后以图片上每一个像素格为圆心B，以AB为半径画圆，记录下圆心坐标和半径。

重复上述过程，知道将圆上的所有点都计算一遍。（我们将圆心坐标和半径当做一个整体来看，你可以将其视为一个元组）

然后统计这些圆心坐标和半径，我们采用投票的方法，这些圆心坐标和半径每出现一次，就相当于得到一票。

这样一定会有一个圆心坐标和半径得票次数最多，那么这个圆心坐标和半径就是我们所要找的。

因为，圆上的每一点都会给真正的圆心和半径投一票。

而，以其他点为圆心所画的圆，通常只能够得到一票或者两票。

显然，这种方法的计算量很大。因此我们在实际应用中并不采用这种方法。

2.2梯度投票法

我们在理解了上述算法的思想之后，就可以考虑如何减少算法的计算量。

我们在学习canny算子的时候，就知道边缘上的每一个像素格都有自己梯度方向。而梯度方向通常和边都是垂直的，因此我们可以知道，沿圆上一点的梯度方向画一条直线，该直线一定会通过圆心。

于是，我们仅让沿边缘上一点的梯度方向的像素格作为圆心画圆，并记录圆心坐标和半径，这样我们同样可以得到真正的圆心坐标和半径。

然而，这种方法仍旧计算量很大，我们考虑设置一个步长。即沿梯度方向，每走一个步长的像素格，我们将其取为圆心，并计算半径。这样就进一步降低了计算量。如下图所示，是以一个像素格为步长所画的圆。黑色的矩形框表示我们的图片的大小，黑色的圆为我们通过canny算子检测得到的圆形边缘。红色的圆为我们所取的圆心和半径。

现在我们需要确定其圆心坐标和半径。

• 1.首先选取圆上任意一点A，获得其梯度方向（这里为了方便画图，选取了一个梯度方向为水平的点），并设定步长。

  梯度方向在我们canny算子检测边缘时就已经获得了。

• 2.我们沿该梯度方向以一定步长选取像素格作为圆心B，并计算AB之间的距离，作为半径。记录圆心坐标，半径。

• 3.对圆上每一点重复上述操作。

• 4.统计所有圆心坐标和半径，获得得票数最多的圆心坐标和半径，作为输出。

3.需要解决的几个问题

3.1非最大化抑制

由于我们取的步长并不一定使我们取到真正的圆心，所以通常统计出来的圆心坐标和半径会出现聚集的现象。

如图所示，在真正的圆心附近，出现几个得票数比较多的像素格（通常是前几名）。因此我们需要对其进行最大化抑制。

但通常情况下，我们并不是单纯的进行非最大化抑制，我们会将得票数最多的像素格，及其周围一定范围内的像素格进行平均操作，得出最准确的真实圆心坐标和半径。

3.2多个圆的输出

我们在处理一个图像的时候，通常并不能够确定这个图像中包含有多少个圆。因此我们通常会设置一个阈值，投票数多于某个阈值的点输出，低于阈值的点舍弃。

4.代码

4.1 canny算子

这里使用的canny代码并非上篇文章中所使用的代码，Python代码实现Canny算法——图像边缘轮廓提取，

为什么不是用之前所使用的代码呢？一是因为之前的代码计算效率低，二是因为我太懒了不想优化了（真实原因是我为了把上篇文章的代码融合到霍夫变换的代码中改了好久的bug，运行结果仍旧不太好），于是就取用北邮鲁鹏老师所给的样例代码￣□￣｜

因此在这里再赘述一下canny算子的代码，虽然代码不同，但算法思想是相同的：

import cv2
import numpy as np

class Canny:

    def __init__(self, Guassian_kernal_size, img, HT_high_threshold, HT_low_threshold):
        '''
        :param Guassian_kernal_size: 高斯滤波器尺寸
        :param img: 输入的图片，在算法过程中改变
        :param HT_high_threshold: 滞后阈值法中的高阈值
        :param HT_low_threshold: 滞后阈值法中的低阈值
        '''
        self.Guassian_kernal_size = Guassian_kernal_size  # 获取高斯核尺寸
        self.img = img  # 获取图片
        self.y, self.x = img.shape[0:2]  # 获取图片的维度
        self.angle = np.zeros([self.y, self.x])  # 定义一个新的矩阵，用来存储梯度方向
        self.img_origin = None
        self.x_kernal = np.array([[-1, 1]])  # x方向偏导卷积核
        self.y_kernal = np.array([[-1], [1]])  # y方向偏导卷积核
        self.HT_high_threshold = HT_high_threshold  # 高门限
        self.HT_low_threshold = HT_low_threshold  # 低门限

    def Get_gradient_img(self):
        '''
        计算梯度图和梯度方向矩阵。
        :return: 生成的梯度图
        '''
        print('Get_gradient_img')

        new_img_x = np.zeros([self.y, self.x], dtype=np.float)  # 定义x偏导存储矩阵
        new_img_y = np.zeros([self.y, self.x], dtype=np.float)  # 定义y偏导存储矩阵
        for i in range(0, self.x):
            for j in range(0, self.y):
                if j == 0:
                    new_img_y[j][i] = 1
                else:
                    new_img_y[j][i] = np.sum(np.array([[self.img[j - 1][i]], [self.img[j][i]]]) * self.y_kernal)
                    # y方向卷积后的灰度图，dy
                if i == 0:
                    new_img_x[j][i] = 1
                else:
                    new_img_x[j][i] = np.sum(np.array([self.img[j][i - 1], self.img[j][i]]) * self.x_kernal)
                    # x方向卷积后的灰度图，dx

        gradient_img, self.angle = cv2.cartToPolar(new_img_x, new_img_y)  # 通过dx，dy求取梯度强度和梯度方向即角度
        self.angle = np.tan(self.angle) 
        self.img = gradient_img.astype(np.uint8)  # 将获得梯度强度转换成无符号八位整形
        return self.img

    def Non_maximum_suppression(self):
        '''
        对生成的梯度图进行非极大化抑制，将tan值的大小与正负结合，确定离散中梯度的方向。
        :return: 生成的非极大化抑制结果图
        '''
        print('Non_maximum_suppression')

        result = np.zeros([self.y, self.x])  # 定义一个新矩阵，用来存储非极大化抑制结果图
        for i in range(1, self.y - 1):
            for j in range(1, self.x - 1):
                if abs(self.img[i][j]) <= 4:
                    result[i][j] = 0
                    continue
                elif abs(self.angle[i][j]) > 1:  # dy>dx
                    gradient2 = self.img[i - 1][j]
                    gradient4 = self.img[i + 1][j]
                    # g1 g2
                    #    C
                    #    g4 g3
                    if self.angle[i][j] > 0:
                        gradient1 = self.img[i - 1][j - 1]
                        gradient3 = self.img[i + 1][j + 1]
                    #    g2 g1
                    #    C
                    # g3 g4
                    else:
                        gradient1 = self.img[i - 1][j + 1]
                        gradient3 = self.img[i + 1][j - 1]
                else:
                    gradient2 = self.img[i][j - 1]
                    gradient4 = self.img[i][j + 1]
                    # g1
                    # g2 C g4
                    #      g3
                    if self.angle[i][j] > 0:
                        gradient1 = self.img[i - 1][j - 1]
                        gradient3 = self.img[i + 1][j + 1]
                    #      g3
                    # g2 C g4
                    # g1
                    else:
                        gradient3 = self.img[i - 1][j + 1]
                        gradient1 = self.img[i + 1][j - 1]

                temp1 = abs(self.angle[i][j]) * gradient1 + (1 - abs(self.angle[i][j])) * gradient2
                temp2 = abs(self.angle[i][j]) * gradient3 + (1 - abs(self.angle[i][j])) * gradient4
                if self.img[i][j] >= temp1 and self.img[i][j] >= temp2:
                    result[i][j] = self.img[i][j]
                else:
                    result[i][j] = 0
        self.img = result
        return self.img

    def Hysteresis_thresholding(self):
        '''
        对生成的非极大化抑制结果图进行滞后阈值法，用强边延伸弱边，这里的延伸方向为梯度的垂直方向，
        将比低阈值大比高阈值小的点置为高阈值大小，方向在离散点上的确定与非极大化抑制相似。
        :return: 滞后阈值法结果图
        '''
        print('Hysteresis_thresholding')

        for i in range(1, self.y - 1):
            for j in range(1, self.x - 1):
                if self.img[i][j] >= self.HT_high_threshold:
                    if abs(self.angle[i][j]) < 1:
                        if self.img_origin[i - 1][j] > self.HT_low_threshold:
                            self.img[i - 1][j] = self.HT_high_threshold
                        if self.img_origin[i + 1][j] > self.HT_low_threshold:
                            self.img[i + 1][j] = self.HT_high_threshold
                        # g1 g2
                        #    C
                        #    g4 g3
                        if self.angle[i][j] < 0:
                            if self.img_origin[i - 1][j - 1] > self.HT_low_threshold:
                                self.img[i - 1][j - 1] = self.HT_high_threshold
                            if self.img_origin[i + 1][j + 1] > self.HT_low_threshold:
                                self.img[i + 1][j + 1] = self.HT_high_threshold
                        #    g2 g1
                        #    C
                        # g3 g4
                        else:
                            if self.img_origin[i - 1][j + 1] > self.HT_low_threshold:
                                self.img[i - 1][j + 1] = self.HT_high_threshold
                            if self.img_origin[i + 1][j - 1] > self.HT_low_threshold:
                                self.img[i + 1][j - 1] = self.HT_high_threshold
                    else:
                        if self.img_origin[i][j - 1] > self.HT_low_threshold:
                            self.img[i][j - 1] = self.HT_high_threshold
                        if self.img_origin[i][j + 1] > self.HT_low_threshold:
                            self.img[i][j + 1] = self.HT_high_threshold
                        # g1
                        # g2 C g4
                        #      g3
                        if self.angle[i][j] < 0:
                            if self.img_origin[i - 1][j - 1] > self.HT_low_threshold:
                                self.img[i - 1][j - 1] = self.HT_high_threshold
                            if self.img_origin[i + 1][j + 1] > self.HT_low_threshold:
                                self.img[i + 1][j + 1] = self.HT_high_threshold
                        #      g3
                        # g2 C g4
                        # g1
                        else:
                            if self.img_origin[i - 1][j + 1] > self.HT_low_threshold:
                                self.img[i + 1][j - 1] = self.HT_high_threshold
                            if self.img_origin[i + 1][j - 1] > self.HT_low_threshold:
                                self.img[i + 1][j - 1] = self.HT_high_threshold
        return self.img

    def canny_algorithm(self):
        '''
        按照顺序和步骤调用以上所有成员函数。
        :return: Canny 算法的结果
        '''
        self.img = cv2.GaussianBlur(self.img, (self.Guassian_kernal_size, self.Guassian_kernal_size), 0)
        self.Get_gradient_img()
        self.img_origin = self.img.copy()
        self.Non_maximum_suppression()
        self.Hysteresis_thresholding()
        return self.img

4.2 霍夫变换检测圆的代码：

import numpy as np
import math


class Hough_transform:
    def __init__(self, img, angle, step=5, threshold=135):
        '''
        img: 输入的边缘图像
        angle: 输入的梯度方向矩阵
        step: Hough 变换步长大小
        threshold: 筛选单元的阈值
        '''
        self.img = img
        self.angle = angle
        self.y, self.x = img.shape[0:2]  
        self.radius = math.ceil(math.sqrt(self.y ** 2 + self.x ** 2))
        self.step = step
        self.vote_matrix = np.zeros(
            [math.ceil(self.y / self.step), math.ceil(self.x / self.step), math.ceil(self.radius / self.step)])
        # 投票矩阵，将原图的宽和高，半径分别除以步长，向上取整，
        self.threshold = threshold
        self.circles = []

    def Hough_transform_algorithm(self):
        '''
        按照 x,y,radius 建立三维空间，根据图片中边上的点沿梯度方向对空间中的所有单元进行投票。每个点投出来结果为一折线。
        return:  投票矩阵
        '''
        print('Hough_transform_algorithm')

        for i in range(1, self.y - 1):
            for j in range(1, self.x - 1):
                if self.img[i][j] > 0:
                    # 沿梯度正方向投票
                    y = i
                    x = j
                    r = 0
                    while y < self.y and x < self.x and y >= 0 and x >= 0: # 保证圆心在图像内
                        self.vote_matrix[math.floor(y / self.step)][math.floor(x / self.step)][
                            math.floor(r / self.step)] += 1
                        # 为了避免 y / self.step 向上取整会超出矩阵索引的情况，这里将该值向下取整
                        y = y + self.step * self.angle[i][j]
                        x = x + self.step
                        r = r + math.sqrt((self.step * self.angle[i][j]) ** 2 + self.step ** 2)
                    # 沿梯度反方向投票
                    y = i - self.step * self.angle[i][j]
                    x = j - self.step
                    r = math.sqrt((self.step * self.angle[i][j]) ** 2 + self.step ** 2)
                    while y < self.y and x < self.x and y >= 0 and x >= 0:
                        self.vote_matrix[math.floor(y / self.step)][math.floor(x / self.step)][
                            math.floor(r / self.step)] += 1
                        y = y - self.step * self.angle[i][j]
                        x = x - self.step
                        r = r + math.sqrt((self.step * self.angle[i][j]) ** 2 + self.step ** 2)
        return self.vote_matrix   # 返回投票矩阵

    def Select_Circle(self):
        '''
        按照阈值从投票矩阵中筛选出合适的圆，并作极大化抑制，这里的非极大化抑制我采
        用的是邻近点结果取平均值的方法，而非单纯的取极大值。
        return: None
        '''
        print('Select_Circle')
		# 挑选投票数大于阈值的圆
        houxuanyuan = []
        for i in range(0, math.ceil(self.y / self.step)):
            for j in range(0, math.ceil(self.x / self.step)):
                for r in range(0, math.ceil(self.radius / self.step)):
                    if self.vote_matrix[i][j][r] >= self.threshold:
                        y = i * self.step + self.step / 2   # 通过投票矩阵中的点，恢复到原图中的点，self.step / 2为补偿值
                        x = j * self.step + self.step / 2
                        r = r * self.step + self.step / 2
                        houxuanyuan.append((math.ceil(x), math.ceil(y), math.ceil(r)))
        if len(houxuanyuan) == 0:
            print("No Circle in this threshold.")
            return
        x, y, r = houxuanyuan[0]
        possible = []
        middle = []
        for circle in houxuanyuan:
            if abs(x - circle[0]) <= 20 and abs(y - circle[1]) <= 20: 
                # 设定一个误差范围（这里设定方圆20个像素以内，属于误差范围），在这个范围内的到圆心视为同一个圆心
                possible.append([circle[0], circle[1], circle[2]])
            else:
                result = np.array(possible).mean(axis=0)  # 对同一范围内的圆心，半径取均值
                middle.append((result[0], result[1], result[2]))
                possible.clear()
                x, y, r = circle
                possible.append([x, y, r])
        result = np.array(possible).mean(axis=0)  # 将最后一组同一范围内的圆心，半径取均值
        middle.append((result[0], result[1], result[2]))  # 误差范围内的圆取均值后，放入其中

        def takeFirst(elem):
            return elem[0]

        middle.sort(key=takeFirst)  # 排序 
        # 重复类似上述取均值的操作，并将圆逐个输出
        x, y, r = middle[0]
        possible = []
        for circle in middle:  
            if abs(x - circle[0]) <= 20 and abs(y - circle[1]) <= 20:
                possible.append([circle[0], circle[1], circle[2]])
            else:
                result = np.array(possible).mean(axis=0)
                print("Circle core: (%f, %f)  Radius: %f" % (result[0], result[1], result[2]))
                self.circles.append((result[0], result[1], result[2]))
                possible.clear()
                x, y, r = circle
                possible.append([x, y, r])
        result = np.array(possible).mean(axis=0)
        print("Circle core: (%f, %f)  Radius: %f" % (result[0], result[1], result[2]))
        self.circles.append((result[0], result[1], result[2]))

    def Calculate(self):
        '''
        按照算法顺序调用以上成员函数
        return: 圆形拟合结果图，圆的坐标及半径集合
        '''
        self.Hough_transform_algorithm()
        self.Select_Circle()
        return self.circles

4.3主函数代码：

import time
import cv2
import math
from my_hough import Hough_transform
from my_canny import Canny

Path = 'C:/Users/Administrator/Desktop/zhongqiu.jpg'   # 图片路径
Save_Path = 'C:/Users/Administrator/Desktop/'        # 结果保存路径
Reduced_ratio = 2    # 为了提高计算效率，将图片进行比例缩放所使用的比例值
Guassian_kernal_size = 3  
HT_high_threshold = 45
HT_low_threshold = 25
Hough_transform_step = 6
Hough_transform_threshold = 110

if __name__ == '__main__':
    start_time = time.time()
    
    img_gray = cv2.imread(Path, cv2.IMREAD_GRAYSCALE)
    img_RGB = cv2.imread(Path)
    y, x = img_gray.shape[0:2]
    img_gray = cv2.resize(img_gray, (int(x / Reduced_ratio), int(y / Reduced_ratio)))   # 图片缩放
    img_RGB = cv2.resize(img_RGB, (int(x / Reduced_ratio), int(y / Reduced_ratio)))
    
    # canny
    print('Canny ...')
    canny = Canny(Guassian_kernal_size, img_gray, HT_high_threshold, HT_low_threshold)
    canny.canny_algorithm()
    cv2.imwrite(Save_Path + "canny_result.jpg", canny.img)

    # hough
    print('Hough ...')
    Hough = Hough_transform(canny.img, canny.angle, Hough_transform_step, Hough_transform_threshold)
    circles = Hough.Calculate()
    for circle in circles:
        cv2.circle(img_RGB, (math.ceil(circle[0]), math.ceil(circle[1])), math.ceil(circle[2]), (132, 135, 239), 2)
    cv2.imwrite(Save_Path + "hough_result.jpg", img_RGB)
    print('Finished!')
    
    end_time = time.time()
    print("running time" + str(end_time - start_time))

4.4运行结果

canny提取到的轮廓
霍夫变换找到的圆（就是月亮！！在月亮上画了红圈圈）