AI 图像傅里叶变换

傅立叶原理表明：

任何连续测量的时序或信号，都可以表示为不同频率的正弦波信号的无限叠加。而根据该原理创立的傅立叶变换算法利用直接测量到的原始信号，以累加方式来计算该信号中不同正弦波信号的频率、振幅和相位。

反变换从本质上说也是一种累加处理，这样就可以将单独改变的正弦波信号转换成一个信号。

傅立叶变换将原来难以处理的时域信号转换成了易于分析的频域信号（信号的频谱），最后还可以利用傅立叶反变换将这些频域信号转换成时域信号。

图像傅立叶变换的物理意义

图像的频率是表征图像中灰度变化剧烈程度的指标，是灰度在平面空间上的梯度。

如：大面积的沙漠在图像中是一片灰度变化缓慢的区域，对应的频率值很低；而对于地表属性变换剧烈的边缘区域在图像中是一片灰度变化剧烈的区域，对应的频率值较高。

傅立叶变换在实际中有非常明显的物理意义，设f是一个能量有限的模拟信号，则其傅立叶变换就表示f的谱。从纯粹的数学意义上看，傅立叶变换是将一个函数转换为一系列周期函数来处理的。

从物理效果看，傅立叶变换是将图像从空间域转换到频率域，其逆变换是将图像从频率域转换到空间域。换句话说，傅立叶变换的物理意义是将图像的灰度分布函数变换为图像的频率分布函数，傅立叶逆变换是将图像的频率分布函数变换为灰度分布函数。

 

傅立叶变换以前，图像（未压缩的位图）是由对在连续空间（现实空间）上的采样得到一系列点的集合，我们习惯用一个二维矩阵表示空间上各点，则图像可由z=f(x,y)来表示。由于空间是三维的，图像是二维的，因此空间中物体在另一个维度上的关系就由梯度来表示，这样我们可以通过观察图像得知物体在三维空间中的对应关系。为什么要提梯度？因为实际上对图像进行二维傅立叶变换得到频谱图，就是图像梯度的分布图，当然频谱图上的各点与图像上各点并不存在一一对应的关系，即使在不移频的情况下也是没有。傅立叶频谱图上我们看到的明暗不一的亮点，实际上图像上某一点与邻域点差异的强弱，即梯度的大小，也即该点的频率的大小（可以这么理解，图像中的低频部分指低梯度的点，高频部分相反）。一般来讲，梯度大则该点的亮度强，否则该点亮度弱。这样通过观察傅立叶变换后的频谱图，也叫功率图，我们首先就可以看出，图像的能量分布，如果频谱图中暗的点数更多，那么实际图像是比较柔和的（因为各点与邻域差异都不大，梯度相对较小），反之，如果频谱图中亮的点数多，那么实际图像一定是尖锐的，边界分明且边界两边像素差异较大的。

 

对频谱移频到原点以后，可以看出图像的频率分布是以原点为圆心，对称分布的。将频谱移频到圆心除了可以清晰地看出图像频率分布以外，还有一个好处，它可以分离出有周期性规律的干扰信号，比如正弦干扰，一副带有正弦干扰，移频到原点的频谱图上可以看出除了中心以外还存在以某一点为中心，对称分布的亮点集合，这个集合就是干扰噪音产生的，这时可以很直观的通过在该位置放置带阻滤波器消除干扰。

 

频率滤波 
在频域中，频率越大说明原始信号变化速度越快；频率越小说明原始信号越平缓。当频率为0时，表示直流信号，没有变化。因此，频率的大小反应了信号的变化快慢。高频分量解释信号的突变部分，而低频分量决定信号的整体形象。

在图像处理中，频域反应了图像在空域灰度变化剧烈程度，也就是图像灰度的变化速度，也就是图像的梯度大小。对图像而言，图像的边缘部分是突变部分，变化较快，因此反应在频域上是高频分量；图像的噪声大部分情况下是高频部分；图像平缓变化部分则为低频分量。也就是说，傅立叶变换提供另外一个角度来观察图像，可以将图像从灰度分布转化到频率分布上来观察图像的特征。书面一点说就是，傅里叶变换提供了一条从空域到频率自由转换的途径。对图像处理而言，以下概念非常的重要： 

图像高频分量：图像突变部分；在某些情况下指图像边缘信息，某些情况下指噪声，更多是两者的混合； 
低频分量：图像变化平缓的部分，也就是图像轮廓信息 高通滤波器：让图像使低频分量抑制，高频分量通过 低通滤波器：与高通相反，让图像使高频分量抑制，低频分量通过 
带通滤波器：使图像在某一部分的频率信息通过，其他过低或过高都抑制 模板运算与卷积定理 
在时域内做模板运算，实际上就是对图像进行卷积。模板运算是图像处理一个很重要的处理过程，很多图像处理过程，比如增强/去噪（这两个分不清楚），边缘检测中普遍用到。

根据卷积定理，时域卷积等价与频域乘积。因此，在时域内对图像做模板运算就等效于在频域内对图像做滤波处理。 
比如说一个均值模板，其频域响应为一个低通滤波器；在时域内对图像作均值滤波就等效于在频域内对图像用均值模板的频域响应对图像的频域响应作一个低通滤波。

在MATLAB中

i=imread('lena2561.bmp');
figure(1);
imshow(i);
colorbar;

j=fft2(i);%完成FFT变换
k=fftshift(j);%%实现居中

%取得傅里叶频谱
figure(2);
l=log(abs(k));
imshow(l,[]);
colorbar;
% s=abs(j);%取得傅里叶频谱
% s=(s-min(min(s)))/(max(max(s))-min(min(s)))*255;
% imshow(s);

% 实现傅里叶逆变换
n=ifft2(j)/255;
figure(3);
imshow(n);
colorbar;

 

OpenCV

#include <stdio.h>
#include <cv.h>
#include <cxcore.h>
#include <highgui.h>
/**************************************************************************
//src IPL_DEPTH_8U
//dst IPL_DEPTH_64F
**************************************************************************/
//傅里叶正变换
void fft2(IplImage *src, IplImage *dst)
{   //实部、虚部
 IplImage *image_Re = 0, *image_Im = 0, *Fourier = 0;
 //   int i, j;
 image_Re = cvCreateImage(cvGetSize(src), IPL_DEPTH_64F, 1);  //实部
 //Imaginary part
 image_Im = cvCreateImage(cvGetSize(src), IPL_DEPTH_64F, 1);  //虚部
 //2 channels (image_Re, image_Im)
 Fourier = cvCreateImage(cvGetSize(src), IPL_DEPTH_64F, 2);
 // Real part conversion from u8 to 64f (double)
 cvConvertScale(src, image_Re, 1, 0);
 // Imaginary part (zeros)
 cvZero(image_Im);
 // Join real and imaginary parts and stock them in Fourier image
 cvMerge(image_Re, image_Im, 0, 0, Fourier);
 // Application of the forward Fourier transform
 cvDFT(Fourier, dst, CV_DXT_FORWARD);
 cvReleaseImage(&image_Re);
 cvReleaseImage(&image_Im);
 cvReleaseImage(&Fourier);
}
/**************************************************************************
//src IPL_DEPTH_64F
//dst IPL_DEPTH_8U
**************************************************************************/
void fft2shift(IplImage *src, IplImage *dst)
{
IplImage *image_Re = 0, *image_Im = 0;
int nRow, nCol, i, j, cy, cx;
double scale, shift, tmp13, tmp24;
image_Re = cvCreateImage(cvGetSize(src), IPL_DEPTH_64F, 1);
//Imaginary part
image_Im = cvCreateImage(cvGetSize(src), IPL_DEPTH_64F, 1);
cvSplit( src, image_Re, image_Im, 0, 0 );
 //具体原理见冈萨雷斯数字图像处理p123
 // Compute the magnitude of the spectrum Mag = sqrt(Re^2 + Im^2)
 //计算傅里叶谱
 cvPow( image_Re, image_Re, 2.0);
 cvPow( image_Im, image_Im, 2.0);
 cvAdd( image_Re, image_Im, image_Re);
 cvPow( image_Re, image_Re, 0.5 );
 //对数变换以增强灰度级细节(这种变换使以窄带低灰度输入图像值映射
 //一宽带输出值，具体可见冈萨雷斯数字图像处理p62)
 // Compute log(1 + Mag);
 cvAddS( image_Re, cvScalar(1.0), image_Re ); // 1 + Mag
 cvLog( image_Re, image_Re ); // log(1 + Mag)
 //Rearrange the quadrants of Fourier image so that the origin is at the image center
 nRow = src->height;
 nCol = src->width;
 cy = nRow/2; // image center
 cx = nCol/2;
//CV_IMAGE_ELEM为OpenCV定义的宏，用来读取图像的像素值，这一部分就是进行中心变换
 for( j = 0; j < cy; j++ ){
  for( i = 0; i < cx; i++ ){
   //中心化，将整体份成四块进行对角交换
   tmp13 = CV_IMAGE_ELEM( image_Re, double, j, i);
   CV_IMAGE_ELEM( image_Re, double, j, i) = CV_IMAGE_ELEM(
    image_Re, double, j+cy, i+cx);
   CV_IMAGE_ELEM( image_Re, double, j+cy, i+cx) = tmp13;
   
   tmp24 = CV_IMAGE_ELEM( image_Re, double, j, i+cx);
   CV_IMAGE_ELEM( image_Re, double, j, i+cx) =
    CV_IMAGE_ELEM( image_Re, double, j+cy, i);
   CV_IMAGE_ELEM( image_Re, double, j+cy, i) = tmp24;
  }
 }
    //归一化处理将矩阵的元素值归一为[0,255]
 //[(f(x,y)-minVal)/(maxVal-minVal)]*255
 double minVal = 0, maxVal = 0;
 // Localize minimum and maximum values
 cvMinMaxLoc( image_Re, &minVal, &maxVal );
 // Normalize image (0 - 255) to be observed as an u8 image
 scale = 255/(maxVal - minVal);
 shift = -minVal * scale;
 cvConvertScale(image_Re, dst, scale, shift);
 cvReleaseImage(&image_Re);
 cvReleaseImage(&image_Im);
}
/***********************************************************************/
int main()
{
IplImage *src;  //源图像
IplImage *Fourier;   //傅里叶系数
IplImage *dst ;
IplImage *ImageRe;
IplImage *ImageIm;
IplImage *Image;
IplImage *ImageDst;
double m,M;
double scale;
double shift;
src = cvLoadImage("D:\main.jpg",0);   //加载源图像，第二个参数表示将输入的图片转为单信道 
Fourier = cvCreateImage(cvGetSize(src),IPL_DEPTH_64F,2);
dst = cvCreateImage(cvGetSize(src),IPL_DEPTH_64F,2);
ImageRe = cvCreateImage(cvGetSize(src),IPL_DEPTH_64F,1);
ImageIm = cvCreateImage(cvGetSize(src),IPL_DEPTH_64F,1);
Image = cvCreateImage(cvGetSize(src),src->depth,src->nChannels);
ImageDst = cvCreateImage(cvGetSize(src),src->depth,src->nChannels);
fft2(src,Fourier);                  //傅里叶变换
fft2shift(Fourier, Image);          //中心化
cvDFT(Fourier,dst,CV_DXT_INV_SCALE);//实现傅里叶逆变换，并对结果进行缩放
cvSplit(dst,ImageRe,ImageIm,0,0);
cvNamedWindow("源图像",0);
cvShowImage("源图像",src);             
//对数组每个元素平方并存储在第二个参数中
cvPow(ImageRe,ImageRe,2);               
cvPow(ImageIm,ImageIm,2);
cvAdd(ImageRe,ImageIm,ImageRe,NULL);
cvPow(ImageRe,ImageRe,0.5);
cvMinMaxLoc(ImageRe,&m,&M,NULL,NULL);
scale = 255/(M - m);
shift = -m * scale;
//将shift加在ImageRe各元素按比例缩放的结果上，存储为ImageDst
cvConvertScale(ImageRe,ImageDst,scale,shift);
cvNamedWindow("傅里叶谱",0);
cvShowImage("傅里叶谱",Image);
cvNamedWindow("傅里叶逆变换",0);
cvShowImage("傅里叶逆变换",ImageDst);
//释放图像
cvWaitKey(10000);
cvReleaseImage(&src);
cvReleaseImage(&Image);
cvReleaseImage(&ImageIm);
cvReleaseImage(&ImageRe);
cvReleaseImage(&Fourier);
cvReleaseImage(&dst);
cvReleaseImage(&ImageDst);
return 0;
}