西电计科大数据计算机视觉作业一sobel算子和canny算子

基于python对Sobel和Canny算子的复现

Sobel算子部分

1. X方向梯度

   **图1.1Sobel x方向卷积核**

通过Sobel的x方向卷积核（如图1.1）与通过opencv读取到的灰度值图像矩阵进行乘法运算卷积运算得到新的图像。

通过循环实现卷积核与该图像的所有像素点都经过计算。

最后过滤掉一些像素值较小的噪音点，达到最终图像，代码如下所示：

kernel\_x = np.array([[-1, 0, 1], [-2, 0, 2], [-1, 0, 1]])  

# sobel x方向卷积核

# x轴方向
def sobel\_x(img, threshold):
      W = np.size(img, 0)
      H = np.size(img, 1)  # 计算图片的长度和宽度
      mag = np.zeros(img.shape)  # 创建一个图片形状的空矩阵
      for i in range(0, W - 2):
          for j in range(0, H - 2):
              v = sum(sum(kernel\_x \* img[i:i + 3, j:j + 3]))  # 进行矩阵卷积运算
              mag[i + 1, j + 1] = v

      for p in range(0, W):  # 过滤掉一些噪音点，让主体突出
          for q in range(0, H):
              if mag[p, q] < threshold:
                  mag[p, q] = 0
      return mag

1. Y方向梯度

Y方向梯度算法与x方向相似，只需将x方向梯度中的x卷积核替换成y方向的卷积核(如图2.1)再进行卷积运算即可。

 **图2.1Sobel y方向卷积核**

其代码如下：

kernel\_y = np.array([[-1, -2, -1], [0, 0, 0], [1, 2, 1]])  # sobel y方向卷积核

# y轴方向
def sobel\_y(img, threshold):
      W = np.size(img, 0)
      H = np.size(img, 1)
      mag = np.zeros(img.shape)
      for i in range(0, W - 2):
          for j in range(0, H - 2):
              h = sum(sum(kernel\_y \* img[i:i + 3, j:j + 3]))
              mag[i + 1, j + 1] = h

      for p in range(0, W):
          for q in range(0, H):
              if mag[p, q] < threshold:
                  mag[p, q] = 0
      return mag

1. 滤波结果

综合前两部的x方向梯度与y方向梯度，我们在算出一个像素点的x方向梯度与y方向梯度之后，通过平方和后取根号的形式，来获得该点的像素值。

kernel\_x = np.array([[-1, 0, 1], [-2, 0, 2], [-1, 0, 1]])  # sobel x方向卷积核
kernel\_y = np.array([[-1, -2, -1], [0, 0, 0], [1, 2, 1]])  # sobel y方向卷积核
def sobel1(img, threshold):
      W = np.size(img, 0)
      H = np.size(img, 1)
      mag = np.zeros(img.shape)
      for i in range(0, W - 2):
          for j in range(0, H - 2):
              v = sum(sum(kernel\_x \* img[i:i + 3, j:j + 3]))  # x轴方向
              h = sum(sum(kernel\_y \* img[i:i + 3, j:j + 3]))  # y轴方向
              z = np.sqrt((v \*\* 2) + (h \*\* 2))  
              if z < threshold:  # 如果小于阈值 则为0
                  mag[i + 1, j + 1] = 0
              else:
                  mag[i + 1, j + 1] = z
      return mag

2. 梯度幅值

		 **图4.1计算幅值时Sobel x方向与y方向卷积核**

计算真实的幅值时我们需要用到计算幅值的卷积核（如图4.1）这与Sobel算子的标准定义不同。其计算方式与上一节提到的计算方式相似。

其代码如下：

kernel\_x = np.array([[-1, 0, 1], [-2, 0, 2], [-1, 0, 1]])  # sobel x方向卷积核
kernel\_y = np.array([[-1, -2, -1], [0, 0, 0], [1, 2, 1]])  # sobel y方向卷积核

a = 1 / 8 \* kernel\_x  # 算幅值时的卷积核
b = 1 / 8 \* kernel\_y

def sobel\_amplitude(img, threshold):
      W = np.size(img, 0)
      H = np.size(img, 1)
      mag = np.zeros(img.shape)
      for i in range(0, W - 2):
          for j in range(0, H - 2):
              v = sum(sum(a \* img[i:i + 3, j:j + 3]))  # x轴方向
              h = sum(sum(b \* img[i:i + 3, j:j + 3]))  # y轴方向
              z = np.sqrt((v \*\* 2) + (h \*\* 2))  # 计算幅值
              if z < threshold:  # 如果幅值小于阈值 则为0
                  mag[i + 1, j + 1] = 0
              else:
                  mag[i + 1, j + 1] = z
      return mag

1. 计算梯度角度

在我们分别计算出gy与gx后，通过求取他们的arctan值来计算出梯度角度。

Python math库中的atan2函数已经考虑到gx等于0的情况所以我们进行分类讨论。

atan2函数的返回值为（-Π，Π），我们通过映射和计算将他转化为（0，360）。

再将这些不同数值的像素通过色彩来进行区分，这里调用了matplotlib的cmap函数，中的rainbow色块如下图所示

							  ** 图5.1rainbow色块展示 **

因cmap rainbow色块支持256个可选值 所以我们再将（0，360）的角度映射到

（0，256）上。如此角度越小的点颜色越偏蓝色，角度越大的点颜色偏红色

其代码如下：

def sobel\_amplitude(img, threshold):
    W = np.size(img, 0)
    H = np.size(img, 1)
    mag = np.zeros(img.shape)
    for i in range(0, W - 2):
        for j in range(0, H - 2):
            v = sum(sum(a \* img[i:i + 3, j:j + 3]))  # x轴方向
            h = sum(sum(b \* img[i:i + 3, j:j + 3]))  # y轴方向
            z = np.sqrt((v \*\* 2) + (h \*\* 2))  # 计算幅值
            if z < threshold:  # 如果幅值小于阈值 则为0
                mag[i + 1, j + 1] = 0
            else:
                # 如果赋值大于阈值 则计算他的角度
                #  这里用了math库的artan2的函数其返回值为（-Π，Π） 通过算数运算将其转化为（0，360）
                #  因matplotlib的cmap库色彩可选值为256个 再经过运算转化为（0，256）
                mag[i + 1, j + 1] = z
    return mag

    plt.imshow(mag\_angle, plt.get\_cmap('rainbow'))

1. Sobel结果展示

6.1Sobel算子各图结果展示

绘图代码如下：

def sobel(image):
    # image = cv2.imread('2.jpg', 0)  # read an image
    mag\_y = sobel\_y(image, 5)
    mag\_x = sobel\_x(image, 5)
    mag\_amplitude = sobel\_amplitude(image, 5)
    mag\_angle = sobel\_angle(image, 50)
    mag\_sobel = sobel1(image, 5)
    plt.figure("Sobel", frameon=False)  # 图像窗口名称
    plt.subplot(2, 3, 5)
    plt.imshow(mag\_x, cmap='gray')
    plt.title("x方向图", fontsize=8)
    plt.xticks([])
    plt.yticks([])
    plt.subplot(2, 3, 6)
    plt.imshow(mag\_y, cmap='gray')
    plt.title("y方向", fontsize=8)
    plt.xticks([])
    plt.yticks([])
    plt.subplot(2, 3, 3)
    #plt.imshow(mag\_angle)
    plt.imshow(mag\_angle, plt.get\_cmap('rainbow'))
    plt.title("角度图", fontsize=8)
    plt.xticks([])
    plt.yticks([])
    plt.subplot(2, 3, 4)
    plt.imshow(mag\_amplitude, cmap='gray')
    plt.title("幅值图", fontsize=8)
    plt.xticks([])
    plt.yticks([])
    plt.subplot(2, 3, 2)
    plt.imshow(mag\_sobel, cmap='gray')
    plt.title("结果图", fontsize=8)
    plt.xticks([])
    plt.yticks([])
    plt.subplot(2, 3, 1)
    plt.imshow(image, cmap='gray')
    plt.title("原图", fontsize=8)
    plt.xticks([])
    plt.yticks([])

7. Sobel算子结果分析

7.1. x方向梯度与y方向梯度对比发现x方向梯度图片在竖直方向和水平方向分别有所空缺

7.2 结果图与幅值图对比

虽然两张图片宏观上看十分相似没有什么区别，我一开始也以为自己是否做错。

但是当我放大两张图片时，发现在微观上，如果卷积核乘1/8后，其边缘将会相比于结果图更加清晰。我们将图片局部放大并将亮度调高后会很明显的看出，在乘1/8后图片跟接近于真实的幅值。

Canny算子部分

1. 高斯滤波

高斯滤波器（kernel）是将高斯函数离散化，将滤波器中对应的横纵坐标索引代入高 斯函数，即可得到对应的值。

(2k+1)x(2k+1) 滤波器的计算公式如右：

常见的高斯滤波器为size=5，其近似值为：

我们依旧用矩阵运算将待测图片进行高斯模糊。

其代码如下：

def smooth(img, sigma=1.4, length=5):
      # 生成高斯核
      k = length // 2
      gaussian = np.zeros([length, length])
      for i in range(length):
          for j in range(length):
              gaussian[i, j] = np.exp(-((i - k) \*\* 2 + (j - k) \*\* 2) / (2 \* sigma \*\* 2))
      gaussian /= 2 \* np.pi \* sigma \*\* 2
      gaussian = gaussian / np.sum(gaussian)

      # 用高斯核进行滤波
      W = np.size(img, 0)
      H = np.size(img, 1)
      new\_image = np.pad(img, ((1, 1), (1, 1)), constant\_values=0)

      for i in range(W - 2 \* k):
          for j in range(H - 2 \* k):
              new\_image[i, j] = np.sum(img[i:i + 5, j:j + 5] \* gaussian)

      return new\_image

1. 计算图片的幅值与角度

其计算步骤与Sobel计算方式相同，此处不再赘述。

其代码如下：

def getGradAngle(image):  # 用sobel核计算图片的幅值和梯度角度
      *""" 
           -1 0 1        -1 -2 -1
      Gx = -2 0 2   Gy =  0  0  0
           -1 0 1         1  2  1
      """*
      Gx = np.array([[-1, 0, 1], [-2, 0, 2], [-1, 0, 1]])
      Gy = np.array([[-1, -2, -1], [0, 0, 0], [1, 2, 1]])

      W = np.size(image, 0)
      H = np.size(image, 1)
      amplitude = np.zeros([W - 2, H - 2])  # 幅值数组
      angle = np.zeros([W - 2, H - 2])  # 角度数组

      for i in range(W - 2):
          for j in range(H - 2):
              dx = np.sum(image[i:i + 3, j:j + 3] \* Gx)
              dy = np.sum(image[i:i + 3, j:j + 3] \* Gy)
              amplitude[i, j] = np.sqrt(dx \*\* 2 + dy \*\* 2)
              angle[i, j] = math.atan2(dy, dx)

      return amplitude, angle

1. 非最大值抑制（NMS）

当我们计算一点C时我们会找到他梯度方向的相邻点dTmp1与dTmp2。如果C点不是这三个点中的最大值时，我们则将C点的像素值置0。

计算dTmp1时我们可能会遇到这两个点并不能直接被获取，这时我们用类似线性插值的方式，用它临近点g1与g2共同来描述该点的像素值。具体的权重通过theta角度算tan值来描述。这里涉及一些数学中角度换算的过程，详细信息如代码所示：

def NMS(amplitude, angle):
      *""" Non-maxima suppression
          非最大值抑制
          遍历梯度方向两个其他节点
          如果有值比本身大，则将本身置为0
      """*
      W = np.size(amplitude, 0)
      H = np.size(angle, 1)
      nms = amplitude.copy()
      # 当梯度不为45的整数倍时 通过同行相邻节点加权算出该点的值
      # 通过角度来计算权重
      for i in range(1, W - 1):
          for j in range(1, H - 1):
              theta = angle[i, j]
              weight = np.tan(theta)
              # 不同角度的权重不同
              if theta > np.pi / 4:
                  d1 = [0, 1]
                  d2 = [1, 1]
                  weight = 1 / weight
              elif theta >= 0:
                  d1 = [1, 0]
                  d2 = [1, 1]
              elif theta >= - np.pi / 4:
                  d1 = [1, 0]
                  d2 = [1, -1]
                  weight \*= -1
              else:
                  d1 = [0, -1]
                  d2 = [1, -1]
                  weight = -1 / weight

              g1 = amplitude[i + d1[0], j + d1[1]]
              g2 = amplitude[i + d2[0], j + d2[1]]
              g3 = amplitude[i - d1[0], j - d1[1]]
              g4 = amplitude[i - d2[0], j - d2[1]]

              grade\_count1 = g1 \* weight + g2 \* (1 - weight)
              grade\_count2 = g3 \* weight + g4 \* (1 - weight)

              if grade\_count1 > amplitude[i, j] or grade\_count2 > amplitude[i, j]:
                  nms[i, j] = 0

      return nms

1. 双阈值抑制

设有阈值T1<T2

这里通过dfs算法扫描所有的大于T2的点，并扫描这些点的边缘中是否有大于T1的点。将所有小于T1的点删除，将大于T2的点和大于T1且与大于T2的点相连的点保留。

注：图片中一个像素点的斜上、下方都是该点的相邻点，所以遍历它相邻点时，方向数组应该有八个值。

 其代码如下：

def double\_threshold(nms, threshold1, threshold2):
      *""" Double Threshold
      通过dfs找出所有强像素点的所有联通点
      """*
      visited = np.zeros\_like(nms)
      output\_image = nms.copy()
      W, H = output\_image.shape

      def dfs(i, j):
          *"""
          当像素值超过第二阈值时直接保留
          像素值小于第一阈值时直接删除
          与第二阈值相连的且大于第一阈值的像素保留
          通过dfs来查询这些与第二阈值相连的大于第一阈值的点
          """*
          #方向数组
          dx = [-1, -1, -1, 0, 0, 1, 1, 1]
          dy = [-1, 0, 1, -1, 1, -1, 0, 1]
          if i >= W or i < 0 or j >= H or j < 0 or visited[i, j] == 1:
              return
          visited[i, j] = 1
          if output\_image[i, j] > threshold1:
              output\_image[i, j] = 255
              for p in range(0, 8):
                  dfs(i + dx[p], j + dy[p])
          else:
              output\_image[i, j] = 0

      for w in range(W):
          for h in range(H):
              if visited[w, h] == 1:
                  continue
              if output\_image[w, h] >= threshold2:
                  dfs(w, h)
              elif output\_image[w, h] <= threshold1:
                  output\_image[w, h] = 0
                  visited[w, h] = 1
      # 将剩余其他不连通的点置为0
      for w in range(W):
          for h in range(H):
              if visited[w, h] == 0:
                  output\_image[w, h] = 0

      return output\_image

1. Canny结果展示

以下三张图片分别为原图、模糊图、与结果图

展示代码如下所示：

def canny(image):
      smoothed\_image = smooth(image)
      amplitude, angle = getGradAngle(smoothed\_image)
      nms = NMS(amplitude, angle)
      output\_image = double\_threshold(nms, 10, 60)
      plt.figure("Canny", frameon=False)
      plt.subplot(1, 3, 1)
      plt.imshow(image)
      plt.title("原图", fontsize=8)
      plt.xticks([])
      plt.yticks([])
      plt.subplot(1, 3, 2)
      plt.imshow(smoothed\_image)
      plt.title("高斯模糊图", fontsize=8)
      plt.xticks([])
      plt.yticks([])
      plt.subplot(1, 3, 3)
      plt.imshow(output\_image)
      plt.title("滤波图", fontsize=8)
      plt.xticks([])
      plt.yticks([])
      # plt.show()

6．Canny结果分析

Canny算子的双阈值抑制算法与最终结果的数值息息相关，要想得到一个非常好的边缘检测结果，那这两个值可能需要精心选择。其结果较为精准，能将轮廓与背景很好的区分。

Sobel与Canny算子对比

相比来说sobel算子处理图片速度更快，但Canny算子更经典，精准度更高，能更好的去除噪音并保留更清晰的线条。Canny边缘检测之所以优秀是因为它在一阶微分算子的基础上，增加了非最大值抑制和双阈值两项改进。利用非极大值抑制不仅可以有效地抑制多响应边缘，而且还可以提高边缘的定位精度；利用双阈值可以有效减少边缘的漏检率。但Canny算子设计更多参数，需要更准确地设置参数，更繁杂一些。

代码地址:

https://gitee.com/orangeinus/xd_-cs_computer_vison_1.git