python实现GPU编程——边缘检测（Prewitt 算子和Sobel算子）_边缘检测模型

1.图像滤波

      “滤波”（平滑）通常又称“模糊”，是一种简单常用的图像处理操作。进行平滑处理的原因有很多，但通常是用来去除噪声和相机失真，平滑在按照一定的原理来降低图像分辨率中也有重要应用。

滤波可理解为滤波器(通常为3×3、5×5矩阵)在图像上进行从上到下（如图1所示），从左到右的遍历，计算滤波器与对应像素的值并根据滤波目的进行数值计算返回值到当前像素点，滤波器对图像进行点积运算并赋值到图像，计算过程如公式（1）所示，其中i、j表示滤波器的坐标，表示第i，j个位子所对应的滤波器的值，表示滤波位子所对应的图像像素值。

                                               公式（1）

若滤波器值全为1，则计算的是当前像素点的平均值，因此滤波完后的图像应表现为模糊的效果，模糊程度取决于滤波器的尺寸(size)，滤波器尺寸越大，模糊效果越明显。

                                                       图1 图像滤波

       根据空间滤波增强目的可分为：平滑滤波和锐化滤波。平滑滤波，能减弱或消除图像中的高频分量，但不影响低频分量。因为高频分量对应图像中的区域边缘等灰度值具有较大、较快变化的部分，平滑滤波将这些分量绿区可减少局部灰度的起伏，使图像变得比较平滑。实际应用中，平滑滤波即可以用来消除噪声，又可以用在提取较大的目标前过滤去除较小的细节或将目标内的小间断连接起来。锐化滤波，能减弱或消除图像中的低频分量，但不影响高频分量。因为低频分量对应图像中灰度值缓慢变化的区域，因而与图像的整体特性如整体对比度和平均灰度值等有关。锐化滤波将这些分量滤去可使图像反差增加，边缘明显。实际应用中，锐化滤波可用于增强图像中被模糊的细节或景物的边缘。

       根据空间滤波的特点可分为：线性滤波和非线性滤波。线性滤波：运算只是对各像素灰度值进行简单处理（如乘一个权值）最后求和。非线性滤波：如果对像素灰度值的复杂运算，而不是最后求和的简单运算。

2.边缘检测

       图像的边缘检测是图像处理的一个基本步骤，是图像处理中的一个基本的研究方向和板块。它的主要原理在于识别出数字图像中那些颜色变化或者亮度变化明显的像素点，这些像素点的显著性变化往往代表图像的这部分属性发生了重要变化，其中包括了深度上的不连续、方向上的不连续及亮度上的不连续等。

       边缘检测算法在对图像的边缘进行检测时，先大概检测出图像轮廓的一些像素电，然后通过一些连接规则将那些像素点连接起来，最后再检测并连接一些之前未被识别的边界点、去除检测到的虚假的像素点和边界点并形成一个整体的边缘。然而在实际的图像中，边缘往往是各种类型的东西或模糊的风景的边缘，同时实际图像中可能存在着噪声，噪声和边缘同属于高频率的信号信息，因此传统中使用的频带过滤方法去检测图像的边缘效果并不好。

       边缘检测本质上就是一种滤波算法，区别在于滤波器的选择，滤波的规则是完全一致的。目前常用的边缘检测模型有很多：一阶的有Roberts算子，Prewitt算子，Sobel算子，Canny算子等；二阶的有Laplacian算子等。图像的边缘检测是基于图像的梯度来实现的，而获得图像的梯度就转化成使用各种算子对图像进行卷积运算来获得的。因此图像的边缘检测算法的核心在于算子。下面分别介绍三个一阶梯度算子。

2.1. Roberts算子

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

                      公式（2）

如公式（2）所示，Roberts算子的模板分为水平方向和垂直方向，从其模板可以看出，Roberts算子能较好的增强正负45度的图像边缘。

Roberts算子常用于垂直边缘明显或具有陡峭的低噪声的图像的边缘检测任务。

2.2. Prewitt 算子

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

           公式（3）

如公式（3）所示，Prewitt算子适合用来识别噪声较多、灰度渐变的图像。

prewitt算子常用于噪声较多、灰度渐变的图像的边缘检测任务。

2.3. Sobel 算子

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

      公式（4）

如公式（4）所示，Sobel算子的边缘定位更准确，常用于噪声较多、灰度渐变的图像。

sobel算子常用于噪声较多，灰度渐变的图像的边缘检测任务，效果更好。

3. 实验

3.1. 图像读取及算子设计代码

# coding=utf-8
import numpy as np
import cv2
from numba import cuda
import math
import os
import time
if __name__ == '__main__':
    # 读取影像并拷贝到device
    img = cv2.imread("D:/Codes/one_two_else/GPU/image1.jpg", cv2.IMREAD_GRAYSCALE)

    # 构造卷积核
    kernel1 = np.array([[-1, 0, 1], [-1, 0, 1], [-1, 0, 1]])
    kernel2 = np.array([[1, 1, 1], [0, 0, 0], [-1, -1, -1]])

#     kernel1 = np.array([[-1, 0, 1], [-2, 0, 2], [-1, 0, 1]])
#     kernel2 = np.array([[-1, -2, -1], [0, 0, 0], [1, 2, 1]])

    # 新建device变量用于存放结果
    result1 = np.zeros((img.shape[0], img.shape[1]))
    result2 = np.zeros((img.shape[0], img.shape[1]))
    start_time = time.time()

    # 执行卷积
    filterImgCPU(img, kernel1, result1)
    filterImgCPU(img, kernel2, result2)
    end_time = time.time()
    print('cost time', end_time - start_time,'s')

    # 直接得到的结果是float类型，而且也不是0到255，所以先将数据转换到合适范围，再转成uint8
    result1_res = np.where(result1 < 0, 0, result1)
    result1 = np.uint8(np.where(result1_res > 255, 255, result1_res))
    result2_res = np.where(result2 < 0, 0, result2)
    result2 = np.uint8(np.where(result2_res > 255, 255, result2_res))
    # 保存各方向的轮廓
 os.makedirs("D:/Codes/one_two_else/GPU/Edgeimg/image1cpu/Prewitt",exist_ok=True)
    cv2.imwrite("D:/Codes/one_two_else/GPU/Edgeimg/image1cpu/Prewitt/contour1.jpg", result1)
cv2.imwrite("D:/Codes/one_two_else/GPU/Edgeimg/image1cpu/Prewitt/contour2.jpg", result2)
    # 将各方向轮廓合并
tmp1 = cv2.bitwise_or(result1, result2)
cv2.imwrite("D:/Codes/one_two_else/GPU/Edgeimg/image1cpu/Prewitt/contour.jpg", tmp1)

3.2. CPU函数实现

def filterImgCPU(input_img, filter, filtered_img):

    filter_h = filter.shape[0]
    filter_w = filter.shape[1]
    half_filter_h = filter_h / 2
    half_filter_w = filter_w / 2
    for x in range(filtered_img.shape[0]):
        for y in range(filtered_img.shape[1]):
            if half_filter_h < x < filtered_img.shape[0] - half_filter_h and half_filter_w < y < filtered_img.shape[1] - half_filter_w:
                patch = input_img[int(x - half_filter_h):int(x + half_filter_h + 1), int(y - half_filter_w):int(y + half_filter_w + 1)]
                tmp = 0.0
                for i in range(filter_h):
                    for j in range(filter_w):
                        tmp += patch[j, i] * filter[j, i]
                filtered_img[x, y] = tmp
            else:
                filtered_img[x, y] = input_img[x, y]

3.3. GPU函数实现

    from numba import cuda

@cuda.jit
def filterImgGPU(input_img, filter, filtered_img):
    filter_h = filter.shape[0]
    filter_w = filter.shape[1]
    half_filter_h = filter_h / 2
    half_filter_w = filter_w / 2

    # x, y = cuda.grid(2)
    x = cuda.blockIdx.x * cuda.blockDim.x + cuda.threadIdx.x
    y = cuda.blockIdx.y * cuda.blockDim.y + cuda.threadIdx.y
    if x < filtered_img.shape[0] and y < filtered_img.shape[1]:
        if half_filter_h < x < filtered_img.shape[0] - half_filter_h and \
                half_filter_w < y < filtered_img.shape[1] - half_filter_w:
            patch = input_img[x - half_filter_h:x + half_filter_h + 1, y - half_filter_w:y + half_filter_w + 1]
            tmp = 0.0
            for i in range(filter.shape[0]):
                for j in range(filter.shape[1]):
                    tmp += patch[j, i] * filter[j, i]
            filtered_img[x, y] = tmp
        else:
            filtered_img[x, y] = input_img[x, y]

3.4. 结果比较

                                                          图2 原始图像

3.4.1. Prewitt 算子边缘检测结果 

如图3所示为GPU实现的Prewitt 算子进行边缘检测所用时间为0.30s，图4所示为CPU实现的Prewitt 算子进行边缘检测所用时间为13.02s。

                                                      图3 Prewitt 算子GPU运行时间

                                                图4 Prewitt 算子CPU运行时间

                              图5  Prewitt算子沿x方向（左）、y方向（右）提取的边缘

                                                          图6 Prewitt算子边缘检测结果

3.4.2. Sobel 算子边缘检测结果

如图7所示为GPU实现的Sobel 算子进行边缘检测所用时间为0.29s，图8所示为CPU实现的Sobel 算子进行边缘检测所用时间为13.14s。

                                                   图7 Sobel 算子GPU运行时间

                                                  图8 Sobel 算子CPU运行时间

                               图9  Sobel 算子沿x方向（左）、y方向（右）提取的边缘

                                                图10 Sobel算子边缘检测结

4. 总结

        本次课程详细介绍了GPU编程，通过本次实验，我学习了如何调用python的numba库进行GPU加速数据运算，学习了如何设置grid、block、thread的大小以及如何遍历每一个线程。通过本次实验证明GPU编程可以在一定程度上提高运算速率，在实现的边缘提取方面，GPU运算时间仅为CPU运算时间的四分之一。