影象邊緣檢測——幾種影象邊緣檢測運算元的學習及python 實現
本文學習利用python學習邊緣檢測的濾波器，首先讀入的圖片程式碼如下：
import cv2
from pylab import *
saber = cv2.imread("construction.jpg")
saber = cv2.cvtColor(saber,cv2.COLOR_BGR2RGB)
plt.imshow(saber)
plt.axis("off")
plt.show()
圖片如下：
邊緣檢測是影象處理和計算機視覺的基本問題，邊緣檢測的目的是標識數字影象中亮度變化明顯的點，影象屬性中的顯著變化通常反映了屬性的重要事件和變化。這些包括：深度上的不連續，表面方向的不連續，物質屬性變化和場景照明變化。邊緣檢測是影象處理和計算機視覺中，尤其是特徵提取中的一個研究領域。影象邊緣檢測大幅度的減少了資料量，並且剔除了可以認為不相關的資訊，保留了影象重要的結構屬性。
在實際的影象分割中，往往只用到一階和二階導數，雖然原理上，可以用更高階的導數，但是因為噪聲的影響，在純粹二階的導數操作中就會出現對噪聲的敏感現象，三階以上的導數資訊往往失去了應用價值。二階導數還可以說明灰度突變的型別。在某些情況下，如灰度變化均勻的影象，只利用一階導數可能找不到邊界，此時二階導數就能提供很有用的資訊。二階導數對噪聲也比較敏感，解決的方法是先對影象進行平滑濾波，消除部分噪聲，再進行邊緣檢測。不過，利用二階導數資訊的演算法是基於過零檢測的，因此得到的邊緣點數比較少，有利於後繼的處理和識別工作。
人類視覺系統認識目標的過程分為兩步:首先，把影象邊緣與背景分離出來;然後，才能知覺到影象的細節，辨認出影象的輪廓。計算機視覺正是模仿人類視覺的這個過程。因此在檢測物體邊緣時，先對其輪廓點進行粗略檢測，然後通過連結規則把原來檢測到的輪廓點連線起來，同時也檢測和連線遺漏的邊界點及去除虛假的邊界點。影象的邊緣是影象的重要特徵，是計算機視覺、模式識別等的基礎，因此邊緣檢測是圖象處理中一個重要的環節。然而，邊緣檢測又是圖象處理中的一個難題，由於實際景物影象的邊緣往往是各種型別的邊緣及它們模糊化後結果的組合，且實際影象訊號存在著噪聲。噪聲和邊緣都屬於高頻訊號，很難用頻帶做取捨。 這就需要邊緣檢測來進行解決的問題了。邊緣檢測的基本方法有很多，一階的有Roberts Cross運算元，Prewitt運算元，Sobel運算元，Canny運算元，Krisch運算元，羅盤運算元；而二階的還有Marr-Hildreth，在梯度方向的二階導數過零點。各種運算元的存在就是對這種導數分割原理進行的例項化計算，是為了在計算過程中直接使用的一種計算單位。在對影象的操作，我們採用模板對原影象進行卷積運算，從而達到我們想要的效果。而獲取一幅影象的梯度就轉化為：模板（Roberts、Prewitt、Sobel、Lapacian運算元）對原影象進行卷積。
一：微分運算元
1.1 Sobel運算元
理論如下：
其主要用於邊緣檢測，在技術上它是以離散型的差分運算元，用來運算影象亮度函式的梯度的近似值， Sobel運算元是典型的基於一階導數的邊緣檢測運算元，由於該運算元中引入了類似區域性平均的運算，因此對噪聲具有平滑作用，能很好的消除噪聲的影響。Sobel運算元對於象素的位置的影響做了加權，與Prewitt運算元、Roberts運算元相比因此效果更好。 Sobel運算元包含兩組3x3的矩陣，分別為橫向及縱向模板，將之與影象作平面卷積，即可分別得出橫向及縱向的亮度差分近似值。實際使用中，常用如下兩個模板來檢測影象邊緣。
檢測水平邊沿 橫向模板 ： 
檢測垂直平邊沿 縱向模板： 
影象的每一個畫素的橫向及縱向梯度近似值可用以下的公式結合，來計算梯度的大小。
然後可用以下公式計算梯度方向。
在以上例子中，如果以上的角度Θ等於零，即代表影象該處擁有縱向邊緣，左方較右方暗。 缺點是Sobel運算元並沒有將影象的主題與背景嚴格地區分開來，換言之就是Sobel運算元並沒有基於影象灰度進行處理，由於Sobel運算元並沒有嚴格地模擬人的視覺生理特徵，所以提取的影象輪廓有時並不能令人滿意。
程式碼及其結果演示：
程式碼1：
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：
#_*_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.2 Isotropic Sobel運算元
Sobel運算元另一種形式是(Isotropic Sobel)運算元，加權平均運算元，權值反比零點與中心店的距離，當沿不同方向檢測邊緣時梯度幅度一致，就是通常所說的各向同性Sobel(Isotropic Sobel)運算元。模板也有兩個，一個是檢測水平邊沿的 ，另一個是檢測垂直平邊沿的 。各向同性Sobel運算元和普通Sobel運算元相比，它的位置加權係數更為準確，在檢測不同方向的邊沿時梯度的幅度一致。
1.3 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更能準確檢測影象邊緣。
程式碼及結果演示：
程式碼：
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.4 Roberts交叉梯度運算元
理論：
Roberts交叉梯度運算元由兩個2×22×2的模版構成，如圖：
對如圖的矩陣分別才用兩個模版相乘，並加在一起
簡單表示即
程式碼及結果演示：
程式碼：
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.5 Sobel運算元，Robert運算元，prewitt運算元的比較
Sobel運算元是濾波運算元的形式來提取邊緣，X，Y方向各用一個模板，兩個模板組合起來構成一個梯度運算元。X方向模板對垂直邊緣影響最大，Y方向模板對水平邊緣影響最大。
Robert運算元是一種梯度運算元，它用交叉的查分表示梯度，是一種利用區域性差分運算元尋找邊緣的運算元，對具有陡峭的低噪聲的影象效果最好。
prewitt運算元是加權平均運算元，對噪聲有抑制作用，但是畫素平均相當於對影象進行的同濾波，所以prewitt運算元對邊緣的定位不如robert運算元。
二：Laplacian運算元
2.1 Laplace運算元理論
Laplace運算元是一種各向同性運算元，二階微分運算元，具有旋轉不變性。在只關心邊緣的位置而不考慮其周圍的象素灰度差值時比較合適。Laplace運算元對孤立象素的響應要比對邊緣或線的響應要更強烈，因此只適用於無噪聲圖象。存在噪聲情況下，使用Laplacian運算元檢測邊緣之前需要先進行低通濾波。所以，通常的分割演算法都是把Laplacian運算元和平滑運算元結合起來生成一個新的模板。
一個二維影象函式的拉普拉斯變換是各向同性的二階導數，定義為：
為了更適合於數字影象處理，將該方程表示為離散形式：
另外，拉普拉斯運算元還可以表示成模板的形式，如下圖所示。從模板形式容易看出，如果在影象中一個較暗的區域中出現了一個亮點，那麼用拉普拉斯運算就會使這個亮點變得更亮。因為影象中的邊緣就是那些灰度發生跳變的區域，所以拉普拉斯銳化模板在邊緣檢測中很有用。一般增強技術對於陡峭的邊緣和緩慢變化的邊緣很難確定其邊緣線的位置。但此運算元卻可用二次微分正峰和負峰之間的過零點來確定，對孤立點或端點更為敏感，因此特別適用於以突出影象中的孤立點、孤立線或線端點為目的的場合。同梯度運算元一樣，拉普拉斯運算元也會增強影象中的噪聲，有時用拉普拉斯運算元進行邊緣檢測時，可將影象先進行平滑處理。
離散拉普拉斯運算元的模板： ，
其擴充套件模板： 。
拉式運算元用來改善因擴散效應的模糊特別有效，因為它符合降制模型。擴散效應是成像過程中經常發生的現象。
Laplacian運算元一般不以其原始形式用於邊緣檢測，因為其作為一個二階導數，Laplacian運算元對噪聲具有無法接受的敏感性；同時其幅值產生算邊緣，這是複雜的分割不希望有的結果；最後Laplacian運算元不能檢測邊緣的方向；所以Laplacian在分割中所起的作用包括：（1）利用它的零交叉性質進行邊緣定位；（2）確定一個畫素是在一條邊緣暗的一面還是亮的一面；一般使用的是高斯型拉普拉斯運算元（Laplacian of a Gaussian,LoG)，由於二階導數是線性運算，利用LoG卷積一幅影象與首先使用高斯型平滑函式卷積改影象，然後計算所得結果的拉普拉斯是一樣的。所以在LoG公式中使用高斯函式的目的就是對影象進行平滑處理，使用Laplacian運算元的目的是提供一幅用零交叉確定邊緣位置的影象；影象的平滑處理減少了噪聲的影響並且它的主要作用還是抵消由Laplacian運算元的二階導數引起的逐漸增加的噪聲影響。
影象銳化處理的作用是使灰度反差增強，從而使模糊影象變得更加清晰。影象模糊的實質就是影象受到平均運算或積分運算，因此可以對影象進行逆運算，如微分運算能夠突出影象細節，使影象變得更為清晰。由於拉普拉斯是一種微分運算元，它的應用可增強影象中灰度突變的區域，減弱灰度的緩慢變化區域。因此，銳化處理可選擇拉普拉斯運算元對原影象進行處理，產生描述灰度突變的影象，再將拉普拉斯影象與原始影象疊加而產生銳化影象。拉普拉斯銳化的基本方法可以由下式表示：
這種簡單的銳化方法既可以產生拉普拉斯銳化處理的效果，同時又能保留背景資訊，將原始影象疊加到拉普拉斯變換的處理結果中去，可以使影象中的各灰度值得到保留，使灰度突變處的對比度得到增強，最終結果是在保留影象背景的前提下，突現出影象中小的細節資訊。
2.2 實驗結果與分析
圖（a）顯示了一幅花朵的圖片，圖（b）顯示了用圖（a）所示的拉普拉斯模板對該影象濾波後的結果。由圖可以看出，將原始影象通過拉普拉斯變換後增強了影象中灰度突變處的對比度，使影象中小的細節部分得到增強並保留了影象的背景色調，使影象的細節比原始影象更加清晰。基於拉普拉斯變換的影象增強已成為影象銳化處理的基本工具。
程式碼：
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運算元實現的區別.
程式碼二：
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擴充套件運算元結果
三： Canny運算元
Canny邊緣檢測是一種比較新的邊緣檢測運算元，具有很好地邊緣檢測效能，該運算元功能比前面幾種都要好，但是它實現起來較為麻煩，Canny運算元是一個具有濾波，增強，檢測的多階段的優化運算元，在進行處理前，Canny運算元先利用高斯平滑濾波器來平滑影象以除去噪聲，Canny分割演算法採用一階偏導的有限差分來計算梯度幅值和方向，在處理過程中，Canny運算元還將經過一個非極大值抑制的過程，最後Canny運算元還採用兩個閾值來連線邊緣。
3.1 Canny邊緣檢測基本原理
(1)圖象邊緣檢測必須滿足兩個條件：一能有效地抑制噪聲；二必須儘量精確確定邊緣的位置。 (2)根據對信噪比與定位乘積進行測度，得到最優化逼近運算元。這就是Canny邊緣檢測運算元。 (3)類似與Marr（LoG）邊緣檢測方法，也屬於先平滑後求導數的方法。
3.2 Canny邊緣檢測演算法
step1: 用高斯濾波器平滑圖象； step2: 用一階偏導的有限差分來計算梯度的幅值和方向； step3: 對梯度幅值進行非極大值抑制 step4: 用雙閾值演算法檢測和連線邊緣
step1 高斯平滑函式
step2 使用一階有限查分計算偏導數陣列P和Q
step3 非極大值抑制
step4 用雙閾值演算法檢測和連線邊緣
對非極大值抑制影象作用兩個閾值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 程式碼及其結果演示
程式碼：
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()
結果演示：
四：降噪後進行邊緣檢測
為了獲得更好的邊緣檢測效果，可以先對影象進行模糊平滑處理，目的是去除影象中的高頻噪聲。
首先用標準差為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()
結果：
參考：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