目錄html
最近在學習ORB的過程當中又仔細學習了Canny,故寫下此篇筆記,以做總結。java
Canny邊緣檢測於1986年由JOHN CANNY首次在論文《A Computational Approach to Edge Detection》中提出,就此拉開了Canny邊緣檢測算法的序幕。算法
Canny邊緣檢測是從不一樣視覺對象中提取有用的結構信息並大大減小要處理的數據量的一種技術,目前已普遍應用於各類計算機視覺系統。Canny發現,在不一樣視覺系統上對邊緣檢測的要求較爲相似,所以,能夠實現一種具備普遍應用意義的邊緣檢測技術。邊緣檢測的通常標準包括:多線程
爲了知足這些要求,Canny使用了變分法。Canny檢測器中的最優函數使用四個指數項的和來描述,它能夠由高斯函數的一階導數來近似。函數
在目前經常使用的邊緣檢測方法中,Canny邊緣檢測算法是具備嚴格定義的,能夠提供良好可靠檢測的方法之一。因爲它具備知足邊緣檢測的三個標準和實現過程簡單的優點,成爲邊緣檢測最流行的算法之一。性能
Canny邊緣檢測算法能夠分爲如下5個步驟:學習
下面詳細介紹每一步的實現思路。優化
高斯濾波是一種線性平滑濾波,適用於消除高斯噪聲,特別是對抑制或消除服從正態分佈的噪聲很是有效。濾波能夠消除或下降圖像中噪聲的影響,使用高斯濾波器主要是基於在濾波降噪的同時也能夠最大限度保留邊緣信息的考慮。this
高斯濾波實現步驟:spa
邊緣檢測是基於對圖像灰度差別運算實現的,因此若是輸入的是RGB彩色圖像,須要先進行灰度圖的轉換。
RGB轉換成灰度圖像的一個經常使用公式是:
\[ Gray = R*0.299 + G*0.587 + B*0.114 \]
注意通常狀況下圖像處理中彩色圖像各份量的排列順序是B、G、R。
RGB原圖像: 轉換後的灰度圖:
Java
代碼調用系統庫實現:
public static Mat RGB2Gray(Mat image) { // Gray = R*0.299 + G*0.587 + B*0.114 Mat gray = new Mat(); Imgproc.cvtColor(image, gray, Imgproc.COLOR_BGR2GRAY); return gray; }
一維高斯函數表述爲:
\[ G(x) = \frac {1}{\sqrt {2\pi}\sigma}\exp(-\frac {(x-\mu_x)^2}{2\sigma^2}) \]
對應圖形:
二維高斯函數表述爲:
\[ G(x) = \frac {1}{ {2\pi}\sigma^2}\exp(-\frac {(x-\mu_x)^2+(y-\mu_y)^2}{2\sigma^2}) \]
對應圖形:
一些重要特性說明:
在圖形或濾波效果上表現爲:\(σ\)越大,曲線越扁平,高斯濾波器的頻帶就越寬,平滑程度就越好,\(σ\)越小,曲線越瘦高,高斯濾波的頻帶就越窄,平滑程度也越弱;
如下對比一下不一樣大小標準差\(σ\)(Sigma)對圖像平滑的影響:
原圖:
卷積核尺寸5*5,σ=0.1:
卷積核尺寸5*5,σ=1:
對比能夠看到,Sigma(σ)越大,平滑效果越明顯。
濾波的主要目的是降噪,通常的圖像處理算法都須要先進行降噪。而高斯濾波主要使圖像變得平滑(模糊),同時也有可能增大了邊緣的寬度。
高斯函數是一個相似與正態分佈的中間大兩邊小的函數。
對於一個位置(m,n)的像素點,其灰度值(這裏只考慮二值圖)爲f(m,n)。
那麼通過高斯濾波後的灰度值將變爲:
\[ g_\sigma(m,n)=\frac 1 {2\pi\sigma^2} \exp(-\frac {m^2+n^2}{2\sigma^2})f(m,n) \]
其中,
\[ \frac 1 {{2\pi\sigma^2}} \exp(-\frac {m^2+n^2}{2\sigma^2}) \]
是二元高斯函數。
爲了儘量減小噪聲對邊緣檢測結果的影響,因此必須濾除噪聲以防止由噪聲引發的錯誤檢測。爲了平滑圖像,使用高斯濾波器與圖像進行卷積,該步驟將平滑圖像,以減小邊緣檢測器上明顯的噪聲影響。大小爲(2k+1)x(2k+1)的高斯濾波器核的生成方程式由下式給出:
下面是一個sigma = 1.4,尺寸爲3x3的高斯卷積核的例子(須要注意歸一化):
若圖像中一個3x3的窗口爲A,要濾波的像素點爲e,則通過高斯濾波以後,像素點e的亮度值爲:
其中*爲卷積符號,sum表示矩陣中全部元素相加求和。
重要的是須要理解,高斯卷積核大小的選擇將影響Canny檢測器的性能。尺寸越大,檢測器對噪聲的敏感度越低,可是邊緣檢測的定位偏差也將略有增長。通常5x5是一個比較不錯的trade off。
public static double[][] getGaussianArray(int size, double sigma) { double[][] array = new double[size][size]; double center_i, center_j; if (size % 2 == 1) { center_i = (double) (size / 2); center_j = (double) (size / 2); } else { center_i = (double) (size / 2) - 0.5; center_j = (double) (size / 2) - 0.5; } double sum = 0.0; for (int i = 0; i < size; i ++) for (int j = 0; j < size; j ++) { array[i][j] = Math.exp((-1.0) * ((i - center_i) * (i - center_i) + (j - center_j) * (j - center_j)) / 2.0 / sigma / sigma); sum += array[i][j]; } for (int i = 0; i < size; i ++) for (int j = 0; j < size; j++) array[i][j] /= sum; return array; }
Sigma爲1時,求得的3*3大小的高斯卷積核參數爲:
Sigma爲1,5*5大小的高斯卷積核參數爲:
在計算出高斯濾波卷積核以後就是用這個卷積核來掃過整張圖像,對每一個像素點進行加權平均。
加入了多線程的優化,代碼實現:
public static Mat greyGaussianFilter(Mat src, double[][] array, int size) throws InterruptedException { Mat temp = src.clone(); final CountDownLatch latch = new CountDownLatch(src.rows()); for (int i = 0; i < src.rows(); i ++) { int finalI = i; new Thread(() -> { for (int j = 0; j < src.cols(); j ++) if (finalI > (size / 2) - 1 && j > (size / 2) - 1 && finalI < src.rows() - (size / 2) && j < src.cols() - (size / 2)) { double sum = 0.0; for (int k = 0; k < size; k++) for (int l = 0; l < size; l++) sum += src.get(finalI - k + size / 2, j - l + size / 2)[0] * array[k][l]; temp.put(finalI, j, sum); } latch.countDown(); }).start(); } latch.await(); return temp; } public static Mat colorGaussianFilter(Mat src, int size, double sigma) throws InterruptedException { // return variable Mat ret = new Mat(); // list for merge and split List<Mat> channels = new ArrayList<>(); List<Mat> new_channels = new ArrayList<>(); Map<Integer, Mat> temp_channels = new TreeMap<>(); // split into 3 channels (r, g, b) Core.split(src, channels); // get gaussion array double[][] array = SmartGaussian.getGaussianArray(size, sigma); // multi-thread final CountDownLatch latch = new CountDownLatch(channels.size()); channels.forEach(mat -> { new Thread(() -> { Mat tmp = new Mat(); try { tmp = SmartGaussian.greyGaussianFilter(mat, array, size); } catch (InterruptedException e) { e.printStackTrace(); } temp_channels.put(channels.indexOf(mat), tmp); latch.countDown(); }).start(); }); latch.await(); for (int i = 0; i < channels.size(); i ++) new_channels.add(temp_channels.get(i)); Core.merge(new_channels, ret); return ret; }
效果圖:(左圖爲原圖,中間爲上面的實現,右邊是OpenCV
實現)
梯度的本意是一個向量(矢量),表示某一函數在該點處的方向導數沿着該方向取得最大值,即函數在該點處沿着該方向(此梯度的方向)變化最快,變化率最大(爲該梯度的模)。
設二元函數:
\[ z=f(x,y) \]
在平面區域D上具備一階連續偏導數,則對於每個點\(P(x,y)\)均可定出一個向量:
\[ \Big\{\frac {∂ f}{∂ x},\frac {∂ f} {∂ y} \Big\} =f_x(x,y)\vec i + f_y(x,y)\vec j \]
該函數就稱爲函數\(z=f(x,y)\)在點\(P(x,y)\)的梯度,記做\(\text{grad}~f(x,y)\)或\(\nabla f(x,y)\),即有:
\[ \text{grad}~f(x,y)=\nabla f(x,y)=\Big\{\frac {∂ f}{∂ x},\frac {∂ f} {∂ y} \Big\} =f_x(x,y)\vec i + f_y(x,y)\vec j \]
其中稱爲(二維的)向量微分算子或Nabla算子。
設 是方向\(l\)上的單位向量,則
因爲當方向\(l\)與梯度方向一致時,有。
因此當\(l\)與梯度方向一致時,方向導數有最大值,且最大值爲梯度的模,即
。
所以說,函數在一點沿梯度方向的變化率最大,最大值爲該梯度的模。
圖像灰度值的梯度可使用最簡單的一階有限差分來進行近似,使用如下圖像在x和y方向上偏導數的兩個矩陣:
計算公式爲:
其中\(f\)爲圖像灰度值,\(P[i,j]\)表明\([i,j]\)在\(X\)方向梯度幅值,\(Q[i,j]\)表明\([i,j]\)在\(Y\)方向的梯度幅值,\(M[i,j]\)是該點幅值,\(\theta[i,j]\)是梯度方向,也就是角度。
圖像中的邊緣能夠指向各個方向,所以Canny算法使用四個算子來檢測圖像中的水平、垂直和對角邊緣。邊緣檢測的算子(如Roberts,Prewitt,Sobel等)返回水平Gx和垂直Gy方向的一階導數值,由此即可以肯定像素點的梯度G和方向theta 。
其中G爲梯度強度, theta表示梯度方向,arctan爲反正切函數。下面以Sobel算子爲例講述如何計算梯度強度和方向。
x和y方向的Sobel算子分別爲:
其中Sx表示x方向的Sobel算子,用於檢測y方向的邊緣; Sy表示y方向的Sobel算子,用於檢測x方向的邊緣(邊緣方向和梯度方向垂直)。在直角座標系中,Sobel算子的方向以下圖所示。
若圖像中一個3x3的窗口爲A,要計算梯度的像素點爲e,則和Sobel算子進行卷積以後,像素點e在x和y方向的梯度值分別爲:
其中*爲卷積符號,sum表示矩陣中全部元素相加求和。根據公式(3-2)即可以計算出像素點e的梯度和方向。
下面是Sobel算子求梯度的java實現:
package edu.sfls.Jeff.JavaDev.CVLib; import org.opencv.core.Core; import org.opencv.core.CvType; import org.opencv.core.Mat; import java.util.concurrent.CountDownLatch; public class SmartSobel { private Mat gradientX = new Mat(), gradientY = new Mat(), gradientXY = new Mat(); private double[][] pointDirection; public SmartSobel() {} public void compute(Mat image) throws InterruptedException { pointDirection = new double[image.rows()][image.cols()]; for (int i = 0; i < image.rows(); i ++) for (int j = 0; j < image.cols(); j ++) pointDirection[i][j] = 0; gradientX = Mat.zeros(image.size(), CvType.CV_32SC1); gradientY = Mat.zeros(image.size(), CvType.CV_32SC1); gradientXY = Mat.zeros(image.size(), CvType.CV_32SC1); final CountDownLatch latch = new CountDownLatch(image.rows() - 2); for (int i = 1; i < image.rows() - 1; i ++) { int finalI = i; new Thread(() -> { for (int j = 1; j < image.cols() - 1; j++) { double gX = (-1) * image.get(finalI - 1, j - 1)[0] + 1 * image.get(finalI - 1, j + 1)[0] + (-2) * image.get(finalI, j - 1)[0] + 2 * image.get(finalI, j + 1)[0] + (-1) * image.get(finalI + 1, j - 1)[0] + 1 * image.get(finalI + 1, j + 1)[0]; double gY = 1 * image.get(finalI - 1, j - 1)[0] + 2 * image.get(finalI - 1, j)[0] + 1 * image.get(finalI - 1, j + 1)[0] + (-1) * image.get(finalI + 1, j - 1)[0] + (-2) * image.get(finalI + 1, j)[0] + (-1) * image.get(finalI + 1, j + 1)[0]; gradientY.put(finalI, j, Math.abs(gY)); gradientX.put(finalI, j, Math.abs(gX)); gradientXY.put(finalI, j, Math.sqrt(gX * gX + gY * gY)); // 防止除以0的狀況發生 if (gX == 0) gX = 0.00000000000000001; pointDirection[finalI][j] = Math.atan(gY / gX); } latch.countDown(); }).start(); } latch.await(); } public void convert() { Core.convertScaleAbs(gradientX, gradientX); Core.convertScaleAbs(gradientY, gradientY); Core.convertScaleAbs(gradientXY, gradientXY); } public Mat getGradientX() { return this.gradientX; } public Mat getGradientY() { return this.gradientY; } public Mat getGradientXY() { return this.gradientXY; } public double[][] getPointDirection() { return this.pointDirection; } }
因爲這裏使用的是\(Math.tan()\),因此最終的角度是映射在\([-\frac \pi 2, \frac \pi 2]\)的範圍以內的。若是使用\(Math.tan2()\)會映射到\([-\pi,\pi]\)的範圍內,而且無需考慮符號影響,更加精確。可是這裏咱們並不關心另外的一個\(\pi\)的狀況,咱們只關心其所在直線(這在後文中會提到,也就是非極大值抑制),因此無需多考慮。
X方向梯度圖: Y方向梯度圖:
X、Y方向梯度融合效果: Opencv Sobel函數效果:
非極大值抑制是一種邊緣稀疏技術,非極大值抑制的做用在於「瘦」邊。對圖像進行梯度計算後,僅僅基於梯度值提取的邊緣仍然很模糊。對於標準3,對邊緣有且應當只有一個準確的響應。而非極大值抑制則能夠幫助將局部最大值以外的全部梯度值抑制爲0,對梯度圖像中每一個像素進行非極大值抑制的算法是:
1) 將當前像素的梯度強度與沿正負梯度方向上的兩個像素進行比較。
2) 若是當前像素的梯度強度與另外兩個像素相比最大,則該像素點保留爲邊緣點,不然該像素點將被抑制。
一般爲了更加精確的計算,在跨越梯度方向的兩個相鄰像素之間使用線性插值來獲得要比較的像素梯度,現舉例以下:
圖3-2 梯度方向分割
如圖3-2所示,將梯度分爲8個方向,分別爲E、NE、N、NW、W、SW、S、SE,其中0表明\(0^\circ\sim45^\circ\),1表明\(45^\circ\sim90^\circ\),2表明\(-90^\circ\sim-45^\circ\),3表明\(-45^\circ\sim0^\circ\)。像素點P的梯度方向爲\(\theta\),則像素點P1和P2的梯度線性插值爲:
\[ \tan \theta = G_y ~/~G_x \\ G_{p1} = (1-\tan\theta)\times E + \tan\theta \times NE \\ G_{p2} = (1-\tan\theta)\times W + \tan\theta \times SW \\ \]
上面也只是圖中的狀況,具體狀況以下:
\[ \theta \in [0, \frac \pi 4]: \begin {cases} G_{p1} = (1-\tan\theta)\times E + \tan\theta \times NE \\ G_{p2} = (1-\tan\theta)\times W + \tan\theta \times SW \end {cases}\\\theta \in [\frac \pi 4, \frac \pi 2]: \begin {cases} G_{p1} = (1-\tan\theta)\times N + \tan\theta \times NE \\ G_{p2} = (1-\tan\theta)\times S + \tan\theta \times SW \end {cases}\\\theta \in [-\frac \pi 4, 0]: \begin {cases} G_{p1} = (1-\tan(-\theta))\times E + \tan(-\theta) \times SE \\ G_{p2} = (1-\tan(-\theta))\times W + \tan(-\theta) \times NW \end {cases}\\\theta \in [-\frac \pi 2, -\frac \pi 4]: \begin {cases} G_{p1} = (1-\tan(-\theta))\times S + \tan(-\theta) \times SE \\ G_{p2} = (1-\tan(-\theta))\times N + \tan(-\theta) \times NW \end {cases}\\ \]
所以非極大值抑制的僞代碼描寫以下:
須要注意的是,如何標誌方向並不重要,重要的是梯度方向的計算要和梯度算子的選取保持一致。
Java實現:
package edu.sfls.Jeff.JavaDev.CVLib; import org.opencv.core.CvType; import org.opencv.core.Mat; import java.util.concurrent.CountDownLatch; public class SmartNMS { public static Mat NMS(Mat gradientImage, double[][] pointDirection) throws InterruptedException { Mat outputImage = gradientImage.clone(); final CountDownLatch latch = new CountDownLatch(gradientImage.rows() - 2); for (int i = 1; i < gradientImage.rows() - 1; i ++) { int finalI = i; new Thread(() -> { for (int j = 1; j < gradientImage.cols() - 1; j ++) { double GP = gradientImage.get(finalI, j)[0], E = gradientImage.get(finalI, j + 1)[0], NE = gradientImage.get(finalI - 1, j + 1)[0], N = gradientImage.get(finalI - 1, j)[0], NW = gradientImage.get(finalI - 1, j - 1)[0], W = gradientImage.get(finalI, j - 1)[0], SW = gradientImage.get(finalI + 1, j - 1)[0], S = gradientImage.get(finalI + 1, j)[0], SE = gradientImage.get(finalI + 1, j + 1)[0]; double GP1 = 0, GP2 = 0; double theta = pointDirection[finalI][j]; if (theta >= 0 && theta <= Math.PI / 4) { GP1 = E * (1 - Math.tan(theta)) + NE * Math.tan(theta); GP2 = W * (1 - Math.tan(theta)) + SW * Math.tan(theta); } else if (theta > Math.PI / 4) { GP1 = N * (1 - 1 / Math.tan(theta)) + NE * 1 / Math.tan(theta); GP2 = S * (1 - 1 / Math.tan(theta)) + SW * 1 / Math.tan(theta); } else if (theta < 0 && theta >= -Math.PI / 4) { GP1 = E * (1 - Math.tan(-theta)) + SE * Math.tan(-theta); GP2 = W * (1 - Math.tan(-theta)) + NW * Math.tan(-theta); } else { GP1 = S * (1 - 1 / Math.tan(-theta)) + SE * 1 / Math.tan(-theta); GP2 = N * (1 - 1 / Math.tan(-theta)) + NW * 1 / Math.tan(-theta); } if (GP < GP1 || GP < GP2) outputImage.put(finalI, j, 0); } latch.countDown(); }).start(); } latch.await(); return outputImage; } }
在施加非極大值抑制以後,剩餘的像素能夠更準確地表示圖像中的實際邊緣。然而,仍然存在因爲噪聲和顏色變化引發的一些邊緣像素。爲了解決這些雜散響應,必須用弱梯度值過濾邊緣像素,並保留具備高梯度值的邊緣像素,能夠經過選擇高低閾值來實現。若是邊緣像素的梯度值高於高閾值,則將其標記爲強邊緣像素;若是邊緣像素的梯度值小於高閾值而且大於低閾值,則將其標記爲弱邊緣像素;若是邊緣像素的梯度值小於低閾值,則會被抑制。閾值的選擇取決於給定輸入圖像的內容。
雙閾值檢測的僞代碼描寫以下:
到目前爲止,被劃分爲強邊緣的像素點已經被肯定爲邊緣,由於它們是從圖像中的真實邊緣中提取出來的。然而,對於弱邊緣像素,將會有一些爭論,由於這些像素能夠從真實邊緣提取也能夠是因噪聲或顏色變化引發的。爲了得到準確的結果,應該抑制由後者引發的弱邊緣。一般,由真實邊緣引發的弱邊緣像素將鏈接到強邊緣像素,而噪聲響應未鏈接。爲了跟蹤邊緣鏈接,經過查看弱邊緣像素及其8個鄰域像素,只要其中一個爲強邊緣像素,則該弱邊緣點就能夠保留爲真實的邊緣。
抑制孤立邊緣點的僞代碼描述以下:
實現:(使用遞歸)
package edu.sfls.Jeff.JavaDev.CVLib; import org.opencv.core.Mat; import java.util.ArrayList; import java.util.List; public class SmartCanny { private List<Integer[]> highPoints = new ArrayList<Integer[]>(); private void DoubleThreshold(Mat image, double lowThreshold, double highThreshold) { for (int i = 1; i < image.rows() - 1; i ++) for (int j = 1; j < image.cols() - 1; j ++) if (image.get(i, j)[0] >= highThreshold) { image.put(i, j, 255); Integer[] p = new Integer[2]; p[0] = i; p[1] = j; highPoints.add(p); } else if (image.get(i, j)[0] < lowThreshold) image.put(i, j, 0); } private void DoubleThresholdLink(Mat image, double lowThreshold) { for (Integer[] p : highPoints) { DoubleThresholdLinkRecurrent(image, lowThreshold, p[0], p[1]); } for (int i = 1; i < image.rows() - 1; i ++) for (int j = 1; j < image.cols() - 1; j ++) if (image.get(i, j)[0] < 255) image.put(i, j, 0); } private void DoubleThresholdLinkRecurrent(Mat image, double lowThreshold, int i, int j) { if (i <= 0 || j <= 0 || i >= image.rows() - 1 || j >= image.cols() - 1) return; if (image.get(i - 1, j - 1)[0] >= lowThreshold && image.get(i - 1, j - 1)[0] < 255) { image.put(i - 1, j - 1, 255); DoubleThresholdLinkRecurrent(image, lowThreshold, i - 1, j - 1); } if (image.get(i - 1, j)[0] >= lowThreshold && image.get(i - 1, j)[0] < 255) { image.put(i - 1, j, 255); DoubleThresholdLinkRecurrent(image, lowThreshold, i - 1, j); } if (image.get(i - 1, j + 1)[0] >= lowThreshold && image.get(i - 1, j + 1)[0] < 255) { image.put(i - 1, j + 1, 255); DoubleThresholdLinkRecurrent(image, lowThreshold, i - 1, j + 1); } if (image.get(i, j - 1)[0] >= lowThreshold && image.get(i, j - 1)[0] < 255) { image.put(i, j - 1, 255); DoubleThresholdLinkRecurrent(image, lowThreshold, i, j - 1); } if (image.get(i, j + 1)[0] >= lowThreshold && image.get(i, j + 1)[0] < 255) { image.put(i, j + 1, 255); DoubleThresholdLinkRecurrent(image, lowThreshold, i, j + 1); } if (image.get(i + 1, j - 1)[0] >= lowThreshold && image.get(i + 1, j - 1)[0] < 255) { image.put(i + 1, j - 1, 255); DoubleThresholdLinkRecurrent(image, lowThreshold, i + 1, j - 1); } if (image.get(i + 1, j)[0] >= lowThreshold && image.get(i + 1, j)[0] < 255) { image.put(i + 1, j, 255); DoubleThresholdLinkRecurrent(image, lowThreshold, i + 1, j); } if (image.get(i + 1, j + 1)[0] >= lowThreshold && image.get(i + 1, j + 1)[0] < 255) { image.put(i + 1, j + 1, 255); DoubleThresholdLinkRecurrent(image, lowThreshold, i + 1, j + 1); } } public Mat Canny(Mat image, int size, double sigma, double lowThreshold, double highThreshold) throws InterruptedException { Mat tmp = SmartConverter.RGB2Gray((SmartGaussian.colorGaussianFilter(image, size, sigma))); SmartSobel ss = new SmartSobel(); ss.compute(tmp); ss.convert(); Mat ret = SmartNMS.NMS(ss.getGradientXY(), ss.getPointDirection()); this.DoubleThreshold(ret, lowThreshold, highThreshold); this.DoubleThresholdLink(ret, lowThreshold); return ret; } }
[1] 高斯濾波及高斯卷積核C++實現 https://blog.csdn.net/dcrmg/article/details/52304446
[2] 邊緣檢測之Canny http://www.javashuo.com/article/p-xoxvvnpc-ko.html
[3] Canny邊緣檢測及C++實現 https://blog.csdn.net/dcrmg/article/details/52344902
[4] Canny邊緣檢測算法 https://zhuanlan.zhihu.com/p/42122107
[5] Sobel算子及C++實現 https://blog.csdn.net/dcrmg/article/details/52280768