学习guide滤波算法

前几天，老大给分陪了新的项目，开始和同事做滤波增强。开始几天看了一点相对好的算法，包括双边滤波，还有一个基于SVM学习的滤波算法（具体名字忘记了），然后就是同事推荐的导向滤波，查了下，确实感觉不错，能够良好的保持边缘特性。

基本思想大致是这样的，通过函数上的某一点与其邻近部分的点成线性关系，一个复杂的函数就可以用很多局部的线性函数来表示，当需求该函数上某一点时，只需要计算包含该店的线性函数的值并做平均即可。这种模型在表示非解析函数上非常有用，同理我们可以认为图像是一个二维函数。

        函数输入输出关系（二维窗口内），满足

q=aI+b  其中q是输出像素的值，I是输入像素的值。   

梯度q=a*梯度I

这个公式解释了输出具有边缘保持效果的原因。

剩下的就是求线性函数系数的回归运算，即希望拟合值与真实值P之间的差距最小

        E（a，b）=sum（（aI+b-p））^2 + e*a^2)

        p是待滤波图像并不像I可以是任意其他的图像。

最终求得系数为：

a=（1/|w| * sum （I*P）-u*P）/（var^2+e）

        b=（mean（p）-a/u）

最终的输出q=mean（a）*mean（I）+mean（b）

  |w|是窗口中像素的数量    var^2：方差  u：均值   mean（p）是滤波图像P在窗口中的均值

大体上就是不断计算，zai计算，求取均值，方差协方差等等吧

应项目需求，参考大神们的足迹，也实现了算法，要学的东西很多，希望毕业前做一个，记住一个，积累吧，还得学习编码。据说是 然要有梦想 万一实现了呢，逗比一下自己。

       

/************************************************************************
*  @Author : YG
*  @  Date : 2015/12/30
*  @ Email : 1851506****@163.com
*  function: 引导滤波
*  introduc：引导滤波是最近比较好的，前沿性的滤波方法，利用了线性模型，不仅能够平滑图像，而且保留了边缘特性，通过函数上一点与其领域部分点的线性关系，一个复杂的函数就可以用很多局部的线性函数来表示，当需求该函数上某一点的值时，只计算包含该点的线性函数的值并做平均即可，这种模型在表示解析函数上非常有效。e是尺度参数，平滑参数，为全局参数，大小影响整个分量的估计结果，尤其是细节比较丰富的图像。
************************************************************************/

#include "stdafx.h"
#include "guilt_me.h"
#include <iostream>
#include <windows.h>

using namespace std;
using namespace cv;


//灰度图处理
/*Mat img_do(Mat &img)
{
	if (!img.data)
	{
		cout << "no data of img1" << endl;
	}
	Mat img_gray(img.size(), CV_8UC1);
	cvtColor(img, img_gray, CV_BGR2GRAY);

	return img_gray;

}*/
/************************************************************************/
/*调节半径r，回调 */
/************************************************************************/
void on_guilt_feliter(int, void*)
{
	vector<Mat> bgr_src, bgr_dst;
	split(img1, bgr_src);
	Mat dst_color;
	double time;
	time = (double)GetTickCount();
	Mat q;
	for (int i = 0; i < 3; i++)
	{
		Mat I = img_Get(bgr_src[i]);
		Mat p = I.clone();
		q = guidedfilter(p, I, r, (double)e / 100);
		bgr_dst.push_back(q);
	}
	merge(bgr_dst, dst_color);

	imwrite("filtered.JPG", q * 255);
	time = 1000 * ((double)getTickCount() - time) / getTickFrequency();

	cout << endl << "Time of guided filter for  runs: " << time << " milliseconds." << endl;

	imshow("效果图", dst_color);
	imwrite("result.jpg", dst_color * 255);

	waitKey(0);

}
/************************************************************************/
/* 调节参数e，进行回调*/
/************************************************************************/
void on_guilt_feliter_e(int, void*)
{
	vector<Mat> bgr_src, bgr_dst;
	split(img1, bgr_src);
	Mat dst_color;
	double time;
	time = (double)GetTickCount();
	Mat q;
	for (int i = 0; i < 3; i++)
	{
		Mat I = img_Get(bgr_src[i]);
		Mat p = I.clone();
		q = guidedfilter(p, I, r, (double)e / 100);
		bgr_dst.push_back(q);
	}
	merge(bgr_dst, dst_color);

	imwrite("filtered.JPG", q * 255);
	time = 1000 * ((double)getTickCount() - time) / getTickFrequency();

	cout << endl << "Time of guided filter for  runs: " << time << " milliseconds." << endl;

	imshow("效果图", dst_color);
	imwrite("result.jpg", dst_color * 255);
	while (char(waitKey(1) != 'o'))
	{
	}

}
/************************************************************************/
/* 无模板，自身为模板*/
/************************************************************************/
void _no_model()
{
	img1 = imread("1.jpg", 1);
	namedWindow("原始图");
	imshow("原始图", img1);
	namedWindow("效果图");
	createTrackbar("半径r：", "效果图", &r, 30, on_guilt_feliter);
	createTrackbar("参数e：", "效果图", &e, 1000, on_guilt_feliter_e);
	on_guilt_feliter(r, 0);
	on_guilt_feliter_e(e, 0);

}
/************************************************************************/
/*主函数入口 */
/************************************************************************/
int _tmain(int argc, _TCHAR* argv[])
{
	_no_model();
	cout << "the programe is over!" << endl;
	return 0;
}

引导滤波的算法在下面

#include <opencv.hpp>
#include <opencv2/core/core.hpp>
#include <opencv2/highgui/highgui.hpp>

using namespace std;
using namespace cv;

//声明全局变量
int r = 16;
int e = 1;

//定义显示图像
Mat img1, img2;
Mat o_img;

//调节滑块
static void on_guilt_feliter(int, void*);
static void on_guilt_feliter_e(int, void*);
/************************************************************************/
/*获取图像，定义64深度图，进行乘积运算，防止过界 */
/************************************************************************/
Mat img_Get(Mat &a)
{
	int h = a.rows;
	int w = a.cols;
	Mat I(h, w, CV_64FC1);

	a.convertTo(I, CV_64FC1);

	for (int i = 0; i < h; i++)
	{
		double *p = I.ptr<double>(i);
		for (int j = 0; j < w; j++)
		{
			p[j] = p[j] / 255.0;
		}
	}
	return I;
}
/************************************************************************/
/* 统计sum*/
/************************************************************************/

Mat calculate_Sum(Mat &img_Src, int rc)
{
	int h = img_Src.rows;
	int w = img_Src.cols;
	Mat img_Sum = img_Src.clone();
	if (rc == 1)
	{
		for (int i = 1; i < h; i++)
		{
			for (int j = 0; j < w; j++)
			{
				img_Sum.at<double>(i, j) += img_Sum.at<double>(i - 1, j);
			}
		}
	}

	if (rc == 2)
	{
		for (int i = 0; i < h; i++)
		{
			for (int j = 1; j < w; j++)
			{
				img_Sum.at<double>(i, j) += img_Sum.at<double>(i, j - 1);
			}
		}
	}
	return img_Sum;
}
/************************************************************************/
/*方框滤波函数调用，开源代码 */
/************************************************************************/
Mat boxfilter(Mat &img_Src, int r)
{
	int h = img_Src.rows;
	int w = img_Src.cols;
	Mat img_dst = Mat::zeros(h, w, CV_64FC1);

	Mat img_Sum = calculate_Sum(img_Src, 1);

	for (int i = 0; i < r + 1; i++)
	{
		for (int j = 0; j < w; j++)
		{
			img_dst.at<double>(i, j) = img_Sum.at<double>(i + r, j);
		}
	}

	for (int i = r + 1; i < h - r; i++)
	{
		for (int j = 0; j < w; j++)
		{
			img_dst.at<double>(i, j) = img_Sum.at<double>(i + r, j) - img_Sum.at<double>(i - r - 1, j);
		}
	}

	for (int i = h - r; i < h; i++)
	{
		for (int j = 0; j < w; j++)
		{
			img_dst.at<double>(i, j) = img_Sum.at<double>(h - 1, j) - img_Sum.at<double>(i - r - 1, j);
		}
	}
	img_Sum = calculate_Sum(img_dst, 2);

	for (int i = 0; i < h; i++)
	{
		for (int j = 0; j < r + 1; j++)
		{
			img_dst.at<double>(i, j) = img_Sum.at<double>(i, j + r);
		}
	}

	for (int i = 0; i < h; i++)
	{
		for (int j = r + 1; j < w - r; j++)
		{
			img_dst.at<double>(i, j) = img_Sum.at<double>(i, j + r) - img_Sum.at<double>(i, j - r - 1);
		}
	}

	for (int i = 0; i < h; i++)
	{
		for (int j = w - r; j < w; j++)
		{
			img_dst.at<double>(i, j) = img_Sum.at<double>(i, w - 1) - img_Sum.at<double>(i, j - r - 1);
		}
	}
	return img_dst;
}
/************************************************************************/
/*引导滤波函数u，实现，参考公式，主要计算均值，方差，协方差等 */
/************************************************************************/
Mat guidedfilter(Mat &I, Mat &p, int r, double eps)
{
	int h = I.rows;
	int w = I.cols;

	Mat one = Mat::ones(h, w, CV_64FC1);
	Mat N = boxfilter(one, r);


	Mat mean_I(h, w, CV_64FC1);
	divide(boxfilter(I, r), N, mean_I);


	Mat mean_p(h, w, CV_64FC1);
	divide(boxfilter(p, r), N, mean_p);


	Mat mul_Ip(h, w, CV_64FC1);
	Mat mean_Ip(h, w, CV_64FC1);
	multiply(I, p, mul_Ip);
	divide(boxfilter(mul_Ip, r), N, mean_Ip);


	Mat mul_mean_Ip(h, w, CV_64FC1);
	Mat cov_Ip(h, w, CV_64FC1);
	multiply(mean_I, mean_p, mul_mean_Ip);
	subtract(mean_Ip, mul_mean_Ip, cov_Ip);


	Mat mul_II(h, w, CV_64FC1);
	Mat mean_II(h, w, CV_64FC1);
	multiply(I, I, mul_II);
	divide(boxfilter(mul_II, r), N, mean_II);


	Mat mul_mean_II(h, w, CV_64FC1);
	Mat var_I(h, w, CV_64FC1);
	multiply(mean_I, mean_I, mul_mean_II);
	subtract(mean_II, mul_mean_II, var_I);


	Mat a(h, w, CV_64FC1);
	for (int i = 0; i < h; i++)
	{
		double *p = var_I.ptr<double>(i);
		for (int j = 0; j < w; j++)
		{
			p[j] = p[j] + eps;
		}
	}
	divide(cov_Ip, var_I, a);


	Mat a_mean_I(h, w, CV_64FC1);
	Mat b(h, w, CV_64FC1);
	multiply(a, mean_I, a_mean_I);
	subtract(mean_p, a_mean_I, b);


	Mat mean_a(h, w, CV_64FC1);
	divide(boxfilter(a, r), N, mean_a);
	Mat mean_b(h, w, CV_64FC1);
	divide(boxfilter(b, r), N, mean_b);

	Mat mean_a_I(h, w, CV_64FC1);
	Mat q(h, w, CV_64FC1);
	multiply(mean_a, I, mean_a_I);
	add(mean_a_I, mean_b, q);

	return q;
}