我们通过这个示例来学习 一个以dft()为核心,对图像求傅里叶变换的过程。
程序示例
如何计算和显示傅里叶变换后的幅度图像。由于数字图像的离散性,我门也仅仅实现的是离散傅里叶变换,通过傅里叶变换得到图像中的几何结构信息。下面以输入图像的单通道灰度图为例。
新函数部分参考上讲的内容: https://editor.csdn.net/md/?articleId=124288761
1.【第一步】载入图像
以灰度图模式读取原始图像,进行是否读取成功的判断,并显示图像。
Mat srcImage = imread("1.jpg", 0);
if(!srcImage.data) {
printf("获取图像失败!\n");
return false;
}
imshow("原始图像", srcImage);
2 .【第二步】将图像扩大到合适的尺寸。
离散傅里叶变换的运行速度与图片的尺寸有很大的关系。当图像尺寸是2,3,5的整数倍时,计算速度最快。因此,为了达到快速计算的目的,通常会添凑新的边缘像素的幅方法来获取最佳的图像尺寸。函数getOptimalDFTSize()
用于返回最佳尺寸,而函数copyMakeBorder()
用于填充边缘像素。
//[2]将输入图像延扩到最佳的尺寸,边界用0补充
int m = getOptimalDFTSize(I.rows);
int n = getOptimalDFTSize(I.cols);
//将添加的像素初始化为0
Mat padded;
copyMakeBorder(I, padded, 0, m-I.rows, 0, n-I.cols, BORDER_CONSTANT, Scalar::all(0));
3 .【第三步】为傅里叶变换的结果(实部和虚部)分配存储空间。
傅里叶变换的结果是复数,对每个图像值会有实和虚部两个图像值。频域值的范围大与时域值范围,因此要将频域值存储在float格式中,所以我们此时就需要将输入图像转为浮点类型,并加一个额为的通道来存储复数部分,
//为傅里叶变换的结果分配(实部和虚部)存储空间
//将palanes数组组合合并成一个多通道的数组complexI
Mat palnes[] = {Mat_<float>(padded), Mat::zeros(padded.size(), CV_32F)};
Mat complexI;
merge(palnes, 2, complexI);
4 .【第四步】进行离散傅里叶变换
图像就地计算模式(in-place),输入输出为同一幅图像。
//进行就地离散傅里叶变换
dft(complexI, complecI);
5 .【第五步】将复数转变为幅值
离散傅里叶变换的结果是复数,对应的幅值:
M
=
R
e
(
D
F
T
(
I
)
)
2
+
I
m
(
D
F
T
(
I
)
)
2
M = \sqrt {Re(DFT(I))^2 + Im(DFT(I))^2}
M=Re(DFT(I))2+Im(DFT(I))2
//将复数转变为幅值,即log(1+ M)
split(complexI, palnes); //将多通道数组complex分离成单通道数组
//palnes[0] = Re(DFT(I)), palnes[1] = Im(DFT(I))
magnitude(palnes[0], palnes[1], palnes[0]); // palnes[0] = magnitude
Mat magnitudeImage = planes[0];
6 .【第六步】进行对数尺度缩放
离散傅里叶的幅度范围不适合在屏幕上显示。高值呈现为白点,低值为黑点,高低值的变化无法有效的分辨,为了在屏幕上显示高低变化的连续性,可用用对数尺度替换线性尺度。
M
1
=
l
o
g
(
1
+
M
)
M_1 = log(1 + M)
M1=log(1+M)
//进行对数尺度的缩放
magnitudeImage += Scalar::all(1);
log(magnitudeImage, magnitudeImage); //求自然对数
7.【第七步】剪切和重分布幅度图象限
在第二步中扩展了图像,将新添加的像素剔除,为了方便显示,可以重新分布幅度图象限的位置(将第五步中的图从中间划开得到四个子图,每张子图就是幅度图的象限,重新分布,即将四个角点重叠到图像中心)。这样原点就位移到了图像的中心。
//剪切和重新分布幅度图象限
//若有奇数行或奇数列,进行频谱裁减
magnitudeImage = magnitudeImage(Rect(0,0,magnitudeImage.cols & -2, magnitudeImage.rows & -1));
//重新排列傅里叶图像中的象限,使得原点位于图像中心
int cv = magnitudeImage.cols/2;
int cy = magnitudeImage.rows/2;
Mat q0(magnitudeImage,Rect(0,0,cx,cy)); //ROI左上角
Mat q1(magnitudeImage, Rect(cx, 0, cx,cy)); // 右下角
Mat q2(magnitudeImage, Rect(0,cy, cx, cy));
Mat q3(magnitudeImage, 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);
8 . 【第八步】归一化
重新分布后的幅度图,其幅度值仍然超过了可示范围[0,1],我们使用normalize()函数将幅度归一化到可显示的范围
//归一化,用0-1之间的浮点值将矩阵变换到可视的图像格式
//normalize(magnitudeImage, magnitudeImage, 0 , 1, CV_MINMAX);//opencv2
normalize(magnitudeImage, magnitudeImage, 0, 1, NORM_MINMAX);//opencv3
9.【第九步】显示效果图
//显示频谱图
imshow("频谱图",magnitudeImage);
分步式的综合示例如下:
#include "opencv2/opencv.hpp"
#include "opencv2/imgproc/imgproc.hpp"
#include "opencv2/highgui/highgui.hpp"
#include<iostream>
using namespace std;
using namespace cv;
int main()
{
//[1]以灰度图读取原始图像并显示
Mat srcImage = imread("../1.jpg", 0);
if(!srcImage.data) {
printf("获取图像失败\n");
return false;
}
imshow("原始图", srcImage);
//[2]将输入图像扩展到最佳尺寸,边界填0
int m= getOptimalDFTSize(srcImage.rows);
int n = getOptimalDFTSize(srcImage.cols);
//将添加的像素初始化为0
Mat padded;
copyMakeBorder(srcImage, padded, 0, m - srcImage.rows, 0, n - srcImage.cols, BORDER_CONSTANT, Scalar::all(0));
//[3]为傅里叶变换的结果分配存储空间(实部和虚部)
//将planes数组合并成一个多通道的数组complexI
Mat planes[] = {Mat_<float>(padded), Mat ::zeros(padded.size(), CV_32F)};
Mat complexI;
merge(planes, 2, complexI);
//[4]进行就地离散傅里叶变换
dft(complexI, complexI);
//[5]将复数转换为幅度值
split(complexI, planes); //将多通道数组分离成几个单通道数组
// planes[0] = Re(DFT(I));
// planes[1] = Im(DFT(I));
magnitude(planes[0], planes[1], planes[0]);
Mat magnitudeImage = planes[0];
//[6]进行对数尺度缩放
magnitudeImage += Scalar::all(1);
log(magnitudeImage, magnitudeImage);
//[7]剪切和重分幅度图象限
//若有奇数行或列,进行频谱裁剪
magnitudeImage = magnitudeImage(Rect(0, 0, magnitudeImage.cols & -2, magnitudeImage.rows & -2));
//重新排列傅里叶变换中的象限,使得原点位于图像中心
int cx = magnitudeImage.cols/2;
int cy = magnitudeImage.rows/2;
Mat q0(magnitudeImage,Rect(0,0,cx,cy)); //ROI左上角
Mat q1(magnitudeImage, Rect(cx, 0, cx,cy)); // 右下角
Mat q2(magnitudeImage, Rect(0,cy, cx, cy));
Mat q3(magnitudeImage, 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);
//[8]归一化,用0-1之间的浮点值将矩阵变换到可视的图像格式
//normalize(magnitudeImage, magnitudeImage, 0 , 1, CV_MINMAX);//opencv2
normalize(magnitudeImage, magnitudeImage, 0, 1, NORM_MINMAX);//opencv3
imshow("频谱图", magnitudeImage);
// waitKey(0);
while(char (waitKey(1)) != 'q'){}
return 0;
}
运行结果: