几种边缘检测算子的python实现

发表于 | 更新于: | 分类于 | 阅读次数: |
字数统计: 2,116字 | 阅读时长 ≈ 12分钟

首先导入必备的库和图像，展示原始图像，如下所示:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
from PIL import Image
import numpy as np
from matplotlib import pyplot as plt

image_ori1 = Image.open("test/circuit.tif")
image_ori1 = np.array(image_ori1)

image_ori2 = Image.open("test/lena.tif")
image_ori2 = np.array(image_ori2)

image_ori3 = Image.open("test/poly.bmp")
image_ori3 = np.array(image_ori3)

plt.figure(figsize=(10,10))
plt.subplot(131)
plt.imshow(image_ori1,cmap ='gray')
plt.subplot(132)
plt.imshow(image_ori2,cmap ='gray')
plt.subplot(133)
plt.imshow(image_ori3,cmap ='gray')
plt.show()

1. Sobel算子

定义sobel算子的核，如下所示：

1
2
3
4
5
def Sobel_Kernel():
    sobel_x = np.array(([-1, 0, 1], [-2, 0, 2], [-1, 0, 1]))
    sobel_y = np.array(([-1, -2, -1], [0, 0, 0], [1, 2, 1]))
    sobel = np.array(([-1, -1, 0], [-1, 0, 1], [0, 1, 1]))
    return sobel,sobel_x,sobel_y

然后定义卷积操作，如下所示：

1
2
3
4
5
6
7
8
def process(img,kernel):
    img_new = np.zeros(shape=(img.shape[0]-kernel.shape[0]//2,img.shape[1]-kernel.shape[0]//2))
    for x in range(0,img.shape[0]-kernel.shape[0]+1):
        for y in range(0,img.shape[1]-kernel.shape[0]+1):
                img_new[x,y] = np.sum(img[x:x+kernel.shape[0],y:y+kernel.shape[0]]*kernel)


    return img_new

然后处理图像，如下所示，x方向，y方向，和综合的sobel边缘提取图如下所示，下面为sobel_x提取y方向的边缘，

1
2
3
4
5
6
7
8
9
10
11
12
13
sobel,sobel_x,sobel_y = Sobel_Kernel()
image1_x = process(image_ori1,kernel=sobel_x)
image2_x = process(image_ori2,kernel=sobel_x)
image3_x = process(image_ori3,kernel=sobel_x)

plt.figure(figsize=(10,10))
plt.subplot(131)
plt.imshow(image1_x,cmap ='gray')
plt.subplot(132)
plt.imshow(image2_x,cmap ='gray')
plt.subplot(133)
plt.imshow(image3_x,cmap ='gray')
plt.show()

x方向边缘如下：

1
2
3
4
5
6
7
8
9
10
11
12
image1_y = process(image_ori1,kernel=sobel_y)
image2_y = process(image_ori2,kernel=sobel_y)
image3_y = process(image_ori3,kernel=sobel_y)

plt.figure(figsize=(10,10))
plt.subplot(131)
plt.imshow(image1_y,cmap ='gray')
plt.subplot(132)
plt.imshow(image2_y,cmap ='gray')
plt.subplot(133)
plt.imshow(image3_y,cmap ='gray')
plt.show()

综合方向如下：

1
2
3
4
5
6
7
8
9
10
11
12
image1_z = process(image_ori1,kernel=sobel)
image2_z = process(image_ori2,kernel=sobel)
image3_z = process(image_ori3,kernel=sobel)

plt.figure(figsize=(10,10))
plt.subplot(131)
plt.imshow(image1_z,cmap ='gray')
plt.subplot(132)
plt.imshow(image2_z,cmap ='gray')
plt.subplot(133)
plt.imshow(image3_z,cmap ='gray')
plt.show()

可以看出来效果还可以。

马尔算子

首先定义一个大小为5的高斯核，如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
def Gaussian_Kernel5(kernel_size=5, sigma=1.0):
    kernel = np.zeros([kernel_size, kernel_size])
    center = kernel_size // 2

    s = 2 * (sigma ** 2)
    sum_val = 0
    for i in range(0, kernel_size):
        for j in range(0, kernel_size):
            x = i - center
            y = j - center
            kernel[i, j] = np.exp(-(x ** 2 + y ** 2) / s)
            sum_val += kernel[i, j]

    sum_val = 1 / sum_val
    return kernel * sum_val

然后定义一个laplacian高斯核，如下所示：

1
2
3
def Laplacian_Kernel():
    laplacain = np.array(([-1, -1, -1], [-1, 8, -1], [-1, -1, -1]))
    return laplacain

判断是否过0点的函数如下，注意很难找到一个正好二阶导为0的像素点，所以这里判断相邻两个像素的成绩是否小于一个足够大的负值，若满足条件，说明其是一个边缘点（梯度足够大且过0），则标记为白色点（255），否则为黑色。

1
2
3
4
5
6
7
8
9
10
def If_Zero(img):
    img_new = np.zeros(shape=img.shape)
    for x in range(1, img.shape[0]-1):
        for y in range(1, img.shape[1]-1):
                if (img[x,y]*img[x+1,y] < -144 or img[x,y]*img[x-1,y] < -144) or (img[x,y]*img[x,y+1] < -144  or img[x,y]*img[x,y-1] < -144):
                    img_new[x,y] = 255
                else:
                    img_new[x,y] = 0

    return img_new

然后是显示最后的图像结果，如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
kerner5_g = Gaussian_Kernel5()
laplacian_k = Laplacian_Kernel()
image1_m = If_Zero(process(process(image_ori1,kerner5_g),laplacian_k))

kerner5_g = Gaussian_Kernel5()
laplacian_k = Laplacian_Kernel()
image2_m = If_Zero(process(process(image_ori2,kerner5_g),laplacian_k))

kerner5_g = Gaussian_Kernel5()
laplacian_k = Laplacian_Kernel()
image3_m = If_Zero(process(process(image_ori3,kerner5_g),laplacian_k))


plt.figure(figsize=(10,10))
plt.subplot(131)
plt.imshow(image1_m,cmap ='gray')
plt.subplot(132)
plt.imshow(image2_m,cmap ='gray')
plt.subplot(133)
plt.imshow(image3_m,cmap ='gray')
plt.show()

可以看出来检测出来的效果不咋地，可能是自己代码实现有问题QAQ。

Canny算子

首先定义一个大小为3的高斯卷积核，如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
def Gaussian_Kernel3(kernel_size=3, sigma=0.00001):
    kernel = np.zeros([kernel_size, kernel_size])
    center = kernel_size // 2

    s = 2 * (sigma ** 2)
    sum_val = 0
    for i in range(0, kernel_size):
        for j in range(0, kernel_size):
            x = i - center
            y = j - center
            kernel[i, j] = np.exp(-(x ** 2 + y ** 2) / s)
            sum_val += kernel[i, j]

    sum_val = 1 / sum_val
    return kernel * sum_val

然后是获得边缘梯度，这里直接使用PPT上获得边缘梯度的公式：

1
2
3
4
5
6
7
8
9
10
11
def Get_Gradient(img):
    img_dx = np.zeros(shape=img.shape)
    img_dy = np.zeros(shape=img.shape)
    img_d = np.zeros(shape=img.shape)
    for i in range(0,img.shape[0]-1):
        for j in range(0,img.shape[1]-1):
            img_dy[i, j] = (img[i+1, j] - img[i, j] + img[i+1, j+1] - img[i, j+1])/2.0
            img_dx[i, j] = (img[i, j+1] - img[i, j] + img[i+1, j+1] - img[i+1, j])/2.0
            img_d[i,j] = np.sqrt(np.square(img_dx[i,j])+np.square(img_dy[i,j]))

    return img_d,img_dx,img_dy

然后是进行非极大值抑制，这里对梯度方向两边的值和梯度值进行比较判断，如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
def Get_NMS(img_d,img_dx,img_dy):
    NMS = img_d
    NMS[0,:] = NMS[img_d.shape[0]-1,:] = NMS[:,0] = NMS[:,img_d.shape[1]-1] = 0
    for i in range(1,img_d.shape[0]-1):
        for j in range(1,img_d.shape[1]-1):
            d_temp1 = d_temp2 = 0
            if img_d[i,j] == 0:
                NMS[i,j] = img_d[i,j]
            else:

                if np.abs(img_dx[i,j])-np.abs(img_dy[i,j]) > 0:
                    d_weight = np.abs(img_dy[i,j]/img_dx[i,j])
                    if img_dy[i,j]*img_dx[i,j]>0:
                        d_temp1 = (1 - d_weight)*img_d[i+1,j+1] + d_weight*img_d[i+1,j]
                        d_temp2 = (1 - d_weight)*img_d[i-1,j-1] + d_weight*img_d[i-1,j]
                    else:
                        d_temp1 = (1 - d_weight) * img_d[i + 1, j - 1] + d_weight * img_d[i + 1, j]
                        d_temp2 = (1 - d_weight) * img_d[i - 1, j + 1] + d_weight * img_d[i - 1, j]
                elif np.abs(img_dx[i,j])-np.abs(img_dy[i,j]) < 0:
                    d_weight = np.abs(img_dx[i, j] / img_dy[i, j])
                    if img_dy[i,j] * img_dx[i,j] > 0:
                        d_temp1 = (1 - d_weight) * img_d[i + 1, j + 1] + d_weight * img_d[i, j+1]
                        d_temp2 = (1 - d_weight) * img_d[i - 1, j - 1] + d_weight * img_d[i, j-1]
                    else:
                        d_temp1 = (1 - d_weight) * img_d[i - 1, j + 1] + d_weight * img_d[i, j + 1]
                        d_temp2 = (1 - d_weight) * img_d[i + 1, j - 1] + d_weight * img_d[i, j - 1]

            if img_d[i,j]>=d_temp2 and img_d[i,j]>=d_temp1:
                NMS[i,j] = img_d[i,j]
            else:
                NMS[i,j] = 0
        return NMS

然后是进行连通性判断，看在阈值之间的区域的八邻域有没有值得连通的点：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
def Connect(img):
    TH = 0.18
    TL = 0.08
    Max = np.max(img)
    Res = img
    for i in range(img.shape[0]-1):
        for j in range(img.shape[1]-1):
            if img[i,j] > TH*Max:
                Res[i,j] = 255
            elif img[i,j] <TL*Max:
                Res[i,j] = 0
            elif (img[i,j+1] > TH*Max or img[i+1,j] > TH*Max or img[i,j-1] > TH*Max or img[i-1,j] > TH*Max or
                img[i-1, j+1] > TH * Max or img[i-1,j-1] > TH*Max or img[i+1,j+1] > TH*Max or img[i+1,j-1] > TH*Max):
                Res[i,j] = 255
    return Res

最后是获得边缘检测的结果，如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
kernel3_g = Gaussian_Kernel3()
img1_d,img1_dx,img1_dy = Get_Gradient(process(image_ori1,kernel3_g))
image1_c = Connect(Get_NMS(img1_d,img1_dx,img1_dy))

img2_d,img2_dx,img2_dy = Get_Gradient(process(image_ori2,kernel3_g))
image2_c = Connect(Get_NMS(img2_d,img2_dx,img2_dy))

img3_d,img3_dx,img3_dy = Get_Gradient(process(image_ori3,kernel3_g))
image3_c = Connect(Get_NMS(img3_d,img3_dx,img3_dy))

plt.figure(figsize=(10,10))
plt.subplot(131)
plt.imshow(image1_c,cmap ='gray')
plt.subplot(132)
plt.imshow(image2_c,cmap ='gray')
plt.subplot(133)
plt.imshow(image3_c,cmap ='gray')
plt.show()

效果一般般，感觉勉强可以接受。

Susan算子

这个地方我都感觉没用到核的定义，直接对USAN区域进行判断就行了，遍历37个值，大于阈值则判定为边缘，如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
def Get_Susan(img):
    T = 17
    new_img = np.zeros(shape=(img.shape[0] - 6,img.shape[1] - 6))
    for i in range(3,img.shape[0] - 3):
        for j in range(3, img.shape[1] - 3):
            num = 0
            for x in range(-1,2):
                for y in [-3,3]:
                    if np.abs(int(img[i+x,j+y]) - int(img[i,j])) < T:
                        num += 1
            for x in range(-2,3):
                for y in [-2,2]:
                    if np.abs(int(img[i + x, j + y]) - int(img[i,j])) < T:
                        num += 1

            for x in range(-3,4):
                for y in range(-1,2):
                    if np.abs(int(img[i+x,j+y])-int(img[i,j]))<T:
                        num += 1

            if 0.75*37-num >= 0:
                new_img[i-3, j-3] = 0.75*37-num
            else:
                new_img[i-3, j-3] = 0

    return new_img

然后是使用Canny里面的NMS和连通性判断，获得最后边缘，结果如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
img1_d,img1_dx,img1_dy = Get_Gradient(Get_Susan(image_ori1))
image1_s = Connect(Get_NMS(img1_d,img1_dx,img1_dy))

img2_d,img2_dx,img2_dy = Get_Gradient(Get_Susan(image_ori2))
image2_s = Connect(Get_NMS(img2_d,img2_dx,img2_dy))

img3_d,img3_dx,img3_dy = Get_Gradient(Get_Susan(image_ori3))
image3_s = Connect(Get_NMS(img3_d,img3_dx,img3_dy))

plt.figure(figsize=(10,10))
plt.subplot(131)
plt.imshow(image1_s,cmap ='gray')
plt.subplot(132)
plt.imshow(image2_s,cmap ='gray')
plt.subplot(133)
plt.imshow(image3_s,cmap ='gray')
plt.show()

emmm,感觉更差了，还不如canny，应该是代码有问题吧…….

0%