c++ opencv数字图像处理：频率域滤波--同态滤波

文章目录
前言
一、同态滤波原理
1.处理原理
二、同态滤波器模板及MATLAB代码
1.同态滤波器
2.代码
前言
数字图像处理c++ opencv（VS2019 opencv4.53）持续更新

一、同态滤波原理
1.处理原理
（1）认为图像f(x,y)由两部分组成：照射分量i(x,y)，反射分量r(x,y):
f ( x , y ) = i ( x , y ) ∗ r ( x , y ) f(x,y)=i(x,y)*r(x,y)
f(x,y)=i(x,y)∗r(x,y)

（2）但上式不能直接用于对两个分量在频率域进行处理，因为乘积的傅里叶变换不等于变换的乘积，因此，定义：
z ( x , y ) = l n f ( x , y ) = l n i ( x , y ) + l n r ( x , y ) z(x,y)=lnf(x,y)=lni(x,y)+lnr(x,y)
z(x,y)=lnf(x,y)=lni(x,y)+lnr(x,y)

（3）经过上面变换后可以进行傅里叶变换，有：
Z ( u , v ) = F i ( u , v ) + F r ( u , v ) Z(u,v)=F_i(u,v)+F_r(u,v)
Z(u,v)=F 
i
​
 (u,v)+F 
r
​
 (u,v)

（4）然后可以设置一个滤波器H(u,v),对Z(u,v)进行滤波：
S ( u , v ) = H ( u , v ) ∗ Z ( u , v ) = H ( u , v ) ∗ F i ( u , v ) + H ( u , v ) ∗ F r ( u , v ) S(u,v)=H(u,v)*Z(u,v)=H(u,v)*F_i(u,v)+H(u,v)*F_r(u,v)
S(u,v)=H(u,v)∗Z(u,v)=H(u,v)∗F 
i
​
 (u,v)+H(u,v)∗F 
r
​
 (u,v)

（5）对S(u,v)使用傅里叶逆变换得到s(x,y)，则滤波后的图像g(x,y):
g ( x , y ) = e s ( x , y ) g(x,y)=e^{s(x,y)}
g(x,y)=e 
s(x,y)
 

在图像中，认为低频成分与照射分量相相联系，高频成分与反射分量相联系。同态滤波就是设置一个滤波器H(x,y)，使用不同的可控方法来影响低频分量和高频分量对图像的影响。

二、同态滤波器模板及MATLAB代码
1.同态滤波器
（1）高斯高通滤波的函数如下：
H ( u , v ) = 1 − e − D 2 ( u , v ) / 2 D 0 2 H(u,v)=1-e^{-D^2(u,v)/2D_0^2}
H(u,v)=1−e 
−D 
2
 (u,v)/2D 
0
2
​
 

高斯高通滤波器剖面图如上，使低频分量处为0，高频分量处保持不变

（2）以高斯高通滤波器为模板改造的同态滤波器

H ( u , v ) = ( γ H − γ L ) [ 1 − e − C [ D 2 ( u , v ) / 2 D 0 2 ] ] + γ L H(u,v)=(γ_H-γ_L)[1-e^{-C[D^2(u,v)/2D_0^2]}]+γ_L
H(u,v)=(γ 
H
​
 −γ 
L
​
 )[1−e 
−C[D 
2
 (u,v)/2D 
0
2
​
 ]
 ]+γ 
L
​
 

其 中 γ H 和 γ L 规 定 频 率 的 范 围 ， c 控 制 中 间 的 偏 斜 度 其中γ_H和γ_L规定频率的范围，c控制中间的偏斜度其中γ 
H
​
 和γ 
L
​
 规定频率的范围，c控制中间的偏斜度

当γ H > 1 γ_H>1γ 
H
​
 >1,γ L < 1 γ_L<1γ 
L
​
 <1时，对应剖面图为：

可以看出这个滤波器减少了低频分量的影响，增加了高频分量的影响。

（3）同态滤波步骤：

2.代码
代码如下：

#include<iostream>
#include<opencv2/opencv.hpp>
#include "MY_DFT.h"

using namespace cv;
using namespace std;

int main()
{
    Mat image, image_gray, image_output, image_transform;   //定义输入图像，灰度图像，输出图像
    image = imread("lena.png");  //读取图像；
    if (image.empty())
    {
        cout << "读取错误" << endl;
        return -1;
    }
    imshow("image", image);

    cvtColor(image, image_gray, COLOR_BGR2GRAY); //转换为灰度图
    imshow("image_gray", image_gray); //显示灰度图

    //1、进行ln变换
    Mat image_gray2(image_gray.size(), CV_32F);
    for (int i = 0; i < image_gray.rows; i++)
    {
        for (int j = 0; j < image_gray.cols; j++)
        {
            image_gray2.at<float>(i, j) = log(image_gray.at<uchar>(i, j) + 0.1);
        }
    }

    //2、傅里叶变换，image_output为可显示的频谱图，image_transform为傅里叶变换的复数结果
    My_DFT(image_gray2, image_output, image_transform);
    imshow("image_output", image_output);

    //3、高斯高通滤波
    Mat planes[] = { Mat_<float>(image_output), Mat::zeros(image_output.size(),CV_32F) };
    split(image_transform, planes);//分离通道，获取实部虚部
    Mat image_transform_real = planes[0];
    Mat image_transform_imag = planes[1];

    int core_x = image_transform_real.rows / 2;//频谱图中心坐标
    int core_y = image_transform_real.cols / 2;
    int d0 = 10;  //滤波半径
    float h;
    //参数：
    float rh = 3;
    float rl = 0.5;
    float c = 5;
    for (int i = 0; i < image_transform_real.rows; i++)
    {
        for (int j = 0; j < image_transform_real.cols; j++)
        {
            h = (rh-rl) * (1 - exp(-c * ((i - core_x) * (i - core_x) + (j - core_y) * (j - core_y)) / (d0 * d0))) + rl;
            image_transform_real.at<float>(i, j) = image_transform_real.at<float>(i, j) * h;
            image_transform_imag.at<float>(i, j) = image_transform_imag.at<float>(i, j) * h;

        }
    }
    planes[0] = image_transform_real;
    planes[1] = image_transform_imag;
    Mat image_transform_ilpf;//定义高斯高通滤波结果
    merge(planes, 2, image_transform_ilpf);

    //4、傅里叶逆变换
    Mat iDft[] = { Mat_<float>(image_output), Mat::zeros(image_output.size(),CV_32F) };
    idft(image_transform_ilpf, image_transform_ilpf);//傅立叶逆变换
    split(image_transform_ilpf, iDft);//分离通道，主要获取0通道
    magnitude(iDft[0], iDft[1], iDft[0]); //计算复数的幅值，保存在iDft[0]
    normalize(iDft[0], iDft[0], 0, 1, NORM_MINMAX);//归一化处理
    exp(iDft[0], iDft[0]);
    normalize(iDft[0], iDft[0], 0, 1, NORM_MINMAX);//归一化处理
    iDft[0].convertTo(iDft[0], CV_8U, 255 / 1.0, 0);
    imshow("idft", iDft[0]);//显示逆变换图像
    waitKey(0);  //暂停，保持图像显示，等待按键结束
    return 0;
}

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
傅里叶变换代码（.h文件）：

#pragma once
#include<iostream>
#include<opencv2/opencv.hpp>
#include<cmath>

using namespace cv;
using namespace std;

void My_DFT(Mat input_image, Mat& output_image, Mat& transform_array);

1
2
3
4
5
6
7
8
9
10
傅里叶变换代码（.cpp文件）：

#include "MY_DFT.h"

//傅里叶变换得到频谱图和复数域结果
void My_DFT(Mat input_image, Mat& output_image, Mat& transform_image)
{
    //1.扩展图像矩阵，为2，3，5的倍数时运算速度快
    int m = getOptimalDFTSize(input_image.rows);
    int n = getOptimalDFTSize(input_image.cols);
    copyMakeBorder(input_image, input_image, 0, m - input_image.rows, 0, n - input_image.cols, BORDER_CONSTANT, Scalar::all(0));

    //2.创建一个双通道矩阵planes，用来储存复数的实部与虚部
    Mat planes[] = { Mat_<float>(input_image), Mat::zeros(input_image.size(), CV_32F) };

    //3.从多个单通道数组中创建一个多通道数组:transform_image。函数Merge将几个数组合并为一个多通道阵列，即输出数组的每个元素将是输入数组元素的级联
    merge(planes, 2, transform_image);

    //4.进行傅立叶变换
    dft(transform_image, transform_image);

    //5.计算复数的幅值，保存在output_image（频谱图）
    split(transform_image, planes); // 将双通道分为两个单通道，一个表示实部，一个表示虚部
    Mat transform_image_real = planes[0];
    Mat transform_image_imag = planes[1];

    magnitude(planes[0], planes[1], output_image); //计算复数的幅值，保存在output_image（频谱图）

    //6.前面得到的频谱图数级过大，不好显示，因此转换
    output_image += Scalar(1);   // 取对数前将所有的像素都加1，防止log0
    log(output_image, output_image);   // 取对数
    normalize(output_image, output_image, 0, 1, NORM_MINMAX); //归一化

    //7.剪切和重分布幅度图像限
    output_image = output_image(Rect(0, 0, output_image.cols & -2, output_image.rows & -2));

    // 重新排列傅里叶图像中的象限，使原点位于图像中心
    int cx = output_image.cols / 2;
    int cy = output_image.rows / 2;
    Mat q0(output_image, Rect(0, 0, cx, cy));   // 左上区域
    Mat q1(output_image, Rect(cx, 0, cx, cy));  // 右上区域
    Mat q2(output_image, Rect(0, cy, cx, cy));  // 左下区域
    Mat q3(output_image, Rect(cx, cy, cx, cy)); // 右下区域

      //交换象限中心化
    Mat tmp;
    q0.copyTo(tmp); q3.copyTo(q0); tmp.copyTo(q3);//左上与右下进行交换
    q1.copyTo(tmp); q2.copyTo(q1); tmp.copyTo(q2);//右上与左下进行交换

    Mat q00(transform_image_real, Rect(0, 0, cx, cy));   // 左上区域
    Mat q01(transform_image_real, Rect(cx, 0, cx, cy));  // 右上区域
    Mat q02(transform_image_real, Rect(0, cy, cx, cy));  // 左下区域
    Mat q03(transform_image_real, Rect(cx, cy, cx, cy)); // 右下区域
    q00.copyTo(tmp); q03.copyTo(q00); tmp.copyTo(q03);//左上与右下进行交换
    q01.copyTo(tmp); q02.copyTo(q01); tmp.copyTo(q02);//右上与左下进行交换

    Mat q10(transform_image_imag, Rect(0, 0, cx, cy));   // 左上区域
    Mat q11(transform_image_imag, Rect(cx, 0, cx, cy));  // 右上区域
    Mat q12(transform_image_imag, Rect(0, cy, cx, cy));  // 左下区域
    Mat q13(transform_image_imag, Rect(cx, cy, cx, cy)); // 右下区域
    q10.copyTo(tmp); q13.copyTo(q10); tmp.copyTo(q13);//左上与右下进行交换
    q11.copyTo(tmp); q12.copyTo(q11); tmp.copyTo(q12);//右上与左下进行交换

    planes[0] = transform_image_real;
    planes[1] = transform_image_imag;
    merge(planes, 2, transform_image);//将傅里叶变换结果中心化
}

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
结果：

文章知识点与官方知识档案匹配，可进一步学习相关知识
————————————————

                            版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
                        
原文链接：https://blog.csdn.net/qq_44785013/article/details/121648185