1、圖像銳化理論基礎
1、銳化的概念
圖像銳化的目的是使模糊的圖像變得清晰起來，主要用於增強圖像的灰度跳變部分，這一點與圖像平滑對灰度跳變的抑制正好相反。而且從算子可以看出來，平滑是基於對圖像領域的加權求和或者說積分運算的，而銳化則是通過其逆運算導數（梯度）或者說有限差分來實現的。
2、圖像的一階微分和二階微分的性質
圖像的銳化也就是增強圖像的突變部分，那么我們也就對圖像的恆定區域中，突變的開始點與結束點（台階和斜坡突變）及沿着灰度斜坡處的微分的性質。微分是對函數局部變化率的一種表示，那么對於一階微分有以下幾個性質：
- 在恆定的灰度區域，圖像的微分值為0.（灰度值沒有發生變換，自然微分為0）
- 在灰度台階或斜坡起點處微分值不為0.（台階是，灰度值的突變變化較大；斜坡則是灰度值變化較緩慢；灰度值發生了變化，微分值不為0）
- 沿着斜坡的微分值不為0.
二階微分，是一階微分的導數，和一階微分相對應，也有以下幾點性質:
- 在恆定區域二階微分值為0
- 在灰度台階或斜坡的起點處微分值不為0
- 沿着斜坡的微分值為0.
從以上圖像灰度的一階和二階微分的性質可以看出，在灰度值變化的地方，一階微分和二階微分的值都不為0；在灰度恆定的地方，微分值都為0.也就是說，不論是使用一階微分還是二階微分都可以得到圖像灰度的變化值。
3、圖像的一階微分和二階微分的計算
圖像可以看着是二維離散函數，對於圖像的一階微分其計算公式如下：
\begin{array}{l}{\frac{\partial f}{\partial x}=f(x+1)-f(x)} \\ {\frac{\partial f}{\partial y}=f(y+1)-f(y)}\end{array}
對於二階微分有:
\begin{array}{l}{\frac{\partial^{2} f}{\partial \tau^{2}}=f(x+1)+f(x-1)-2 f(x)} \\ {\frac{\partial^{2} f}{\partial y^{2}}=f(y+1)+f(y-1)-2 f(y)}\end{array}
在圖像處理中的一階微分通常使用梯度的幅值來實現，f在(x,y)處梯度的列向量為：
$\nabla f=\operatorname{grad}(f)=\left[ \begin{array}{l}{g_{x}} \\ {g_{y}}\end{array}\right]=\left[ \begin{array}{l}{\frac{\partial f}{\partial x}} \\ {\frac{\partial f}{\partial y}}\end{array}\right]$
該向量表示圖像中的像素在點(x,y)處灰度變化率最大的方向。
這個列向量的幅值就是圖像f(x,y)的梯度圖：
$M(x, y)=\operatorname{mag}(\nabla f)=\sqrt{g_{x}^{2}+g_{y}^{2}}$
由於求平方的根運算比較費時，通常可以使用絕對值的和來近似
$$
M(x, y) \approx\left|g_{x}\right|+\left|g_{y}\right|
$$
對於二階微分的計算，則有：
$$
\nabla^{2} f=\frac{\partial^{2} f}{\partial x^{2}}+\frac{\partial^{2} f}{\partial y^{2}}
$$
其中：
$$
\begin{aligned} \frac{\partial^{2} f}{\partial x^{2}} &=f(x+1, y)+f(x-1, y)-2 f(x, y) \\ \frac{\partial^{2} f}{\partial y^{2}} &=f(x, y+1)+f(x, y-1)-2 f(x, y) \end{aligned}
$$
2、基於一階導的梯度算子
1、Robort算子
根據梯度的定義
$\begin{aligned} g_{x} &=f(x+1, y)-f(x, y) \\ g_{y} &=f(x, y+1)-f(x, y) \end{aligned}$
可以得到模板$\left[ \begin{array}{ll}{-1} & {1}\end{array}\right]$和$\left[ \begin{array}{c}{-1} \\ {1}\end{array}\right]$
可以看出，得到的邊緣一部分是在內邊界，一部分是外邊界，並且，黃色像素點並未有計算結果，也就是，邊緣候選點丟失了一個。且使用該方法計算的圖像的梯度只是考慮單個像素的差值，並沒有利用到圖像的像素的鄰域特性。考慮到每個像素的某個鄰域的灰度變化。因此，通常不會簡單的利用梯度的定義進行梯度的計算，而是在像素的某個鄰域內設置梯度算子。
考慮，3X3區域的像素，使用如下矩陣表示：
$$
\left[ \begin{array}{lll}{z_{1}} & {z_{2}} & {z_{3}} \\ {z_{4}} & {z_{5}} & {z_{6}} \\ {z_{7}} & {z_{8}} & {z_{9}}\end{array}\right]
$$
令中心點表示圖像中任一像素，那么根據梯度的定義，其在x和y方向的梯度分別為：$g_{x}=z_{9}-z_{5}$和$g_{y}=z_{8}-z_{6}$，梯度圖像
$$
M(x, y) \approx\left|z_{9}-z_{5}\right|+\left|z_{8}-z_{6}\right|
$$
根據上述公式，Robert在1965年提出的Robert交叉算子
$$
\left[ \begin{array}{cc}{-1} & {0} \\ {0} & {1}\end{array}\right] \text { and } \left[ \begin{array}{cc}{0} & {-1} \\ {1} & {0}\end{array}\right]
$$
以下是一組Robort算子計算邊界的示例，可以幫助你對Robort算子進行理解：
示例:
//實現直方圖的反向投影
#include "stdafx.h"
#include <math.h>
#include <random>
#include <iostream>
#include <opencv2\core\core.hpp>
#include <opencv2\highgui\highgui.hpp>
#include <opencv2\imgproc\imgproc.hpp>
using namespace cv;
using namespace std;
void robert_grad(const Mat& src, Mat &dst)
{
Mat grad_x, grad_y;
Mat kernel_x = (Mat_<float>(2, 2) << -1, 0, 0, 1);
Mat kernel_y = (Mat_<float>(2, 2) << 0, -1, 1, 0);
filter2D(src, grad_x, CV_32F, kernel_x);
filter2D(src, grad_y, CV_32F, kernel_y);
//convertScaleAbs(grad_x, grad_x);
//convertScaleAbs(grad_y, grad_y);
//addWeighted(grad_x, 1, grad_y, 1, 0, dst);
magnitude(grad_x, grad_y, dst);
convertScaleAbs(dst, dst);
}
int main()
{
Mat srcImage = imread("111.jpg",0);
if (!srcImage.data)
{
cout << "圖像打開失敗！" << endl;
return -1;
}
Mat robort_Img;
robert_grad(srcImage, robort_Img);
imshow("原圖", srcImage);
imshow("Robort算子處理后", robort_Img);
waitKey();
return 0;
}
2、Sobel算子
Robert交叉算子的尺寸是偶數，偶數尺寸濾波器沒有對稱中心計算效率較低，所以通常濾波器的模板尺寸是奇數。
則有：
$$
\begin{aligned} g_{x} &=\left(z_{7}+2 z_{8}+z_{9}\right)-\left(z_{1}+2 z_{2}+z_{3}\right) \\ g_{y} &=\left(z_{3}+2 z_{0}+z_{9}\right)-\left(z_{1}+2 z_{4}+z_{7}\right) \end{aligned}
$$
利用上述公式可以得到如下兩個卷積模板，分別計算圖像在x和y方向的梯度
$$
\left[ \begin{array}{ccc}{-1} & {-2} & {-1} \\ {0} & {0} & {0} \\ {1} & {2} & {1}\end{array}\right]
$$
和
$$
\left[ \begin{array}{rrr}{-1} & {0} & {1} \\ {-2} & {0} & {2} \\ {-1} & {0} & {1}\end{array}\right]
$$
算法的相關實現如下所示：
bool Sobel(const Mat& image, Mat& result, int TYPE)
{
if (image.channels() != 1)
return false;
// 系數設置
int kx(0);
int ky(0);
if (TYPE == 0) {
kx = 0; ky = 1;
}
else if (TYPE == 1) {
kx = 1; ky = 0;
}
else if (TYPE == 2) {
kx = 1; ky = 1;
}
else
return false;
// 設置mask
float mask[3][3] = { {1,2,1},{0,0,0},{-1,-2,-1} };
Mat y_mask = Mat(3, 3, CV_32F, mask) / 8;
Mat x_mask = y_mask.t(); // 轉置
// 計算x方向和y方向上的濾波
Mat sobelX, sobelY;
filter2D(image, sobelX, CV_32F, x_mask);
filter2D(image, sobelY, CV_32F, y_mask);
sobelX = abs(sobelX);
sobelY = abs(sobelY);
// 梯度圖
Mat gradient = kx * sobelX.mul(sobelX) + ky * sobelY.mul(sobelY);
// 計算閾值
int scale = 4;
double cutoff = scale * mean(gradient)[0];
result.create(image.size(), image.type());
result.setTo(0);
for (int i = 1; i < image.rows - 1; i++)
{
float* sbxPtr = sobelX.ptr<float>(i);
float* sbyPtr = sobelY.ptr<float>(i);
float* prePtr = gradient.ptr<float>(i - 1);
float* curPtr = gradient.ptr<float>(i);
float* lstPtr = gradient.ptr<float>(i + 1);
uchar* rstPtr = result.ptr<uchar>(i);
// 閾值化和極大值抑制
for (int j = 1; j < image.cols - 1; j++)
{
if (curPtr[j] > cutoff && (
(sbxPtr[j] > kx*sbyPtr[j] && curPtr[j] > curPtr[j - 1] && curPtr[j] > curPtr[j + 1]) ||
(sbyPtr[j] > ky*sbxPtr[j] && curPtr[j] > prePtr[j] && curPtr[j] > lstPtr[j])))
rstPtr[j] = 255;
}
}
return true;
}
這里的非極大抑制Non-maximum suppression，這一步的目的是將模糊（blurred）的邊界變得清晰（sharp）。通俗的講，就是保留了每個像素點上梯度強度的極大值，而刪掉其他的值。
示例:
//實現直方圖的反向投影
#include "stdafx.h"
#include <math.h>
#include <random>
#include <iostream>
#include <opencv2\core\core.hpp>
#include <opencv2\highgui\highgui.hpp>
#include <opencv2\imgproc\imgproc.hpp>
using namespace cv;
using namespace std;
bool Sobel(const Mat& image, Mat& result, int TYPE)
{
if (image.channels() != 1)
return false;
// 系數設置
int kx(0);
int ky(0);
if (TYPE == 0) {
kx = 0; ky = 1;
}
else if (TYPE == 1) {
kx = 1; ky = 0;
}
else if (TYPE == 2) {
kx = 1; ky = 1;
}
else
return false;
// 設置mask
float mask[3][3] = { {1,2,1},{0,0,0},{-1,-2,-1} };
Mat y_mask = Mat(3, 3, CV_32F, mask) / 8;
Mat x_mask = y_mask.t(); // 轉置
// 計算x方向和y方向上的濾波
Mat sobelX, sobelY;
filter2D(image, sobelX, CV_32F, x_mask);
filter2D(image, sobelY, CV_32F, y_mask);
sobelX = abs(sobelX);
sobelY = abs(sobelY);
// 梯度圖
Mat gradient = kx * sobelX.mul(sobelX) + ky * sobelY.mul(sobelY);
// 計算閾值
int scale = 4;
double cutoff = scale * mean(gradient)[0];
result.create(image.size(), image.type());
result.setTo(0);
for (int i = 1; i < image.rows - 1; i++)
{
float* sbxPtr = sobelX.ptr<float>(i);
float* sbyPtr = sobelY.ptr<float>(i);
float* prePtr = gradient.ptr<float>(i - 1);
float* curPtr = gradient.ptr<float>(i);
float* lstPtr = gradient.ptr<float>(i + 1);
uchar* rstPtr = result.ptr<uchar>(i);
// 閾值化和極大值抑制
for (int j = 1; j < image.cols - 1; j++)
{
if (curPtr[j] > cutoff && (
(sbxPtr[j] > kx*sbyPtr[j] && curPtr[j] > curPtr[j - 1] && curPtr[j] > curPtr[j + 1]) ||
(sbyPtr[j] > ky*sbxPtr[j] && curPtr[j] > prePtr[j] && curPtr[j] > lstPtr[j])))
rstPtr[j] = 255;
}
}
return true;
}
int main()
{
Mat srcImage = imread("111.jpg",0);
if (!srcImage.data)
{
cout << "圖像打開失敗！" << endl;
return -1;
}
Mat sobel_x_Img, sobel_y_Img, sobel_xy_Img;
blur(srcImage, srcImage, Size(5, 5), Point(-1, -1));//由於直接濾波效果不好，這里用均值濾波濾除噪點
Sobel(srcImage, sobel_x_Img,0);
Sobel(srcImage, sobel_y_Img, 1);
Sobel(srcImage, sobel_xy_Img, 2);
imshow("原圖", srcImage);
imshow("sobel_x算子處理后", sobel_x_Img);
imshow("sobel_y算子處理后", sobel_y_Img);
imshow("sobel_xy算子處理后", sobel_xy_Img);
waitKey();
return 0;
}
程序運行效果如下：
同樣的，OpenCV提供了Sobel算子的API函數
Sobel(
InputArray Src // 輸入圖像
OutputArray dst// 輸出圖像，大小與輸入圖像一致
int depth // 輸出圖像深度.
Int dx. // X方向，幾階導數
int dy // Y方向，幾階導數.
int ksize, SOBEL算子kernel大小，必須是1、3、5、7、
double scale = 1
double delta = 0
int borderType = BORDER_DEFAULT
)
示例如下：
#include "stdafx.h"
#include<opencv2\opencv.hpp>
#include<opencv2\highgui\highgui.hpp>
using namespace std;
using namespace cv;
//邊緣檢測
int main()
{
Mat img = imread("111.jpg");
imshow("原始圖", img);
Mat grad_x, grad_y;
Mat abs_grad_x, abs_grad_y, dst;
//求x方向梯度
Sobel(img, grad_x, CV_16S, 1, 0, 3, 1, 1, BORDER_DEFAULT);
convertScaleAbs(grad_x, abs_grad_x);
imshow("x方向soble", abs_grad_x);
//求y方向梯度
Sobel(img, grad_y, CV_16S, 0, 1, 3, 1, 1, BORDER_DEFAULT);
convertScaleAbs(grad_y, abs_grad_y);
imshow("y向soble", abs_grad_y);
//合並梯度
addWeighted(abs_grad_x, 0.5, abs_grad_y, 0.5, 0, dst);
imshow("整體方向soble", dst);
waitKey(0);
}
程序運行結果如下：
3、基於二階微分的算子
1、Laplacian算子
二階微分算子的代表就是拉普拉斯算子，其定義如下：
$$
\nabla^{2} f=\frac{\partial^{2} f}{\partial x^{2}}+\frac{\partial^{2} f}{\partial y^{2}}
$$
其中：
$$
\begin{aligned} \frac{\partial^{2} f}{\partial x^{2}} &=f(x+1, y)+f(x-1, y)-2 f(x, y) \\ \frac{\partial^{2} f}{\partial y^{2}} &=f(x, y+1)+f(x, y-1)-2 f(x, y) \end{aligned}
$$
對於上述的3X3區域，則有
$$
\nabla^{2} f=z_{2}+z_{4}+z_{6}+z_{8}-4 z_{5}
$$
其得到的模板如下：
$$
\left[ \begin{array}{ccc}{0} & {1} & {0} \\ {1} & {-4} & {1} \\ {0} & {1} & {0}\end{array}\right]
$$
因為在銳化增強中，絕對值相同的正值和負值實際上表示相同的響應，故也等同於使用模板
$$
\left[ \begin{array}{ccc}{0} & {-1} & {0} \\ {-1} & {4} & {-1} \\ {0} & {-1} & {0}\end{array}\right]
$$
分析模板結構，可知模板對於90°的旋轉是各向同性的。所謂對於某角度各向同性是指把原圖像旋轉該角度后在進行濾波與對原圖像濾波在旋轉該角度的結果相同。這說明laplacian算子對於接近水平和接近豎直放i想的邊緣都有很好的增強，從而避免了在使用梯度算子時要進行多次濾波的麻煩。更進一步，我們還可以得到45°旋轉各向同性的濾波器：
$$
W_{3}=\left[ \begin{array}{ccc}{1} & {1} & {1} \\ {1} & {-8} & {1} \\ {1} & {1} & {1}\end{array}\right]
$$
和
$$
W_{4}=\left[ \begin{array}{ccc}{-1} & {-1} & {-1} \\ {-1} & {8} & {-1} \\ {-1} & {-1} & {-1}\end{array}\right]
$$
沿用高斯平滑的思想，根據到中心點的距離給模板周邊的點賦予不同的權重，還可得到如下模板：
$$
W_{5}=\left[ \begin{array}{ccc}{1} & {4} & {1} \\ {4} & {-20} & {4} \\ {1} & {4} & {1}\end{array}\right]
$$
在OpenCV中默認的計算方式如下，假設有一個5*5的小圖像，原始值依次為1,2,…25,如下圖紅色部分，首先將5*5映射到(5+3-1)*(5+3-1)大小，然后和3*3的kernel做累積和，因為計算結果有可能超出有效值[0, 255]范圍，因此在最終還需要判斷使其在有效范圍內，即小於0為0，大於255為255：
如坐標為(0,0)的計算過程為：12 = 7*0+6*1+7*0+2*1+1*(-4)+2*1+7*0+6*1+7*0.
OpenCV提供的Laplacian的函數API為：
Laplacian( src_gray, dst, ddepth, kernel_size, scale, delta, BORDER_DEFAULT );
參數意義為，
- src_gray，輸入圖像
- dst，Laplace操作結果
- ddepth，輸出圖像深度，因為輸入圖像一般為CV_8U，為了避免數據溢出，輸出圖像深度應該設置為CV_16S
- 濾波器孔徑尺寸；
- 比例因子；
- 表示結果存入目標圖
示例
#include "stdafx.h"
#include<opencv2\opencv.hpp>
#include<opencv2\highgui\highgui.hpp>
using namespace std;
using namespace cv;
//邊緣檢測
int main()
{
Mat img = imread("111.jpg");
imshow("原始圖", img);
Mat gray, dst, abs_dst;
//轉換為灰度圖
cvtColor(img, gray, COLOR_RGB2GRAY);
//使用Laplace函數
//第三個參數：目標圖像深度；第四個參數：濾波器孔徑尺寸；第五個參數：比例因子；第六個參數：表示結果存入目標圖
Laplacian(gray, dst, CV_16S, 3, 1, 0, BORDER_DEFAULT);
//計算絕對值，並將結果轉為8位
convertScaleAbs(dst, abs_dst);
waitKey(0);
}
2、Log算子
1980年，Marr和Hildreth提出將Laplace算子與高斯低通濾波相結合，提出了LOG（Laplace and Guassian）算子。
高斯函數和一級、二階導數如下圖所示：
步驟如下：
1.對圖像先進性高斯濾波（G × f），再進行Laplace算子運算Δ（G × f）；
2.保留一階導數峰值的位置，從中尋找Laplace過零點；
3.對過零點的精確位置進行插值估計。尋找Laplace零點的原因是如果圖像中有噪聲，噪聲在一階導數處也會取得極大值從而被當作邊緣。然而求解這個極大值也不方便，采用二階導數后，極大值點就為0了，因此值為0的地方就是邊界。
由上圖可以看出，高斯濾波之后邊緣信息才顯現出來。
LOG算子如下：
$$
\nabla^{2} G=\frac{\partial^{2} G}{\partial x^{2}}+\frac{\partial^{2} G}{\partial y^{2}}=\frac{-2 \sigma^{2}+x^{2}+y^{2}}{2 \pi \sigma^{6}} e^{-\left(x^{2}+y^{2}\right) / 2 \sigma^{2}}
$$
常用模板如下：
示例：
在OpenCV中實現Log算子的方法如下：
#include "stdafx.h"
#include <opencv2/highgui/highgui.hpp>
#include "opencv2/imgproc/imgproc.hpp"
#include <iostream>
int main()
{
cv::Mat srcImage =
cv::imread("111.jpg", 0);
if (!srcImage.data)
return -1;
// 高斯平滑
GaussianBlur(srcImage, srcImage, cv::Size(3, 3),
0, 0, cv::BORDER_DEFAULT);
cv::Mat dstImage;
// 拉普拉斯變換
Laplacian(srcImage, dstImage, CV_16S, 3);
convertScaleAbs(dstImage, dstImage);
cv::imshow("srcImage", srcImage);
cv::imshow("dstImage", dstImage);
cv::waitKey(0);
return 0;
}
4、導向濾波
導向圖濾波主要是通過一張引導圖G，對目標圖像P（輸入圖像）進行濾波處理，使得最后的輸出圖像大體上與目標圖像P相似，但是紋理部分與引導圖G相似。導向濾波不僅能夠實現雙邊濾波的邊緣平滑，而且在檢測到邊緣附近有很好的表現，可以應用在圖像增強、HDR壓縮、圖像摳圖以及圖像去霧等場景中。
導向濾波的實現原理如下：
對於普通的線性變換濾波器，輸入圖像是I，輸出圖像是S，導向函數為T，導向濾波定義在像素點j處的濾波結果可以表示為下式：
$$
S_{j}=\sum_{i} W_{i j}(T) I_{j}
$$
其中，i，j是圖像像素的下標，濾波器核函數Wij是圖像的加權系數，景點的雙邊濾波器核給定如下：
$$
W_{i j}^{b f}(I)=\frac{1}{K_{i}} \exp \left(-\frac{\left|x_{i}-x_{j}\right|^{2}}{\sigma_{s}^{2}}\right) \exp \left(-\frac{\left|I_{i}-I_{j}\right|^{2}}{\sigma_{r}^{2}}\right)
$$
其中x是像素坐標Ki是一個歸一化參數
$$
\sum_{i} W_{i j}^{b f}=1
$$
參數σs和σr代表空間相似度以及顏色范圍相似度的靈敏性
假設導向濾波器存在線性模型如下：
$$
s_{i}=a_{k} I_{i}+b_{k} \quad i \in W_{k}
$$
其中窗口函數Wk是S以I為中心在像素k位置形成的線性變換，敘述ak和bk滿足相同的線性系數，為確定線性系數，定義的輸出相應q應滿足下式：
$$
q_{i}=I_{i}-n_{i}
$$
其中ni表示噪聲，局部線性模型確定輸入與輸出相應的邊緣，就可以最小化q與I之間的差異，同事保持原線性變換。
最小化窗口函數Wk：
$$
E\left(a_{k}, b_{k}\right)=\sum_{i \in w k}\left(\left(a_{k} I_{i}+b_{k}-p_{i}\right)^{2}+\psi a_{k}^{2}\right)
$$
其中Ψ是一個正則化參數，根據線性變換得到系數ak和bk：
$$
\begin{array}{l}{a_{k}=\frac{\frac{1}{|w|} \sum_{i \in w k} T_{i} I_{i}-\mu_{i} \overline{I_{k}}}{\sigma_{k}^{2}+\psi}} \\ {b_{k}=\overline{I_{k}}-a_{k} \mu_{k}}\end{array}
$$
其中μk與σk是導向圖像I的窗口wk的均值和方差，|w|是像素的總個數，
$$
\overline{I_{k}}=\frac{1}{|w|} \sum_{i \in w k} I_{i}
$$
導向濾波實現步驟如下：
（1.利用boxFilter濾波器完成均值計算，其中均值包括導向均值、原始均值、互相關均值以及自相關均值
（2.根據均值計算相關系數參數，包括自相關與互相關方差。
（3.計算窗口線性變換參數系數a、b。
（4.根據公式計算參數a、b的均值。
（5.利用參數得到導向濾波輸出矩陣S。
用程序表示如下所示：
Mat guidedfilter(Mat &srcImage, Mat &srcClone, int r, double eps)
{
//轉換源圖像信息
srcImage.convertTo(srcImage, CV_64FC1);
srcClone.convertTo(srcClone, CV_64FC1);
int NumRows = srcImage.rows;
int NumCols = srcImage.cols;
Mat boxResult;
//下面按照步驟進行導向濾波操作
/////////////////////////////////////////////////////////////
//步驟一：計算均值
boxFilter(Mat::ones(NumRows, NumCols, srcImage.type()),
boxResult, CV_64FC1, Size(r, r));
//生成導向均值mean_I
Mat mean_I;
boxFilter(srcImage, mean_I, CV_64FC1, Size(r, r));
//生成原始均值mean_P
Mat mean_P;
boxFilter(srcClone, mean_P, CV_64FC1, Size(r, r));
//生成互相關均值mean_IP
Mat mean_IP;
boxFilter(srcImage.mul(srcClone), mean_IP,
CV_64FC1, Size(r, r));
Mat cov_IP = mean_IP - mean_I.mul(mean_P);
//生成自相關均值mean_II
Mat mean_II;
//應用盒濾波計算相關均值
boxFilter(srcImage.mul(srcImage), mean_II, CV_64FC1, Size(r, r));
//步驟二：計算相關系數
Mat var_I = mean_II - mean_I.mul(mean_I);
Mat var_IP = mean_IP - mean_I.mul(mean_P);
//步驟三:計算參數系數a,b
Mat a = cov_IP / (var_I + eps);
Mat b = mean_P = a.mul(mean_I);
//步驟四：計算系數a,b的均值
Mat mean_a;
boxFilter(a, mean_a, CV_64FC1, Size(r, r));
mean_a = mean_a / boxResult;
Mat mean_b;
boxFilter(b, mean_b, CV_64FC1, Size(r, r));
mean_b = mean_b / boxResult;
//步驟五：生成輸出矩陣
Mat resultMat = mean_a.mul(srcImage) + mean_b;
return resultMat;
}
示例:
#include "stdafx.h"
#include <iostream>
#include <opencv2\core\core.hpp>
#include <opencv2\highgui\highgui.hpp>
#include <opencv2\imgproc\imgproc.hpp>
using namespace cv;
using namespace std;
//導向濾波器
Mat guidedfilter(Mat &srcImage, Mat &srcClone, int r, double eps);
int main()
{
Mat srcImage = imread("111.jpg");
if (srcImage.empty())
{
cout << "讀入圖片錯誤！" << endl;
system("pause");
return -1;
}
//進行通道分離
vector<Mat>vSrcImage, vResultImage;
split(srcImage, vSrcImage);
Mat resultMat;
for (int i = 0; i < 3; i++)
{
//分通道轉換成浮點型數據
Mat tempImage;
vSrcImage[i].convertTo(tempImage, CV_64FC1, 1.0 / 255.0);
Mat p = tempImage.clone();
//分別進行導向濾波
Mat resultImage = guidedfilter(tempImage, p, 4, 0.01);
vResultImage.push_back(resultImage);
}
//通道結果合並
merge(vResultImage, resultMat);
imshow("原圖像", srcImage);
imshow("導向濾波后圖像", resultMat);
waitKey(0);
return 0;
}
Mat guidedfilter(Mat &srcImage, Mat &srcClone, int r, double eps)
{
//轉換源圖像信息
srcImage.convertTo(srcImage, CV_64FC1);
srcClone.convertTo(srcClone, CV_64FC1);
int NumRows = srcImage.rows;
int NumCols = srcImage.cols;
Mat boxResult;
//下面按照步驟進行導向濾波操作
/////////////////////////////////////////////////////////////
//步驟一：計算均值
boxFilter(Mat::ones(NumRows, NumCols, srcImage.type()),
boxResult, CV_64FC1, Size(r, r));
//生成導向均值mean_I
Mat mean_I;
boxFilter(srcImage, mean_I, CV_64FC1, Size(r, r));
//生成原始均值mean_P
Mat mean_P;
boxFilter(srcClone, mean_P, CV_64FC1, Size(r, r));
//生成互相關均值mean_IP
Mat mean_IP;
boxFilter(srcImage.mul(srcClone), mean_IP,
CV_64FC1, Size(r, r));
Mat cov_IP = mean_IP - mean_I.mul(mean_P);
//生成自相關均值mean_II
Mat mean_II;
//應用盒濾波計算相關均值
boxFilter(srcImage.mul(srcImage), mean_II, CV_64FC1, Size(r, r));
//步驟二：計算相關系數
Mat var_I = mean_II - mean_I.mul(mean_I);
Mat var_IP = mean_IP - mean_I.mul(mean_P);
//步驟三:計算參數系數a,b
Mat a = cov_IP / (var_I + eps);
Mat b = mean_P = a.mul(mean_I);
//步驟四：計算系數a,b的均值
Mat mean_a;
boxFilter(a, mean_a, CV_64FC1, Size(r, r));
mean_a = mean_a / boxResult;
Mat mean_b;
boxFilter(b, mean_b, CV_64FC1, Size(r, r));
mean_b = mean_b / boxResult;
//步驟五：生成輸出矩陣
Mat resultMat = mean_a.mul(srcImage) + mean_b;
return resultMat;
}
5、canny算子——邊緣檢測方法的提出
Canny算子是John F. Canny於 1986 年開發出來的一個多級邊緣檢測算法。Canny 的目標是找到一個最優的邊緣檢測算法，最優邊緣檢測的含義是：
-
應用高斯濾波來平滑圖像，目的是去除噪聲，任何邊緣檢測算法都不可能在未經處理的原始數據上很好地工作，所以第一步是對原始數據與高斯 mask 作 卷積，得到的圖像與原始圖像相比有些輕微的模糊（blurred）。
-
找尋圖像的強度梯度（intensity gradients），即找尋一幅圖像中灰度強度變化最強的位置。
-
應用非最大抑制（non-maximum suppression）技術來消除邊誤檢（本來不是但檢測出來是）。就是保留了每個像素點上梯度強度的極大值，而刪掉其他的值。
-
應用雙閾值的方法來決定可能的（潛在的）邊界
-
利用滯后技術來跟蹤邊界
非最大抑制
非最大抑制的目的是將模糊（blurred）的邊界變得清晰（sharp）。通俗的講，就是保留了每個像素點上梯度強度的極大值，而刪掉其他的值。對於每個像素點，進行如下操作：
圖中的數字代表了像素點的梯度強度，箭頭方向代表了梯度方向。以第二排第三個像素點為例，由於梯度方向向上，則將這一點的強度（7）與其上下兩個像素點的強度（5和4）比較，由於這一點強度最大，則保留。
雙閾值
雙閾值的技術即設定一個閾值上界和閾值下界（opencv中通常由人為指定的），圖像中的像素點如果大於閾值上界則認為必然是邊界（稱為強邊界，strong edge），小於閾值下界則認為必然不是邊界，兩者之間的則認為是候選項（稱為弱邊界，weak edge），需進行進一步處理。
滯后技術
滯后技術即將和強邊界相連的弱邊界認為是邊界，其他的弱邊界則被抑制。
OpenCV提供了對應的Canny算子的API函數，函數如下所示：
void Canny(inputArray,outputArray,double threshold1,double threshold2,int apertureSize=3,bool L2gradient=false) *第一個參數，輸入圖像，且需為單通道8位圖像。 *第二個參數，輸出的邊緣圖。 *第三個參數，第一個滯后性閾值。用於邊緣連接。 *第四個參數，第二個滯后性閾值。用於控制強邊緣的初始段，高低閾值比在2:1到3:1之間。 *第五個參數，表明應用sobel算子的孔徑大小，默認值為3。 *第六個參數，bool類型L2gradient，一個計算圖像梯度幅值的標識，默認值false。
示例：
#include "stdafx.h"
#include <opencv2/core/core.hpp>
#include <opencv2/highgui/highgui.hpp>
#include<opencv2/imgproc/imgproc.hpp>
#include<iostream>
#include"math.h"
using namespace std;
using namespace cv;
Mat Img_in, Img_gray, Img_out,Img_canny;
Mat scr;
int main()
{
Img_in = imread("111.jpg");
int rows = Img_in.rows;
int cols = Img_in.cols;//獲取圖像尺寸
cvtColor(Img_in, Img_gray, CV_BGR2GRAY);
imshow("【灰度圖】", Img_gray);//轉化為灰度圖
//step1:高斯平滑
Mat img_filt;
GaussianBlur(Img_gray, Img_out, Size(3, 3), 0, 0);
Canny(Img_out, Img_canny, 50, 150, 3);
imshow("自帶函數結果", Img_canny);
//adaptiveThreshold(img_filt , Img_out , 255 ,ADAPTIVE_THRESH_MEAN_C , THRESH_BINARY,min(rows,cols), 0);
//imshow("【二值圖】",Img_out );
Img_out.convertTo(Img_out, CV_32FC1); //將圖像轉換為float或double型，否則算梯度會報錯
/*step2:計算梯度（幅度和方向）
選擇一階差分卷積模板：
dx=[-1,-1;1,1] dy=[1,-1;1,-1]
*/
Mat gy = (Mat_<char>(2, 2) << 1, -1,
1, -1);
//定義一階差分卷積梯度模板
Mat gx = (Mat_<char>(2, 2) << -1, -1,
1, 1); //定義一階差分卷積梯度模板
Mat img_gx, img_gy, img_g;//定義矩陣
Mat img_dir = Mat::zeros(rows, cols, CV_32FC1);//定義梯度方向矩陣,計算角度為float型
filter2D(Img_out, img_gx, Img_out.depth(), gx); //獲取x方向的梯度圖像.使用梯度模板進行二維卷積,結果與原圖像大小相同
filter2D(Img_out, img_gy, Img_out.depth(), gy); //獲取x方向的梯度圖像.使用梯度模板進行二維卷積,結果與原圖像大小相同
img_gx = img_gx.mul(img_gx);//點乘（每個像素值平方）
img_gy = img_gy.mul(img_gy);//點乘（每個像素值平方）
img_g = img_gx + img_gy;
sqrt(img_g, img_g); //梯度幅值圖像
imshow("梯度圖", img_g);
//求梯度方向圖像
for (int i = 0; i < rows; i++)
{
for (int j = 0; j < cols; j++)
{
img_dir.at<float>(i, j) = fastAtan2(img_gy.at<float>(i, j), img_gx.at<float>(i, j));//求角度
}
}
/* step3:對梯度幅值進行非極大值抑制
首先將角度划分成四個方向范圍:水平(0°)、45°、垂直(90°)、135°
*/
Mat Nms = Mat::zeros(rows, cols, CV_32FC1);//定義一個非極大值抑制圖像，float型
for (int i = 0; i < rows; i++)
{
for (int j = 0; j < cols; j++)
{
if (img_dir.at<float>(i, j) <= 22.5 && img_dir.at<float>(i, j) >= 0 || img_dir.at<float>(i, j) >= 157.5 && img_dir.at<float>(i, j) <= 202.5
|| img_dir.at<float>(i, j) >= 337.5 && img_dir.at<float>(i, j) <= 360)
img_dir.at<float>(i, j) = 0;
else if (img_dir.at<float>(i, j) > 22.5 && img_dir.at<float>(i, j) <= 67.5 || img_dir.at<float>(i, j) > 202.5 && img_dir.at<float>(i, j) <= 247.5)
img_dir.at<float>(i, j) = 45;
else if (img_dir.at<float>(i, j) > 67.5 && img_dir.at<float>(i, j) <= 112.5 || img_dir.at<float>(i, j) > 247.5 && img_dir.at<float>(i, j) <= 292.5)
img_dir.at<float>(i, j) = 90;
else if (img_dir.at<float>(i, j) > 112.5 && img_dir.at<float>(i, j) < 157.5 || img_dir.at<float>(i, j) > 292.5 && img_dir.at<float>(i, j) < 337.5)
img_dir.at<float>(i, j) = 135;
}
}
for (int i = 1; i < rows - 1; i++)
{
for (int j = 1; j < cols - 1; j++)
{
if (img_dir.at<float>(i, j) == 90 && img_g.at<float>(i, j) == max(img_g.at<float>(i, j), max(img_g.at<float>(i, j + 1), img_g.at<float>(i, j - 1))))
Nms.at<float>(i, j) = img_g.at<float>(i, j);
else if (img_dir.at<float>(i, j) == 45 && img_g.at<float>(i, j) == max(img_g.at<float>(i, j), max(img_g.at<float>(i - 1, j + 1), img_g.at<float>(i + 1, j - 1))))
Nms.at<float>(i, j) = img_g.at<float>(i, j);
else if (img_dir.at<float>(i, j) == 0 && img_g.at<float>(i, j) == max(img_g.at<float>(i, j), max(img_g.at<float>(i - 1, j), img_g.at<float>(i + 1, j))))
Nms.at<float>(i, j) = img_g.at<float>(i, j);
else if (img_dir.at<float>(i, j) == 135 && img_g.at<float>(i, j) == max(img_g.at<float>(i, j), max(img_g.at<float>(i - 1, j - 1), img_g.at<float>(i + 1, j + 1))))
Nms.at<float>(i, j) = img_g.at<float>(i, j);
}
}
/*step4:雙閾值檢測和連接邊緣
*/
Mat img_dst = Mat::zeros(rows, cols, CV_32FC1);//定義一個雙閾值圖像，float型
double TH, TL;
double maxVal = 0;//必須為double類型，且必須賦初值，否則報錯
Nms.convertTo(Nms, CV_64FC1); //為了計算，將非極大值抑制圖像轉為double型
minMaxLoc(Nms, NULL, &maxVal, NULL, NULL); //求矩陣 Nms最大值
TH = 0.5*maxVal;//高閾值
TL = 0.3*maxVal;//低閾值
for (int i = 0; i < rows; i++)
{
for (int j = 0; j < cols; j++)
{
if (Nms.at<double>(i, j) < TL)
img_dst.at<float>(i, j) = 0;
else if (Nms.at<double>(i, j) > TH)
img_dst.at<float>(i, j) = 1;
else if (Nms.at<double>(i - 1, j - 1) < TL || Nms.at<double>(i - 1, j) < TL || Nms.at<double>(i - 1, j + 1) < TL ||
Nms.at<double>(i, j - 1) < TL || Nms.at<double>(i, j + 1) < TL || Nms.at<double>(i + 1, j - 1) < TL ||
Nms.at<double>(i + 1, j) < TL || Nms.at<double>(i + 1, j + 1) < TL)
img_dst.at<float>(i, j) = 1;
}
}
imshow("非極大值抑制圖", Nms);
imshow(" 邊緣檢測圖", img_dst);
imwrite(" 邊緣檢測效果圖.jpg", img_dst);//保存圖像
waitKey(0);
return 0;
}
