C++ opencv 卷积 傅里叶变换
程序
废话不多说,直接上代码。
因为两个函数卷积的傅里叶变换等于两个函数傅里叶变换的点乘,所以在下面的程序中,我们就直接进行点乘,不再进行中心化处理。
//该函数只能处理灰度图像
Mat FrequencyFilter(Mat& src, int Sx, int Sy)
{
// src的宽高
int R = src.rows;
int C = src.cols;
// 卷积核的宽高
int r = 2 * Sx + 1;
int c = 2 * Sy + 1;
//将mask的深度作改变,提高精度
Mat kernel = cv::Mat::ones(r, c, CV_64FC1);
/* 1.卷积边界扩充 */
Mat src_padded;
copyMakeBorder(src, src_padded, Sx, Sx, Sy, Sy, BORDER_REPLICATE);
//将src.padded的深度作改变,提高精度
src_padded.convertTo(src_padded, CV_64FC1);
/* 2.补0以满足快速傅里叶变换的行数和列数 */
// 获取满足二维快速傅里叶变换的行数和列数,并且进行了避免缠绕错误的处理
int rows = getOptimalDFTSize(src_padded.rows + r - 1);
int cols = getOptimalDFTSize(src_padded.cols + c - 1);
//输出适合快速傅里叶变化的行数和列数
std::cout << cols << " " << rows << std::endl;
Mat src_padded_zeros, kernel_zeros;
//填充参数不一样,输出的图像边缘也会有很大的差别,在这里我们采用BORDER_REPLICATE,其会使输出的图像边界比直接补0更加平缓,大家可以改变边界填充方式,观察输出结果的边界形式
copyMakeBorder(src_padded, src_padded_zeros, 0, rows - src_padded.rows, 0, cols - src_padded.cols, BORDER_REPLICATE);
copyMakeBorder(kernel, kernel_zeros, 0, rows - r, 0, cols - c, BORDER_CONSTANT, Scalar(0, 0, 0));
/* 3.快速傅里叶变换 */
Mat fft_src, fft_kernel;
//设定输出结果为复数
dft(src_padded_zeros, fft_src, DFT_COMPLEX_OUTPUT);
dft(kernel_zeros, fft_kernel, DFT_COMPLEX_OUTPUT);
/* 4. 两个傅里叶变换点乘 */
Mat Dot_Multiplication;
mulSpectrums(fft_src, fft_kernel, Dot_Multiplication, DFT_ROWS);
/* 5.1傅里叶逆变换,只取实部 */
//个人认为会存在一些误差,但在实际图像效果上来看与取幅值基本无差别,即人眼看不出什么差别
/*
Mat dst;
dft(Dot_Multiplication, dst, DFT_INVERSE + DFT_SCALE + DFT_REAL_OUTPUT);
normalize(dst, dst, 0, 255, NORM_MINMAX);
Mat result = dst(Rect(Sy * 2, Sx * 2, C, R)).clone();
result.convertTo(result, CV_8UC1);
return result;
*/
/* 5.2傅里叶逆变换,取实部和虚部的幅值 */
Mat dst;
dft(Dot_Multiplication, dst, DFT_INVERSE + DFT_SCALE + DFT_COMPLEX_OUTPUT);
Mat plane[2];
split(dst, plane);
magnitude(plane[0], plane[1], plane[0]);
//打印一下幅度值,看一下是幅度的范围,大家可以用这段程序去查看各个阶段变换后的数据
for (int i = 0; i < 60; i++)
{
for (int j = 0; j < 60; j++)
{
std::cout << plane[0].at<double>(i, j) << std::endl;
}
}
//对幅度进行归一化,理论上来说是不需要的,欢迎大家在评论区给出自己的想法
//2021.09.08更新,需要归一化的原因在于做傅里叶变换时的mask没有进行归一化
//所以在此处把灰度值缩放到0~255
normalize(plane[0], plane[0], 0, 255, NORM_MINMAX);
//* 6.裁剪,与所输入的图像矩阵的尺寸相同
//注意坐标系,故写作以下形式
//此处的从(Sx*2,Sy*2)开始裁剪的原因尚且未知,欢迎大家在评论区给出自己的想法
Mat result = plane[0](Rect(Sy * 2, Sx * 2, C, R)).clone();
//为了显示图像,作深度变换,不作变换时,实际打印数据为doulble类型
result.convertTo(result, CV_8UC1);
return result;
}
输出
原灰度图
时间域平滑结果
频率域平滑结果
分析
假设图像大小为m × n,mask大小为Sx × Sy
----- | 时间域平滑 | 积分图 | 频率域平滑 |
---|---|---|---|
输出效果 | 彼此一致 | 彼此一致 | 彼此一致 |
输出速度 | 慢 | 快 | 中等 |
时间复杂度 | m * n * Sx * Sy | m * n | (m * n) * log(m * n) |