查看“︁索伯算子”︁的源代码

{{DISPLAYTITLE:索伯算子}}
{{FeatureDetectionCompVisNavbox}}
'''索伯算子'''（Sobel operator）是[[圖像處理]]中的[[算子]]之一，有時又稱為'''索伯-費爾德曼算子'''或'''索伯濾波器'''，在[[影像處理]]及[[電腦視覺領域]]中常被用來做[[邊緣檢測]]。索伯算子最早是由美國計算機科學家[[艾爾文·索伯]]及蓋瑞·費爾德曼（Gary Feldman）於1968年在[[史丹佛大學]]的人工智慧實驗室（SAIL）所提出，因此為了表揚他們的貢獻，而用他們的名字命名。在技術上，它是一離散性[[差分]]算子，用來運算圖像亮度函數的梯度之近似值。在圖像的任何一點使用此算子，索伯算子的運算將會產生對應的梯度向量或是其範數。概念上，索伯算子就是一個小且是整數的濾波器對整張影像在水平及垂直方向上做捲積，因此它所需的運算資源相對較少，另一方面，對於影像中的頻率變化較高的地方，它所得的梯度之近似值也比較粗糙。

==公式==
該算子包含兩組3x3的矩陣，分別為橫向及縱向，將之與圖像作平面[[卷積]]，即按所示方程式計算：

<math>
\begin{align}
& \left( 
\begin{bmatrix}
a & b & c \\
d & e & f \\
g & h & i
\end{bmatrix}
*
\begin{bmatrix}
1 & 2 & 3 \\
4 & 5 & 6 \\
7 & 8 & 9
\end{bmatrix} \right) [2,2]
=  \\
& (i \cdot 1)+(h \cdot 2)+(g \cdot 3)+ \\
& (f \cdot 4)+(e \cdot 5)+(d \cdot 6)+ \\
& (c \cdot 7)+(b \cdot 8)+(a \cdot 9)
\end{align}
</math>

即可分別得出橫向及縱向的亮度差分近似值。如果以<math>\mathbf{A}</math>代表原始圖像，<math>\mathbf{G_x}</math>及<math>\mathbf{G_y}</math>分別代表經橫向及縱向邊緣檢測的圖像，其公式如下：

:<math>
\mathbf{G_x} = \begin{bmatrix} 
+1 & 0 & -1 \\
+2 & 0 & -2 \\
+1 & 0 & -1
\end{bmatrix} * \mathbf{A}
\quad \mbox{and} \quad 
\mathbf{G_y} = \begin{bmatrix} 
+1 & +2 & +1 \\
0 & 0 & 0 \\
-1 & -2 & -1 
\end{bmatrix} * \mathbf{A}
</math>

圖像的每一個像素的橫向及縱向梯度近似值可用以下的公式結合，來計算梯度的大小。
:<math>\mathbf{G} = \sqrt{ \mathbf{G_x}^2 + \mathbf{G_y}^2 }</math>

然後可用以下公式計算梯度方向。
:<math>\mathbf{\Theta} = \operatorname{arctan}\left({ \mathbf{G_y} \over \mathbf{G_x} }\right)</math>

以縱向邊緣為例，如果角度<math>\Theta</math>等於零，即代表圖像的縱向邊緣右方较亮；若为 <math>\pi</math>，则左方较亮。

==更正式來說==
雖然索伯算子代表著相對沒那麼精準的影像梯度近似，但其品質也足夠在很多應用派上用場。更精確來說，對於每個點，它只用其亮度的值在周圍3x3的區域去近似該點的梯度，而且在近似的過程，它也只用整數去代表影像亮度的權重。

==索伯算子在其他維度的延伸==
索伯算子包含著兩個分開的運算:

* <span class="">利用三角形函數濾波器來達到平滑化</span> : <math> h(-1) = 1,  h(0) = 2,  h(1) = 1 </math>
* <span class="">一次微分的中央差值</span> : <math> h'(-1) = 1,  h'(0) = 0,  h'(1) = -1 </math>

索伯算子在不同維度的表達形式，<math>x,y,z,t \in \left\{0, -1, 1\right\} </math>:

1D: <math> h_x'(x) = h'(x); </math>

2D: <math> h_x'(x,y) = h'(x)h(y) </math>

3D: <math> h_x'(x,y,z) = h'(x)h(y)h(z)</math>

4D: <math> h_x'(x,y,z,t) = h'(x)h(y)h(z)h(t)</math>

下面是一個三維的索伯算子在z軸的實例:

: <math>
h_z'(:,:,-1) = 
\begin{bmatrix} 
+1 & +2 & +1  \\
+2 & +4 & +2 \\
+1 & +2 & +1 
\end{bmatrix}
\quad
h_z'(:,:,0) = 
\begin{bmatrix} 
0 & 0 & 0  \\
0 & 0 & 0 \\
0 & 0 & 0 
\end{bmatrix}
\quad
h_z'(:,:,1) = 
\begin{bmatrix} 
-1 & -2 & -1  \\
-2 & -4 & -2 \\
-1 & -2 & -1 
\end{bmatrix}

</math>

==技術細節==
索伯算子在實作上可以用很簡單的硬體和軟體: 只有該點的周圍八個點需要被計算，同時在梯度向量的近似上也都只牽涉到整數的計算，而且上述的兩個離散濾波器是可分解的:

: <math>\begin{bmatrix} 
1 & 0 & -1 \\
2 & 0 & -2 \\
1 & 0 & -1 
\end{bmatrix} = \begin{bmatrix} 
1 \\
2 \\
1  
\end{bmatrix} \begin{bmatrix} 
1 & 0 & -1
\end{bmatrix} = \begin{bmatrix} 
1 \\
1  
\end{bmatrix} * \begin{bmatrix} 
1 \\
1  
\end{bmatrix} \begin{bmatrix} 
1 & -1
\end{bmatrix} * \begin{bmatrix} 
1 & 1
\end{bmatrix}</math>
: <math>\begin{bmatrix} 
\ \ 1 & \ \ 2 & \ \ 1 \\
\ \ 0 & \ \ 0 & \ \ 0 \\
-1 & -2 & -1 
\end{bmatrix} = \begin{bmatrix} 
\ \ 1 \\
\ \ 0 \\
-1  
\end{bmatrix} \begin{bmatrix} 
1 & 2 & 1
\end{bmatrix} = \begin{bmatrix} 
1 \\
1  
\end{bmatrix} * \begin{bmatrix} 
\ \ 1 \\
-1  
\end{bmatrix} \begin{bmatrix} 
1 & 1
\end{bmatrix} * \begin{bmatrix} 
1 & 1
\end{bmatrix}
</math>

因此 '''G'''<sub>''x''</sub> 和 '''G'''<sub>''y''</sub> 可推導成

: <math>
\mathbf{G}_x = \begin{bmatrix} 
1 \\
2 \\
1
\end{bmatrix} * \left ( \begin{bmatrix} 
1 & 0 & -1  
\end{bmatrix} * \mathbf{A} \right )
\quad \mbox{and} \quad 
\mathbf{G}_y = \begin{bmatrix} 
\ \ 1 \\
\ \ 0 \\
-1  
\end{bmatrix} * \left ( \begin{bmatrix} 
1 & 2 & 1
\end{bmatrix} * \mathbf{A} \right )
</math>

在特別的實作上，可分解的運算可能會有優勢因為在每個點的運算上會有少一些的算術運算。

應用捲積核K在像素群P上的虛擬碼可以表示成:

N(x,y) =  Sum of { K(i,j).P(x-i,y-j)}, for i,j running from -1 to 1.

P代表著像素矩陣，所以N(x,y) 代表的應用捲積核K在像素群P後產生的矩陣

== 實例 ==
經過索伯算子所產生的結果是一張二維的圖，代表的每個點的梯度。下圖中高梯度的地方會以白色的線表示。為了說明這一點，下方展示了在一张簡單的圖片上使用索伯算子的結果。
{|
|[[File:Bikesgray.jpg|缩略图|200x200像素|測試的灰階影像]]
|[[File:Bikesgraysobel.jpg|缩略图|200x200像素|經過標準化後的索伯算子(梯度)結果]]
|-
|[[File:Bikesgraygv.jpg|缩略图|200x200像素|經過標準化後的x軸方向梯度結果]]
|[[File:Bikesgraygh.jpg|缩略图|200x200像素|經過標準化後的y軸方向梯度結果]]
|}
下方的灰階影像說明了梯度的方向。紅色和黃色代表梯度角度是正的，相對地，藍色和青色代表梯度角度是負的。在圓的最左邊和最右邊，垂直的梯度角度為零是因為沒有區域的變化，而在圓的最頂部和底部，水平的梯度角度為零是因為沒有區域的變化。在圓的最頂部和底部，水平梯度角度變化分別是 −π/2 和π/2是因為沒有區域的變化。頂部的負梯度角度代表這條邊是從亮到暗的區域界線，而底部的正梯度角度代表這條邊是從暗到亮的區域界線。其他顯示為黑色的像素是因為附近完全沒有任何變化。值得注意的是，因為角度是對應到像素值的函數，因此可能一些微小的變化也可能會導致很大的角度變化。因此，些微的雜訊可能會對角度響應變化，而這通常是不樂見的情況。由上述可知，使用角度梯度在影像處理的應用時，需要在影像去噪多加著墨以減少錯誤。

{|
|[[File:Black_Circle.jpg|缩略图|150x150像素|黑色圓圈與白底的灰階影像]]
|[[File:Sobel_Operator_Gradient_Angle.JPG|缩略图|180x180像素|經過索伯算子得到的梯度方向結果]]
|}

== 替代的運算子 ==
索伯算子在轉動上沒有完美的對稱。因此沙爾想要改進這個特性。它曾提出了更大的5x5核心，不過後來最常被使用的是:

: <math>
\begin{bmatrix} 
+3 & 0 & -3 \\
+10 & 0 & -10 \\
+3 & 0 & -3 
\end{bmatrix}
\ \ \ \ \ \ \ \ \  
\begin{bmatrix} 
+3 & +10 & +3 \\
 0 &   0 &  0 \\
-3 & -10 & -3 
\end{bmatrix}
</math>

沙爾算子的概念源自於試圖在頻域上最小化加權的均方角度差，這個最佳化問題需要考量濾波器在數值上的一致性，因此它們是微分核。 

另一個相似的改進方法是被法立德和西蒙切利所提出<ref>H. Farid and E. P. Simoncelli, [http://www.cs.dartmouth.edu/~farid/downloads/publications/caip97.pdf ''Optimally Rotation-Equivariant Directional Derivative Kernels''] {{Wayback|url=http://www.cs.dartmouth.edu/~farid/downloads/publications/caip97.pdf |date=20150923211417 }}, Int'l Conf Computer Analysis of Images and Patterns, pp. 207–214, Sep 1997.</ref><ref>H. Farid and E. P. Simoncelli, [http://www.cns.nyu.edu/pub/lcv/farid03-reprint.pdf ''Differentiation of discrete multi-dimensional signals''] {{Wayback|url=http://www.cns.nyu.edu/pub/lcv/farid03-reprint.pdf |date=20191116012714 }}, IEEE Trans Image Processing, vol.13(4), pp. 496–508, Apr 2004.</ref>，它們的研究是針對更高次方的微分，相較於沙爾算子，這些濾波器沒有考量到數值的一致性。

除此之外，克龍也提出了利用數值方法的最佳化來設計微分濾波器<ref>D. Kroon, 2009, Short Paper University Twente, [http://www.k-zone.nl/Kroon_DerivativePaper.pdf ''Numerical Optimization of Kernel  Based Image Derivatives''] {{Wayback|url=http://www.k-zone.nl/Kroon_DerivativePaper.pdf |date=20180705003712 }} .</ref>。

2014年，哈斯特也利用任意的三次樣條曲線去設計微分濾波器<ref>A. Hast., [http://www.sciencedirect.com/science/article/pii/S0167865514000282 "Simple filter design for first and second order derivatives by a double filtering approach"] {{Wayback|url=http://www.sciencedirect.com/science/article/pii/S0167865514000282 |date=20150924164047 }}, Pattern Recognition Letters, Vol. 42, no.1 June, pp. 65–71. 2014.</ref>，他的研究指出可以利用長度為7的濾波器作雙重濾波在三次樣條曲線及三角樣條上已取得精準的一次微分和二次微分。

另外一個相似的算子是卡雅利算子<ref>{{Cite journal|title=Alternative Approach for Satellite Cloud Classification: Edge Gradient Application|url=https://www.hindawi.com/journals/amete/2013/584816/|last=Dim|first=Jules R.|last2=Takamura|first2=Tamio|date=2013-12-11|journal=Advances in Meteorology|doi=10.1155/2013/584816|volume=2013|pages=1–8|language=en|issn=1687-9309|access-date=2018-07-04|archive-date=2019-08-20|archive-url=https://web.archive.org/web/20190820161509/https://www.hindawi.com/journals/amete/2013/584816/|dead-url=no}}</ref>，概念也是來自於索伯算子，但它是個具有完美的旋轉對稱的3x3捲積濾波器。

在那之後的工作，沙爾又針對了光流的預測上提出了更好的濾波器<ref>Scharr, Hanno, [http://www.springerlink.com/content/w527h8j3348w5075/ ''Optimal Filters for Extended Optical Flow '']{{Dead link|date=2020年3月 |bot=InternetArchiveBot |fix-attempted=yes }} In: Jähne, B., Mester, R., Barth, E., Scharr, H. (eds.) IWCM 2004. LNCS, vol. 3417, pp. 14–29. Springer, Heidelberg (2007).</ref>，同一年，他也提出了更好的二次微分濾波器應用在動作預測上<ref>Scharr, Hanno, [http://www.eurasip.org/Proceedings/Eusipco/Eusipco2007/Papers/a2l-g03.pdf ''OPTIMAL SECOND ORDER DERIVATIVE FILTER FAMILIES FOR TRANSPARENT MOTION ESTIMATION''] {{Wayback|url=http://www.eurasip.org/Proceedings/Eusipco/Eusipco2007/Papers/a2l-g03.pdf |date=20210118213146 }} 15th European Signal Processing Conference (EUSIPCO 2007), Poznan, Poland, September 3–7, 2007.</ref>。結果顯示使用的核心越大，就越有可能去逼近高斯濾波的微分。

== 結果比較 ==
下圖是四張不同的梯度算子所預測出的測試圖的梯度強度

{|
|[[File:Bikesgray.jpg|缩略图|200x200像素|灰階原圖]]
|[[File:Bikesgray_sobel.JPG|缩略图|200x200像素|索伯算子的結果]]
|[[File:Bikesgray-scharr.png|缩略图|200x200像素|沙爾算子的結果]]
|}
{|
|[[File:Bikesgray_roberts.JPG|缩略图|200x200像素|羅伯茨克羅斯算子的結果]]
|[[File:Bikesgray_prewitt.JPG|缩略图|200x200像素|普魯伊特算子的結果]]
|}

== 伪代码 ==
 function sobel(A : as two dimensional image array)
 	Gx=[-1 0 1; -2 0 2; -1 0 1]
 	Gy=[-1 -2 -1; 0 0 0; 1 2 1]
 	
 	rows = size(A,1)
 	columns = size(A,2)
 	mag=zeros(A)
 
 	for i=1:rows-2
 		for j=1:columns-2
 			S1=sum(sum(Gx.*A(i:i+2,j:j+2)))
 			S2=sum(sum(Gy.*A(i:i+2,j:j+2)))
 
 			mag(i+1,j+1)=sqrt(S1.^2+S2.^2)
 		end for
 	end for
 	
 	threshold = 70 %varies for application [0 255]
 	output_image = max(mag,threshold)
 	output_image(output_image==round(threshold))=0;
 	return output_image
 end function

== 参考文献 ==
{{Reflist}}

== 参见 ==
{{Portal box|计算机科学}}
* [[數字圖像處理]]
* [[邊緣檢測]]
* [[圖像梯度]]

[[Category:特征检测_(计算机视觉)]]