频率域滤波

最近有个朋友创业，找我实现个算法，用到离散傅里叶变换。实现后，正好把以前相关的笔记整理完。

频域滤波

空间域和频率域为我们提供了不同的视角。在空域中，数字图像f(x,y)为一个定义在二维空间中的矩形区域上的离散函数；换一个角度，，如果将f(x,y)视为幅值变化的二维信号，则可以通过某些变换手段（如傅里叶变换、离散余弦变换、沃尔什变换和小波变换等）在频域下对它进行分析。离散傅里叶变化情况下，相位谱决定图像结构，幅度谱反映了图像整体各个方向的频率分量的相对强度。

在多数情况下，频率域滤波和空间域滤波可以视为同一个图像增强问题的殊途同归的两种解决方式。傅里叶变换提供了一种变换到频率域的手段。关于傅里叶变换，这里就不累赘介绍了。最早我是从本科接触傅里叶变换，先是《信号与系统》，后来是《数字信号处理》。

为了得到合适的空间域滤波器，我们可以首先设计频域滤波器H(u, v)，然后根据卷积定理，将H(u, v)反变换到空间域得到空间域中滤波使用的卷积模板h(x,y)。

频率滤波的基本步骤

1. 计算原图像f(x,y)的DFT，得到F(u,v)。
2. 将频谱F(u,v)的零频点移到中心位置。
3. 计算滤波器函数H(u,v)与F(u,v)的乘积G(u, v)。
4. 将频谱G(u, v)的零频点移回频谱图的左上角位置。
5. 计算上步计算结果的IDFT，得到g(x,y)。
6. 取g(x,y)的实部作为最终滤波后的结果图像。

示例演示

我们按照上面的基本步骤实现一个高斯低通滤波来去噪。高斯低通滤波如下

代码实现如下

void MainWindow::gaussianBlur(const cv::Mat &image)
{
    Mat input;
    image.convertTo(input, CV_32F);
    std::vector<cv::Mat> channels;
    split(image, channels);  //分离图像的RGB通道, OpenCV:BGR
    for(auto& splitedimage : channels){
        //第一步
        int w = getOptimalDFTSize(splitedimage.cols);
        int h = getOptimalDFTSize(splitedimage.rows);
        Mat padded;
        copyMakeBorder(splitedimage, padded, 0, h - splitedimage.rows, 0, w - splitedimage.cols, BORDER_CONSTANT, Scalar::all(0));
        //imshow("padded", padded);
        Mat plane[] = {cv::Mat_<float>(padded), Mat::zeros(padded.size(),CV_32F)};
        Mat compleximage;
        merge(plane, 2, compleximage);
        dft(compleximage, compleximage);
        //第二步
        split(compleximage, plane);//以下的操作是移动图像  (零频移到中心)
        int cx = plane[0].cols / 2;
        int cy = plane[0].rows / 2;
        Mat part1r(plane[0],Rect(0,0,cx,cy));  //元素坐标表示为(cx,cy)
        Mat part2r(plane[0],Rect(cx,0,cx,cy));
        Mat part3r(plane[0],Rect(0,cy,cx,cy));
        Mat part4r(plane[0],Rect(cx,cy,cx,cy));
        Mat temp;
        part1r.copyTo(temp);  //左上与右下交换位置(实部)
        part4r.copyTo(part1r);
        temp.copyTo(part4r);
        part2r.copyTo(temp);  //右上与左下交换位置(实部)
        part3r.copyTo(part2r);
        temp.copyTo(part3r);
        Mat part1i(plane[1],Rect(0,0,cx,cy));  //元素坐标(cx,cy)
        Mat part2i(plane[1],Rect(cx,0,cx,cy));
        Mat part3i(plane[1],Rect(0,cy,cx,cy));
        Mat part4i(plane[1],Rect(cx,cy,cx,cy));
        part1i.copyTo(temp);  //左上与右下交换位置(虚部)
        part4i.copyTo(part1i);
        temp.copyTo(part4i);
        part2i.copyTo(temp);  //右上与左下交换位置(虚部)
        part3i.copyTo(part2i);
        temp.copyTo(part3i);
        //第三步
        Mat gaussianblur(padded.size(), CV_32FC1);
        float sigma = 45.0;
        float d0 = 2* sigma* sigma;
        for(int i = 0;i < padded.rows; i++ ) {
            for(int j = 0; j < padded.cols; j++ ) {
              float d = pow(float(i - padded.rows / 2), 2) + pow(float(j - padded.cols / 2 ), 2);//分子,计算pow必须为float型
              gaussianblur.at<float>(i,j) = expf(-d / d0);//expf为以e为底求幂（必须为float型）
            }
        }
        imshow("GausianBlur", gaussianblur);
        multiply(plane[0], gaussianblur, plane[0]); //滤波（实部与滤波器模板对应元素相乘）
        multiply(plane[1], gaussianblur, plane[1]);//滤波（虚部与滤波器模板对应元素相乘）
        //第四步
        //以下的操作是移动图像  (零频移到左上角)
        Mat rpart1r(plane[0],Rect(0,0,cx,cy));  //元素坐标表示为(cx,cy)
        Mat rpart2r(plane[0],Rect(cx,0,cx,cy));
        Mat rpart3r(plane[0],Rect(0,cy,cx,cy));
        Mat rpart4r(plane[0],Rect(cx,cy,cx,cy));
        rpart1r.copyTo(temp);  //左上与右下交换位置(实部)
        rpart4r.copyTo(rpart1r);
        temp.copyTo(rpart4r);
        rpart2r.copyTo(temp);  //右上与左下交换位置(实部)
        rpart3r.copyTo(rpart2r);
        temp.copyTo(rpart3r);
        Mat rparti(plane[1],Rect(0,0,cx,cy));  //元素坐标(cx,cy)
        Mat rpart2i(plane[1],Rect(cx,0,cx,cy));
        Mat rpart3i(plane[1],Rect(0,cy,cx,cy));
        Mat rpart4i(plane[1],Rect(cx,cy,cx,cy));
        rparti.copyTo(temp);  //左上与右下交换位置(虚部)
        rpart4i.copyTo(rparti);
        temp.copyTo(rpart4i);
        rpart2i.copyTo(temp);  //右上与左下交换位置(虚部)
        rpart3i.copyTo(rpart2i);
        temp.copyTo(rpart3i);
        Mat blur_result;
        merge(plane, 2, blur_result);//实部与虚部合并
        //第五步
        idft(blur_result, blur_result);	//idft结果也为复数
        //第六步
        split(blur_result, plane);//分离通道，主要获取通道
        magnitude(plane[0],plane[1], plane[0]);  //求幅值(模)
        normalize(plane[0],plane[0], 1,0, CV_MINMAX);  //归一化便于显示
        plane[0].copyTo(splitedimage);
    }
    Mat reuslut;
    merge(channels, reuslut);
    imshow("Result", reuslut);
}

运行结果