图像处理 | 最常用的边缘检测详解与代码(Robert, Sober, Prewitt, Canny, Kirsch, Laplacian, LOG, DOG算子)

边缘检测

目的:找到图像中亮度变化剧烈的像素点构成的集合,即表现出来往往是图像的轮廓,即边缘。
传统边缘检测的步骤:①滤波去除噪声;②增强边缘的特征;③提取边缘,完成检测。

类别常见算子
一阶Roberts Cross算子,Prewitt算子,Sobel算子,Canny算子,罗盘算子,krisch算子
二阶Laplacian算子,LoG,DoG,Marr-Hildreth,在梯度方向的二阶导数过零点。

在这里插入图片描述

一阶


Roberts Cross 罗伯茨交叉算子

Roberts算子是一种斜向偏差分的梯度计算方法,梯度的大小代表边缘的强度,梯度的方向与边缘的走向垂直。它是2X2算子模板。2个卷积核形成了Roberts算子。图象中的每一个点都用这2个核做卷积。更详细计算过程查看 博主苏源流

优点:计算简单,边缘定位较准,但对噪声极敏感。适用于边缘明显且噪声较少的图像分割。

Reboerts算子是一种利用局部差分来寻找边缘,采用对角方向相邻两像素值之差,算子计算形式如下:

G x = f ( i , j ) − f ( i − 1 , j − 1 ) ; G y = f ( i − 1 , j ) − f ( i , j − 1 ) Gx = f(i,j) - f(i-1,j-1); Gy = f(i-1,j) - f(i,j-1) Gx=f(i,j)f(i1,j1)Gy=f(i1,j)f(i,j1)
∣ G ( x , y ) ∣ = s p r t ( G x 2 − G y 2 ) |G(x,y)| = sprt(Gx^2-Gy^2) G(x,y)=sprt(Gx2Gy2)
即 G ( x , y ) = a b s ( f ( x , y ) − f ( x + 1 , y + 1 ) ) + a b s ( f ( x , y + 1 ) − f ( x + 1 , y ) ) 即 G(x,y)=abs(f(x,y)-f(x+1,y+1))+abs(f(x,y+1)-f(x+1,y)) G(x,y)=abs(f(x,y)f(x+1,y+1))+abs(f(x,y+1)f(x+1,y))
图片来自 残月飞雪,感谢

# python
# 边缘检测 Robertes cross算子=[[-1,-1],[1,1]] 图像增强锐化
# 一种斜向偏差分的梯度计算,梯度的大小代表边缘的强度,梯度的方向与边缘的走向垂直
# @Zivid 2021/8/16
import cv2
import numpy as np
from skimage import filters


# 方法一:自定义
def Robert01(pic):
    pic_f = np.copy(pic)
    pic_f = pic_f.astype("float")
    Roberts = np.zeros((row, column))
    for x in range(row - 1):
        for y in range(column - 1):
            gx = abs(pic_f[x + 1, y + 1] - pic_f[x, y])
            gy = abs(pic_f[x + 1, y] - pic_f[x, y + 1])
            Roberts[x, y] = gx + gy

    sharp = pic_f + Roberts
    sharp = np.where(sharp < 0, 0, np.where(sharp < 0, 0, np.where(sharp > 255, 255, sharp))
    sharp = sharp.astype("uint8")
    return sharp

# 方法二:自定义 这个比较好 ↓
def Robert02(pic):
    ro = [[-1, -1], [1, 1]]
    for i in range(row):
        for j in range(column):
            if(j+2<column) and (i+2<=row):
                process_img = pic[i:i+2, j:j+2]
                list_robert = ro * process_img
                pic[i,j] = abs(list_robert.sum())
    return pic

# 方法三:自定义
def Robert03(pic):
    x_kernel = np.array([[-1,0],[0,1]], dtype=int)
    y_kernel = np.array([[0,-1],[1,0]], dtype=int)
    x = cv2.filter2D(pic, cv2.CV_16S, x_kernel)
    y = cv2.filter2D(pic, cv2.CV_16S, y_kernel)
    absX = cv2.convertScaleAbs(x)
    absY = cv2.convertScaleAbs(y)
    Prewitt = cv2.addWeighted(absX, 0.5, absY, 0.5, 0)
    return Prewitt



# 方法四 直接调用 skimage库内的函数
def sys_robret(pic):
    edge_pic = filters.roberts(pic)
    return edge_pic):
    edge_pic = filters.roberts(pic)
    return edge_pic
    
if __name__ == "__main__":
    pic = cv2.imread('H:\Zivid\pho\e.jpg', 0)
    row, column = pic.shape
    cv2.imshow("original", pic)
    img01 = Robert01(pic)
    cv2.imshow("Robertd-01", img01)
    img02 = Robert02(pic)
    cv2.imshow("Roborts-02", img02)
    img03 = Robert03(pic)
    cv2.imshow("Roborts-03", img03)
    img04 = sys_robret(pic)
    cv2.imshow("Roborts-04", img04)
    cv2.waitKey(0)


Sobel 索贝尔算子

Sobel算子根据像素点上下、左右邻点灰度加权差,在边缘处达到极值这一现象检测边缘。它结合高斯平滑微分求导,用来计算图像灰度函数的近似梯度。在图像的任何一点使用此算子,将会产生对应的梯度矢量或是其法矢量。

它对噪声具有平滑作用,提供较为精确的边缘方向信息,边缘定位精度不够高。当对精度要求不是很高时,是一种较为常用的边缘检测方法。

Sobel算子实现水平边缘检测、垂直边缘检测;45度、135度角边缘检测。

在这里插入图片描述

算法步骤:
该算子包含两组3x3的矩阵,分别为横向及纵向,将之与图像作平面卷积,即可分别得出横向及纵向的亮度差分近似值。如果以A代表原始图像,Gx及Gy分别代表经横向及纵向边缘检测的图像灰度值,其公式如下:

在这里插入图片描述

图像的每一个像素的横向及纵向灰度值通过以下公式结合,来计算该点灰度的大小:

在这里插入图片描述

通常,为了提高效率 使用不开平方的近似值:

在这里插入图片描述

如果梯度G大于某一阀值 则认为该点(x,y)为边缘点。然后可用以下公式计算梯度方向:

在这里插入图片描述

Sobel算子根据像素点上下、左右邻点灰度加权差,在边缘处达到极值这一现象检测边缘。对噪声具有平滑作用,提供较为精确的边缘方向信息,边缘定位精度不够高。当对精度要求不是很高时,是一种较为常用的边缘检测方法。

# python
# 边缘检测 Sober索贝尔算子
# 对噪声具有平滑作用,提供较为精确的边缘方向信息,边缘定位精度不够高。
# @Zivid 2021/8/16

import cv2
import numpy as np
from skimage import filters

# 方法一 自定义
def def_sober(gray_img):
    # 算子
    x_sobel = np.array([[-1, 0, 1],
                        [-2, 0, 2],
                        [-1, 0, 1]])
    y_sobel = np.array([[-1, -2, -1],
                        [0, 0, 0],
                        [1, 2, 1]])

    h, w = gray_img.shape
    img = np.zeros([h + 2, w + 2], np.uint8)
    img[2:h + 2, 2:w + 2] = gray_img[0:h, 0:w]
    x_edge_img = cv2.filter2D(img, -1, x_sobel)
    y_edge_img = cv2.filter2D(img, -1, y_sobel)
    edge_img = np.zeros([h, w])
    for i in range(h):
        for j in range(w):
            edge_img[i][j] = np.sqrt(x_edge_img[i][j] ** 2 + y_edge_img[i][j] ** 2) / (np.sqrt(2))
    return edge_img

# 方法二 系统自带opencv
def cv2_sober(img):
    # 算子
    gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)  # 灰度化
    sobel_x = cv2.Sobel(gray, cv2.CV_8U, 1, 0)  # 通过控制dx,dy的值来控制梯度的方向,这里是求x方向的梯度
    sobel_y = cv2.Sobel(gray, cv2.CV_8U, 0, 1)  # 这里是求y方向的梯度
    sobel = cv2.Sobel(gray, cv2.CV_8U, 1, 1)  # 这里是求x,y方向综合的梯度
    return [sobel_x, sobel_y, sobel]

# 方法三 直接调用系统自带 skimage库内的函数
def skimage_sobel(img):
    edge_pic = filters.sobel(img)
    return edge_pic


if __name__ == '__main__':
    # 读取图像
    img = cv2.imread('H:\Zivid\pho\gray1.jpg', cv2.COLOR_BGR2GRAY)
    rgb_img = cv2.cvtColor(img, cv2.COLOR_BGR2RGB)  # 先转化成RGB
    # 灰度化处理图像
    grayImage = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
    cv2.imshow('original', grayImage)
    pic1 = def_sober(grayImage)
    cv2.imshow('def_sober', pic1)
    pic2 = skimage_sobel(grayImage)
    cv2.imshow('skimage_sobel', pic2)
    sobel_x,sobel_y,sobel = cv2_sober(img)
    cv2.imshow("cv2_Sobel_x", sobel_x)
    cv2.imshow("cv2_Sobel_y", sobel_y)
    cv2.imshow("cv2_Sobel", sobel)
    cv2.waitKey(0)

Prewitt 普利维特算子

Prewitt算子的计算步骤如sober算子,将方向的差分运算和局部平均相结合的方法,也是取水平垂直两个卷积核来分别对图像中各个像素点做卷积运算,所不同的是,Sobel 算子是先做加权平均然后再微分,Prewitt 算子是先平均后求微分

Prewitt 算子通过对图像上的每个像素点的八方向邻域的灰度加权差之和来进行检测边缘,对噪声有一定抑制作用,抗噪性较好,但由于采用了局部灰度平均,容易检测出伪边缘,且边缘定位精度较低。

在这里插入图片描述

# python 
# 边缘检测 Prewitt普利维特算子
# 利用像素点上下、左右邻点灰度差,在边缘处达到极值检测边缘。对噪声具有平滑作用
# @Zivid 2021/8/16
import cv2
import numpy as np
from skimage import filters


# 方法一 自定义
def prewitt01(pic):
    # 算子
    x_kernel = np.array([[1, 1, 1], [0, 0, 0], [-1, -1, -1]], dtype=int)
    y_kernel = np.array([[-1, 0, 1], [-1, 0, 1], [-1, 0, 1]], dtype=int)
    x = cv2.filter2D(pic, cv2.CV_16S, x_kernel)
    y = cv2.filter2D(pic, cv2.CV_16S, y_kernel)
    absX = cv2.convertScaleAbs(x)
    absY = cv2.convertScaleAbs(y)
    Prewitt1 = cv2.addWeighted(absX, 0.5, absY, 0.5, 0)
    return Prewitt1


# 方法二 直接调用系统自带 skimage库内的函数
def sys_prewitt(pic):
    edge_img = filters.prewitt(pic)
    return edge_img


if __name__ == '__main__':
    # 读取图像
    img = cv2.imread('H:\Zivid\pho\gray1.jpg', cv2.COLOR_BGR2GRAY)
    rgb_img = cv2.cvtColor(img, cv2.COLOR_BGR2RGB)  # 先转化成RGB
    # 灰度化处理图像
    grayImage = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
    cv2.imshow('original', grayImage)
    pic1 = prewitt01(grayImage)
    cv2.imshow('Prewitt01', pic1)
    pic2 = sys_prewitt(grayImage)
    cv2.imshow('Prewitt02', pic2)
    cv2.waitKey(0)


Canny 算子

Canny 边缘检测算法是被业界公认的性能最为优良的边缘检测算法之一。

优点:结合了高斯平滑和微分求导,用来计算图像灰度函数的近似梯度
缺点:
(1)为了得到较好的边缘检测结果,它通常需要使用较大的滤波尺度,这样容易丢失一些细节
(2)Canny算子的双阈值要人为的选取,不能够自适应

算法步骤(图片来自作者 拾牙慧者,更详细可见作者 有点方):

来自 拾牙慧者

# python
# 边缘检测 Canny算法
# 步骤:
# (1)灰度化(通常灰度化采用的公式是:Gray=0.299R+0.587G+0.114B;)
# (2)高斯滤波
# (3)计算图像的梯度和梯度方向,(本文使用Sobel算子)
# (4)非极大值抑制(上一步得到的边缘较粗,这里会细化边缘)
# (5)双阈值筛选边缘(二值化显示)
# @Zivid 2021/8/18

import cv2
import numpy as np
# from skimage import filters


# 方法一 自定义
def def_canny(img, top=1, buttom=1, left=1, right=1):
    # prewitt算子。可采用不同的算子,如sober,prewitt
    m1 = np.array([[-1, 0, 1], [-1, 0, 1], [-1, 0, 1]])
    m2 = np.array([[-1, -1, -1], [0, 0, 0], [1, 1, 1]])

    # 第一步 完成高斯平滑
    img = cv2.GaussianBlur(img, (3, 3), 2)  # 高斯滤波

    # 第二步 一阶差分,极端没一点的梯度幅值和方向
    img1 = np.zeros(img.shape, dtype='uint8')  # 与原图大小相同的图片,用于存放梯度值
    theta = np.zeros(img.shape, dtype='float')  # 方向矩阵与原图像大小相同,存放梯度方向
    img = cv2.copyMakeBorder(img, top, buttom, left, right, borderType=cv2.BORDER_CONSTANT, value=0)
    # ↑ 表示增加的像素值,类似补零操作,这里的上下左右都增加了一行补充方式问问诶cv2.BORDER_REPLICATE
    rows, cols = img.shape
    # 卷积操作,以核与图像的卷积代替核之间对应的那个值
    for i in range(1, rows - 1):
        for j in range(1, cols - 1):
            # Gy, Gx相乘再相加
            Gy = (np.dot(np.array([1, 1, 1]), (m1 * img[i - 1:i + 2, j - 1:j + 2]))).dot(np.array([[1], [1], [1]]))
            Gx = (np.dot(np.array([1, 1, 1]), (m2 * img[i - 1:i + 2, j - 1:j + 2]))).dot(np.array([[1], [1], [1]]))
            # 将所有的角度转换到-180到180之间
            if Gx[0] == 0:
                theta[i - 1, j - 1] = 90
                continue
            else:
                # 求梯度方向
                temp = (np.arctan(Gy[0] / Gx[0])) * 180 / np.pi
            if Gx[0] * Gy[0] > 0:
                if Gx[0] > 0:
                    theta[i - 1, j - 1] = np.abs(temp)
                else:
                    theta[i - 1, j - 1] = (np.abs(temp) - 180)
            if Gx[0] * Gy[0] < 0:
                if Gx[0] > 0:
                    theta[i - 1, j - 1] = -1 * np.abs(temp)
                else:
                    theta[i - 1, j - 1] = (180 - np.abs(temp))
            img1[i - 1, j - 1] = (np.sqrt(Gx ** 2 + Gy ** 2))

    # 对梯度方向进行划分,将梯度全部划分为0 90 45 -45四个方向
    for i in range(1, rows - 2):
        for j in range(1, cols - 2):
            if (((theta[i, j] >= -22.5) and (theta[i, j] < 22.5)) or
                    ((theta[i, j] >= 157.5) and (theta[i, j] < 180)) or
                    ((theta[i, j] <= -157.5) and (theta[i, j] >= -180))):
                theta[i, j] = 0.0
            elif (((theta[i, j] >= 22.5) and (theta[i, j] < 67.5)) or
                  ((theta[i, j] <= -112.5) and (theta[i, j] >= -157.5))):
                theta[i, j] = 45.0
            elif (((theta[i, j] >= 67.5) and (theta[i, j] < 112.5)) or
                  ((theta[i, j] <= -67.5) and (theta[i, j] >= -112.5))):
                theta[i, j] = 90.0
            elif (((theta[i, j] >= 122.5) and (theta[i, j] < 157.5)) or
                  ((theta[i, j] <= -22.5) and (theta[i, j] >= -67.5))):
                theta[i, j] = -45.0

    # 第三步 进行非极大值抑制计算
    # 寻找极值点:像素点的值是与该像素点的方向垂直方向上相邻三个像素最大的,则该像素点定义为机制,即边缘像素
    img2 = np.zeros(img1.shape)
    for i in range(1, img2.shape[0] - 1):
        for j in range(1, img2.shape[1] - 1):
            if (theta[i, j] == 0.0) and (img1[i, j] == np.max([img1[i, j], img1[i + 1, j], img1[i - 1, j]])):
                img2[i, j] = img1[i, j]
            if (theta[i, j] == -45.0) and (img1[i, j] == np.max([img1[i, j], img1[i - 1, j - 1], img1[i + 1, j + 1]])):
                img2[i, j] = img1[i, j]
            if (theta[i, j] == 90.0) and (img1[i, j] == np.max([img1[i, j], img1[i, j + 1], img1[i, j - 1]])):
                img2[i, j] = img1[i, j]
            if (theta[i, j] == 45.0) and (img1[i, j] == np.max([img1[i, j], img1[i - 1, j + 1], img1[i + 1, j - 1]])):
                img2[i, j] = img1[i, j]

    # 第四步 双阈值检测和边缘链接
    img3 = np.zeros(img2.shape)  # 定义双阈值图像
    TL = 50
    TH = 100
    # 小于低阈值的设为0,大于高阈值的设为255
    for i in range(1, img3.shape[0] - 1):
        for j in range(1, img3.shape[1] - 1):
            # 将明确的灰度值的像素点绘制
            if img2[i, j] < TL:
                img3[i, j] = 0
            elif img2[i, j] > TH:
                img3[i, j] = 255
            # 如果一个像素点的值在二者之间,且这个像素的临近八个像素点有一个像素值是小于高阈值,则为255
            # 将剩余像素点都设置为255,实现边缘链接
            elif ((img2[i + 1, j] < TH) or (img2[i - 1, j] < TH) or (img2[i, j + 1] < TH) or
                  (img2[i, j - 1] < TH) or (img2[i - 1, j - 1] < TH) or (img2[i - 1, j + 1] < TH) or
                  (img2[i + 1, j + 1] < TH) or (img2[i + 1, j - 1] < TH)):
                img3[i, j] = 255
    return [img1, img2, img3, theta]


# 方法2 opencv自带函数
def cv_canny(img):
    img = cv2.GaussianBlur(img, (3, 3), 2)
    edge_img = cv2.Canny(img, 50, 100)
    return edge_img


if __name__ == '__main__':
    img = cv2.imread('H:\Zivid\pho\ground.jpg', 0)  # 彩色变成灰白
    img1, img2, img3, theta = def_canny(img)
    img_cv = cv_canny(img)

    cv2.imshow('original', img)
    cv2.imshow('img1', img1)
    cv2.imshow('img2', img2)
    cv2.imshow('img3', img3)
    cv2.imshow('theta', theta)
    cv2.imshow('opencv', img_cv)
    cv2.waitKey(0)

Kirsch 算子

Kirsch 算子是一种 3×3 的非线性方向算子。其基本思想是希望改进取平均值的过程,从而尽量使边缘两侧的像素各自与自己同类的像素取平均值,然后再求平均值之差,来减小由于取平均值所造成的边缘细节丢失。通常采用八方向 Kirsch 模板的方法进行检测,取其中最大的值作为边缘强度,而将与之对应的方向作为边缘方向。常用的八方向 Kirsch 模板如下所示,详细可见博主 飘云之下

在这里插入图片描述
对每个像素点都用 这8个模板进行进行卷积(注意,每个卷积值都应取绝对值),在利用 max{ temp0,temp1,temp2,temp3,temp4,temp5,temp6,temp7} 求出该点的最大卷积值。

# 转自 博主 [飘云之下]
import cv2
import numpy as np
from matplotlib import pyplot as plt
#计算Kirsch 边沿检测算子

#定义Kirsch 卷积模板
m1 = np.array([[5, 5, 5],[-3, 0, -3],[-3, -3, -3]])
m2 = np.array([[-3, 5, 5],[-3, 0, 5],[-3, -3, -3]])
m3 = np.array([[-3, -3, 5],[-3, 0, 5],[-3, -3, 5]])
m4 = np.array([[-3, -3, -3],[-3, 0, 5],[-3, 5, 5]])
m5 = np.array([[-3, -3, -3],[-3, 0, -3],[5, 5, 5]])
m6 = np.array([[-3, -3, -3],[5, 0, -3],[5, 5, -3]])
m7 = np.array([[5, -3, -3],[5, 0, -3],[5, -3, -3]])
m8 = np.array([[5, 5, -3],[5, 0, -3],[-3, -3, -3]])
img = cv2.imread("H:\Zivid\pho\gray1.jpg",0)
# img = cv2.GaussianBlur(img,(3,3),5)
#周围填充一圈
#卷积时,必须在原图周围填充一个像素
img = cv2.copyMakeBorder(img,1,1,1,1,borderType=cv2.BORDER_REPLICATE
temp = list(range(8))
img1 = np.zeros(img.shape) #复制空间  此处必须的重新复制一块和原图像矩阵一样大小的矩阵,以保存计算后的结果
for i in range(1,img.shape[0]-1):
    for j in range(1,img.shape[1]-1):
        temp[0] = np.abs( ( np.dot( np.array([1,1,1]) , ( m1*img[i-1:i+2,j-1:j+2]) ) ).dot(np.array([[1],[1],[1]])) )
			#利用矩阵的二次型表达,可以计算出矩阵的各个元素之和
        temp[1] = np.abs(
            (np.dot(np.array([1, 1, 1]), (m2 * img[i - 1:i + 2, j - 1:j + 2]))).dot(np.array([[1],[1],[1]])) )
        temp[2] = np.abs( ( np.dot( np.array([1,1,1]) , ( m1*img[i-1:i+2,j-1:j+2]) ) ).dot(np.array([[1],[1],[1]])) )
        temp[3] = np.abs(
            (np.dot(np.array([1, 1, 1]), (m3 * img[i - 1:i + 2, j - 1:j + 2]))).dot(np.array([[1],[1],[1]])) )
        temp[4] = np.abs(
            (np.dot(np.array([1, 1, 1]), (m4 * img[i - 1:i + 2, j - 1:j + 2]))).dot(np.array([[1],[1],[1]])) )
        temp[5] = np.abs(
            (np.dot(np.array([1, 1, 1]), (m5 * img[i - 1:i + 2, j - 1:j + 2]))).dot(np.array([[1],[1],[1]])) )
        temp[6] = np.abs(
            (np.dot(np.array([1, 1, 1]), (m6 * img[i - 1:i + 2, j - 1:j + 2]))).dot(np.array([[1],[1],[1]])) )
        temp[7] = np.abs(
            (np.dot(np.array([1, 1, 1]), (m7 * img[i - 1:i + 2, j - 1:j + 2]))).dot(np.array([[1],[1],[1]])) )
        img1[i,j] = np.max(temp)
        if img1[i, j] > 255:  #此处的阈值一般写255,根据实际情况选择0~255之间的值
            img1[i, j] = 255
        else:
            img1[i, j] = 0
cv2.imshow("1",img1)
print(img.shape)
cv2.waitKey(0)

二阶

Laplacian 拉普拉斯算子

拉普拉斯算子是 n 维欧几里德空间中的一个二阶微分算子,常用于图像增强领域和边缘提取。它通过灰度差分计算邻域内的像素,查看中间值与周围像素值的差别,如果相差大则判断为边界。

它不依赖于边缘方向的二阶微分算子,对图像中的阶跃型边缘点定位准确,该算子对噪声非常敏感,它使噪声成分得到加强,这两个特性使得该算子容易丢失一部分边缘的方向信息,造成一些不连续的检测边缘,同时抗噪声能力比较差,由于其算法可能会出现双像素边界,常用来判断边缘像素位于图像的明区或暗区。算法详细可见博主 yanghan742915081
在这里插入图片描述

# 边缘检测 二阶微分算子 laplacian算法
# 它通过灰度差分计算邻域内的像素。
# 1)判断图像中心像素灰度值与它周围其他像素的灰度值,如果中心像素的灰度更高,
# 则提升中心像素的灰度;反之降低中心像素的灰度,从而实现图像锐化操作;
# 2)在算法实现过程中,Laplacian算子通过对邻域中心像素的四方向或八方向求梯度,
# 再将梯度相加起来判断中心像素灰度与邻域内其他像素灰度的关系;
# 3)最后通过梯度运算的结果对像素灰度进行调整。
# @Zivid 2021/8/18
import cv2
import numpy as np

# 方法1 自定义函数
def def_laplacian(img):
    rows, cols = img.shape
    image1 = np.zeros(img.shape, dtype='uint8')
    image2 = image1
        lap4_filter = np.array([[0,-1,0],
                            [-1,4,-1],
                            [0,-1,0]])
    for i in range(rows -2):
        for j in range(cols -2):
            image1[i + 1, j + 1] = abs(np.sum(img[i:i + 3, j:j + 3] * lap4_filter))
            if image1[i,j] < 0:

    lap8_filter = np.array([[-1, -1, -1],
                            [-1, 8, -1],
                            [-1, -1, -1]])
    for i in range(rows -2):
        for j in range(cols -2):
            image2[i + 1, j + 1] = abs(np.sum(img[i:i + 3, j:j + 3] * lap8_filter))
            if image2[i,j] < 0:
                image2[i,j] = 0

    image1 = cv2.convertScaleAbs(image1)
    image2 = cv2.convertScaleAbs(image2)
    return image1, image2

# 方法2 opencv系统函数
def cv_laplacian(img):
    lap = cv2.Laplacian(img, cv2.CV_16S, ksize=3)  # 算子的大小,必须为1、3、5、7
    laplacian = cv2.convertScaleAbs(lap)  # 转回uint8
    return laplacian

if __name__ == '__main__':
    img = cv2.imread('H:\Zivid\pho\gray1.jpg',0)  # 彩色变成灰白
    img1, img2 = def_laplacian(img)
    lap = cv_laplacian(img)
    cv2.imshow('original', img)
    cv2.imshow('4阶', img1)
    cv2.imshow('8阶', img2)
    cv2.imshow('laplacian', lap)
    cv2.waitKey(0)
    

当Laplace算子输出出现过零点时就表明有边缘存在,其中忽略无意义的过零点(均匀零区)。原则上,过零点的位置精度可以通过线性内插方法精确到子像素分辨率。
拉普拉斯算子在图像边缘检测中并不常用,主要原因有:

  1. 任何包含有二阶导数的算子比只包含有一阶导数的算子更易受噪声的影响,一阶导数很小的局部峰值也能导致二阶导数过零点,所以Laplace算子对噪声具有无法接受的敏感性
  2. Laplace算子的幅值产生双边元,这是复杂的分割不希望有结果;
  3. Laplace算子不能检测边缘的方向。为了避免噪声的影响,必须采用特别有效的滤波方法。

所以,人们提出了改进的LOG算子来处理Laplace不能处理的离散点和噪声。


LoG 高斯拉普拉斯算子

高斯拉普拉斯算子Lapacian of gaussian先对图像进行高斯平滑滤波处理,然后再与Laplacian算子进行卷积。其处理过程:灰度-高斯-拉普拉斯-负值为0。

更详细的计算过程请查看博主 ingydreamguard 【很详细】

滤波器的选择要考虑以下两个因素:其一是滤波器在空间上要求平稳,即要求空间位置误差 Δ x要小;其二是平滑滤波器本身要求是带通滤波器,并且在有限的带通内是平稳的,即要求频域误差 Δω 要小。

根据信号处理中的测不准原理, Δx 和 Δ ω是相互矛盾的,而达到测不准下限的滤波器就是高斯滤波器。Marr 和 Hildreth 提出的这种差分算子是各向同性的拉普拉斯二阶差分算子。

算法步骤:
Laplace算子通过对图像求取二阶导数的零交叉点(zero-cross)来进行边缘检测,其计算公式如下:

∇ 2 f ( x , y ) = ∂ 2 f ∂ x 2 + ∂ 2 f ∂ y 2 ∇^2f(x,y)= \frac{∂^2f}{∂x^2}+\frac{∂^2f}{∂y^2} 2f(x,y)=x22f+y22f

由于微分运算对噪声比较敏感,可以先对图像进行高斯平滑滤波,再使用Laplace算子进行边缘检测,以降低噪声的影响。由此便形成了用于极值点检测的LoG算子。常用的二维高斯函数如下:

G σ ( x , y ) = 1 2 π σ 2 e x p ( − x 2 + y 2 2 σ 2 ) G_σ(x,y)=\frac{1}{\sqrt{2πσ^2}}exp(−\frac{x^2+y^2}{2σ^2}) Gσ(x,y)=2πσ2 1exp(2σ2x2+y2)

原图像与高斯核函数卷积后再做laplace运算

Δ [ G σ ( x , y ) ∗ f ( x , y ) ] = [ Δ G σ ( x , y ) ] ∗ f ( x , y ) Δ[Gσ(x,y)∗f(x,y)]=[ΔGσ(x,y)]∗f(x,y) Δ[Gσ(x,y)f(x,y)]=[ΔGσ(x,y)]f(x,y)
L o G = Δ G σ ( x , y ) = ∂ 2 G σ ( x , y ) ∂ x 2 + ∂ 2 G σ ( x , y ) ∂ y 2 = x 2 + y 2 − 2 σ 2 σ 4 e − ( x 2 + y 2 ) / 2 σ 2 LoG=ΔGσ(x,y)=\frac{∂^2Gσ(x,y)}{∂x^2}+\frac{∂^2Gσ(x,y)}{∂y^2}=\frac{x^2+y^2−2σ^2}{σ^4} e^{−(x^2+y^2)/2σ^2} LoG=ΔGσ(x,y)=x22Gσ(x,y)+y22Gσ(x,y)=σ4x2+y22σ2e(x2+y2)/2σ2

该边缘检测器的基本特征是:
(1) 所用的平滑滤波器是高斯滤波器
(2) 增强步骤采用的是二阶导数(即二维拉普拉斯函数)
(3) 边缘检测的判据是二阶导数过零点并且对应一阶导数的极大值
该方法的特点是先用高斯滤波器与图像进行卷积,既平滑了图像又降低了噪声,使孤立的噪声点和较小的结构组织被滤除。然而由于对图像的平滑会导致边缘的延展,因此只考虑那些具有局部梯度极大值的点作为边缘点,这可以用二阶导数的零交叉来实现。

# python
# 边缘检测 二阶微分算子 LoG算法
# 先对图像进行高斯平滑处理,然后再与Laplacian算子进行卷积
# @Zivid 2021/8/18
import cv2
import numpy as np


def def_log(img):
    rows, cols = img.shape
    # 第一步:高斯滤波
    gau_filter = np.array([[0, 0, 1, 0, 0],
                           [0, 1, 2, 1, 0],
                           [1, 2, -16, 2, 1],
                           [0, 1, 2, 1, 0],
                           [0, 0, 1, 0, 0]])
    # gau_filter = np.array([[1,2,1],
    #                        [2,4,2],
    #                        [1,2,1]])
    ex_pad1 = np.pad(img, ((2, 2), (2, 2)), 'constant')  # 扩展操作
    # 逐个像素进行高斯滤波操作
    for i in range(rows - 4):
        for j in range(cols - 4):
            ex_pad1[i, j] = np.sum(ex_pad1[i:i + 5, j:j + 5] * gau_filter)

    # 第二步:计算laplacian的二阶导数,
    # 四邻域
    lap4_filter = np.array([[0, -1, 0],
                            [-1, 4, -1],
                            [0, -1, 0]])
    lap4_pad = np.pad(ex_pad1, ((1, 1), (1, 1)), 'constant')
    image1 = np.zeros(img.shape, dtype='uint8')
    for i in range(rows - 2):
        for j in range(cols - 2):
            image1[i, j] = np.sum(lap4_pad[i:i + 3, j:j + 3] * lap4_filter)
            if image1[i, j] < 0:
                image1[i, j] = 0
    # 八邻域
    lap8_filter = np.array([[-1, -1, -1],
                            [-1, 8, -1],
                            [-1, -1, -1]])
    lap8_pad = np.pad(ex_pad1, ((1, 1), (1, 1)), 'constant')
    image2 = np.zeros(img.shape, dtype='uint8')
    for i in range(1, rows - 1):
        for j in range(1, cols - 1):
            image2[i, j] = np.sum(lap8_pad[i - 1:i + 2, j - 1:j + 2] * lap8_filter)
            if image2[i, j] < 0:
                image2[i, j] = 0

    return ex_pad1, image1, image2

def sys_log(img):
    # 以核为3*3,方差为0的高斯函数进行高斯滤波
    image1 = cv2.GaussianBlur(img, (3, 3), sigmaX=0)
    sys_img = cv2.Laplacian(image1, cv2.CV_16S, ksize=3)
    sys_img = cv2.convertScaleAbs(sys_img)
    return sys_img


if __name__ == '__main__':
    img = cv2.imread('H:\Zivid\pho\gray1.jpg', 0)  # 彩色变成灰白
    gau, img1, img2 = def_log(img)
    lap = sys_log(img)
    cv2.imshow('original', img)
    cv2.imshow('gause', gau)
    cv2.imshow('4阶', img1)
    cv2.imshow('8阶', img2)
    cv2.imshow('laplacian of gaussian', lap)
    cv2.waitKey(0)
                      

DoG 高斯函数差分

Difference of Gaussian(DOG)是高斯函数的差分。它是可以通过将图像与高斯函数进行卷积得到一幅图像的低通滤波结果,即去噪过程,这里的Gaussian和高斯低通滤波器的高斯一样,是一个函数,即为正态分布函数。同时,它对高斯拉普拉斯LoG的近似,在某一尺度上的特征检测可以通过对两个相邻高斯尺度空间的图像相减,得到DoG的响应值图像。

常见的二维高斯函数公式为:

G σ 1 ( x , y ) = 1 2 π σ 1 2 e x p ( − x 2 + y 2 2 σ 1 2 ) G_{σ1}(x,y)=\frac{1}{\sqrt{2πσ_1^2}}exp(−\frac{x^2+y^2}{2σ_1^2}) Gσ1(x,y)=2πσ12 1exp(2σ12x2+y2)

其次,两幅图像的高斯滤波表示为:

Δ g 1 ( x , y ) = G σ 1 ( x , y ) ∗ f ( x , y ) ; Δ g 2 ( x , y ) = G σ 2 ( x , y ) ∗ f ( x , y ) Δg_1(x,y)=Gσ_1(x,y)∗f(x,y) ; Δg_2(x,y)=Gσ_2(x,y)∗f(x,y) Δg1(x,y)=Gσ1(x,y)f(x,y)Δg2(x,y)=Gσ2(x,y)f(x,y)

最后,将滤波得到的两幅图g1和g2相减得到:

Δ g 1 ( x , y ) − Δ g 2 ( x , y ) = G σ 1 ( x , y ) ∗ f ( x , y ) − G σ 2 ( x , y ) ∗ f ( x , y ) = ( G σ 1 − G σ 2 ) ∗ f ( x , y ) = D o G ∗ f ( x , y ) Δg_1(x,y)-Δg_2(x,y)=Gσ_1(x,y)∗f(x,y)-Gσ_2(x,y)∗f(x,y)=(Gσ_1-Gσ_2)*f(x,y)=DoG*f(x,y) Δg1(x,y)Δg2(x,y)=Gσ1(x,y)f(x,y)Gσ2(x,y)f(x,y)=(Gσ1Gσ2)f(x,y)=DoGf(x,y)

则由上式得到:

D o G = G σ 1 − G σ 2 = 1 2 π ( 1 σ 1 e − ( x 2 + y 2 ) / 2 σ 1 2 − 1 σ 2 e − ( x 2 + y 2 ) / 2 σ 2 2 ) DoG = G_{σ1}-G_{σ2}=\frac{1}{\sqrt{2π}}(\frac{1}{σ_1}e^{-(x^2+y^2)/2σ_1^2}-\frac{1}{σ_2}e^{-(x^2+y^2)/2σ_2^2}) DoG=Gσ1Gσ2=2π 1(σ11e(x2+y2)/2σ12σ21e(x2+y2)/2σ22)

DoG和LoG对比:
在这里插入图片描述

# python
# 边缘检测 二阶微分算子 dog算法
# 确定σ值来对一幅图像进行滤波,所得零交叉边缘图与仅为全部图形保留的公共边缘相结合
# @Zivid 2021/8/18

from skimage import filters
import cv2
import matplotlib.pyplot as plt

# 方法:采用cv或skimage库自带高斯函数
def sys_dog(img):
    ## 先对图像进行两次高斯滤波运算
    # 1.采用skimage进行高斯滤波
    gimg1 = filters.gaussian(img, sigma=6)
    gimg2 = filters.gaussian(img, sigma=3)
    # 2.采用cv进行高斯滤波
    # gimg1 = cv2.GaussianBlur(img, (3, 3), sigmaX=6)
    # gimg2 = cv2.GaussianBlur(img, (3, 3), sigmaX=3)
    # 两个高斯运算的差分
    dimg = gimg1 - gimg2
    return gimg1, gimg2, dimg


if __name__ == '__main__':
    gray_img = cv2.imread('H:\Zivid\pho\gray1.jpg', 0)
    gimg1, gimg2, dimg = sys_dog(gray_img)
    cv2.imshow('ori', gray_img)
    cv2.imshow('gaus1', gimg1)
    cv2.imshow('gaus2', gimg2)
    cv2.imshow('gaus1-2', dimg)

    cv2.waitKey(0)


感谢以下博主,都是大佬:
菜菜粥brycezou有情怀的机械男zhulu_20TechArtisan6疯狂的Little_BeeSenit_Co ,等等等

  • 23
    点赞
  • 114
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
高斯差分算子(Difference of Gaussian,DoG)是一种常用边缘检测算法。它是通过对进行高斯滤波后,计算不同尺度下的高斯滤波结果之间的差异来实现的。具体步骤如下: 1. 对原始进行高斯滤波,得到不同尺度下的高斯模糊像; 2. 计算相邻两个尺度下的高斯模糊像之差,得到一组差分像; 3. 对每个差分进行非极大值抑制,得到一组非极大值抑制像; 4. 对所有非极大值抑制进行二值化处理,得到一组二值化像; 5. 将所有二值化进行叠加,得到最终的边缘检测结果。 以下是Python实现高斯差分边缘检测算子代码: ```python import cv2 import numpy as np # 读取像并加上高斯噪声 img = cv2.imread('lena.jpg', 0) img = cv2.GaussianBlur(img, (5, 5), 0) img = img + np.random.normal(0, 25, img.shape) # 定义高斯差分函数 def DoG(img, ksize, sigma1, sigma2): g1 = cv2.GaussianBlur(img, ksize, sigma1) g2 = cv2.GaussianBlur(img, ksize, sigma2) return g1 - g2 # 计算高斯差分dog1 = DoG(img, (5, 5), 1, 2) dog2 = DoG(img, (5, 5), 2, 4) dog3 = DoG(img, (5, 5), 4, 8) # 非极大值抑制 def non_max_suppression(img): h, w = img.shape out = np.zeros((h, w), dtype=np.float32) for y in range(1, h-1): for x in range(1, w-1): dx = img[y, x+1] - img[y, x-1] dy = img[y+1, x] - img[y-1, x] gradient = np.sqrt(dx**2 + dy**2) if gradient == 0: out[y, x] = 0 else: angle = np.rad2deg(np.arctan(dy/dx)) if angle < 0: angle += 180 if (angle <= 22.5) or (angle > 157.5): if (img[y, x] >= img[y, x+1]) and (img[y, x] >= img[y, x-1]): out[y, x] = img[y, x] elif (22.5 < angle <= 67.5): if (img[y, x] >= img[y-1, x+1]) and (img[y, x] >= img[y+1, x-1]): out[y, x] = img[y, x] elif (67.5 < angle <= 112.5): if (img[y, x] >= img[y-1, x]) and (img[y, x] >= img[y+1, x]): out[y, x] = img[y, x] else: if (img[y, x] >= img[y-1, x-1]) and (img[y, x] >= img[y+1, x+1]): out[y, x] = img[y, x] return out # 非极大值抑制像 nms1 = non_max_suppression(dog1) nms2 = non_max_suppression(dog2) nms3 = non_max_suppression(dog3) # 二值化处理 th1 = cv2.threshold(nms1, 0, 255, cv2.THRESH_BINARY+cv2.THRESH_OTSU)[1] th2 = cv2.threshold(nms2, 0, 255, cv2.THRESH_BINARY+cv2.THRESH_OTSU)[1] th3 = cv2.threshold(nms3, 0, 255, cv2.THRESH_BINARY+cv2.THRESH_OTSU)[1] # 叠加所有二值化像 result = th1 + th2 + th3 # 显示结果 cv2.imshow('result', result) cv2.waitKey(0) cv2.destroyAllWindows() ```

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值