边缘检测应用之钱币检测

看了北京邮电大学鲁鹏老师讲授的计算机视觉课程，在看完前4节课后尝试做了留下的作业，可以说是对讲授内容的一次复习及应用，受益匪浅，就以此篇博客来记录过程中学到的知识与遇到的问题。

任务说明

编写一个钱币定位系统，其不仅能够检测出输入图像中各个钱币的边缘，同时，还能给出各个钱币的圆心坐标与半径。

解决思路

步骤1. 使用Canny算法提取图像边缘：

1. 为了防止边缘附近的噪声对之后进行的检测产生影响，先对图像进行高斯滤波，消除噪声
2. 计算图像的梯度图并获取梯度方向，得到粗略的边缘图像
3. 对梯度图进行非极大化抑制，得到不那么粗略的边缘图像
4. 使用双阈值法增加边缘图像的细节部分，得到最终的边缘图像

步骤2. 在边缘图上利用Hough变换计算圆心与半径

1. 根据图像空间进行变换得到以圆心的x，y坐标及圆半径r为轴的三维参数空间(之后简称参数空间)
2. 根据边缘点的梯度方向对参数空间投票
3. 依据预设置的阈值筛选出初步结果
4. 对筛选出的结果进行非极大化抑制，得到更精确的结果

具体实现

大体流程 (main.py)

main.py 负责提供大体框架

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

Path = "picture_source/picture.jpg"     #图片路径
Save_Path = "picture_result/"           #结果存放路径
Reduced_ratio = 2                       #图像缩放比例
Guassian_kernal_size = 3                #高斯滤波器的核的大小
HT_high_threshold = 45                  #双阈值法的高阈值，若检测结果图不明显，可以降低适当该值
HT_low_threshold = 25                   #双阈值法的低阈值，若检测结果图不明显，可以降低适当该值
Hough_transform_step = 7                #图像空间转参数空间的变化比例
Hough_transform_threshold = 40         #Hough变换的门限值，若检测出的圆数少于实际值，可以适当减小该值

if __name__ == '__main__':
    img_gray = cv2.imread(Path, cv2.IMREAD_GRAYSCALE)   #Canny检测输入为灰度图
    img_RGB = cv2.imread(Path)
    y, x = img_gray.shape[0:2]      #y为图像的行数，x为图像的列数
    img_gray = cv2.resize(img_gray, (int(x / Reduced_ratio), int(y / Reduced_ratio)))       #将图像大小变为原先的1/2
    img_RGB = cv2.resize(img_RGB, (int(x / Reduced_ratio), int(y / Reduced_ratio)))
    # canny takes about 40 seconds
    print ('Canny ...')     
    canny = Canny(Guassian_kernal_size, img_gray, HT_high_threshold, HT_low_threshold)      
    canny.canny_algorithm()     #进行Canny边缘检测
    cv2.imwrite(Save_Path + "MYcanny_result.jpg", canny.img)
    
    # hough takes about 30 seconds
    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 + "MYhough_result_an.jpg", img_RGB)
    print ('Finished!')

几个变量中 Reduced_ratio , HT_high_threshold , HT_low_threshold , Hough_transform_step , Hough_transform_threshold ,可以根据自己实际得到图像结果进行更改，以实现更好的结果。

ps:其实初始值我也不太清楚如何设定的，也没找到具体的说法，可能属于经验值吧

边缘检测 (my_canny.py)

求梯度图以及梯度方向矩阵 (Get_gradient_img)

Get_gradient_img(self):该函数用于计算梯度图和梯度方向矩阵，返回生成的梯度图

#self.x_kernal = np.array([[-1, 1]])
#self.y_kernal = np.array([[-1], [1]])

def Get_gradient_img(self):

        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.y):
            for j in range(0, self.x):
                if j == 0:
                    new_img_x[i][j] = 1
                else:
                    new_img_x[i][j] = np.sum(np.array([self.img[i][j - 1], self.img[i][j]]) * self.x_kernal)          # 获得x方向的梯度图
                if i == 0:
                    new_img_y[i][j] = 1
                else:
                    new_img_y[i][j] = np.sum(np.array([[self.img[i - 1][j]], [self.img[i][j]]]) * self.y_kernal)      # 获得y方向的梯度图

        gradient_img, self.angle = cv2.cartToPolar(new_img_x, new_img_y)        # cartToPolar(...) 计算梯度幅值及方向, 生成梯度图
        self.angle = np.tan(self.angle)
        self.img = gradient_img.astype(np.uint8)
    
        return self.img

该函数中：

new_img_x[i][j] = np.sum(np.array([self.img[i][j - 1], self.img[i][j]]) * self.x_kernal)

new_img_y[i][j] = np.sum(np.array([[self.img[i - 1][j]], [self.img[i][j]]]) * self.y_kernal)

两句体现了算法的核心。
要理解这个算法，就先要理解什么是边缘。

边缘，就是图像中发生变化的边界

例如，下图中沿红线方向，图像颜色由白变黑再变白，其中黑色线条与白色背景交界处就是图像变化的边界，即为边缘。

知道了边缘其实表示的是图像的变化，现在要做的就是在数学世界中找到类似可以表示变化的工具。回想这么多年的学习历程，大家应该都可以想到这个工具，他的名字就是导数。
我们所熟知的导数公式如下：

然而在计算机这个离散的（0，1）的世界中，想要表示极限很难。所以我们可以将上述式子进行小小的变化：

这样就可以通过求相邻的两像素的差值来近似求导。

这也是代码段中一开始 self.x_kernal，self.y_kernal，的用处，即利用卷积的方法计算x，y的导数。这是我觉得最妙的地方！

之后便是求图像的梯度图，梯度方向矩阵，数学公式为：

可以利用cv2.cartToPolar函数求得。

ps:一些使用的函数的用法会在本文最后整理

到此为止，我们就完成了边缘检测的第一步。

非极大化抑制 (Non_maximum_suppression)

Non_maximum_suppression(self)：该函数对图像进行非极大化抑制，返回经过非极大化抑制处理后的结果图

def Non_maximum_suppression (self):

        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:  #初定的阈值，也是经验值
                    self.img[i][j] = 0
                    continue
                elif abs(self.angle[i][j]) > 1:
                    gradient2 = self.img[i - 1][j]
                    gradient3 = self.img[i + 1][j]
                    if self.angle[i][j] > 0:            # 注意，此时的y轴正方向是向下的
                        gradient1 = self.img[i - 1][j - 1]
                        gradient4 = self.img[i + 1][j + 1]
                    else:
                        gradient1 = self.img[i - 1][j + 1]
                        gradient4 = self.img[i + 1][j - 1]
                else:
                    gradient2 = self.img[i][j - 1]
                    gradient3 = self.img[i][j + 1]
                    if self.angle[i][j] > 0:
                        gradient1 = self.img[i - 1][j - 1]
                        gradient4 = self.img[i + 1][j + 1]
                    else:
                        gradient1 = self.img[i + 1][j - 1]
                        gradient4 = self.img[i - 1][j + 1]
                dTemp1 = abs(self.angle[i][j]) * gradient1 + (1 - abs(self.angle[i][j])) * gradient2
                dTemp2 = abs(self.angle[i][j]) * gradient4 + (1 - abs(self.angle[i][j])) * gradient3
        
                if self.img[i][j] >= dTemp1 and self.img[i][j] >= dTemp2:
                    result[i][j] = self.img[i][j]
                else:
                    result[i][j] = 0
        self.img = result

        return self.img

这里就不得不说说为什么要进行这一步了。按照视频中鲁鹏老师的讲授，虽然我们通过初定的阈值过滤掉了部分值，但剩余的值并不都是我们需要的，我们只需要一段中的极大值。非极大化抑制可以帮助我们做到这一点。按我的理解就是可以使图片更加‘精细’。

了解了为什么要进行这一步后我们再来看看如何实现这一步。思路很简单，就是比大小。将一个像素与他沿梯度方向正向以及负向相邻的像素比较像素值的大小，找出其中的最大值，其余全部赋0。

上述方法中，难点在于如何获取梯度方向的像素值。由于计算机世界的离散性，可能在梯度方向上并不存在值，这时候就需要使用插值法，利用存在的店来近似出需要的点。具体的方法可以参考这篇博客，相信你看完之后就能明白了。

再记录一个困扰了我的问题，就是角度的理解。OpenCV打开图片，是以图片的左上方点为原点，沿图片水平向右为x轴正方向，沿图片水平向下为y轴正方向。所以，我常规理解的角度是逆时针旋转的角度，而在这里角度其实是顺时针旋转得到的。理解了这个再去想函数中的角度问题就没什么困难了。

滞后阈值法 (Hysteresis_thresholding)

Hysteresis_thresholding(self)：该函数用于增强图像细节，弥补非极大化抑制后的缺失

def Hysteresis_thresholding(self):

        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:   # 注意此时使用的是img_origin,因为非极大化抑制后的图中与极大值点相邻的点已经为0了，所以需要使用未经过非极大化处理的图像来判断
                            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
                        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
                        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
                        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
                        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

这时候可能会想：既然这个函数会弥补非极大化抑制造成的缺失，那么为什么要进行这两步呢？

不妨尝试下如果不进行这两步会发生什么。

这是原图：

这是经过非极大化抑制，滞后阈值处理的图像：

这是未经处理的图像：

看上去甚至是不处理的图像更加清晰!

但别着急，我们还要找出图中圆的位置

这是经过非极大化抑制，滞后阈值处理的最终结果：

这是未处理的最终结果：

那这是为什么呢？明明未处理的图像更清晰啊，为什么会是这样的结果呢？

我认为，正是这更清晰的原因导致了最终的结果。

我们将两张边缘图放大看看。

这是经过非极大化抑制，滞后阈值处理的图像放大：

这是未经处理的图像放大：

可以明显看到未处理的图像线条更粗，而这个更粗的线条中含有我们不需要的，冗余的信息，这些信息导致了最终的失败。边缘检测，边缘检测，就是为了检测边缘，这么粗的线条，他又怎么能叫做边缘呢？

所以，经过我们的实践，非极大化抑制与双阈值法还是边缘检测中不可少的部分！

拟合 (my_hough.py)

参数空间建立及投票 (Hough_transform_algorithm)

Hough_transform_algorithm(self)：建立三维空间，根据图片中边上的点沿梯度方向对空间中的所有单元进行投票，返回投票矩阵。

def Hough_transform_algorithm(self):

        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，y是边缘线上的点
                    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     # 图像空间与转换后的坐标空间之间的变换关系：图像空间中每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

有关hough变换的详细内容请参考鲁鹏老师的计算机视觉课程

选择圆 (Select_Circle)

Select_Circle(self)：该函数用于筛选出图像中的圆。无返回值。

def Select_Circle(self):
        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       # 逆变换回几个点中的一个代表点
                        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:      # 将每个点周围一定范围内的点的坐标求均值记录在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)        # 几个距离邻近的圆心的均值
                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]))

至于为什么要进行两遍选择的操作，是因为一次选择后可能还有冗余信息。而两次选择可以去除这些冗余。

具体如下：

这是两次选择的图像：

这是一次选择的图像：

可以看出一次选择的图中还是存在冗余的圆。

到这里我们的任务就算是圆满结束了！！！

函数解析

cv2.cartToPolar

#用于计算二维向量的大小和角度
CV_EXPORTS_W void cartToPolar(InputArray x, InputArray y,
                              OutputArray magnitude, OutputArray angle,
                              bool angleInDegrees = false);

• x：x轴坐标数组，必须是单精度或双精度浮点数组
• y：y轴坐标数组，大小和类型必须与x相同
• magnitude：输出与x相同大小和类型的大小数组
• angle：与x具有相同大小和类型的角度的输出数组；角度以弧度（从0到2*Pi）或度（0到360度）度量
• angleInDegrees：标志，指示结果是以弧度（默认情况下是以弧度）还是以度度量

numpy.tan

#对数组中每一个元素求正切值
numpy.tan(x, /, out=None, *, where=True, casting='same_kind', order='K', dtype=None, subok=True[, signature, extobj]) = <ufunc 'tan'>

• x：输入数组
• out：指定结果存储的位置。若提供此参数，其维度必须与输入数组广播后的维度一致。若不提供此参数或参数值为None，将返回新开辟的数组。若此参数为元组，其长度必须和返回值的个数保持一致。
• where：True用于标记进行函数计算的位置，False用于标记此位置不进行函数计算，直接将输入值原样返回，通常用默认值即可。

math.ceil | math.floor

#将x向上舍入到最接近的整数
math.ceil(x)
#将x向下舍入到最接近的整数
math.floor(x)

• x：必需，数字。如果 x 不是一个数字，返回 TypeError。

参考资料

计算机视觉（本科） 北京邮电大学 鲁鹏 清晰完整合集

边缘检测中非极大值抑制简单解释