#include <opencv2\highgui\highgui.hpp>
#include <iostream>
#include<vector>
using namespace  cv;
using namespace  std;

void GetGaussianKernel(double*& gaus_1, const int size, const double sigma_s);
void gaussianFilter2(const vector<uchar>& corrupted, vector<uchar> &smooth, double*& templates,const int width,const int height, const double sigma_r, const int size_m);


int main() {

	const double sigma_r = 30;	//值域的sigma
	const double PI = 4.0*atan(1.0); //圓周率π賦值
	const int size_m = 5;	//模板大小
	const double sigma_s = 10;	//空間域的sigma
	double *templates;

	templates = new double[size_m*size_m];
	GetGaussianKernel(templates, size_m, sigma_s);
	Mat img = imread("123.jpg", 3);

//	namedWindow("MyWindow");
//	imshow("MyWindow", img);

	vector<uchar> array(img.rows*img.cols*3);
	if (img.isContinuous()) { array.assign(img.datastart, img.dataend); }

	vector<uchar> no(img.rows*img.cols*3);

	gaussianFilter2(array, no, templates ,int(img.cols)*3, img.rows,sigma_r,size_m);

	Mat now((int)img.rows, (int)img.cols, CV_8UC3 );

	for (int i = 0; i < img.rows; i++)
		for (int j = 0; j < img.cols; j++) {
			now.at<Vec3b>(i, j)[0] = no[i*img.cols*3 + j*3];
			now.at<Vec3b>(i, j)[1] = no[i*img.cols*3 + j*3 + 1];
			now.at<Vec3b>(i, j)[2] =no[i*img.cols*3 + j*3 + 2];
		}

	imwrite("1123.jpg", now);
	namedWindow("MyWindow1");
	imshow("MyWindow1", now);
	waitKey(0);
	return 0;
}


void GetGaussianKernel(double*& gaus_1, const int size, const double sigma_s)
{
	double **gaus = new double*[size];
	for (int i = 0; i<size; i++)gaus[i] = new double[size];

	double sum = 0;
	for (int i = -size / 2; i<size / 2 + 1; i++) {
		for (int j = -size / 2; j<size / 2 + 1; j++) {
			gaus[i + size / 2][j + size / 2] = exp(-((i*i) + (j*j)) / (2 * sigma_s*sigma_s));
			sum += gaus[i + size / 2][j + size / 2];
		}
	}

	for (int i = 0; i<size; i++) {
		for (int j = 0; j<size; j++) {
			gaus[i][j] /= sum;
			gaus_1[i*size + j] = gaus[i][j];	//使用一維更簡單
		}
	}
	return;
}



void gaussianFilter2(const vector<uchar>& corrupted, vector<uchar> &smooth,double*& templates,const int width,const int height,const double sigma_r, const int size_m) {
	int len = size_m / 2;
	smooth = corrupted;  //復制像素

	for (int j = 0; j<height ; j++) {  //邊緣不處理
		for (int i = 0; i<width ; i++) {
			double sum = 0;
			int index = 0;
			double sum_c = 0;
			double temp = 0;
			for (int m = j - len; m<j + len + 1; m++) {
				for (int n = i - 3 * len; n<i + 3 * len + 1; n += 3) {
					if (m<0 || n<0 || m>height - 1 || n>width - 1)continue;　　//邊緣處理
					temp = templates[index++] * exp(-(corrupted[m*width + n] - corrupted[j*width + i])*(corrupted[m*width + n] - corrupted[j*width + i]) / (2.0*sigma_r*sigma_r));
					sum += corrupted[m*width + n] * temp;
					sum_c += temp;
				}
			}
			sum /= sum_c;
			if (sum > 255)
				sum = 255;
			if (sum < 0)
				sum = 0;
			smooth[j*width + i] = sum;
		}
	}
}

vs2015+opencv3.3.1 實現 c++ 雙邊濾波器