利⽤幅度谱和相位谱重构图像
⼀、概要
图像经过傅⾥叶变换后，将图像在空域中的信息映射⾄频域空间中。图像的频域空间包含幅度谱以及相位谱，其中幅度谱反映的是图像的灰度信息，相位谱反应的是图像的位置信息，如轮廓。
本博⽂将基于傅⾥叶分析理论演⽰利⽤幅度谱与相位谱重构图像的过程，并验证幅度谱与相位谱在图像重构过程中的作⽤，最后在本⽂末给出全部代码。
⼆、图像的傅⾥叶变换及其可视化
⾸先将图像进⾏傅⾥叶变换得到图像的频谱并使其可视化。在这⾥应当注意傅⾥叶变换对是复数，因此需要构建双通道的Mat类存储复数的实部与虚部。
代码的主函数部分：
int main()
{
Mat src =imread("T.jpg");
Mat gray_src;
cvtColor(src, gray_src, CV_BGR2GRAY);
imshow("灰度化输⼊图像", gray_src);
int m =getOptimalDFTSize(ws);
int n =getOptimalDFTSize(ls);
Mat padded;
copyMakeBorder(gray_src, padded,0, m - ws,0, n - ls, BORDER_CONSTANT, Scalar::all(0));
//中⼼化
for(int row =0; row < ws; row++)
{
for(int col =0; col < ls; col++)
{
padded.at<float>(row, col)*=pow(-1, row + col);
}
}
Mat planes[]={ Mat_<float>(padded),Mat::zeros(padded.size(),CV_32F)};
Mat complexI;
merge(planes,2, complexI);
dft(complexI, complexI);
split(complexI, planes);
/
/显⽰幅度谱和相位谱
Show_Spectrum(planes[0], planes[1]);
//利⽤相位谱重构图像
Phase_Spectrum_Reconstruction(planes[0], planes[1]);
//利⽤幅度谱重构图像
Amplitude_Spectrum_Reconstruction(planes[0], planes[1]);
//利⽤幅度谱和相位谱重构图像
Amplitude_Phase_Spectrum_Reconstruction(planes[0], planes[1]);
waitKey(0);
system("pause");
return0;
}
对图像求得傅⾥叶变换后，将图像的频谱可视化。该部分的⼦函数如下所⽰：
//根据幅度谱的定义计算傅⾥叶变换的幅度谱
Mat Magnitude(Mat Re, Mat Im)
{
Mat rst = Mat::ws, Re.cols, CV_32F);
for(int row =0; row < Re.rows; row++)
{
for(int col =0; col < Re.cols; col++)
{
rst.at<float>(row, col)=sqrt(pow(Re.at<float>(row, col),2)+pow(Im.at<float>(row, col),2));
}
}
return rst;
}
//根据相位谱的定义计算傅⾥叶变换的相位谱
Mat Phase(Mat Re, Mat Im)
{
Mat rst = Mat::ws, Re.cols, CV_32F);
for(int row =0; row < Re.rows; row++)
{
for(int col =0; col < Re.cols; col++)
{
rst.at<float>(row, col)=atan2(Im.at<float>(row, col), Re.at<float>(row, col));
}
}
return rst;
}
//显⽰图像的相位谱和幅度谱
void Show_Spectrum(Mat Re, Mat Im)
{
Mat amplitude =Magnitude(Re, Im);
Mat angle =Phase(Re, Im);
amplitude += Scalar::all(1);
log(amplitude, amplitude);
normalize(amplitude, amplitude,0,1, NORM_MINMAX);
normalize(angle, angle,0,1, NORM_MINMAX);
imshow("图像傅⾥叶变换的幅度谱", amplitude);
imshow("图像傅⾥叶变换的相位谱", angle);
}
本例使⽤泰勒斯威夫特的图⽚进⾏演⽰，输⼊图像如下：
将上图灰度化后得到单通道的灰度图像，对该灰度图像进⾏傅⾥叶变换后得到的幅度谱如下图所⽰：
相位谱如下图所⽰：
三、利⽤相位谱重构图像
经过上⾯的傅⾥叶变换容易得到图像的相位谱与幅度谱，此时图像的傅⾥叶谱为|F(u,v)|*exp(jφ(u,v))。令|F(u,v)|=1，得到图像的相位谱exp(jφ(u,v))。最后利⽤欧拉公式得到： exp(jφ(u,v))=cos(φ(u,v))+jsin(φ(u,v))。
代码实现时创建⼀个双通道的Mat类存储相位谱的实部与虚部，然后对相位谱进⾏傅⾥叶反变换即可。
该部分的⼦函数如下所⽰：
void Phase_Spectrum_Reconstruction(Mat Re, Mat Im)
{
Mat angle =Phase(Re, Im);
Mat Phase_Spectrum = Mat::ws, ls, CV_32FC2);//定义⼀个双通道的Mat类分
别存储相位谱的实部与虚部
for(int row =0; row < ws; row++)
{
for(int col =0; col < ls; col++)
{
Phase_Spectrum.at<Vec2f>(row, col)[0]=cos(angle.at<float>(row, col));//实部
Phase_Spectrum.at<Vec2f>(row, col)[1]=sin(angle.at<float>(row, col));//虚部
}
}
idft(Phase_Spectrum, Phase_Spectrum);
Mat temp[]={ Mat_<float>(angle),Mat::zeros(angle.size(),CV_32F)};
split(Phase_Spectrum, temp);
magnitude(temp[0], temp[1], temp[0]);
normalize(temp[0], temp[0],0,1, CV_MINMAX);
imshow("相位谱重构图像", temp[0]);
}
利⽤相位谱重构的图像如下图所⽰：
利⽤相位谱很好的将图像的轮廓信息重构了出来，由此容易验证图像的相位谱存储的是图像的位置信息。但是由于缺乏幅度谱，因此重构后的图像缺少像素值上的变化。
四、利⽤幅度谱重构图像
令图像的傅⾥叶谱|F(u,v)|*exp(jφ(u,v))中的相位谱φ(u,v)等于0，得到幅度谱|F(u,v)|，然后对幅度谱进⾏傅⾥叶变换重构图像。
该部分的⼦函数如下所⽰：
void Amplitude_Spectrum_Reconstruction(Mat Re, Mat Im)
{
Mat amplitude =Magnitude(Re, Im);
Mat Amplitude_Spectrum = Mat::ws, ls, CV_32FC2);//定义⼀个双通道的Mat类分别存储幅度谱谱的实部与虚部
for(int row =0; row < ws; row++)
{
for(int col =0; col < ls; col++)
{
Amplitude_Spectrum.at<Vec2f>(row, col)[0]= amplitude.at<float>(row, col)*pow(-1, row + col);//实部
Amplitude_Spectrum.at<Vec2f>(row, col)[1]=0;//虚部
}
}
idft(Amplitude_Spectrum, Amplitude_Spectrum);
Mat temp[]={ Mat_<float>(amplitude),Mat::zeros(amplitude.size(),CV_32F)};
split(Amplitude_Spectrum, temp);
magnitude(temp[0], temp[1], temp[0]);
normalize(temp[0], temp[0],0,1, CV_MINMAX);
imshow("幅度谱重构图像", temp[0]);
}
利⽤幅度谱重构图像得到的结果如下图所⽰：
由幅度谱⽆法重构原来的图像，这也验证了幅度谱只包含图像的灰度信息，对图像内容起决定性作⽤的是图像的相位谱。
五、利⽤幅度谱和相位谱重构图像
利⽤欧拉公式得到：|F(u,v)|*exp(jφ(u,v))=|F(u,v)|*cos(φ(u,v))+j|F(u,v)|*sin(φ(u,v))。
代码实现时创建⼀个双通道的Mat类存储上述傅⾥叶谱的实部与虚部，然后对其进⾏傅⾥叶反变换即可。
该部分的⼦函数如下所⽰：
void Amplitude_Phase_Spectrum_Reconstruction(Mat Re, Mat Im)
{
Mat angle =Phase(Re, Im);
Mat amplitude =Magnitude(Re, Im);
Mat Amplitude_Phase_Spectrum = Mat::ws, ls, CV_32FC2);//定义⼀个双通道
的Mat类分别存储相位谱的实部与虚部for(int row =0; row < Amplitude_ws; row++)
{
for(int col =0; col < Amplitude_ls; col++)
{
Amplitude_Phase_Spectrum.at<Vec2f>(row, col)[0]= amplitude.at<float>(row, col)*cos(angle.at<float>(row, col));//实部
Amplitude_Phase_Spectrum.at<Vec2f>(row, col)[1]= amplitude.at<float>(row, col)*sin(angle.at<float>(row, col));//虚部
}
}
idft(Amplitude_Phase_Spectrum, Amplitude_Phase_Spectrum);
Mat temp[]={ Mat_<float>(angle),Mat::zeros(angle.size(),CV_32F)};
split(Amplitude_Phase_Spectrum, temp);
magnitude(temp[0], temp[1], temp[0]);
normalize(temp[0], temp[0],0,1, CV_MINMAX);
imshow("利⽤幅度谱和相位谱重构图像", temp[0]);
}
利⽤图像的幅度谱和相位谱重构的图像如下图所⽰：
很显然利⽤图像的幅度谱与相位谱完美的重构了图像。
六、完整代码
#include<iostream>
#include<opencv2/opencv.hpp>
#include<math.h>
#include<complex.h>
#include<complex>
using namespace std;
using namespace cv;
Mat Magnitude(Mat Re, Mat Im);
Mat Phase(Mat Re, Mat Im);
void Show_Spectrum(Mat Re, Mat Im);
void Phase_Spectrum_Reconstruction(Mat Re, Mat Im);
void Amplitude_Spectrum_Reconstruction(Mat Re, Mat Im);
void Amplitude_Phase_Spectrum_Reconstruction(Mat Re, Mat Im);
int main()
{
Mat src =imread("T.jpg");
Mat gray_src;
cvtColor(src, gray_src, CV_BGR2GRAY);
imshow("灰度化输⼊图像", gray_src);
int m =getOptimalDFTSize(ws);
int n =getOptimalDFTSize(ls);
Mat padded;
copyMakeBorder(gray_src, padded,0, m - ws,0, n - ls, BORDER_CONSTANT, Scalar::all(0));
//中⼼化
matlab傅里叶变换的幅度谱和相位谱for(int row =0; row < ws; row++)
{
for(int col =0; col < ls; col++)
{
padded.at<float>(row, col)*=pow(-1, row + col);
}
}
Mat planes[]={ Mat_<float>(padded),Mat::zeros(padded.size(),CV_32F)};
Mat complexI;
merge(planes,2, complexI);
dft(complexI, complexI);
split(complexI, planes);
split(complexI, planes);
//显⽰幅度谱和相位谱
Show_Spectrum(planes[0], planes[1]);
//利⽤相位谱重构图像
Phase_Spectrum_Reconstruction(planes[0], planes[1]);
//利⽤幅度谱重构图像
Amplitude_Spectrum_Reconstruction(planes[0], planes[1]);
//利⽤幅度谱和相位谱重构图像
Amplitude_Phase_Spectrum_Reconstruction(planes[0], planes[1]);
waitKey(0);
system("pause");
return0;
}
Mat Magnitude(Mat Re, Mat Im)
{
Mat rst = Mat::ws, Re.cols, CV_32F);
for(int row =0; row < Re.rows; row++)
{
for(int col =0; col < Re.cols; col++)
{
rst.at<float>(row, col)=sqrt(pow(Re.at<float>(row, col),2)+pow(Im.at<float>(row, col),2));
}
}
return rst;
}
Mat Phase(Mat Re, Mat Im)
{
Mat rst = Mat::ws, Re.cols, CV_32F);
for(int row =0; row < Re.rows; row++)
{
for(int col =0; col < Re.cols; col++)
{
rst.at<float>(row, col)=atan2(Im.at<float>(row, col), Re.at<float>(row, col));
}
}
return rst;
}
void Show_Spectrum(Mat Re, Mat Im)
{
Mat amplitude =Magnitude(Re, Im);
Mat angle =Phase(Re, Im);
amplitude += Scalar::all(1);
log(amplitude, amplitude);
normalize(amplitude, amplitude,0,1, NORM_MINMAX);
normalize(angle, angle,0,1, NORM_MINMAX);
imshow("图像傅⾥叶变换的幅度谱", amplitude);
imshow("图像傅⾥叶变换的相位谱", angle);
}
void Phase_Spectrum_Reconstruction(Mat Re, Mat Im)
{
Mat angle =Phase(Re, Im);
Mat Phase_Spectrum = Mat::ws, ls, CV_32FC2);//定义⼀个双通道的Mat类分别存储相位谱的实部与虚部for(int row =0; row < ws; row++)
{
for(int col =0; col < ls; col++)
{
Phase_Spectrum.at<Vec2f>(row, col)[0]=cos(angle.at<float>(row, col));//实部
Phase_Spectrum.at<Vec2f>(row, col)[1]=sin(angle.at<float>(row, col));//虚部
}
}
idft(Phase_Spectrum, Phase_Spectrum);
Mat temp[]={ Mat_<float>(angle),Mat::zeros(angle.size(),CV_32F)};
split(Phase_Spectrum, temp);
magnitude(temp[0], temp[1], temp[0]);
normalize(temp[0], temp[0],0,1, CV_MINMAX);