OpenCV計算機視覺學習（6）——圖像梯度計算&邊緣檢測（Sobel算子，Laplacian算子，Canny算子）

　　本文學習利用python學習邊緣檢測的濾波器，首先讀入的圖片代碼如下：

import cv2
from pylab import *

img = cv2.imread("construction.jpg")
img = cv2.cvtColor(img,cv2.COLOR_BGR2RGB)
plt.imshow(img)
plt.axis("off")
plt.show()

 　　圖片如下：

　　下面列一下本節學習的內容目錄（這里於2020.6.9修改）：

　　前言：邊緣檢測的定義和類型

        1，一階微分算子
            1.1  Roberts交叉梯度算子
            1.2  Prewitt算子
            1.3  Sobel算子
            1.4  Isotropic Sobel算子
            1.5  Scharr算子
            1.6  Sobel算子，Roberts算子，Prewitt算子的比較

        2，二階微分算子
            2.1  Laplacian算子
            2.2  LOG算子

        3，非微分邊緣檢測算子——Canny算子
            3.1  Canny算子邊緣檢測基本原理
            3.2  Canny算子算法步驟
            3.3  Canny算子python實現

        4，降噪后進行邊緣檢測
            4.1   邊緣檢測示例1
            4.2   邊緣檢測示例2
            4.3   邊緣檢測示例3

前言：邊緣檢測的定義和類型

　　邊緣檢測是圖像處理和計算機視覺的基本問題，邊緣檢測的目的是標識數字圖像中亮度變化明顯的點，圖像屬性中的顯著變化通常反映了屬性的重要事件和變化。這些包括：深度上的不連續，表面方向的不連續，物質屬性變化和場景照明變化。邊緣檢測是圖像處理和計算機視覺中，尤其是特征提取中的一個研究領域。圖像邊緣檢測大幅度的減少了數據量，並且剔除了可以認為不相關的信息，保留了圖像重要的結構屬性。

　　在實際的圖像分割中，往往只用到一階和二階導數，雖然原理上，可以用更高階的導數，但是因為噪聲的影響，在純粹二階的導數操作中就會出現對噪聲的敏感現象，三階以上的導數信息往往失去了應用價值。二階導數還可以說明灰度突變的類型。在某些情況下，如灰度變化均勻的圖像，只利用一階導數可能找不到邊界，此時二階導數就能提供很有用的信息。二階導數對噪聲也比較敏感，解決的方法是先對圖像進行平滑濾波，消除部分噪聲，再進行邊緣檢測。不過，利用二階導數信息的算法是基於過零檢測的，因此得到的邊緣點數比較少，有利於后繼的處理和識別工作。

　　邊緣類型：簡單分為四種類型，階躍型，屋脊型，斜坡型，脈沖型，其中階躍型和斜坡型是類似的，只是變化的快慢不同。

　　人類視覺系統認識目標的過程分為兩步:首先，把圖像邊緣與背景分離出來;然后，才能知覺到圖像的細節，辨認出圖像的輪廓。計算機視覺正是模仿人類視覺的這個過程。因此在檢測物體邊緣時，先對其輪廓點進行粗略檢測，然后通過鏈接規則把原來檢測到的輪廓點連接起來，同時也檢測和連接遺漏的邊界點及去除虛假的邊界點。圖像的邊緣是圖像的重要特征，是計算機視覺、模式識別等的基礎，因此邊緣檢測是圖象處理中一個重要的環節。然而，邊緣檢測又是圖象處理中的一個難題，由於實際景物圖像的邊緣往往是各種類型的邊緣及它們模糊化后結果的組合，且實際圖像信號存在着噪聲。噪聲和邊緣都屬於高頻信號，很難用頻帶做取舍。
　　這就需要邊緣檢測來進行解決的問題了。邊緣檢測的基本方法有很多，一階的有Roberts Cross算子，Prewitt算子，Sobel算子，Krisch算子，羅盤算子；而二階的還有Marr-Hildreth算子（又稱為LOG算子），在梯度方向的二階導數過零點，而Canny算子屬於非微分邊緣檢測算子。各種算子的存在就是對這種導數分割原理進行的實例化計算，是為了在計算過程中直接使用的一種計算單位。在對圖像的操作，我們采用模板對原圖像進行卷積運算，從而達到我們想要的效果。而獲取一幅圖像的梯度就轉化為：模板（Roberts、Prewitt、Sobel、Lapacian算子）對原圖像進行卷積。

　　所以我們可以用一幅圖來對這些算子進行比較。

一：一階微分算子

1.1 Roberts交叉梯度算子

　　Roberts 算子又稱為交叉微分算子，它是基於交叉差分的梯度算法，通過局部差分計算檢測邊緣線條。常用來處理具有陡峭的低噪聲圖像，當圖像邊緣接近於正 45 度或負 45 度時，該算法處理效果更理想。其缺點是對邊緣的定位不太准確，提取的邊緣線條較粗。

1.1.1  Roberts交叉微分算子原理：

　　Roberts交叉梯度算子的模板分為水平方向和垂直方向，由兩個2*2的模版構成，如圖：

　　從上面模板中可以看出，Roberts算子能較好的增強正負 45度 的圖像邊緣。

　　詳細計算公式如下所示：

　　對於圖像來說，是一個二維的離散型數集，通過推廣二維連續型求函數偏導的方法，來求得圖像的偏導數，即在（x, y）處的最大變換率，也就是這里的梯度：

　　梯度是一個矢量，則（x,  y）處的梯度表示為：

 　　其大小為：

 　　平方和平方根需要大量的計算開銷，所以使用絕對值來近似梯度幅值：

　　方向與 α (x,  y) 正交：

 　　對應的模板為：

　　上圖是圖像的垂直和水平梯度，但是我們有時候也需要對角線方向的梯度，定義如下：

 　　對應的模板為：

 　　2*2 大小的模板在概念上很簡單，但是他們對於用關於中心點對稱的模板來計算邊緣方向不是很有用，其最小模板大小為3*3， 3*3模板考慮了中心點對段數據的性質，並攜帶有關於邊緣方向的更多信息。

1.1.2  Roberts交叉微分算子Python實現

　　在Python中，Roberts算子主要通過 Numpy 定義模板，再調用 OpenCV的 filter2D() 函數實現邊緣提取。該函數主要是利用內核實現對圖像的卷積運算，其函數原型如下所示：

dst = filter2D(src, ddepth, kernel[, dst[, anchor[, delta[, borderType]]]])

變量解釋：
    src表示輸入圖像

    dst表示輸出的邊緣圖，其大小和通道數與輸入圖像相同

    ddepth表示目標圖像所需的深度

    kernel表示卷積核，一個單通道浮點型矩陣

    anchor表示內核的基准點，其默認值為（-1，-1），位於中心位置

    delta表示在儲存目標圖像前可選的添加到像素的值，默認值為0

    borderType表示邊框模式

　　Python代碼實現：

# -*- coding: utf-8 -*-
import cv2  
import numpy as np  
import matplotlib.pyplot as plt
 
#讀取圖像
img = cv2.imread('lena.png')
lenna_img = cv2.cvtColor(img,cv2.COLOR_BGR2RGB)

#灰度化處理圖像
grayImage = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
 
#Roberts算子
kernelx = np.array([[-1,0],[0,1]], dtype=int)
kernely = np.array([[0,-1],[1,0]], dtype=int)
x = cv2.filter2D(grayImage, cv2.CV_16S, kernelx)
y = cv2.filter2D(grayImage, cv2.CV_16S, kernely)
#轉uint8 
absX = cv2.convertScaleAbs(x)      
absY = cv2.convertScaleAbs(y)    
Roberts = cv2.addWeighted(absX,0.5,absY,0.5,0)

#用來正常顯示中文標簽
plt.rcParams['font.sans-serif']=['SimHei']

#顯示圖形
titles = [u'原始圖像', u'Roberts算子']  
images = [lenna_img, Roberts]  
for i in range(2):  
   plt.subplot(1,2,i+1), plt.imshow(images[i], 'gray')  
   plt.title(titles[i])  
   plt.xticks([]),plt.yticks([])  
plt.show()

　　不調庫的代碼實現：

import cv2
from pylab import *

saber  = cv2.imread("construction.jpg")
# 首先將原圖像進行邊界擴展，並將其轉換為灰度圖。
gray_saber = cv2.cvtColor(saber,cv2.COLOR_RGB2GRAY)
gray_saber = cv2.resize(gray_saber,(200,200))

def RobertsOperator(roi):
    operator_first = np.array([[-1,0],[0,1]])
    operator_second = np.array([[0,-1],[1,0]])
    return np.abs(np.sum(roi[1:,1:]*operator_first))+np.abs(np.sum(roi[1:,1:]*operator_second))

def RobertsAlogrithm(image):
    image = cv2.copyMakeBorder(image,1,1,1,1,cv2.BORDER_DEFAULT)
    for i in range(1,image.shape[0]):
        for j in range(1,image.shape[1]):
            image[i,j] = RobertsOperator(image[i-1:i+2,j-1:j+2])
    return image[1:image.shape[0],1:image.shape[1]]

Robert_saber = RobertsAlogrithm(gray_saber)
plt.imshow(Robert_saber,cmap="binary")
plt.axis("off")
plt.show()

 　　結果演示：

1.2  Prewitt算子

　　Prewitt 是一種圖像邊緣檢測的微分算子，其原理是利用特定區域內像素灰度值產生的差分實現邊緣檢測。由於 Prewitt 算子采用 3*3 模板對區域內的像素值進行計算，而Robert算子的模板是 2*2，故 Prewitt 算子的邊緣檢測結果在水平方向和垂直方向均比 Robert 算子更加明顯，Prewitt算子適合用來識別噪聲較多，灰度漸變的圖像。

1.2.1  Prewitt算子原理

       Prewitt算子是一種一階微分算子的邊緣檢測，利用像素點上下、左右鄰點的灰度差，在邊緣處達到極值檢測邊緣，去掉部分偽邊緣，對噪聲具有平滑作用 。其原理是在圖像空間利用兩個方向模板與圖像進行鄰域卷積來完成的，這兩個方向模板一個檢測水平邊緣，一個檢測垂直邊緣。

 　　對數字圖像f(x，y)，Prewitt算子的定義如下：

　　　　　　G(i)=|[f(i-1,j-1)+f(i-1,j)+f(i-1，j+1)]-[f(i+1,j-1)+f(i+1，j)+f(i+1，j+1)]|

　　　　　　G(j)=|[f(i-1,j+1)+f(i,j+1)+f(i+1，j+1)]-[f(i-1,j-1)+f(i,j-1)+f(i+1，j-1)]|

 

　　則 P(i,j)=max[G(i),G(j)]或 P(i,j)=G(i)+G(j)

 

　　經典Prewitt算子認為：凡灰度新值大於或等於閾值的像素點都是邊緣點。即選擇適當的閾值T，若P(i,j)≥T，則(i,j)為邊緣點，P(i,j)為邊緣圖像。這種判定是欠合理的，會造成邊緣點的誤判，因為許多噪聲點的灰度值也很大，而且對於幅值較小的邊緣點，其邊緣反而丟失了。

　　Prewitt算子對噪聲有抑制作用，抑制噪聲的原理是通過像素平均，但是像素平均相當於對圖像的低通濾波，所以Prewitt算子對邊緣的定位不如Roberts算子。

 　　計算公式為：

　　因為平均能減少或消除噪聲，Prewitt梯度算子法就是先求平均，再求差分來求梯度。水平和垂直梯度模板分別為：

　　　　檢測水平邊沿 橫向模板                          

 

　　　　檢測垂直平邊沿 縱向模板：

　　對於如圖的矩陣起始值：

　　就是以下兩個式子：

　　該算子與Sobel算子類似，只是權值有所變化，但兩者實現起來功能還是有差距的，據經驗得知Sobel要比Prewitt更能准確檢測圖像邊緣。

1.2.2  Prewitt算子Python實現

　　在Python中，Prewitt 算子的實現過程與 Roberts 算子比較相似。通過 Numpy定義模板，再調用OpenCV的 filter2D() 函數對圖像的卷積運算，最終通過  convertScaleAbs() 和 addWeighted() 函數實現邊緣提出，代碼如下所示：

# -*- coding: utf-8 -*-
import cv2  
import numpy as np  
import matplotlib.pyplot as plt
 
#讀取圖像
img = cv2.imread('lena.png')
lenna_img = cv2.cvtColor(img,cv2.COLOR_BGR2RGB)

#灰度化處理圖像
grayImage = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
 
#Prewitt算子
kernelx = np.array([[1,1,1],[0,0,0],[-1,-1,-1]],dtype=int)
kernely = np.array([[-1,0,1],[-1,0,1],[-1,0,1]],dtype=int)
x = cv2.filter2D(grayImage, cv2.CV_16S, kernelx)
y = cv2.filter2D(grayImage, cv2.CV_16S, kernely)
#轉uint8
absX = cv2.convertScaleAbs(x)       
absY = cv2.convertScaleAbs(y)    
Prewitt = cv2.addWeighted(absX,0.5,absY,0.5,0)

#用來正常顯示中文標簽
plt.rcParams['font.sans-serif']=['SimHei']

#顯示圖形
titles = [u'原始圖像', u'Prewitt算子']  
images = [lenna_img, Prewitt]  
for i in xrange(2):  
   plt.subplot(1,2,i+1), plt.imshow(images[i], 'gray')  
   plt.title(titles[i])  
   plt.xticks([]),plt.yticks([])  
plt.show()

　　不調用庫的代碼實現：

import cv2
from pylab import *

saber  = cv2.imread("construction.jpg")
# 首先將原圖像進行邊界擴展，並將其轉換為灰度圖。
gray_saber = cv2.cvtColor(saber,cv2.COLOR_RGB2GRAY)
gray_saber = cv2.resize(gray_saber,(200,200))


def PreWittOperator(roi, operator_type):
    if operator_type == "horizontal":
        prewitt_operator = np.array([[-1, -1, -1], [0, 0, 0], [1, 1, 1]])
    elif operator_type == "vertical":
        prewitt_operator = np.array([[-1, 0, 1], [-1, 0, 1], [-1, 0, 1]])
    else:
        raise ("type Error")
    result = np.abs(np.sum(roi * prewitt_operator))
    return result


def PreWittAlogrithm(image, operator_type):
    new_image = np.zeros(image.shape)
    image = cv2.copyMakeBorder(image, 1, 1, 1, 1, cv2.BORDER_DEFAULT)
    for i in range(1, image.shape[0] - 1):
        for j in range(1, image.shape[1] - 1):
            new_image[i - 1, j - 1] = PreWittOperator(image[i - 1:i + 2, j - 1:j + 2], operator_type)
    new_image = new_image * (255 / np.max(image))
    return new_image.astype(np.uint8)


plt.subplot(121)
plt.title("horizontal")
plt.imshow(PreWittAlogrithm(gray_saber,"horizontal"),cmap="binary")
plt.axis("off")
plt.subplot(122)
plt.title("vertical")
plt.imshow(PreWittAlogrithm(gray_saber,"vertical"),cmap="binary")
plt.axis("off")
plt.show()

 　　結果演示：

　　下面試一下Prewitt對噪聲的敏感性

　　代碼：

import cv2
from pylab import *

saber  = cv2.imread("construction.jpg")
# 首先將原圖像進行邊界擴展，並將其轉換為灰度圖。
gray_saber = cv2.cvtColor(saber,cv2.COLOR_RGB2GRAY)
gray_saber = cv2.resize(gray_saber,(200,200))


def PreWittOperator(roi, operator_type):
    if operator_type == "horizontal":
        prewitt_operator = np.array([[-1, -1, -1], [0, 0, 0], [1, 1, 1]])
    elif operator_type == "vertical":
        prewitt_operator = np.array([[-1, 0, 1], [-1, 0, 1], [-1, 0, 1]])
    else:
        raise ("type Error")
    result = np.abs(np.sum(roi * prewitt_operator))
    return result


def PreWittAlogrithm(image, operator_type):
    new_image = np.zeros(image.shape)
    image = cv2.copyMakeBorder(image, 1, 1, 1, 1, cv2.BORDER_DEFAULT)
    for i in range(1, image.shape[0] - 1):
        for j in range(1, image.shape[1] - 1):
            new_image[i - 1, j - 1] = PreWittOperator(image[i - 1:i + 2, j - 1:j + 2], operator_type)
    new_image = new_image * (255 / np.max(image))
    return new_image.astype(np.uint8)

def noisy(noise_typ,image):
    if noise_typ == "gauss":
        row,col,ch= image.shape
        mean = 0
        var = 0.1
        sigma = var**0.5
        gauss = np.random.normal(mean,sigma,(row,col,ch))
        gauss = gauss.reshape(row,col,ch)
        noisy = image + gauss
        return noisy
    elif noise_typ == "s&p":
        row,col,ch = image.shape
        s_vs_p = 0.5
        amount = 0.004
        out = np.copy(image)
        num_salt = np.ceil(amount * image.size * s_vs_p)
        coords = [np.random.randint(0, i - 1, int(num_salt))
              for i in image.shape]
        out[coords] = 1
        num_pepper = np.ceil(amount* image.size * (1. - s_vs_p))
        coords = [np.random.randint(0, i - 1, int(num_pepper)) for i in image.shape]
        out[coords] = 0
        return out
    elif noise_typ == "poisson":
        vals = len(np.unique(image))
        vals = 2 ** np.ceil(np.log2(vals))
        noisy = np.random.poisson(image * vals) / float(vals)
        return noisy
    elif noise_typ =="speckle":
        row,col,ch = image.shape
        gauss = np.random.randn(row,col,ch)
        gauss = gauss.reshape(row,col,ch)
        noisy = image + image * gauss
        return noisy

dst = noisy("s&p",saber)
plt.subplot(131)
plt.title("add noise")
plt.axis("off")
plt.imshow(dst)

plt.subplot(132)
plt.title("Prewitt Process horizontal")
plt.axis("off")
plt.imshow(PreWittAlogrithm(gray_saber,"horizontal"),cmap="binary")

plt.subplot(133)
plt.title("Prewitt Process vertical")
plt.axis("off")
plt.imshow(PreWittAlogrithm(gray_saber,"vertical"),cmap="binary")
plt.show()

 　　結果演示：

　　選擇水平梯度或垂直梯度從上圖可以看出對於邊緣的影響還是相當大的.

1.3  Sobel算子

　　Sobel算子是一種用於邊緣檢測的離散微分算子，它結合了高斯平滑和微分求導。該算子用於計算圖像明暗程度近似值。根據圖像邊緣旁邊明暗程度把該區域內超過某個數的特定點記為邊緣。Sobel 算子在Prewitt算子的基礎上增加了權重的概念，認為相鄰點的距離遠近對當前像素點的影響是不同的，距離越近的像素點對應當前像素的影響越大，從而實現圖像銳化並突出邊緣輪廓。

　　Sobel算子的邊緣定位更准確，常用於噪聲較多，灰度漸變的圖像。

1.3.1  Sobel算子原理

　　Sobel算子主要用於邊緣檢測，在技術上它是以離散型的差分算子，用來運算圖像亮度函數的梯度的近似值， Sobel算子是典型的基於一階導數的邊緣檢測算子，由於該算子中引入了類似局部平均的運算，因此對噪聲具有平滑作用，能很好的消除噪聲的影響。Sobel算子是在Prewitt算子的基礎上改進的，在中心系數上使用一個權值，與Prewitt算子、Roberts算子相比因此效果更好，能較好的抑制（平滑）噪聲。 
　　Sobel算子包含兩組3x3的矩陣，分別為橫向及縱向模板，將之與圖像作平面卷積，即可分別得出橫向及縱向的亮度差分近似值。

　　計算公式為：

　　實際使用中，常用如下兩個模板來檢測圖像邊緣。

　　檢測水平邊沿 橫向模板 ：        

　　檢測垂直平邊沿 縱向模板：       

　　Sobel 算子根據像素點上下，左右鄰點灰度加權差，在邊緣處達到極值這一現象檢測邊緣，對噪聲具有平滑作用，提供較為精確的邊緣方向信息。因為Sobel算子結合了高斯平滑和微分求導（分化），因此結果會具有更多的抗噪性，當對精度要求不是很高時，Sobel 算子是一種較為常用的邊緣 檢測方法。

　　圖像的每一個像素的橫向及縱向梯度近似值可用以下的公式結合，來計算梯度的大小。

　　然后可用以下公式計算梯度方向。 

　　在以上例子中，如果以上的角度Θ等於零，即代表圖像該處擁有縱向邊緣，左方較右方暗。 
　　缺點是Sobel算子並沒有將圖像的主題與背景嚴格地區分開來，換言之就是Sobel算子並沒有基於圖像灰度進行處理，由於Sobel算子並沒有嚴格地模擬人的視覺生理特征，所以提取的圖像輪廓有時並不能令人滿意。

1.3.2  Sobel算子的Python實現

　　Sobel算子依然是一種過濾器，只是其是帶有方向的。在OpenCV-Python中，使用Sobel的算子的函數原型如下：

dst = cv2.Sobel(src, ddepth, dx, dy[, dst[, ksize[, scale[, delta[, borderType]]]]])

　　參數解釋：

前四個是必須的參數：

• dst 表示輸出的邊緣圖，其大小和通道數與輸入圖像相同
• src 表示需要處理的圖像；
• ddepth 表示圖像的深度，-1表示采用的是與原圖像相同的深度。目標圖像的深度必須大於等於原圖像的深度；
• dx和dy表示的是求導的階數，dx 表示x方向上的差分階數，取值為1或者0，dy表示y方向上的差分階數，取值為1或0，0表示這個方向上沒有求導，一般為0、1。

其后是可選的參數：

• ksize是Sobel算子的大小，其值必須是正數和奇數，通常為1、3、5、7。
• scale是縮放導數的比例常數，默認情況下沒有伸縮系數；
• delta是一個可選的增量，將會加到最終的dst中，同樣，默認情況下沒有額外的值加到dst中；
• borderType是判斷圖像邊界的模式。這個參數默認值為cv2.BORDER_DEFAULT。

　　注意：在進行Sobel算子處理之后，還需要調用convertScaleAbs() 函數計算絕對值，並將圖像轉換為8位圖像進行顯示。原因是sobel算子求導的話，白到黑是正數，但是黑到白就是負數了，所有的負數都會被截斷為0，所以要取絕對值。

　　下面看一下 convertScaleAbs()函數原型：

dst = convertScaleAbs(src[, dst[, alpha[, beta]]])

    src表示原數組

    dst表示輸出數組，深度為8位

    alpha表示比例因子

    beta表示原數組元素按比例縮放后添加的值

　　在OpenCV-Python中，Sobel函數的使用如下：

#coding=utf-8
import cv2
import numpy as np  
 
img = cv2.imread("D:/lion.jpg", 0)
 
x = cv2.Sobel(img,cv2.CV_16S,1,0)
y = cv2.Sobel(img,cv2.CV_16S,0,1)
 
absX = cv2.convertScaleAbs(x)   # 轉回uint8
absY = cv2.convertScaleAbs(y)
 
dst = cv2.addWeighted(absX,0.5,absY,0.5,0)
 
cv2.imshow("absX", absX)
cv2.imshow("absY", absY)
 
cv2.imshow("Result", dst)
 
cv2.waitKey(0)
cv2.destroyAllWindows()

　　解釋：

　　在Sobel函數的第二個參數這里使用了cv2.CV_16S。因為OpenCV文檔中對Sobel算子的介紹中有這么一句：“in the case of 8-bit input images it will result in truncated derivatives”。即Sobel函數求完導數后會有負值，還有會大於255的值。而原圖像是uint8，即8位無符號數，所以Sobel建立的圖像位數不夠，會有截斷。因此要使用16位有符號的數據類型，即cv2.CV_16S。

　　在經過處理后，別忘了用convertScaleAbs()函數將其轉回原來的uint8形式。否則將無法顯示圖像，而只是一副灰色的窗口。convertScaleAbs()的原型為：

dst = cv2.convertScaleAbs(src[, dst[, alpha[, beta]]])

　　其中可選參數alpha是伸縮系數，beta是加到結果上的一個值。結果返回uint8類型的圖片。

由於Sobel算子是在兩個方向計算的，最后還需要用cv2.addWeighted(...)函數將其組合起來。其函數原型為：

dst = cv2.addWeighted(src1, alpha, src2, beta, gamma[, dst[, dtype]])

　　其中alpha是第一幅圖片中元素的權重，beta是第二個的權重，gamma是加到最后結果上的一個值。

　　下面有兩個自己實現的 addWeighted（）的方法，效果和直接使用函數一樣：

def tes(filename1, filename2):
    img1 = cv2.imread(filename1, 0)
    img2 = cv2.imread(filename2, 0)
    # (352, 642) (221, 405)
    # print(img1.shape, img2.shape)
    img1 = cv2.resize(img1, (405, 221))
    # print(img1.shape,img2.shape)
    # (221, 405) (221, 405)

    res1 = cv2.addWeighted(img1, 0.5, img2, 0.5, 0)
    res2 = img1 * 0.5 + img2 * 0.5
    res2 = np.uint8(res2)

    res3 = np.zeros(img1.shape)
    for i in range(img1.shape[0]):
        for j in range(img1.shape[1]):
            res3[i, j] = int(img1[i, j] * 0.5 + img2[i, j] * 0.5)
    res3 = np.uint8(res3)


    cv2.imshow('res1', res1)
    cv2.imshow('res2', res2)
    cv2.imshow('res3', res3)
    cv2.waitKey(0)
    cv2.destroyAllWindows()

 

　　代碼1實現Sobel算子：

import cv2
from pylab import *

saber  = cv2.imread("construction.jpg")
gray_saber = cv2.cvtColor(saber,cv2.COLOR_RGB2GRAY)
gray_saber = cv2.resize(gray_saber,(200,200))

def SobelOperator(roi, operator_type):
    if operator_type == "horizontal":
        sobel_operator = np.array([[-1, -2, -1], [0, 0, 0], [1, 2, 1]])
    elif operator_type == "vertical":
        sobel_operator = np.array([[-1, 0, 1], [-2, 0, 2], [-1, 0, 1]])
    else:
        raise ("type Error")
    result = np.abs(np.sum(roi * sobel_operator))
    return result


def SobelAlogrithm(image, operator_type):
    new_image = np.zeros(image.shape)
    image = cv2.copyMakeBorder(image, 1, 1, 1, 1, cv2.BORDER_DEFAULT)
    for i in range(1, image.shape[0] - 1):
        for j in range(1, image.shape[1] - 1):
            new_image[i - 1, j - 1] = SobelOperator(image[i - 1:i + 2, j - 1:j + 2], operator_type)
    new_image = new_image * (255 / np.max(image))
    return new_image.astype(np.uint8)

plt.subplot(121)
plt.title("horizontal")
plt.imshow(SobelAlogrithm(gray_saber,"horizontal"),cmap="binary")
plt.axis("off")
plt.subplot(122)
plt.title("vertical")
plt.imshow(SobelAlogrithm(gray_saber,"vertical"),cmap="binary")
plt.axis("off")
plt.show()

 　　結果1展示：

　　代碼2實現Sobel算子：

#_*_coding:utf-8_*_
from PIL import Image
from PIL import ImageEnhance
from numpy import *
from pylab import *
from scipy.ndimage import filters

image1 = Image.open('construction.jpg').convert('L')
im = array(image1)
#soble 導數濾波器  使用 Sobel 濾波器來計算 x 和 y 的方向導數，
imx = zeros(im.shape)
# print(imx)
filters.sobel(im,1,imx)

imy = zeros(im.shape)
filters.sobel(im,0,imy)

magnitude = sqrt(imx**2 + imy**2)
# print(magnitude)
def deal_with(a):
    for i in range(len(a)):
        if a[i] <50:
            a[i] =0
        elif a[i] >200:
            a[i] =255
        # else:
        #     a[i] = 155
    return a
a = np.apply_along_axis(deal_with,1,magnitude)
result =  contour(magnitude, origin='image')
axis('equal')
axis('off')
figure()
hist(magnitude.flatten(),128)
show()

 　　結果2展示：

1.4  Isotropic Sobel算子

　　Sobel算子另一種形式是(Isotropic Sobel)算子，加權平均算子，權值反比零點與中心店的距離，當沿不同方向檢測邊緣時梯度幅度一致，就是通常所說的各向同性Sobel(Isotropic Sobel)算子。模板也有兩個，一個是檢測水平邊沿的 ，另一個是檢測垂直平邊沿的 。各向同性Sobel算子和普通Sobel算子相比，它的位置加權系數更為准確，在檢測不同方向的邊沿時梯度的幅度一致。

1.5  Scharr算子

1.5.1  Scharr算子原理

　　Scharr算子與Sobel算子的不同點是在平滑部分，這里所用的平滑算子是 1/16 *[3, 10, 3]，相比於 1/4*[1, 2, 1]，中心元素占的權重更重，這可能是相對於圖像這種隨機性較強的信號，領域相關性不大，所以鄰域平滑應該使用相對較小的標准差的高斯函數，也就是更瘦高的模板。

　　由於Sobel算子在計算相對較小的核的時候，其近似計算導數的精度比較低，比如一個3*3的Sobel算子，當梯度角度接近水平或垂直方向時，其不精確性就越發明顯。Scharr算子同Sobel算子的速度一樣快，但是准確率更高，尤其是計算較小核的情景，所以利用3*3濾波器實現圖像邊緣提取更推薦使用Scharr算子。

　　Scharr算子又稱為Scharr濾波器，也是計算x或y方向上的圖像差分，在OpenCV中主要配合Sobel算子的運算而存在的，下面對比一下Sobel算子和scharr算子的核函數對比：

 

 1.5.2  Scharr算子Python實現

　　Scharr算子和Sobel算子類似，這里簡單說一下其函數用法。在OpenCV-Python中，使用Scharrl的算子的函數原型如下：

dst = cv2.Scharr(src, ddepth, dx, dy[, dst[, ksize[, scale[, delta[, borderType]]]]])

　　參數解釋：

前四個是必須的參數：

• dst 表示輸出的邊緣圖，其大小和通道數與輸入圖像相同
• src 表示需要處理的圖像；
• ddepth 表示圖像的深度，-1表示采用的是與原圖像相同的深度。目標圖像的深度必須大於等於原圖像的深度；
• dx和dy表示的是求導的階數，dx 表示x方向上的差分階數，取值為1或者0，dy表示y方向上的差分階數，取值為1或0，0表示這個方向上沒有求導，一般為0、1。

其后是可選的參數：

• ksize是Sobel算子的大小，其值必須是正數和奇數，通常為1、3、5、7。
• scale是縮放導數的比例常數，默認情況下沒有伸縮系數；
• delta是一個可選的增量，將會加到最終的dst中，同樣，默認情況下沒有額外的值加到dst中；
• borderType是判斷圖像邊界的模式。這個參數默認值為cv2.BORDER_DEFAULT。

　　注意：在進行Scharr算子處理之后，也需要調用convertScaleAbs() 函數計算絕對值，並將圖像轉換為8位圖像進行顯示。原因是sobel算子求導的話，白到黑是正數，但是黑到白就是負數了，所有的負數都會被截斷為0，所以要取絕對值。

　　這里不再贅述convertScaleAbs()函數了，直接看一個Scharr算子的實現，並將其與Sobel算子對比，我們看看效果。

　　代碼如下：

# coding=utf-8
import cv2
import numpy as np

img = cv2.imread("durant.jpg", 0)
img = cv2.resize(img, (0, 0), fx=0.5, fy=0.5)

sobel_x = cv2.Sobel(img, cv2.CV_16S, 1, 0)
sobel_y = cv2.Sobel(img, cv2.CV_16S, 0, 1)
scharr_x = cv2.Scharr(img, cv2.CV_16S, 1, 0)
scharr_y = cv2.Scharr(img, cv2.CV_16S, 0, 1)

sobel_absX = cv2.convertScaleAbs(sobel_x)  # 轉回uint8
sobel_absY = cv2.convertScaleAbs(sobel_y)
scharr_absX = cv2.convertScaleAbs(scharr_x)  # 轉回uint8
scharr_absY = cv2.convertScaleAbs(scharr_y)

Sobel_dst = cv2.addWeighted(sobel_absX, 0.5, sobel_absY, 0.5, 0)
Scharr_dst = cv2.addWeighted(scharr_absX, 0.5, scharr_absY, 0.5, 0)

# cv2.imshow("absX", scharr_absX)
# cv2.imshow("absY", scharr_absY)
# cv2.imshow("Result", dst)

sobel_image = np.hstack((sobel_absX, sobel_absY, Sobel_dst))
scharr_image = np.hstack((scharr_absX, scharr_absY, Scharr_dst))
all_image = np.vstack((sobel_image, scharr_image))
cv2.imshow("Result", all_image)


cv2.waitKey(0)
cv2.destroyAllWindows()

　　原圖如下（其實這張圖可以很明顯的看出求x軸和y軸方向梯度的差異）：

 　　效果如下：

 　　這里還是用我最喜歡的球星做一下，哈哈哈哈哈哈哈：

 　　上面三張圖分別是Sobel算子的x軸，y軸，和對x軸和y軸加起來的效果，下面三張圖分別是Scharr算子的x軸，y軸，和對x軸和y軸加起來的效果。但是上面最后面兩個效果圖的對比，並不能說明Scharr算子比Sobel算子效果更好，可能在具體的實際場景中才會出現那個更好。

1.6  Sobel算子，Robert算子，prewitt算子的比較

　　Sobel算子是濾波算子的形式來提取邊緣，X，Y方向各用一個模板，兩個模板組合起來構成一個梯度算子。X方向模板對垂直邊緣影響最大，Y方向模板對水平邊緣影響最大。

　　Robert算子是一種梯度算子，它用交叉的查分表示梯度，是一種利用局部差分算子尋找邊緣的算子，對具有陡峭的低噪聲的圖像效果最好。

　　prewitt算子是加權平均算子，對噪聲有抑制作用，但是像素平均相當於對圖像進行的同濾波，所以prewitt算子對邊緣的定位不如robert算子。

二：二階微分算子

　　拉普拉斯（Laplacian）算子是n維歐幾里德空間中的一個二階微分算子，常用於圖像增強領域和邊緣提取。它通過灰度差分計算鄰域內的像素，基本流程是：判斷圖像中心像素灰度值與它周圍其他像素的灰度值，如果中心像素的灰度更高，則提升中心像素的灰度；反之降低中心像素的灰度，從而實現圖像銳化操作。在算法實現過程中，Laplacian算子通過對鄰域中心像素的四方向或八方向求梯度，再將梯度相加起來判斷中心像素灰度與鄰域內其他像素灰度的關系，最后通過梯度運算的結果對像素灰度進行調整。

2.1 Laplacian算子

　　拉普拉斯（Laplacian）算子是 n 維歐幾里得空間中的一個二階微分算子，常用於圖像增強領域和邊緣提取，它通過灰度差分計算領域內的像素。

2.1.1 Laplacian原理

　　Laplacian算子的基本流程是：判斷圖像中心像素灰度值與它周圍其他像素的灰度值，如果中心像素的灰度更高，則提升中心像素的灰度；反之降低中心像素的灰度，從而實現圖像銳化操作。在算法實現過程中，Laplacian算子通過對鄰域中心像素的四方向或八方向求梯度，再將梯度相加起來判斷中心像素灰度與鄰域內其他像素灰度的關系，最后通過梯度運算的結果對像素灰度進行調整。

         Laplace算子是一種各向同性算子，二階微分算子，具有旋轉不變性。在只關心邊緣的位置而不考慮其周圍的象素灰度差值時比較合適。Laplace算子對孤立象素的響應要比對邊緣或線的響應要更強烈，因此只適用於無噪聲圖象。存在噪聲情況下，使用Laplacian算子檢測邊緣之前需要先進行低通濾波。所以，通常的分割算法都是把Laplacian算子和平滑算子結合起來生成一個新的模板。

　　一階導數為：

 　　二階導數為：

 　　我們這里需要的是關於x的二階導數，故將上式中的變量減去1后，得到：

 　　在圖像處理中，通過拉普拉斯模板求二階導數，其定義如下：

　　為了更適合於數字圖像處理，將該方程表示為離散形式：

　　　　   另外，拉普拉斯算子還可以表示成模板的形式，如下圖所示。從模板形式容易看出，如果在圖像中一個較暗的區域中出現了一個亮點，那么用拉普拉斯運算就會使這個亮點變得更亮。因為圖像中的邊緣就是那些灰度發生跳變的區域，所以拉普拉斯銳化模板在邊緣檢測中很有用。一般增強技術對於陡峭的邊緣和緩慢變化的邊緣很難確定其邊緣線的位置。但此算子卻可用二次微分正峰和負峰之間的過零點來確定，對孤立點或端點更為敏感，因此特別適用於以突出圖像中的孤立點、孤立線或線端點為目的的場合。同梯度算子一樣，拉普拉斯算子也會增強圖像中的噪聲，有時用拉普拉斯算子進行邊緣檢測時，可將圖像先進行平滑處理。

 　　Laplacian算子分為四鄰域和八鄰域，四鄰域是對鄰域中心像素的四方向求梯度，八鄰域是對八方向求梯度。其中四鄰域模板如公式所示：

　　　　離散拉普拉斯算子的模板：　　　　，

　　通過模板可以發現，當鄰域內像素灰度相同時，模板的卷積運算結果為0；當中心像素灰度高於鄰域內其他像素的平均灰度時，模板的卷積運算結果為正數；當中心像素的灰度低於鄰域內其他像素的平均灰度時，模板的卷積為負數。對卷積運算的結果用適當的衰弱因子處理並加在原中心像素上，就可以實現圖像的銳化處理。

　　Laplacian算子的八鄰域模板如下：

　　　　其擴展模板：　　　　　　　　　　 。

      拉式算子用來改善因擴散效應的模糊特別有效，因為它符合降制模型。擴散效應是成像過程中經常發生的現象。

      Laplacian算子一般不以其原始形式用於邊緣檢測，因為其作為一個二階導數，Laplacian算子對噪聲具有無法接受的敏感性；同時其幅值產生算邊緣，這是復雜的分割不希望有的結果；最后Laplacian算子不能檢測邊緣的方向；所以Laplacian在分割中所起的作用包括：（1）利用它的零交叉性質進行邊緣定位；（2）確定一個像素是在一條邊緣暗的一面還是亮的一面；一般使用的是高斯型拉普拉斯算子（Laplacian of a Gaussian,LoG)，由於二階導數是線性運算，利用LoG卷積一幅圖像與首先使用高斯型平滑函數卷積改圖像，然后計算所得結果的拉普拉斯是一樣的。所以在LoG公式中使用高斯函數的目的就是對圖像進行平滑處理，使用Laplacian算子的目的是提供一幅用零交叉確定邊緣位置的圖像；圖像的平滑處理減少了噪聲的影響並且它的主要作用還是抵消由Laplacian算子的二階導數引起的逐漸增加的噪聲影響。

　　  圖像銳化處理的作用是使灰度反差增強，從而使模糊圖像變得更加清晰。圖像模糊的實質就是圖像受到平均運算或積分運算，因此可以對圖像進行逆運算，如微分運算能夠突出圖像細節，使圖像變得更為清晰。由於拉普拉斯是一種微分算子，它的應用可增強圖像中灰度突變的區域，減弱灰度的緩慢變化區域。因此，銳化處理可選擇拉普拉斯算子對原圖像進行處理，產生描述灰度突變的圖像，再將拉普拉斯圖像與原始圖像疊加而產生銳化圖像。拉普拉斯銳化的基本方法可以由下式表示：

　　這種簡單的銳化方法既可以產生拉普拉斯銳化處理的效果，同時又能保留背景信息，將原始圖像疊加到拉普拉斯變換的處理結果中去，可以使圖像中的各灰度值得到保留，使灰度突變處的對比度得到增強，最終結果是在保留圖像背景的前提下，突現出圖像中小的細節信息。

2.1.2 Laplacian算子Python實現

   　　圖（a）顯示了一幅花朵的圖片，圖（b）顯示了用圖（a）所示的拉普拉斯模板對該圖像濾波后的結果。由圖可以看出，將原始圖像通過拉普拉斯變換后增強了圖像中灰度突變處的對比度，使圖像中小的細節部分得到增強並保留了圖像的背景色調，使圖像的細節比原始圖像更加清晰。基於拉普拉斯變換的圖像增強已成為圖像銳化處理的基本工具。

　　Python和OpenCV 將Laplacian算子封裝在Laplacian()函數中，其函數原型如下所示：

dst = Laplacian(src, ddepth[, dst[, ksize[, scale[, delta[, borderType]]]]])

    src表示輸入圖像


    dst表示輸出的邊緣圖，其大小和通道數與輸入圖像相同
    ddepth表示目標圖像所需的深度

    ksize表示用於計算二階導數的濾波器的孔徑大小，其值必須是正數和奇數，
且默認值為1，更多詳細信息查閱getDerivKernels

    scale表示計算拉普拉斯算子值的可選比例因子。默認值為1，更多詳細信息查閱getDerivKernels

    delta表示將結果存入目標圖像之前，添加到結果中的可選增量值，默認值為0

    borderType表示邊框模式，更多詳細信息查閱BorderTypes

    注意，Laplacian算子其實主要是利用Sobel算子的運算，通過加上Sobel算子運算出的
圖像x方向和y方向上的導數，得到輸入圖像的圖像銳化結果。同時，在進行Laplacian算子
處理之后，還需要調用convertScaleAbs()函數計算絕對值，並將圖像轉換為8位圖進行顯示。

　　當ksize=1時，Laplacian()函數采用3×3的孔徑（四鄰域模板）進行變換處理。下面的代碼是采用ksize=3的Laplacian算子進行圖像銳化處理，其代碼如下：

# -*- coding: utf-8 -*-
import cv2  
import numpy as np  
import matplotlib.pyplot as plt
 
#讀取圖像
img = cv2.imread('lena.png')
lenna_img = cv2.cvtColor(img,cv2.COLOR_BGR2RGB)

#灰度化處理圖像
grayImage = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
 
#拉普拉斯算法
dst = cv2.Laplacian(grayImage, cv2.CV_16S, ksize = 3)
Laplacian = cv2.convertScaleAbs(dst) 

#用來正常顯示中文標簽
plt.rcParams['font.sans-serif']=['SimHei']

#顯示圖形
titles = [u'原始圖像', u'Laplacian算子']  
images = [lenna_img, Laplacian]  
for i in range(2):  
   plt.subplot(1,2,i+1), plt.imshow(images[i], 'gray')  
   plt.title(titles[i])  
   plt.xticks([]),plt.yticks([])  
plt.show()

　　不調庫實現的代碼：

import cv2
from pylab import *

saber  = cv2.imread("construction.jpg")
# 首先將原圖像進行邊界擴展，並將其轉換為灰度圖。
gray_saber = cv2.cvtColor(saber,cv2.COLOR_RGB2GRAY)
gray_saber = cv2.resize(gray_saber,(200,200))


def LaplaceOperator(roi, operator_type):
    if operator_type == "fourfields":
        laplace_operator = np.array([[0, 1, 0], [1, -4, 1], [0, 1, 0]])
    elif operator_type == "eightfields":
        laplace_operator = np.array([[1, 1, 1], [1, -8, 1], [1, 1, 1]])
    else:
        raise ("type Error")
    result = np.abs(np.sum(roi * laplace_operator))
    return result


def LaplaceAlogrithm(image, operator_type):
    new_image = np.zeros(image.shape)
    image = cv2.copyMakeBorder(image, 1, 1, 1, 1, cv2.BORDER_DEFAULT)
    for i in range(1, image.shape[0] - 1):
        for j in range(1, image.shape[1] - 1):
            new_image[i - 1, j - 1] = LaplaceOperator(image[i - 1:i + 2, j - 1:j + 2], operator_type)
    new_image = new_image * (255 / np.max(image))
    return new_image.astype(np.uint8)

def noisy(noise_typ,image):
    if noise_typ == "gauss":
        row,col,ch= image.shape
        mean = 0
        var = 0.1
        sigma = var**0.5
        gauss = np.random.normal(mean,sigma,(row,col,ch))
        gauss = gauss.reshape(row,col,ch)
        noisy = image + gauss
        return noisy
    elif noise_typ == "s&p":
        row,col,ch = image.shape
        s_vs_p = 0.5
        amount = 0.004
        out = np.copy(image)
        num_salt = np.ceil(amount * image.size * s_vs_p)
        coords = [np.random.randint(0, i - 1, int(num_salt))
              for i in image.shape]
        out[coords] = 1
        num_pepper = np.ceil(amount* image.size * (1. - s_vs_p))
        coords = [np.random.randint(0, i - 1, int(num_pepper)) for i in image.shape]
        out[coords] = 0
        return out
    elif noise_typ == "poisson":
        vals = len(np.unique(image))
        vals = 2 ** np.ceil(np.log2(vals))
        noisy = np.random.poisson(image * vals) / float(vals)
        return noisy
    elif noise_typ =="speckle":
        row,col,ch = image.shape
        gauss = np.random.randn(row,col,ch)
        gauss = gauss.reshape(row,col,ch)
        noisy = image + image * gauss
        return noisy

plt.subplot(121)
plt.title("fourfields")
plt.imshow(LaplaceAlogrithm(gray_saber,"fourfields"),cmap="binary")
plt.axis("off")
plt.subplot(122)
plt.title("eightfields")
plt.imshow(LaplaceAlogrithm(gray_saber,"eightfields"),cmap="binary")
plt.axis("off")
plt.show()

 

　　結果演示：

　　上圖為比較取值不同的laplace算子實現的區別。

　　PIL庫實現的代碼：

import numpy as np
from PIL import Image
import matplotlib.pyplot as plt
import matplotlib.cm as cm
import scipy.signal as signal     # 導入sicpy的signal模塊

# Laplace算子
suanzi1 = np.array([[0, 1, 0],
                    [1,-4, 1],
                    [0, 1, 0]])

# Laplace擴展算子
suanzi2 = np.array([[1, 1, 1],
                    [1,-8, 1],
                    [1, 1, 1]])

# 打開圖像並轉化成灰度圖像
image = Image.open("construction.jpg").convert("L")
image_array = np.array(image)

# 利用signal的convolve計算卷積
image_suanzi1 = signal.convolve2d(image_array,suanzi1,mode="same")
image_suanzi2 = signal.convolve2d(image_array,suanzi2,mode="same")

# 將卷積結果轉化成0~255
image_suanzi1 = (image_suanzi1/float(image_suanzi1.max()))*255
image_suanzi2 = (image_suanzi2/float(image_suanzi2.max()))*255

# 為了使看清邊緣檢測結果，將大於灰度平均值的灰度變成255(白色)
image_suanzi1[image_suanzi1>image_suanzi1.mean()] = 255
image_suanzi2[image_suanzi2>image_suanzi2.mean()] = 255

# 顯示圖像
plt.subplot(2,1,1)
plt.imshow(image_array,cmap=cm.gray)
plt.axis("off")
plt.subplot(2,2,3)
plt.imshow(image_suanzi1,cmap=cm.gray)
plt.axis("off")
plt.subplot(2,2,4)
plt.imshow(image_suanzi2,cmap=cm.gray)
plt.axis("off")
plt.show()

 　　結果：

　　其中上方為原圖像，下方：左邊為Laplace算子結果，右邊為Laplace擴展算子結果

2.2  LOG算子

　　LOG（Laplacian of Gaussian）邊緣檢測算子是David Courtnay Marr和 Ellen Hildreth在 1980年共同提出，也稱為 Marr & Hildreth算子。它根據圖像的信噪比來求檢測邊緣的最優濾波器。下面學習一下其原理和應用。

2.2.1  LOG算子原理

　　LOG算子首先對圖像做高斯濾波，然后再求其拉普拉斯（Laplacian）二階導數，根據二階導數的鍋零點來檢測圖像的邊界，即通過檢測濾波結果的零交叉（Zero  crossings）來獲得圖像或物體的邊緣。

　　LOG算子綜合考慮了對噪聲的抑制和對邊緣的檢測兩個方向，並且把Gauss平滑濾波器和Laplacian銳化濾波器結合了起來，先平滑掉噪聲，再進行邊緣檢測，所以效果會更好。該算子與視覺生理中的數學模型相似，因此在圖像處理領域中得到了廣泛的應用。它具有抗干擾能力強，邊界定位精度高，邊緣連續性好，能有效提取對比度弱的邊界等特點。

　　常見的LOG算子是5*5模板，如下圖所示：

 　　由於LOG算子到中心的距離與位置加權系數的關系曲線像墨西哥草帽的剖面，所以LOG算子也叫墨西哥草帽濾波器，如下圖所示：

2.2.2  LOG算子Python實現

　　LOG算子的邊緣提取實現代碼如下：

# coding=utf-8
import cv2
import numpy as np
import matplotlib.pyplot as plt

img = cv2.imread("durant.jpg")
rgb_img = cv2.cvtColor(img, cv2.COLOR_BGR2RGB)
img = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)

# 先通過高斯濾波降噪
gaussian = cv2.GaussianBlur(img, (3, 3), 0)

# 再通過拉普拉斯算子做邊緣檢測
dst = cv2.Laplacian(gaussian, cv2.CV_16S, ksize=3)
LOG = cv2.convertScaleAbs(dst)

# 用來正常顯示中文標簽
plt.rcParams['font.sans-serif'] = ['SimHei']

# 顯示圖形
titles = [u'原始圖像', u'LOG算子']
images = [rgb_img, LOG]
for i in range(2):
    plt.subplot(1, 2, i + 1), plt.imshow(images[i], 'gray')
    plt.title(titles[i])
    plt.xticks([]), plt.yticks([])
plt.show()

 　　結果如下：

三： 非微分邊緣檢測算子——Canny算子

　　John F.Canny 於 1986年發明了一個多級邊緣檢測算法——Canny邊緣檢測算子，並創立了邊緣檢測計算理論（Computational  theory of edge detection），該理論有效的解釋了這項技術的工作理論。

3.1  Canny算子邊緣檢測基本原理

　　Canny邊緣檢測是一種比較新的邊緣檢測算子，具有很好地邊緣檢測性能，該算子功能比前面幾種都要好，但是它實現起來較為麻煩，Canny算子是一個具有濾波，增強，檢測的多階段的優化算子，在進行處理前，Canny算子先利用高斯平滑濾波器來平滑圖像以除去噪聲，Canny分割算法采用一階偏導的有限差分來計算梯度幅值和方向，在處理過程中，Canny算子還將經過一個非極大值抑制的過程，最后Canny算子還采用兩個閾值來連接邊緣（高低閾值輸出二值圖像）。

　　邊緣檢測通常是在保留原有圖像屬性的情況下，對圖像數據規模進行縮減，提取圖像邊緣輪廓的處理方式。

　　(1)圖象邊緣檢測必須滿足兩個條件：一能有效地抑制噪聲；二必須盡量精確確定邊緣的位置。

　　(2)根據對信噪比與定位乘積進行測度，得到最優化逼近算子。這就是Canny邊緣檢測算子。

　　(3)類似與Marr（LoG）邊緣檢測方法，也屬於先平滑后求導數的方法。

 

3.2  Canny算子的算法步驟

　　Canny算法是一種被廣泛應用於邊緣檢測的標准算法，其目標是找到一個最優的邊緣檢測解或找尋一幅圖像中灰度強度變換最強的位置。最優邊緣檢測主要通過低錯誤率，高定位性和最小響應三個標准進行評價。Canny算子的實現步驟如下：

• step1: 用高斯濾波器平滑圖象；
• step2: 計算圖像中每個像素點的梯度強度和方向（用一階偏導的有限差分來計算梯度的幅值和方向）；
• step3: 對梯度幅值進行非極大值抑制（Non-Maximum Suppression），以消除邊緣檢測帶來的雜散響應；
• step4: 用雙閾值算法(Double-Threshold)檢測來確定真實和潛在的邊緣，通過抑制孤立的弱邊緣最終完成邊緣檢測；

 

3.2.1   step1：使用高斯平滑函數去除噪聲

 

　　 3*3的高斯核函數：

 　　5*5 的高斯核函數：

3.2.2  step2：使用一階有限查分計算偏導數陣列P和Q

　　這里需要按照Sobel濾波器步驟計算梯度幅值和方向，尋找圖像的強度梯度。先將卷積模板分別作用 x 和 y 方向，再計算梯度幅值和方向，其公式如下所示。梯度方向一般取 0 度，45 度，90度和 135度四個方向。

 

 

3.2.3  step3： 非極大值抑制

　　通過非極大值抑制（Non - maximum  Suppression）過濾掉非邊緣像素，將模糊的邊界變得清晰。該過程保留了每隔像素點上梯度強度的極大值，過濾掉其他的值。

　　對於每個像素點，它進行如下操作：

• 1，將其梯度方向近似為以下值中的一個，包括0，45，90，135，180，225，270和315，即表示上下左右和45度方向。
• 2，比較該像素點和其梯度正負方向的像素點的梯度強度，如果該像素點梯度強度最大則保留，否則抑制（刪除，即置為零）。

 

 

 

3.2.3  step4： 用雙閾值算法檢測和連接邊緣　

　　利用雙閾值方法來確定潛在的邊界。經過非極大抑制后圖像中仍然有很多噪聲點，此時需要通過雙閾值技術處理，即設定一個閾值上界和閾值下界。圖像中的像素點如果大於閾值上界則認為必然是邊界（稱為強邊界，strong  edge），小於閾值下界則認為必然不是邊界，兩者之間的則認為是候選項（稱為弱邊界，weak edge）。

　　對非極大值抑制圖像作用兩個閾值th1和th2，兩者關系th1=0.4th2。我們把梯度值小於th1的像素的灰度值設為0，得到圖像1。然后把梯度值小於th2的像素的灰度值設為0，得到圖像2。由於圖像2的閾值較高，去除大部分噪音，但同時也損失了有用的邊緣信息。而圖像1的閾值較低，保留了較多的信息，我們可以以圖像2為基礎，以圖像1為補充來連結圖像的邊緣。

　　鏈接邊緣的具體步驟如下：

　　　　對圖像2進行掃描，當遇到一個非零灰度的像素p(x,y)時，跟蹤以p(x,y)為開始點的輪廓線，直到輪廓線的終點q(x,y)。

　　　　考察圖像1中與圖像2中q(x,y)點位置對應的點s(x,y)的8鄰近區域。如果在s(x,y)點的8鄰近區域中有非零像素s(x,y)存在，則將其包括到圖像2中，作為r(x,y)點。從r(x,y)開始，重復第一步，直到我們在圖像1和圖像2中都無法繼續為止。

　　　　當完成對包含p(x,y)的輪廓線的連結之后，將這條輪廓線標記為已經訪問。回到第一步，尋找下一條輪廓線。重復第一步、第二步、第三步，直到圖像2中找不到新輪廓線為止。

　　至此，完成canny算子的邊緣檢測。

 

3.3 Canny算子Python實現

　　在OpenCV中，canny() 函數原型如下所示：

edges = Canny(image, threshold1, threshold2[, edges[, apertureSize[, L2gradient]]])

 　　參數含義：

• image：表示輸入圖像（必須是單通道圖像）
• edges：表示輸出的邊緣圖，其大小和類型與輸入圖像相同
• threshold1：表示第一個滯后性閾值
• threshold2：表示第二個滯后性閾值
• apertureSize：表示應用Sobel算子的孔徑大小，其默認值為3
• L2gradient：表示一個計算圖像梯度幅值的標識，默認值為FALSE　　

　　　　　　true：表示使用更精確的L2范數進行計算（即兩個方向的倒數的平方和再開方）

　　　　　　false：表示使用L1范數（直接將兩個方向導數的絕對值相加）

 　　Canny算子的邊緣提取實現代碼如下：

# Canny邊緣提取
import cv2 as cv


def edge_demo(image):
    blurred = cv.GaussianBlur(image, (3, 3), 0)
    gray = cv.cvtColor(blurred, cv.COLOR_RGB2GRAY)
    # xgrad = cv.Sobel(gray, cv.CV_16SC1, 1, 0) #x方向梯度
    # ygrad = cv.Sobel(gray, cv.CV_16SC1, 0, 1) #y方向梯度
    # edge_output = cv.Canny(xgrad, ygrad, 50, 150)
    edge_output = cv.Canny(gray, 50, 150)
    cv.imshow("Canny Edge", edge_output)
    dst = cv.bitwise_and(image, image, mask=edge_output)
    cv.imshow("Color Edge", dst)


src = cv.imread('logo1.jpg')
# 設置為WINDOW_NORMAL可以任意縮放
# cv.namedWindow('input_image', cv.WINDOW_NORMAL)
cv.imshow('input_image', src)
edge_demo(src)
cv.waitKey(0)
cv.destroyAllWindows()

 　　結果如下：

　　高斯模糊代碼：

import cv2
from pylab import *

saber  = cv2.imread("construction.jpg")
# 首先將原圖像進行邊界擴展，並將其轉換為灰度圖。
gray_saber = cv2.cvtColor(saber,cv2.COLOR_RGB2GRAY)
gray_saber = cv2.resize(gray_saber,(200,200))


def GaussianOperator(roi):
    GaussianKernel = np.array([[1,2,1],[2,4,2],[1,2,1]])
    result = np.sum(roi*GaussianKernel/16)
    return result

def GaussianSmooth(image):
    new_image = np.zeros(image.shape)
    image = cv2.copyMakeBorder(image,1,1,1,1,cv2.BORDER_DEFAULT)
    for i in range(1,image.shape[0]-1):
        for j in range(1,image.shape[1]-1):
            new_image[i-1,j-1] =GaussianOperator(image[i-1:i+2,j-1:j+2])
    return new_image.astype(np.uint8)

smooth_saber = GaussianSmooth(gray_saber)
plt.subplot(121)
plt.title("Origin Image")
plt.axis("off")
plt.imshow(gray_saber,cmap="gray")
plt.subplot(122)
plt.title("GaussianSmooth Image")
plt.axis("off")
plt.imshow(smooth_saber,cmap="gray")
plt.show()

　　結果演示：

 四：降噪后進行邊緣檢測

4.1 邊緣檢測示例1

 　　為了獲得更好的邊緣檢測效果，可以先對圖像進行模糊平滑處理，目的是去除圖像中的高頻噪聲。

　　首先用標准差為5的5*5高斯算子對圖像進行平滑處理，然后利用Laplace的擴展算子對圖像進行邊緣檢測。

　　代碼：

import numpy as np
from PIL import Image
import matplotlib.pyplot as plt
import matplotlib.cm as cm
import scipy.signal as signal

# 生成高斯算子的函數
def func(x,y,sigma=1):
    return 100*(1/(2*np.pi*sigma))*np.exp(-((x-2)**2+(y-2)**2)/(2.0*sigma**2))

# 生成標准差為5的5*5高斯算子
suanzi1 = np.fromfunction(func,(5,5),sigma=5)

# Laplace擴展算子
suanzi2 = np.array([[1, 1, 1],
                    [1,-8, 1],
                    [1, 1, 1]])

# 打開圖像並轉化成灰度圖像
image = Image.open("4.JPG").convert("L")
image_array = np.array(image)

# 利用生成的高斯算子與原圖像進行卷積對圖像進行平滑處理
image_blur = signal.convolve2d(image_array, suanzi1, mode="same")

# 對平滑后的圖像進行邊緣檢測
image2 = signal.convolve2d(image_blur, suanzi2, mode="same")

# 結果轉化到0-255
image2 = (image2/float(image2.max()))*255

# 將大於灰度平均值的灰度值變成255（白色），便於觀察邊緣
image2[image2>image2.mean()] = 255

# 顯示圖像
plt.subplot(2,1,1)
plt.imshow(image_array,cmap=cm.gray)
plt.axis("off")
plt.subplot(2,1,2)
plt.imshow(image2,cmap=cm.gray)
plt.axis("off")
plt.show()

　　結果：

　　於 2020.1.18 號，再次修改，加深度邊緣檢測算子的學習。

4.2  邊緣檢測示例2

　　邊緣檢測算法主要是基於圖像強度的一階和二階導數，但是導數通常對噪聲很敏感，因此需要采用濾波器來過濾噪聲，並調用圖像增強或閾值化算法進行處理，最后再進行邊緣檢測。下面是采用高斯濾波去噪和閾值化處理之后，再進行邊緣檢測的過程，並對比了三種常見的邊緣提取算法。

 　　代碼如下：

# -*- coding: utf-8 -*-
import cv2
import numpy as np
import matplotlib.pyplot as plt

# 讀取圖像
img = cv2.imread('durant.jpg')
lenna_img = cv2.cvtColor(img, cv2.COLOR_BGR2RGB)

# 灰度化處理圖像
grayImage = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)

# 高斯濾波
gaussianBlur = cv2.GaussianBlur(grayImage, (3, 3), 0)

# 閾值處理
ret, binary = cv2.threshold(gaussianBlur, 127, 255, cv2.THRESH_BINARY)

# Scharr算子
x = cv2.Scharr(grayImage, cv2.CV_32F, 1, 0)  # X方向
y = cv2.Scharr(grayImage, cv2.CV_32F, 0, 1)  # Y方向
absX = cv2.convertScaleAbs(x)
absY = cv2.convertScaleAbs(y)
Scharr = cv2.addWeighted(absX, 0.5, absY, 0.5, 0)

# Canny算子
gaussian = cv2.GaussianBlur(grayImage, (3, 3), 0)  # 高斯濾波降噪
Canny = cv2.Canny(gaussian, 50, 150)

# LOG算子
gaussian = cv2.GaussianBlur(grayImage, (3, 3), 0)  # 先通過高斯濾波降噪
dst = cv2.Laplacian(gaussian, cv2.CV_16S, ksize=3)  # 再通過拉普拉斯算子做邊緣檢測
LOG = cv2.convertScaleAbs(dst)

# 效果圖
titles = ['Source Image', 'Gray Image', 'Binary Image',
          'Scharr Image', 'Canny Image', 'LOG Image']
images = [lenna_img, grayImage, binary, Scharr, Canny, LOG]
for i in np.arange(6):
    plt.subplot(2, 3, i + 1), plt.imshow(images[i], 'gray')
    plt.title(titles[i])
    plt.xticks([]), plt.yticks([])
plt.show()

 　　結果如下：

 4.3 邊緣檢測示例3

　　下面是采用高斯濾波去噪和閾值化處理之后，再進行邊緣檢測的過程，並對比了四種常見的邊緣提取算法。

　　代碼如下：

# -*- coding: utf-8 -*-
import cv2
import numpy as np
import matplotlib.pyplot as plt

# 讀取圖像
img = cv2.imread('durant.jpg')
lenna_img = cv2.cvtColor(img, cv2.COLOR_BGR2RGB)

# 灰度化處理圖像
grayImage = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)

# 高斯濾波
gaussianBlur = cv2.GaussianBlur(grayImage, (3, 3), 0)

# 閾值處理
ret, binary = cv2.threshold(gaussianBlur, 127, 255, cv2.THRESH_BINARY)

# Roberts算子
kernelx = np.array([[-1, 0], [0, 1]], dtype=int)
kernely = np.array([[0, -1], [1, 0]], dtype=int)
x = cv2.filter2D(binary, cv2.CV_16S, kernelx)
y = cv2.filter2D(binary, cv2.CV_16S, kernely)
absX = cv2.convertScaleAbs(x)
absY = cv2.convertScaleAbs(y)
Roberts = cv2.addWeighted(absX, 0.5, absY, 0.5, 0)

# Prewitt算子
kernelx = np.array([[1, 1, 1], [0, 0, 0], [-1, -1, -1]], dtype=int)
kernely = np.array([[-1, 0, 1], [-1, 0, 1], [-1, 0, 1]], dtype=int)
x = cv2.filter2D(binary, cv2.CV_16S, kernelx)
y = cv2.filter2D(binary, cv2.CV_16S, kernely)
absX = cv2.convertScaleAbs(x)
absY = cv2.convertScaleAbs(y)
Prewitt = cv2.addWeighted(absX, 0.5, absY, 0.5, 0)

# Sobel算子
x = cv2.Sobel(binary, cv2.CV_16S, 1, 0)
y = cv2.Sobel(binary, cv2.CV_16S, 0, 1)
absX = cv2.convertScaleAbs(x)
absY = cv2.convertScaleAbs(y)
Sobel = cv2.addWeighted(absX, 0.5, absY, 0.5, 0)

# 拉普拉斯算法
dst = cv2.Laplacian(binary, cv2.CV_16S, ksize=3)
Laplacian = cv2.convertScaleAbs(dst)

# 效果圖
titles = ['Source Image', 'Binary Image', 'Roberts Image',
          'Prewitt Image', 'Sobel Image', 'Laplacian Image']
images = [lenna_img, binary, Roberts, Prewitt, Sobel, Laplacian]
for i in np.arange(6):
    plt.subplot(2, 3, i + 1), plt.imshow(images[i], 'gray')
    plt.title(titles[i])
    plt.xticks([]), plt.yticks([])
plt.show()

 　　結果如下：

 

 

 

 

 

 

 

 

參考：https://www.cnblogs.com/cfantaisie/archive/2011/06/05/2073151.html

https://blog.csdn.net/sinat_32974931/article/details/51125516

https://www.cnblogs.com/lynsyklate/p/7881300.html

https://blog.csdn.net/Eastmount/article/details/89001702

https://blog.csdn.net/wsp_1138886114/article/details/82935839

https://blog.csdn.net/JimmyFu0055/article/details/83719438

https://blog.csdn.net/Eastmount/article/details/89056240