python數字圖像處理（三）邊緣檢測經常使用算子

時間  2019-11-15
標籤 python 數字 圖像處理 邊緣 檢測 經常 使用 算子 欄目 Python 简体版
原文   原文鏈接

在該文將介紹基本的幾種應用於邊緣檢測的濾波器,首先咱們讀入saber用來作爲示例的圖像html

#讀入圖像代碼,在此以前應當引入必要的opencv matplotlib numpy
saber = cv2.imread("saber.png")
saber = cv2.cvtColor(saber,cv2.COLOR_BGR2RGB)
plt.imshow(saber)
plt.axis("off")
plt.show()

使用上圖做爲濾波器使用的圖形python

Roberts 交叉梯度算子

Roberts交叉梯度算子由兩個\(2\times2\)的模版構成，如圖。
\[ \begin{bmatrix} -1 & 0 \\ 0 & 1 \end{bmatrix} \]算法

\[ \begin{bmatrix} 0 & -1 \\ 1 & 0 \end{bmatrix} \]編程

對如圖的矩陣分別才用兩個模版相乘，並加在一塊兒
\[ \begin{bmatrix} Z_1 & Z_2 & Z_3 \\ Z_4 & Z_5 & Z_6 \\ Z_7 & Z_8 & Z_9 \\ \end{bmatrix} \]
簡單表示即
\(|z_9-z_5|+|z_8-z_6|\)。
能夠簡單用編程實現app

首先將原圖像進行邊界擴展，並將其轉換爲灰度圖。dom

gray_saber = cv2.cvtColor(saber,cv2.COLOR_RGB2GRAY)
gray_saber = cv2.resize(gray_saber,(200,200))

在接下來的的代碼中咱們對圖像都會進行邊界擴展，是由於將圖像增長寬度爲1的圖像可以使在計算3*3的算子時不用考慮邊界狀況函數

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()

python的在這種純計算方面實在有些慢，固然這也是由於我過多的使用了for循環而沒有利用好numpy的緣由優化

咱們能夠看出Roberts算子在檢測邊緣效果已經至關不錯了，可是它對噪聲至關敏感，咱們給它加上噪聲在看看它的效果
噪聲代碼來源於stacksoverflowui

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(121)
plt.title("add noise")
plt.axis("off")
plt.imshow(dst)
plt.subplot(122)
dst = cv2.cvtColor(dst,cv2.COLOR_RGB2GRAY)
plt.title("Robert Process")
plt.axis("off")
dst = RobertsAlogrithm(dst)
plt.imshow(dst,cmap="binary")
plt.show()

這裏咱們能夠看出Robert算子對於噪聲過於敏感,在加了噪聲後,人的邊緣變得很是不清晰es5

Prewitt算子

Prewitt算子是一種\(3\times3\)模板的算子,它有兩種形式,分別表示水平和垂直的梯度.

垂直梯度
\[ \begin{bmatrix} -1 & -1 & -1 \\ 0 & 0 & 0 \\ 1 & 1 & 1 \\ \end{bmatrix} \]
水平梯度
\[ \begin{bmatrix} -1 & 0 & 1 \\ -1 & 0 & 1 \\ -1 & 0 & 1 \\ \end{bmatrix} \]

對於如圖的矩陣起始值
\[ \begin{bmatrix} Z_1 & Z_2 & Z_3 \\ Z_4 & Z_5 & Z_6 \\ Z_7 & Z_8 & Z_9 \\ \end{bmatrix} \]
就是如下兩個式子
\[ |(z_7+z_8+z_9)-(z_1+z_2+z_3)| \]

\[ |(z_1+z_4+z_7)-(z_3+z_6+z_9)| \]
1與-1更換位置對結果並不會產生影響

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)

在PreWitt算子執行中要考慮大於255的狀況,使用下面代碼

new_image = new_image*(255/np.max(image))

將其放縮到0-255之間,並使用astype將其轉換到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對噪聲的敏感性

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()

選擇水平梯度或垂直梯度從上圖能夠看出對於邊緣的影響仍是至關大的.

Sobel算子

Sobel算子根據像素點上下、左右鄰點灰度加權差，在邊緣處達到極值這一現象檢測邊緣。對噪聲具備平滑做用，提供較爲精確的邊緣方向信息，邊緣定位精度不夠高。當對精度要求不是很高時，是一種較爲經常使用的邊緣檢測方法。相比Prewitt他周邊像素對於判斷邊緣的貢獻是不一樣的.
垂直梯度
\[ \begin{bmatrix} -1 & -2 & -1 \\ 0 & 0 & 0 \\ 1 & 2 & 1 \\ \end{bmatrix} \]
水平梯度
\[ \begin{bmatrix} -1 & 0 & 1 \\ -2 & 0 & 2 \\ -1 & 0 & 1 \\ \end{bmatrix} \]

對上面的代碼稍微修改算子便可

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()

不管是Sobel仍是Prewitt實際上是求了單個方向顏色變化的極值

在某個方向上顏色的變化如圖(此圖來源於opencv官方文檔)

對其求導能夠獲得這樣一條曲線

由此咱們能夠得知邊緣能夠經過定位梯度值大於鄰域的相素的方法找到(或者推廣到大於一個閥值).
而Sobel就是利用了這必定理

Laplace算子

不管是Sobel仍是Prewitt都是求單方向上顏色變化的一階導(Sobel相比Prewitt增長了權重),那麼如何才能反映空間上的顏色變化的一個極值呢?

咱們試着在對邊緣求一次導

在一階導數的極值位置，二階導數爲0。因此咱們也能夠用這個特色來做爲檢測圖像邊緣的方法。 可是， 二階導數的0值不只僅出如今邊緣(它們也可能出如今無心義的位置),可是咱們能夠過濾掉這些點。
在實際狀況中咱們不可能用0做爲判斷條件,而是設定一個閾值,這就引出來Laplace算子
\begin{bmatrix}
0 & 1 & 0 \
1 & -4 & 1 \
0 & 1 & 0 \
\end{bmatrix}

\begin{bmatrix}
1 & 1 & 1 \
1 & -8 & 1 \
1 & 1 & 1 \
\end{bmatrix}

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)
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算子實現的區別.

Canny算子

Canny邊緣檢測算子是一種多級檢測算法。1986年由John F. Canny提出，同時提出了邊緣檢測的三大準則：

  • 低錯誤率的邊緣檢測：檢測算法應該精確地找到圖像中的儘量多的邊緣，儘量的減小漏檢和誤檢。
  • 最優定位：檢測的邊緣點應該精確地定位於邊緣的中心。
  • 圖像中的任意邊緣應該只被標記一次，同時圖像噪聲不該產生僞邊緣。

Canny邊緣檢測算法的實現較爲複雜,主要分爲如下步驟,

  1. 高斯模糊
  2. 計算梯度幅值和方向
  3. 非極大值 抑制
  4. 滯後閾值

高斯模糊

在前面的章節中已經實現過均值模糊,模糊在圖像處理中常常用來去除一些無關的細節.

均值模糊引出了另外一個思考,周邊每一個像素點的權重設爲同樣是否合適?
(周邊像素點的權重是不少算法的區別,好比Sobel和Prewitt)

所以機率中很是熟悉的正態分佈函數便被引入.正太分佈函數的密度函數就是高斯函數,使用高斯函數來給周邊像素點分配權重.

更詳細的能夠參考高斯模糊的算法

在opencv的官方文檔中它介紹了這樣一個\(5\times5\)的高斯算子

相比直接的使用密度函數來分配權重,Opencv介紹的高斯算子是各像素點乘以一個大於1的比重,最後乘以一個\(1 \frac {159}\),其實就是歸一化的密度函數.這有效避免了計算機計算的精度問題

咱們在此處直接使用\(3\times3\)的算子

咱們將其歸一化而後代入上面的代碼便可

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()

計算梯度幅值和方向

咱們知道邊緣實質是顏色變換的極值,可能出如今各個方向上,所以可否在計算時考慮多個方向的梯度,是算法效果優異的體現之一

在以前的Sobel等算子中咱們將水平和垂直方向是分開考慮的.在Canny算子中利用下面公式將多個方向都考慮到

Gx = SobelAlogrithm(smooth_saber,"horizontal")
Gy = SobelAlogrithm(smooth_saber,"vertical")
G = np.sqrt(np.square(Gx.astype(np.float64))+np.square(Gy.astype(np.float64)))
cita = np.arctan2(Gy.astype(np.float64),Gx.astype(np.float64))

cita求出了圖像每個點的梯度,在以後的非極大值抑制中會有做用

plt.imshow(G.astype(np.uint8),cmap="gray")
plt.axis("off")
plt.show()

咱們發現對於圖像使用sobel算子使得得到的邊緣較粗,下一步咱們將使圖像的線條變細

非最大值抑制

在得到梯度的方向和大小以後，應該對整幅圖想作一個掃描，出去那些非邊界上的點。對每個像素進行檢查，看這個點的梯度是否是周圍具備相同梯度方向的點中最大的。

(圖像來源於opencv官方文檔)

非極大值的算法實現很是簡單

  • 比較當前點的梯度強度和正負梯度方向點的梯度強度。
  • 若是當前點的梯度強度和同方向的其餘點的梯度強度相比較是最大，保留其值。不然抑制，即設爲0。好比當前點的方向指向正上方90°方向，那它須要和垂直方向，它的正上方和正下方的像素比較。
def NonmaximumSuppression(image,cita):
    keep = np.zeros(cita.shape)
    cita = np.abs(cv2.copyMakeBorder(cita,1,1,1,1,cv2.BORDER_DEFAULT))
    for i in range(1,cita.shape[0]-1):
        for j in range(1,cita.shape[1]-1):
            if cita[i][j]>cita[i-1][j] and cita[i][j]>cita[i+1][j]:
                keep[i-1][j-1] = 1
            elif cita[i][j]>cita[i][j+1] and cita[i][j]>cita[i][j-1]:
                keep[i-1][j-1] = 1
            elif cita[i][j]>cita[i+1][j+1] and cita[i][j]>cita[i-1][j-1]:
                keep[i-1][j-1] = 1
            elif cita[i][j]>cita[i-1][j+1] and cita[i][j]>cita[i+1][j-1]:
                keep[i-1][j-1] = 1
            else:
                keep[i-1][j-1] = 0
    return keep*image

這裏的代碼咱們先建立了一個keep矩陣,凡是在cita中梯度大於周邊的點設爲1,不然設爲0,而後與原圖相乘.(這是我想得比較巧妙的一個方法)

nms_image = NonmaximumSuppression(G,cita)
nms_image = (nms_image*(255/np.max(nms_image))).astype(np.uint8)
plt.imshow(nms_image,cmap="gray")
plt.axis("off")
plt.show()

在上面的極大值抑制中咱們使邊緣的線變得很是細,可是仍然有至關多的小點,得使用一種算法將他們去除

滯後閾值

滯後閾值的算法以下
Canny 使用了滯後閾值，滯後閾值須要兩個閾值(高閾值和低閾值):

  • 若是某一像素位置的幅值超過 高 閾值, 該像素被保留爲邊緣像素。
  • 若是某一像素位置的幅值小於 低 閾值, 該像素被排除。
  • 若是某一像素位置的幅值在兩個閾值之間,該像素僅僅在鏈接到一個高於 高 閾值的像素時被保留。

在opencv官方文檔中推薦咱們用2:1或3:1的高低閾值.

咱們使用3:1的閾值來進行滯後閾值算法.

滯後閾值的算法能夠分紅兩步首先使用最大閾值,最小閾值排除掉部分點.

MAXThreshold = np.max(nms_image)/4*3
MINThreshold = np.max(nms_image)/4
usemap = np.zeros(nms_image.shape)
high_list = []
for i in range(nms_image.shape[0]):
    for j in range(nms_image.shape[1]):
        if nms_image[i,j]>MAXThreshold:
            high_list.append((i,j))
direct = [(0,1),(1,1),(-1,1),(-1,-1),(1,-1),(1,0),(-1,0),(0,-1)]
def DFS(stepmap,start):
    route = [start]
    while route:
        now = route.pop()
        if usemap[now] == 1:
            break
        usemap[now] = 1
        for dic in direct:
            next_coodinate = (now[0]+dic[0],now[1]+dic[1])
            if not usemap[next_coodinate] and nms_image[next_coodinate]>MINThreshold \
            and next_coodinate[0]<stepmap.shape[0]-1 and next_coodinate[0]>=0 \
            and next_coodinate[1]<stepmap.shape[1]-1 and next_coodinate[1]>=0:
                route.append(next_coodinate)
for i in high_list:
    DFS(nms_image,i)
plt.imshow(nms_image*usemap,cmap="gray")
plt.axis("off")
plt.show()

能夠發現此次的效果並很差,這是因爲最大最小閾值設定不是很好致使的,咱們將上面代碼整理,將最大閾值,最小閾值改成可變參數,調參來優化效果.

Canny邊緣檢測代碼

def CannyAlogrithm(image,MINThreshold,MAXThreshold):
    image = GaussianSmooth(image)
    Gx = SobelAlogrithm(image,"horizontal")
    Gy = SobelAlogrithm(image,"vertical")
    G = np.sqrt(np.square(Gx.astype(np.float64))+np.square(Gy.astype(np.float64)))
    G = G*(255/np.max(G)).astype(np.uint8)
    cita = np.arctan2(Gy.astype(np.float64),Gx.astype(np.float64))
    nms_image = NonmaximumSuppression(G,cita)
    nms_image = (nms_image*(255/np.max(nms_image))).astype(np.uint8)
    usemap = np.zeros(nms_image.shape)
    high_list = []
    for i in range(nms_image.shape[0]):
        for j in range(nms_image.shape[1]):
            if nms_image[i,j]>MAXThreshold:
                high_list.append((i,j))
                
    direct = [(0,1),(1,1),(-1,1),(-1,-1),(1,-1),(1,0),(-1,0),(0,-1)]
    def DFS(stepmap,start):
        route = [start]
        while route:
            now = route.pop()
            if usemap[now] == 1:
                break
            usemap[now] = 1
            for dic in direct:
                next_coodinate = (now[0]+dic[0],now[1]+dic[1])
                if next_coodinate[0]<stepmap.shape[0]-1 and next_coodinate[0]>=0 \
                and next_coodinate[1]<stepmap.shape[1]-1 and next_coodinate[1]>=0 \
                and not usemap[next_coodinate] and nms_image[next_coodinate]>MINThreshold:
                    route.append(next_coodinate)
    for i in high_list:
        DFS(nms_image,i)
    return nms_image*usemap

咱們試着調整閾值,來對其餘圖片進行檢測

經過調整閾值能夠減小毛點的存在,控制線的清晰

(PS:在使用canny過程當中爲了減小python運算的時間,我將圖片進行了縮放,這也致使了效果不夠明顯,canny邊緣檢測在高分辨率狀況下效果會更好)

參考資料

Opencv文檔
stackoverflow Numpy設計的是如此出色,以致於在實現算法時我可以忽略如此之多無關的細節

相關文章
相關標籤/搜索
每日一句
    每一个你不满意的现在，都有一个你没有努力的曾经。
本站公眾號
   歡迎關注本站公眾號,獲取更多信息
聯繫我們 最近搜索 最新文章 沪ICP备13005482号-10 MyBatis教程 SQL 教程 MySQL教程 Java 教程 Thymeleaf 教程 Hibernate教程 Spring教程 Redis教程