一、概要
图像经过傅里叶变换后，将图像在空域中的信息映射至频域空间中。图像的频域空间包含幅度谱以及相位谱，其中幅度谱反映的是图像的灰度信息，相位谱反应的是图像的位置信息，如轮廓。
本博文将基于傅里叶分析理论演示利用幅度谱与相位谱重构图像的过程，并验证幅度谱与相位谱在图像重构过程中的作用，最后在本文末给出全部代码。
二、图像的傅里叶变换及其可视化
首先将图像进行傅里叶变换得到图像的频谱并使其可视化。在这里应当注意傅里叶变换对是复数，因此需要构建双通道的Mat类存储复数的实部与虚部。
代码的主函数部分：
int main()
{Mat src = imread("T.jpg");Mat gray_src;cvtColor(src, gray_src, CV_BGR2GRAY);imshow("灰度化输入图像", gray_src);int m = getOptimalDFTSize(gray_src.rows);int n = getOptimalDFTSize(gray_src.cols);Mat padded;copyMakeBorder(gray_src, padded, 0, m - gray_src.rows, 0, n - gray_src.cols, BORDER_CONSTANT, Scalar::all(0));padded.convertTo(padded, CV_32FC1);//中心化for (int row = 0; row < padded.rows; row++){for (int col = 0; col < padded.cols; 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");return 0;
}
对图像求得傅里叶变换后，将图像的频谱可视化。该部分的子函数如下所示：
//根据幅度谱的定义计算傅里叶变换的幅度谱
Mat Magnitude(Mat Re, Mat Im)
{Mat rst = Mat::zeros(Re.rows, 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::zeros(Re.rows, 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::zeros(angle.rows, angle.cols, CV_32FC2);//定义一个双通道的Mat类分别存储相位谱的实部与虚部for (int row = 0; row < Phase_Spectrum.rows; row++){for (int col = 0; col < Phase_Spectrum.cols; 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::zeros(amplitude.rows, amplitude.cols, CV_32FC2);//定义一个双通道的Mat类分别存储幅度谱谱的实部与虚部for (int row = 0; row < Amplitude_Spectrum.rows; row++){for (int col = 0; col < Amplitude_Spectrum.cols; 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::zeros(angle.rows, angle.cols, CV_32FC2);//定义一个双通道的Mat类分别存储相位谱的实部与虚部for (int row = 0; row < Amplitude_Phase_Spectrum.rows; row++){for (int col = 0; col < Amplitude_Phase_Spectrum.cols; 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(gray_src.rows);int n = getOptimalDFTSize(gray_src.cols);Mat padded;copyMakeBorder(gray_src, padded, 0, m - gray_src.rows, 0, n - gray_src.cols, BORDER_CONSTANT, Scalar::all(0));padded.convertTo(padded, CV_32FC1);//中心化for (int row = 0; row < padded.rows; row++){for (int col = 0; col < padded.cols; 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");return 0;
}
Mat Magnitude(Mat Re, Mat Im)
{Mat rst = Mat::zeros(Re.rows, 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::zeros(Re.rows, 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::zeros(angle.rows, angle.cols, CV_32FC2);//定义一个双通道的Mat类分别存储相位谱的实部与虚部for (int row = 0; row < Phase_Spectrum.rows; row++){for (int col = 0; col < Phase_Spectrum.cols; 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]);
}
void Amplitude_Spectrum_Reconstruction(Mat Re, Mat Im)
{Mat amplitude = Magnitude(Re, Im);Mat Amplitude_Spectrum = Mat::zeros(amplitude.rows, amplitude.cols, CV_32FC2);//定义一个双通道的Mat类分别存储幅度谱谱的实部与虚部for (int row = 0; row < Amplitude_Spectrum.rows; row++){for (int col = 0; col < Amplitude_Spectrum.cols; 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]);
}
void Amplitude_Phase_Spectrum_Reconstruction(Mat Re, Mat Im)
{Mat angle = Phase(Re, Im);Mat amplitude = Magnitude(Re, Im);Mat Amplitude_Phase_Spectrum = Mat::zeros(angle.rows, angle.cols, CV_32FC2);//定义一个双通道的Mat类分别存储相位谱的实部与虚部for (int row = 0; row < Amplitude_Phase_Spectrum.rows; row++){for (int col = 0; col < Amplitude_Phase_Spectrum.cols; 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]);
}
