高级图像处理基础实验与图像处理基础理论知识

实验1：图像灰度变换

实验一内容

1、利用OpenCV读取图像
具体内容：用OpenCV打开图像，并在窗口中显示。
2、灰度图像二值化处理
具体内容：设置并调整阈值对图像进行二值化处理。
3、灰度图像的对数变换
具体内容：设置并调整r值对图像进行对数变换。
4、灰度图像的伽马变换
具体内容：设置并调整γ值对图像进行伽马变换。
5、彩色图像的补色变换
具体内容：对彩色图像进行补色变换。

实验一理论知识

灰度值计算公式:
$$Gray=0.299R+0.587G+0.114B$$
对于灰度值，0是黑色，255是白色。
对数变换公式:
$$s=\frac{c\log (1+vr)}{log(v+1)}$$
r的灰度值归一化为[0,1]之间的值。对数变换可以拉伸范围较窄的低灰度值，同时压缩范围较宽的高灰度值。可以用来扩展图像中的暗像素值，同时压缩亮像素值。常用值：c = 1.0, v = 5.0。
伽马变换公式:
$$s=c{r}^{\gamma }$$
当γ>1时将较窄范围的低灰度值映射为较宽范围的灰度值，同时将较宽范围的高灰度值映射为较窄范围的灰度值；当γ<1时，情况相反，与反对数变换类似。
当γ<1时，γ的值越小，对图像低灰度值的扩展越明显；当γ>1时，γ的值越大，对图像高灰度值部分的扩展越明显。这样就能够显示更多的图像的低灰度或者高灰度细节。
反色、补色变换公式:
$$C_{New}=255-C_{Old}$$
$$C_{New}=Max(R,G,B)+Min(R,G,B)-C_{Old}$$
一种颜色反色或补色后，颜色的色相值与反相前相差180度，这是反色和补色的共同点。不同点是，反色除了要改变色相以外，纯度和明度也做出了相应的改变，而补色不改变纯度和明度。例如，黑色的反色是白色，黑色的补色还是黑色。

VS2017中配置opencv3.4.1环境

从这里:https://opencv.org/releases/page/2/ 下载opencv3.4.1，选择windows版，如果下载时速度太慢请打开vpn。
下载完成后得到一个自解压压缩包，解压，将得到的opencv341文件夹放在如C:\Program Files位置，并在环境变量中添加C:\Program Files\opencv341\build\x64\vc15\bin。
新建一个VS的C++ Win32 Console Application工程，点击项目->属性->VC++目录，平台选择x64，配置选择活动（debug）。
包含目录项添加下面的路径:

C:\Program Files\opencv341\build\include
C:\Program Files\opencv341\build\include\opencv
C:\Program Files\opencv341\build\include\opencv2

库目录项添加下面的路径:

C:\Program Files\opencv341\build\x64\vc15\lib

点击链接器->输入->附加依赖项，添加下面的路径:

opencv_world341d.lib //d表示debug版本

这样我们的opencv环境就配置好了。我们还可以把这个配置保存下来，在菜单栏中选择：视图-> 其他窗口->属性管理器，找到其中的Debug|x64项，点击添加新属性表，命名为opencv341.props，然后将这个属性表中属性按上面修改一遍，保存下来，这样在新建其他项目时只需要在这个属性管理器中导入opencv341.props属性表即可。

实验一代码

#include "pch.h"
#include <string>
#include <iostream>
#include<opencv2/opencv.hpp>
#include <opencv2/core/core.hpp> 
#include <opencv2/highgui/highgui.hpp> //包含imread, imshow等
#include <opencv2/imgproc/imgproc.hpp> //包含cvtColor等

using namespace std;
using namespace cv;

// 读取、显示、保存图片
void loadAndSaveImage(string load_path, string save_path) {
	Mat pScr;
	// 1为加载图像的颜色类型,1为原图返回,0为灰度返回
	pScr = imread(load_path, 1);
	imshow("picture", pScr);
	imwrite(save_path, pScr);
	// 等待按下某键(ASCII码值),若为0则按下任意键后程序退出,若写3000就是3秒后程序自动退出
	waitKey(0);
	destroyAllWindows();
}

// 图像二值化处理,使用整幅图像的灰度平均值作为二值化的阈值
void binaryImage(string load_path, string save_path_gray, string save_path_gray_binary) {
	Mat pScr;
	// 1为加载图像的颜色类型,1为原图返回,0为灰度返回
	pScr = imread(load_path, 1);
	// 图像行列数
	int row = pScr.rows;
	int col = pScr.cols;
	Mat graypScr(row, col, CV_8UC1);
	for (int i = 0; i < row; i++) {
		for (int j = 0; j < col; j++)
			// 0、1、2位置分别为b、g、r通道
			graypScr.at<uchar>(i, j) = int(0.2989*int(pScr.at<Vec3b>(i, j)[0]) + 0.5870*int(pScr.at<Vec3b>(i, j)[1]) + 0.1140*int(pScr.at<Vec3b>(i, j)[2]));
	}
	// imshow("gray_picture", graypScr);
	// waitKey(0);
	// destroyAllWindows();
	imwrite(save_path_gray, graypScr);
	int sum = 0, threshold = 0;
	// 计算阈值为所有灰度值的平均值
	for (int i = 0; i < row; i++) {
		for (int j = 0; j < col; j++)
			sum = sum + int(graypScr.at<uchar>(i, j));
	}
	threshold = int(sum / (row*col));
	// 灰度图像二值化,大于等于阈值置为255,否则置为0
	for (int i = 0; i < row; i++) {
		for (int j = 0; j < col; j++) {
			if (int(graypScr.at<uchar>(i, j)) >= threshold)
				graypScr.at<uchar>(i, j) = 255;
			else
				graypScr.at<uchar>(i, j) = 0;
		}
	}
	// imshow("binary_gray_picture", graypScr);
	// waitKey(0);
	// destroyAllWindows();
	imwrite(save_path_gray_binary, graypScr);;
}

// 灰度图像的对数变换s = c * log(1+r*v)/log(v),更突出暗部细节
// 对数变换将源图像中范围较窄的低灰度值映射到范围较宽的灰度区间,同时将范围较宽的高灰度值区间映射为较窄的灰度区间
// 从而扩展了暗像素的值，压缩了高灰度的值，能够对图像中低灰度区域进行增强
void LogImage(string load_path, string save_path_gray_log) {
	Mat graypScr;
	// 1为加载图像的颜色类型,1为原图返回,0为灰度返回
	graypScr = imread(load_path, 0);
	// imshow("picture", graypScr);
	// waitKey(0);
	// destroyAllWindows();
	// 图像行列数
	int row = graypScr.rows;
	int col = graypScr.cols;
	double c = 1.0, v = 5.0;
	// 灰度图像的对数变换
	for (int i = 0; i < row; i++) {
		for (int j = 0; j < col; j++) {
			double r = graypScr.at<uchar>(i, j) / 255.0;
			double s = c * log(1 + r * v) / log(v + 1);
			int d = int(s * 255);
			if (d < 0)
				d = 0;
			else if (d > 255)
				d = 255;
			graypScr.at<uchar>(i, j) = d;
		}
	}
	// imshow("gray_log_picture", graypScr);
	// waitKey(0);
	// destroyAllWindows();
	imwrite(save_path_gray_log, graypScr);
}

// 灰度图像的伽马变换s = c * pow(r, gama)
// (255,255,255)是白色，灰度值最高
// γ > 1时，会将低于某个灰度值K的灰度区域压缩到较小的灰度区间,而将高于K的灰度区域扩展到较大灰度区间。令L为灰度的最大值,k = 3 / 4L
// [0, 3 / 4L]的灰度区域映射到为[0, 1 / 8L]的输出,而将[3 / 4L, L]这部分高灰度区域映射到[1 / 8L, L]区间
// 这样变换的结果就是，低于k的灰度区域被压缩到更低灰度区间,而较亮的高灰度区域的灰度值被扩展到较大的灰度区间变的不那么亮,整体的效果就是图像的对比度增加,但是由于亮度区域被扩展,也就不那么亮了。
// γ < 1时,会将灰度值较小的低灰度区域扩展到较宽的灰度区间,而将较宽的高灰度区域压缩到较小的灰度区间。这样低灰度区域扩展开来变亮,而宽的高灰度区域被压缩到较窄的区间,也变亮,故变换后的整体效果是变亮了。
void gamaImage(string load_path, string save_path_gray_gamma) {
	Mat graypScr;
	// 1为加载图像的颜色类型,1为原图返回,0为灰度返回
	graypScr = imread(load_path, 0);
	// 图像行列数
	int row = graypScr.rows;
	int col = graypScr.cols;
	double c = 1.0, gama = 3.0;
	for (int i = 0; i < row; i++) {
		for (int j = 0; j < col; j++) {
			double r = graypScr.at<uchar>(i, j) / 255.0;
			double s = c * pow(r, gama);
			int d = int(s * 255);
			if (d < 0)
				d = 0;
			else if (d > 255)
				d = 255;
			graypScr.at<uchar>(i, j) = d;
		}
	}
	// imshow("gray_gamma_picture", graypScr);
	// waitKey(0);
	// destroyAllWindows();
	imwrite(save_path_gray_gamma, graypScr);
}

// 彩色图像的补色变换(反色差)
void imcomplementImage(string load_path, string save_path_color) {
	Mat pScr;
	// 1为加载图像的颜色类型,1为原图返回,0为灰度返回
	pScr = imread(load_path, 1);
	// 图像行列数
	int row = pScr.rows;
	int col = pScr.cols;
	for (int i = 0; i < row; i++) {
		for (int j = 0; j < col; j++) {
			pScr.at<Vec3b>(i, j)[0] = 255 - pScr.at<Vec3b>(i, j)[0];
			pScr.at<Vec3b>(i, j)[1] = 255 - pScr.at<Vec3b>(i, j)[1];
			pScr.at<Vec3b>(i, j)[2] = 255 - pScr.at<Vec3b>(i, j)[2];
		}
	}
	// imshow("color_picture", graypScr);
	// waitKey(0);
	// destroyAllWindows();
	imwrite(save_path_color, pScr);
}


int main(int argc, char * argv[]) {
	string load_path = "C:/Users/zgcr6/Desktop/高图实验/1/save/lenna.bmp";
	string save_path = "C:/Users/zgcr6/Desktop/高图实验/1/save/lenna_save.bmp";
	string save_path_gray = "C:/Users/zgcr6/Desktop/高图实验/1/save/lenna_gray.bmp";
	string save_path_gray_binary = "C:/Users/zgcr6/Desktop/高图实验/1/save/lenna_binary_gray.bmp";
	string save_path_gray_log = "C:/Users/zgcr6/Desktop/高图实验/1/save/lenna_gray_log.bmp";
	string save_path_gray_gamma = "C:/Users/zgcr6/Desktop/高图实验/1/save/lenna_gray_gamma.bmp";
	string save_path_color = "C:/Users/zgcr6/Desktop/高图实验/1/save/lenna_color.bmp";
	loadAndSaveImage(load_path, save_path);
	binaryImage(load_path, save_path_gray, save_path_gray_binary);
	LogImage(load_path, save_path_gray_log);
	gamaImage(load_path, save_path_gray_gamma);
	imcomplementImage(load_path, save_path_color);
	return 0;
}

实验二：直方图均衡

实验二内容

1、计算灰度图像的归一化直方图
具体内容：利用OpenCV对图像像素进行操作，计算归一化直方图，并在窗口中以图形的方式显示出来。
2、灰度图像直方图均衡处理
具体内容：通过计算归一化直方图，设计算法实现直方图均衡化处理。
3、彩色图像直方图均衡处理
具体内容： 在灰度图像直方图均衡处理的基础上实现彩色直方图均衡处理。

实验二理论知识

归一化直方图公式:
先计算每个灰度值的像素个数
$$\mathrm{h}\left(r_k\right)=n_k$$
然后都除以全图像素个数，得到灰度值分布概率函数P（rk），我们就可以得到归一化的灰度直方图。
$$\mathrm{P}\left(r_k\right)=\frac{n_k}{n}$$
直方图均衡化公式:
用上面得到的灰度值分布概率函数P（rk）计算灰度值累积分布函数Sk:
$$S_k=\sum _{i=0}^{k-1}P\left(r_k\right)=\sum _{i=0}^{k-1}\frac{n_k}{n}$$
映射后的灰度值为:
$$G_k=(L-1)S_k$$
注意L-1=255，累积分布函数Sk和输入图像的灰度级相乘，不一定得到一个整数，但是输出的灰度级要求整数，于是我们这里四舍五入取整数灰度级。最后我们只要遍历所有像素点将原灰度值换成映射后的灰度值即可。对于彩色图像的直方图均衡，在R、G、B三个通道上分别做直方图均衡即可。

实验二代码

#include "pch.h"
#include <string>
#include <iostream>
#include <opencv2/opencv.hpp>
#include <opencv2/core/core.hpp> 
#include <opencv2/highgui/highgui.hpp> //包含imread, imshow等
#include <opencv2/imgproc/imgproc.hpp> //包含cvtColor等

using namespace std;
using namespace cv;

void grayImageGetHistImage(string load_path, string save_path_hist_image) {
	Mat graypScr;
	// 1为加载图像的颜色类型,1为原图返回,0为灰度返回
	graypScr = imread(load_path, 0);
	// imshow("gray_picture", graypScr);
	// waitKey(0);
	// destroyAllWindows();
	int row = graypScr.rows;
	int col = graypScr.cols;
	// 用一个数组记录0-255的灰度值的数量,数组下标就是灰度值
	int grayNum[256] = { 0 };
	// 遍历所有像素点,记录每一个灰度值的数量
	for (int i = 0; i < row; i++) {
		for (int j = 0; j < col; j++) {
			grayNum[int(graypScr.at<uchar>(i, j))] = grayNum[int(graypScr.at<uchar>(i, j))] + 1;
		}
	}
	// 找出灰度值数量的最大值
	double maxValue = 0;
	for (int i = 0; i < 256; i++) {
		if (grayNum[i] > maxValue)
			maxValue = grayNum[i];
	}
	// 划分成256块
	int bin_num = 256;
	// 生成白色画布(255)
	Mat histImg(bin_num, bin_num, CV_8U, Scalar(255));
	// 每个块是一条垂线	
	int hpt = static_cast<int>(0.9*bin_num);
	// 对每一个块,用其灰度值数量乘以hpt再除以灰度值数量最大值，这样就转化成程度
	for (int i = 0; i < bin_num; i++) {
		int intensity = static_cast<int>(grayNum[i] * hpt / maxValue);
		//直方图底部点与直方图顶部点之间绘制一条线		
		line(histImg, Point(i, bin_num), Point(i, bin_num - intensity), Scalar::all(0));
	}
	// imshow("hist_image", histImg);
	// waitKey(0);
	// destroyAllWindows();
	imwrite(save_path_hist_image, histImg);
}

// 灰度直方图均衡化步骤:
// 统计原始图像各灰度值的数量;图像中灰度值为i的像素的出现概率是:px(i)=p(x=i)=ni/n,n是图像中所有的像素数,px(i)实际上是灰度值为i的图像的直方图,归一化到[0, 1];
// px的累积分布函数,是图像的累计归一化直方图:cdfx(i)=j从0到i求和px(j)
// 直方图均衡化计算公式,cdfmin为累积分布函数最小值,M和N分别代表了图像的长宽像素个数,L是灰度值的个数(如256),v为原始图像中为v的像素值
// h(v)=round((cdf(v)-cdfmin)*(L-1)/(M*N-cdfmin))
void grayImageHistEqualization(string load_path, string save_path_gray, string save_path_gray_image_equalization, string save_path_hist_image_equalization) {
	Mat graypScr;
	// 1为加载图像的颜色类型,1为原图返回,0为灰度返回
	graypScr = imread(load_path, 0);
	// imshow("gray_picture", graypScr);
	// waitKey(0);
	// destroyAllWindows();
	imwrite(save_path_gray, graypScr);
	int row = graypScr.rows;
	int col = graypScr.cols;
	// 用一个数组记录0-255的灰度值的数量,数组下标就是灰度值
	double grayProbability[256] = { 0.0 };
	// 遍历所有像素点,记录每一个灰度值的数量
	for (int i = 0; i < row; i++) {
		for (int j = 0; j < col; j++) {
			grayProbability[int(graypScr.at<uchar>(i, j))] = grayProbability[int(graypScr.at<uchar>(i, j))] + 1.0;
		}
	}
	// 求图像中灰度值为i的像素出现概率
	for (int i = 0; i < 256; i++)
		grayProbability[i] = grayProbability[i] / (row*col*1.0);
	double cdf[256] = { 0.0 };
	// 求px的累积分布函数,每个i的累积分布函数是灰度值从0到i的像素出现概率的累加之和
	for (int i = 0; i < 256; i++) {
		for (int j = 0; j <= i; j++)
			cdf[i] += grayProbability[j];
	}
	int L = 255;
	// 该数组记录了均衡后原始灰度值(数组下标)对应的均衡后的灰度值,计算式为D[i]=L*cdf[i]
	int D[256] = { 0 };
	for (int i = 0; i < 256; i++) {
		D[i] = int(L*cdf[i]);
	}
	// 对灰度图像做均衡化处理,每个像素点原有灰度值替换为均衡后的灰度值
	for (int i = 0; i < row; i++) {
		for (int j = 0; j < col; j++) {
			int grayValue = int(graypScr.at<uchar>(i, j));
			graypScr.at<uchar>(i, j) = D[grayValue];
		}
	}
	// imshow("gray_equalization_picture", graypScr);
	// waitKey(0);
	// destroyAllWindows();
	imwrite(save_path_gray_image_equalization, graypScr);
	grayImageGetHistImage(save_path_gray_image_equalization, save_path_hist_image_equalization);
}

void colorImageHistEqualization(string load_path, string save_path_color_image_equalization) {
	Mat pScr;
	// 1为加载图像的颜色类型,1为原图返回,0为灰度返回
	pScr = imread(load_path, 1);
	// imshow("gray_picture", graypScr);
	// waitKey(0);
	// destroyAllWindows();
	int row = pScr.rows;
	int col = pScr.cols;
	for (int k = 0; k < 3; k++) {
		Mat singleChannel(row, col, CV_8UC1);
		for (int i = 0; i < row; i++) {
			for (int j = 0; j < col; j++) {
				// 0、1、2位置分别为b、g、r通道,每次提取一个通道的值存在一个单通道的Mat矩阵中,作为灰度图像处理
				singleChannel.at<uchar>(i, j) = int(pScr.at<Vec3b>(i, j)[k]);
			}
		}
		// 用一个数组记录0-255的灰度值的数量,数组下标就是灰度值
		double grayProbability[256] = { 0.0 };
		// 遍历所有像素点,记录每一个灰度值的数量
		for (int i = 0; i < row; i++) {
			for (int j = 0; j < col; j++) {
				grayProbability[int(singleChannel.at<uchar>(i, j))] = grayProbability[int(singleChannel.at<uchar>(i, j))] + 1.0;
			}
		}
		// 求图像中灰度值为i的像素出现概率
		for (int i = 0; i < 256; i++)
			grayProbability[i] = grayProbability[i] / (row*col*1.0);
		double cdf[256] = { 0.0 };
		// 求px的累积分布函数,每个i的累积分布函数是灰度值从0到i的像素出现概率的累加之和
		for (int i = 0; i < 256; i++) {
			for (int j = 0; j <= i; j++)
				cdf[i] += grayProbability[j];
		}
		int L = 255;
		// 该数组记录了均衡后原始灰度值(数组下标)对应的均衡后的灰度值,计算式为D[i]=L*cdf[i]
		int D[256] = { 0 };
		for (int i = 0; i < 256; i++) {
			D[i] = int(L*cdf[i]);
		}
		// 对灰度图像做均衡化处理,每个像素点原有灰度值替换为均衡后的灰度值
		for (int i = 0; i < row; i++) {
			for (int j = 0; j < col; j++) {
				int grayValue = int(singleChannel.at<uchar>(i, j));
				singleChannel.at<uchar>(i, j) = D[grayValue];
			}
		}
		// imshow("single_channel_equalization_picture", singleChannel);
		// waitKey(0);
		// destroyAllWindows();
		for (int i = 0; i < row; i++) {
			for (int j = 0; j < col; j++) {
				// 将均衡化后的灰度值分别写回b、g、r通道
				pScr.at<Vec3b>(i, j)[k] = int(singleChannel.at<uchar>(i, j));
			}
		}
	}
	// imshow("color_image_equalization_picture", pScr);
	// waitKey(0);
	// destroyAllWindows();
	imwrite(save_path_color_image_equalization, pScr);
}

int main(int argc, char * argv[]) {
	string load_path = "C:/Users/zgcr6/Desktop/高图实验/2/save/lenna.bmp";
	string save_path_hist_image = "C:/Users/zgcr6/Desktop/高图实验/2/save/lenna_hist_image.bmp";
	string save_path_gray = "C:/Users/zgcr6/Desktop/高图实验/2/save/lenna_gray.bmp";
	string save_path_gray_image_equalization = "C:/Users/zgcr6/Desktop/高图实验/2/save/lenna_gray_image_equalization.bmp";
	string save_path_hist_image_equalization = "C:/Users/zgcr6/Desktop/高图实验/2/save/lenna_hist_image_equalization.bmp";
	string save_path_color_image_equalization = "C:/Users/zgcr6/Desktop/高图实验/2/save/lenna_color_image_equalization.bmp";
	grayImageGetHistImage(load_path, save_path_hist_image);
	grayImageHistEqualization(load_path, save_path_gray, save_path_gray_image_equalization, save_path_hist_image_equalization);
	colorImageHistEqualization(load_path, save_path_color_image_equalization);
	return 0;
}

实验三：空域滤波

实验三内容

1、利用均值模板平滑灰度图像
具体内容：利用OpenCV对图像像素进行操作，分别利用3x3、5x5和9x9尺寸的均值模板平滑灰度图像。
2、利用高斯模板平滑灰度图像
具体内容：利用OpenCV对图像像素进行操作，分别利用3x3、5x5和9x9尺寸的高斯模板平滑灰度图像。
3、利用Laplacian、Robert、Sobel模板锐化灰度图像
具体内容：利用OpenCV对图像像素进行操作，分别利用Laplacian、Robert、Sobel模板锐化灰度图像。
4、利用高提升滤波算法增强灰度图像
具体内容：利用OpenCV对图像像素进行操作，设计高提升滤波算法增强图像。
5、利用均值模板平滑彩色图像
具体内容：利用OpenCV分别对图像像素的RGB三个通道进行操作，利用3x3、5x5和9x9尺寸的均值模板平滑彩色图像。
6、利用高斯模板平滑彩色图像
具体内容：利用OpenCV分别对图像像素的RGB三个通道进行操作，分别利用3x3、5x5和9x9尺寸的高斯模板平滑彩色图像。
7、利用Laplacian、Robert、Sobel模板锐化灰度图像
具体内容：利用OpenCV分别对图像像素的RGB三个通道进行操作，分别利用Laplacian、Robert、Sobel模板锐化彩色图像。

实验三理论知识

算数平均值滤波公式:
$$f(x,y)=\frac{1}{mn}\sum _{(x,y)\in S_{xy}}g(s,t)$$
Sxy表示中心点在（x,y）处，大小为m×n的滤波器窗口。g（s,t）表示原始图像灰度值，f（x,y）表示均值滤波后得到的图像灰度值。均值滤波就是把窗口内的像素点上灰度值全部按权重1相加求和，再除以像素点个数，得到的灰度值就是窗口正中心像素点的灰度值。对于彩色图像，对其三个通道分别执行均值滤波即可。
高斯滤波公式:
先计算mxn大小的窗口中每个像素点的权重:
$$h(x,y)=e^{-\frac{(x-x_0{)}^2+(y-y_0{)}^2}{2{\sigma }^2}}$$
其中（x0，y0）是窗口正中心的像素点坐标。
然后我们要将所有权重归一化，即上面的h（x，y）都乘以下面的系数:
$$\frac{1}{{\sum }_{(i,j)\in w}{w}_{i,j}}$$
最后将窗口中每个像素点灰度值乘以对应的权重并求和，得到的灰度值就是窗口正中心像素点的灰度值。对于彩色图像，对其三个通道分别执行高斯滤波即可。
Laplacian算子:
差分就是离散函数连续两项之差。
我们用图像灰度值梯度的变化来检测图像中物体的边缘，边缘处通常灰度值变化会很剧烈，相应的灰度值梯度会较大。但是在图像上各个像素点是离散点，因此我们用对x,y两个方向的二阶导数的差分来近似表示梯度。这就是Laplacian算子。
一维一阶差分公式:
$$\frac{\partial f}{\partial x}=f(x+1)-f(x)$$
一维二阶差分公式:
$$\frac{{\partial }^{2}f}{\partial x^2}=f(x+1)+f(x-1)-2f(x)$$
在一个二维函数f（x,y）中，x,y两个方向的二阶差分分别为:
$$\frac{{\partial }^{2}f}{\partial x^2}=f(x+1,y)+f(x-1,y)-2f(x,y)$$
$$\frac{{\partial }^{2}f}{\partial y^2}=f(x,y+1)+f(x,y-1)-2f(x,y)$$
Laplace算子就是将上面两个二阶差分的结果相加:
$${\nabla }^{2}f(x,y)=f(x+1,y)+f(x-1,y)+f(x,y+1)+f(x,y-1)-4f(x,y)$$
写成模板形式为:

0    1    0
1    -4   1
0    1    0

当然还有扩展的模板形式:

1    1    1        0    -1    0        -1    -1    -1
1    -8   1       -1     4   -1        -1    8     -1        
1    1    1        0    -1    0        -1    -1    -1

对于彩色图像，对其三个通道分别执行Laplacian算子即可。
Robert算子:
使用一个对角线之差的绝对值之和来近似代表梯度。
$$M(x,y)\approx \left|z_9-z_5\right|+\left|z_8-z_6\right|$$
模板形式:

z1    z2    z3
z4    z5    z6        -1    0        0    -1
z7    z8    z9         0    1        1    0

对于彩色图像，对其三个通道分别执行Robert算子即可。
Sobel算子:
对Robert算子的改进。
$$M(x,y)\approx \left|\left(z_7+2z_8+z_9\right)-\left(z_1+2z_2+z_3\right)\right|+\left|\left(z_3+2z_6+z_9\right)-\left(z_1+2z_4+z_7\right)\right|$$
模板形式:

-1    0    1        -1    -2    -1
-2    0    2        0     0     0
-1    0    1        1     2     1

对于彩色图像，对其三个通道分别执行Sobel算子即可。
Laplacian算子、Robert算子、Sobel算子都对噪声敏感，因此一般提取边缘前都需要先平滑。
高提升滤波公式:
$$g=f+kmk\{\begin{array}{l}k>1高提升滤波\\ \\ k=1非锐化掩蔽\\ \\ 0<k<1不强调非锐化模板的贡献\end{array}$$
先平滑原图像f，得到s；从原图像f中减去模糊图像s，产生的差值图像称为m=f-s。系数k越大对细节增强越明显；平滑时减弱的边缘，锐化后增强的更加明显。

实验三代码

#include "pch.h"
#include <iostream>
#include <string>
#include <math.h>
#include <opencv2/opencv.hpp>
#include <opencv2/core/core.hpp> 
#include <opencv2/highgui/highgui.hpp> //包含imread, imshow等
#include <opencv2/imgproc/imgproc.hpp> //包含cvtColor等

using namespace std;
using namespace cv;

void gray_image(string load_path, string save_path_gray_image) {
	Mat graypScr;
	// 1为加载图像的颜色类型,1为原图返回,0为灰度返回
	graypScr = imread(load_path, 0);
	//imshow("gray_picture", graypScr);
	//waitKey(0);
	//destroyAllWindows();
	imwrite(save_path_gray_image, graypScr);
}

// 均值滤波,既可以对灰度图像,也可以对彩色图像
void mean_filter(string load_path, int mean_filter_size, string save_path_mean_filter) {
	// 滤波窗口尺寸必须为奇数
	Mat pScr;
	// 1为加载图像的颜色类型,1为原图返回,0为灰度返回
	pScr = imread(load_path, 1);
	int row = pScr.rows;
	int col = pScr.cols;
	int channel = pScr.channels();
	//imshow("picture", pScr);
	//waitKey(0);
	//destroyAllWindows();
	// 滤波窗口中心
	int pos = (mean_filter_size - 1) / 2;
	// 从图片矩阵左上角开始从左到右从上到下使用滤波窗口处理图片矩阵,每次提取filter_size*filter_size的数据进入滤波窗口,计算平均值,然后替代滤波窗口中心值
	// 图片按每个通道来处理
	for (int k = 0; k < channel; k++) {
		// 每个滤波窗口中心位置都计算滤波窗口内通道值的平均值
		for (int m = pos; m < row - pos; m++) {
			for (int n = pos; n < col - pos; n++) {
				int centerpointValue = 0;
				for (int i = m - (mean_filter_size - 1) / 2; i <= m + (mean_filter_size - 1) / 2; i++)
					for (int j = n - (mean_filter_size - 1) / 2; j <= n + (mean_filter_size - 1) / 2; j++)
						centerpointValue += int(pScr.at<Vec3b>(i, j)[k]);
				centerpointValue = int(centerpointValue / (mean_filter_size*mean_filter_size));
				pScr.at<Vec3b>(m, n)[k] = centerpointValue;
			}
		}
	}
	// 对边界进行处理,边界上的像素值取最邻近行的滤波窗口中心位置像素值
	for (int k = 0; k < channel; k++)
		for (int i = 0; i < pos; i++)
			for (int j = pos; j < col - pos; j++)
				pScr.at<Vec3b>(i, j)[k] = pScr.at<Vec3b>(pos, j)[k];
	for (int k = 0; k < channel; k++)
		for (int i = row - pos; i < row; i++)
			for (int j = pos; j < col - pos; j++)
				pScr.at<Vec3b>(i, j)[k] = pScr.at<Vec3b>(row - pos - 1, j)[k];
	for (int k = 0; k < channel; k++)
		for (int i = 0; i < row; i++)
			for (int j = 0; j < pos; j++)
				pScr.at<Vec3b>(i, j)[k] = pScr.at<Vec3b>(i, pos)[k];
	for (int k = 0; k < channel; k++)
		for (int i = 0; i < row; i++)
			for (int j = col - pos; j < col; j++)
				pScr.at<Vec3b>(i, j)[k] = pScr.at<Vec3b>(i, col - pos - 1)[k];
	imwrite(save_path_mean_filter, pScr);
}

// 二维高斯滤波,既可以对灰度图像,也可以对彩色图像
void gaussian_filter(string load_path, int gaussian_filter_size, double sigma, string save_path_gaussian_filter) {
	double PI = 3.1415926;
	// 滤波窗口尺寸必须为奇数
	Mat pScr;
	// 1为加载图像的颜色类型,1为原图返回,0为灰度返回
	pScr = imread(load_path, 1);
	int row = pScr.rows;
	int col = pScr.cols;
	int channel = pScr.channels();
	//imshow("picture", pScr);
	//waitKey(0);
	//destroyAllWindows();
	// 滤波窗口中心
	int pos = (gaussian_filter_size - 1) / 2;
	// 记录滤波窗口所有像素点的权值
	vector<double> weight(gaussian_filter_size*gaussian_filter_size, 0.0);
	double weight_sum = 0;
	for (int i = 0; i < gaussian_filter_size; i++) {
		for (int j = 0; j < gaussian_filter_size; j++) {
			weight[i*gaussian_filter_size + j] = exp(-(pow(i - pos, 2) + pow(j - pos, 2)) / (2 * sigma * sigma)) / (2 * PI * sigma*sigma);
			weight_sum += weight[i*gaussian_filter_size + j];
		}
	}
	// 权值归一化
	for (int index = 0; index < gaussian_filter_size*gaussian_filter_size; index++)
		weight[index] /= weight_sum;
	// 这样我们就得到了归一化后的高斯模板权值
	// 从图片矩阵左上角开始从左到右从上到下使用滤波窗口处理图片矩阵,每次提取filter_size*filter_size的数据进入滤波窗口,利用二维高斯函数计算权值,然后将滤波窗口权值归一化,用权值计算出用来替代滤波窗口中心的像素值
	// 图片按每个通道来处理
	for (int k = 0; k < channel; k++) {
		// 每个滤波窗口中心位置都计算滤波窗口内通道值的高斯加权平均值
		for (int m = pos; m < row - pos; m++) {
			for (int n = pos; n < col - pos; n++) {
				// 使用归一化后的权重计算最终加权得到的像素值,用centerpointValue用来记录
				double centerpointValue = 0;
				// 将高斯滤波窗口中每个点的加权值和其像素值相乘再相加,得到的和就是滤波窗口中心点的像素值
				for (int i = m - (gaussian_filter_size - 1) / 2; i <= m + (gaussian_filter_size - 1) / 2; i++) {
					for (int j = n - (gaussian_filter_size - 1) / 2; j <= n + (gaussian_filter_size - 1) / 2; j++) {
						int index = (i - (m - (gaussian_filter_size - 1) / 2))*gaussian_filter_size + (j - (n - (gaussian_filter_size - 1) / 2));
						centerpointValue += weight[index] * int(pScr.at<Vec3b>(i, j)[k]);
					}
				}
				pScr.at<Vec3b>(m, n)[k] = int(centerpointValue);
			}
		}
	}
	imwrite(save_path_gaussian_filter, pScr);
}

// 先对x方向进行高斯滤波,再对y方向进行高斯滤波,既可以对灰度图像,也可以对彩色图像
void separate_gaussian_filter(string load_path, int gaussian_filter_size, double sigma, string save_path_separate_gaussian_filter) {
	const double PI = 3.1415926;
	// 滤波窗口尺寸必须为奇数
	Mat pScr;
	// 1为加载图像的颜色类型,1为原图返回,0为灰度返回
	pScr = imread(load_path, 1);
	int row = pScr.rows;
	int col = pScr.cols;
	int channel = pScr.channels();
	//imshow("picture", pScr);
	//waitKey(0);
	//destroyAllWindows();
	// 滤波窗口中心
	int pos = (gaussian_filter_size - 1) / 2;
	// 计算高斯滤波窗口中每个位置的权值
	vector<double> weight(gaussian_filter_size*gaussian_filter_size, 0.0);
	double weight_sum = 0;
	// 每个滤波窗口记录权重，并记录权值和用来归一化
	for (int i = 0; i < gaussian_filter_size; i++) {
		weight[i] = exp(-pow(i - pos, 2) / (2 * sigma * sigma)) / (sigma*sqrt(2 * PI));
		weight_sum += weight[i];
	}
	// 权值归一化
	for (int i = 0; i < gaussian_filter_size; i++)
		weight[i] /= weight_sum;
	// 先对水平方向进行高斯滤波
	// 图片按每个通道来处理
	for (int k = 0; k < channel; k++) {
		// 对每行上进行高斯滤波
		for (int m = pos; m < col - pos; m++) {
			for (int n = 0; n < row; n++) {
				double centerpointValue = 0.0;
				// 计算每个中心点的像素值
				for (int i = m - (gaussian_filter_size - 1) / 2; i <= m + (gaussian_filter_size - 1) / 2; i++) {
					centerpointValue += weight[i - (m - (gaussian_filter_size - 1) / 2)] * int(pScr.at<Vec3b>(n, i)[k]);
				}
				pScr.at<Vec3b>(n, m)[k] = int(centerpointValue);
			}
		}
	}
	// 再对竖直方向进行高斯滤波
	for (int k = 0; k < channel; k++) {
		// 对每列上进行高斯滤波
		for (int n = 0; n < col; n++) {
			for (int m = pos; m < row - pos; m++) {
				double centerpointValue = 0.0;
				// 计算每个中心点的像素值
				for (int i = m - (gaussian_filter_size - 1) / 2; i <= m + (gaussian_filter_size - 1) / 2; i++) {
					centerpointValue += weight[i - (m - (gaussian_filter_size - 1) / 2)] * int(pScr.at<Vec3b>(i, n)[k]);
				}
				pScr.at<Vec3b>(m, n)[k] = int(centerpointValue);
			}
		}
	}
	imwrite(save_path_separate_gaussian_filter, pScr);
}

// 一阶差分就是离散函数中连续两项之差,当自变量从x变到x+1时,函数y(x)在点x的一阶差分∆y(x)=y(x+1)-y(x)
// 二阶差分就是一阶差分的连续两项之差,当自变量从x变到x+1时，函数y=y(x)一阶差分的差分为∆(∆y(x))=∆y(x+1)-∆y(x)=(y(x+2)-y(x+1))-(y(x+1)-y(x))=y(x+2)-2y(x+1)+y(x)
// 分别对Laplace算子x,y两个方向的二阶导数进行差分就得到了离散函数的拉普拉斯算子(Laplacian):f(x+1,y)+f(x-1,y)+f(x,y+1)+f(x,y-1)-4f(x,y)
// 这是最基本的拉普拉斯模板a,即像素点的上下左右系数为1,再减去4倍的本像素点像素值;扩展模板b是像素点周围8个点系数为1,再减去8倍的本像素点像素值
// 还有模板c像素点的上下左右系数为-1,再加上4倍的本像素点像素值;模板d是像素点周围8个点系数为-1,再加上8倍的本像素点像素值
// 模板得到的是二阶差分结果,再用原图像像素值根据掩膜中心系数为1还是-1加上一个掩膜中心系数乘以二阶差分的结果
// 拉普拉斯算子可增强图像中灰度突变的区域，减弱灰度的缓慢变化区域
// 选择拉普拉斯算子对原图像进行处理,可产生描述灰度突变的图像,再将拉普拉斯图像与原始图像叠加而产生锐化图像
// 注意拉普拉斯结果的像素值范围可能落在了[0,255]之外,而计算机在显示的时候将赋值全部置为0,大于255的像素全部显示成255
// 因为图像中的边缘就是那些灰度发生跳变的区域，所以拉普拉斯锐化模板在边缘检测中很有用
// 简单来说,拉普拉斯算子可以使黑区域的亮点更亮,白区域的暗点更暗
// 同梯度算子一样，拉普拉斯算子也会增强图像中的噪声，有时用拉普拉斯算子进行边缘检测时，可将图像先进行平滑处理
// 拉普拉斯算子,既可以对灰度图像,也可以对彩色图像
void laplacian(string load_path, string save_path_laplacian) {
	Mat pScr, res;
	// 1为加载图像的颜色类型,1为原图返回,0为灰度返回
	pScr = imread(load_path, 1);
	int row = pScr.rows;
	int col = pScr.cols;
	int channel = pScr.channels();
	//imshow("picture", pScr);
	//waitKey(0);
	//destroyAllWindows();
	// 锐化后的图像像素值g(x,y)=f(x,y)+或-二阶差分值(即filter2D)计算后的结果,这里是-
	Mat kernel = (Mat_<float>(3, 3) << 0, 1, 0, 1, -4, 1, 0, 1, 0);
	// 使用上面的核中的权重对图像中每个像素点的像素值进行卷积计算时与下面拉普拉斯函数效果相同
	filter2D(pScr, res, pScr.depth(), kernel);
	// 使用拉普拉斯函数进行锐化,ksize=1时拉普拉斯算子是0,1,0,1,-4,1,0,1,0
	//Laplacian(pScr, res, CV_16S, 1);
	// 将CV_16S型图像转为CV_8U型图像
	//convertScaleAbs(res, res);
	imwrite(save_path_laplacian, res);
}

// 任意一对互相垂直方向上(即斜对角)的差分可以看成求梯度的近似方法,Robert算子是利用这种原理,采用对角方向相邻两像素值之差代替该梯度值
// 假设当前点(x,y)在x方向的梯度为gx,在y方向上的梯度为gy
// 那么当前点的总梯度长度为根号下(gx的平方+gy的平方)
// 但是上面的式子计算比较复杂,我们将总梯度长度简化为|gx|+|gy|
// 我们用gx=f(x,y)-f(x+1,y+1),gy=f(x+1,y)-f(x,y+1)来近似代表x方向和y方向上的梯度
// 所以Robert模板的权重为0,0,0,0,-1,0,0,0,1和0,0,0,0,0,-1,0,1,0
// robert算子,既可以对灰度图像,也可以对彩色图像
void robert(string load_path, string save_path_robert) {
	Mat pScr;
	// 1为加载图像的颜色类型,1为原图返回,0为灰度返回
	pScr = imread(load_path, 1);
	int row = pScr.rows;
	int col = pScr.cols;
	int channel = pScr.channels();
	//imshow("picture", pScr);
	//waitKey(0);
	//destroyAllWindows();
	// 滤波窗口中心
	int pos = 1;
	vector<int> kernel_1 = { 0,0,0,0,-1,0,0,0,1 };
	vector<int> kernel_2 = { 0,0,0,0,0,-1,0,1,0 };
	// 从图片矩阵左上角开始从左到右从上到下使用滤波窗口处理图片矩阵,每次提取filter_size*filter_size的数据进入滤波窗口,计算平均值,然后替代滤波窗口中心值
	// 图片按每个通道来处理
	for (int k = 0; k < channel; k++) {
		// 每个滤波窗口中心位置都计算滤波窗口内通道值的平均值
		for (int m = pos; m < row - pos; m++) {
			for (int n = pos; n < col - pos; n++) {
				int centerpointValue_1 = 0;
				int centerpointValue_2 = 0;
				for (int i = m - 1; i <= m + 1; i++)
					for (int j = n - 1; j <= n + 1; j++)
						centerpointValue_1 += kernel_1[(i - (m - 1)) * 3 + (j - (n - 1))] * int(pScr.at<Vec3b>(i, j)[k]);
				centerpointValue_1 = abs(centerpointValue_1);
				for (int i = m - 1; i <= m + 1; i++)
					for (int j = n - 1; j <= n + 1; j++)
						centerpointValue_2 += kernel_2[(i - (m - 1)) * 3 + (j - (n - 1))] * int(pScr.at<Vec3b>(i, j)[k]);
				centerpointValue_2 = abs(centerpointValue_2);
				pScr.at<Vec3b>(m, n)[k] = centerpointValue_1 + centerpointValue_2;
				if (pScr.at<Vec3b>(m, n)[k] > 255)
					pScr.at<Vec3b>(m, n)[k] = 255;
				else if (pScr.at<Vec3b>(m, n)[k] < 0)
					pScr.at<Vec3b>(m, n)[k] = 0;
			}
		}
	}
	imwrite(save_path_robert, pScr);
}

// Sobel算子和Robert算子计算方法类似,也使用gx和gy的近似值代替x方向和y方向上的梯度,总梯度也使用|gx|+|gy|
// gx模板-1,0,1,-2,0,2,-1,0,1,gy模板-1,-2,-1,0,0,0,1,2,1
// 这两个模板作卷积矩阵前要旋转180度，即gx=1,0,-1,2,0,-2,1,0,-1,gy=1,2,1,0,0,0,-1,-2,-1
// sobel算子,既可以对灰度图像,也可以对彩色图像 
void sobel(string load_path, string save_path_sobel) {
	Mat pScr, grad_x, grad_y, abs_grad_x, abs_grad_y, grad;
	// 1为加载图像的颜色类型,1为原图返回,0为灰度返回
	pScr = imread(load_path, 1);
	int row = pScr.rows;
	int col = pScr.cols;
	int channel = pScr.channels();
	//imshow("picture", pScr);
	//waitKey(0);
	//destroyAllWindows();
	// 求x方向梯度
	// 第一个1,0表示在x方向求导,不在y方向求导,ksize=3是sobel算子大小,1是scale,0是delta
	Sobel(pScr, grad_x, CV_16S, 1, 0, 3, 1, 0, BORDER_DEFAULT);
	// 转成CV_8U型图像
	convertScaleAbs(grad_x, abs_grad_x);
	// 求y方向梯度
	Sobel(pScr, grad_y, CV_16S, 0, 1, 3, 1, 0, BORDER_DEFAULT);
	convertScaleAbs(grad_y, abs_grad_y);
	/// 合并梯度(近似)
	addWeighted(abs_grad_x, 0.5, abs_grad_y, 0.5, 0, grad);
	imwrite(save_path_sobel, grad);
}

// 锐化图像得s(x,y)。 
// 从原图像中减去锐化图像，产生的差值图像称为模板G(x,y)=f(x,y)-s(x,y)
// 将模板加到原图像中G(x,y)=p*f(x,y)-k*G(x,y)=(p-1)*f(x,y)+k*s(x,y)
// 当K=1时，为非锐化掩蔽;当k > 1时，高提升滤波，系数越大对细节的增强越明显;当k < 1时，不强调非锐化模板的贡献

//p是高提升滤波中原图像的混合比例,enhance_filter_size是高提升滤波模板的尺寸(长宽相等),k是模板的混合比例
void enhance_filter(string load_path, string save_path_enhance_filter, double p, int enhance_filter_size, double k, Mat kernel) {
	Mat pScr, res;
	// 1为加载图像的颜色类型,1为原图返回,0为灰度返回
	pScr = imread(load_path, 1);
	int row = pScr.rows;
	int col = pScr.cols;
	int channel = pScr.channels();
	//imshow("picture", pScr);
	//waitKey(0);
	//destroyAllWindows();
	// 使用上面的核中的权重对图像中每个像素点的像素值进行卷积计算时与下面拉普拉斯函数效果相同
	filter2D(pScr, res, pScr.depth(), kernel);
	// 滤波窗口中心
	int pos = (enhance_filter_size - 1) / 2;
	for (int k = 0; k < channel; k++) {
		for (int i = pos; i < row - pos; i++) {
			for (int j = pos; j < col - pos; j++) {
				res.at<Vec3b>(i, j)[k] = int((p - 1) * pScr.at<Vec3b>(i, j)[k] + k * res.at<Vec3b>(i, j)[k]);
			}
		}
	}
	normalize(res, res, 0, 255, NORM_MINMAX);
	imwrite(save_path_enhance_filter, res);
}

int main() {
	string load_path = "C:/Users/zgcr6/Desktop/高图实验/3/save/lenna.bmp";
	string save_path_gray_image = "C:/Users/zgcr6/Desktop/高图实验/3/save/lenna_gray.bmp";
	//gray_image(load_path, save_path_gray_image);
	string load_gray_path = save_path_gray_image;
	string save_path_mean_filter_gray_3 = "C:/Users/zgcr6/Desktop/高图实验/3/save/lenna_mean_filter_gray_3.bmp";
	string save_path_mean_filter_color_3 = "C:/Users/zgcr6/Desktop/高图实验/3/save/lenna_mean_filter_color_3.bmp";
	string save_path_mean_filter_gray_5 = "C:/Users/zgcr6/Desktop/高图实验/3/save/lenna_mean_filter_gray_5.bmp";
	string save_path_mean_filter_color_5 = "C:/Users/zgcr6/Desktop/高图实验/3/save/lenna_mean_filter_color_5.bmp";
	string save_path_mean_filter_gray_9 = "C:/Users/zgcr6/Desktop/高图实验/3/save/lenna_mean_filter_gray_9.bmp";
	string save_path_mean_filter_color_9 = "C:/Users/zgcr6/Desktop/高图实验/3/save/lenna_mean_filter_color_9.bmp";
	//mean_filter(load_gray_path, 3, save_path_mean_filter_gray_3);
	//mean_filter(load_path, 3, save_path_mean_filter_color_3);
	//mean_filter(load_gray_path, 5, save_path_mean_filter_gray_5);
	//mean_filter(load_path, 5, save_path_mean_filter_color_5);
	//mean_filter(load_gray_path, 9, save_path_mean_filter_gray_9);
	//mean_filter(load_path, 9, save_path_mean_filter_color_9);
	double sigma = 1.5;
	string save_path_gaussian_filter_gray_3 = "C:/Users/zgcr6/Desktop/高图实验/3/save/lenna_gaussian_filter_gray_3.bmp";
	string save_path_gaussian_filter_color_3 = "C:/Users/zgcr6/Desktop/高图实验/3/save/lenna_gaussian_filter_color_3.bmp";
	string save_path_gaussian_filter_gray_5 = "C:/Users/zgcr6/Desktop/高图实验/3/save/lenna_gaussian_filter_gray_5.bmp";
	string save_path_gaussian_filter_color_5 = "C:/Users/zgcr6/Desktop/高图实验/3/save/lenna_gaussian_filter_color_5.bmp";
	string save_path_gaussian_filter_gray_9 = "C:/Users/zgcr6/Desktop/高图实验/3/save/lenna_gaussian_filter_gray_9.bmp";
	string save_path_gaussian_filter_color_9 = "C:/Users/zgcr6/Desktop/高图实验/3/save/lenna_gaussian_filter_color_9.bmp";
	gaussian_filter(load_gray_path, 3, sigma, save_path_gaussian_filter_gray_3);
	gaussian_filter(load_path, 3, sigma, save_path_gaussian_filter_color_3);
	//separate_gaussian_filter(load_gray_path, 3, sigma, save_path_gaussian_filter_gray_3);
	//separate_gaussian_filter(load_path, 3, sigma, save_path_gaussian_filter_color_3);
	//separate_gaussian_filter(load_gray_path, 5, sigma, save_path_gaussian_filter_gray_5);
	//separate_gaussian_filter(load_path, 5, sigma, save_path_gaussian_filter_color_5);
	//separate_gaussian_filter(load_gray_path, 9, sigma, save_path_gaussian_filter_gray_9);
	//separate_gaussian_filter(load_path, 9, sigma, save_path_gaussian_filter_color_9);
	string save_path_laplacian_gray = "C:/Users/zgcr6/Desktop/高图实验/3/save/lenna_laplacian_gray.bmp";
	string save_path_laplacian_color = "C:/Users/zgcr6/Desktop/高图实验/3/save/lenna_laplacian_color.bmp";
	//laplacian(load_gray_path, save_path_laplacian_gray);
	//laplacian(load_path, save_path_laplacian_color);
	string save_path_robert_gray = "C:/Users/zgcr6/Desktop/高图实验/3/save/lenna_robert_gray.bmp";
	string save_path_robert_color = "C:/Users/zgcr6/Desktop/高图实验/3/save/lenna_robert_color.bmp";
	//robert(load_gray_path, save_path_robert_gray);
	//robert(load_path, save_path_robert_color);
	string save_path_sobel_gray = "C:/Users/zgcr6/Desktop/高图实验/3/save/lenna_sobel_gray.bmp";
	string save_path_sobel_color = "C:/Users/zgcr6/Desktop/高图实验/3/save/lenna_sobel_color.bmp";
	sobel(load_gray_path, save_path_sobel_gray);
	sobel(load_path, save_path_sobel_color);
	Mat kernel = (Mat_<float>(3, 3) << 1, 1, 1, 1, -8, 1, 1, 1, 1);
	string save_path_enhance_filter_gray = "C:/Users/zgcr6/Desktop/高图实验/3/save/lenna_enhance_filter_gray.bmp";
	string save_path_enhance_filter_color = "C:/Users/zgcr6/Desktop/高图实验/3/save/lenna_enhance_filter_color.bmp";
	//enhance_filter(load_gray_path, save_path_enhance_filter_gray, 2.0, 3, 2.0, kernel);
	//enhance_filter(load_path, save_path_enhance_filter_color, 1.0, 3, 2.0, kernel);
	return 0;
}

实验四：图像去噪

实验四内容

1、均值滤波
具体内容：利用OpenCV对灰度图像像素进行操作，分别利用算术均值滤波器、几何均值滤波器、谐波和逆谐波均值滤波器进行图像去噪。模板大小为5x5。（注：请分别为图像添加高斯噪声、胡椒噪声、盐噪声和椒盐噪声，并观察滤波效果）。
2、中值滤波
具体内容：利用OpenCV对灰度图像像素进行操作，分别利用5x5和9x9尺寸的模板对图像进行中值滤波。（注：请分别为图像添加胡椒噪声、盐噪声和椒盐噪声，并观察滤波效果）。
3、自适应均值滤波。
具体内容：利用OpenCV对灰度图像像素进行操作，设计自适应局部降低噪声滤波器去噪算法。模板大小7x7（对比该算法的效果和均值滤波器的效果）。
4、自适应中值滤波
具体内容：利用OpenCV对灰度图像像素进行操作，设计自适应中值滤波算法对椒盐图像进行去噪。模板大小7x7（对比中值滤波器的效果）。
5、彩色图像均值滤波
具体内容：利用OpenCV对彩色图像RGB三个通道的像素进行操作，利用算术均值滤波器和几何均值滤波器进行彩色图像去噪。模板大小为5x5。

实验四理论知识

椒噪声、盐噪声、椒盐噪声:
椒噪声即确定一个百分比，然后从图像中所有像素点中随机取出这个百分比数量的像素点，令其灰度值为0（黑色）。
盐噪声即确定一个百分比，然后从图像中所有像素点中随机取出这个百分比数量的像素点，令其灰度值为255（白色）。
椒盐噪声即确定一个百分比，然后从图像中所有像素点中随机取出这个百分比数量的像素点，令其灰度值随机为0（黑色）或255（白色）中的一个值。
高斯噪声:
Box–Muller变换:
基于这种变换，我们便可以得到一个从均匀分布中得到正态分布采样的算法。选取两个服从[0, 1]上均匀分布的随机变量U1和U2，有:
$$X=\sqrt{-2\ln U_1}\cos \left(2\pi U_2\right)$$
$$Y=\sqrt{-2\ln U_1}\sin \left(2\pi U_2\right)$$
则X与Y服从均值为0，方差为1的高斯分布。然后我们输入σ和μ，求Yxσ+μ，结果就是一个服从均值为μ，方差为σ2的高斯分布的噪声。原图每个像素点的像素值加上高斯噪声即为加噪后的图象。
算术均值滤波公式:
$$f(x,y)=\frac{1}{mn}\sum _{(x,y)\in S_{xy}}g(s,t)$$
Sxy表示中心点在（x,y）处，大小为m×n的滤波器窗口。g（s,t）表示原始图像灰度值，f（x,y）表示均值滤波后得到的图像灰度值。均值滤波就是把窗口内的像素点上灰度值全部按权重1相加求和，再除以像素点个数，得到的灰度值就是窗口正中心像素点的灰度值。对于彩色图像，对其三个通道分别执行算术均值滤波即可。
几何均值滤波公式:
$$f(x,y)={\left[\prod _{(s,t)\in S_{xy}g(s,t)}\right]}^{\frac{1}{mn}}$$
和算术均值滤波器相比，几何均值滤波器能够更好的取出高斯噪声，并且能够更多的保留图像的边缘信息。但，其对0值是非常敏感的，在滤波器的窗口内只要有一个像素的灰度值为0，就会造成滤波器的输出结果为0。
实际实现时如果像素点灰度值为0则不乘，开方的幂次也要相应减1。对于彩色图像，对其三个通道分别执行几何均值滤波即可。
谐波均值滤波公式:
$$f(x,y)=\frac{mn}{{\sum }_{(x,y)\in S_{xy}}\frac{1}{g(s,t)}}$$
谐波均值滤波器对盐粒噪声（白噪声）效果较好，不适用于胡椒噪声，比较适合处理高斯噪声。
逆谐波均值滤波公式:
$$f(x,y)=\frac{{\sum }_{(x,y)\in S_{xy}}g(s,t{)}^{Q+1}}{{\sum }_{(x,y)\in S_{xy}}g(s,t{)}^Q}$$
Q称为滤波器的阶数，该滤波器可以用来消除椒盐噪声。但是需要不同同时处理盐粒噪声和胡椒噪声，当Q为正时，可以消除胡椒噪声；当Q为负时，消除盐粒噪声。当Q=0时，该滤波器退化为算术均值滤波器；Q=-1时，退化为谐波均值滤波器。
中值滤波公式:
对mxn大小窗口内所有像素点灰度值按从小到大排序，取中值为窗口正中心的像素点的灰度值。
中值滤波对脉冲噪声（如椒盐噪声）有良好的滤除作用，特别是在滤除噪声的同时，能够保护信号的边缘，使之不被模糊。
自适应均值滤波公式:
先计算4个量:
g（x,y）带噪图像在点（x,y）上的值；
噪声的方差σ2；
Sxy区域内,像素点的局部均值mL；
Sxy区域内,像素点的局部方差S2；
则自适应均值计算公式为:
$$f(x,y)=g(x,y)-\frac{{\sigma }^2}{L^2}\left[g(x,y)-m_L\right]$$
自适应中值滤波公式:
自适应中值滤波器运行两个进程A,B；
过程A确定当前窗口内得到的中值Zmed是否是噪声。如果Zmin<Zmed<Zmax，则中值Zmed不是噪声，这时转到过程B测试，当前窗口的中心位置的像素Zxy是否是一个噪声点。如果Zmin<Zxy<Zmax，则Zxy不是一个噪声，此时滤波器输出Zxy；如果不满足上述条件，则可判定Zxy是噪声，这是输出中值Zmed（在A中已经判断出Zmed不是噪声）。
如果在过程A中，得到则Zmed不符合条件Zmin<Zmed<Zmax，则可判断得到的中值Zmed是一个噪声。在这种情况下，需要增大滤波器的窗口尺寸，在一个更大的范围内寻找一个非噪声点的中值，直到找到一个非噪声的中值，跳转到B；或者，窗口的尺寸达到了最大值，这时返回找到的中值，退出。

实验四代码

#include "pch.h"
#include <iostream>
#include <string>
#include <math.h>
#include <algorithm>
#include <opencv2/opencv.hpp>
#include <opencv2/core/core.hpp> 
#include <opencv2/highgui/highgui.hpp> //包含imread, imshow等
#include <opencv2/imgproc/imgproc.hpp> //包含cvtColor等

using namespace std;
using namespace cv;

void gray_image(string load_path, string save_path_gray_image) {
	Mat graypScr;
	// 1为加载图像的颜色类型,1为原图返回,0为灰度返回
	graypScr = imread(load_path, 0);
	//imshow("gray_picture", graypScr);
	//waitKey(0);
	//destroyAllWindows();
	imwrite(save_path_gray_image, graypScr);
}

// 盐噪声其实就是使图像的一些随机的像素为白色(255),胡椒噪声就是使图像的一些随机的像素为黑色(0)
// 盐椒噪声就是在使图像中一些随机的像素为黑色(0)或白色(255)
// 盐噪声
void salt(string load_path, string save_path_salt, double percent) {
	Mat pScr;
	// 1为加载图像的颜色类型,1为原图返回,0为灰度返回
	pScr = imread(load_path, 1);
	int row = pScr.rows;
	int col = pScr.cols;
	int channel = pScr.channels();
	int num = int(percent*row*col);
	for (int n = 0; n < num; n++) {
		// i,j为像素点坐标
		int i = rand() % row;
		int j = rand() % col;
		for (int k = 0; k < channel; k++) {
			pScr.at<Vec3b>(i, j)[k] = 255;
		}
	}
	imwrite(save_path_salt, pScr);
}

// 胡椒噪声
void pepper(string load_path, string save_path_pepper, double percent) {
	Mat pScr;
	// 1为加载图像的颜色类型,1为原图返回,0为灰度返回
	pScr = imread(load_path, 1);
	int row = pScr.rows;
	int col = pScr.cols;
	int channel = pScr.channels();
	int num = int(percent*row*col);
	for (int n = 0; n < num; n++) {
		// i,j为像素点坐标
		int i = rand() % row;
		int j = rand() % col;
		for (int k = 0; k < channel; k++) {
			pScr.at<Vec3b>(i, j)[k] = 0;
		}
	}
	imwrite(save_path_pepper, pScr);
}

// 椒盐噪声
void peppersalt(string load_path, string save_path_peppersalt, double percent) {
	Mat pScr;
	// 1为加载图像的颜色类型,1为原图返回,0为灰度返回
	pScr = imread(load_path, 1);
	int row = pScr.rows;
	int col = pScr.cols;
	int channel = pScr.channels();
	int num = int(percent*row*col);
	for (int n = 0; n < num; n++) {
		// i,j为像素点坐标
		int i = rand() % row;
		int j = rand() % col;
		int flag = rand() % 2;
		for (int k = 0; k < channel; k++) {
			if (flag == 0)
				pScr.at<Vec3b>(i, j)[k] = 0;
			else if (flag == 1)
				pScr.at<Vec3b>(i, j)[k] = 255;
		}
	}
	imwrite(save_path_peppersalt, pScr);
}

// Box–Muller变换是一种从一个均匀分布中得到正态分布采样的算法
// 选取两个服从[0, 1]上均匀分布的随机变量U1和U2,X=cos(2πU1)*sqrt(−2lnU2),Y=sin(2πU1)*sqrt(−2lnU2)
// 则X与Y服从均值为0，方差为1的高斯分布
// 生成高斯噪声,mu是均值,sigma是方差,X服从N(0,1)分布 
int generategaussiannoise(double mu, double sigma) {
	//定义一个特别小的值
	const double epsilon = numeric_limits<double>::min();//返回目标数据类型能表示的最逼近1的正数和1的差的绝对值
	static double z0, z1;
	static bool flag = false;
	flag = !flag;
	//flag为假，构造高斯随机变量
	if (!flag)
		return int(z1 * sigma + mu);
	double u1, u2;
	//构造随机变量
	do {
		u1 = rand()*(1.0 / RAND_MAX);
		u2 = rand()*(1.0 / RAND_MAX);
	} while (u1 <= epsilon);
	//flag为真构造高斯随机变量X
	z0 = sqrt(-2.0*log(u1))*cos(2 * CV_PI * u2);
	z1 = sqrt(-2.0*log(u1))*sin(2 * CV_PI * u2);
	return int(z1 * sigma + mu);
}

// 原图每个像素点的像素值加上高斯噪声即为加噪后的图象
// k可取16,32,64,128,256
// sigma可取0.1,0.5,1.0,2.0,5.0
// k为高斯噪声的系数,系数越大,高斯噪声越强;噪声服从高斯分布,所以方差越大,数据越分散,噪声也就越多
void addgaussiannoise(string load_path, string save_path_gaussiannoise, double mu, double sigma, int k) {
	Mat pScr;
	// 1为加载图像的颜色类型,1为原图返回,0为灰度返回
	pScr = imread(load_path, 1);
	int row = pScr.rows;
	int col = pScr.cols;
	int channel = pScr.channels();
	for (int k = 0; k < channel; k++) {
		for (int i = 0; i < row; i++) {
			for (int j = 0; j < col; j++) {
				pScr.at<Vec3b>(i, j)[k] = pScr.at<Vec3b>(i, j)[k] + generategaussiannoise(mu, sigma) * k;
				if (pScr.at<Vec3b>(i, j)[k] > 255)
					pScr.at<Vec3b>(i, j)[k] = 255;
				else if (pScr.at<Vec3b>(i, j)[k] < 0)
					pScr.at<Vec3b>(i, j)[k] = 0;
			}
		}
	}
	imwrite(save_path_gaussiannoise, pScr);
}

// 算术平均值滤波,既可以对灰度图像,也可以对彩色图像
void arithmetic_mean_filter(string load_path, int arithmetic_mean_filter_size, string save_path_arithmetic_mean_filter) {
	// 滤波窗口尺寸必须为奇数
	Mat pScr;
	// 1为加载图像的颜色类型,1为原图返回,0为灰度返回
	pScr = imread(load_path, 1);
	int row = pScr.rows;
	int col = pScr.cols;
	int channel = pScr.channels();
	//imshow("picture", pScr);
	//waitKey(0);
	//destroyAllWindows();
	// 滤波窗口中心
	int pos = (arithmetic_mean_filter_size - 1) / 2;
	// 从图片矩阵左上角开始从左到右从上到下使用滤波窗口处理图片矩阵,每次提取filter_size*filter_size的数据进入滤波窗口,计算平均值,然后替代滤波窗口中心值
	// 图片按每个通道来处理
	for (int k = 0; k < channel; k++) {
		// 每个滤波窗口中心位置都计算滤波窗口内通道值的平均值
		for (int m = pos; m < row - pos; m++) {
			for (int n = pos; n < col - pos; n++) {
				int centerpointValue = 0;
				for (int i = m - (arithmetic_mean_filter_size - 1) / 2; i <= m + (arithmetic_mean_filter_size - 1) / 2; i++)
					for (int j = n - (arithmetic_mean_filter_size - 1) / 2; j <= n + (arithmetic_mean_filter_size - 1) / 2; j++)
						centerpointValue += int(pScr.at<Vec3b>(i, j)[k]);
				centerpointValue = int(centerpointValue / (arithmetic_mean_filter_size*arithmetic_mean_filter_size));
				pScr.at<Vec3b>(m, n)[k] = centerpointValue;
			}
		}
	}
	imwrite(save_path_arithmetic_mean_filter, pScr);
}

// 几何平均值滤波,既可以对灰度图像,也可以对彩色图像
void geometric_mean_filter(string load_path, int geometric_mean_filter_size, string save_path_geometric_mean_filter) {
	// 滤波窗口尺寸必须为奇数
	Mat pScr;
	// 1为加载图像的颜色类型,1为原图返回,0为灰度返回
	pScr = imread(load_path, 1);
	int row = pScr.rows;
	int col = pScr.cols;
	int channel = pScr.channels();
	//imshow("picture", pScr);
	//waitKey(0);
	//destroyAllWindows();
	// 滤波窗口中心
	int pos = (geometric_mean_filter_size - 1) / 2;
	// 从图片矩阵左上角开始从左到右从上到下使用滤波窗口处理图片矩阵,每次提取filter_size*filter_size的数据进入滤波窗口,计算平均值,然后替代滤波窗口中心值
	// 图片按每个通道来处理
	for (int k = 0; k < channel; k++) {
		// 每个滤波窗口中心位置都计算滤波窗口内通道值的平均值
		for (int m = pos; m < row - pos; m++) {
			for (int n = pos; n < col - pos; n++) {
				double centerpointValue = 1.0;
				int count = 0;
				for (int i = m - (geometric_mean_filter_size - 1) / 2; i <= m + (geometric_mean_filter_size - 1) / 2; i++)
					for (int j = n - (geometric_mean_filter_size - 1) / 2; j <= n + (geometric_mean_filter_size - 1) / 2; j++) {
						if (pScr.at<Vec3b>(i, j)[k] != 0) {
							centerpointValue *= int(pScr.at<Vec3b>(i, j)[k]);
							count++;
						}
					}
				centerpointValue = pow(centerpointValue, 1.0 / count);
				pScr.at<Vec3b>(m, n)[k] = int(centerpointValue);
			}
		}
	}
	imwrite(save_path_geometric_mean_filter, pScr);
}

// 谐波均值滤波
void harmonic_mean_filter(string load_path, int harmonic_mean_filter_size, string save_path_harmonic_mean_filter) {
	// 滤波窗口尺寸必须为奇数
	Mat pScr;
	// 1为加载图像的颜色类型,1为原图返回,0为灰度返回
	pScr = imread(load_path, 1);
	int row = pScr.rows;
	int col = pScr.cols;
	int channel = pScr.channels();
	int pos = (harmonic_mean_filter_size - 1) / 2;
	// 从图片矩阵左上角开始从左到右从上到下使用滤波窗口处理图片矩阵,每次提取filter_size*filter_size的数据进入滤波窗口,计算平均值,然后替代滤波窗口中心值
	// 图片按每个通道来处理
	for (int k = 0; k < channel; k++) {
		// 每个滤波窗口中心位置都计算滤波窗口内通道值的平均值
		for (int m = pos; m < row - pos; m++) {
			for (int n = pos; n < col - pos; n++) {
				double centerpointValue = 0.0;
				int count = 0;
				for (int i = m - (harmonic_mean_filter_size - 1) / 2; i <= m + (harmonic_mean_filter_size - 1) / 2; i++)
					for (int j = n - (harmonic_mean_filter_size - 1) / 2; j <= n + (harmonic_mean_filter_size - 1) / 2; j++) {
						if (pScr.at<Vec3b>(i, j)[k] != 0) {
							centerpointValue += 1.0 / int(pScr.at<Vec3b>(i, j)[k]);
							count++;
						}
					}
				centerpointValue = count / centerpointValue;
				pScr.at<Vec3b>(m, n)[k] = int(centerpointValue);
			}
		}
	}
	imwrite(save_path_harmonic_mean_filter, pScr);
}

// 逆谐波均值滤波
// Q为滤波器的阶数,当Q为正时，可以消除胡椒噪声;当Q为负时,消除盐粒噪声;当Q=0时，该滤波器退化为算术均值滤波器;Q = -1时,退化为谐波均值滤波器
void inverse_harmonic_mean_filter(string load_path, int inverse_harmonic_mean_filter_size, string save_path_inverse_harmonic_mean_filter, double q) {
	// 滤波窗口尺寸必须为奇数
	Mat pScr;
	// 1为加载图像的颜色类型,1为原图返回,0为灰度返回
	pScr = imread(load_path, 1);
	int row = pScr.rows;
	int col = pScr.cols;
	int channel = pScr.channels();
	int pos = (inverse_harmonic_mean_filter_size - 1) / 2;
	// 从图片矩阵左上角开始从左到右从上到下使用滤波窗口处理图片矩阵,每次提取filter_size*filter_size的数据进入滤波窗口,计算平均值,然后替代滤波窗口中心值
	// 图片按每个通道来处理
	for (int k = 0; k < channel; k++) {
		// 每个滤波窗口中心位置都计算滤波窗口内通道值的平均值
		for (int m = pos; m < row - pos; m++) {
			for (int n = pos; n < col - pos; n++) {
				double centerpointValue_1 = 0.0;
				double centerpointValue_2 = 0.0;
				for (int i = m - (inverse_harmonic_mean_filter_size - 1) / 2; i <= m + (inverse_harmonic_mean_filter_size - 1) / 2; i++)
					for (int j = n - (inverse_harmonic_mean_filter_size - 1) / 2; j <= n + (inverse_harmonic_mean_filter_size - 1) / 2; j++) {
						centerpointValue_1 += pow(pScr.at<Vec3b>(i, j)[k], q + 1.0);
						if (pScr.at<Vec3b>(i, j)[k] != 0)
							centerpointValue_2 += pow(pScr.at<Vec3b>(i, j)[k], q);
					}
				pScr.at<Vec3b>(m, n)[k] = int(centerpointValue_1 / centerpointValue_2);
			}
		}
	}
	imwrite(save_path_inverse_harmonic_mean_filter, pScr);
}

//中值滤波
void median_filter(string load_path, string save_path_median_filter, int median_filter_size) {
	// 滤波窗口尺寸必须为奇数
	Mat pScr;
	// 1为加载图像的颜色类型,1为原图返回,0为灰度返回
	pScr = imread(load_path, 1);
	int row = pScr.rows;
	int col = pScr.cols;
	int channel = pScr.channels();
	int pos = (median_filter_size - 1) / 2;
	// 从图片矩阵左上角开始从左到右从上到下使用滤波窗口处理图片矩阵,每次提取filter_size*filter_size的数据进入滤波窗口,计算平均值,然后替代滤波窗口中心值
	// 图片按每个通道来处理
	for (int k = 0; k < channel; k++) {
		// 每个滤波窗口中心位置都计算滤波窗口内通道值的平均值
		for (int m = pos; m < row - pos; m++) {
			for (int n = pos; n < col - pos; n++) {
				int centerpointValue = 0;
				vector<int> filter_value;
				for (int i = m - (median_filter_size - 1) / 2; i <= m + (median_filter_size - 1) / 2; i++)
					for (int j = n - (median_filter_size - 1) / 2; j <= n + (median_filter_size - 1) / 2; j++)
						filter_value.push_back(int(pScr.at<Vec3b>(i, j)[k]));
				// 快速排序
				sort(filter_value.begin(), filter_value.end());
				pScr.at<Vec3b>(m, n)[k] = filter_value[int(median_filter_size*median_filter_size / 2)];
			}
		}
	}
	imwrite(save_path_median_filter, pScr);
}

// 自适应均值滤波
void self_adapt_mean_filter(string load_path, string save_path_self_adapt_mean_filter, int self_adapt_mean_filter_size) {
	// 滤波窗口尺寸必须为奇数
	Mat pScr, dst;
	// 1为加载图像的颜色类型,1为原图返回,0为灰度返回
	pScr = imread(load_path, 1);
	int row = pScr.rows;
	int col = pScr.cols;
	int channel = pScr.channels();
	blur(pScr, dst, Size(7, 7));
	int pos = (self_adapt_mean_filter_size - 1) / 2;
	double Sn = 100.0;
	for (int k = 0; k < channel; k++) {
		// 每个滤波窗口中心位置都计算滤波窗口内通道值的平均值
		for (int m = pos; m < row - pos; m++) {
			for (int n = pos; n < col - pos; n++) {
				int Zxy = int(pScr.at<Vec3b>(m, n)[k]);
				int Zmed = int(dst.at<Vec3b>(m, n)[k]);
				double Sl = 0.0;
				int count = 0;
				for (int i = m - (self_adapt_mean_filter_size - 1) / 2; i <= m + (self_adapt_mean_filter_size - 1) / 2; i++) {
					for (int j = n - (self_adapt_mean_filter_size - 1) / 2; j <= n + (self_adapt_mean_filter_size - 1) / 2; j++) {
						int Sxy = pScr.at<Vec3b>(i, j)[k];
						Sl += pow(Sxy - Zmed, 2);
						count++;
					}
				}
				Sl = Sl / count;
				pScr.at<Vec3b>(m, n)[k] = (int)(Zxy - Sn / Sl * (Zxy - Zmed));

			}
		}
	}
	imwrite(save_path_self_adapt_mean_filter, pScr);
}

// 过程A的目的是确定当前窗口内得到的中值Zmed是否是噪声。如果Zmin<Zmed<Zmax,则中值Zmed不是噪声,此时转到过程B测试当前窗口的中心位置的像素Zxy是否是一个噪声点
// 如果Zmin<Zxy<Zmax,则Zxy不是一个噪声，此时滤波器输出Zxy;如果不满足上述条件，则可判定Zxy是噪声,这时输出中值Zmed(在A中已经判断出Zmed不是噪声)。
// 如果在过程A中，得到Zmed不符合条件Zmin<Zmed<Zmax,则可判断得到的中值Zmed是一个噪声
// 在这种情况下,需要增大滤波器的窗口尺寸,在一个更大的范围内寻找一个非噪声点的中值，直到找到一个非噪声的中值，跳转到B
// 或者，窗口的尺寸达到了最大值，这时返回找到的中值,退出
// 自适应中值滤波
void self_adapt_median_filter(string load_path, string save_path_self_adapt_median_filter, int self_adapt_median_filter_size) {
	// 滤波窗口尺寸必须为奇数
	Mat pScr;
	// 1为加载图像的颜色类型,1为原图返回,0为灰度返回
	pScr = imread(load_path, 1);
	int row = pScr.rows;
	int col = pScr.cols;
	int channel = pScr.channels();
	// 滤波器窗口的起始尺寸
	int pos = (self_adapt_median_filter_size - 1) / 2;
	for (int k = 0; k < channel; k++) {
		// 每个滤波窗口中心位置都计算滤波窗口内通道值的平均值
		for (int m = pos; m < row - pos; m++) {
			for (int n = pos; n < col - pos; n++) {
				int minValue = 255;
				int maxValue = 0;
				int medValue = 0;
				vector<int> filter_value;
				for (int i = m - (self_adapt_median_filter_size - 1) / 2; i <= m + (self_adapt_median_filter_size - 1) / 2; i++) {
					for (int j = n - (self_adapt_median_filter_size - 1) / 2; j <= n + (self_adapt_median_filter_size - 1) / 2; j++) {
						filter_value.push_back(int(pScr.at<Vec3b>(i, j)[k]));
						if (int(pScr.at<Vec3b>(i, j)[k]) > maxValue)
							maxValue = int(pScr.at<Vec3b>(i, j)[k]);
						if (int(pScr.at<Vec3b>(i, j)[k]) < minValue)
							minValue = int(pScr.at<Vec3b>(i, j)[k]);
					}
				}
				sort(filter_value.begin(), filter_value.end());
				medValue = filter_value[int(self_adapt_median_filter_size*self_adapt_median_filter_size / 2)];
				int fxy = int(pScr.at<Vec3b>(m, n)[k]);
				if (medValue > minValue && medValue < maxValue)
					if (fxy <= minValue || fxy >= maxValue)
						pScr.at<Vec3b>(m, n)[k] = medValue;
					else
						pScr.at<Vec3b>(m, n)[k] = medValue;
			}
		}
	}
	imwrite(save_path_self_adapt_median_filter, pScr);
}

int main() {
	string load_path = "C:/Users/zgcr6/Desktop/高图实验/4/save/lenna.bmp";
	string save_path_gray_image = "C:/Users/zgcr6/Desktop/高图实验/4/save/lenna_gray.bmp";
	gray_image(load_path, save_path_gray_image);
	string load_gray_path = save_path_gray_image;
	string save_path_gray_salt = "C:/Users/zgcr6/Desktop/高图实验/4/save/lenna_gray_salt.bmp";
	string save_path_color_salt = "C:/Users/zgcr6/Desktop/高图实验/4/save/lenna_color_salt.bmp";
	salt(load_gray_path, save_path_gray_salt, 0.15);
	salt(load_path, save_path_color_salt, 0.15);
	string save_path_gray_pepper = "C:/Users/zgcr6/Desktop/高图实验/4/save/lenna_gray_pepper.bmp";
	string save_path_color_pepper = "C:/Users/zgcr6/Desktop/高图实验/4/save/lenna_color_pepper.bmp";
	pepper(load_gray_path, save_path_gray_pepper, 0.15);
	pepper(load_path, save_path_color_pepper, 0.15);
	string save_path_gray_peppersalt = "C:/Users/zgcr6/Desktop/高图实验/4/save/lenna_gray_peppersalt.bmp";
	string save_path_color_peppersalt = "C:/Users/zgcr6/Desktop/高图实验/4/save/lenna_color_peppersalt.bmp";
	peppersalt(load_gray_path, save_path_gray_peppersalt, 0.15);
	peppersalt(load_path, save_path_color_peppersalt, 0.15);
	string save_path_gray_gaussiannoise = "C:/Users/zgcr6/Desktop/高图实验/4/save/lenna_gray_gaussiannoise.bmp";
	string save_path_color_gaussiannoise = "C:/Users/zgcr6/Desktop/高图实验/4/save/lenna_color_gaussiannoise.bmp";
	addgaussiannoise(load_gray_path, save_path_gray_gaussiannoise, 1.0, 20, 32);
	addgaussiannoise(load_path, save_path_color_gaussiannoise, 1.0, 20, 32);
	string save_path_arithmetic_mean_filter_salt_gray = "C:/Users/zgcr6/Desktop/高图实验/4/save/lenna_arithmetic_mean_filter_salt_gray.bmp";
	string save_path_arithmetic_mean_filter_pepper_gray = "C:/Users/zgcr6/Desktop/高图实验/4/save/lenna_arithmetic_mean_filter_pepper_gray.bmp";
	string save_path_arithmetic_mean_filter_peppersalt_gray = "C:/Users/zgcr6/Desktop/高图实验/4/save/lenna_arithmetic_mean_filter_peppersalt_gray.bmp";
	string save_path_arithmetic_mean_filter_gaussiannoise_gray = "C:/Users/zgcr6/Desktop/高图实验/4/save/lenna_arithmetic_mean_filter_gaussiannoise_gray.bmp";
	arithmetic_mean_filter(save_path_gray_salt, 5, save_path_arithmetic_mean_filter_salt_gray);
	arithmetic_mean_filter(save_path_gray_pepper, 5, save_path_arithmetic_mean_filter_pepper_gray);
	arithmetic_mean_filter(save_path_gray_peppersalt, 5, save_path_arithmetic_mean_filter_peppersalt_gray);
	arithmetic_mean_filter(save_path_gray_gaussiannoise, 5, save_path_arithmetic_mean_filter_gaussiannoise_gray);
	string save_path_geometric_mean_filter_salt_gray = "C:/Users/zgcr6/Desktop/高图实验/4/save/lenna_geometric_mean_filter_salt_gray.bmp";
	string save_path_geometric_mean_filter_pepper_gray = "C:/Users/zgcr6/Desktop/高图实验/4/save/lenna_geometric_mean_filter_pepper_gray.bmp";
	string save_path_geometric_mean_filter_peppersalt_gray = "C:/Users/zgcr6/Desktop/高图实验/4/save/lenna_geometric_mean_filter_peppersalt_gray.bmp";
	string save_path_geometric_mean_filter_gaussiannoise_gray = "C:/Users/zgcr6/Desktop/高图实验/4/save/lenna_geometric_mean_filter_gaussiannoise_gray.bmp";
	geometric_mean_filter(save_path_gray_salt, 5, save_path_geometric_mean_filter_salt_gray);
	geometric_mean_filter(save_path_gray_pepper, 5, save_path_geometric_mean_filter_pepper_gray);
	geometric_mean_filter(save_path_gray_peppersalt, 5, save_path_geometric_mean_filter_peppersalt_gray);
	geometric_mean_filter(save_path_gray_gaussiannoise, 5, save_path_geometric_mean_filter_gaussiannoise_gray);
	string save_path_harmonic_mean_filter_salt_gray = "C:/Users/zgcr6/Desktop/高图实验/4/save/lenna_harmonic_mean_filter_salt_gray.bmp";
	string save_path_harmonic_mean_filter_pepper_gray = "C:/Users/zgcr6/Desktop/高图实验/4/save/lenna_harmonic_mean_filter_pepper_gray.bmp";
	string save_path_harmonic_mean_filter_peppersalt_gray = "C:/Users/zgcr6/Desktop/高图实验/4/save/lenna_harmonic_mean_filter_peppersalt_gray.bmp";
	string save_path_harmonic_mean_filter_gaussiannoise_gray = "C:/Users/zgcr6/Desktop/高图实验/4/save/lenna_harmonic_mean_filter_gaussiannoise_gray.bmp";
	harmonic_mean_filter(save_path_gray_salt, 5, save_path_harmonic_mean_filter_salt_gray);
	harmonic_mean_filter(save_path_gray_pepper, 5, save_path_harmonic_mean_filter_pepper_gray);
	harmonic_mean_filter(save_path_gray_peppersalt, 5, save_path_harmonic_mean_filter_peppersalt_gray);
	harmonic_mean_filter(save_path_gray_gaussiannoise, 5, save_path_harmonic_mean_filter_gaussiannoise_gray);
	string save_path_inverse_harmonic_mean_filter_salt_gray = "C:/Users/zgcr6/Desktop/高图实验/4/save/lenna_inverse_harmonic_mean_filter_salt_gray.bmp";
	string save_path_inverse_harmonic_mean_filter_pepper_gray = "C:/Users/zgcr6/Desktop/高图实验/4/save/lenna_inverse_harmonic_mean_filter_pepper_gray.bmp";
	string save_path_inverse_harmonic_mean_filter_peppersalt_gray = "C:/Users/zgcr6/Desktop/高图实验/4/save/lenna_inverse_harmonic_mean_filter_peppersalt_gray.bmp";
	string save_path_inverse_harmonic_mean_filter_gaussiannoise_gray = "C:/Users/zgcr6/Desktop/高图实验/4/save/lenna_inverse_harmonic_mean_filter_gaussiannoise_gray.bmp";
	inverse_harmonic_mean_filter(save_path_gray_salt, 5, save_path_inverse_harmonic_mean_filter_salt_gray, -0.5);
	inverse_harmonic_mean_filter(save_path_gray_pepper, 5, save_path_inverse_harmonic_mean_filter_pepper_gray, 1.0);
	inverse_harmonic_mean_filter(save_path_gray_peppersalt, 5, save_path_inverse_harmonic_mean_filter_peppersalt_gray, 0);
	inverse_harmonic_mean_filter(save_path_gray_gaussiannoise, 5, save_path_inverse_harmonic_mean_filter_gaussiannoise_gray, -1);
	string save_path_median_filter_salt_gray_5 = "C:/Users/zgcr6/Desktop/高图实验/4/save/lenna_median_filter_salt_gray_5.bmp";
	string save_path_median_filter_pepper_gray_5 = "C:/Users/zgcr6/Desktop/高图实验/4/save/lenna_median_filter_pepper_gray_5.bmp";
	string save_path_median_filter_peppersalt_gray_5 = "C:/Users/zgcr6/Desktop/高图实验/4/save/lenna_median_filter_peppersalt_gray_5.bmp";
	median_filter(save_path_gray_salt, save_path_median_filter_salt_gray_5, 5);
	median_filter(save_path_gray_pepper, save_path_median_filter_pepper_gray_5, 5);
	median_filter(save_path_gray_peppersalt, save_path_median_filter_peppersalt_gray_5, 5);
	string save_path_median_filter_salt_gray_9 = "C:/Users/zgcr6/Desktop/高图实验/4/save/lenna_median_filter_salt_gray_9.bmp";
	string save_path_median_filter_pepper_gray_9 = "C:/Users/zgcr6/Desktop/高图实验/4/save/lenna_median_filter_pepper_gray_9.bmp";
	string save_path_median_filter_peppersalt_gray_9 = "C:/Users/zgcr6/Desktop/高图实验/4/save/lenna_median_filter_peppersalt_gray_9.bmp";
	median_filter(save_path_gray_salt, save_path_median_filter_salt_gray_9, 9);
	median_filter(save_path_gray_pepper, save_path_median_filter_pepper_gray_9, 9);
	median_filter(save_path_gray_peppersalt, save_path_median_filter_peppersalt_gray_9, 9);
	string save_path_self_adapt_mean_filter_gray_7 = "C:/Users/zgcr6/Desktop/高图实验/4/save/lenna_self_adapt_mean_filter_gray_7.bmp";
	arithmetic_mean_filter(load_gray_path, 7, save_path_arithmetic_mean_filter_salt_gray);
	self_adapt_mean_filter(load_gray_path, save_path_self_adapt_mean_filter_gray_7, 7);
	string save_path_median_filter_peppersalt_gray_7 = "C:/Users/zgcr6/Desktop/高图实验/4/save/lenna_median_filter_peppersalt_gray_7.bmp";
	string save_path_self_adapt_median_filter_gray_7 = "C:/Users/zgcr6/Desktop/高图实验/4/save/lenna_self_adapt_median_filter_gray_7.bmp";
	median_filter(save_path_gray_peppersalt, save_path_median_filter_peppersalt_gray_7, 7);
	self_adapt_median_filter(save_path_gray_peppersalt, save_path_self_adapt_median_filter_gray_7, 7);
	string save_path_arithmetic_mean_filter_salt_color = "C:/Users/zgcr6/Desktop/高图实验/4/save/lenna_arithmetic_mean_filter_salt_color.bmp";
	string save_path_arithmetic_mean_filter_pepper_color = "C:/Users/zgcr6/Desktop/高图实验/4/save/lenna_arithmetic_mean_filter_pepper_color.bmp";
	string save_path_arithmetic_mean_filter_peppersalt_color = "C:/Users/zgcr6/Desktop/高图实验/4/save/lenna_arithmetic_mean_filter_peppersalt_color.bmp";
	string save_path_arithmetic_mean_filter_gaussiannoise_color = "C:/Users/zgcr6/Desktop/高图实验/4/save/lenna_arithmetic_mean_filter_gaussiannoise_color.bmp";
	salt(load_path, save_path_color_salt, 0.15);
	pepper(load_path, save_path_color_pepper, 0.15);
	peppersalt(load_path, save_path_color_peppersalt, 0.15);
	addgaussiannoise(load_path, save_path_color_gaussiannoise, 1.0, 20, 32);
	//arithmetic_mean_filter(save_path_color_salt, 5, save_path_arithmetic_mean_filter_salt_color);
	//arithmetic_mean_filter(save_path_color_pepper, 5, save_path_arithmetic_mean_filter_pepper_color);
	//arithmetic_mean_filter(save_path_color_peppersalt, 5, save_path_arithmetic_mean_filter_peppersalt_color);
	//arithmetic_mean_filter(save_path_color_gaussiannoise, 5, save_path_arithmetic_mean_filter_gaussiannoise_color);
	string save_path_geometric_mean_filter_salt_color = "C:/Users/zgcr6/Desktop/高图实验/4/save/lenna_geometric_mean_filter_salt_color.bmp";
	string save_path_geometric_mean_filter_pepper_color = "C:/Users/zgcr6/Desktop/高图实验/4/save/lenna_geometric_mean_filter_pepper_color.bmp";
	string save_path_geometric_mean_filter_peppersalt_color = "C:/Users/zgcr6/Desktop/高图实验/4/save/lenna_geometric_mean_filter_peppersalt_color.bmp";
	string save_path_geometric_mean_filter_gaussiannoise_color = "C:/Users/zgcr6/Desktop/高图实验/4/save/lenna_geometric_mean_filter_gaussiannoise_color.bmp";
	//geometric_mean_filter(save_path_color_salt, 5, save_path_geometric_mean_filter_salt_color);
	//geometric_mean_filter(save_path_color_pepper, 5, save_path_geometric_mean_filter_pepper_color);
	//geometric_mean_filter(save_path_color_peppersalt, 5, save_path_geometric_mean_filter_peppersalt_color);
	//geometric_mean_filter(save_path_color_gaussiannoise, 5, save_path_geometric_mean_filter_gaussiannoise_color);
	return 0;
}

实验五：频域滤波

实验五内容

1、 灰度图像的DFT和IDFT
具体内容：利用OpenCV提供的cvDFT函数对图像进行DFT和IDFT变换。
2、利用理想高通和低通滤波器对灰度图像进行频域滤波
具体内容：利用cvDFT函数实现DFT，在频域上利用理想高通和低通滤波器进行滤波，并把滤波过后的图像显示在屏幕上（观察振铃现象），要求截止频率可输入。
3、利用布特沃斯高通和低通滤波器对灰度图像进行频域滤波
具体内容：利用cvDFT函数实现DFT ，在频域上进行利用布特沃斯高通和低通滤波器进行滤波，并把滤波过后的图像显示在屏幕上（观察振铃现象），要求截止频率和n可输入。

实验五理论知识

二维离散傅里叶变换（DFT）与逆变换（IDFT）

二维DFT公式:
$$F(u,v)=\sum _{x=0}^{M-1}\sum _{y=0}^{N-1}f(x,y)e^{-j2\pi \left(\frac{\mathrm{u}\mathrm{x}}{\mathrm{M}}+\frac{vy}{N}\right)}$$
二维IDFT公式:
$$f(x,y)=\frac{1}{MN}\sum _{u=0}^{M-1}\sum _{v=0}^{N-1}F(u,v)e^{j2\pi \left(\frac{\mathrm{ux}}{M}+\frac{vy}{N}\right)}$$
注意上面两个公式来自数字图像处理（第三版）这本书，网上有不少博客把这两个公式写反了。
其中（x，y）为时域图上坐标（即图像上像素点坐标），f（x,y）为（x，y）处的灰度值，（u，v）为频域图上坐标，F（u，v）为（u，v）处的幅值。M和N分别为图像长和宽上像素点个数。j是虚数。
二维DFT变换可分离为两次一维DFT变换，这是由二维DFT变换的可分离性决定的:
$$F(u,v)=\sum _{x=0}^{M-1}\left[\sum _{y=0}^{N-1}f(x,y)e^{-j2\pi vy/N}\right]e^{-j2\pi ux/M}$$
先对行（y变量）做变换:
$$F(x,v)=\sum _{y=0}^{N-1}f(x,y)e^{-j2\pi ux/N}$$
然后对列（x变量）进行变换:
$$F(u,v)=\sum _{x=0}^{M-1}F(x,v)e^{-j2\pi ux/M}$$
同理，二维IDFT变换也可分离为两次一维IDFT变换:
$$f(x,y)=\frac{1}{MN}\sum _{u=0}^{M-1}\left[\sum _{v=0}^{N-1}F(u,v)e^{j2\pi vy/N}\right]e^{j2\pi ux/M}$$

一维离散傅里叶变换（DFT）与逆变换（IDFT）

DFT和IDFT既可以指一维的离散傅里叶变换和逆变换，也可以指二维的离散傅里叶变换和逆变换。
一维DFT公式:
$$F(u)=\sum _{n=0}^{N-1}f(n)e^{-j2\pi un/N}$$
一维IDFT公式:
$$f(n)=\frac{1}{N}\sum _{n=0}^{N-1}F(u)e^{j2\pi un/N}$$
其中f(n)是一维的时间序列信号，长度为N，u为频率，F（u）为频率u对应的幅值。
注意上面两个公式来自数字图像处理（第三版）这本书，网上有不少博客把这两个公式写反了。
j是一个虚数单位:
$$j^2=-1$$
一个复数z可以表示为z=a+bi（a,b∈R），其中a为实部，b为虚部，i为虚数单位。
欧拉公式:
$$e^{j\theta }=\cos \theta +j\sin \theta$$
$$e^{-j\theta }=\cos \theta -j\sin \theta$$
由上面欧拉公式，一维DFT公式可写为:
$$F(u)=\sum _{n=0}^{N-1}f(n)\left(\cos \frac{2\pi un}{N}-j\sin \frac{2\pi un}{N}\right)$$

快速一维傅里叶变换（FFT）

一维DFT公式为:
$$F(u)=\sum _{n=0}^{N-1}f(n)e^{-j2\pi un/N}$$
注意从这里开始所有的N都默认为2的整数次幂。可以看出，上面的公式时间复杂度为
$$O\left(n^2\right)$$
当N较大时，这个计算量是很大的。
令
$$W_N=e^{-j\frac{2\pi }{N}}$$
则一维DFT公式写为:
$$F(u)=\sum _{n=0}^{N-1}f(n)W_N^{nu}$$
一维傅里叶变换的周期性:
离散傅里叶变换DFT和它的里变换都以傅里叶变换的点数N为周期的。即:
$$F(u)=F(u\pm kN)$$
我们很容易证明下列结论:
$$W_N^{nu}=W_N^{(n+rN)u}=W_N^{(u+rN)n},\text{r为整数}$$
$$W_N^{nu}=W_{mN}^{mnu},W_N^{nu}=W_{N/m}^{nu/m}$$
我们还可以求出一些特殊值:
$$W_N^{N/2}=-1$$
$$W_N^{k+\frac{N}{2}}=-W_N^k$$
$$W_N^{(N-k)n}=W_N^{(N-n)k}=W_N^{-nk}$$
假设N是2的整数次幂，FFT的思想是每次都把序列分为奇序列与偶序列。x1(n)和x2(n)的长度都是N/2，x1(n)是偶数序列，x2(n)是奇数序列，
$$\{\begin{array}{l}x(2r)=x_1(r)\\ \\ x(2r+1)=x_2(r)\end{array}$$
$$r=0,1,\dots ,\frac{N}{2}-1$$
那么原序列可以表示为奇序列与偶序列之和:
$$x(n)=x_1(r)+x_2(r)$$
则有:
$$F(u)=\sum _{n=0(n为偶数)}^{N-1}x(2r)W_N^{2nu}+\sum _{n=0(n为奇数)}^{N-1}x(2r+1)W_N^{(2n+1)u}$$
$$F(u)=\sum _{r=0}^{\frac{N}{2}-1}{x}_1(r)W_N^{2ru}+\sum _{r=0}^{\frac{N}{2}-1}{x}_2(r)W_N^{(2r+1)u}$$
由上面周期性推导出的几个公式，可得:
$$F(u)=\sum _{r=0}^{\frac{N}{2}-1}{x}_1(r)W_{N/2}^{ru}+W_N^u\sum _{r=0}^{\frac{N}{2}-1}{x}_2(r)W_{N/2}^{ru}$$
令
$$F_1(u)=\sum _{r=0}^{\frac{N}{2}-1}{x}_1(r)W_{N/2}^{ru}$$
$$F_2(u)=\sum _{r=0}^{\frac{N}{2}-1}{x}_2(r)W_{N/2}^{ru}$$
则
$$F(u)=F_1(u)+W_N^u{F}_2(u)$$
其中
$$r=0,1,\dots ,\frac{N}{2}-1，u=0,1,\dots ,\frac{N}{2}-1$$
F1（u），F2（u）分别是x1（r）和x2（r）的N/2点DFT。这样我们就可以把一个N点DFT分解为两个N/2点DFT，并通过上式组合成N点DFT。式中的r:
$$r=0,1,\dots ,\frac{N}{2}-1，u=0,1,\dots ,\frac{N}{2}-1$$
也就是说，这只是F（u）的前一半序列。我们可以应用周期性求解F（u）的后一半序列。由前面周期性推导的几个公式，可以很容易地证明:
$$X_1\left(\frac{N}{2}+k\right)=X_1(k)$$
$$X_2\left(\frac{N}{2}+k\right)=X_2(k)$$
$$W_N^{k+\frac{N}{2}}=-W_N^k$$
于是有
$$F\left(u+\frac{N}{2}\right)=X_1(u)-W_N^u{X}_2(u),k=0,1,\dots ,\frac{N}{2}-1$$
这样我们就通过两个N/2点DFT序列组合出一个完整的N点DFT序列。
利用FFT求解DFT时， 复数乘法N2log2N次，复数加法Nlog2N次，故总时间复杂度为:
$$O\left(nlo{g}_{2}n\right)$$

一维傅里叶变换（DFT)与逆变换（IDFT）计算实例

令
$$W_N=e^{-j\frac{2\pi }{N}}$$
则一维DFT公式写为:
$$F(u)=\sum _{n=0}^{N-1}f(n)W_N^{nu}$$
假设有一个序列长度N=4，具体的序列x（n）=1，2，-1，3，n=0，1，2，3。
由N=4，得:
$$W=e^{-j\frac{2\pi }{4}}=\cos \left(\frac{\pi }{2}\right)-j\sin \left(\frac{\pi }{2}\right)=-j$$
于是有:
$$\left[\begin{array}{c}F(0)\\ \\ F(1)\\ \\ F(2)\\ \\ F(3)\end{array}\right]=\left[\begin{array}{cccc}{W}^0& W^0& W^0& W^0\\ & & & \\ W^0& W^1& W^2& W^3\\ & & & \\ W^0& W^2& W^4& W^6\\ & & & \\ W^0& W^3& W^6& W^9\end{array}\right]\left[\begin{array}{c}f(0)\\ \\ f(1)\\ \\ f(2)\\ \\ f(3)\end{array}\right]$$
$$=\left[\begin{array}{cccc}1& 1& 1& 1\\ & & & \\ 1& -j& -1& j\\ & & & \\ 1& -1& 1& -1\\ & & & \\ 1& j& -1& -j\end{array}\right]\left[\begin{array}{c}1\\ \\ 2\\ \\ -1\\ \\ 3\end{array}\right]=\left[\begin{array}{c}5\\ \\ 2+j\\ \\ -5\\ \\ 2-j\end{array}\right]$$
逆变换则为:
$$f(n)=\frac{1}{4}\left[\begin{array}{cccc}{W}^0& W^0& W^0& W^0\\ & & & \\ W^0& W^{-1}& W^{-2}& W^{-3}\\ & & & \\ W^0& W^{-2}& W^{-4}& W^{-6}\\ & & & \\ W^0& W^{-3}& W^{-6}& W^{-9}\end{array}\right]\left[\begin{array}{c}5\\ \\ 2+j\\ \\ -5\\ \\ 2-j\end{array}\right]$$
$$=\frac{1}{4}\left[\begin{array}{cccc}1& 1& 1& 1\\ & & & \\ 1& j& -1& -j\\ & & & \\ 1& -1& 1& -1\\ & & & \\ 1& 1& -j& j\end{array}\right]\left[\begin{array}{c}5\\ \\ 2+j\\ \\ -5\\ \\ 2-j\end{array}\right]=\left[\begin{array}{c}1\\ \\ 2\\ \\ -1\\ \\ 3\end{array}\right]$$

对一维傅里叶变换的含义解释

一维傅里叶变换就是将一个一维的时间序列信号分解成成无数个同时间序列的x同向、频率u不同的三角函数（即一维正交基），原信号f（x）可视作是e的复指数F（u）乘积的线性组合。
对于频谱图上的每一个频率值u都对应一个三角函数。每一个频率值u对应一个F（u）。F（u）是一个复数a+bi的形式，这个复数的模就是频率为u的三角函数的幅值（三角函数的最大/最小值），复数的角度arctan（b/a）就是三角函数的相位（相位应该在0到2pai之间）。
想象一下，现在我们输入一连串的时间序列f（x），其下标即时间点x，每个时间点x的值为f（x），以x为x轴，f（x）为y轴画一个直角坐标系图。我们的频谱图的频率长度和时间序列长度一样，但其x轴的含义是频率u。每一个频率u都可以确定一个三角函数，这个三角函数也可以画在时间序列的x和y轴上，三角函数的x轴也是时间点，其y值由两个方面确定：F（u）是一个复数a+bi，这个复数的模确定三角函数的幅值，同时arctan（b/a）就是三角函数的相位（相位应该在0到2pai之间）。频谱图上每一个u对应的三角函数以e的复指数与F（u）乘积的线性组合的形式叠加起来，就是我们的时间序列f（x）。

对二维傅里叶变换的含义解释

二维傅里叶变换就是将一个二维的时间序列信号（时间点为（x，y））分解成无数个三角平面波之和。一个三角平面波也可以看成由两个三角函数所构成的乘积，因此我们也可以说一个二维时间序列信号可以表示为无数个x向正余弦函数与y向正余弦函数乘积（即二维正交基）的线性组合。
三角平面波就是一个在二维直角坐标系（x，y）的三角函数向z轴方向伸展的结果。确定一个三角平面波需要四个参数，除了频率，幅度，相位之外，此外还有一个方向参数。我们使用一个二维的频域图（u，v，F（u，v））来保存二维傅里叶变换的结果。频域图上的点（u，v） 就代表这个点上的频率对应的三角平面波的法向量n，这个法向量的模代表这个平面波的频率，这个点里面保存的F（u，v）复数a+bi的模就是三角平面波的幅度，arctan（b/a）就是三角平面波的相位。
在一维FT变换后，频域呈现对称性，即前半段代表（0，fs/2）频率，而后半段代表（-fs/2，0）频率。在二维中的u和v也是如此。因此为了方便，一般会将图像1和4象限、2和3象限进行互换，这样就将0频率移动到了频域图中心。
对于互换后的频域图，中心的低频（图像上亮白部分）贡献了图像的主体，周围高频（黑色部分）提供图像的细节和边缘。

小波变换

傅里叶变换只能获取一段时间的信号总体上包含哪些频率的成分，但是对各成分出现的时刻不知道。两段完全不相同的不规律信号可能经过傅里叶变换后会有相同的频谱图。所以，对于不规律信号，傅里叶变换有很大局限性。
对于不规律信号，只知道包含哪些频率成分是不够的，我们还想知道各个频率成分出现的时间、频率随时间变化的情况，各个时刻的瞬时频率及其幅值，这就是时频分析。小波变换得到的正是一段时间信号的时频图。
傅里叶变换把无限长的三角函数作为基函数，而小波变换使用有限长的会衰减的小波作为基函数，这个小波只在一小段时间内有突变值，其余时间都为0值。
傅里叶变换变量只有频率ω，小波变换有两个变量：尺度a（scale）和平移量τ（translation）。尺度a控制小波函数的伸缩，平移量τ控制小波函数的平移。尺度就对应于频率（反比），平移量τ就对应于时间。
每个小波变换都会有一个mother wavelet，我们称之为母小波，同时还有一个scaling function，中文是尺度函数，也被成为父小波。任何小波变换的basis函数，其实就是对这个母小波和父小波缩放和平移后的集合。
当伸缩、平移到时间信号与某个频率的小波突变处重合时，会相乘得到一个大的值。这样我们就可以不仅可以知道频率，还知道该频率在时域上的具体时间点。

理想低通滤波

理想低通滤波器在以频域图上以中心点（p，q）为圆心、D0（又叫截止频率）为半径的圆内，通过所有的频率，而在圆外截断所有的频率。即:
$$D(u,v)=\sqrt{{\left(u-\frac{p}{2}\right)}^2+{\left(v-\frac{q}{2}\right)}^2}$$
D（u, v）表示频域图上某个点（u，v）距离频域图圆心（p，q）的距离。
$$H(u,v)=\{\begin{array}{l}1,D(u,v)\le D_0\\ \\ 0,D(u,v)>D_0\end{array}$$
D0为截止频率。即先做傅里叶变换得到频域图，频域图上每个点的F（u，v）乘以对应的H（u，v）后得到新的F（u，v）。然后再做逆傅里叶变换得到原图像。理想低通滤波器可以达到模糊/平滑图像的效果。
振铃现象:
由于理想低通滤波在频域图的圆内外边缘的F（u，v）急剧变化，其高频信息全部丢失（图像的细节和边缘丢失），这样做逆傅里叶变换得到的图像就会产生振铃现象，即图像的灰度剧烈变化处产生的震荡，就好像钟被敲击后产生的空气震荡。

理想高通滤波

理想高通滤波器与理想低通滤波器正好相反，其在以频域图上以中心点（p，q）为圆心、D0（又叫截止频率）为半径的圆内，截断所有的频率，而在圆外通过所有的频率。即:
$$D(u,v)=\sqrt{{\left(u-\frac{p}{2}\right)}^2+{\left(v-\frac{q}{2}\right)}^2}$$
D（u, v）表示频域图上某个点（u，v）距离频域图圆心（p，q）的距离。
$$H(u,v)=\{\begin{array}{l}0,D(u,v)\le D_0\\ \\ 1,D(u,v)>D_0\end{array}$$
D0为截止频率。即先做傅里叶变换得到频域图，频域图上每个点的F（u，v）乘以对应的H（u，v）后得到新的F（u，v）。然后再做逆傅里叶变换得到原图像。理想高通滤波器可以提取物体边缘。

Butterworth低通滤波

$$H(u,v)=\frac{1}{1+{\left(\frac{D(u,v)}{D_0}\right)}^{2n}}$$
其中n称为阶数。D0为截止频率。1阶Butterworth低通滤波没有出现振铃现象，随着阶数增大，振铃现象越来越明显。
先做傅里叶变换得到频域图，频域图上每个点的F（u，v）乘以对应的H（u，v）后得到新的F（u，v）。然后再做逆傅里叶变换得到原图像。Butterworth低通滤波是对理想低通滤波的改进。

Butterworth高通滤波

$$H(u,v)=\frac{1}{1+{\left(\frac{D_0}{D(u,v)}\right)}^{2n}}$$
其中n称为阶数。D0为截止频率。先做傅里叶变换得到频域图，频域图上每个点的F（u，v）乘以对应的H（u，v）后得到新的F（u，v）。然后再做逆傅里叶变换得到原图像。Butterworth高通滤波是对理想高通滤波的改进。

实验五代码

#include "pch.h"
#include <iostream>
#include <string>
#include <math.h>
#include <algorithm>
#include <opencv2/opencv.hpp>
#include <opencv2/core/core.hpp> 
#include <opencv2/highgui/highgui.hpp> //包含imread, imshow等
#include <opencv2/imgproc/imgproc.hpp> //包含cvtColor等

using namespace std;
using namespace cv;

void gray_image(string load_path, string save_path_gray_image) {
	Mat graypScr;
	// 1为加载图像的颜色类型,1为原图返回,0为灰度返回
	graypScr = imread(load_path, 0);
	//imshow("gray_picture", graypScr);
	//waitKey(0);
	//destroyAllWindows();
	imwrite(save_path_gray_image, graypScr);
}

// F(u,v)即求频域函数,u,v代表空间频率，即像素值在x方向和y方向上的梯度,图像尺寸为M×N,傅里叶变换是从时域变换到频域
// 对于图像,其中物体的边界处频率较高,因为差异大
// 频谱图中图像最明亮的像素放到中央,然后逐渐变暗,在边缘上的像素最暗,这样可以发现图像中亮、暗像素的百分比，即为频域中的振幅A的强度
// 对二维傅里叶变换公式,t=(ux/m+vy/n),F(u,v)=f(x,y)exp(2*pi*t),x从0到m,y从0到n遍历
// 两个求和号对图像进行遍历,f(x,y)取出原像素的数值,当x=0,(横轴不动),对y进行遍历时
// t的后一项表示变换前像素的位置比例与变换后的位置相乘,映射到新的位置,且能够反映像素延y方向距离的差异,越靠后的像素(y越大)t值越大,即t能够反映出不同位置(纵轴)像素之间的差异,前一项含义为保留像素相对位置(横轴)的信息(遍历y时为常数),2π为修正参数
// f(x,y)代表当前像素点(x,y)像素值
// 在图像处理领域,通过DFT可以将图像转换到频域,实现高通和低通滤波,频域相当于是图像的特征
// 频谱图中心看到更多的白色区域,表示低振幅的波占到了多数,高振幅的波占少数,而且说明了波的整体频率偏低
// 傅立叶频域图上我们看到的明暗不一的亮点,实际上是图像上某一点与邻域点差异的强弱,即梯度的大小,也即该点的频率的大小
// 可以这么理解,图像中的低频部分指低梯度的点，高频部分相反
// 梯度大则该点的亮度强,否则该点亮度弱。这样通过观察傅立叶变换后的频域图，也叫功率图
// 如果频域图中暗的点数更多,那么实际图像是比较柔和的(因为各点与邻域差异都不大，梯度相对较小)
// 如果频谱图中亮的点数多,那么实际图像一定是尖锐的,边界分明且边界两边像素差异较大的
// 频域变换的用途:压缩和去噪
// 压缩的原理就是在频域中,大部分频域的值为0(或接近0,可以进行有损压缩,如jpeg图像),只要压缩频域中的少数非0值即可达到图片压缩的目的
// 去噪则是通过频域的滤波实现,因为噪声大部分情况下体现为高频信号,使用低通滤波器即可滤除高频噪声(当然,也会带来损失,那就是边缘会变得模糊,边缘也是高频信号)
// 二维的傅里叶变换(DFT)实际上是计算两次一维的FFT变换实现的,根据FFT变换的计算要求,需要图像的行数、列数均满足2的n次方,如果不满足,在计算FFT之前先要对图像补零以满足2的n次
// FFT是离散傅氏变换的快速算法,即为快速傅氏变换,FFT其实是一个用O(nlog⁡2n)的时间将一个用系数表示的多项式转换成它的点值表示的算法
// n次多项式的系数为f(x)={a0,a1,...,an-1},这样表示出来的多项式称为系数表示法
// 点值表示法即把n个不同的(x,y)点代入，把n条式子联立起来成为一个有n条方程的n元方程组，每一项的系数都可以解出来,这n个点唯一确定该多项式
// 多项式可用f(x)={(x0,y0),(x1,y1),...,(xn-1,yn-1)}表示,这就是点值表示法
// 系数表示法做多项式乘法时间复杂度O(n2),但两个用点值表示的多项式相乘，只需要O(n)的时间
// 设两个多项式为f(x)={(x0,f0),(x1,f1),...,(xn-1,fn-1)}和g(x)={(x0,g0),(x1,g1),...,(xn-1,gn-1)}
// 乘积h(x)={(x0,f0*g0),(x1,f1*g1),...,(xn-1,fn-1*gn-1)}
// FFT变换:设一个多项式A(x),按A(x)下标的奇偶性把A(x)分成两半，右边再提一个x,即左右分别设为A1(x)和A2(x),A(x)=A1(x2)+xA2(x2),然后可以用分治法来作
// 设一个序列为1,2,-1,3,N=4,若W=exp(-j*(2*pi/4),欧拉公式exp(-it)=cos(t)-i*sin(t),exp(it)=cos(t)+i*sin(t),故W=cos(pi/2)-j*sin(pi/2)=-j
// 一维DFT变换公式X(k)=x(n)W的n*k次方(W即上面计算出的-j),n从0到n-1求和;逆变换x(n)=X(k)W的-n*k次方(w=-j),n从0到n-1求和后乘以1/N,原始计算为exp(j*2*pi/N)
// 如X(0)=x(0)*W的0*0次方+x(1)*W的1*0次方+...+x(3)*W的3*0次方,X(3)=x(0)*W的0*3次方+...+x(3)*W的3*3次方,W=-j,W的0次方为1,W的平方为-1
void DFT(string load_path, string save_path_DFT) {
	Mat pScr;
	// 1为加载图像的颜色类型,1为原图返回,0为灰度返回,注意这里得到的pScr必须为灰度图格式,否则后面会报错,原因是如果读入RGB图像,padded也是三通道的RGB格式,那么Mat planes[]就会报错
	pScr = imread(load_path, 0);
	int row = pScr.rows;
	int col = pScr.cols;
	int channel = pScr.channels();
	//以0填充输入图像矩阵
	Mat padded;
	int m = getOptimalDFTSize(pScr.rows);
	int n = getOptimalDFTSize(pScr.cols);
	//扩充原图的边缘,将图像变大,然后以各种外插方式自动填充图像边界,比如均值滤波或者中值滤波中,使用copyMakeBorder将原图稍微放大,然后我们就可以处理边界的情况
	//填充输入图像I，输入矩阵为padded,上方和左方不做填充处理,padded与pScr的差集部分填充为0,其他部分同pScr,padded为我们的目标图像
	copyMakeBorder(pScr, padded, 0, m - pScr.rows, 0, n - pScr.cols, BORDER_CONSTANT, Scalar::all(0));
	// CV_32S等同于CV_32SC1,是单通道的,CV_32F取值范围为0-1.0，imshow的时候会把图像乘以255后再显示
	Mat planes[] = { Mat_<float>(padded), Mat::zeros(padded.size(),CV_32F) };
	Mat complexI;
	// 将planes融合合并成一个多通道数组complexI,2是需要合并的矩阵的个数
	merge(planes, 2, complexI);
	//进行傅里叶变换
	dft(complexI, complexI);
	//计算幅值,转换到对数尺度(logarithmic scale)
	//=> log(1 + sqrt(Re(DFT(I))^2 + Im(DFT(I))^2))
	// planes[0] = Re(DFT(I),planes[1] = Im(DFT(I)),即planes[0]为实部,planes[1]为虚部
	// 振幅和相位可以用实部和虚部表示出来
	split(complexI, planes);
	// 计算二维矢量的幅值,第一个参数表示矢量的浮点型X坐标值,也就是实部;第二个参数表示矢量的浮点型Y坐标值,也就是虚部;第三个参数是输出的幅值
	magnitude(planes[0], planes[1], planes[0]);
	Mat magI = planes[0];
	// 幅值都加上1再转换成对数尺度
	magI += Scalar::all(1);
	log(magI, magI);
	//如果有奇数行或列，则对频谱进行裁剪
	magI = magI(Rect(0, 0, magI.cols&-2, magI.rows&-2));
	//重新排列傅里叶图像中的象限，使得原点位于图像中心
	int cx = magI.cols / 2;
	int cy = magI.rows / 2;
	// 根据中心点划定象限
	Mat q0(magI, Rect(0, 0, cx, cy));       //左上角图像划定ROI区域
	Mat q1(magI, Rect(cx, 0, cx, cy));      //右上角图像
	Mat q2(magI, Rect(0, cy, cx, cy));      //左下角图像
	Mat q3(magI, Rect(cx, cy, cx, cy));     //右下角图像
	//变换左上角和右下角象限,即交换q0和q3的幅值
	Mat tmp;
	q0.copyTo(tmp);
	q3.copyTo(q0);
	tmp.copyTo(q3);
	//变换右上角和左下角象限,即交换q1和q2的幅值
	q1.copyTo(tmp);
	q2.copyTo(q1);
	tmp.copyTo(q2);
	// 最后频谱图应当归一化到0-1
	normalize(magI, magI, 0, 1, CV_MINMAX);
	//imshow("magnitude", magI);
	//waitKey(0);
	//destroyAllWindows();
	normalize(magI, magI, 0, 255, CV_MINMAX);
	imwrite(save_path_DFT, magI);
}

// 傅里叶变换后再逆变换得到原图
void DFT_IDFT(string load_path, string save_path_DFT_IDFT_1, string save_path_DFT_IDFT_2) {
	Mat pScr;
	// 1为加载图像的颜色类型,1为原图返回,0为灰度返回,注意这里得到的pScr必须为灰度图格式,否则后面会报错,原因是如果读入RGB图像,padded也是三通道的RGB格式,那么Mat planes[]就会报错
	pScr = imread(load_path, 0);
	int row = pScr.rows;
	int col = pScr.cols;
	int channel = pScr.channels();
	//以0填充输入图像矩阵
	Mat padded;
	int m = getOptimalDFTSize(pScr.cols);
	int n = getOptimalDFTSize(pScr.rows);//获取最佳尺寸，快速傅立叶变换要求尺寸为2的n次方
	//扩充原图的边缘,将图像变大,然后以各种外插方式自动填充图像边界,比如均值滤波或者中值滤波中,使用copyMakeBorder将原图稍微放大,然后我们就可以处理边界的情况
	//填充输入图像I，输入矩阵为padded,上方和左方不做填充处理,padded与pScr的差集部分填充为0,其他部分同pScr,padded为我们的目标图像
	copyMakeBorder(pScr, padded, 0, n - pScr.rows, 0, m - pScr.cols, BORDER_CONSTANT, Scalar::all(0));
	// CV_32S等同于CV_32SC1,是单通道的,CV_32F取值范围为0-1.0，imshow的时候会把图像乘以255后再显示
	Mat planes[] = { Mat_<float>(padded),Mat::zeros(padded.size(),CV_32F) };
	Mat complexIm;
	// 将planes融合合并成一个多通道数组complexI,2是需要合并的矩阵的个数
	merge(planes, 2, complexIm);
	//进行傅立叶变换，结果保存在自身
	dft(complexIm, complexIm);
	// 将多通道的Mat矩阵分离成单通道的序列
	split(complexIm, planes);
	//获取幅度图像，0通道为实数通道，1为虚数通道，因为二维傅立叶变换结果是复数
	// 计算二维矢量的幅值,第一个参数表示矢量的浮点型X坐标值,也就是实部;第二个参数表示矢量的浮点型Y坐标值,也就是虚部;第三个参数是输出的幅值
	magnitude(planes[0], planes[1], planes[0]);
	Mat magI = planes[0];
	// 幅值都加上1再转换成对数尺度
	magI += Scalar::all(1);
	log(magI, magI);
	//重新排列傅里叶图像中的象限，使得原点位于图像中心
	int cx = padded.cols / 2;
	int cy = padded.rows / 2;
	// 根据中心点划定象限
	Mat part1(planes[0], Rect(0, 0, cx, cy));
	Mat part2(planes[0], Rect(cx, 0, cx, cy));
	Mat part3(planes[0], Rect(0, cy, cx, cy));
	Mat part4(planes[0], Rect(cx, cy, cx, cy));
	Mat temp;
	//变换左上角和右下角象限,即交换p1和p4的幅值
	part1.copyTo(temp);
	part4.copyTo(part1);
	temp.copyTo(part4);
	//变换右上角和左下角象限,即交换p2和p3的幅值
	part2.copyTo(temp);
	part3.copyTo(part2);
	temp.copyTo(part3);
	// 最后频谱图应当归一化到0-1,然后显示频谱图
	normalize(magI, magI, 1, 0, CV_MINMAX);
	//imshow("magnitude", magI);
	//waitKey(0);
	//destroyAllWindows();
	normalize(magI, magI, 0, 255, CV_MINMAX);
	imwrite(save_path_DFT_IDFT_1, magI);

	// 频谱图逆变换回原图
	Mat _complexim;
	// 把complexIm的内容复制粘贴到_complexim
	complexIm.copyTo(_complexim);
	//创建两个通道，类型为float，大小为填充后的尺寸,0通道为实数通道，1为虚数通道，因为二维傅立叶变换结果是复数
	Mat iDft[] = { Mat::zeros(planes[0].size(),CV_32F),Mat::zeros(planes[0].size(),CV_32F) };
	//傅立叶逆变换,求出原来数值
	idft(_complexim, _complexim);
	// 分离成单通道序列
	split(_complexim, iDft);
	// 我们主要使用0通道的值,即实部值就是我们原有的图像的像素值
	magnitude(iDft[0], iDft[1], iDft[0]);
	normalize(iDft[0], iDft[0], 0, 1, CV_MINMAX);//归一化处理，float类型的显示范围为0-1
	// 显示逆变换后图片
	//imshow("idft_picture", iDft[0]);
	//waitKey(0);
	//destroyAllWindows();
	normalize(iDft[0], iDft[0], 0, 255, CV_MINMAX);
	imwrite(save_path_DFT_IDFT_2, iDft[0]);
}

Mat freqfilt(Mat &scr, Mat &blur) {
	//创建通道，存储dft后的实部与虚部（CV_32F，必须为单通道数）
	Mat plane[] = { scr, Mat::zeros(scr.size() , CV_32FC1) };
	Mat complexIm;
	//合并通道(把两个矩阵合并为一个2通道的Mat类容器)
	merge(plane, 2, complexIm);
	//进行傅立叶变换，结果保存在自身
	dft(complexIm, complexIm);
	//分离通道(数组分离)
	split(complexIm, plane);
	plane[0] = plane[0](Rect(0, 0, plane[0].cols & -2, plane[0].rows & -2));
	int cx = plane[0].cols / 2;
	int cy = plane[0].rows / 2;
	Mat part1_r(plane[0], Rect(0, 0, cx, cy));
	Mat part2_r(plane[0], Rect(cx, 0, cx, cy));
	Mat part3_r(plane[0], Rect(0, cy, cx, cy));
	Mat part4_r(plane[0], Rect(cx, cy, cx, cy));
	Mat temp;
	//左上与右下交换位置(实部)
	part1_r.copyTo(temp);
	part4_r.copyTo(part1_r);
	temp.copyTo(part4_r);
	//右上与左下交换位置(实部)
	part2_r.copyTo(temp);
	part3_r.copyTo(part2_r);
	temp.copyTo(part3_r);
	//元素坐标(cx,cy)
	Mat part1_i(plane[1], Rect(0, 0, cx, cy));
	Mat part2_i(plane[1], Rect(cx, 0, cx, cy));
	Mat part3_i(plane[1], Rect(0, cy, cx, cy));
	Mat part4_i(plane[1], Rect(cx, cy, cx, cy));
	//左上与右下交换位置(虚部)
	part1_i.copyTo(temp);
	part4_i.copyTo(part1_i);
	temp.copyTo(part4_i);
	//右上与左下交换位置(虚部)
	part2_i.copyTo(temp);
	part3_i.copyTo(part2_i);
	temp.copyTo(part3_i);
	// 滤波器函数与DFT结果的乘积
	Mat blur_r, blur_i, BLUR;
	//滤波(实部与滤波器模板对应元素相乘),plane[0]与blur逐个元素相乘
	multiply(plane[0], blur, blur_r);
	//滤波(虚部与滤波器模板对应元素相乘),plane[1]与blur逐个元素相乘
	multiply(plane[1], blur, blur_i);
	Mat plane1[] = { blur_r, blur_i };
	//实部与虚部合并
	merge(plane1, 2, BLUR);
	//得到原图频谱图
	//获取幅度图像，0通道为实部通道，1为虚部，因为二维傅立叶变换结果是复数
	magnitude(plane[0], plane[1], plane[0]);
	//傅立叶变换后的图片不好分析，进行对数处理，结果比较好看
	plane[0] += Scalar::all(1);
	log(plane[0], plane[0]);
	//归一化便于显示
	normalize(plane[0], plane[0], 0, 1, CV_MINMAX);
	//imshow("原图像频谱图",plane[0]);
	//waitKey(0);
	//destroyAllWindows();
	//idft结果也为复数
	idft(BLUR, BLUR);
	//分离通道，主要获取通道
	split(BLUR, plane);
	//求幅值(模)
	magnitude(plane[0], plane[1], plane[0]);
	//归一化便于显示
	normalize(plane[0], plane[0], 0, 1, CV_MINMAX);
	//返回参数
	return plane[0];
}

// 理想低通滤波模板
Mat ideal_low_pass_kernel(Mat &scr, float D0, string save_path_ideal_low_pass_kernel) {
	Mat ideal_low_pass(scr.size(), CV_32FC1);
	for (int i = 0; i < scr.rows; i++) {
		for (int j = 0; j < scr.cols; j++) {
			double d = sqrt(pow((i - scr.rows / 2), 2) + pow((j - scr.cols / 2), 2));
			if (d <= D0) {
				ideal_low_pass.at<float>(i, j) = 1.0;
			}
			else
				ideal_low_pass.at<float>(i, j) = 0.0;
		}
	}
	//imshow("ideal_low_pass_kernel", ideal_low_pass);
	//waitKey(0);
	//destroyAllWindows();
	// 归一化成255范围保存的图片和上面显示的效果一样
	normalize(ideal_low_pass, ideal_low_pass, 0, 255, CV_MINMAX);
	imwrite(save_path_ideal_low_pass_kernel, ideal_low_pass);
	return ideal_low_pass;
}

// 理想低通滤波,D0表示截止频率,D0是一个半径
// 半径D0越小,模糊越大;半径D0越大,模糊越小
// 中心点取频率矩形的中心,从点(u,v)到中心的距离为欧式距离D(u,v),如果距离小于等于D0,则H(u,v)=1,否则为0
// 振铃现象就是指输出图像的灰度剧烈变化处产生的震荡,就好像钟被敲击后产生的空气震荡
// 理想低通滤波器
void ideal_low_pass_filter(string load_path, float d0, string save_path_ideal_low_pass_kernel, string save_path_ideal_low_pass_filter) {
	Mat pScr;
	// 1为加载图像的颜色类型,1为原图返回,0为灰度返回,注意这里得到的pScr必须为灰度图格式,否则后面会报错,原因是如果读入RGB图像,padded也是三通道的RGB格式,那么Mat planes[]就会报错
	pScr = imread(load_path, 0);
	int row = pScr.rows;
	int col = pScr.cols;
	int channel = pScr.channels();
	int M = getOptimalDFTSize(pScr.rows);
	int N = getOptimalDFTSize(pScr.cols);
	Mat padded;
	copyMakeBorder(pScr, padded, 0, M - pScr.rows, 0, N - pScr.cols, BORDER_CONSTANT, Scalar::all(0));
	//将图像转换为float型
	padded.convertTo(padded, CV_32FC1);
	Mat ideal_low_kernel = ideal_low_pass_kernel(padded, d0, save_path_ideal_low_pass_kernel);//理想低通滤波器
	Mat result = freqfilt(padded, ideal_low_kernel);
	//imshow("ideal_low_pass_kernel_picture", result);
	//waitKey(0);
	//destroyAllWindows();
	normalize(result, result, 0, 255, CV_MINMAX);
	imwrite(save_path_ideal_low_pass_filter, result);
}

// 理想高通滤波模板
Mat ideal_high_pass_kernel(Mat &scr, float D0, string save_path_ideal_high_pass_kernel) {
	Mat ideal_high_pass(scr.size(), CV_32FC1);
	for (int i = 0; i < scr.rows; i++) {
		for (int j = 0; j < scr.cols; j++) {
			double d = sqrt(pow((i - scr.rows / 2), 2) + pow((j - scr.cols / 2), 2));
			if (d <= D0) {
				ideal_high_pass.at<float>(i, j) = 0.0;
			}
			else
				ideal_high_pass.at<float>(i, j) = 1.0;
		}
	}
	//imshow("ideal_high_pass_kernel", ideal_high_pass);
	//waitKey(0);
	//destroyAllWindows();
	// 归一化成255范围保存的图片和上面显示的效果一样
	normalize(ideal_high_pass, ideal_high_pass, 0, 255, CV_MINMAX);
	imwrite(save_path_ideal_high_pass_kernel, ideal_high_pass);
	return ideal_high_pass;
}

// 理想高通滤波,和理想低通滤波正好相反,中心点取频率矩形的中心,从点(u,v)到中心的距离为欧式距离D(u,v),如果距离小于等于D0,则H(u,v)=0,否则为1
// 高通滤波后图像得到锐化处理,图像中仅剩下边缘
// 理想高通滤波器
void ideal_high_pass_filter(string load_path, float d0, string save_path_ideal_high_pass_kernel, string save_path_ideal_high_pass_filter) {
	Mat pScr;
	// 1为加载图像的颜色类型,1为原图返回,0为灰度返回,注意这里得到的pScr必须为灰度图格式,否则后面会报错,原因是如果读入RGB图像,padded也是三通道的RGB格式,那么Mat planes[]就会报错
	pScr = imread(load_path, 0);
	int row = pScr.rows;
	int col = pScr.cols;
	int channel = pScr.channels();
	int M = getOptimalDFTSize(pScr.rows);
	int N = getOptimalDFTSize(pScr.cols);
	Mat padded;
	copyMakeBorder(pScr, padded, 0, M - pScr.rows, 0, N - pScr.cols, BORDER_CONSTANT, Scalar::all(0));
	//将图像转换为float型
	padded.convertTo(padded, CV_32FC1);
	Mat ideal_high_kernel = ideal_high_pass_kernel(padded, d0, save_path_ideal_high_pass_kernel);//理想低通滤波器
	Mat result = freqfilt(padded, ideal_high_kernel);
	//imshow("ideal_high_pass_kernel_picture", result);
	//waitKey(0);
	//destroyAllWindows();
	normalize(result, result, 0, 255, CV_MINMAX);
	imwrite(save_path_ideal_high_pass_filter, result);
}

// 巴特沃斯低通滤波模板
// 阶数n=1 无振铃和负值
// 阶数n=2 轻微振铃和负值
// 阶数n=5 明显振铃和负值
// 半径D0越小,模糊越大;半径D0越大,模糊越小
Mat butterworth_low_pass_kernel(Mat &scr, float D0, int n, string save_path_butterworth_low_pass_kernel) {
	Mat butterworth_low_pass(scr.size(), CV_32FC1);
	for (int i = 0; i < scr.rows; i++) {
		for (int j = 0; j < scr.cols; j++) {
			//计算pow必须为float型
			double d = sqrt(pow((i - scr.rows / 2), 2) + pow((j - scr.cols / 2), 2));
			butterworth_low_pass.at<float>(i, j) = 1.0 / (1 + pow(d / D0, 2 * n));
		}
	}
	//imshow("butterworth_low_path_kernel", butterworth_low_pass);
	//waitKey(0);
	//destroyAllWindows();
	// 归一化成255范围保存的图片和上面显示的效果一样
	normalize(butterworth_low_pass, butterworth_low_pass, 0, 255, CV_MINMAX);
	imwrite(save_path_butterworth_low_pass_kernel, butterworth_low_pass);
	return butterworth_low_pass;
}

// 巴特沃斯低通滤波器
void butterworth_low_pass_filter(string load_path, float d0, int n, string save_path_butterworth_low_pass_kernel, string save_path_butterworth_low_pass_filter) {
	Mat pScr;
	// 1为加载图像的颜色类型,1为原图返回,0为灰度返回,注意这里得到的pScr必须为灰度图格式,否则后面会报错,原因是如果读入RGB图像,padded也是三通道的RGB格式,那么Mat planes[]就会报错
	pScr = imread(load_path, 0);
	int row = pScr.rows;
	int col = pScr.cols;
	int channel = pScr.channels();
	int M = getOptimalDFTSize(pScr.rows);
	int N = getOptimalDFTSize(pScr.cols);
	Mat padded;
	copyMakeBorder(pScr, padded, 0, M - pScr.rows, 0, N - pScr.cols, BORDER_CONSTANT, Scalar::all(0));
	//将图像转换为float型
	padded.convertTo(padded, CV_32FC1);
	Mat butterworth_kernel = butterworth_low_pass_kernel(padded, d0, n, save_path_butterworth_low_pass_kernel);//理想低通滤波器
	Mat result = freqfilt(padded, butterworth_kernel);
	//imshow("butterworth_low_pass_filter_picture", result);
	//waitKey(0);
	//destroyAllWindows();
	normalize(result, result, 0, 255, CV_MINMAX);
	imwrite(save_path_butterworth_low_pass_filter, result);
}

// 巴特沃斯高通滤波模板
Mat butterworth_high_pass_kernel(Mat &scr, float D0, int n, string save_path_butterworth_high_pass_kernel) {
	Mat butterworth_high_pass(scr.size(), CV_32FC1);
	for (int i = 0; i < scr.rows; i++) {
		for (int j = 0; j < scr.cols; j++) {
			//计算pow必须为float型
			double d = sqrt(pow((i - scr.rows / 2), 2) + pow((j - scr.cols / 2), 2));
			butterworth_high_pass.at<float>(i, j) = 1.0 / (1 + pow(D0 / d, 2 * n));
		}
	}
	//imshow("butterworth_high_path_kernel", butterworth_high_pass);
	//waitKey(0);
	//destroyAllWindows();
	// 归一化成255范围保存的图片和上面显示的效果一样
	normalize(butterworth_high_pass, butterworth_high_pass, 0, 255, CV_MINMAX);
	imwrite(save_path_butterworth_high_pass_kernel, butterworth_high_pass);
	return butterworth_high_pass;
}

// 巴特沃斯高通滤波器
void butterworth_high_pass_filter(string load_path, float d0, int n, string save_path_butterworth_high_pass_kernel, string save_path_butterworth_high_pass_filter) {
	Mat pScr;
	// 1为加载图像的颜色类型,1为原图返回,0为灰度返回,注意这里得到的pScr必须为灰度图格式,否则后面会报错,原因是如果读入RGB图像,padded也是三通道的RGB格式,那么Mat planes[]就会报错
	pScr = imread(load_path, 0);
	int row = pScr.rows;
	int col = pScr.cols;
	int channel = pScr.channels();
	int M = getOptimalDFTSize(pScr.rows);
	int N = getOptimalDFTSize(pScr.cols);
	Mat padded;
	copyMakeBorder(pScr, padded, 0, M - pScr.rows, 0, N - pScr.cols, BORDER_CONSTANT, Scalar::all(0));
	//将图像转换为float型
	padded.convertTo(padded, CV_32FC1);
	Mat butterworth_kernel = butterworth_high_pass_kernel(padded, d0, n, save_path_butterworth_high_pass_kernel);//理想低通滤波器
	Mat result = freqfilt(padded, butterworth_kernel);
	//imshow("butterworth_high_pass_filter_picture", result);
	//waitKey(0);
	//destroyAllWindows();
	normalize(result, result, 0, 255, CV_MINMAX);
	imwrite(save_path_butterworth_high_pass_filter, result);
}

int main() {
	string load_path = "C:/Users/zgcr6/Desktop/高图实验/5/save/lenna.bmp";
	string save_path_gray_image = "C:/Users/zgcr6/Desktop/高图实验/5/save/lenna_gray.bmp";
	gray_image(load_path, save_path_gray_image);
	string load_gray_path = save_path_gray_image;
	string save_path_DFT_gray = "C:/Users/zgcr6/Desktop/高图实验/5/save/lenna_DFT_gray.bmp";
	DFT(load_gray_path, save_path_DFT_gray);
	string save_path_DFT_IDFT_gray_1 = "C:/Users/zgcr6/Desktop/高图实验/5/save/lenna_DFT_IDFT_gray_1 .bmp";
	string save_path_DFT_IDFT_gray_2 = "C:/Users/zgcr6/Desktop/高图实验/5/save/lenna_DFT_IDFT_gray_2 .bmp";
	DFT_IDFT(load_gray_path, save_path_DFT_IDFT_gray_1, save_path_DFT_IDFT_gray_2);
	string save_path_ideal_low_pass_kernel = "C:/Users/zgcr6/Desktop/高图实验/5/save/lenna_ideal_low_pass_kernel .bmp";
	string save_path_ideal_low_pass_filter_gray = "C:/Users/zgcr6/Desktop/高图实验/5/save/lenna_ideal_low_pass_filter_gray .bmp";
	ideal_low_pass_filter(load_gray_path, 60.0, save_path_ideal_low_pass_kernel, save_path_ideal_low_pass_filter_gray);
	string save_path_ideal_high_pass_kernel = "C:/Users/zgcr6/Desktop/高图实验/5/save/lenna_ideal_high_pass_kernel.bmp";
	string save_path_ideal_high_pass_filter_gray = "C:/Users/zgcr6/Desktop/高图实验/5/save/lenna_ideal_high_pass_filter_gray .bmp";
	ideal_high_pass_filter(load_gray_path, 30.0, save_path_ideal_high_pass_kernel, save_path_ideal_high_pass_filter_gray);
	string save_path_butterworth_low_path_kernel = "C:/Users/zgcr6/Desktop/高图实验/5/save/lenna_butterworth_low_path_kernel_gray .bmp";
	string save_path_butterworth_low_pass_filter_gray = "C:/Users/zgcr6/Desktop/高图实验/5/save/lenna_butterworth_low_pass_filter_gray .bmp";
	butterworth_low_pass_filter(load_gray_path, 60.0, 2, save_path_butterworth_low_path_kernel,save_path_butterworth_low_pass_filter_gray);
	string save_path_butterworth_high_path_kernel = "C:/Users/zgcr6/Desktop/高图实验/5/save/lenna_butterworth_high_path_kernel_gray .bmp";
	string save_path_butterworth_high_pass_filter_gray = "C:/Users/zgcr6/Desktop/高图实验/5/save/lenna_butterworth_high_pass_filter_gray .bmp";
	butterworth_high_pass_filter(load_gray_path, 60.0, 2, save_path_butterworth_high_path_kernel, save_path_butterworth_high_pass_filter_gray);
	return 0;
}