opencv 傅里叶变换（低通滤波 + 高通滤波）

1、傅里叶变换

时域分析：以时间作为参照物，世间万物都是随着时间变化而变化，并且不会停止
频域分析：认为世间万物都是静止的，永恒不变的

通过以下制作饮料的过程可以很好的理解傅里叶变换。

1、从时域分析：就是六点零一放了1块冰糖，3颗红豆，2颗绿豆，4块西红柿，1杯纯净水，六点零二放了1块冰糖。。。。随着时间的变化一直在变化

2、从频域角度分析：不在是以时间为参照物了，而是这个事情的频率，1分钟放1块冰糖，2分钟放3粒红豆，3分钟放2粒绿豆，4分钟放4块西红柿，5分钟放1杯水。

下面这两个图都是描述同一个事情，可以更明显看出，两者的区别。

在两个角度去看周期函数的变化

任何连续的周期信号，可以由一组适当的正𫠊曲线组成

相为：三个开始起点不一致的余𫠊函数，组成了这个曲线

2、通过numpy实现

通过将原图进行傅里叶变换，得到频域图像，获得高频和低频，对高频和低频进行操作之后，进行逆变换回原图像达到对图像进行特色操作：图像增强、图像去噪、边缘检测、特征提取、压缩、加密等
1、低通滤波器：只保留低频信息，去掉高频信息，会去掉边缘特征信息，会让图像变模糊
2、高频滤波器：只保留高频信息，去掉低频信息，会增强图像的边缘和特征信息，但是会失去一些细节信息

def test9():img = cv2.imread("1.jpg", 0)# 执行傅里叶变换，转化为频域f = np.fft.fft2(img)# 低频在左上角，为了方便，将其移到中心位置（带负数的数组）fshift = np.fft.fftshift(f)# 通过将其转换为（0-255）中result = 20 * np.log(np.abs(fshift))# 原图显示# 创建窗口一行两列，第一咧plt.subplot(121)plt.imshow(img, cmap='gray')plt.title('orginal')# 不用坐标系plt.axis('off')# 傅里叶变换之后的图plt.subplot(122)plt.imshow(result, cmap='gray')plt.title('result')plt.axis('off')plt.show()

从频域转换会原图像（因为没有做任何操作，所以输出图像不会发生改变）

def test10():img = cv2.imread("1.jpg", 0)# 执行傅里叶变换，转化为频域f = np.fft.fft2(img)# 低频在左上角，为了方便，将其移到中心位置（带负数的数组）fshift = np.fft.fftshift(f)# 低频谱从中心移动到左上角（相当于又移回去）ishift = np.fft.ifftshift(fshift)# 从频域转换回原图像iimg = np.fft.ifft2(ishift)iimg = np.abs(iimg)# 原图显示# 创建窗口一行两列，第一咧plt.subplot(121)plt.imshow(img, cmap='gray')plt.title('orginal')# 不用坐标系plt.axis('off')# 傅里叶变换之后的图plt.subplot(122)plt.imshow(iimg, cmap='gray')plt.title('result')plt.axis('off')plt.show()

3、高通滤波器

前面已经说过，高通滤波器是去除所有低频信息，达到图像增强的效果但是会失去一些细节信息。

def test11():img = cv2.imread("1.jpg", 0)# 执行傅里叶变换，转化为频域f = np.fft.fft2(img)# 低频在左上角，为了方便，将其移到中心位置（带负数的数组）fshift = np.fft.fftshift(f)# 高宽 创建高通滤波器rows, cols = img.shapecrows, ccols = int(rows / 2), int(cols / 2)# 前面我们将低频信息移动到图像中间来了# 这里将低频信息全部设置为0，达到去掉低频信息目的fshift[crows - 30:crows + 30, ccols - 30:ccols + 30] = 0ifshift = np.fft.ifftshift(fshift)# 将频谱逆变换到图像iimg = np.fft.ifft2(ifshift)iimg = np.abs(iimg)plt.subplot(121)plt.imshow(img, cmap='gray')plt.title('orginal')# 不用坐标系plt.axis('off')# 傅里叶变换之后的图plt.subplot(122)plt.imshow(iimg, cmap='gray')plt.title('result')plt.axis('off')plt.show()

5、通过opencv实现傅里叶变换

使用opencv中的函数实现

def test12():img = cv2.imread("1.jpg", 0)# dft返回的是两个通道的频域，0是频域实部分，1是频域图像虚部分# DFT_COMPLEX_OUTPUT输出复数dft = cv2.dft(np.float32(img), flags=cv2.DFT_COMPLEX_OUTPUT)dftshift = np.fft.fftshift(dft)idftshift = np.fft.ifftshift(dftshift)# 傅里叶逆变换iimg = cv2.idft(idftshift)# magnitude函数频域图像的幅度谱# result = 20 * np.log(cv2.magnitude(dftshift[:, :, 0], dftshift[:, :, 1]))result = cv2.magnitude(iimg[:, :, 0], iimg[:, :, 1])plt.subplot(121)plt.imshow(img, cmap='gray')plt.title('orginal')# 不用坐标系plt.axis('off')# 傅里叶变换之后的图plt.subplot(122)plt.imshow(result, cmap='gray')plt.title('result')plt.axis('off')plt.show()

6、低通滤波器

创建一个mask，将中心位置设为1，其他位置设为0，然后和频谱图像相乘之后就只保留了低频信息了

def test13():img = cv2.imread("1.jpg", 0)# # 执行傅里叶变换dft = cv2.dft(np.float32(img), flags=cv2.DFT_COMPLEX_OUTPUT)ifshift = np.fft.fftshift(dft)# 创建低通滤波器rows, cols = img.shapemask = np.zeros((rows, cols, 2), np.int8)crows, ccols = int(rows / 2), int(cols / 2)mask[crows - 30:crows + 30, ccols - 30:ccols + 30] = 1fshift = ifshift * maskishift = np.fft.ifftshift(fshift)io = cv2.idft(ishift)result = cv2.magnitude(io[:, :, 0], io[:, :, 1])plt.subplot(121)plt.imshow(img, cmap='gray')plt.title('orginal')# 不用坐标系plt.axis('off')# 傅里叶变换之后的图plt.subplot(122)plt.imshow(result, cmap='gray')plt.title('result')plt.axis('off')plt.show()

7、C++实现傅里叶变换

傅里叶变换：cv::dft()
执行一维或二维浮点数组的正向或反向离散傅里叶变换。

#include <opencv2/core.hpp>
函数说明：void cv::dft( InputArray src, OutputArray dst, int flags = 0, int nonzeroRows = 0 );
输入参数：（1）src					输入数组。可以是实数也可以是复数。（2）dst					输出数组。其大小和类型取决于flags。（3）flags = 0				转换标志。cv::DFT_INVERSE 				执行1D或2D逆变换，而不是默认的正转换。cv::DFT_SCALE 					缩放比例标识符，输出的结果会以1/N进行缩放。N=数组元素的数量cv::DFT_ROWS 					对输入矩阵的每一行执行正变换或逆变换。能够同时处理多个向量，并减少3D和高维变换的开销。cv::DFT_COMPLEX_OUTPUT 			一维或二维实数数组正变换。cv::DFT_REAL_OUTPUT 			一维或二维复数数组逆变换。cv::DFT_COMPLEX_INPUT 			指定输入为实数输入。输入必须有2个通道。cv::DCT_INVERSE					执行逆1D或2D转换，而不是默认的正向转换。cv::DCT_ROWS 					对输入矩阵的每一行执行正变换或逆变换。能够同时处理多个向量，并减少3D和高维变换的开销。（4）nonzeroRows = 0		默认值为0。若设为非零值，dft函数会将该值作为非零行的有效区间长度，只对非零行进行处理，提高计算效率。

傅里叶反变换：cv::idft()
计算一维或二维阵列的离散傅里叶反变换。
默认情况下，dft和idft都不会缩放结果。因此，您应该显式地将DFT_SCALE传递给dft或idft中的一个，以使这些变换相互逆。

#include <opencv2/core.hpp>
函数说明：void cv::idft( InputArray src, OutputArray dst, int flags = 0, int nonzeroRows = 0 );
输入参数：（1）src					输入数组。可以是实数也可以是复数。（2）dst					输出数组。其大小和类型取决于flags。（3）flags = 0				转换标志。cv::DFT_INVERSE 				执行1D或2D逆变换，而不是默认的正转换。cv::DFT_SCALE 					缩放比例标识符，输出的结果会以1/N进行缩放。N=数组元素的数量cv::DFT_ROWS 					对输入矩阵的每一行执行正变换或逆变换。能够同时处理多个向量，并减少3D和高维变换的开销。cv::DFT_COMPLEX_OUTPUT 			一维或二维实数数组正变换。cv::DFT_REAL_OUTPUT 			一维或二维复数数组逆变换。cv::DFT_COMPLEX_INPUT 			指定输入为实数输入。输入必须有2个通道。cv::DCT_INVERSE					执行逆1D或2D转换，而不是默认的正向转换。cv::DCT_ROWS 					对输入矩阵的每一行执行正变换或逆变换。能够同时处理多个向量，并减少3D和高维变换的开销。（4）nonzeroRows = 0		默认值为0。若设为非零值，dft函数会将该值作为非零行的有效区间长度，只对非零行进行处理，提高计算效率。

计算相位谱：cv::phase()
计算由x和y的对应元素组成的每个2D矢量的旋转角度。计算公式：angle(I)=atan2(y(I), x(I))
角度估计精度约为0.3度。当x(I)=y(I)=0时，对应角(I)设为0

#include <opencv2/core.hpp>
函数说明：void cv::phase( InputArray x, InputArray y, OutputArray angle, bool angleInDegrees = false );
输入参数：（1）x					输入二维矢量的x坐标的浮点数组。（2）y					输入二维矢量的y坐标数组。它必须和x有相同的大小和类型。（3）angle				角度。输出与x大小和类型相同的数组。（4）angleInDegrees = false		当为true时，该函数以度数计算角度，否则以弧度计算。

计算幅度谱：cv::magnitude()
计算由x和y数组的相应元素组成的2D向量的大小。计算公式：dst(I) = sqrt(x(I)^2 + y(I)^2)

#include <opencv2/core.hpp>
函数说明：void cv::magnitude( InputArray x, InputArray y, OutputArray magnitude );
输入参数：（1）x					输入二维矢量的x坐标的浮点数组。（2）y					输入二维矢量的y坐标数组。它必须和x有相同的大小和类型。（3）magnitude 			幅值。输出与x大小和类型相同的数组。

计算x和y的坐标：cv::polarToCart()
通过二维矢量的大小和角度来计算x和y的坐标。估计坐标的相对精度约为1e-6。
计算公式：x(I) = magnitude(I) * cos(angle(I))
计算公式：y(I) = magnitude(I) * sin(angle(I))

#include <opencv2/core.hpp>
函数说明：void cv::polarToCart( InputArray magnitude, InputArray angle, OutputArray x, OutputArray y, bool angleInDegrees = false );
输入参数：（1）magnitude			输入二维矢量（大小）的浮点数组。若为空矩阵，则假设所有的大小都是=1;若不为空，则必须具有与angle相同的大小和类型。（2）angle				输入二维矢量（角度）的浮点数组。（3）x					二维矢量的x坐标输出数组。它有相同的尺寸和类型的角度。（4）y					二维矢量的y坐标输出数组。它有相同的尺寸和类型的角度。（5）angleInDegrees = false  	当为true时，输入角以度数表示，否则以弧度表示。

获取最适合傅里叶正变换的宽 / 高：cv::getOptimalDFTSize()
DFT（傅里叶正变换）性能不是向量大小的单调函数。因此，当您计算两个数组的卷积或执行数组的频谱分析时，通常有必要在输入数据中填充零，以获得比原始数组转换速度快得多的更大的数组。大小为2的幂(2,4,8,16,32，…)的数组处理速度最快。但是，数组的大小是2、3和5的乘积(例如，300 = 55322)的处理效率也很高。

#include <opencv2/core.hpp>
函数说明：int cv::getOptimalDFTSize( int vecsize);
输入参数：		vecsize					给定向量。如果vecsize太大(非常接近INT_MAX)，则返回一个负数。
返回值：			N						返回大于或等于vecsize的最小数 N。

// 傅里叶变换
// 高频：变化剧烈的灰度分量，列如边界
// 低频：变化缓慢的灰度分量，例如一边大海
// 高频滤波器：只保留高频，会使得图像细节增强
// 低频滤波器：只保留低频，会使图像变模糊
void test_f()
{Mat img = imread("H:\\数据集\\已标注\\images\\datas\\all_datas\\1616003650999.jpg", IMREAD_GRAYSCALE);img.convertTo(img, CV_32F);// 数据准备// 离散傅里叶变换的运行速度与图像的大小有很大的关系，当图像的尺寸使2，3，5的整数倍时，计算速度最快// 为了达到快速计算的目的，经常通过添加新的边缘像素的方法获取最佳图像尺寸int w1 = getOptimalDFTSize(img.rows);int h1 = getOptimalDFTSize(img.cols);Mat padding;copyMakeBorder(img, padding, 0, w1-img.rows, 0, h1-img.cols, BORDER_CONSTANT,Scalar::all(0));// 执行傅里叶变换// 为傅立叶变换的结果分配存储空间// 将plannes数组组合成一个多通道的数组，两个同搭配，分别保存实部和虚部// 傅里叶变换的结果使复数，这就是说对于每个图像原像素值，会有两个图像值// 此外，频域值范围远远超过图象值范围，因此至少将频域储存在float中// 所以我们将输入图像转换成浮点型，并且多加一个额外通道来存储复数部分Mat planes[] = { Mat_<float>(padding),Mat::zeros(padding.size(),CV_32F) };Mat complexI;// planes[0]是实部，planes[1]是虚部merge(planes, 2,complexI);cout << complexI.size() << endl;cout << planes->size() << endl;dft(complexI, complexI, DFT_SCALE | DFT_COMPLEX_OUTPUT);split(complexI, planes);// 计算幅度普和相位谱Mat ph,mag,idft;phase(planes[0], planes[1], ph);magnitude(planes[0], planes[1], mag);// 重新排列傅里叶图像中的象限，使得原点位于图像中心int cx = mag.cols/2;int cy = mag.rows/2;Mat q1 = mag(Rect(0, 0, cx, cy));Mat q2 = mag(Rect(cx, 0, cx, cy));Mat q3 = mag(Rect(0, cy, cx, cy));Mat q4 = mag(Rect(cx, cy, cx, cy));// 变换左上角和右下角象限Mat tmp;q1.copyTo(tmp);q4.copyTo(q1);tmp.copyTo(q4);// 变换右上角和左下角q2.copyTo(tmp);q3.copyTo(q2);tmp.copyTo(q3);imshow("mag", mag);// 对频域进行滤波操作// 对频域进行滤波操作 高频滤波器,也可以直接这句代替// mag(Rect(cx - 30, cy - 30, 30 * 2, 30 * 2)) = 0;// 如果像素距离图像中心的水平偏差（abs(i - mag.cols / 2)）或垂直偏差（abs(j - mag.rows / 2)）大于图像尺寸的十分之一（mag.cols / 10 或 mag.rows / 10），则满足条件。for (int i = 0; i < mag.cols; i++) {for (int j = 0; j < mag.rows; j++) {if (abs(i - mag.cols / 2) > mag.cols / 10 || abs(j - mag.rows / 2) > mag.rows / 10)// 0表示低通、1表示高通mag.at<float>(j, i) = 1;}}imshow("mag1", mag);//3.4、变换左上角和右下角象限q1.copyTo(tmp);q4.copyTo(q1);tmp.copyTo(q4);//3.5、变换右上角和左下角象限q2.copyTo(tmp);q3.copyTo(q2);tmp.copyTo(q3);//（4）傅里叶逆变换polarToCart(mag, ph, planes[0], planes[1]);//由幅度谱mag和相位谱ph恢复实部planes[0]和虚部planes[1]merge(planes, 2, idft);dft(idft, idft, cv::DFT_INVERSE | cv::DFT_REAL_OUTPUT);//  & -2 操作用于将宽度和高度向下取整为偶数img = idft(cv::Rect(0, 0, img.cols & -2, img.rows & -2));img.convertTo(img, CV_8U);imshow("3", img);waitKey(0);
}