OpenCV C++（十一）----频率域滤波

肉松饼饼

IP属地: 广东

0.1 2020.03.08 17:42 字数 1440

11.1概述和原理

低频指的是图像的傅里叶变换“中心位置”附近的区域。低频信息表示图像中灰度值缓慢变化的区域。
高频随着到“中心位置”距离的增加而增加，即傅里叶变换中心位置的外围区域，这里的“中心位置”指的是傅里叶变换所对应的幅度谱最大值的位置。高频信息则正好相反，表示灰度值变化迅速的部分，如边缘。

频率域滤波器在程序或者数学运算中的呈现可以理解为一个矩阵。下面所涉及的常用的低通、高通、带通、带阻等滤波的关键步骤，就是通过一定的准则构造该矩阵的。

image.png

算法步骤：

第一步：输入图像矩阵I。
第二步：图像矩阵的每一个像素值乘以（-1） r+c得到矩阵I′， I′ =I.*（-1） r+c ，其中r和c代表当前像素值在矩阵中的位置索引。
第三步：因为图像矩阵的宽和高均为7，为了利用傅里叶变换的快速算法，对I′补0，使用命令getOptimalDFTSize（7） 得到一个不小于7且可以分解为2p ×3q ×5r的最小整数，计算结果为8。所以在矩阵I′的右侧和下侧各补一行0，记为f
第四步：利用傅里叶变换的快速算法得到复数矩阵F。
第五步：构建频率域滤波器Filter。频率域滤波器本质上是一个和第四步得到的快速傅里叶变换矩阵F 具有相同行数、列数的复数矩阵，一般情况下为实数矩阵
第六步：将第四步得到的快速傅里叶变换矩阵F 和第五步得到的频率域滤波器Filter 的对应位置相乘（矩阵的点乘）。当然，如果滤波器是一个实数矩阵，那么在代码实现中，将傅里叶变换的实部和虚部分别与频率域滤波器进行点乘即可，即F filter =F.*Filter，因为这里构造的滤波器是一个全是1的矩阵，所以F filter =F。
第七步：对第六步得到的点乘矩阵F filter进行傅里叶逆变换，得到复数矩阵F′。
第八步：取复数矩阵F′的实部。
第九步：与第二步类似，将第八步得到的矩阵乘以（-1） r+c。
第十步：进行裁剪，取该实部矩阵的左上角，尺寸和原图相同。裁剪得到的结果，即为频率域滤波的结果。

11.2、低通滤波器和高通滤波器

11.2.1、低通滤波器

1、理想低通滤波器

image.png

2、巴特沃斯低通滤波器

image.png

3、高斯低通滤波器

image.png

//创建三种常见的低通滤波器
enum FILTER_TYPE {
    ILP_FILTER = 0,
    BLP_FILTER = 1,
    GLP_FILTER = 2
};

Mat creatLPFilter(Size size, Point center, float radius, FILTER_TYPE type, int n)
{
    Mat lpFilter = Mat::zeros(size, CV_32FC1);
    int rows = size.height;
    int cols = size.width;
    if (radius <= 0)
    {
        return lpFilter;
    }
    //构建理想低通滤波器
    if (type == ILP_FILTER)
    {
        for (int r = 0; r < rows; r++)
        {
            for (int c = 0; c < cols; c++)
            {
                float norm2 = pow(abs(float(r - center.y)), 2) + pow(abs(float(c - center.x)), 2);
                if (sqrt(norm2) < radius)
                {
                    lpFilter.at<float>(r,c) = 1;
                }
                else
                {
                    lpFilter.at<float>(r, c) = 0;
                }
            }
        }
    }
    //构建巴特沃斯低通滤波器
    if (type == BLP_FILTER)
    {
        for (int r = 0; r < rows; r++)
        {
            for (int c = 0; c < cols; c++)
            {
                lpFilter.at<float>(r, c) = float(1.0 / (1.0 + pow(sqrt(pow(r - center.y, 2.0) + pow(c - center.x, 2.0)) / radius, 2.0*n)));
            }
        }
    }
    //构建高斯低通滤波器
    if (type == GLP_FILTER)
    {
        for (int r = 0; r < rows; r++)
        {
            for (int c = 0; c < cols; c++)
            {
                lpFilter.at<float>(r, c) = float(exp(-(pow(r - center.y, 2.0) + pow(c - center.x, 2.0)) / (2 * pow(radius, 2.0))));
            }
        }
    }
    return lpFilter;
}

Mat I;//输入的图像矩阵
Mat F;//图像的快速傅里叶变换
Point maxLoc;//傅里叶谱的最大值的坐标
int radius = 20;//截断频率
const int Max_RADIUS = 100;//设置的最大的截断频率
Mat lpFilter;//低通滤波器
int lpType = 0;//低通滤波器的类型
const int MAX_LPTYPE = 2;
Mat F_lpFilter;//低通傅里叶变换
Mat flpSpectrum;//低通傅里叶变换的傅里叶谱的灰度级
Mat result;//低通滤波后的效果
string lpFilterspectrum = "低通傅里叶谱";//显示窗口的名称
void callback_lpFilter(int, void*);
int main()
{
    //step1：读入图像矩阵
    Mat I = imread("Koala.jpg", CV_LOAD_IMAGE_GRAYSCALE);
    if (!I.data)
    {
        return -1;
    }
    //数据类型转换为浮点型
    Mat fI;
    I.convertTo(fI, CV_32FC1, 1.0, 0.0);
    //step2:每一个数乘于（-1）^（r+c)
    for (int r = 0; r < fI.rows; r++)
    {
        for (int c = 0; c < fI.cols; c++)
        {
            if ((r + c) % 2)
            {
                fI.at<float>(r, c) *= -1;
            }
        }
    }
    //step3:补零；step4：快速傅里叶变换
    fft2Image(fI, F);
    //傅里叶谱
    Mat amplSpec;
    amplitudeSpectrum(F, amplSpec);
    //傅里叶谱的灰度级显示
    Mat spectrum = graySpectrum(amplSpec);
    //找到傅里叶谱的最大值的坐标
    minMaxLoc(spectrum, NULL, NULL, NULL, &maxLoc);
    //低通滤波
    namedWindow(lpFilterspectrum, WINDOW_AUTOSIZE);
    createTrackbar("低通类型", lpFilterspectrum, &lpType, MAX_LPTYPE, callback_lpFilter);
    createTrackbar("半径", lpFilterspectrum, &radius, Max_RADIUS, callback_lpFilter);
    callback_lpFilter(0, 0);
    waitKey(0);
    return 0;
}


void callback_lpFilter(int, void*)
{
    //step5：构建低通滤波器
    lpFilter = creatLPFilter(F.size(), maxLoc, radius, lpType, 2);
    //step6:低通滤波器和图像的快速傅里叶变换点乘
    F_lpFilter.create(F.size(), F.type());
    for (int r = 0; r < F_lpFilter.rows; r++)
    {
        for (int c = 0; c < F_lpFilter.cols; c++)
        {
            //分别取出当前位置的快速傅里叶变换和理想低通滤波器的值
            Vec2f F_rc = F.at<Vec2f>(r, c);
            float lpFilter_rc = lpFilter.at<float>(r, c);
            //低通滤波器和图像的快速傅里叶变换的对应位置相乘
            F_lpFilter.at<Vec2f>(r, c) = F_rc * lpFilter_rc;
        }
    }
    //低通傅里叶变换的傅里叶谱
    amplitudeSpectrum(F_lpFilter, flpSpectrum);
    //低通傅里叶谱的灰度级显示
    flpSpectrum = graySpectrum(flpSpectrum);
    //step7、8：对低通滤波器变换执行傅里叶逆变换，并只取实部
    dft(F_lpFilter, result, DFT_SCALE + DFT_INVERSE + DFT_REAL_OUTPUT);
    //step9：同乘以(-1)^(x+y)
    for (int r = 0; r < result.rows; r++)
    {
        for (int c = 0; c < result.cols; c++)
        {
            if ((r + c) % 2)
            {
                result.at<float>(r, c)  *= -1;
            }
        }
    }
    //注意将结果转换为CV_8U类型
    result.convertTo(result, CV_8UC1, 1.0, 0);
    //step10：截取左上部分，其大小和输入图像的大小相同
    result = result(Rect(0, 0, I.cols, I.rows)).clone();
}

11.2.2、高通滤波器

1、理想高通滤波器

image.png

2、巴特沃斯高通滤波器

image.png

3、高斯高通滤波器

image.png

显然，高通滤波器和低通滤波器满足这样的关系：

hpFilter=1-lpFilter

即1减去通过createLPFilter得到的矩阵就可以得到高通滤波器。

11.3、带通和带阻滤波器

带通滤波是指只保留某一范围区域的频率带。

11.3.1、带通滤波器

1、理想带通滤波器

image.png

2、巴特沃斯带通滤波器

image.png

3、高斯带通滤波器

image.png

11.3.2、带阻滤波器

与带通滤波相反，带阻滤波是指撤销或者消弱指定范围区域的频率带。

1、理想带阻滤波器

image.png

2、巴特沃斯带阻滤波器

image.png

3、高斯带阻滤波器

image.png

显然， 1减去带通滤波就可以得到相应的带阻滤波。

11.4、自定义滤波器

有一种通过交互式的方式，构建自定义滤波，便于消除指定的频率。自定义滤波器通常用于消除结构化噪声或者目标。

11.5、同态滤波

同态滤波是一种减少低频增加高频，从而减少光照变化并锐化边缘或者细节的滤波方法。

对于一幅图像f(x,y)，可以表示为照射分量i(x,y)和反射分量r(x,y)的乘积。
其中0<i(x,y)<∞，0<r(x,y)<1。

i(x,y)描述景物的照明，变化缓慢，处于低频成分。

r(x,y)描述景物的细节，变化较快，处于高频成分。

因为该性质是乘性的，所以不能直接使用傅里叶变换对i(x,y)和r(x,y)进行控制，因此可以先对f(x,y)取对数，分离i(x,y)和r(x,y)。令z(x,y) = ln f(x,y) = ln i(x,y) + ln r(x,y)。由于f(x,y)的取值范围为[0, L-1]，为了避免出现ln(0)的情况，故采用ln ( f(x,y) + 1 ) 来计算。