1四、OpenCV實現圖像的空間濾波——圖像銳化及邊緣檢測 圖像處理基礎(6):銳化空間濾波器 OpenCV探索之路(六):邊緣檢測(canny、sobel、laplacian) OpenCV-跟我一

一、圖像銳化理論基礎

一、銳化的概念

    圖像銳化的目的是使模糊的圖像變得清晰起來,主要用於加強圖像的灰度跳變部分,這一點與圖像平滑對灰度跳變的抑制正好相反。並且從算子能夠看出來,平滑是基於對圖像領域的加權求和或者說積分運算的,而銳化則是經過其逆運算導數(梯度)或者說有限差分來實現的。html

二、圖像的一階微分和二階微分的性質

圖像的銳化也就是加強圖像的突變部分,那麼咱們也就對圖像的恆定區域中,突變的開始點與結束點(臺階和斜坡突變)及沿着灰度斜坡處的微分的性質。微分是對函數局部變化率的一種表示,那麼對於一階微分有如下幾個性質:ios

  • 在恆定的灰度區域,圖像的微分值爲0.(灰度值沒有發生變換,天然微分爲0)
  • 在灰度臺階或斜坡起點處微分值不爲0.(臺階是,灰度值的突變變化較大;斜坡則是灰度值變化較緩慢;灰度值發生了變化,微分值不爲0)
  • 沿着斜坡的微分值不爲0.

二階微分,是一階微分的導數,和一階微分相對應,也有如下幾點性質:算法

  • 在恆定區域二階微分值爲0
  • 在灰度臺階或斜坡的起點處微分值不爲0
  • 沿着斜坡的微分值爲0.

從以上圖像灰度的一階和二階微分的性質能夠看出,在灰度值變化的地方,一階微分和二階微分的值都不爲0;在灰度恆定的地方,微分值都爲0.也就是說,不管是使用一階微分仍是二階微分均可以獲得圖像灰度的變化值。app

三、圖像的一階微分和二階微分的計算

圖像能夠看着是二維離散函數,對於圖像的一階微分其計算公式以下:dom

\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}ide

對於二階微分有:函數

\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}post

 在圖像處理中的一階微分一般使用梯度的幅值來實現,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}
$$

二、基於一階導的梯度算子

 一、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;
}

 

二、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大小,必須是一、三、五、七、

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

}

  程序運行結果以下:

 

 

三、基於二階微分的算子

一、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 );

參數意義爲,

  1. src_gray,輸入圖像
  2. dst,Laplace操做結果
  3. ddepth,輸出圖像深度,由於輸入圖像通常爲CV_8U,爲了不數據溢出,輸出圖像深度應該設置爲CV_16S
  4. 濾波器孔徑尺寸;
  5. 比例因子;
  6. 表示結果存入目標圖

 示例

#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);

}

  

 

二、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;

}
  

 

四、導向濾波

  導向圖濾波主要是經過一張引導圖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;
}

 

 五、canny算子——邊緣檢測方法的提出

   Canny算子是John F. Canny於 1986 年開發出來的一個多級邊緣檢測算法。Canny 的目標是找到一個最優的邊緣檢測算法,最優邊緣檢測的含義是:

(1)最優檢測:算法可以儘量多地標識出圖像中的實際邊緣,漏檢真實邊緣的機率和誤檢非邊緣的機率都儘量小;
(2)最優定位準則:檢測到的邊緣點的位置距離實際邊緣點的位置最近,或者是因爲噪聲影響引發檢測出的邊緣偏離物體的真實邊緣的程度最小;
(3)檢測點與邊緣點一一對應:算子檢測的邊緣點與實際邊緣點應該是一一對應。
爲了知足這些要求 Canny 使用了變分法(calculus of variations),這是一種尋找優化特定功能的函數的方法。最優檢測使用四個指數函數項表示,可是它很是近似於高斯函數的一階導數。
 
Canny邊緣檢測算法能夠分爲如下5個步驟:
  1. 應用高斯濾波來平滑圖像,目的是去除噪聲,任何邊緣檢測算法都不可能在未經處理的原始數據上很好地工做,因此第一步是對原始數據與高斯 mask 做 卷積,獲得的圖像與原始圖像相比有些輕微的模糊(blurred)。
  2. 找尋圖像的強度梯度(intensity gradients),即找尋一幅圖像中灰度強度變化最強的位置。
  3. 應用非最大抑制(non-maximum suppression)技術來消除邊誤檢(原本不是但檢測出來是)。就是保留了每一個像素點上梯度強度的極大值,而刪掉其餘的值。
  4. 應用雙閾值的方法來決定可能的(潛在的)邊界
  5. 利用滯後技術來跟蹤邊界

 非最大抑制

  非最大抑制的目的是將模糊(blurred)的邊界變得清晰(sharp)。通俗的講,就是保留了每一個像素點上梯度強度的極大值,而刪掉其餘的值。對於每一個像素點,進行以下操做:

a) 將其梯度方向近似爲如下值中的一個(0,45,90,135,180,225,270,315)(即上下左右和45度方向)
b) 比較該像素點,和其梯度方向正負方向的像素點的梯度強度
c) 若是該像素點梯度強度最大則保留,不然抑制(刪除,即置爲0)
爲了更好的解釋這個概念,看下圖。

圖中的數字表明瞭像素點的梯度強度,箭頭方向表明了梯度方向。以第二排第三個像素點爲例,因爲梯度方向向上,則將這一點的強度(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;
}

  

 

參考資料:

圖像處理基礎(6):銳化空間濾波器

邊緣檢測之Robert算子

【數字圖像處理】5.8:灰度圖像-圖像加強 Robert算子、Sobel算子

 

 

圖像銳化(加強)和邊緣檢測

OpenCV探索之路(六):邊緣檢測(canny、sobel、laplacian)

圖像邊緣檢測之拉普拉斯(Laplacian)C++實現

OpenCV-跟我一塊兒學數字圖像處理之拉普拉斯算子

【OpenCV圖像處理入門學習教程四】基於LoG算子的圖像邊緣檢測

LOG邊緣檢測--Marr-Hildreth邊緣檢測算法

圖像邊緣檢測——二階微分算子(上)Laplace算子、LOG算子、DOG算子(Matlab實現)

圖像邊緣檢測——二階微分算子(下)Canny算子(Matlab實現)

邊緣檢測之Canny

canny邊緣檢測算法原理與C語言實現

Canny算子邊緣檢測詳細原理(OpenCV+MATLAB實現)

【OpenCV圖像處理】十7、圖像的導向濾波

相關文章
相關標籤/搜索