影像處理中的微分算子原理與實作
一階微分邊緣算子：經典算子比如：Roberts(羅伯特)、Prewitt（普魯伊特）、Sobel(索貝爾)，Canny（坎尼）等，
二階微分邊緣算子：Laplacian算子，LoG（ Laplace of Gaussian function）邊緣檢測算子和DoG（Difference of Gaussian）高斯差分算子，
-
Roberts算子
1963年，Roberts算子，又稱羅伯茨算子，是一種最簡單的算子，是一種利用區域差分算子尋找邊緣的算子，他采用對角線方向相鄰兩象素之差近似梯度幅值檢測邊緣，檢測垂直邊緣的效果好于斜向邊緣，定位精度高，對噪聲敏感,無法抑制噪聲的影響，
Roberts算子的模板分為水平方向和垂直方向，如下式所示，從其模板可以看出，Roberts算子能較好的增強正負45度的影像邊緣，
void quick_opencv::roberts_Demo(Mat &img)
{
// Roberts算子邊緣檢測
cvtColor(img,img,COLOR_BGR2GRAY);
Mat grad;
grad.create(img.size(), CV_8UC1);
for (int i = 1; i < img.rows - 1; i++)
{
for (int j = 1; j < img.cols - 1; j++)
{
grad.at<uchar>(i, j) = saturate_cast<uchar>(fabs(img.at<uchar>(i, j) - img.at<uchar>(i - 1, j - 1)) + fabs(img.at<uchar>(i, j - 1) - img.at<uchar>(i - 1, j)));
}
}
imshow("Roberts算子", grad);
}
-
Prewitt算子
1970年，Prewitt算子來自J.M.S. Prewitt “Object Enhancement and Extraction” in “Picture processing and Psychopictorics”, Academic Press
Prewitt算子是一種影像邊緣檢測的微分算子，其原理是利用特定區域內像素灰度值產生的差分實作邊緣檢測，由于Prewitt算子采用 33 模板對區域內的像素值進行計算，而Robert算子的模板為 22，故Prewitt算子的邊緣檢測結果在水平方向和垂直方向均比Robert算子更加明顯，Prewitt算子適合用來識別噪聲較多、灰度漸變的影像，其計算公式如下所示：
void getPrewitt_oper(Mat& getPrewitt_horizontal, Mat& getPrewitt_vertical, Mat& getPrewitt_Diagonal1, Mat& getPrewitt_Diagonal2) {
//水平方向
getPrewitt_horizontal = (Mat_<float>(3, 3) << -1, -1, -1, 0, 0, 0, 1, 1, 1);
//垂直方向
getPrewitt_vertical = (Mat_<float>(3, 3) << -1, 0, 1, -1, 0, 1, -1, 0, 1);
//對角135°
getPrewitt_Diagonal1 = (Mat_<float>(3, 3) << 0, 1, 1, -1, 0, 1, -1, -1, 0);
//對角45°
getPrewitt_Diagonal2 = (Mat_<float>(3, 3) << -1, -1, 0, -1, 0, 1, 0, 1, 1);
//逆時針反轉180°得到卷積核
flip(getPrewitt_horizontal, getPrewitt_horizontal, -1);
flip(getPrewitt_vertical, getPrewitt_vertical, -1);
flip(getPrewitt_Diagonal1, getPrewitt_Diagonal1, -1);
flip(getPrewitt_Diagonal2, getPrewitt_Diagonal2, -1);
}
void edge_Prewitt(Mat& src, Mat& dst1, Mat& dst2, Mat& dst3, Mat& dst4, Mat& dst, int ddepth, double delta = 0, int borderType = BORDER_DEFAULT) {
//獲取Prewitt算子
Mat getPrewitt_horizontal;
Mat getPrewitt_vertical;
Mat getPrewitt_Diagonal1;
Mat getPrewitt_Diagonal2;
getPrewitt_oper(getPrewitt_horizontal, getPrewitt_vertical, getPrewitt_Diagonal1, getPrewitt_Diagonal2);
//卷積得到水平方向邊緣
filter2D(src, dst1, ddepth, getPrewitt_horizontal, Point(-1, -1), delta, borderType);
//卷積得到4垂直方向邊緣
filter2D(src, dst2, ddepth, getPrewitt_vertical, Point(-1, -1), delta, borderType);
//卷積得到45°方向邊緣
filter2D(src, dst3, ddepth, getPrewitt_Diagonal1, Point(-1, -1), delta, borderType);
//卷積得到135°方向邊緣
filter2D(src, dst4, ddepth, getPrewitt_Diagonal2, Point(-1, -1), delta, borderType);
//邊緣強度（近似）
convertScaleAbs(dst1, dst1); //求絕對值并轉為無符號8位圖
convertScaleAbs(dst2, dst2);
convertScaleAbs(dst3, dst3); //求絕對值并轉為無符號8位圖
convertScaleAbs(dst4, dst4);
dst = dst1 + dst2;
}
void quick_opencv::prewitt_Demo(Mat &img)
{
Mat dst, dst1, dst2, dst3, dst4;
cvtColor(img, img, COLOR_BGR2GRAY);
//注意：要采用CV_32F，因為有些地方卷積后為負數，若用8位無符號，則會導致這些地方為0
edge_Prewitt(img, dst1, dst2, dst3, dst4, dst, CV_32F);
namedWindow("水平邊緣", WINDOW_NORMAL);
imshow("水平邊緣", dst1);
namedWindow("垂直邊緣", WINDOW_NORMAL);
imshow("垂直邊緣", dst2);
namedWindow("45°邊緣", WINDOW_NORMAL);
imshow("45°邊緣", dst3);
namedWindow("135°邊緣", WINDOW_NORMAL);
imshow("135°邊緣", dst4);
namedWindow("邊緣強度", WINDOW_NORMAL);
imshow("邊緣強度", dst);
}
-
Sobel算子
1973年，Sobel邊緣算子，當年作者并沒有公開發表過論文，僅僅是在一次博士生課題討論會(1968)上提出(“A 3x3 Isotropic Gradient Operator for Image Processing”)，后在1973年出版的一本專著（“Pattern Classification and Scene Analysis”）的腳注里作為注釋出現和公開的，
Sobel算子是一種用于邊緣檢測的離散微分算子，它結合了高斯平滑和微分求導，該算子用于計算影像明暗程度近似值，根據影像邊緣旁邊明暗程度把該區域內超過某個數的特定點記為邊緣，Sobel算子在Prewitt算子的基礎上增加了權重的概念，認為相鄰點的距離遠近對當前像素點的影響是不同的，距離越近的像素點對應當前像素的影響越大，從而實作影像銳化并突出邊緣輪廓，
Sobel算子根據像素點上下、左右鄰點灰度加權差，在邊緣處達到極值這一現象檢測邊緣，對噪聲具有平滑作用，提供較為精確的邊緣方向資訊，因為Sobel算子結合了高斯平滑和微分求導（分化），因此結果會具有更多的抗噪性，當對精度要求不是很高時，Sobel算子是一種較為常用的邊緣檢測方法，
Sobel算子的邊緣定位更準確，常用于噪聲較多、灰度漸變的影像，其演算法模板如下面的公式所示，其中 表示水平方向， 表示垂直方向，
void quick_opencv::sobel_Demo(Mat &img)
{
cvtColor(img, img, COLOR_BGR2GRAY);
Mat imageX = Mat::zeros(img.size(), CV_16SC1);
Mat imageY = Mat::zeros(img.size(), CV_16SC1);
Mat imageXY = Mat::zeros(img.size(), CV_16SC1);
Mat imageX8UC;
Mat imageY8UC;
Mat imageXY8UC;
GaussianBlur(img, img, Size(3, 3), 0); //高斯濾波器（模糊/平滑/近似）消除噪點
uchar *P = img.data;
uchar *PX = imageX.data;
uchar *PY = imageY.data;
int step = img.step;
int stepXY = imageX.step;
for (int i = 1;i < img.rows - 1;i++)
{
for (int j = 1;j < img.cols - 1;j++)
{
// 通過指標遍歷影像上每一個像素
// 求出X,Y方向的導數（梯度） sobel算子加權的結果
PX[i*imageX.step + j * (stepXY / step)] = abs(P[(i - 1)*step + j + 1] + P[i*step + j + 1] * 2 + P[(i + 1)*step + j + 1] - P[(i - 1)*step + j - 1] - P[i*step + j - 1] * 2 - P[(i + 1)*step + j - 1]);
PY[i*imageX.step + j * (stepXY / step)] = abs(P[(i + 1)*step + j - 1] + P[(i + 1)*step + j] * 2 + P[(i + 1)*step + j + 1] - P[(i - 1)*step + j - 1] - P[(i - 1)*step + j] * 2 - P[(i - 1)*step + j + 1]);
}
}
addWeighted(imageX, 0.5, imageY, 0.5, 0, imageXY);//融合X、Y方向的梯度
convertScaleAbs(imageX, imageX8UC);
convertScaleAbs(imageY, imageY8UC);
convertScaleAbs(imageXY, imageXY8UC); //轉換為8bit影像
Mat imageSobel;
Mat x_grad, y_grad;
Sobel(img, x_grad, CV_16S, 1, 0, 3);
Sobel(img, y_grad, CV_16S, 0, 1, 3);
convertScaleAbs(x_grad, x_grad);
convertScaleAbs(y_grad, y_grad);
addWeighted(x_grad, 0.5, y_grad, 0.5, 0, imageSobel);
imshow("Source Image", img);
imshow("X Direction", imageX8UC);
imshow("Y Direction", imageY8UC);
imshow("XY Direction", imageXY8UC);
imshow("Opencv Soble", imageSobel);
}
-
Canny算子
1986年，Canny邊緣檢測算子是John F. Canny開發出來的一個多級邊緣檢測演算法，更為重要的是Canny創立了“邊緣檢測計算理論”（computational theory of edge detection）解釋這項技術如何作業，到今天已經30多年過去了，但Canny演算法仍然是影像邊緣檢測演算法中最經典、先進的演算法之一，
1、高斯平滑
2、計算梯度幅度和方向
可選用的模板：soble算子、Prewitt算子、Roberts模板等等;
一般采用soble算子，OpenCV也是如此，利用soble水平和垂直算子與輸入影像卷積計算dx、dy
3、根據角度對幅值進行非極大值抑制
沿著梯度方向對幅值進行非極大值抑制，而非邊緣方向
在每一點上，領域中心 x 與沿著其對應的梯度方向的兩個像素相比，若中心像素為最大值，則保留，否則中心置0，這樣可以抑制非極大值，保留區域梯度最大的點，以得到細化的邊緣，
4、用雙閾值演算法檢測和連接邊緣
選取系數TH和TL，比率為2:1或3:1，（一般取TH=0.3或0.2,TL=0.1）；
將小于低閾值的點拋棄，賦0；將大于高閾值的點立即標記（這些點為確定邊緣點），賦1或255；
將小于高閾值，大于低閾值的點使用8連通區域確定（即：只有與TH像素連接時才會被接受，成為邊緣點，賦 1或255）
void quick_opencv::canny_Demo(Mat &img)
{
ConvertRGB2GRAY(img, imageGray); //RGB轉換為灰度圖
imshow("Gray Image", imageGray);
int size = 5; //定義卷積核大小
double **gaus = new double *[size]; //卷積核陣列
for (int i = 0;i < size;i++)
{
gaus[i] = new double[size]; //動態生成矩陣
}
GetGaussianKernel(gaus, 5, 1); //生成5*5 大小高斯卷積核，Sigma=1；
imageGaussian = Mat::zeros(imageGray.size(), CV_8UC1);
GaussianFilter(imageGray, imageGaussian, gaus, 5); //高斯濾波
imshow("Gaussian Image", imageGaussian);
Mat imageSobelY;
Mat imageSobelX;
double *pointDirection = new double[(imageSobelX.cols - 1)*(imageSobelX.rows - 1)]; //定義梯度方向角陣列
SobelGradDirction(imageGaussian, imageSobelX, imageSobelY, pointDirection); //計算X、Y方向梯度和方向角
imshow("Sobel Y", imageSobelY);
imshow("Sobel X", imageSobelX);
Mat SobelGradAmpl;
SobelAmplitude(imageSobelX, imageSobelY, SobelGradAmpl); //計算X、Y方向梯度融合幅值
imshow("Soble XYRange", SobelGradAmpl);
Mat imageLocalMax;
LocalMaxValue(SobelGradAmpl, imageLocalMax, pointDirection); //區域非極大值抑制
imshow("Non-Maximum Image", imageLocalMax);
Mat cannyImage;
cannyImage = Mat::zeros(imageLocalMax.size(), CV_8UC1);
DoubleThreshold(imageLocalMax, 90, 160); //雙閾值處理
imshow("Double Threshold Image", imageLocalMax);
DoubleThresholdLink(imageLocalMax, 90, 160); //雙閾值中間閾值濾除及連接
imshow("Canny Image", imageLocalMax);
}
-
Laplace算子
1812年，拉普拉斯（Laplace, 1749 – 1827）發表了重要的《概率分析理論》一書，在該書中總結了當時整個概率論的研究，論述了概率在選舉審判調查、氣象等方面的應用，匯入「拉普拉斯變換」等，
拉普拉斯(Laplacian) 算子是 維歐幾里德空間中的一個二階微分算子，常用于影像增強領域和邊緣提取，它通過灰度差分計算鄰域內的像素，
演算法基本流程
1）判斷影像中心像素灰度值與它周圍其他像素的灰度值，如果中心像素的灰度更高，則提升中心像素的灰度；反之降低中心像素的灰度，從而實作影像銳化操作；
2）在演算法實作程序中，Laplacian算子通過對鄰域中心像素的四方向或八方向求梯度，再將梯度相加起來判斷中心像素灰度與鄰域內其他像素灰度的關系；
3）最后通過梯度運算的結果對像素灰度進行調整，
Laplacian算子分為四鄰域和八鄰域，四鄰域是對鄰域中心像素的四個方向求梯度，八鄰域是對八個方向求梯度，
其中，Laplacian算子四鄰域模板如下所示：
Laplacian算子的八鄰域模板如下所示：
通過Laplacian算子的模板可以發現：
1）當鄰域內像素灰度相同時，模板的卷積運算結果為0；
2）當中心像素灰度高于鄰域內其他像素的平均灰度時，模板的卷積運算結果為正數；
3）當中心像素的灰度低于鄰域內其他像素的平均灰度時，模板的卷積為負數，對卷積運算的結果用適當的衰弱因子處理并加在原中心像素上，就可以實作影像的銳化處理，
void quick_opencv::laplacian_Demo(Mat &img)
{
Mat src_gray;
int kernel_size = 3;
const char* window_name = "Laplacian Demo";
cvtColor(img, src_gray, COLOR_RGB2GRAY);
namedWindow(window_name, WINDOW_AUTOSIZE);
Mat dst, abs_dst;
Laplacian(src_gray, dst, CV_16S, kernel_size);
convertScaleAbs(dst, abs_dst);
imshow(window_name, abs_dst);
}
-
LoG算子
1980年，LoG邊緣檢測算子是David Courtnay Marr和Ellen Hildreth共同提出的，因此，也稱為Marr & Hildreth 邊緣檢測演算法或Marr & Hildreth算子，該演算法首先對影像做高斯濾波，然后再求其拉普拉斯（Laplacian）二階導數，即影像與 Laplacian of the Gaussian function 進行濾波運算，最后，可以通過檢測濾波結果的零交叉（Zero crossings）獲得影像或物體的邊緣，因而，也被業界簡稱為Laplacian-of-Gaussian (LoG)算子，
Log算子的運算式如下：
常用的卷積模板是5*5的模板：
void quick_opencv::LOG_Demo(Mat &img)
{
//高斯-拉普拉斯算子 二階微分 用于邊緣檢測
Mat src_gray;
int kernel_size = 3;
const char* window_name = "LOG Demo";
GaussianBlur(img, img, Size(3, 3), 0, 0, BORDER_DEFAULT);
cvtColor(img, src_gray, COLOR_RGB2GRAY);
namedWindow(window_name, WINDOW_AUTOSIZE);
Mat dst, abs_dst;
Laplacian(src_gray, dst, CV_16S, kernel_size);
convertScaleAbs(dst, abs_dst);
imshow(window_name, abs_dst);
}
-
DoG算子
1980年，高斯差分(DoG)算子，Marr and Hildreth指出，使用高斯差分(DOG)來近似Log算子是可能的：
LoG算子和DoG算子的函式波形對比如下圖所示，由于高斯差分的計算更加簡單，因此可用DoG算子近似替代LoG算子：
//x，y方向聯合實作獲取高斯模板
void generateGaussMask(Mat& Mask, Size wsize, double sigma) {
Mask.create(wsize, CV_64F);
int h = wsize.height;
int w = wsize.width;
int center_h = (h - 1) / 2;
int center_w = (w - 1) / 2;
double sum = 0.0;
double x, y;
for (int i = 0; i < h; ++i) {
y = pow(i - center_h, 2);
for (int j = 0; j < w; ++j) {
x = pow(j - center_w, 2);
//因為最后都要歸一化的，常數部分可以不計算，也減少了運算量
double g = exp(-(x + y) / (2 * sigma*sigma));
Mask.at<double>(i, j) = g;
sum += g;
}
}
Mask = Mask / sum;
}
//按二維高斯函式實作高斯濾波
void GaussianFilter(Mat& src, Mat& dst, Mat window) {
int hh = (window.rows - 1) / 2;
int hw = (window.cols - 1) / 2;
dst = Mat::zeros(src.size(), src.type());
//邊界填充
Mat Newsrc;
copyMakeBorder(src, Newsrc, hh, hh, hw, hw, BORDER_REPLICATE);//邊界復制
//高斯濾波
for (int i = hh; i < src.rows + hh; ++i) {
for (int j = hw; j < src.cols + hw; ++j) {
double sum[3] = { 0 };
for (int r = -hh; r <= hh; ++r) {
for (int c = -hw; c <= hw; ++c) {
if (src.channels() == 1) {
sum[0] = sum[0] + Newsrc.at<uchar>(i + r, j + c) * window.at<double>(r + hh, c + hw);
}
else if (src.channels() == 3) {
Vec3b rgb = Newsrc.at<Vec3b>(i + r, j + c);
sum[0] = sum[0] + rgb[0] * window.at<double>(r + hh, c + hw);//B
sum[1] = sum[1] + rgb[1] * window.at<double>(r + hh, c + hw);//G
sum[2] = sum[2] + rgb[2] * window.at<double>(r + hh, c + hw);//R
}
}
}
for (int k = 0; k < src.channels(); ++k) {
if (sum[k] < 0)
sum[k] = 0;
else if (sum[k] > 255)
sum[k] = 255;
}
if (src.channels() == 1)
{
dst.at<uchar>(i - hh, j - hw) = static_cast<uchar>(sum[0]);
}
else if (src.channels() == 3)
{
Vec3b rgb = { static_cast<uchar>(sum[0]), static_cast<uchar>(sum[1]), static_cast<uchar>(sum[2]) };
dst.at<Vec3b>(i - hh, j - hw) = rgb;
}
}
}
}
//DOG高斯差分
void DOG1(Mat &src, Mat &dst, Size wsize, double sigma, double k = 1.6) {
Mat Mask1, Mask2, gaussian_dst1, gaussian_dst2;
generateGaussMask(Mask1, wsize, k*sigma);//獲取二維高斯濾波模板1
generateGaussMask(Mask2, wsize, sigma);//獲取二維高斯濾波模板2
//高斯濾波
GaussianFilter(src, gaussian_dst1, Mask1);
GaussianFilter(src, gaussian_dst2, Mask2);
dst = gaussian_dst1 - gaussian_dst2 - 1;
threshold(dst, dst, 0, 255, THRESH_BINARY);
}
//DOG高斯差分--使用opencv的GaussianBlur
void DOG2(Mat &src, Mat &dst, Size wsize, double sigma, double k = 1.6) {
Mat gaussian_dst1, gaussian_dst2;
//高斯濾波
GaussianBlur(src, gaussian_dst1, wsize, k*sigma);
GaussianBlur(src, gaussian_dst2, wsize, sigma);
dst = gaussian_dst1 - gaussian_dst2;
threshold(dst, dst, 0, 255, cv::THRESH_BINARY);
}
// 實作
void quick_opencv::DOG_Demo(Mat &img)
{
if (img.channels() > 1) cv::cvtColor(img, img, COLOR_BGR2GRAY);
Mat edge1, edge2;
DOG1(img, edge1, cv::Size(7, 7), 2);
DOG2(img, edge2, cv::Size(7, 7), 2);
namedWindow("My_DOG", WINDOW_NORMAL);
imshow("My_DOG", edge1);
namedWindow("Opencv_DOG", WINDOW_NORMAL);
imshow("Opencv_DOG", edge2);
}
轉載請註明出處，本文鏈接：https://www.uj5u.com/qita/295285.html
標籤：其他
上一篇：【深度好文】Python影像處理二值影像投影量的計算和應用
下一篇：百度影像識別（3）垃圾分類應用