本文共 2958 字,大约阅读时间需要 9 分钟。
最近由于作业原因,试着用OpenCV实现频率域滤波,但是OpenCV中并没有像MATLAB中fftshift这样的中心化操作,所以我写了一个频率域滤波的函数,以后用频率域滤波的时候在主函数中调用即可。当然,水平有限,编写的代码并不优美,有问题请大家留言批评指正。
在这里我不介绍傅里叶变换,频率域滤波和高斯低通滤波器的原理,想必大家已经有了大概了解,本文关注于OpenCV中的代码实现。
废话少说,先上一张例图:
这幅图像的噪声为周期性噪声,带用阻滤波器或者陷波滤波器去噪较好,这里用高斯低通去噪,效果一般,周期性没有完全去除。
下面是频率域滤波的流程:
1、计算原图像的DFT,得到F(U,V);
2、将频谱F(U,V)的零频点移到中心位置;
3、计算滤波器函数H(U,V)与上步结果F(U,V)的乘积G(U,V);
4、将上步的IDFT,得到g(X,Y);
5、取g(X,Y)的模作为结果。
下面根据上述步骤编写代码如下:
#include#include #include using namespace std;using namespace cv;Mat gaussianlbrf(Mat scr, float sigma) { Mat gaussianBlur(scr.size(), CV_32FC1); float d0 = 2 * sigma * sigma; for (int i = 0; i < scr.rows; ++i) { for (int j = 0; j < scr.cols; ++j) { float d = powf((i - scr.rows / 2) * (i - scr.rows / 2) + (j - scr.cols / 2) * (j - scr.cols / 2)); gaussianBlur.at (i, j) = expf(-d / d0); } } return gaussianBlur;}Mat freqfilt(Mat scr, Mat blur) { Mat complexIm; split(scr, complexIm); dft(complexIm, complexIm); // 将零频点移到中心 int cx = complexIm.cols / 2; int cy = complexIm.rows / 2; // 处理实部 Mat part1_r(complexIm, Rect(0, 0, cx, cy)); Mat part2_r(complexIm, Rect(cx, 0, cx, cy)); Mat part3_r(complexIm, Rect(0, cy, cx, cy)); Mat part4_r(complexIm, 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); // 处理虚部 Mat part1_i(complexIm, Rect(0, 0, cx, cy)); Mat part2_i(complexIm, Rect(cx, 0, cx, cy)); Mat part3_i(complexIm, Rect(0, cy, cx, cy)); Mat part4_i(complexIm, 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); // 滤波 Mat blur_r, blur_i; multiply(complexIm, blur, blur_r); multiply(complexIm, blur, blur_i); Mat BLUR; merge({blur_r, blur_i}, BLUR); magnitude(complexIm, complexIm, complexIm); complexIm += Scalar::all(1); log(complexIm, complexIm); normalize(complexIm, complexIm, 1, 0, CV_MINMAX); idft(BLUR, BLUR); magnitude(BLUR, BUR, BUR); normalize(BUR, BUR, 1, 0, CV_MINMAX); return BUR;}
int main() { Mat input = imread("imageHill.jpg", CV_LOAD_IMAGE_GRAYSCALE); imshow("原图", input); int w = getOptimalDFTSize(input.cols); int h = getOptimalDFTSize(input.rows); Mat padded; copyMakeBorder(input, padded, 0, h - input.rows, 0, w - input.cols, BORDER_CONSTANT, Scalar::all(0)); padded.convertTo(padded, CV_32FC1); Mat gaussianKernel = gaussianlbrf(padded, 45); Mat out = freqfilt(padded, gaussianKernel); imshow("结果图 sigma=40", out); waitKey(0); return 0;}
注:1,结果的黑边为填充效果。
2,务必注意矩阵数据类型,做DFT,必须为浮点型,计算滤波器函数H(U,V)与DFT结果的乘积G(U,V)时,数据类型务必一致,包括通道数,单通道类型不能与多通道类型相乘,关于数据类型,可以参看。
参考:
转载地址:http://bwgfk.baihongyu.com/