常见边缘检测算子极其实现(基于Python和OpenCV)

Posted Hello World

tags:

篇首语:本文由小常识网(cha138.com)小编为大家整理,主要介绍了常见边缘检测算子极其实现(基于Python和OpenCV)相关的知识,希望对你有一定的参考价值。

1.边缘检测类型和基本原理

在图像处理中,图像边缘常包括三种模型

1).台阶模型:相邻两个像素的灰度值快速变化;如:在离散灰度图像中,灰度值为...0,0,255,255,...就可视为台阶模型;

2).斜坡模型:图像中从亮到暗(或暗到亮)呈现一个类似于斜坡,如在离散灰度图像灰度:...,0,50,100,150,200,255,255...;

3).屋顶模型:...0,0,80,160,255,160,80,0,0,0...

如下图所示

基本原理

边缘检测简单来说,就是对图像进行按模板进行卷积;边缘在图像上一般就表示为明暗交接的地方,从上面的图中看出,若我们将图像看为一个连续可导的函数,那么,对边缘检测,就是求突变的位置,在数学上,我们对这样连续可导图像是通过求导,因此,边缘也是通过同样的原理。但由于图像是离散矩阵形式,求导采用微分的形式。边缘检测中常见的数学方法如下图所示

2.常见的边缘检测算子

2.1Roberts算子

当对对角线方向的边缘感兴趣时,需要一个二维模板,Roberts交叉梯度算子是最早尝试使用这样的具有对角优势的二维算子,具体的函数实现如下

    #传入灰度图像
    #返回处理后图像
    def robert(self,img):
        #获取边界模板
        roberX=roberY=np.zeros((2,2))
        roberX[0][0]=-1;roberX[1][1]=1
        roberY[0][1]=-1;roberY[1][0]=1
        #计算图像大小
        m,n=np.shape(img)
        gx=0
        gy=0
        M=img[:]
        for i in range(m):
            for j in range(n):
                gx=np.sum(img[i:i+2,j:j+2]*roberX)
                gy=np.sum(img[i:i+2,j:j+2]*roberY)
                M[i][j]=abs(gx)+abs(gy)
        return M
View Code

2.2Prewitt算子

    #传入灰度图像
    #返回处理后图像
    def prewitt(self,img):
        prewitX=[[-1,-1,-1],
                 [0,0,0],
                 [1,1,1]]
        prewitY=[[-1,0,1],
                 [-1,0,1],
                 [-1,0,1]]
        m,n=np.shape(img)
        M=img[:]
        gx=0
        gy=0
        for i in range(m-3):
            for j in range(n-3):
                gx=np.sum(img[i:i+3,j:j+3]*prewitX)
                gy=np.sum(img[i:i+3,j:j+3]*prewitY)
                M[i][j]=abs(gx)+abs(gy)
        return M
View Code

2.3Sobel算子

可以看出两个算子计算相差并不大,

    #传入灰度图像
    #返回处理后图像
    def sobel(self,img):
        #计算边界模板
        sobelX=np.array([[-1,-2,-1],
            [0,0,0],
            [1,2,1]])
        sobelY=np.array([[-1,0,1],
                          [-2,0,2],
                          [-1,0,1]])
        m,n=np.shape(img)
        gx=0
        gy=0
        #初始化与原图像相同数组
        M=np.zeros_like(img)
        for i in range(1,m-1):
            for j in range(1,n-1):
                gx=np.sum(img[i-1:i+2,j-1:j+2]*sobelX)
                gy=np.sum(img[i-1:i+2,j-1:j+2]*sobelY)
                M[i][j]=abs(gx)+abs(gy)
        return M
View Code

3.LOG算子

LOG边缘检测步骤:

1.先用n*n的高斯进行平滑处理;

2.对平滑后的图像再使用拉普拉斯算子进行边缘检测;

3.再找出图像的零交叉点,进行边界连接。

传统的边缘检测算法,实质就是一些模板(又称为算子)在图像上滑动,LOG算子也不例外,其模板计算公式如下图所示:

对于参数δ是要输入的参数,因此,对于模板大小n>=6δ原则,以上函数图像如下,因长的像草帽,又叫墨西哥草帽算子。

图像零交叉点的查找是该方法的核心部分之一,从上述公式可以看出,这是对原始图像进行二阶导数处理过程。若我们考虑上面提到的边缘检测的斜坡模型,求导数的图像如下所示

最左边的为原始图像,中间部分为一阶导数图像,经过一阶导后,斜坡部分为一个常数,最右边为二级求导后图像,二阶处理后斜坡部分变为0,但其进入斜坡的两个端点处,却不为零,而且方向相反(一正一负)。变化部分如图像中箭头所示。因此。我们可以基于该原理,进行边界提取。即通过3*3模板进行边界查找,判断经过上面模板处理后图像为0像素的上下、左右,两个对角4个方向是否相反,若相反,则标记为边缘。对于台阶模型边缘只需其上下、左右4个方向上符号相反即可。

#MarrHildreth边界检测
class marrHildreth():
    def edgesMarrHildreth(self,img,sigma):
        #动态确定滤波器大小,大于6sigma的奇数
        size=int(2*(np.ceil(3*sigma))+1)
        x,y=np.meshgrid(np.arange(-size/2-1,size/2+1),np.arange(-size/2-1,size/2+1))
        # normal=1/(2.0*np.pi*sigma**2)
        #计算LoG滤波器
        kernel=((x**2+y**2-(2.0*sigma**2))/sigma**4)*np.exp(-(x**2+y**2)/(2.0*sigma**2))
        kern_size=kernel.shape[0]

        m,n=np.shape(img)
        log=np.zeros_like(img,dtype=float)

        t=int(kern_size/2)#定义起始点

        for i in range(t,m-t):
            for j in range(t,n-t):
                window=img[i-t:i+(t+1),j-t:j+(t+1)]*kernel
                log[i,j]=np.sum(window)
        log=log.astype(np.int64,copy=False)
        zero_crossing=np.zeros_like(log)
        #寻找零交叉点
        for i in range(1,log.shape[0]-1):
            for j in range(1,log.shape[1]-1):
                #斜坡型,使用3*3模板,零交叉点处,至少相邻俩个方向二阶导数符号相反,
                #上/下;左/右,两个对角方向
                if log[i][j]==0:
                    if(log[i][j-1] < 0 and log[i][j+1] > 0) \\
                            or (log[i][j-1] > 0 and log[i][j+1] < 0) \\
                            or (log[i-1][j] < 0 and log[i+1][j] > 0) \\
                            or (log[i-1][j] > 0 and log[i+1][j] < 0) \\
                            or ((log[i - 1][i - 1] > 0) and (log[i + 1][j + 1] < 0)) \\
                            or ((log[i - 1][i - 1] < 0) and (log[i + 1][j + 1] > 0)):
                        zero_crossing[i][j]=255
                #阶梯型,只要两个相邻(4领域)中符号相反
                if log[i][j]<0:
                    if(log[i][j-1]>0)or(log[i][j+1]>0)or(log[i-1][j]>0)or(log[i+1][j]>0):
                        zero_crossing[i][j]=255
        return zero_crossing
View Code

可以看出LOG结果并不理想,这应该是程序设计有问题!

4.Canny算子,最经典的一种算子之一

基本步骤

1.用二维高斯平滑图像;由于二维高斯过程计算比较复杂,通常可以将其近似看为沿x,y两个方向上计算;

2.计算图像梯度;

3.非最大抑制;

4.双阈值进行边界提取

高斯平滑、计算梯度和角度的公式如上所示。这里还有一个重要的性质是梯度的方向和边缘的方向相互垂直。

    #将RGB图像转化为灰度图像
    def gray(self,img_path):
        img=plt.imread(img_path)
        img_rgb=cv2.cvtColor(img,cv2.COLOR_BGR2RGB)
        img_gray=np.dot(img_rgb[...,:3],[0.299,0.587,0.114])
        return img_gray

    #高斯滤波平滑图像
    def smooth(self,img_gray,sigmal1=1.4,sigmal2=1.4):
        sigmal1=sigmal1
        sigmal2=sigmal2
        gau_sum=0
        gaussian=np.zeros([5,5])
        for i in range(5):
            for j in range(5):
                gaussian[i,j]=math.exp((-1/(2*sigmal2*sigmal1))*(np.square(i-3)+
                                                               np.square(j-3)))/(2*math.pi*sigmal1*sigmal2)
                gau_sum=gau_sum+gaussian[i,j]
            gaussian=gaussian/gau_sum
            W,H=img_gray.shape
            new_gray=np.zeros([W-5,H-5])
            for i in range(W-5):
                for j in range(H-5):
                    new_gray[i,j]=np.sum(img_gray[i:i+5,j:j+5]*gaussian)
        return new_gray
    #计算dx,dy,M,thana
    def gradients(self,new_gray):
        W,H=new_gray.shape
        dx=np.zeros([W-1,H-1])
        dy=np.ones([W-1,H-1])
        M=np.zeros([W-1,H-1])
        theta=np.zeros([W-1,H-1])
        for i in range(W-1):
            for j in range(H-1):
                dx[i,j]=new_gray[i+1,j]-new_gray[i,j]
                dy[i,j]=new_gray[i,j+1]-new_gray[i,j]
                M[i,j]=np.sqrt(np.square(dx[i,j])+np.square(dy[i,j]))
                theta[i,j]=math.atan(dx[i,j]/(dy[i,j]+0.0000001))
        return dx,dy,M,theta
View Code

非最大值抑制:

核心:梯度是有方向矢量,通常在进行非最大抑制时使用3*3模板进行查找该方向线上的两个向量,要求给定像元梯度要大于这两个方向线上的梯度值。

下面介绍两种非最大值抑制的方法

1)基于灰度图像插值的方法

从上图中可以看到,如前面介绍,梯度方向与边缘相互垂直,即,图中短线代表边缘,长的带箭头的表示梯度方向,在图中所示,上面部分在b中,但离a较近;下面部分在g中,但其离h较近,这时,我们不能简单的将上面部分看为b,下半部视为g,因该通过插值的方法通过a、b计算上面部分值,同样通过下半部分g、h计算下面值。

 

将以上过程简化为上图所示,图中虚线代表梯度方向,中间两个图均有一个共同的性质,即沿y轴方向导数大于沿x轴方向导数。同理,最下面两个图沿x轴方向导数大于沿y方向导数。首先我们要明确我们的目标是插值计算出图中g,h值。以中间的两个图为例,前面介绍了中间两个图沿y轴方向导数大于x轴;那么我们就可以明确在该条件下首先可以获取到上下两个值。之后还要进一步判断,中间两个图中左图:两个方向(x,y)上的梯度相乘大于0(判断方法:梯度有方向,首先假设梯度方向,再看x、y方向);中间两个图中右图,两个方向相乘小于零。这样就可以判断出要进行插值的另外两个点。

    #非最大抑制方案2
    def NMS(self,M,dx,dy):
        #M:梯度
        #dx:x方向导数
        #dy:y方向导数
        d=np.copy(M)
        W,H=M.shape
        NMS=np.zeros((W,H))
        NMS[0,:]=NMS[W-1,:]=NMS[:,0]=NMS[:,H-1]=0
        for i in range(1,W-1):
            for j in range(1,H-1):
                #找到可能边缘点
                if M[i,j]!=0:
                    #x方向导数
                    gradx=dx[i,j]
                    #y方向导数
                    grady=dy[i,j]
                    #当前梯度
                    gradTemp=M[i,j]
                    #如果 x 方向梯度值比较大,说明导数方向趋向于 x 分量
                    if np.abs(gradx)>np.abs(grady):
                        #权重
                        weight=np.abs(grady)/np.abs(gradx)
                        grad2=M[i-1,j]
                        grad4=M[i+1,j]
                        if gradx*grady>0:
                            # 如果 x, y 方向导数符号一致
                            # 像素点位置关系
                            # g1
                            # g2  c  g4
                            #        g3
                            grad1=M[i-1,j-1]
                            grad3=M[i+1,j+1]
                        else:
                            # 如果 x,y 方向导数符号相反
                            # 像素点位置关系
                            #        g3
                            # g2  c  g4
                            # g1
                            grad1=M[i-1,j+1]
                            grad3=M[i+1,j-1]
                    else:
                        weight=np.abs(gradx)/np.abs(grady)
                        grad2=M[i,j-1]
                        grad4=M[i,j+1]
                        # 如果 x, y 方向导数符号相反
                        # 像素点位置关系
                        #    g2 g1
                        #    c
                        # g3 g4
                        if gradx*grady<0:
                            grad1=M[i+1,j-1]
                            grad3=M[i-1,j+1]
                        # 如果 x,y 方向导数符号一致
                        # 像素点位置关系
                        # g1 g2
                        #    c
                        #    g4  g3
                        else:
                            grad1=M[i-1,j-1]
                            grad3=M[i+1,j+1]
                    #分别计算两个方向的灰度差值
                    gradTemp1=weight*grad2+(1-weight)*grad1
                    gradTemp2=weight*grad4+(1-weight)*grad3
                    #若该梯度均大于差值的灰度值,则保留,否则赋值为零,抑制
                    if gradTemp>=gradTemp1 and gradTemp>=gradTemp2:
                        NMS[i,j]=gradTemp
        return NMS
View Code

2)基于角度非最大值抑制方法

前面提到图像边界与梯度方向相互垂直。

基本流程:

1).令d1、d2、d3、d4分别表示水平、垂直、-45°、+45°角度范围。

2).同样采用3*3模板进行非极大值抑制,但现在使用的是角度。

3).首先根据角度判断是在什么方向。

4).中间像素要大于相邻俩个方向上的值。

    #非最大抑制方案1
    def NMS1(self,M,dx,dy,theta):
        #定义常数,将弧度转化为度
        gree=180.0/math.pi
        W,H=np.shape(M)
        #获取角度值
        d = np.copy(theta)
        NMS=np.zeros((W,H))
        NMS[0,:]=NMS[W-1,:]=NMS[:,0]=NMS[:,H-1]=0
        for i in range(1,W-1):
            for j in range(1,H-1):
                if M[i,j]!=0:
                    gradTemp=d[i][j]*gree
                    #水平边缘
                    if((-157.5<=gradTemp<=157.5) or (-22.5<=gradTemp<=22.5)):
                        if(M[i][j]>M[i-1][j] and M[i][j]>M[i+1][j]):
                            NMS[i][j]=M[i][j]
                    #45°边缘
                    if((112.5<=gradTemp<=157.5) or (-67.5<=gradTemp<=-22.5)):
                        if(M[i][j]>M[i+1][j-1] and(M[i][j]>M[i-1][j+1])):
                            NMS[i][j]=M[i][j]
                    #垂直边缘
                    if((-112.5<=gradTemp<=-67.5) or (67.5<=gradTemp<=122.5)):
                        if(M[i][j]>M[i-1][j]and(M[i][j]>M[i][j+1])):
                            NMS[i][j]=M[i][j]
                    #-45°边缘
                    if((-157.5<=gradTemp<=112.5) or (22.5<=gradTemp<=67.5)):
                        if(M[i][j]>M[i-1][j-1]and(M[i][j]>M[i+1][j+1])):
                            NMS[i][j]=M[i][j]
        return NMS
View Code

双阈值边界提取

1).Canny算子采用的是高低俩个阈值来进行边界连接,通常高低阈值取值比为2:1或3:1,即TH:TL=2:1(或3:1);

2).高于TH高阈值的点认为是边界点,而低于低阈值点认为不可能为边界点,

3).高于低阈值,低于高阈值间的点,要进行判断(这里采用3*3模板)在该点的8个邻域中是否有高于高阈值(TH)点,若有,则将该点视为边界点

#双阈值处理
    def double_threshold(self,NMS):
        W,H=NMS.shape
        DT=np.zeros([W,H])
        #高低阈值比为:1:3
        TL=0.1*np.max(NMS)
        TH=0.3*np.max(NMS)

        for i in range(1,W-1):
            for j in range(1,H-1):
                #将低于最小阈值都设为0
                if(NMS[i,j]<TL):
                    DT[i,j]=0
                #将大于最高阈值都设为1
                elif(NMS[i,j]>=TH):
                    DT[i,j]=255

                #用8连通方法得到高低阈值间的边缘像素
                elif((NMS[i-1,:]>=TH).any()) or ((NMS[i+1,:]>=TH).any()) or ((NMS[i,[j-1,j+1]]>=TH).any()):
                    DT[i,j]=255
                # elif((NMS[i-1,[j-1,j+1]>TH].any())or((NMS[i+1,[j-1,j+1]]>TH).any())or((NMS[i,[j-1,j+1]]>TH).any())):
                #     DT[i,j]=255
        return DT
View Code

Canny算子总代码

#Canny边缘检测
class canny():
    #将RGB图像转化为灰度图像
    def gray(self,img_path):
        img=plt.imread(img_path)
        img_rgb=cv2.cvtColor(img,cv2.COLOR_BGR2RGB)
        img_gray=np.dot(img_rgb[...,:3],[0.299,0.587,0.114])
        return img_gray

    #高斯滤波平滑图像
    def smooth(self,img_gray,sigmal1=1.4,sigmal2=1.4):
        sigmal1=sigmal1
        sigmal2=sigmal2
        gau_sum=0
        gaussian=np.zeros([5,5])
        for i in range(5):
            for j in range(5):
                gaussian[i,j]=math.exp((-1/(2*sigmal2*sigmal1))*(np.square(i-3)+
                                                               np.square(j-3)))/(2*math.pi*sigmal1*sigmal2)
                gau_sum=gau_sum+gaussian[i,j]
            gaussian=gaussian/gau_sum
            W,H=img_gray.shape
            new_gray=np.zeros([W-5,H-5])
            for i in range(W-5):
                for j in range(H-5):
                    new_gray[i,j]=np.sum(img_gray[i:i+5,j:j+5]*gaussian)
        return new_gray
    #计算dx,dy,M,thana
    def gradients(self,new_gray):
        W,H=new_gray.shape
        dx=np.zeros([W-1,H-1])
        dy=np.ones([W-1,H-1])
        M=np.zeros([W-1,H-1])
        theta=np.zeros([W-1,H-1])
        for i in range(W-1):
            for j in range(H-1):
                dx[i,j]=new_gray[i+1,j]-new_gray[i,j]
                dy[i,j]=new_gray[i,j+1]-new_gray[i,j]
                M[i,j]=np.sqrt(np.square(dx[i,j])+np.square(dy[i,j]))
                theta[i,j]=math.atan(dx[i,j]/(dy[i,j]+0.0000001))
        return dx,dy,M,theta

    #非最大抑制方案1
    def NMS1(self,M,dx,dy,theta):
        #定义常数,将弧度转化为度
        gree=180.0/math.pi
        W,H=np.shape(M)
        #获取角度值
        d = np.copy(theta)
        NMS=np.zeros((W,H))
        NMS[0,:]=NMS[W-1,:]=NMS[:,0]=NMS[:,H-1]=0
        for i in range(1,W-1):
            for j in range(1,H-1):
                if M[i,j]!=0:
                    gradTemp=d[i][j]*gree
                    #水平边缘
                    if((-157.5<=gradTemp<=157.5) or (-22.5<=gradTemp<=22.5)):
                        if(M[i][j]>M[i-1][j] and M[i][j]>M[i+1][j]):
                            NMS[i][j]=M[i][j]
                    #45°边缘
                    if((112.5<=gradTemp<=157.5) or (-67.5<=gradTemp<=-22.5)):
                        if(M[i][j]>M[i+1][j-1] and(M[i][j]>M[i-1][j+1])):
                            NMS[i][j]=M[i][j]
                    #垂直边缘
                    if((-112.5<=gradTemp<=-67.5) or (67.5<=gradTemp<=122.5)):
                        if(M[i][j]>M[i-1][j]and(M[i][j]>M[i][j+1])):
                            NMS[i][j]=M[i][j]
                    #-45°边缘
                    if((-157.5<=gradTemp<=112.5) or (22.5<=gradTemp<=67.5)):
                        if(M[i][j]>M[i-1][j-1]and(M[i][j]>M[i+1][j+1])):
                            NMS[i][j]=M[i][j]
        return NMS

    #非最大抑制方案2
    def NMS(self,M,dx,dy):
        #M:梯度
        #dx:x方向导数
        #dy:y方向导数
        d=np.copy(M)
        W,H=M.shape
        NMS=np.zeros((W,H))
        NMS[0,:]=NMS[W-1,:]=NMS[:,0]=NMS[:,H-1]=0
        for i in range(1,W-1):
            for j in range(1,H-1):
                #找到可能边缘点
                if M[i,j]!=0:
                    #x方向导数
                    gradx=dx[i,j]
                    #y方向导数
                    grady=dy[i,j]
                    #当前梯度
                    gradTemp=M[i,j]
                    #如果 x 方向梯度值比较大,说明导数方向趋向于 x 分量
                    if np.abs(gradx)>np.abs(grady):
                        #权重
                        weight=np.abs(grady)/np.abs(gradx)
                        grad2=M[i-1,j]
                        grad4=M[i+1,j]
                        if gradx*grady>0:
                            # 如果 x, y 方向导数符号一致
                            # 像素点位置关系
                            # g1
                            # g2  c  g4
                            #        g3
                            grad1=M[i-1,j-1]
                            grad3=M[i+1,j+1]
                        else:
                            # 如果 x,y 方向导数符号相反
                            # 像素点位置关系
                            #        g3
                            # g2  c  g4
                            # g1
                            grad1=M[i-1,j+1]
                            grad3=M[i+1,j-1]
                    else:
                        weight=np.abs(gradx)/np.abs(grady)
                        grad2=M[i,j-1]
                        grad4=M[i,j+1]
                        # 如果 x, y 方向导数符号相反
                        # 像素点位置关系
                        #    g2 g1
                        #    c
                        # g3 g4
                        if gradx*grady<0:
                            grad1=M[i+1,j-1]
                            grad3=M[i-1,j+1]
                        # 如果 x,y 方向导数符号一致
                        # 像素点位置关系
                        # g1 g2
                        #    c
                        #    g4  g3
                        else:
                            grad1=M[i-1,j-1]
                            grad3=M[i+1,j+1]
                    #分别计算两个方向的灰度差值
                    gradTemp1=weight*grad2+(1-weight)*grad1
                    gradTemp2=weight*grad4+(1-weight)*grad3
                    #若该梯度均大于差值的灰度值,则保留,否则赋值为零,抑制
                    if gradTemp>=gradTemp1 and gradTemp>=gradTemp2:
                        NMS[i,j]=gradTemp
        return NMS

    #双阈值处理
    def double_threshold(self,NMS):
        W,H=NMS.shape
        DT=np.zeros([W,H])
        #高低阈值比为:1:3
        TL=0.1*np.max(NMS)
        TH=0.3*np.max(NMS)

        for i in range(1,W-1):
            for j in range(1,H-1):
                #将低于最小阈值都设为0
                if(NMS[i,j]<TL):
                    DT[i,j]=0
                #将大于最高阈值都设为1
                elif(NMS[i,j]>=TH):
                    DT[i,j]=255

                #用8连通方法得到高低阈值间的边缘像素
                elif((NMS[i-1,:]>=TH).any()) or ((NMS[i+1,:]>=TH).any()) or ((NMS[i,[j-1,j+1]]>=TH).any()):
                    DT[i,j]=255
                # elif((NMS[i-1,[j-1,j+1]>TH].any())or((NMS[i+1,[j-1,j+1]]>TH).any())or((NMS[i,[j-1,j+1]]>TH).any())):
                #     DT[i,j]=255
        return DT
View Code

 

以上是关于常见边缘检测算子极其实现(基于Python和OpenCV)的主要内容,如果未能解决你的问题,请参考以下文章

常见边缘检测对比(Roberts算子Prewitt算子Sobel算子Laplacian算子Canny算子)

常见边缘检测对比(Roberts算子Prewitt算子Sobel算子Laplacian算子Canny算子)

常见边缘检测对比(Roberts算子Prewitt算子Sobel算子Laplacian算子Canny算子)

计算机视觉教程3-1:全面详解图像边缘检测算法(附Python实战)

边缘检测算子中抗噪性能最好的是?

图像处理之Canny边缘检测