13、OpenCV實現圖像的空間濾波——圖像平滑

1、空間濾波基礎概念

1、空間濾波基礎

空間濾波一詞中濾波取自數字信號處理，指接受或拒絕一定的頻率成分，但是空間濾波學習內容實際上和通過傅里葉變換實現的頻域的濾波是等效的，故而也稱為濾波。空間濾波主要直接基於領域（空間域）對圖像中的像素執行計算，用濾波器（也成為空間掩膜、核、模板和窗口）直接作用於圖像本身完成類似的平滑。

2、空間濾波機理

對空間域中的每一點(x,y)，重復如下操作：

• 對預先定義的以（x,y)為中心的領域內的像素進行預定義運算。
• 將（1）中運算的結果作為(x,y)點新的響應。

上述過程就稱為鄰域處理或空間域濾波。

3、空間濾波分類

根據預定義的操作，可以將濾波器分為：

• 線性濾波器
• 非線性濾波器

而根據濾波器最終對圖像造成的影響，可以將濾波器分為：

• 平滑濾波器 ，通常用於模糊圖像或者去除圖像中的噪聲
• 銳化濾波器，突出圖像中的邊緣細節部分

4、空間相關和卷積

 　在執行線性空間vlbo時，有兩個相近的概念需要理解，相關和卷積。相關是濾波器模板移過圖像並計算每個位置乘積之和的處理。卷積機理類似，但是濾波器首先需要旋轉180°。關於卷積的詳細說明，可以參考文章

卷積在圖像處理中的應用（轉自https://medium.com/@irhumshafkat/intuitively-understanding-convolutions-for-deep-learning-1f6f42faee1）

一些算子及不同類型的卷積

這里不再做重復說明。

2、為圖片添加噪聲

使用攝像頭拍攝圖片，噪聲是難以避免的，這里重點提下兩種噪聲：

• 椒鹽噪聲：噪聲的幅值基本上相同，但是噪聲出現的位置是隨機的；（中值濾波效果好）
• 高斯噪聲：每一點都存在噪聲，但噪聲的幅值是隨機分布的。

但是，我們網上找到的圖片大多是經過處理的，或者說噪聲不是明顯的，這里為了研究，為圖片添加噪聲。

椒鹽噪聲表現在圖片上是圖片中離散分布的黑點和白點

// 添加椒鹽噪聲
void addSaltNoise(Mat &m, int num)
{
    // 隨機數產生器
    std::random_device rd; //種子
    std::mt19937 gen(rd()); // 隨機數引擎

    auto cols = m.cols * m.channels();

    for (int i = 0; i < num; i++)
    {
        auto row = static_cast<int>(gen() % m.rows);
        auto col = static_cast<int>(gen() % cols);

        auto p = m.ptr<uchar>(row);
        p[col++] = 255;
        p[col++] = 255;
        p[col] = 255;
    }
}

高斯噪聲是一種加性噪聲，為圖像添加高斯噪聲的代碼如下：

// 添加Gussia噪聲
// 使用指針訪問
void addGaussianNoise(Mat &m, int mu, int sigma)
{
    // 產生高斯分布隨機數發生器
    std::random_device rd;
    std::mt19937 gen(rd());

    std::normal_distribution<> d(mu, sigma);

    auto rows = m.rows; // 行數
    auto cols = m.cols * m.channels(); // 列數

    for (int i = 0; i < rows; i++)
    {
        auto p = m.ptr<uchar>(i); // 取得行首指針
        for (int j = 0; j < cols; j++)
        {
            auto tmp = p[j] + d(gen);
            tmp = tmp > 255 ? 255 : tmp;
            tmp = tmp < 0 ? 0 : tmp;
            p[j] = tmp;
        }
    }
}

示例：

 　　　為圖像添加兩種噪聲的程序如下：

//實現直方圖的反向投影
#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 addSaltNoise(Mat &m, int num)
{
	// 隨機數產生器
	std::random_device rd; //種子
	std::mt19937 gen(rd()); // 隨機數引擎

	auto cols = m.cols * m.channels();

	for (int i = 0; i < num; i++)
	{
		auto row = static_cast<int>(gen() % m.rows);
		auto col = static_cast<int>(gen() % cols);

		auto p = m.ptr<uchar>(row);
		p[col++] = 255;
		p[col++] = 255;
		p[col] = 255;
	}
}
void addGaussianNoise(Mat &m, int mu, int sigma)
{
	// 產生高斯分布隨機數發生器
	std::random_device rd;
	std::mt19937 gen(rd());

	std::normal_distribution<> d(mu, sigma);

	auto rows = m.rows; // 行數
	auto cols = m.cols * m.channels(); // 列數

	for (int i = 0; i < rows; i++)
	{
		auto p = m.ptr<uchar>(i); // 取得行首指針
		for (int j = 0; j < cols; j++)
		{
			auto tmp = p[j] + d(gen);
			tmp = tmp > 255 ? 255 : tmp;
			tmp = tmp < 0 ? 0 : tmp;
			p[j] = tmp;
		}
	}
}

int main()
{
	Mat srcImage = imread("111.jpg");
	Mat src_salt, src_gaussian;
	src_salt = imread("111.jpg");
	src_gaussian = imread("111.jpg");
	if (!srcImage.data)
	{
		cout << "圖像打開失敗！" << endl;
		return -1;
	}
	addSaltNoise(src_salt, 100);
	addGaussianNoise(src_gaussian,100,1);
	imshow("原圖", srcImage);
	imshow("添加椒鹽噪聲后", src_salt);
	imshow("添加高斯噪聲后", src_gaussian);
	waitKey();
	return 0;
}

添加完噪聲后的圖像如下：

 

3、線性平滑空間濾波

平滑濾波器主要用於模糊處理和降低噪音。

1、平滑線性濾波器(又稱均值濾波）

　　首先需要說明的是圖像平滑是一種減少和抑制圖像噪聲的實用數字圖像處理技術，一般來說，圖像具有局部連續性質，即相鄰像素的數值相近，而噪聲的存在使得在噪聲點處產生灰度跳躍，但一般可以合理的假設偶爾出現的噪聲並沒有改變圖像局部連續的性質。其中，線性濾波是現實將空域模板對應鄰域內像素的灰度值加權之和作為鄰域內中心像素的相應輸出。線性平滑模板的權系數全為正值而且系數之和等於1，如下圖所示：

 

　　因此這種方法不會增加圖像中總體的灰度程度。也就是說在灰度一致的區域，線性平滑濾波的相應輸出不變。此外，可以知道，線性平滑濾波等效於低通濾波。但是同樣的我們可以看出，線性平滑模板在模糊和降噪的同時，圖像中的邊緣和細節的銳度也都丟失了，也就是說，在平滑的過程中，使得圖像的一部分細節信息丟失。

 

在OpenCV中，函數blur表示使用該模板的均值濾波器，其聲明如下：

void blur( InputArray src, OutputArray dst, Size ksize, Point anchor = Point(-1,-1), int borderType = BORDER_DEFAULT );

src是輸入圖像，dst為輸出圖像；ksize是濾波器模板窗口的大小；后兩個參數分別表示，待處理像素在模板窗口的位置，默認值是窗口的中心位置，所以窗口的大小一般為奇數，最后一個參數表示對編解類型的處理，使用默認值即可。其調用示例blur(src,dst,Size(5,5)，模板窗口的大小為5×55×5。

示例：

下面通過程序對兩種含噪聲的圖片進行線性平滑濾波處理

//實現直方圖的反向投影
#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 addSaltNoise(Mat &m, int num)
{
	// 隨機數產生器
	std::random_device rd; //種子
	std::mt19937 gen(rd()); // 隨機數引擎

	auto cols = m.cols * m.channels();

	for (int i = 0; i < num; i++)
	{
		auto row = static_cast<int>(gen() % m.rows);
		auto col = static_cast<int>(gen() % cols);

		auto p = m.ptr<uchar>(row);
		p[col++] = 255;
		p[col++] = 255;
		p[col] = 255;
	}
}
void addGaussianNoise(Mat &m, int mu, int sigma)
{
	// 產生高斯分布隨機數發生器
	std::random_device rd;
	std::mt19937 gen(rd());

	std::normal_distribution<> d(mu, sigma);

	auto rows = m.rows; // 行數
	auto cols = m.cols * m.channels(); // 列數

	for (int i = 0; i < rows; i++)
	{
		auto p = m.ptr<uchar>(i); // 取得行首指針
		for (int j = 0; j < cols; j++)
		{
			auto tmp = p[j] + d(gen);
			tmp = tmp > 255 ? 255 : tmp;
			tmp = tmp < 0 ? 0 : tmp;
			p[j] = tmp;
		}
	}
}

int main()
{
	Mat srcImage = imread("111.jpg");
	Mat src_salt, src_gaussian;
	src_salt = imread("111.jpg");
	src_gaussian = imread("111.jpg");
	if (!srcImage.data)
	{
		cout << "圖像打開失敗！" << endl;
		return -1;
	}
	addSaltNoise(src_salt, 100);
	addGaussianNoise(src_gaussian,100,1);

    //平滑濾波
	Mat aft_mean_filter_salt, aft_mean_filter_gaussian;
	blur(src_salt, aft_mean_filter_salt, Size(5, 5), Point(-1, -1));
	blur(src_gaussian, aft_mean_filter_gaussian, Size(5, 5), Point(-1, -1));


	//增強圖像對比度
	//Mat kernel = (Mat_<int>(3, 3) << 0, -1, 0, -1, 5, -1, 0, -1, 0);
	//filter2D(src_salt, src_salt, src_salt.depth(), kernel, Point(-1, -1));

	imshow("原圖", srcImage);
	imshow("添加椒鹽噪聲后", src_salt);
	imshow("添加高斯噪聲后", src_gaussian);
	imshow("經過平滑濾波后含椒鹽噪聲的的結果", aft_mean_filter_salt);
	imshow("經過平滑濾波后含高斯噪聲的的結果", aft_mean_filter_gaussian);
	waitKey();
	return 0;
}

　　

 

2、高斯平滑

　　均值平滑對於領域內的像素一視同仁，為了減少平滑處理中的模糊，得到更自然的平滑效果，我們可以很自然的想到適當加大模板中心點的權重，隨着距離中心點的距離增加，權重迅速減小，從而可以確保中心點看起來更接近於與它距離更近的點，基於這樣的考慮得到的模板即為高斯模板。常用的3X3高斯模板如下所示：

　　高斯模板名字的由來是高斯函數，即二維正態分布密度函數，一個二維的高斯函數如下,

　　$$
　　h(x, y)=e^{-\frac{x^{2}+y^{2}}{2 \sigma^{2}}}
　　$$

 公式中(x,y)是點坐標，在圖像處理中可認為是整數，σ是標准差。

  對於窗口模板的大小為 (2K+1)X(2K+1)的矩陣，其(i，j）位置的元素值可如下確定：

高斯模板的生成方法如下所示：

void generateGaussianTemplate(double window[][11], int ksize, double sigma)
{
    static const double pi = 3.1415926;
    int center = ksize / 2; // 模板的中心位置，也就是坐標的原點
    double x2, y2;
    for (int i = 0; i < ksize; i++)
    {
        x2 = pow(i - center, 2);
        for (int j = 0; j < ksize; j++)
        {
            y2 = pow(j - center, 2);
            double g = exp(-(x2 + y2) / (2 * sigma * sigma));
            g /= 2 * pi * sigma;
            window[i][j] = g;
        }
    }
    double k = 1 / window[0][0]; // 將左上角的系數歸一化為1
    for (int i = 0; i < ksize; i++)
    {
        for (int j = 0; j < ksize; j++)
        {
            window[i][j] *= k;
        }
    }
}

　　需要一個二維數組，存放生成的系數（這里假設模板的最大尺寸不會超過11）；第二個參數是模板的大小（不要超過11）；第三個參數就比較重要了，是高斯分布的標准差。

生成的過程，首先根據模板的大小，找到模板的中心位置ksize/2。 然后就是遍歷，根據高斯分布的函數，計算模板中每個系數的值。

需要注意的是，最后歸一化的過程，使用模板左上角的系數的倒數作為歸一化的系數（左上角的系數值被歸一化為1），模板中的每個系數都乘以該值（左上角系數的倒數），然后將得到的值取整，就得到了整數型的高斯濾波器模板。

至於小數形式的生成也比較簡單，去掉歸一化的過程，並且在求解過程后，模板的每個系數要除以所有系數的和。具體代碼如下：

void generateGaussianTemplate(double window[][11], int ksize, double sigma)
{
    static const double pi = 3.1415926;
    int center = ksize / 2; // 模板的中心位置，也就是坐標的原點
    double x2, y2;
    double sum = 0;
    for (int i = 0; i < ksize; i++)
    {
        x2 = pow(i - center, 2);
        for (int j = 0; j < ksize; j++)
        {
            y2 = pow(j - center, 2);
            double g = exp(-(x2 + y2) / (2 * sigma * sigma));
            g /= 2 * pi * sigma;
            sum += g;
            window[i][j] = g;
        }
    }
    //double k = 1 / window[0][0]; // 將左上角的系數歸一化為1
    for (int i = 0; i < ksize; i++)
    {
        for (int j = 0; j < ksize; j++)
        {
            window[i][j] /= sum;
        }
    }
}

σ選取技巧

　　當標准差σ選取不同的值是，二維高斯函數形狀會發生很大的變化，

• 如果σ過小，偏離中心的所有像素權杖重將會非常小，相當於加權和響應基本不考慮領域像素的作用，這樣無法起到濾波平滑的效果。
• 如果σ過大，而領域相對較小，這樣在領域內高斯模板則退化為平均模板。

在matlab中，σ的默認值為0.5，在實際應用中，通常3X3的模板取σ為0.8左右，更大的模板可以適當增加σ的值。

在OpenCV中使用GuassianBlur()函數來實現高斯平滑，函數的聲明如下：

void GaussianBlur( InputArray src, OutputArray dst, Size ksize, double sigmaX, double sigmaY=0, int borderType=BORDER_DEFAULT );

　這里由於高斯函數的可分離性，尺寸較大的高斯濾波器可以分成兩步進行：首先將圖像在水平（豎直）方向與一維高斯函數進行卷積；然后將卷積后的結果在豎直（水平）方向使用相同的一維高斯函數得到的模板進行卷積運算。

你也可以用如下源碼進行實現：

void separateGaussianFilter(const Mat &src, Mat &dst, int ksize, double sigma)
{
    CV_Assert(src.channels()==1 || src.channels() == 3); // 只處理單通道或者三通道圖像
    // 生成一維的高斯濾波模板
    double *matrix = new double[ksize];
    double sum = 0;
    int origin = ksize / 2;
    for (int i = 0; i < ksize; i++)
    {
        // 高斯函數前的常數可以不用計算，會在歸一化的過程中給消去
        double g = exp(-(i - origin) * (i - origin) / (2 * sigma * sigma));
        sum += g;
        matrix[i] = g;
    }
    // 歸一化
    for (int i = 0; i < ksize; i++)
        matrix[i] /= sum;
    // 將模板應用到圖像中
    int border = ksize / 2;
    copyMakeBorder(src, dst, border, border, border, border, BorderTypes::BORDER_REFLECT);
    int channels = dst.channels();
    int rows = dst.rows - border;
    int cols = dst.cols - border;
    // 水平方向
    for (int i = border; i < rows; i++)
    {
        for (int j = border; j < cols; j++)
        {
            double sum[3] = { 0 };
            for (int k = -border; k <= border; k++)
            {
                if (channels == 1)
                {
                    sum[0] += matrix[border + k] * dst.at<uchar>(i, j + k); // 行不變，列變化；先做水平方向的卷積
                }
                else if (channels == 3)
                {
                    Vec3b rgb = dst.at<Vec3b>(i, j + k);
                    sum[0] += matrix[border + k] * rgb[0];
                    sum[1] += matrix[border + k] * rgb[1];
                    sum[2] += matrix[border + k] * rgb[2];
                }
            }
            for (int k = 0; k < channels; k++)
            {
                if (sum[k] < 0)
                    sum[k] = 0;
                else if (sum[k] > 255)
                    sum[k] = 255;
            }
            if (channels == 1)
                dst.at<uchar>(i, j) = static_cast<uchar>(sum[0]);
            else if (channels == 3)
            {
                Vec3b rgb = { static_cast<uchar>(sum[0]), static_cast<uchar>(sum[1]), static_cast<uchar>(sum[2]) };
                dst.at<Vec3b>(i, j) = rgb;
            }
        }
    }
    // 豎直方向
    for (int i = border; i < rows; i++)
    {
        for (int j = border; j < cols; j++)
        {
            double sum[3] = { 0 };
            for (int k = -border; k <= border; k++)
            {
                if (channels == 1)
                {
                    sum[0] += matrix[border + k] * dst.at<uchar>(i + k, j); // 列不變，行變化；豎直方向的卷積
                }
                else if (channels == 3)
                {
                    Vec3b rgb = dst.at<Vec3b>(i + k, j);
                    sum[0] += matrix[border + k] * rgb[0];
                    sum[1] += matrix[border + k] * rgb[1];
                    sum[2] += matrix[border + k] * rgb[2];
                }
            }
            for (int k = 0; k < channels; k++)
            {
                if (sum[k] < 0)
                    sum[k] = 0;
                else if (sum[k] > 255)
                    sum[k] = 255;
            }
            if (channels == 1)
                dst.at<uchar>(i, j) = static_cast<uchar>(sum[0]);
            else if (channels == 3)
            {
                Vec3b rgb = { static_cast<uchar>(sum[0]), static_cast<uchar>(sum[1]), static_cast<uchar>(sum[2]) };
                dst.at<Vec3b>(i, j) = rgb;
            }
        }
    }
    delete[] matrix;
}

示例：

//實現直方圖的反向投影
#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 addSaltNoise(Mat &m, int num)
{
	// 隨機數產生器
	std::random_device rd; //種子
	std::mt19937 gen(rd()); // 隨機數引擎

	auto cols = m.cols * m.channels();

	for (int i = 0; i < num; i++)
	{
		auto row = static_cast<int>(gen() % m.rows);
		auto col = static_cast<int>(gen() % cols);

		auto p = m.ptr<uchar>(row);
		p[col++] = 255;
		p[col++] = 255;
		p[col] = 255;
	}
}
void addGaussianNoise(Mat &m, int mu, int sigma)
{
	// 產生高斯分布隨機數發生器
	std::random_device rd;
	std::mt19937 gen(rd());

	std::normal_distribution<> d(mu, sigma);

	auto rows = m.rows; // 行數
	auto cols = m.cols * m.channels(); // 列數

	for (int i = 0; i < rows; i++)
	{
		auto p = m.ptr<uchar>(i); // 取得行首指針
		for (int j = 0; j < cols; j++)
		{
			auto tmp = p[j] + d(gen);
			tmp = tmp > 255 ? 255 : tmp;
			tmp = tmp < 0 ? 0 : tmp;
			p[j] = tmp;
		}
	}
}

int main()
{
	Mat srcImage = imread("111.jpg");
	Mat src_salt, src_gaussian;
	src_salt = imread("111.jpg");
	src_gaussian = imread("111.jpg");
	if (!srcImage.data)
	{
		cout << "圖像打開失敗！" << endl;
		return -1;
	}
	addSaltNoise(src_salt, 100);
	addGaussianNoise(src_gaussian,100,1);

    //平滑濾波
	Mat aft_mean_filter_salt, aft_mean_filter_gaussian;
	blur(src_salt, aft_mean_filter_salt, Size(5, 5), Point(-1, -1));
	blur(src_gaussian, aft_mean_filter_gaussian, Size(5, 5), Point(-1, -1));


	//增強圖像對比度
	//Mat kernel = (Mat_<int>(3, 3) << 0, -1, 0, -1, 5, -1, 0, -1, 0);
	//filter2D(src_salt, src_salt, src_salt.depth(), kernel, Point(-1, -1));

	imshow("原圖", srcImage);
	imshow("添加椒鹽噪聲后", src_salt);
	imshow("添加高斯噪聲后", src_gaussian);
	imshow("經過平滑濾波后含椒鹽噪聲的的結果", aft_mean_filter_salt);
	imshow("經過平滑濾波后含高斯噪聲的的結果", aft_mean_filter_gaussian);
	waitKey();
	return 0;
}

　

 

4、非線性平滑空間濾波

1、中值濾波

　　中值濾波本質上是一種統計排序濾波器，對於原圖像中某點(i,j)，中值濾波以該點為中心的領域內的所有像素的統計排序中值作為(i,j)的響應，即它將每一像素點的灰度值設置為該點某鄰域窗口內的所有像素點灰度值的中值，也就是將中心像素的值用所有像素值的中間值（不是平均值）替換。中值濾波不同於均值，是指排序隊列中位於中間位置的元素的值。

　　中值濾波對於某些類型的隨機噪聲具有非常理想的降噪能力，對於線性平滑濾波而言，在處理的像素領域之內包含噪聲點時，噪聲的存在總會或多或少的影響該點的像素值的計算(對於高斯平滑，影響程度同噪聲點到中點的距離成正比），但在中值濾波中，噪聲點則常常是直接忽略掉的；而且在同線性平滑濾波器相比，中值濾波在降噪的同時引起的模糊效應較低。

　　OpenCV提供了中值濾波的API函數如下：

void medianBlur(InputArray src,OutputArray dst,int ksize)

　　參數解釋：原圖像，目標圖像，模板尺寸。

示例：

//實現直方圖的反向投影
#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 addSaltNoise(Mat &m, int num)
{
	// 隨機數產生器
	std::random_device rd; //種子
	std::mt19937 gen(rd()); // 隨機數引擎

	auto cols = m.cols * m.channels();

	for (int i = 0; i < num; i++)
	{
		auto row = static_cast<int>(gen() % m.rows);
		auto col = static_cast<int>(gen() % cols);

		auto p = m.ptr<uchar>(row);
		p[col++] = 255;
		p[col++] = 255;
		p[col] = 255;
	}
}
void addGaussianNoise(Mat &m, int mu, int sigma)
{
	// 產生高斯分布隨機數發生器
	std::random_device rd;
	std::mt19937 gen(rd());

	std::normal_distribution<> d(mu, sigma);

	auto rows = m.rows; // 行數
	auto cols = m.cols * m.channels(); // 列數

	for (int i = 0; i < rows; i++)
	{
		auto p = m.ptr<uchar>(i); // 取得行首指針
		for (int j = 0; j < cols; j++)
		{
			auto tmp = p[j] + d(gen);
			tmp = tmp > 255 ? 255 : tmp;
			tmp = tmp < 0 ? 0 : tmp;
			p[j] = tmp;
		}
	}
}

int main()
{
	Mat srcImage = imread("111.jpg");
	Mat src_salt, src_gaussian;
	src_salt = imread("111.jpg");
	src_gaussian = imread("111.jpg");
	if (!srcImage.data)
	{
		cout << "圖像打開失敗！" << endl;
		return -1;
	}
	addSaltNoise(src_salt, 100);
	addGaussianNoise(src_gaussian,100,1);

    //平滑濾波
	Mat aft_median_filter_salt, aft_median_filter_gaussian;
	medianBlur(src_salt, aft_median_filter_salt, 3);
	medianBlur(src_gaussian, aft_median_filter_gaussian, 3);
	//增強圖像對比度
	//Mat kernel = (Mat_<int>(3, 3) << 0, -1, 0, -1, 5, -1, 0, -1, 0);
	//filter2D(src_salt, src_salt, src_salt.depth(), kernel, Point(-1, -1));

	imshow("原圖", srcImage);
	imshow("添加椒鹽噪聲后", src_salt);
	imshow("添加高斯噪聲后", src_gaussian);
	imshow("經過平滑濾波后含椒鹽噪聲的的結果", aft_median_filter_salt);
	imshow("經過平滑濾波后含高斯噪聲的的結果", aft_median_filter_gaussian);
	waitKey();
	return 0;
}

　　

2、自適應中值濾波

 　　中值濾波的效果依賴於濾波窗口的大小，太大會使圖像邊緣模糊，太小則去噪效果不好。因為噪聲點和邊緣點同樣是灰度變化較為劇烈的像素，普通中值濾波在改變噪聲點灰度值的時候，會一定程度上改變邊緣像素的灰度值。但是噪聲點幾乎就是鄰近像素的極值，而邊緣往往不是，我們可以通過這個特性來限制中值濾波。

　　具體的改進方法如下：追行掃描圖像， 當處理每一個像素時，判斷該像素是否是濾波窗口所覆蓋下領域像素的極大或者極小值，如果是，則采用正常的中值濾波處理該像素，如果不是，則不予處理。這種方法能夠非常有效的去除突發噪聲點，尤其是椒鹽噪聲，而且幾乎不影響邊緣。這種方法又稱為自適應中值濾波。

　　算法實現如下所示：

uchar adaptiveProcess(const Mat &im, int row, int col, int kernelSize, int maxSize)
{
	vector<uchar> pixels;
	for (int a = -kernelSize / 2; a <= kernelSize / 2; a++)
		for (int b = -kernelSize / 2; b <= kernelSize / 2; b++)
		{
			pixels.push_back(im.at<uchar>(row + a, col + b));
		}
	sort(pixels.begin(), pixels.end());
	auto min = pixels[0];
	auto max = pixels[kernelSize * kernelSize - 1];
	auto med = pixels[kernelSize * kernelSize / 2];
	auto zxy = im.at<uchar>(row, col);
	if (med > min && med < max)
	{
		// to B
		if (zxy > min && zxy < max)
			return zxy;
		else
			return med;
	}
	else
	{
		kernelSize += 2;
		if (kernelSize <= maxSize)
			return adaptiveProcess(im, row, col, kernelSize, maxSize); // 增大窗口尺寸，繼續A過程。
		else
			return med;
	}
}

　示例：

//實現直方圖的反向投影
#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;
uchar adaptiveProcess(const Mat &im, int row, int col, int kernelSize, int maxSize);
// 添加椒鹽噪聲
void addSaltNoise(Mat &m, int num)
{
	// 隨機數產生器
	std::random_device rd; //種子
	std::mt19937 gen(rd()); // 隨機數引擎

	auto cols = m.cols * m.channels();

	for (int i = 0; i < num; i++)
	{
		auto row = static_cast<int>(gen() % m.rows);
		auto col = static_cast<int>(gen() % cols);

		auto p = m.ptr<uchar>(row);
		p[col++] = 255;
		p[col++] = 255;
		p[col] = 255;
	}
}
void addGaussianNoise(Mat &m, int mu, int sigma)
{
	// 產生高斯分布隨機數發生器
	std::random_device rd;
	std::mt19937 gen(rd());

	std::normal_distribution<> d(mu, sigma);

	auto rows = m.rows; // 行數
	auto cols = m.cols * m.channels(); // 列數

	for (int i = 0; i < rows; i++)
	{
		auto p = m.ptr<uchar>(i); // 取得行首指針
		for (int j = 0; j < cols; j++)
		{
			auto tmp = p[j] + d(gen);
			tmp = tmp > 255 ? 255 : tmp;
			tmp = tmp < 0 ? 0 : tmp;
			p[j] = tmp;
		}
	}
}

uchar adaptiveProcess(const Mat &im, int row, int col, int kernelSize, int maxSize)
{
	vector<uchar> pixels;
	for (int a = -kernelSize / 2; a <= kernelSize / 2; a++)
		for (int b = -kernelSize / 2; b <= kernelSize / 2; b++)
		{
			pixels.push_back(im.at<uchar>(row + a, col + b));
		}
	sort(pixels.begin(), pixels.end());
	auto min = pixels[0];
	auto max = pixels[kernelSize * kernelSize - 1];
	auto med = pixels[kernelSize * kernelSize / 2];
	auto zxy = im.at<uchar>(row, col);
	if (med > min && med < max)
	{
		// to B
		if (zxy > min && zxy < max)
			return zxy;
		else
			return med;
	}
	else
	{
		kernelSize += 2;
		if (kernelSize <= maxSize)
			return adaptiveProcess(im, row, col, kernelSize, maxSize); // 增大窗口尺寸，繼續A過程。
		else
			return med;
	}
}

int main()
{
	Mat srcImage = imread("111.jpg",0);
	Mat src_salt, src_gaussian;
	src_salt = imread("111.jpg",0);
	src_gaussian = imread("111.jpg",0);
	if (!srcImage.data)
	{
		cout << "圖像打開失敗！" << endl;
		return -1;
	}
	addSaltNoise(src_salt, 100);
	addGaussianNoise(src_gaussian,100,1);

	int minSize = 3; // 濾波器窗口的起始尺寸
	int maxSize = 7; // 濾波器窗口的最大尺寸
	//平滑濾波
	Mat aft_admedian_filter_salt, aft_admedian_filter_gaussian;
	// 擴展圖像的邊界
	copyMakeBorder(src_salt, aft_admedian_filter_salt, maxSize / 2, maxSize / 2, maxSize / 2, maxSize / 2, BorderTypes::BORDER_REFLECT);
	for (int j = maxSize / 2; j < aft_admedian_filter_salt.rows - maxSize / 2; j++)
	{
		for (int i = maxSize / 2; i < aft_admedian_filter_salt.cols * aft_admedian_filter_salt.channels() - maxSize / 2; i++)
		{
			aft_admedian_filter_salt.at<uchar>(j, i) = adaptiveProcess(aft_admedian_filter_salt, j, i, minSize, maxSize);
		}
	}

	copyMakeBorder(src_salt, aft_admedian_filter_gaussian, maxSize / 2, maxSize / 2, maxSize / 2, maxSize / 2, BorderTypes::BORDER_REFLECT);
	for (int j = maxSize / 2; j < aft_admedian_filter_gaussian.rows - maxSize / 2; j++)
	{
		for (int i = maxSize / 2; i < aft_admedian_filter_gaussian.cols * aft_admedian_filter_gaussian.channels() - maxSize / 2; i++)
		{
			aft_admedian_filter_gaussian.at<uchar>(j, i) = adaptiveProcess(aft_admedian_filter_gaussian, j, i, minSize, maxSize);
		}
	}

	medianBlur(src_gaussian, aft_admedian_filter_gaussian, 3);
	//增強圖像對比度
	//Mat kernel = (Mat_<int>(3, 3) << 0, -1, 0, -1, 5, -1, 0, -1, 0);
	//filter2D(src_salt, src_salt, src_salt.depth(), kernel, Point(-1, -1));

	imshow("原圖", srcImage);
	imshow("添加椒鹽噪聲后", src_salt);
	imshow("添加高斯噪聲后", src_gaussian);
	imshow("經過平滑濾波后含椒鹽噪聲的的結果", aft_admedian_filter_salt);
	imshow("經過平滑濾波后含高斯噪聲的的結果", aft_admedian_filter_gaussian);
	waitKey();
	return 0;
}

 

3、雙邊濾波

 　　雙邊濾波（Bilateral filter）是一種非線性的濾波方法，是結合圖像的空間鄰近度和像素值相似度的一種折衷處理，同時考慮空域信息和灰度相似性，達到保邊去噪的目的。雙邊濾波能夠提供一種不會將邊緣平滑掉的方法，但作為代價，需要更多的處理時間。與高斯濾波類似，雙邊濾波會依據每個像素及其領域構造一個加權平均值，加權計算包括兩個部分，其中第一部分加權方式與高斯平滑中相同，第二部分也屬於高斯加權，但不是基於中心像素點與其他像素點的空間距離之上的加權，而是基於其他像素與中心像素的亮度差值的加權。可以將雙邊濾波視為高斯平滑，對相似的像素賦予較高的權重，不相似的像素賦予較小的權重，也可用於圖像分割之中。

　　雙邊濾波器之所以能夠做到在平滑去噪的同時還能夠很好的保存邊緣（Edge Preserve），是由於其濾波器的核由兩個函數生成：

• 一個函數由像素歐式距離決定濾波器模板的系數
• 另一個函數由像素的灰度差值決定濾波器的系數　　

　　其綜合了高斯濾波器（Gaussian Filter）和α-截尾均值濾波器（Alpha-Trimmed mean Filter）的特點。高斯濾波器只考慮像素間的歐式距離，其使用的模板系數隨着和窗口中心的距離增大而減小；α-截尾均值濾波器則只考慮了像素灰度值之間的差值，去掉α%的最小值和最大值后再計算均值。

雙邊濾波器使用二維高斯函數生成距離模板，使用一維高斯函數生成值域模板。
距離模板系數的生成公式如下：

$$
d(i, j, k, l)=\exp \left(-\frac{(i-k)^{2}+(j-l)^{2}}{2 \sigma_{d}^{2}}\right)
$$

其中,(k,l)為模板窗口的中心坐標;(i,j)為模板窗口的其他系數的坐標；σ_d 為高斯函數的標准差。 使用該公式生成的濾波器模板和高斯濾波器使用的模板是沒有區別的。

值域模板系數的生成公式如下：

$$
r(i, j, k, l)=\exp \left(-\frac{\|f(i, j)-f(k, l)\|^{2}}{2 \sigma_{r}^{2}}\right)
$$

其中，函數f(x,y)表示要處理的圖像，f(x,y)表示圖像在點(x,y)處的像素值；(k,l)為模板窗口的中心坐標；(i,j)為模板窗口的其他系數的坐標；σ_r為高斯函數的標准差。

將上述兩個模板相乘就得到了雙邊濾波器的模板

$$
w(i, j, k, l)=d(i, j, k, l) * r(i, j, k, l)=\exp \left(-\frac{(i-k)^{2}+(j-l)^{2}}{2 \sigma_{d}^{2}}-\frac{\|f(i, j)-f(k, l)\|^{2}}{2 \sigma_{r}^{2}}\right)
$$

其中定義域濾波和值域濾波可以用下圖形象表示：

                                    

 　　代碼實現如下：

void myBilateralFilter(const Mat &src, Mat &dst, int ksize, double space_sigma, double color_sigma)
{
    int channels = src.channels();
    CV_Assert(channels == 1 || channels == 3);
    double space_coeff = -0.5 / (space_sigma * space_sigma);
    double color_coeff = -0.5 / (color_sigma * color_sigma);
    int radius = ksize / 2;
    Mat temp;
    copyMakeBorder(src, temp, radius, radius, radius, radius, BorderTypes::BORDER_REFLECT);
    vector<double> _color_weight(channels * 256); // 存放差值的平方
    vector<double> _space_weight(ksize * ksize); // 空間模板系數
    vector<int> _space_ofs(ksize * ksize); // 模板窗口的坐標
    double *color_weight = &_color_weight[0];
    double *space_weight = &_space_weight[0];
    int    *space_ofs = &_space_ofs[0];
    for (int i = 0; i < channels * 256; i++)
        color_weight[i] = exp(i * i * color_coeff);
    // 生成空間模板
    int maxk = 0;
    for (int i = -radius; i <= radius; i++)
    {
        for (int j = -radius; j <= radius; j++)
        {
            double r = sqrt(i*i + j * j);
            if (r > radius)
                continue;
            space_weight[maxk] = exp(r * r * space_coeff); // 存放模板系數
            space_ofs[maxk++] = i * temp.step + j * channels; // 存放模板的位置，和模板系數相對應
        }
    }
    // 濾波過程
    for (int i = 0; i < src.rows; i++)
    {
        const uchar *sptr = temp.data + (i + radius) * temp.step + radius * channels;
        uchar *dptr = dst.data + i * dst.step;
        if (channels == 1)
        {
            for (int j = 0; j < src.cols; j++)
            {
                double sum = 0, wsum = 0;
                int val0 = sptr[j]; // 模板中心位置的像素
                for (int k = 0; k < maxk; k++)
                {
                    int val = sptr[j + space_ofs[k]];
                    double w = space_weight[k] * color_weight[abs(val - val0)]; // 模板系數 = 空間系數 * 灰度值系數
                    sum += val * w;
                    wsum += w;
                }
                dptr[j] = (uchar)cvRound(sum / wsum);
            }
        }
        else if (channels == 3)
        {
            for (int j = 0; j < src.cols * 3; j+=3)
            {
                double sum_b = 0, sum_g = 0, sum_r = 0, wsum = 0;
                int b0 = sptr[j];
                int g0 = sptr[j + 1];
                int r0 = sptr[j + 2];
                for (int k = 0; k < maxk; k++)
                {
                    const uchar *sptr_k = sptr + j + space_ofs[k];
                    int b = sptr_k[0];
                    int g = sptr_k[1];
                    int r = sptr_k[2];
                    double w = space_weight[k] * color_weight[abs(b - b0) + abs(g - g0) + abs(r - r0)];
                    sum_b += b * w;
                    sum_g += g * w;
                    sum_r += r * w;
                    wsum += w;
                }
                wsum = 1.0f / wsum;
                b0 = cvRound(sum_b * wsum);
                g0 = cvRound(sum_g * wsum);
                r0 = cvRound(sum_r * wsum);
                dptr[j] = (uchar)b0;
                dptr[j + 1] = (uchar)g0;
                dptr[j + 2] = (uchar)r0;
            }
        }
    }
}

　　OpenCV提供了對應雙邊濾波的API接口

void bilateralFilter(InputArray src,OutputArray dst,int d,double sigmaColor,double sigmaSpace,int borderType=BORDER_DEFAULT)

第一個參數，InputArray類型的src，輸入圖像，即源圖像，需要為8為或者浮點型單通道、三通道的圖像；
第二個參數，OutputArray類型的dst，即目標圖像，需要和源圖像有一樣的尺寸和類型；
第三個參數，int類型的d，表示在過濾過程中每個像素鄰域的直徑。如果這個參數被設置為負值，那么OpenCV會從第五個參數sigmaSpace來計算出它；
第四個參數，double類型的sigmaColor，顏色空間濾波器的sigma值，這個參數的值越大，就表明該像素鄰域內有越寬廣的顏色會被混合到一起，產生較大的半相等顏色區域；
第五個參數，double類型的 sigmaSpace，坐標空間中的sigma值，坐標空間的標注方差。它的數值越大，意味着越遠的像素會相互影響，從而使更大區域中足夠相似的顏色獲取相同的顏色，當d>0時，d制定了鄰域大小且與sigmaSpace無關否則，d正比於sigmaSpace；
第六個參數，int類型的borderType，用於推斷圖像外部像素的某種邊界模式，有默認值BORDER_DEFAULT。

示例：

　

//實現直方圖的反向投影
#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 addSaltNoise(Mat &m, int num)
{
	// 隨機數產生器
	std::random_device rd; //種子
	std::mt19937 gen(rd()); // 隨機數引擎

	auto cols = m.cols * m.channels();

	for (int i = 0; i < num; i++)
	{
		auto row = static_cast<int>(gen() % m.rows);
		auto col = static_cast<int>(gen() % cols);

		auto p = m.ptr<uchar>(row);
		p[col++] = 255;
		p[col++] = 255;
		p[col] = 255;
	}
}
void addGaussianNoise(Mat &m, int mu, int sigma)
{
	// 產生高斯分布隨機數發生器
	std::random_device rd;
	std::mt19937 gen(rd());

	std::normal_distribution<> d(mu, sigma);

	auto rows = m.rows; // 行數
	auto cols = m.cols * m.channels(); // 列數

	for (int i = 0; i < rows; i++)
	{
		auto p = m.ptr<uchar>(i); // 取得行首指針
		for (int j = 0; j < cols; j++)
		{
			auto tmp = p[j] + d(gen);
			tmp = tmp > 255 ? 255 : tmp;
			tmp = tmp < 0 ? 0 : tmp;
			p[j] = tmp;
		}
	}
}
void myBilateralFilter(const Mat &src, Mat &dst, int ksize, double space_sigma, double color_sigma)
{
	int channels = src.channels();
	CV_Assert(channels == 1 || channels == 3);
	double space_coeff = -0.5 / (space_sigma * space_sigma);
	double color_coeff = -0.5 / (color_sigma * color_sigma);
	int radius = ksize / 2;
	Mat temp;
	copyMakeBorder(src, temp, radius, radius, radius, radius, BorderTypes::BORDER_REFLECT);
	vector<double> _color_weight(channels * 256); // 存放差值的平方
	vector<double> _space_weight(ksize * ksize); // 空間模板系數
	vector<int> _space_ofs(ksize * ksize); // 模板窗口的坐標
	double *color_weight = &_color_weight[0];
	double *space_weight = &_space_weight[0];
	int    *space_ofs = &_space_ofs[0];
	for (int i = 0; i < channels * 256; i++)
		color_weight[i] = exp(i * i * color_coeff);
	// 生成空間模板
	int maxk = 0;
	for (int i = -radius; i <= radius; i++)
	{
		for (int j = -radius; j <= radius; j++)
		{
			double r = sqrt(i*i + j * j);
			if (r > radius)
				continue;
			space_weight[maxk] = exp(r * r * space_coeff); // 存放模板系數
			space_ofs[maxk++] = i * temp.step + j * channels; // 存放模板的位置，和模板系數相對應
		}
	}
	// 濾波過程
	for (int i = 0; i < src.rows; i++)
	{
		const uchar *sptr = temp.data + (i + radius) * temp.step + radius * channels;
		uchar *dptr = dst.data + i * dst.step;
		if (channels == 1)
		{
			for (int j = 0; j < src.cols; j++)
			{
				double sum = 0, wsum = 0;
				int val0 = sptr[j]; // 模板中心位置的像素
				for (int k = 0; k < maxk; k++)
				{
					int val = sptr[j + space_ofs[k]];
					double w = space_weight[k] * color_weight[abs(val - val0)]; // 模板系數 = 空間系數 * 灰度值系數
					sum += val * w;
					wsum += w;
				}
				dptr[j] = (uchar)cvRound(sum / wsum);
			}
		}
		else if (channels == 3)
		{
			for (int j = 0; j < src.cols * 3; j += 3)
			{
				double sum_b = 0, sum_g = 0, sum_r = 0, wsum = 0;
				int b0 = sptr[j];
				int g0 = sptr[j + 1];
				int r0 = sptr[j + 2];
				for (int k = 0; k < maxk; k++)
				{
					const uchar *sptr_k = sptr + j + space_ofs[k];
					int b = sptr_k[0];
					int g = sptr_k[1];
					int r = sptr_k[2];
					double w = space_weight[k] * color_weight[abs(b - b0) + abs(g - g0) + abs(r - r0)];
					sum_b += b * w;
					sum_g += g * w;
					sum_r += r * w;
					wsum += w;
				}
				wsum = 1.0f / wsum;
				b0 = cvRound(sum_b * wsum);
				g0 = cvRound(sum_g * wsum);
				r0 = cvRound(sum_r * wsum);
				dptr[j] = (uchar)b0;
				dptr[j + 1] = (uchar)g0;
				dptr[j + 2] = (uchar)r0;
			}
		}
	}
}
int main()
{
	Mat srcImage = imread("111.jpg");
	Mat src_salt, src_gaussian;
	src_salt = imread("111.jpg");
	src_gaussian = imread("111.jpg");
	if (!srcImage.data)
	{
		cout << "圖像打開失敗！" << endl;
		return -1;
	}
	addSaltNoise(src_salt, 100);
	addGaussianNoise(src_gaussian,100,1);

	//平滑濾波
	Mat aft_bilateral_filter_salt, aft_bilateral_filter_gaussian, srcImage1;

	bilateralFilter(srcImage, srcImage1, 5, 100, 3);
	bilateralFilter(src_salt, aft_bilateral_filter_salt, 5, 100, 3);
	bilateralFilter(src_gaussian, aft_bilateral_filter_gaussian, 5, 100, 3);
	//增強圖像對比度
	Mat kernel = (Mat_<int>(3, 3) << 0, -1, 0, -1, 5, -1, 0, -1, 0);
	filter2D(srcImage1, srcImage1, srcImage1.depth(), kernel, Point(-1, -1));

	imshow("原圖", srcImage);
	imshow("雙邊濾波處理過后的原圖", srcImage1);
	imshow("添加椒鹽噪聲后", src_salt);
	imshow("添加高斯噪聲后", src_gaussian);
	imshow("經過平滑濾波后含椒鹽噪聲的的結果", aft_bilateral_filter_salt);
	imshow("經過平滑濾波后含高斯噪聲的的結果", aft_bilateral_filter_gaussian);
	waitKey();
	return 0;
}

　　程序運行結果如下，雙邊濾波很好的對原圖像進行了優化，但是對應椒鹽噪聲之類的高頻噪聲沒有濾波效果。

 

參考資料：

圖像處理基礎(1)：噪聲的添加和過濾

圖像處理基礎(2)：自適應中值濾波器(基於OpenCV實現)

圖像處理基礎(3)：均值濾波器及其變種

圖像處理基礎(4)：高斯濾波器詳解

圖像處理基礎(5)：雙邊濾波器

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

OpenCV--Python 圖像平滑之高斯平滑、均值平滑

【OpenCV圖像處理】十五、圖像空域濾波（上）

【OpenCV圖像處理】十六、圖像空域濾波（下）

【OpenCV】鄰域濾波：方框、高斯、中值、雙邊濾波

【OpenCV學習筆記】之圖像平滑（線性/非線性濾波器）

[Python圖像處理] 四.圖像平滑之均值濾波、方框濾波、高斯濾波及中值濾波

C++圖像處理 -- 表面模糊

選擇性模糊及其算法的實現。