常見邊緣檢測算子極其實現(基於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

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

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

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

可以看出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

非最大值抑制：

核心：梯度是有方向矢量，通常在進行非最大抑制時使用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

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

雙閾值邊界提取

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

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