一、概述
图像的傅里叶变换及其两个重要的度量:幅度谱和相位谱。了解两个重要的概念:低频和高频。低频指的是图
的傅里叶变换
“
中心位置
”
附近的区域。注意,如无特殊说明,后面所提到的图像的傅里叶变换都是中心化后的。高频随着到“
中心位置
”
距离的增加而增加,即傅里叶变换中心位置的外围区域,这里的“
中心位置
”
指的是傅里叶变换所对应的幅度谱最大值的位置。
频率域滤波器在程序或者数学运算中的呈现可以理解为一个矩阵,该矩阵的宽、高和图像的傅里叶变换的宽、高是相同的,下面所涉及的常用的低通、高通、带通、带阻 等滤波的关键步骤,就是通过一定的准则构造该矩阵的。
二、步骤
步骤如下:
![](https://i-blog.csdnimg.cn/blog_migrate/f8145f46ab559e6c12b688e5f24bafe4.png)
下面通过一个简单的例子来详细解释频率域滤波的整个步骤。
第一步:输入图像矩阵
I
。假设为:
第二步:图像矩阵的每一个像素值乘以
(-1)r+c
得到矩阵
I′
,
I′ =I.*
(
-1
)
r+c
,其
中
r
和
c
代表当前像素值在矩阵中的位置索引。
![](https://i-blog.csdnimg.cn/blog_migrate/39433e456fdec278bf3523714e3fcb5d.png)
![](https://i-blog.csdnimg.cn/blog_migrate/24a0a878cf3bf69c15708772f6c6edfd.png)
第三步:因为图像矩阵的宽和高均为
7
,为了利用傅里叶变换的快速算法,对
I′
补
0
,使用命令getOptimalDFTSize
(
7
)得到一个不小于
7
且可以分解为
2p ×3q ×5r
的最小整数,计算结果为8
。所以在矩阵
I′
的右侧和下侧各补一行
0
,记为
f
:
![](https://i-blog.csdnimg.cn/blog_migrate/fb359a647c85def6770a261015a43069.png)
第四步:利用傅里叶变换的快速算法得到复数矩阵F。
OpenCV是将复数矩阵按照双通道存储的,即第一通道存储的是复数矩阵的实部,第二通道存储的是复数矩阵的虚部。
第五步:构建频率域滤波器
Filter
。频率域滤波器本质上是一个和第四步得到的快速傅里叶变换矩阵F 具有相同行数、列数的复数矩阵,一般情况下为实数矩阵,这里假设是 一个全是1的矩阵:
![](https://i-blog.csdnimg.cn/blog_migrate/eab3ec9f9ee7c2e588772a9caea520ac.png)
本章提到的频率域滤波,如低通滤波、高通滤波、自定义滤波等,其关键步骤就是 通过一定的标准构造该矩阵以完成图像在频率域上的滤波。
第六步:将第四步得到的快速傅里叶变换矩阵
F
和第五步得到的频率域滤波器
Filter 的对应位置相乘(矩阵的点乘)。当然,如果滤波器是一个实数矩阵,那么在代码实现 中,将傅里叶变换的实部和虚部分别与频率域滤波器进行点乘即可,即Ffilter =F.*Filter,因为这里构造的滤波器是一个全是1
的矩阵,所以
F filter =F
。
第七步:对第六步得到的点乘矩阵
Ffilter
进行傅里叶逆变换,得到复数矩阵
F′
。
第八步:
取复数矩阵F′的实部
。
第九步:与第二步类似,将第八步得到的矩阵乘以(
-1
)
r+c
。
第十步:因为在快速傅里叶变换的步骤中进行了补
0
操作,所以第九步得到的实部矩 阵的尺寸有可能比原图大,所以要进行裁剪,取该实部矩阵的左上角,尺寸和原图相同。裁剪得到的结果,即为频率域滤波的结果。在该示例中,因为滤波器是一个全是1
的
矩阵,相当于对原图没有做任何处理,即最后滤波的结果和原图是一样的。
频率域滤波算法均是按照上述十个步骤完成的,接下来详细介绍常用滤波器的构建
方法、代码实现及其效果。
三、低通滤波
针对图像的傅里叶变换,低频信息表示图像中灰度值缓慢变化的区域;而高频信息则正好相反,表示灰度值变化迅速的部分,如边缘。低通滤波,顾名思义,保留傅里叶变换的低频信息;或者削弱傅里叶变换的高频信息;而高通滤波则正好相反,保留傅里叶变换的高频信息,移除或者削弱傅里叶变换的低频信息。
三种常用的低通滤波器:
H、
W 分别代表图像快速傅里叶变换的高、宽,
傅里叶谱的最大值(中心点)的位置在maxR,maxC
,
radius
代表截断频率,
D
(
r
,
c
)代表到中心位置的距离
1、理想低通滤波:ilpFilter=[ilpFilter(r,c)]H*w,
![](https://i-blog.csdnimg.cn/blog_migrate/c54d85f504f9d284beeb3d8b8b907965.png)
2、巴特沃斯低通滤波器
第二种是巴特沃斯低通滤波器,记
blpFilter=[blpFilter
(
r
,
c
)
]H×W
![](https://i-blog.csdnimg.cn/blog_migrate/61e2cf1b7e5c9915981ad02f750d92b4.png)
3、高斯低通滤波器
记glpFilter=[glpFilter(r,c)]H×W;
作用:
滤波器越靠近中心点位置的值越接近于1,越远离中心位置的值就越小于1,与傅里叶变换相乘后,相当于保留了低频信息,消弱或者移除了高频信息
四、低通滤波的代码实现
// 快速傅里叶变换
void fft2Image(InputArray image, OutputArray F)
{
//的到Mat 的类型
Mat Image = image.getMat();
int rows = Image.rows;
int cols = Image.cols;
// 获得最佳两列
int r = getOptimalDFTSize(rows);
int c = getOptimalDFTSize(cols);
// 在左侧和下面进10
Mat f;
copyMakeBorder(Image, f, 0, r - rows, 0, c - cols, BORDER_CONSTANT,Scalar::all(0));
// 两个通道,第一个通道用存储实部 第二个通道是存放虚部的
dft(f, F, DFT_COMPLEX_OUTPUT);
//傅里叶逆变换
//dft(f, F, DFT_INVERSE+DFT_REAL_OUTPUT+DFT_SCALE);//只返回实部
}
// 幅度谱
void amplitudeSpectrum(InputArray _srcfft, OutputArray _dstSpectrum)
{
if (_srcfft.channels() != 2) return ;
// 如果是两个通道则开始分离通道
vector<Mat> FFT2Channels;
split(_srcfft, FFT2Channels);
// 计算傅里叶变换幅度值 sqrt(pow(r,2)+pow(L,2))
magnitude(FFT2Channels[0], FFT2Channels[1], _dstSpectrum);
}
// 傅里叶谱的灰度级显示
Mat graySpectrum(Mat s)
{
Mat dst;
log(s + 1, dst);
// 归一化
normalize(dst,dst,0,1,NORM_MINMAX);
// 为了显示 灰度级
dst.convertTo(dst,CV_8UC1,255,0);
return dst;
}
// 计算傅里叶变换的相位角==phase 算子
Mat phaseS(Mat sfft)
{
// 相位谱
Mat phase;
phase.create(sfft.size(), CV_64FC1);
// 分离通道
vector<Mat> fft2Channels;
split(sfft,fft2Channels);
// 开始计算相位谱
for (int r = 0; r < phase.rows;r++) {
for (int c = 0; c < phase.cols; c++) {
// 实部 虚部
double real = fft2Channels[0].at<double>(r,c);
double img = fft2Channels[1].at<double>(r,c);
// atan2 的返回值就是【0,180】,【-180,0】
phase.at<double>(r, c) = atan2(img,real);
}
}
return phase;
}
// -----------------------------------------------------------
#if 1 //低通滤波
enum typeFilter
{
ILPFILTER,
BLPFILTER,
GLPFILTER,
};
Mat createFilter(Size size, Point center,float radius, int type, int n = 2)
{
Mat ilpFilter = Mat::zeros(size,CV_32FC1); // 定义一个size大小的滤波器
int rows = size.height;
int cols = size.width;
if (radius <= 0)return ilpFilter;
// 构建理想的低通滤波器
if (type== ILPFILTER)
{
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(r - center.y)), 2);
if (norm2 < radius)
{
ilpFilter.at<float>(r, c) = 1;
}
else
{
ilpFilter.at<float>(r, c) = 0;
}
}
}
}
//
else if(type == BLPFILTER)
{
for (int r = 0; r < rows; r++)
{
for (int c = 0; c < cols; c++)
{
/* float m = sqrt(pow(r - center.y, 2.0) + pow(c - center.x, 2.0));
float n = m / radius;
float o = pow(n,2*n);
float s = 1.0 + o;
float i = 1.0 / s;*/
ilpFilter.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 * n))));
}
}
}
// 高斯滤波
else if (type == GLPFILTER)
{
for (int r = 0; r < rows; r++)
{
for (int c = 0; c < cols; c++)
{
ilpFilter.at<float>(r, c) = float(exp(-(pow(c-center.x,2.0)+pow(r-center.y,2.0))/(2*pow(radius,2.0))));
}
}
}
return ilpFilter;
}
Mat src;// 输入的图像
Mat F;// 图像的快速傅里叶变换
Point maxLocation;//傅里叶谱的最大值的坐标
int radius = 20; // 截断频率
const int max_Radius = 100;// 最大截断频率
Mat ilpFilter;//低通滤波器
int filterType=0;//
int max_FilterType = 2;
Mat F_ilpFilter;// 低通傅里叶变换
Mat FlpSpetrum;// 低通傅里叶变换的傅里叶谱的灰度级
Mat result;// 滤波后的图像
string lpFilterSpectrum = "低通傅里叶谱";//
void callback_lpFilter(int ,void *);
int main()
{
src = imread("C:\\Users\\19473\\Desktop\\opencv_images\\611.png",CV_LOAD_IMAGE_GRAYSCALE);
if (!src.data)
{
cout << "原始图读取失败" << endl;
return -1;
}
imshow("原图",src);
// 数据转换 将数据转换为double
Mat f1;
src.convertTo(f1,CV_32FC1,1.0,0.0);
// 2 每个数 x pow(-1,c+r);
for (int r = 0; r < f1.rows; r++)
{
for (int c = 0; c < f1.cols; c++)
{
// 判断奇偶
if (r+c%2)
{
f1.at<float>(r, c) *= -1;
}
}
}
// 开始傅里叶变换
fft2Image(f1,F);// F 就是傅里叶变换之后的
Mat as;
//得到傅里叶谱
amplitudeSpectrum(F,as);
//傅里叶谱的灰度级显示
Mat s = graySpectrum(as);
imshow("原傅里叶谱的灰度级显示",s);
// 找到傅里叶谱的最大值的坐标
minMaxLoc(s,NULL,NULL,NULL,&maxLocation);
//------------------------------------ 低通滤波----------------------------------------------
namedWindow(lpFilterSpectrum,WINDOW_AUTOSIZE);
createTrackbar("低通滤波:", lpFilterSpectrum,&filterType, max_FilterType, callback_lpFilter);
createTrackbar("半径:", lpFilterSpectrum, &radius, max_Radius, callback_lpFilter);
callback_lpFilter(0,0);
waitKey(0);
return 0;
}
// 回调函数
void callback_lpFilter(int, void*)
{
// ----------------- 构建低通滤波器----------------------------------
ilpFilter = createFilter(F.size(),maxLocation,radius, filterType,2);
// -------------------- 低通滤波器和图像的傅里叶变换开始点乘=====================================
F_ilpFilter.create(F.size(),F.type());
for (int r = 0; r < F_ilpFilter.rows; r++)
{
for (int c = 0; c < F_ilpFilter.cols; c++)
{
// 分别取出当前位置的快速傅里叶变换和理想低通滤波的值
Vec2f f_rc = F.at<Vec2f>(r,c);
float lpFilter_rc = ilpFilter.at<float>(r,c);
// 低通滤波器和图像的快速傅里叶变换的对应位置相乘
F_ilpFilter.at<Vec2f>(r, c) = f_rc * lpFilter_rc;
}
}
// 低通傅里叶变换的傅里叶谱
amplitudeSpectrum(F_ilpFilter, FlpSpetrum);
FlpSpetrum = graySpectrum(FlpSpetrum);
imshow(lpFilterSpectrum, FlpSpetrum);
// 对傅里叶变换逆变换 并且只要实部
dft(F_ilpFilter,result,DFT_SCALE+DFT_INVERSE+DFT_REAL_OUTPUT);
// 乘以 (-1)的(r+c)次方
for (int r = 0; r < result.rows; r++)
{
for (int c = 0; c < result.cols; c++)
{
if ((c + r) % 2) result.at<float>(r, c) *= -1;
}
}
// 很重要的一步 将结果转换成 CV_8u
result.convertTo(result, CV_8UC1, 1.0, 0);
result = result(Rect(0,0,src.cols,src.rows)).clone();
imshow("经过低通滤波后的图像",result);
}
#endif
五、高通滤波
maxR
,
maxC
代表图像傅里叶谱的最大值的位置,D(r,c):
![](https://i-blog.csdnimg.cn/blog_migrate/7d66e5bd488caea6553d8d53a0a21eb8.png)
1、 理想高通滤波
![](https://i-blog.csdnimg.cn/blog_migrate/dbd2f0a6db732581f62be542253fc1c4.png)
2. 巴特沃斯高通滤波器
3. 高斯高通滤波器
六、带通 带阻滤波器
带通滤波是指只保留某一范围区域的频率带。下面介绍三种常用的带通滤波器。假 设BW
代表带宽,
D0
代表带宽的径向中心,其他符号与低通、高通滤波相同。
1. 理想带通滤波器和理想带阻滤波器
![](https://i-blog.csdnimg.cn/blog_migrate/ebd0757f54d885a46fa214e68b7c6a4f.png)
![](https://i-blog.csdnimg.cn/blog_migrate/d7858bb9341c760601cab7b0de9e73a3.png)
2. 巴特沃斯带通滤波器和巴特沃斯带阻滤波器
3. 高斯带通滤波器和高斯带阻滤波器
带阻滤波器和这个上述的恰好相反
七、自定义滤波器
#if 1 //自定义滤波器
Mat image;// 输入的图像
Mat ImageFFT;// 图像的快速傅里叶变换
Point maxLocation;//傅里叶谱的最大值的坐标
Mat ilpFilter;//低通滤波器
Mat F_ImageSpetrum;// 傅里叶变换的傅里叶谱
Mat result;// 滤波后的图像
string windowName = "傅里叶幅度谱的灰度级";//
bool drawing_box = false;// 鼠标事件
Point drawPoint;
Rect rectFilter;
bool gotRectFilter = false;
void mouseRectHandler(int event,int x,int y,int,void* )
{
switch (event)
{
case CV_EVENT_LBUTTONDOWN:// 鼠标按下
drawing_box = true;
// 记录起点
drawPoint = Point(x, y);
break;
case CV_EVENT_MOUSEMOVE:
if (drawing_box)
{
// 将鼠标移动到右下角
if (x> drawPoint.x&&y>= drawPoint.y)
{
rectFilter.x = drawPoint.x;
rectFilter.y = drawPoint.y;
rectFilter.width = x - drawPoint.x;
rectFilter.height= y - drawPoint.y;
}
// 将鼠标移动到有右上角
if (x >= drawPoint.x && y <= drawPoint.y)
{
rectFilter.x = drawPoint.x;
rectFilter.y = y;
rectFilter.width = x - drawPoint.x;
rectFilter.height = drawPoint.y- y ;
}
// 将鼠标移动到有左上角
if (x >= drawPoint.x && y <= drawPoint.y)
{
rectFilter.x = x;
rectFilter.y = y;
rectFilter.width = drawPoint.x - x;
rectFilter.height = drawPoint.y - y;
}
// 将鼠标移动到有左下角
if (x >= drawPoint.x && y <= drawPoint.y)
{
rectFilter.x = x;
rectFilter.y = drawPoint.y;
rectFilter.width = drawPoint.x-x;
rectFilter.height = y- drawPoint.y ;
}
}
break;
case CV_EVENT_LBUTTONUP:
drawing_box = false;
gotRectFilter = true;
break;
default:
break;
}
}
int main()
{
image = imread("C:\\Users\\19473\\Desktop\\opencv_images\\612.png",CV_LOAD_IMAGE_GRAYSCALE);
if (!image.data)
{
cout << "原始图读取失败" << endl;
return -1;
}
imshow("原图", image);
// 数据转换 将数据转换为double
Mat f1;
image.convertTo(f1,CV_32FC1,1.0,0.0);
// 2 每个数 x pow(-1,c+r);
for (int r = 0; r < f1.rows; r++)
{
for (int c = 0; c < f1.cols; c++)
{
// 判断奇偶
if (r+c%2)
{
f1.at<float>(r, c) *= -1;
}
}
}
// 开始傅里叶变换
fft2Image(f1,ImageFFT);// F 就是傅里叶变换之后的
Mat amplSpec;
//得到傅里叶谱
amplitudeSpectrum(ImageFFT, amplSpec);
//傅里叶谱的灰度级显示
Mat spectrum = graySpectrum(amplSpec);
// 找到傅里叶谱的最大值的坐标
minMaxLoc(amplSpec,NULL,NULL,NULL,&maxLocation);
//------------------------------------ 自定义 滤波----------------------------------------------
namedWindow(windowName,WINDOW_AUTOSIZE);
setMouseCallback(windowName, mouseRectHandler,NULL);
for (;;)
{
spectrum(rectFilter).setTo(0);
// 自定义 滤波器和傅里叶变换点乘
ImageFFT(rectFilter).setTo(Scalar::all(0));
imshow(windowName, spectrum);
// 按esc 见退出编辑
if (waitKey(10)==27)
{
break;
}
}
Mat result;
// 对傅里叶变换逆变换 并且只要实部
dft(ImageFFT, result, DFT_SCALE + DFT_INVERSE + DFT_REAL_OUTPUT);
// 乘以 (-1)的(r+c)次方
for (int r = 0; r < result.rows; r++)
{
for (int c = 0; c < result.cols; c++)
{
if ((c + r) % 2) result.at<float>(r, c) *= -1;
}
}
// 很重要的一步 将结果转换成 CV_8u
result.convertTo(result, CV_8UC1, 1.0, 0);
result = result(Rect(0, 0, image.cols, image.rows)).clone();
imshow("经过自定义滤波后的图像", result);
waitKey(0);
return 0;
}
// 回调函数
#endif