OpenCV计算机视觉学习（6）——图像梯度计算&边缘检测（Sobel算子，scharr算法，Laplacian算子，Canny算子）

如果需要处理的原图及代码，请移步小编的GitHub地址

　　传送门：请点击我

　　如果点击有误：https://github.com/LeBron-Jian/ComputerVisionPractice

　　本文学习利用python学习边缘检测的滤波器，首先读入的图片代码如下：

import cv2
from pylab import *

img = cv2.imread("construction.jpg")
img = cv2.cvtColor(img,cv2.COLOR_BGR2RGB)
plt.imshow(img)
plt.axis("off")
plt.show()

 　　图片如下：

　　下面列一下本节学习的内容目录（这里于2020.6.9修改）：

　　前言：边缘检测的定义和类型

        1，一阶微分算子
            1.1  Roberts交叉梯度算子
            1.2  Prewitt算子
            1.3  Sobel算子
            1.4  Isotropic Sobel算子
            1.5  Scharr算子
            1.6  Sobel算子，Roberts算子，Prewitt算子的比较

        2，二阶微分算子
            2.1  Laplacian算子
            2.2  LOG算子

        3，非微分边缘检测算子——Canny算子
            3.1  Canny算子边缘检测基本原理
            3.2  Canny算子算法步骤
            3.3  Canny算子python实现

        4，降噪后进行边缘检测
            4.1   边缘检测示例1
            4.2   边缘检测示例2
            4.3   边缘检测示例3

前言：边缘检测的定义和类型

　　边缘检测是图像处理和计算机视觉的基本问题，边缘检测的目的是标识数字图像中亮度变化明显的点，图像属性中的显著变化通常反映了属性的重要事件和变化。这些包括：深度上的不连续，表面方向的不连续，物质属性变化和场景照明变化。边缘检测是图像处理和计算机视觉中，尤其是特征提取中的一个研究领域。图像边缘检测大幅度的减少了数据量，并且剔除了可以认为不相关的信息，保留了图像重要的结构属性。

　　在实际的图像分割中，往往只用到一阶和二阶导数，虽然原理上，可以用更高阶的导数，但是因为噪声的影响，在纯粹二阶的导数操作中就会出现对噪声的敏感现象，三阶以上的导数信息往往失去了应用价值。二阶导数还可以说明灰度突变的类型。在某些情况下，如灰度变化均匀的图像，只利用一阶导数可能找不到边界，此时二阶导数就能提供很有用的信息。二阶导数对噪声也比较敏感，解决的方法是先对图像进行平滑滤波，消除部分噪声，再进行边缘检测。不过，利用二阶导数信息的算法是基于过零检测的，因此得到的边缘点数比较少，有利于后继的处理和识别工作。

　　边缘类型：简单分为四种类型，阶跃型，屋脊型，斜坡型，脉冲型，其中阶跃型和斜坡型是类似的，只是变化的快慢不同。

　　人类视觉系统认识目标的过程分为两步:首先，把图像边缘与背景分离出来;然后，才能知觉到图像的细节，辨认出图像的轮廓。计算机视觉正是模仿人类视觉的这个过程。因此在检测物体边缘时，先对其轮廓点进行粗略检测，然后通过链接规则把原来检测到的轮廓点连接起来，同时也检测和连接遗漏的边界点及去除虚假的边界点。图像的边缘是图像的重要特征，是计算机视觉、模式识别等的基础，因此边缘检测是图象处理中一个重要的环节。然而，边缘检测又是图象处理中的一个难题，由于实际景物图像的边缘往往是各种类型的边缘及它们模糊化后结果的组合，且实际图像信号存在着噪声。噪声和边缘都属于高频信号，很难用频带做取舍。
　　这就需要边缘检测来进行解决的问题了。边缘检测的基本方法有很多，一阶的有Roberts Cross算子，Prewitt算子，Sobel算子，Krisch算子，罗盘算子；而二阶的还有Marr-Hildreth算子（又称为LOG算子），在梯度方向的二阶导数过零点，而Canny算子属于非微分边缘检测算子。各种算子的存在就是对这种导数分割原理进行的实例化计算，是为了在计算过程中直接使用的一种计算单位。在对图像的操作，我们采用模板对原图像进行卷积运算，从而达到我们想要的效果。而获取一幅图像的梯度就转化为：模板（Roberts、Prewitt、Sobel、Lapacian算子）对原图像进行卷积。

　　所以我们可以用一幅图来对这些算子进行比较。

一：一阶微分算子

1.1 Roberts交叉梯度算子

　　Roberts 算子又称为交叉微分算子，它是基于交叉差分的梯度算法，通过局部差分计算检测边缘线条。常用来处理具有陡峭的低噪声图像，当图像边缘接近于正 45 度或负 45 度时，该算法处理效果更理想。其缺点是对边缘的定位不太准确，提取的边缘线条较粗。

1.1.1  Roberts交叉微分算子原理：

　　Roberts交叉梯度算子的模板分为水平方向和垂直方向，由两个2*2的模版构成，如图：

　　从上面模板中可以看出，Roberts算子能较好的增强正负 45度 的图像边缘。

　　详细计算公式如下所示：

　　对于图像来说，是一个二维的离散型数集，通过推广二维连续型求函数偏导的方法，来求得图像的偏导数，即在（x, y）处的最大变换率，也就是这里的梯度：

　　梯度是一个矢量，则（x,  y）处的梯度表示为：

 　　其大小为：

 　　平方和平方根需要大量的计算开销，所以使用绝对值来近似梯度幅值：

　　方向与 α (x,  y) 正交：

 　　对应的模板为：

　　上图是图像的垂直和水平梯度，但是我们有时候也需要对角线方向的梯度，定义如下：

 　　对应的模板为：

 　　2*2 大小的模板在概念上很简单，但是他们对于用关于中心点对称的模板来计算边缘方向不是很有用，其最小模板大小为3*3， 3*3模板考虑了中心点对段数据的性质，并携带有关于边缘方向的更多信息。

1.1.2  Roberts交叉微分算子Python实现

　　在Python中，Roberts算子主要通过 Numpy 定义模板，再调用 OpenCV的 filter2D() 函数实现边缘提取。该函数主要是利用内核实现对图像的卷积运算，其函数原型如下所示：

dst = filter2D(src, ddepth, kernel[, dst[, anchor[, delta[, borderType]]]])

变量解释：
    src表示输入图像

    dst表示输出的边缘图，其大小和通道数与输入图像相同

    ddepth表示目标图像所需的深度

    kernel表示卷积核，一个单通道浮点型矩阵

    anchor表示内核的基准点，其默认值为（-1，-1），位于中心位置

    delta表示在储存目标图像前可选的添加到像素的值，默认值为0

    borderType表示边框模式

　　Python代码实现：

# -*- coding: utf-8 -*-
import cv2 
import numpy as np 
import matplotlib.pyplot as plt
  
#读取图像
img = cv2.imread('lena.jpg')
lenna_img = cv2.cvtColor(img,cv2.COLOR_BGR2RGB)
 
#灰度化处理图像
grayImage = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
  
#Roberts算子
kernelx = np.array([[-1,0],[0,1]], dtype=int)
kernely = np.array([[0,-1],[1,0]], dtype=int)
x = cv2.filter2D(grayImage, cv2.CV_16S, kernelx)
y = cv2.filter2D(grayImage, cv2.CV_16S, kernely)
#转uint8
absX = cv2.convertScaleAbs(x)     
absY = cv2.convertScaleAbs(y)   
Roberts = cv2.addWeighted(absX,0.5,absY,0.5,0)
 
#用来正常显示中文标签
plt.rcParams['font.sans-serif']=['SimHei']
 
#显示图形
titles = [u'原始图像', u'Roberts算子'] 
images = [lenna_img, Roberts] 
for i in range(2): 
   plt.subplot(1,2,i+1), plt.imshow(images[i], 'gray') 
   plt.title(titles[i]) 
   plt.xticks([]),plt.yticks([]) 
plt.show()

　　效果如下：

　　不调库的代码实现：

import cv2
from pylab import *

saber  = cv2.imread("construction.jpg")
# 首先将原图像进行边界扩展，并将其转换为灰度图。
gray_saber = cv2.cvtColor(saber,cv2.COLOR_RGB2GRAY)
gray_saber = cv2.resize(gray_saber,(200,200))

def RobertsOperator(roi):
    operator_first = np.array([[-1,0],[0,1]])
    operator_second = np.array([[0,-1],[1,0]])
    return np.abs(np.sum(roi[1:,1:]*operator_first))+np.abs(np.sum(roi[1:,1:]*operator_second))

def RobertsAlogrithm(image):
    image = cv2.copyMakeBorder(image,1,1,1,1,cv2.BORDER_DEFAULT)
    for i in range(1,image.shape[0]):
        for j in range(1,image.shape[1]):
            image[i,j] = RobertsOperator(image[i-1:i+2,j-1:j+2])
    return image[1:image.shape[0],1:image.shape[1]]

Robert_saber = RobertsAlogrithm(gray_saber)
plt.imshow(Robert_saber,cmap="binary")
plt.axis("off")
plt.show()

 　　结果演示：

1.2  Prewitt算子

　　Prewitt 是一种图像边缘检测的微分算子，其原理是利用特定区域内像素灰度值产生的差分实现边缘检测。由于 Prewitt 算子采用 3*3 模板对区域内的像素值进行计算，而Robert算子的模板是 2*2，故 Prewitt 算子的边缘检测结果在水平方向和垂直方向均比 Robert 算子更加明显，Prewitt算子适合用来识别噪声较多，灰度渐变的图像。

1.2.1  Prewitt算子原理

       Prewitt算子是一种一阶微分算子的边缘检测，利用像素点上下、左右邻点的灰度差，在边缘处达到极值检测边缘，去掉部分伪边缘，对噪声具有平滑作用 。其原理是在图像空间利用两个方向模板与图像进行邻域卷积来完成的，这两个方向模板一个检测水平边缘，一个检测垂直边缘。

 　　对数字图像f(x，y)，Prewitt算子的定义如下：

　　　　　　G(i)=|[f(i-1,j-1)+f(i-1,j)+f(i-1，j+1)]-[f(i+1,j-1)+f(i+1，j)+f(i+1，j+1)]|

　　　　　　G(j)=|[f(i-1,j+1)+f(i,j+1)+f(i+1，j+1)]-[f(i-1,j-1)+f(i,j-1)+f(i+1，j-1)]|

 

　　则 P(i,j)=max[G(i),G(j)]或 P(i,j)=G(i)+G(j)

 

　　经典Prewitt算子认为：凡灰度新值大于或等于阈值的像素点都是边缘点。即选择适当的阈值T，若P(i,j)≥T，则(i,j)为边缘点，P(i,j)为边缘图像。这种判定是欠合理的，会造成边缘点的误判，因为许多噪声点的灰度值也很大，而且对于幅值较小的边缘点，其边缘反而丢失了。

　　Prewitt算子对噪声有抑制作用，抑制噪声的原理是通过像素平均，但是像素平均相当于对图像的低通滤波，所以Prewitt算子对边缘的定位不如Roberts算子。

 　　计算公式为：

　　因为平均能减少或消除噪声，Prewitt梯度算子法就是先求平均，再求差分来求梯度。水平和垂直梯度模板分别为：

　　　　检测水平边沿 横向模板                          

 

　　　　检测垂直平边沿 纵向模板：

　　对于如图的矩阵起始值：

　　就是以下两个式子：

　　该算子与Sobel算子类似，只是权值有所变化，但两者实现起来功能还是有差距的，据经验得知Sobel要比Prewitt更能准确检测图像边缘。

1.2.2  Prewitt算子Python实现

　　在Python中，Prewitt 算子的实现过程与 Roberts 算子比较相似。通过 Numpy定义模板，再调用OpenCV的 filter2D() 函数对图像的卷积运算，最终通过  convertScaleAbs() 和 addWeighted() 函数实现边缘提出.

　　代码如下所示：

# -*- coding: utf-8 -*-
import cv2  
import numpy as np  
import matplotlib.pyplot as plt
 
#读取图像
img = cv2.imread('lena.png')
lenna_img = cv2.cvtColor(img,cv2.COLOR_BGR2RGB)

#灰度化处理图像
grayImage = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
 
#Prewitt算子
kernelx = np.array([[1,1,1],[0,0,0],[-1,-1,-1]],dtype=int)
kernely = np.array([[-1,0,1],[-1,0,1],[-1,0,1]],dtype=int)
x = cv2.filter2D(grayImage, cv2.CV_16S, kernelx)
y = cv2.filter2D(grayImage, cv2.CV_16S, kernely)
#转uint8
absX = cv2.convertScaleAbs(x)       
absY = cv2.convertScaleAbs(y)    
Prewitt = cv2.addWeighted(absX,0.5,absY,0.5,0)

#用来正常显示中文标签
plt.rcParams['font.sans-serif']=['SimHei']

#显示图形
titles = [u'原始图像', u'Prewitt算子']  
images = [lenna_img, Prewitt]  
for i in xrange(2):  
   plt.subplot(1,2,i+1), plt.imshow(images[i], 'gray')  
   plt.title(titles[i])  
   plt.xticks([]),plt.yticks([])  
plt.show()

 　　效果如下：

　　不调用库的代码实现：

import cv2
from pylab import *

saber  = cv2.imread("construction.jpg")
# 首先将原图像进行边界扩展，并将其转换为灰度图。
gray_saber = cv2.cvtColor(saber,cv2.COLOR_RGB2GRAY)
gray_saber = cv2.resize(gray_saber,(200,200))


def PreWittOperator(roi, operator_type):
    if operator_type == "horizontal":
        prewitt_operator = np.array([[-1, -1, -1], [0, 0, 0], [1, 1, 1]])
    elif operator_type == "vertical":
        prewitt_operator = np.array([[-1, 0, 1], [-1, 0, 1], [-1, 0, 1]])
    else:
        raise ("type Error")
    result = np.abs(np.sum(roi * prewitt_operator))
    return result


def PreWittAlogrithm(image, operator_type):
    new_image = np.zeros(image.shape)
    image = cv2.copyMakeBorder(image, 1, 1, 1, 1, cv2.BORDER_DEFAULT)
    for i in range(1, image.shape[0] - 1):
        for j in range(1, image.shape[1] - 1):
            new_image[i - 1, j - 1] = PreWittOperator(image[i - 1:i + 2, j - 1:j + 2], operator_type)
    new_image = new_image * (255 / np.max(image))
    return new_image.astype(np.uint8)


plt.subplot(121)
plt.title("horizontal")
plt.imshow(PreWittAlogrithm(gray_saber,"horizontal"),cmap="binary")
plt.axis("off")
plt.subplot(122)
plt.title("vertical")
plt.imshow(PreWittAlogrithm(gray_saber,"vertical"),cmap="binary")
plt.axis("off")
plt.show()

 　　结果演示：

　　下面试一下Prewitt对噪声的敏感性

　　代码：

import cv2
from pylab import *

saber  = cv2.imread("construction.jpg")
# 首先将原图像进行边界扩展，并将其转换为灰度图。
gray_saber = cv2.cvtColor(saber,cv2.COLOR_RGB2GRAY)
gray_saber = cv2.resize(gray_saber,(200,200))


def PreWittOperator(roi, operator_type):
    if operator_type == "horizontal":
        prewitt_operator = np.array([[-1, -1, -1], [0, 0, 0], [1, 1, 1]])
    elif operator_type == "vertical":
        prewitt_operator = np.array([[-1, 0, 1], [-1, 0, 1], [-1, 0, 1]])
    else:
        raise ("type Error")
    result = np.abs(np.sum(roi * prewitt_operator))
    return result


def PreWittAlogrithm(image, operator_type):
    new_image = np.zeros(image.shape)
    image = cv2.copyMakeBorder(image, 1, 1, 1, 1, cv2.BORDER_DEFAULT)
    for i in range(1, image.shape[0] - 1):
        for j in range(1, image.shape[1] - 1):
            new_image[i - 1, j - 1] = PreWittOperator(image[i - 1:i + 2, j - 1:j + 2], operator_type)
    new_image = new_image * (255 / np.max(image))
    return new_image.astype(np.uint8)

def noisy(noise_typ,image):
    if noise_typ == "gauss":
        row,col,ch= image.shape
        mean = 0
        var = 0.1
        sigma = var**0.5
        gauss = np.random.normal(mean,sigma,(row,col,ch))
        gauss = gauss.reshape(row,col,ch)
        noisy = image + gauss
        return noisy
    elif noise_typ == "s&p":
        row,col,ch = image.shape
        s_vs_p = 0.5
        amount = 0.004
        out = np.copy(image)
        num_salt = np.ceil(amount * image.size * s_vs_p)
        coords = [np.random.randint(0, i - 1, int(num_salt))
              for i in image.shape]
        out[coords] = 1
        num_pepper = np.ceil(amount* image.size * (1. - s_vs_p))
        coords = [np.random.randint(0, i - 1, int(num_pepper)) for i in image.shape]
        out[coords] = 0
        return out
    elif noise_typ == "poisson":
        vals = len(np.unique(image))
        vals = 2 ** np.ceil(np.log2(vals))
        noisy = np.random.poisson(image * vals) / float(vals)
        return noisy
    elif noise_typ =="speckle":
        row,col,ch = image.shape
        gauss = np.random.randn(row,col,ch)
        gauss = gauss.reshape(row,col,ch)
        noisy = image + image * gauss
        return noisy

dst = noisy("s&p",saber)
plt.subplot(131)
plt.title("add noise")
plt.axis("off")
plt.imshow(dst)

plt.subplot(132)
plt.title("Prewitt Process horizontal")
plt.axis("off")
plt.imshow(PreWittAlogrithm(gray_saber,"horizontal"),cmap="binary")

plt.subplot(133)
plt.title("Prewitt Process vertical")
plt.axis("off")
plt.imshow(PreWittAlogrithm(gray_saber,"vertical"),cmap="binary")
plt.show()

 　　结果演示：

　　选择水平梯度或垂直梯度从上图可以看出对于边缘的影响还是相当大的.

1.3  Sobel算子

　　Sobel算子是一种用于边缘检测的离散微分算子，它结合了高斯平滑和微分求导。该算子用于计算图像明暗程度近似值。根据图像边缘旁边明暗程度把该区域内超过某个数的特定点记为边缘。Sobel 算子在Prewitt算子的基础上增加了权重的概念，认为相邻点的距离远近对当前像素点的影响是不同的，距离越近的像素点对应当前像素的影响越大，从而实现图像锐化并突出边缘轮廓。

　　Sobel算子的边缘定位更准确，常用于噪声较多，灰度渐变的图像。

1.3.1  Sobel算子原理

　　Sobel算子主要用于边缘检测，在技术上它是以离散型的差分算子，用来运算图像亮度函数的梯度的近似值， Sobel算子是典型的基于一阶导数的边缘检测算子，由于该算子中引入了类似局部平均的运算，因此对噪声具有平滑作用，能很好的消除噪声的影响。Sobel算子是在Prewitt算子的基础上改进的，在中心系数上使用一个权值，与Prewitt算子、Roberts算子相比因此效果更好，能较好的抑制（平滑）噪声。 
　　Sobel算子包含两组3x3的矩阵，分别为横向及纵向模板，将之与图像作平面卷积，即可分别得出横向及纵向的亮度差分近似值。

　　计算公式为：

　　实际使用中，常用如下两个模板来检测图像边缘。

　　检测水平边沿 横向模板 ：        

　　检测垂直平边沿 纵向模板：       

　　Sobel 算子根据像素点上下，左右邻点灰度加权差，在边缘处达到极值这一现象检测边缘，对噪声具有平滑作用，提供较为精确的边缘方向信息。因为Sobel算子结合了高斯平滑和微分求导（分化），因此结果会具有更多的抗噪性，当对精度要求不是很高时，Sobel 算子是一种较为常用的边缘 检测方法。

　　图像的每一个像素的横向及纵向梯度近似值可用以下的公式结合，来计算梯度的大小。

　　然后可用以下公式计算梯度方向。 

　　在以上例子中，如果以上的角度Θ等于零，即代表图像该处拥有纵向边缘，左方较右方暗。 
　　缺点是Sobel算子并没有将图像的主题与背景严格地区分开来，换言之就是Sobel算子并没有基于图像灰度进行处理，由于Sobel算子并没有严格地模拟人的视觉生理特征，所以提取的图像轮廓有时并不能令人满意。

1.3.2  Sobel算子的Python实现

　　Sobel算子依然是一种过滤器，只是其是带有方向的。在OpenCV-Python中，使用Sobel的算子的函数原型如下：

dst = cv2.Sobel(src, ddepth, dx, dy[, dst[, ksize[, scale[, delta[, borderType]]]]])

　　参数解释：

前四个是必须的参数：

• dst 表示输出的边缘图，其大小和通道数与输入图像相同
• src 表示需要处理的图像；
• ddepth 表示图像的深度，-1表示采用的是与原图像相同的深度。目标图像的深度必须大于等于原图像的深度；
• dx和dy表示的是求导的阶数，dx 表示x方向上的差分阶数，取值为1或者0，dy表示y方向上的差分阶数，取值为1或0，0表示这个方向上没有求导，一般为0、1。

其后是可选的参数：

• ksize是Sobel算子的大小，其值必须是正数和奇数，通常为1、3、5、7。
• scale是缩放导数的比例常数，默认情况下没有伸缩系数；
• delta是一个可选的增量，将会加到最终的dst中，同样，默认情况下没有额外的值加到dst中；
• borderType是判断图像边界的模式。这个参数默认值为cv2.BORDER_DEFAULT。

　　注意：在进行Sobel算子处理之后，还需要调用convertScaleAbs() 函数计算绝对值，并将图像转换为8位图像进行显示。原因是sobel算子求导的话，白到黑是正数，但是黑到白就是负数了，所有的负数都会被截断为0，所以要取绝对值。

　　下面看一下 convertScaleAbs()函数原型：

dst = convertScaleAbs(src[, dst[, alpha[, beta]]])

    src表示原数组

    dst表示输出数组，深度为8位

    alpha表示比例因子

    beta表示原数组元素按比例缩放后添加的值

　　在OpenCV-Python中，Sobel函数的使用如下：

#coding=utf-8
import cv2
import numpy as np  
 
img = cv2.imread("D:/lena.jpg", 0)
 
x = cv2.Sobel(img,cv2.CV_16S,1,0)
y = cv2.Sobel(img,cv2.CV_16S,0,1)
 
absX = cv2.convertScaleAbs(x)   # 转回uint8
absY = cv2.convertScaleAbs(y)
 
dst = cv2.addWeighted(absX,0.5,absY,0.5,0)
 
cv2.imshow("absX", absX)
cv2.imshow("absY", absY)
 
cv2.imshow("Result", dst)
 
cv2.waitKey(0)
cv2.destroyAllWindows()

　　解释：

　　在Sobel函数的第二个参数这里使用了cv2.CV_16S。因为OpenCV文档中对Sobel算子的介绍中有这么一句：“in the case of 8-bit input images it will result in truncated derivatives”。即Sobel函数求完导数后会有负值，还有会大于255的值。而原图像是uint8，即8位无符号数，所以Sobel建立的图像位数不够，会有截断。因此要使用16位有符号的数据类型，即cv2.CV_16S。

　　在经过处理后，别忘了用convertScaleAbs()函数将其转回原来的uint8形式。否则将无法显示图像，而只是一副灰色的窗口。convertScaleAbs()的原型为：

dst = cv2.convertScaleAbs(src[, dst[, alpha[, beta]]])

　　其中可选参数alpha是伸缩系数，beta是加到结果上的一个值。结果返回uint8类型的图片。

由于Sobel算子是在两个方向计算的，最后还需要用cv2.addWeighted(...)函数将其组合起来。其函数原型为：

dst = cv2.addWeighted(src1, alpha, src2, beta, gamma[, dst[, dtype]])

　　其中alpha是第一幅图片中元素的权重，beta是第二个的权重，gamma是加到最后结果上的一个值。

　　下面举一个例子：

　　分析分别计算图像的x和y方向的sobel算子，再求和与直接计算x和y方向的sobel算子的区别，代码如下：

import cv2
import numpy as np
import matplotlib.pyplot as plt
 
 
img = cv2.imread('circle.jpg', cv2.IMREAD_GRAYSCALE)

# 白到黑是正数，黑到白就是负数，所有的负数都会被截断成0，所以要取绝对值
sobelx = cv2.Sobel(img, cv2.CV_64F, 1, 0, ksize=3)
sobelxAbs = cv2.convertScaleAbs(sobelx)

sobely = cv2.Sobel(img, cv2.CV_64F, 0, 1, ksize=3)
sobelyAbs = cv2.convertScaleAbs(sobely)

# 分别计算x和y，再求和
slbelxy = cv2.addWeighted(sobelx, 0.5, sobely, 0.5, 0)


# 不建议直接计算，这里尝试直接计算，看看效果
sobelxy1 = cv2.Sobel(img, cv2.CV_64F, 1, 1, ksize=3)
sobelxy1 = cv2.convertScaleAbs(sobelxy1)



 
result = [img, sobelxAbs, sobelyAbs, slbelxy, sobelxy1]
titles = ['origin img', 'sobel x img', 'sobel y img', 'sobel x-y img', 'sobel xy img']
for i in range(5):
    plt.subplot(2, 3, i+1), plt.imshow(result[i])
    plt.title(titles[i])
    plt.xticks([]), plt.yticks([])
plt.show()

 　　效果如下：

 　　从图中可以看到，如果分开计算效果会比直接计算效果好。（如果想要看黑白图效果，可以直接使用opencv的show函数）。也可以看看lena：

　　代码2实现Sobel算子：

import cv2
from pylab import *

saber  = cv2.imread("construction.jpg")
gray_saber = cv2.cvtColor(saber,cv2.COLOR_RGB2GRAY)
gray_saber = cv2.resize(gray_saber,(200,200))

def SobelOperator(roi, operator_type):
    if operator_type == "horizontal":
        sobel_operator = np.array([[-1, -2, -1], [0, 0, 0], [1, 2, 1]])
    elif operator_type == "vertical":
        sobel_operator = np.array([[-1, 0, 1], [-2, 0, 2], [-1, 0, 1]])
    else:
        raise ("type Error")
    result = np.abs(np.sum(roi * sobel_operator))
    return result


def SobelAlogrithm(image, operator_type):
    new_image = np.zeros(image.shape)
    image = cv2.copyMakeBorder(image, 1, 1, 1, 1, cv2.BORDER_DEFAULT)
    for i in range(1, image.shape[0] - 1):
        for j in range(1, image.shape[1] - 1):
            new_image[i - 1, j - 1] = SobelOperator(image[i - 1:i + 2, j - 1:j + 2], operator_type)
    new_image = new_image * (255 / np.max(image))
    return new_image.astype(np.uint8)

plt.subplot(121)
plt.title("horizontal")
plt.imshow(SobelAlogrithm(gray_saber,"horizontal"),cmap="binary")
plt.axis("off")
plt.subplot(122)
plt.title("vertical")
plt.imshow(SobelAlogrithm(gray_saber,"vertical"),cmap="binary")
plt.axis("off")
plt.show()

 　　结果2展示：

　　代码3实现Sobel算子：

#_*_coding:utf-8_*_
from PIL import Image
from PIL import ImageEnhance
from numpy import *
from pylab import *
from scipy.ndimage import filters

image1 = Image.open('construction.jpg').convert('L')
im = array(image1)
#soble 导数滤波器  使用 Sobel 滤波器来计算 x 和 y 的方向导数，
imx = zeros(im.shape)
# print(imx)
filters.sobel(im,1,imx)

imy = zeros(im.shape)
filters.sobel(im,0,imy)

magnitude = sqrt(imx**2 + imy**2)
# print(magnitude)
def deal_with(a):
    for i in range(len(a)):
        if a[i] <50:
            a[i] =0
        elif a[i] >200:
            a[i] =255
        # else:
        #     a[i] = 155
    return a
a = np.apply_along_axis(deal_with,1,magnitude)
result =  contour(magnitude, origin='image')
axis('equal')
axis('off')
figure()
hist(magnitude.flatten(),128)
show()

 　　结果3展示：

1.4  Isotropic Sobel算子

　　Sobel算子另一种形式是(Isotropic Sobel)算子，加权平均算子，权值反比零点与中心店的距离，当沿不同方向检测边缘时梯度幅度一致，就是通常所说的各向同性Sobel(Isotropic Sobel)算子。模板也有两个，一个是检测水平边沿的 ，另一个是检测垂直平边沿的 。各向同性Sobel算子和普通Sobel算子相比，它的位置加权系数更为准确，在检测不同方向的边沿时梯度的幅度一致。

1.5  Scharr算子

1.5.1  Scharr算子原理

　　Scharr算子与Sobel算子的不同点是在平滑部分，这里所用的平滑算子是 1/16 *[3, 10, 3]，相比于 1/4*[1, 2, 1]，中心元素占的权重更重，这可能是相对于图像这种随机性较强的信号，领域相关性不大，所以邻域平滑应该使用相对较小的标准差的高斯函数，也就是更瘦高的模板。

　　由于Sobel算子在计算相对较小的核的时候，其近似计算导数的精度比较低，比如一个3*3的Sobel算子，当梯度角度接近水平或垂直方向时，其不精确性就越发明显。Scharr算子同Sobel算子的速度一样快，但是准确率更高，尤其是计算较小核的情景，所以利用3*3滤波器实现图像边缘提取更推荐使用Scharr算子。

　　Scharr算子又称为Scharr滤波器，也是计算x或y方向上的图像差分，在OpenCV中主要配合Sobel算子的运算而存在的，下面对比一下Sobel算子和scharr算子的核函数对比：

 1.5.2  Scharr算子Python实现

　　Scharr算子和Sobel算子类似，这里简单说一下其函数用法。在OpenCV-Python中，使用Scharrl的算子的函数原型如下：

dst = cv2.Scharr(src, ddepth, dx, dy[, dst[, ksize[, scale[, delta[, borderType]]]]])

　　参数解释：

前四个是必须的参数：

• dst 表示输出的边缘图，其大小和通道数与输入图像相同
• src 表示需要处理的图像；
• ddepth 表示图像的深度，-1表示采用的是与原图像相同的深度。目标图像的深度必须大于等于原图像的深度；
• dx和dy表示的是求导的阶数，dx 表示x方向上的差分阶数，取值为1或者0，dy表示y方向上的差分阶数，取值为1或0，0表示这个方向上没有求导，一般为0、1。

其后是可选的参数：

• ksize是Sobel算子的大小，其值必须是正数和奇数，通常为1、3、5、7。
• scale是缩放导数的比例常数，默认情况下没有伸缩系数；
• delta是一个可选的增量，将会加到最终的dst中，同样，默认情况下没有额外的值加到dst中；
• borderType是判断图像边界的模式。这个参数默认值为cv2.BORDER_DEFAULT。

　　注意：在进行Scharr算子处理之后，也需要调用convertScaleAbs() 函数计算绝对值，并将图像转换为8位图像进行显示。原因是sobel算子求导的话，白到黑是正数，但是黑到白就是负数了，所有的负数都会被截断为0，所以要取绝对值。

　　这里不再赘述convertScaleAbs()函数了，直接看一个Scharr算子的实现，并将其与Sobel算子对比，我们看看效果。

　　代码如下：

# coding=utf-8
import cv2
import numpy as np

img = cv2.imread("durant.jpg", 0)
img = cv2.resize(img, (0, 0), fx=0.5, fy=0.5)

sobel_x = cv2.Sobel(img, cv2.CV_16S, 1, 0)
sobel_y = cv2.Sobel(img, cv2.CV_16S, 0, 1)
scharr_x = cv2.Scharr(img, cv2.CV_16S, 1, 0)
scharr_y = cv2.Scharr(img, cv2.CV_16S, 0, 1)

sobel_absX = cv2.convertScaleAbs(sobel_x)  # 转回uint8
sobel_absY = cv2.convertScaleAbs(sobel_y)
scharr_absX = cv2.convertScaleAbs(scharr_x)  # 转回uint8
scharr_absY = cv2.convertScaleAbs(scharr_y)

Sobel_dst = cv2.addWeighted(sobel_absX, 0.5, sobel_absY, 0.5, 0)
Scharr_dst = cv2.addWeighted(scharr_absX, 0.5, scharr_absY, 0.5, 0)

# cv2.imshow("absX", scharr_absX)
# cv2.imshow("absY", scharr_absY)
# cv2.imshow("Result", dst)

sobel_image = np.hstack((sobel_absX, sobel_absY, Sobel_dst))
scharr_image = np.hstack((scharr_absX, scharr_absY, Scharr_dst))
all_image = np.vstack((sobel_image, scharr_image))
cv2.imshow("Result", all_image)


cv2.waitKey(0)
cv2.destroyAllWindows()

　　原图如下（其实这张图可以很明显的看出求x轴和y轴方向梯度的差异）：

 　　效果如下：

 　　这里还是用我最喜欢的球星做一下，哈哈哈哈哈哈哈：

 　　上面三张图分别是Sobel算子的x轴，y轴，和对x轴和y轴加起来的效果，下面三张图分别是Scharr算子的x轴，y轴，和对x轴和y轴加起来的效果。但是上面最后面两个效果图的对比，并不能说明Scharr算子比Sobel算子效果更好，可能在具体的实际场景中才会出现那个更好。

1.6  Sobel算子，Robert算子，prewitt算子的比较

　　Sobel算子是滤波算子的形式来提取边缘，X，Y方向各用一个模板，两个模板组合起来构成一个梯度算子。X方向模板对垂直边缘影响最大，Y方向模板对水平边缘影响最大。

　　Robert算子是一种梯度算子，它用交叉的差分表示梯度，是一种利用局部差分算子寻找边缘的算子，对具有陡峭的低噪声的图像效果最好。

　　prewitt算子是加权平均算子，对噪声有抑制作用，但是像素平均相当于对图像进行的同滤波，所以prewitt算子对边缘的定位不如robert算子。

二：二阶微分算子

　　拉普拉斯（Laplacian）算子是n维欧几里德空间中的一个二阶微分算子，常用于图像增强领域和边缘提取。它通过灰度差分计算邻域内的像素，基本流程是：判断图像中心像素灰度值与它周围其他像素的灰度值，如果中心像素的灰度更高，则提升中心像素的灰度；反之降低中心像素的灰度，从而实现图像锐化操作。在算法实现过程中，Laplacian算子通过对邻域中心像素的四方向或八方向求梯度，再将梯度相加起来判断中心像素灰度与邻域内其他像素灰度的关系，最后通过梯度运算的结果对像素灰度进行调整。

2.1 Laplacian算子

　　拉普拉斯（Laplacian）算子是 n 维欧几里得空间中的一个二阶微分算子，常用于图像增强领域和边缘提取，它通过灰度差分计算领域内的像素。

2.1.1 Laplacian原理

　　Laplacian算子的基本流程是：判断图像中心像素灰度值与它周围其他像素的灰度值，如果中心像素的灰度更高，则提升中心像素的灰度；反之降低中心像素的灰度，从而实现图像锐化操作。在算法实现过程中，Laplacian算子通过对邻域中心像素的四方向或八方向求梯度，再将梯度相加起来判断中心像素灰度与邻域内其他像素灰度的关系，最后通过梯度运算的结果对像素灰度进行调整。

         Laplace算子是一种各向同性算子，二阶微分算子，具有旋转不变性。在只关心边缘的位置而不考虑其周围的象素灰度差值时比较合适。Laplace算子对孤立象素的响应要比对边缘或线的响应要更强烈，因此只适用于无噪声图象。存在噪声情况下，使用Laplacian算子检测边缘之前需要先进行低通滤波。所以，通常的分割算法都是把Laplacian算子和平滑算子结合起来生成一个新的模板。

　　一阶导数为：

 　　二阶导数为：

 　　我们这里需要的是关于x的二阶导数，故将上式中的变量减去1后，得到：

 　　在图像处理中，通过拉普拉斯模板求二阶导数，其定义如下：

　　为了更适合于数字图像处理，将该方程表示为离散形式：

　　　　   另外，拉普拉斯算子还可以表示成模板的形式，如下图所示。从模板形式容易看出，如果在图像中一个较暗的区域中出现了一个亮点，那么用拉普拉斯运算就会使这个亮点变得更亮。因为图像中的边缘就是那些灰度发生跳变的区域，所以拉普拉斯锐化模板在边缘检测中很有用。一般增强技术对于陡峭的边缘和缓慢变化的边缘很难确定其边缘线的位置。但此算子却可用二次微分正峰和负峰之间的过零点来确定，对孤立点或端点更为敏感，因此特别适用于以突出图像中的孤立点、孤立线或线端点为目的的场合。同梯度算子一样，拉普拉斯算子也会增强图像中的噪声，有时用拉普拉斯算子进行边缘检测时，可将图像先进行平滑处理。

 　　Laplacian算子分为四邻域和八邻域，四邻域是对邻域中心像素的四方向求梯度，八邻域是对八方向求梯度。其中四邻域模板如公式所示：

　　　　离散拉普拉斯算子的模板：　　　　，

　　通过模板可以发现，当邻域内像素灰度相同时，模板的卷积运算结果为0；当中心像素灰度高于邻域内其他像素的平均灰度时，模板的卷积运算结果为负数；当中心像素的灰度低于邻域内其他像素的平均灰度时，模板的卷积为正数。对卷积运算的结果用适当的衰弱因子处理并加在原中心像素上，就可以实现图像的锐化处理。

　　Laplacian算子的八邻域模板如下：

　　　　其扩展模板：　　　　　　　　　　 。

      拉式算子用来改善因扩散效应的模糊特别有效，因为它符合降制模型。扩散效应是成像过程中经常发生的现象。

      Laplacian算子一般不以其原始形式用于边缘检测，因为其作为一个二阶导数，Laplacian算子对噪声具有无法接受的敏感性；同时其幅值产生算边缘，这是复杂的分割不希望有的结果；最后Laplacian算子不能检测边缘的方向；所以Laplacian在分割中所起的作用包括：（1）利用它的零交叉性质进行边缘定位；（2）确定一个像素是在一条边缘暗的一面还是亮的一面；一般使用的是高斯型拉普拉斯算子（Laplacian of a Gaussian,LoG)，由于二阶导数是线性运算，利用LoG卷积一幅图像与首先使用高斯型平滑函数卷积改图像，然后计算所得结果的拉普拉斯是一样的。所以在LoG公式中使用高斯函数的目的就是对图像进行平滑处理，使用Laplacian算子的目的是提供一幅用零交叉确定边缘位置的图像；图像的平滑处理减少了噪声的影响并且它的主要作用还是抵消由Laplacian算子的二阶导数引起的逐渐增加的噪声影响。

　　  图像锐化处理的作用是使灰度反差增强，从而使模糊图像变得更加清晰。图像模糊的实质就是图像受到平均运算或积分运算，因此可以对图像进行逆运算，如微分运算能够突出图像细节，使图像变得更为清晰。由于拉普拉斯是一种微分算子，它的应用可增强图像中灰度突变的区域，减弱灰度的缓慢变化区域。因此，锐化处理可选择拉普拉斯算子对原图像进行处理，产生描述灰度突变的图像，再将拉普拉斯图像与原始图像叠加而产生锐化图像。拉普拉斯锐化的基本方法可以由下式表示：

　　这种简单的锐化方法既可以产生拉普拉斯锐化处理的效果，同时又能保留背景信息，将原始图像叠加到拉普拉斯变换的处理结果中去，可以使图像中的各灰度值得到保留，使灰度突变处的对比度得到增强，最终结果是在保留图像背景的前提下，突现出图像中小的细节信息。

2.1.2 Laplacian算子Python实现

   　　图（a）显示了一幅花朵的图片，图（b）显示了用图（a）所示的拉普拉斯模板对该图像滤波后的结果。由图可以看出，将原始图像通过拉普拉斯变换后增强了图像中灰度突变处的对比度，使图像中小的细节部分得到增强并保留了图像的背景色调，使图像的细节比原始图像更加清晰。基于拉普拉斯变换的图像增强已成为图像锐化处理的基本工具。

　　Python和OpenCV 将Laplacian算子封装在Laplacian()函数中，其函数原型如下所示：

dst = Laplacian(src, ddepth[, dst[, ksize[, scale[, delta[, borderType]]]]])

    src表示输入图像


    dst表示输出的边缘图，其大小和通道数与输入图像相同
    ddepth表示目标图像所需的深度

    ksize表示用于计算二阶导数的滤波器的孔径大小，其值必须是正数和奇数，
且默认值为1，更多详细信息查阅getDerivKernels

    scale表示计算拉普拉斯算子值的可选比例因子。默认值为1，更多详细信息查阅getDerivKernels

    delta表示将结果存入目标图像之前，添加到结果中的可选增量值，默认值为0

    borderType表示边框模式，更多详细信息查阅BorderTypes

    注意，Laplacian算子其实主要是利用Sobel算子的运算，通过加上Sobel算子运算出的
图像x方向和y方向上的导数，得到输入图像的图像锐化结果。同时，在进行Laplacian算子
处理之后，还需要调用convertScaleAbs()函数计算绝对值，并将图像转换为8位图进行显示。

　　当ksize=1时，Laplacian()函数采用3×3的孔径（四邻域模板）进行变换处理。下面的代码是采用ksize=3的Laplacian算子进行图像锐化处理，其代码如下

# -*- coding: utf-8 -*-
import cv2  
import numpy as np  
import matplotlib.pyplot as plt
 
#读取图像
img = cv2.imread('lena.jpg')
lenna_img = cv2.cvtColor(img,cv2.COLOR_BGR2RGB)

#灰度化处理图像
grayImage = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
 
#拉普拉斯算法
dst = cv2.Laplacian(grayImage, cv2.CV_16S, ksize = 3)
Laplacian = cv2.convertScaleAbs(dst) 

#用来正常显示中文标签
plt.rcParams['font.sans-serif']=['SimHei']

#显示图形
titles = [u'原始图像', u'Laplacian算子']  
images = [lenna_img, Laplacian]  
for i in range(2):  
   plt.subplot(1,2,i+1), plt.imshow(images[i], 'gray')  
   plt.title(titles[i])  
   plt.xticks([]),plt.yticks([])  
plt.show()

　　效果如下：

　　不调库实现的代码：

import cv2
from pylab import *

saber  = cv2.imread("construction.jpg")
# 首先将原图像进行边界扩展，并将其转换为灰度图。
gray_saber = cv2.cvtColor(saber,cv2.COLOR_RGB2GRAY)
gray_saber = cv2.resize(gray_saber,(200,200))


def LaplaceOperator(roi, operator_type):
    if operator_type == "fourfields":
        laplace_operator = np.array([[0, 1, 0], [1, -4, 1], [0, 1, 0]])
    elif operator_type == "eightfields":
        laplace_operator = np.array([[1, 1, 1], [1, -8, 1], [1, 1, 1]])
    else:
        raise ("type Error")
    result = np.abs(np.sum(roi * laplace_operator))
    return result


def LaplaceAlogrithm(image, operator_type):
    new_image = np.zeros(image.shape)
    image = cv2.copyMakeBorder(image, 1, 1, 1, 1, cv2.BORDER_DEFAULT)
    for i in range(1, image.shape[0] - 1):
        for j in range(1, image.shape[1] - 1):
            new_image[i - 1, j - 1] = LaplaceOperator(image[i - 1:i + 2, j - 1:j + 2], operator_type)
    new_image = new_image * (255 / np.max(image))
    return new_image.astype(np.uint8)

def noisy(noise_typ,image):
    if noise_typ == "gauss":
        row,col,ch= image.shape
        mean = 0
        var = 0.1
        sigma = var**0.5
        gauss = np.random.normal(mean,sigma,(row,col,ch))
        gauss = gauss.reshape(row,col,ch)
        noisy = image + gauss
        return noisy
    elif noise_typ == "s&p":
        row,col,ch = image.shape
        s_vs_p = 0.5
        amount = 0.004
        out = np.copy(image)
        num_salt = np.ceil(amount * image.size * s_vs_p)
        coords = [np.random.randint(0, i - 1, int(num_salt))
              for i in image.shape]
        out[coords] = 1
        num_pepper = np.ceil(amount* image.size * (1. - s_vs_p))
        coords = [np.random.randint(0, i - 1, int(num_pepper)) for i in image.shape]
        out[coords] = 0
        return out
    elif noise_typ == "poisson":
        vals = len(np.unique(image))
        vals = 2 ** np.ceil(np.log2(vals))
        noisy = np.random.poisson(image * vals) / float(vals)
        return noisy
    elif noise_typ =="speckle":
        row,col,ch = image.shape
        gauss = np.random.randn(row,col,ch)
        gauss = gauss.reshape(row,col,ch)
        noisy = image + image * gauss
        return noisy

plt.subplot(121)
plt.title("fourfields")
plt.imshow(LaplaceAlogrithm(gray_saber,"fourfields"),cmap="binary")
plt.axis("off")
plt.subplot(122)
plt.title("eightfields")
plt.imshow(LaplaceAlogrithm(gray_saber,"eightfields"),cmap="binary")
plt.axis("off")
plt.show()

　　结果演示：

　　上图为比较取值不同的laplace算子实现的区别。

　　PIL库实现的代码：

import numpy as np
from PIL import Image
import matplotlib.pyplot as plt
import matplotlib.cm as cm
import scipy.signal as signal     # 导入sicpy的signal模块

# Laplace算子
suanzi1 = np.array([[0, 1, 0],
                    [1,-4, 1],
                    [0, 1, 0]])

# Laplace扩展算子
suanzi2 = np.array([[1, 1, 1],
                    [1,-8, 1],
                    [1, 1, 1]])

# 打开图像并转化成灰度图像
image = Image.open("construction.jpg").convert("L")
image_array = np.array(image)

# 利用signal的convolve计算卷积
image_suanzi1 = signal.convolve2d(image_array,suanzi1,mode="same")
image_suanzi2 = signal.convolve2d(image_array,suanzi2,mode="same")

# 将卷积结果转化成0~255
image_suanzi1 = (image_suanzi1/float(image_suanzi1.max()))*255
image_suanzi2 = (image_suanzi2/float(image_suanzi2.max()))*255

# 为了使看清边缘检测结果，将大于灰度平均值的灰度变成255(白色)
image_suanzi1[image_suanzi1>image_suanzi1.mean()] = 255
image_suanzi2[image_suanzi2>image_suanzi2.mean()] = 255

# 显示图像
plt.subplot(2,1,1)
plt.imshow(image_array,cmap=cm.gray)
plt.axis("off")
plt.subplot(2,2,3)
plt.imshow(image_suanzi1,cmap=cm.gray)
plt.axis("off")
plt.subplot(2,2,4)
plt.imshow(image_suanzi2,cmap=cm.gray)
plt.axis("off")
plt.show()

 　　结果：

　　其中上方为原图像，下方：左边为Laplace算子结果，右边为Laplace扩展算子结果

2.2  LOG算子

　　LOG（Laplacian of Gaussian）边缘检测算子是David Courtnay Marr和 Ellen Hildreth在 1980年共同提出，也称为 Marr & Hildreth算子。它根据图像的信噪比来求检测边缘的最优滤波器。下面学习一下其原理和应用。

2.2.1  LOG算子原理

　　LOG算子首先对图像做高斯滤波，然后再求其拉普拉斯（Laplacian）二阶导数，根据二阶导数的锅零点来检测图像的边界，即通过检测滤波结果的零交叉（Zero  crossings）来获得图像或物体的边缘。

　　LOG算子综合考虑了对噪声的抑制和对边缘的检测两个方向，并且把Gauss平滑滤波器和Laplacian锐化滤波器结合了起来，先平滑掉噪声，再进行边缘检测，所以效果会更好。该算子与视觉生理中的数学模型相似，因此在图像处理领域中得到了广泛的应用。它具有抗干扰能力强，边界定位精度高，边缘连续性好，能有效提取对比度弱的边界等特点。

　　常见的LOG算子是5*5模板，如下图所示：

 　　由于LOG算子到中心的距离与位置加权系数的关系曲线像墨西哥草帽的剖面，所以LOG算子也叫墨西哥草帽滤波器，如下图所示：

2.2.2  LOG算子Python实现

　　LOG算子的边缘提取实现代码如下：

# coding=utf-8
import cv2
import numpy as np
import matplotlib.pyplot as plt

img = cv2.imread("durant.jpg")
rgb_img = cv2.cvtColor(img, cv2.COLOR_BGR2RGB)
img = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)

# 先通过高斯滤波降噪
gaussian = cv2.GaussianBlur(img, (3, 3), 0)

# 再通过拉普拉斯算子做边缘检测
dst = cv2.Laplacian(gaussian, cv2.CV_16S, ksize=3)
LOG = cv2.convertScaleAbs(dst)

# 用来正常显示中文标签
plt.rcParams['font.sans-serif'] = ['SimHei']

# 显示图形
titles = [u'原始图像', u'LOG算子']
images = [rgb_img, LOG]
for i in range(2):
    plt.subplot(1, 2, i + 1), plt.imshow(images[i], 'gray')
    plt.title(titles[i])
    plt.xticks([]), plt.yticks([])
plt.show()

 　　结果如下：

三： 非微分边缘检测算子——Canny算子

　　John F.Canny 于 1986年发明了一个多级边缘检测算法——Canny边缘检测算子，并创立了边缘检测计算理论（Computational  theory of edge detection），该理论有效的解释了这项技术的工作理论。

3.1  Canny算子边缘检测基本原理

　　Canny边缘检测是一种比较新的边缘检测算子，具有很好地边缘检测性能，该算子功能比前面几种都要好，但是它实现起来较为麻烦，Canny算子是一个具有滤波，增强，检测的多阶段的优化算子，在进行处理前，Canny算子先利用高斯平滑滤波器来平滑图像以除去噪声，Canny分割算法采用一阶偏导的有限差分来计算梯度幅值和方向，在处理过程中，Canny算子还将经过一个非极大值抑制的过程，最后Canny算子还采用两个阈值来连接边缘（高低阈值输出二值图像）。

　　边缘检测通常是在保留原有图像属性的情况下，对图像数据规模进行缩减，提取图像边缘轮廓的处理方式。

　　(1)图象边缘检测必须满足两个条件：一能有效地抑制噪声；二必须尽量精确确定边缘的位置。

　　(2)根据对信噪比与定位乘积进行测度，得到最优化逼近算子。这就是Canny边缘检测算子。

　　(3)类似与Marr（LoG）边缘检测方法，也属于先平滑后求导数的方法。

 

3.2  Canny算子的算法步骤

　　Canny算法是一种被广泛应用于边缘检测的标准算法，其目标是找到一个最优的边缘检测解或找寻一幅图像中灰度强度变换最强的位置。最优边缘检测主要通过低错误率，高定位性和最小响应三个标准进行评价。Canny算子的实现步骤如下：

• step1: 用高斯滤波器平滑图象；
• step2: 计算图像中每个像素点的梯度强度和方向（用一阶偏导的有限差分来计算梯度的幅值和方向）；
• step3: 对梯度幅值进行非极大值抑制（Non-Maximum Suppression），以消除边缘检测带来的杂散响应；
• step4: 用双阈值算法(Double-Threshold)检测来确定真实和潜在的边缘，通过抑制孤立的弱边缘最终完成边缘检测；

 

3.2.1   step1：使用高斯平滑函数去除噪声

 

　　 3*3的高斯核函数：

 　　5*5 的高斯核函数：

3.2.2  step2：使用一阶有限查分计算偏导数阵列P和Q

　　这里需要按照Sobel滤波器步骤计算梯度幅值和方向，寻找图像的强度梯度。先将卷积模板分别作用 x 和 y 方向，再计算梯度幅值和方向，其公式如下所示。梯度方向一般取 0 度，45 度，90度和 135度四个方向。

 

 

3.2.3  step3： 非极大值抑制

　　通过非极大值抑制（Non - maximum  Suppression）过滤掉非边缘像素，将模糊的边界变得清晰。该过程保留了每隔像素点上梯度强度的极大值，过滤掉其他的值。

　　对于每个像素点，它进行如下操作：

• 1，将其梯度方向近似为以下值中的一个，包括0，45，90，135，180，225，270和315，即表示上下左右和45度方向。
• 2，比较该像素点和其梯度正负方向的像素点的梯度强度，如果该像素点梯度强度最大则保留，否则抑制（删除，即置为零）。

 

 

 

3.2.3  step4： 用双阈值算法检测和连接边缘　

　　利用双阈值方法来确定潜在的边界。经过非极大抑制后图像中仍然有很多噪声点，此时需要通过双阈值技术处理，即设定一个阈值上界和阈值下界。图像中的像素点如果大于阈值上界则认为必然是边界（称为强边界，strong  edge），小于阈值下界则认为必然不是边界，两者之间的则认为是候选项（称为弱边界，weak edge）。

　　对非极大值抑制图像作用两个阈值th1和th2，两者关系th1=0.4th2。我们把梯度值小于th1的像素的灰度值设为0，得到图像1。然后把梯度值小于th2的像素的灰度值设为0，得到图像2。由于图像2的阈值较高，去除大部分噪音，但同时也损失了有用的边缘信息。而图像1的阈值较低，保留了较多的信息，我们可以以图像2为基础，以图像1为补充来连结图像的边缘。

　　链接边缘的具体步骤如下：

　　　　对图像2进行扫描，当遇到一个非零灰度的像素p(x,y)时，跟踪以p(x,y)为开始点的轮廓线，直到轮廓线的终点q(x,y)。

　　　　考察图像1中与图像2中q(x,y)点位置对应的点s(x,y)的8邻近区域。如果在s(x,y)点的8邻近区域中有非零像素s(x,y)存在，则将其包括到图像2中，作为r(x,y)点。从r(x,y)开始，重复第一步，直到我们在图像1和图像2中都无法继续为止。

　　　　当完成对包含p(x,y)的轮廓线的连结之后，将这条轮廓线标记为已经访问。回到第一步，寻找下一条轮廓线。重复第一步、第二步、第三步，直到图像2中找不到新轮廓线为止。

　　至此，完成canny算子的边缘检测。

 

3.3 Canny算子Python实现

　　在OpenCV中，canny() 函数原型如下所示：

edges = Canny(image, threshold1, threshold2[, edges[, apertureSize[, L2gradient]]])

 　　参数含义：

• image：表示输入图像（必须是单通道图像）
• edges：表示输出的边缘图，其大小和类型与输入图像相同
• threshold1：表示第一个滞后性阈值
• threshold2：表示第二个滞后性阈值
• apertureSize：表示应用Sobel算子的孔径大小，其默认值为3
• L2gradient：表示一个计算图像梯度幅值的标识，默认值为FALSE　　

　　　　　　true：表示使用更精确的L2范数进行计算（即两个方向的倒数的平方和再开方）

　　　　　　false：表示使用L1范数（直接将两个方向导数的绝对值相加）

 　　Canny算子的边缘提取实现代码如下：

# Canny边缘提取
import cv2 as cv
import numpy as np
 
 
def edge_demo(image):
    blurred = cv.GaussianBlur(image, (3, 3), 0)
    gray = cv.cvtColor(blurred, cv.COLOR_RGB2GRAY)
    # xgrad = cv.Sobel(gray, cv.CV_16SC1, 1, 0) #x方向梯度
    # ygrad = cv.Sobel(gray, cv.CV_16SC1, 0, 1) #y方向梯度
    # edge_output = cv.Canny(xgrad, ygrad, 50, 150)
    edge_output = cv.Canny(gray, 50, 150)
    cv.imshow("Canny Edge", edge_output)
    dst = cv.bitwise_and(image, image, mask=edge_output)
    cv.imshow("Color Edge", dst)
 
 
# src = cv.imread('logo1.jpg')
# # 设置为WINDOW_NORMAL可以任意缩放
# # cv.namedWindow('input_image', cv.WINDOW_NORMAL)
# cv.imshow('input_image', src)
# edge_demo(src)
# cv.waitKey(0)
# cv.destroyAllWindows()

img = cv.imread('lena.jpg', cv.IMREAD_GRAYSCALE)

v1 = cv.Canny(img, 80, 150)
v2 = cv.Canny(img, 50, 100)

res = np.hstack((v1, img, v2))
cv.imshow('img', res)
cv.waitKey(0)
cv.destroyAllWindows()

　　结果如下：

　　lena 效果图如下：

 

　　高斯模糊代码：

import cv2
from pylab import *

saber  = cv2.imread("construction.jpg")
# 首先将原图像进行边界扩展，并将其转换为灰度图。
gray_saber = cv2.cvtColor(saber,cv2.COLOR_RGB2GRAY)
gray_saber = cv2.resize(gray_saber,(200,200))


def GaussianOperator(roi):
    GaussianKernel = np.array([[1,2,1],[2,4,2],[1,2,1]])
    result = np.sum(roi*GaussianKernel/16)
    return result

def GaussianSmooth(image):
    new_image = np.zeros(image.shape)
    image = cv2.copyMakeBorder(image,1,1,1,1,cv2.BORDER_DEFAULT)
    for i in range(1,image.shape[0]-1):
        for j in range(1,image.shape[1]-1):
            new_image[i-1,j-1] =GaussianOperator(image[i-1:i+2,j-1:j+2])
    return new_image.astype(np.uint8)

smooth_saber = GaussianSmooth(gray_saber)
plt.subplot(121)
plt.title("Origin Image")
plt.axis("off")
plt.imshow(gray_saber,cmap="gray")
plt.subplot(122)
plt.title("GaussianSmooth Image")
plt.axis("off")
plt.imshow(smooth_saber,cmap="gray")
plt.show()

　　结果演示：

 四：降噪后进行边缘检测

4.1 边缘检测示例1

 　　为了获得更好的边缘检测效果，可以先对图像进行模糊平滑处理，目的是去除图像中的高频噪声。

　　首先用标准差为5的5*5高斯算子对图像进行平滑处理，然后利用Laplace的扩展算子对图像进行边缘检测。

　　代码：

import numpy as np
from PIL import Image
import matplotlib.pyplot as plt
import matplotlib.cm as cm
import scipy.signal as signal

# 生成高斯算子的函数
def func(x,y,sigma=1):
    return 100*(1/(2*np.pi*sigma))*np.exp(-((x-2)**2+(y-2)**2)/(2.0*sigma**2))

# 生成标准差为5的5*5高斯算子
suanzi1 = np.fromfunction(func,(5,5),sigma=5)

# Laplace扩展算子
suanzi2 = np.array([[1, 1, 1],
                    [1,-8, 1],
                    [1, 1, 1]])

# 打开图像并转化成灰度图像
image = Image.open("4.JPG").convert("L")
image_array = np.array(image)

# 利用生成的高斯算子与原图像进行卷积对图像进行平滑处理
image_blur = signal.convolve2d(image_array, suanzi1, mode="same")

# 对平滑后的图像进行边缘检测
image2 = signal.convolve2d(image_blur, suanzi2, mode="same")

# 结果转化到0-255
image2 = (image2/float(image2.max()))*255

# 将大于灰度平均值的灰度值变成255（白色），便于观察边缘
image2[image2>image2.mean()] = 255

# 显示图像
plt.subplot(2,1,1)
plt.imshow(image_array,cmap=cm.gray)
plt.axis("off")
plt.subplot(2,1,2)
plt.imshow(image2,cmap=cm.gray)
plt.axis("off")
plt.show()

　　结果：

　　于 2020.1.18 号，再次修改，加深度边缘检测算子的学习。

4.2  边缘检测示例2

　　边缘检测算法主要是基于图像强度的一阶和二阶导数，但是导数通常对噪声很敏感，因此需要采用滤波器来过滤噪声，并调用图像增强或阈值化算法进行处理，最后再进行边缘检测。下面是采用高斯滤波去噪和阈值化处理之后，再进行边缘检测的过程，并对比了三种常见的边缘提取算法。

 　　代码如下：

# -*- coding: utf-8 -*-
import cv2
import numpy as np
import matplotlib.pyplot as plt

# 读取图像
img = cv2.imread('durant.jpg')
lenna_img = cv2.cvtColor(img, cv2.COLOR_BGR2RGB)

# 灰度化处理图像
grayImage = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)

# 高斯滤波
gaussianBlur = cv2.GaussianBlur(grayImage, (3, 3), 0)

# 阈值处理
ret, binary = cv2.threshold(gaussianBlur, 127, 255, cv2.THRESH_BINARY)

# Scharr算子
x = cv2.Scharr(grayImage, cv2.CV_32F, 1, 0)  # X方向
y = cv2.Scharr(grayImage, cv2.CV_32F, 0, 1)  # Y方向
absX = cv2.convertScaleAbs(x)
absY = cv2.convertScaleAbs(y)
Scharr = cv2.addWeighted(absX, 0.5, absY, 0.5, 0)

# Canny算子
gaussian = cv2.GaussianBlur(grayImage, (3, 3), 0)  # 高斯滤波降噪
Canny = cv2.Canny(gaussian, 50, 150)

# LOG算子
gaussian = cv2.GaussianBlur(grayImage, (3, 3), 0)  # 先通过高斯滤波降噪
dst = cv2.Laplacian(gaussian, cv2.CV_16S, ksize=3)  # 再通过拉普拉斯算子做边缘检测
LOG = cv2.convertScaleAbs(dst)

# 效果图
titles = ['Source Image', 'Gray Image', 'Binary Image',
          'Scharr Image', 'Canny Image', 'LOG Image']
images = [lenna_img, grayImage, binary, Scharr, Canny, LOG]
for i in np.arange(6):
    plt.subplot(2, 3, i + 1), plt.imshow(images[i], 'gray')
    plt.title(titles[i])
    plt.xticks([]), plt.yticks([])
plt.show()

 　　结果如下：

 4.3 边缘检测示例3

　　下面是采用高斯滤波去噪和阈值化处理之后，再进行边缘检测的过程，并对比了四种常见的边缘提取算法。

　　代码如下：

# -*- coding: utf-8 -*-
import cv2
import numpy as np
import matplotlib.pyplot as plt

# 读取图像
img = cv2.imread('durant.jpg')
lenna_img = cv2.cvtColor(img, cv2.COLOR_BGR2RGB)

# 灰度化处理图像
grayImage = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)

# 高斯滤波
gaussianBlur = cv2.GaussianBlur(grayImage, (3, 3), 0)

# 阈值处理
ret, binary = cv2.threshold(gaussianBlur, 127, 255, cv2.THRESH_BINARY)

# Roberts算子
kernelx = np.array([[-1, 0], [0, 1]], dtype=int)
kernely = np.array([[0, -1], [1, 0]], dtype=int)
x = cv2.filter2D(binary, cv2.CV_16S, kernelx)
y = cv2.filter2D(binary, cv2.CV_16S, kernely)
absX = cv2.convertScaleAbs(x)
absY = cv2.convertScaleAbs(y)
Roberts = cv2.addWeighted(absX, 0.5, absY, 0.5, 0)

# Prewitt算子
kernelx = np.array([[1, 1, 1], [0, 0, 0], [-1, -1, -1]], dtype=int)
kernely = np.array([[-1, 0, 1], [-1, 0, 1], [-1, 0, 1]], dtype=int)
x = cv2.filter2D(binary, cv2.CV_16S, kernelx)
y = cv2.filter2D(binary, cv2.CV_16S, kernely)
absX = cv2.convertScaleAbs(x)
absY = cv2.convertScaleAbs(y)
Prewitt = cv2.addWeighted(absX, 0.5, absY, 0.5, 0)

# Sobel算子
x = cv2.Sobel(binary, cv2.CV_16S, 1, 0)
y = cv2.Sobel(binary, cv2.CV_16S, 0, 1)
absX = cv2.convertScaleAbs(x)
absY = cv2.convertScaleAbs(y)
Sobel = cv2.addWeighted(absX, 0.5, absY, 0.5, 0)

# 拉普拉斯算法
dst = cv2.Laplacian(binary, cv2.CV_16S, ksize=3)
Laplacian = cv2.convertScaleAbs(dst)

# 效果图
titles = ['Source Image', 'Binary Image', 'Roberts Image',
          'Prewitt Image', 'Sobel Image', 'Laplacian Image']
images = [lenna_img, binary, Roberts, Prewitt, Sobel, Laplacian]
for i in np.arange(6):
    plt.subplot(2, 3, i + 1), plt.imshow(images[i], 'gray')
    plt.title(titles[i])
    plt.xticks([]), plt.yticks([])
plt.show()

 　　结果如下：

 

 　　于2020.10.21 再度学习，并修改内容。

参考：https://www.cnblogs.com/cfantaisie/archive/2011/06/05/2073151.html

https://blog.csdn.net/sinat_32974931/article/details/51125516

https://www.cnblogs.com/lynsyklate/p/7881300.html

https://blog.csdn.net/Eastmount/article/details/89001702

https://blog.csdn.net/wsp_1138886114/article/details/82935839

https://blog.csdn.net/JimmyFu0055/article/details/83719438

https://blog.csdn.net/Eastmount/article/details/89056240