实现图像的二维傅里叶变换的频谱图
图像的频谱分析有助于我们理解图像的二维傅里叶变换（离散非周期），并且以直观的方式来展现图像的低通或高通滤波，然而如何获得图像的频谱呢？在matlab中只要短短的几行代码，就可以利用库中的函数轻松地做到。
可是，在此过程中都发生了哪些数学变换，以及得到的频谱图的含义，以及二维傅里叶变换的执行过程…都被蒙在了暗盒里。此篇文章将以应用的角度来加深对二维傅里叶变换执行过程的理解。不过，谈不上实际应用，因为普通的二维傅里叶变换实在是…太！慢！了！
由于我只熟悉java，所以用java来写。如果你熟悉的语言不是java并不会影响你对本文的理解，我会一步步解释，保证你用C/C++也能实现它。
1. 图像矩阵f(x,y)装入二维数组
这个简单，在双重循环里获取图像（i，j）坐标处的像素值（这里为了减少运算量，只获取了一个颜色通道），存到二维数组中待下一步操作
for (int i = 0;i<width;i++) {
for(int j=0;j<height;j++){
pixel[i][j]=image.getRGB(i,j)<<24>>24&0xff;
}
}
2.编写二维傅里叶变换类
F ( k , l ) = ∑ x = 0 w i d t h − 1 ∑ y = 0 h e i g h t − 1 f ( x , y ) e − i 2 π ( k x w i d t h + l y h e i g h t ) F(k,l)=\sum_{x=0}^{width-1}\sum_{y=0}^{height-1} f(x,y) e^{-i2\pi(\frac{kx}{width}+\frac{ly}{height})} F(k,l)=x=0∑width−1y=0∑height−1f(x,y)e−i2π(widthkx+heightly)
由欧拉公式：
e − i x = c o s ( x ) − i s i n ( x ) e^{-ix}=cos(x)-isin(x) e−ix=cos(x)−isin(x)
得：
F ( k , l ) = ∑ x = 0 w i d t h − 1 ∑ y = 0 h e i g h t − 1 f ( x , y ) [ c o s ( 2 π ( k x w i d t h + l y h e i g h t ) ) − i s i n ( 2 π ( k x w i d t h + l y h e i g h t ) ) ] F(k,l)=\sum_{x=0}^{width-1}\sum_{y=0}^{height-1} f(x,y)[cos(2\pi(\frac{kx}{width}+\frac{ly}{height}))-isin(2\pi(\frac{kx}{width}+\frac{ly}{height}))] F(k,l)=x=0∑width−1y=0∑height−1f(x,y)[cos(2π(widthkx+heightly))−isin(2π(widthkx+heightly))]
即最终变换出的结果应该是复数的形式(re指实部，im指虚部)
F ( k , l ) = r e + i ∗ i m F(k,l)=re+i*im F(k,l)=re+i∗im
所以在类中定义一个求实部re的方法和求虚部im的方法
static double doDFTre(int pixel[][],int width,int height,int u,int v){
double sum=0;
for(int i=0;i<width;i++){
for(int j=0;j<height;j++){
sum=sum+(pixel[i][j]*pow(-1,i+j)*cos(2*PI*u*i/width+2*PI*v*j/height));
}
}
return sum;
}//Real part
static double doDFTim(int pixel[][],int width,int height,int u,int v){
double sum=0;
for(int i=0;i<width;i++){
for(int j=0;j<height;j++){
sum=sum+(pixel[i][j]*pow(-1,i+j)*sin(2*PI*u*i/width+2*PI*v*j/height));
}
}
return sum;
}//Imaginary part
图像的频谱图中每个像素点的值表示的是以幅值的对数来刻画的（可类比地震级、亮度级）
求幅值Am
A m = r e 2 + i m 2 Am=\sqrt{re^2+im^2} Am=re2+im2
进行对数变换
方 法 的 返 回 值 是 l g ( 1 + A m ) 2 方法的返回值是lg(1+Am)^2 方法的返回值是lg(1+Am)2
static int doDFTam(int pixel[][],int width,int height,int i,int j){
Thread thread1=new Thread(()->{
re=doDFTre(pixel,width,height,i,j);
});
Thread thread2=new Thread(()->{
im=doDFTim(pixel,width,height,i,j);
});
thread1.start();
thread2.start();
Am=sqrt(pow(re,2)+pow(im,2));
return (int)pow(log(1+Am),2.3)
在多核处理器上，为了计算更快，可以分多个线程并发计算
这里设了两个线程，一个求实部，一个求虚部
3.将低频部分移动到频谱图中心
以上得出的频谱图低频部分分散在四周，高频部分在中央，为了集中低频于中央，可在步骤1中对原图f(i,j)作如下变换得到u(i,j)再对u(i,j)作二维傅里叶变换（这部分我也不理解，只是套公式）
u ( i , j ) = f ( i ， j ) ( − 1 ) i + j u(i,j)=f(i，j)(-1)^{i+j} u(i,j)=f(i，j)(−1)i+j
4.运行
可以再写一个守护线程用来计时
5.完整的代码
import javax.imageio.ImageIO;
import java.awt.image.BufferedImage;
import java.io.File;
import java.io.IOException;
import static java.lang.Math.*;
public class Fourier2D {
static int logam=0;
public static void main(String[] args) {
Thread timethread=new Thread(()->{
int i=0;
while(true){
System.out.println("Time Consumption "+i+" sec");
i=i+1;
try {
Thread.sleep(1000);
} catch (InterruptedException e) {
e.printStackTrace();
}
}
});//计时器
timethread.setDaemon(true);
timethread.setPriority(1);
timethread.start();
try {
BufferedImage image= ImageIO.read(new File("C:\\Users\\Unicode718A\\Desktop\\学习笔记\\IO\\in.png"));
int width=image.getWidth();
int height=image.getHeight();
BufferedImage outimage=new BufferedImage(width,height,2);
int pixel[][]=new int[width][height];
for (int i = 0;i<width;i++) {
for(int j=0;j<height;j++){
pixel[i][j]=image.getRGB(i,j)<<24>>24&0xff;
}
}
for(int i=0;i<width;i++){
for(int j=0;j<width;j++){
logam=DFT.doDFTam(pixel,width,height,i,j);
if(logam>255){
logam=255;
}//8-bit png图片亮度等级最大只有255
if(logam<0){
logam=0;
}
outimage.setRGB(i,j,255<<24|logam<<16|logam<<8|logam);
}
}
ImageIO.write(outimage,"png",new File("C:\\Users\\Unicode718A\\Desktop\\学习笔记\\IO\\out.png"));
} catch (IOException e) {
e.printStackTrace();
}
}
}
class DFT{
static double re;
static double im;
static double Am;
static double doDFTre(int pixel[][],int width,int height,int u,int v){
double sum=0;
for(int i=0;i<width;i++){
for(int j=0;j<height;j++){
sum=sum+(pixel[i][j]*pow(-1,i+j)*cos(2*PI*u*i/width+2*PI*v*j/height));
}
}
return sum;
}//Real part
static double doDFTim(int pixel[][],int width,int height,int u,int v){
double sum=0;
for(int i=0;i<width;i++){
for(int j=0;j<height;j++){
sum=sum+(pixel[i][j]*pow(-1,i+j)*sin(2*PI*u*i/width+2*PI*v*j/height));
}
}
return sum;
}//Imaginary part
static int doDFTam(int pixel[][],int width,int height,int i,int j){
Thread thread1=new Thread(()->{
re=doDFTre(pixel,width,height,i,j);
});//线程1
Thread thread2=new Thread(()->{
im=doDFTim(pixel,width,height,i,j);
});//线程2
thread1.start();
thread2.start();
try {
thread1.join();
thread2.join();
} catch (InterruptedException e) {
e.printStackTrace();
}
Am=sqrt(pow(re,2)+pow(im,2));
return (int)pow(log(1+Am),2.3);//Log Transformation of Amplitude
}
}
6.结果展示
100 x 100像素原图
对应频谱图
256 x 256像素的原图
对应频谱图
7.普通二维傅里叶变换运算量巨大
256 x 256 像素转换的运算时间竟然达到了10分钟，而且还是双线程并发的
而100 x 100像素转换的时间”只要”19秒
为什么会相差…这么巨大？
算算处理器执行循环的总次数就大略明白了：
100 x 100的情况：共2x100 x100 x100 x100=2 x 10^8次
256 x 256的情况：共2x256 x256 x256 x256=8.6 x10^9此
工作量增加了43倍，运算时间也会增加相应的倍数
要是图像再大点，耗时更是无法想象，所以实际应用中的傅里叶变换依赖更高效的算法
以上仅为我个人的理解，如有错误之处，还望大佬们指正。
参考：
1. http://fourier.eng.hmc.edu/e101/lectures/Image_processing/node6.html
2. https://plus.maths.org/content/fourier-transforms-images
3. 数字图像处理 冈萨雷斯 4th Edition