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

		}

	}

}