图像 快速傅里叶变换 及 频率域滤波 java 实现

首先感谢中山大学12级软件学院计算机应用方向和我同班的乔勃大God，以及软件学院副院长、数图seisei朝老师的帮助！马屁还是要拍的o(*￣▽￣*)ブ

//-------------------------------------------------OK 正题---------------------------------------------------------

很多学过数图的小伙伴都经历过傅里叶变换这一坑，如果使用中文书的话错误就更多了。

So，今天小渣我就用java实现一个快速傅里叶变换和频率域的滤波，函数都是纯手写（用来弄懂原理）

废话不多说，先看看频率域滤波的大致步骤：

1、给定一副A*B的图像，对其进行扩充（补0），使其长和宽变为大于或等于自身长度的最小2的整数次幂。

2、对扩充后的图像进行移动中心操作，即乘以（-1）^(x+y)。

3、计算第二步的FFT（利用一维快速傅里叶变换实现二维）【生成一个大小和扩充图像一致的二维复数数组F（u，v）】

4、生成一个频率域的滤波器H（u，v）【注意这个是复数，大小和第三步的结果一样】

5、点乘，令g（u，v）= H（u，v）·F（u，v）

6、对g做快速傅里叶逆变换IFFT

7、对第6步的结果取实部

8、队第7步的结果乘以（-1）^(x+y)再次移动中心。

9、剪裁获得与原图大小一致的结果。

sorry，今天好累，先开个头，把代码放上去，明天继续。

/*   这个文件里包含两个DFT函数，一个是主DFT，在这里完成一个二维图片的
扩充、中心移动、傅里叶变换、滤波（包括滤波器的处理等）、逆变换、取实部、
移回中心等操作并返回一个处理完后的二维图片。另一个DFT是指普通的（慢速）
傅里叶变换，输入一个二维数组、u、v并返回一个值F（u、v）。
除此之外，本文件还提供：
·get2PowerEdge：用来计算一个图像边的最小2的整数次幂扩展；
·showFouriedImage：返回一个傅里叶变换后的频谱图（取模、取log再量化）
·fft：一维快速傅里叶变换
·ifft：一维快速傅里叶逆变换
·convolve：一维卷积。*/

import java.awt.Image;
import java.awt.image.BufferedImage;
import java.io.File;
import java.io.IOException;
import java.lang.reflect.Array;
import java.util.Scanner;
import java.lang.Math;

import javax.imageio.ImageIO;

public class ImageFourier {

public static void main (String[] args) throws IOException {
BufferedImage srcimg = ImageIO.read(new File("./84.png"));
BufferedImage destimg = DFT(srcimg);
ImageIO.write(destimg, "png", new File("G:/pic/NormalDFT84.png"));
System.out.println("finished!");
}

public static BufferedImage DFT(BufferedImage img) throws IOException {

int w = img.getWidth(null);
int h = img.getHeight(null);
int m = get2PowerEdge(w); // 获得2的整数次幂
//		System.out.println(m);
int n = get2PowerEdge(h);
//		System.out.println(n);
int[][] last = new int[m]
;
Complex[][] next = new Complex[m]
;
int pixel, alpha = -1, newred, newgreen, newblue, newrgb;

BufferedImage destimg = new BufferedImage(w, h, BufferedImage.TYPE_INT_ARGB);
//----------------------------------------------------------------------
//  first: Image Padding and move it to center 填充图像至2的整数次幂并乘以（-1）^(x+y)
//  use 2-D array last to store
for (int i = 0; i < m; i++) {
for (int j = 0; j < n; j++) {
if (i < w && j < h) {
pixel = img.getRGB(i, j);
if ((i+j)%2==0) {
newred = pixel&0x00ff0000>>16;
}
else {
newred = -(pixel&0x00ff0000>>16);
}
last[i][j] = newred;
}
else {
last[i][j] = 0;
}
}
}

//----------------------------------------------------------------------
// second: Fourier Transform 离散傅里叶变换
// u-width v-height x-width y-height

//------------Normal-DFT-------------
//			for (int u = 0; u < m; u++){
//				for (int v = 0; v <n; v++) {
//					next[u][v] = DFT(last, u, v);
//					System.out.println("U: "+u+"---v: "+v);
//				}
//			}
//			if (true) { // 生成DFT图片,记得修改图片大小
//			destimg = showFourierImage(next);
//			return destimg;
//		}

//---------------FFT-----------------
// 先把所有的行都做一维傅里叶变换，再放回去
Complex[] temp1 = new Complex
;
for (int x = 0; x < m; x++) {
for (int y = 0; y < n; y++) {
Complex c = new Complex(last[x][y],0);
temp1[y] = c;
}
next[x] = fft(temp1);
}

// 再把所有的列（已经被行的一维傅里叶变换所替代）都做一维傅里叶变换
Complex[] temp2 = new Complex[m];
for (int y = 0; y < n; y++) {
for (int x = 0; x < m; x++) {
Complex c = next[x][y];
temp2[x] = c;
}
temp2 = fft(temp2);
for (int i = 0; i < m; i++) {
next[i][y] = temp2[i];
}
}

//		if (true) { // 生成DFT图片,记得修改图片大小
//			destimg = showFourierImage(next);
//			return destimg;
//		}

//----------------------------------------------------------------------
// third: Generate the frequency filter and filter the image in frequency domain 生成频率域滤波器并滤波

// 构造原始滤波函数
Complex[][] filter = new Complex[m]
;
//这个是11X11均值滤波
//		for (int x = 0; x < m; x++) {
//			for (int y = 0; y < n; y++) {
//				if (x < 11 && y < 11) {
//					if ((x+y)%2==0)
//						filter[x][y] = new Complex(1/121d, 0); // double 后面赋值数字记得加d！！！！！！！
//					else
//						filter[x][y] = new Complex(-1/121d, 0);
//				}
//				else {
//					filter[x][y] = new Complex(0, 0);
//				}
//			}
//		}

//下面这个是拉普拉斯滤波
filter[0][0] = new Complex(0, 0);
filter[0][1] = new Complex(-1, 0);
filter[0][2] = new Complex(0, 0);
filter[1][0] = new Complex(-1, 0);
filter[1][1] = new Complex(4, 0);
filter[1][2] = new Complex(-1, 0);
filter[2][0] = new Complex(0, 0);
filter[2][1] = new Complex(-1, 0);
filter[2][2] = new Complex(0, 0);
for (int x = 0; x < m; x++) {
for (int y = 0; y < n; y++) {
if (x < 3 && y < 3) {/*上面已经写好了*/}
else {
filter[x][y] = new Complex(0, 0);
}
}
}

// 傅里叶变换 转换为频率域
for (int x = 0; x < m; x++) {
for (int y = 0; y < n; y++) {
Complex c = new Complex(filter[x][y].getR(), filter[x][y].getI());
temp1[y] = c;
}
filter[x] = fft(temp1);
}

for (int y = 0; y < n; y++) {
for (int x = 0; x < m; x++) {
Complex c = new Complex(filter[x][y].getR(), filter[x][y].getI());
//				Complex c = filter[x][y];
temp2[x] = c;
}
temp2 = fft(temp2);
for (int i = 0; i < m; i++) {
filter[i][y] = temp2[i];
}
}

//		if (true) {
//			destimg = showFourierImage(filter);
//			return destimg;
//		}

// point-wise multiply
Complex[][] g = new Complex[m]
;
for (int x = 0; x < m; x++) {
for (int y = 0; y < n; y++) {
g[x][y] = filter[x][y].times(next[x][y]);
//				System.out.println("g: "+g[x][y].getR()+"  "+g[x][y].getI());
}
}
//----------------------------------------------------------------------
// fourth: use IDFT to get the image 傅里叶逆变换
for (int x = 0; x < m; x++) {
for (int y = 0; y < n; y++) {
Complex c = new Complex(g[x][y].getR(), g[x][y].getI());
temp1[y] = c;
}
g[x] = ifft(temp1);
}

//		for (int x = 0; x < m; x++) {
//			for (int y = 0; y < n; y++) {
//				System.out.println("gifft-g: "+g[x][y].getR()+"  "+g[x][y].getI());
//			}
//		}

for (int y = 0; y < n; y++) {
for (int x = 0; x < m; x++) {
Complex c = g[x][y];
temp2[x] = c;
}
temp2 = ifft(temp2);
for (int i = 0; i < m; i++) {
g[i][y] = temp2[i];
}
}

for (int x = 0; x < m; x++) {
for (int y = 0; y < n; y++) {
//				System.out.println("ifft-g: "+g[x][y].getR()+"  "+g[x][y].getI());
}
}
//----------------------------------------------------------------------
// fifth：取实部
for (int x = 0; x < m; x++) {
for (int y = 0; y < n; y++) {
last[x][y] = (int)g[x][y].getR();
//				System.out.println(last[x][y]);
}
}
//----------------------------------------------------------------------
// sixth: move the image back and cut the image 乘以(-1)^(x+y)再剪裁图像
//		int srcpixel, srcred;
int newalpha = (-1) << 24;
for (int i = 0; i < w; i++) {
for (int j = 0; j < h; j++) {
//	srcpixel = img.getRGB(i, j);
//	srcred = srcpixel&0x00ff0000>>16;
newred = last[i][j];
if ((i+j)%2!=0)
newred = -newred;
//	newred = srcred-newred;
newblue = newred; // 先写这个 ，如果先改变newred的值，newblue也会变成改过后的newred！
newgreen = newred << 8; // 这个也一样，反正不能放到newred改变自己之前！
newred = newred << 16;
newrgb = newalpha | newred | newgreen | newblue;
destimg.setRGB(i, j, newrgb);
//				System.out.println("R: "+newred+"---G: "+newgreen+"---B: "+newblue);
}
}
//----------------------------------------------------------------------
return destimg;
}

//---------------------other functions--------------------------

// 根据图像的长获得2的整数次幂
public static int get2PowerEdge(int e) {
if (e == 1)
return 1;
int cur = 1;
while(true) {
if (e > cur && e <= 2 * cur)
return 2*cur;
else
cur *= 2;
}
}

// 返回傅里叶频谱图
public static BufferedImage showFourierImage (Complex[][] f) {
int w = f.length;
int h = f[0].length;
double max = 0;
double min = 0;
BufferedImage destimg = new BufferedImage(w, h, BufferedImage.TYPE_INT_ARGB);
//-------------------First get abs(取模)--------------------------
double[][] abs = new double[w][h];
for (int i = 0; i < w; i++) {
for (int j = 0; j < h; j++) {
abs[i][j] = f[i][j].abs();
//				System.out.println(f[i][j].getR()+"  "+f[i][j].getI());
}
}
//-------------------Second get log(取log + 1)-------------------
for (int i = 0; i < w; i++) {
for (int j = 0; j < h; j++) {
abs[i][j] = Math.log(abs[i][j]+1);
}
}
//-------------------Third quantization(量化)---------------------
max = abs[0][0];
min = abs[0][0];
for (int i = 0; i < w; i++) {
for (int j = 0; j < h; j++) {
if (abs[i][j] > max)
max = abs[i][j];
if (abs[i][j] < min)
min = abs[i][j];
}
}
int level = 255;
double interval = (max - min) / level;
for (int i = 0; i < w; i++) {
for (int j = 0; j < h; j++) {
for (int k = 0; k <= level; k++) {
if (abs[i][j] >= k * interval && abs[i][j] < (k + 1) * interval) {
abs[i][j] = (k * interval / (max - min)) * level;
break;
}
}
}
}
//-------------------Fourth setImage----------------------------
int newalpha = (-1) << 24;
int newred;
int newblue;
int newgreen;
int newrgb;
for (int i = 0; i < w; i++) {
for (int j = 0; j < h; j++) {
newred = (int)abs[i][j] << 16;
newgreen = (int)abs[i][j] << 8;
newblue = (int)abs[i][j];
newrgb = newalpha | newred | newgreen | newblue;
destimg.setRGB(i, j, newrgb);
}
}
return destimg;
}

// normal 2-D DFT
public static Complex DFT(int[][] f, int u, int v) {
int M = f.length;
int N = f[0].length;
Complex c = new Complex(0, 0);
for (int x = 0; x < M; x++) {
for (int y = 0; y < N; y++) {
Complex temp = new Complex(0, -2*Math.PI*(u*x/M+v*y/N));
c = c.plus(temp.exp().times(f[x][y]));
}
}
return c;
}

// 快速一维傅里叶变换
public static Complex[] fft (Complex[] x) { //传入的全都是206
int N = x.length;
//		if (N == 256) {
//			for (int i = 0; i < N; i++)
//			System.out.println(i+"---"+x[i].getR()+"  "+x[i].getI());
//		}
//		System.out.println(N);

// base case
if (N == 1) {
//			System.out.println(x[0].getR()+"  "+x[0].getI()); // !!!!ERROR
//			return new Complex[] {x[0]};
return x;
}

// radix 2 Cooley-Turkey FFT
if (N % 2 != 0) {
throw new RuntimeException("N is not a power of 2");
}

// fft of even terms
Complex[] even = new Complex[N/2];
for (int k = 0; k < N/2; k++) {
even[k] = x[2*k];
}
Complex[] q = fft(even);

// fft of odd terms
Complex[] odd = new Complex[N/2];
for (int k = 0; k < N/2; k++) {
odd[k] = x[2*k+1]; //DEBUG  之前这里忘记+1  差点搞死我
}
Complex[] r = fft(odd);

// combine
Complex[] y = new Complex
;
for (int k = 0; k < N/2; k++) {
double kth = -2 * k * Math.PI / N;
Complex wk = new Complex(Math.cos(kth), Math.sin(kth)); // all small number not 0
y[k] = q[k].plus(wk.times(r[k]));
y[k + N/2] = q[k].minus(wk.times(r[k]));
//			System.out.println("wk: "+N+"---"+wk.getR()+"  "+wk.getI());
//			System.out.println("q[k]: "+N+"---"+q[k].getR()+"  "+q[k].getI());
//			System.out.println("r[k]: "+N+"---"+r[k].getR()+"  "+r[k].getI());
//			System.out.println("wk.times(r[k]): "+N+"---"+wk.times(r[k]).getR()+"  "+wk.times(r[k]).getI());
}

return y;
}

// 快速一维傅里叶逆变换
public static Complex[] ifft(Complex[] x) {
int N = x.length;
Complex[] y = new Complex
;

// take conjugate
for (int i = 0; i < N; i++) {
y[i] = x[i].conjugate();
}

// compute forward fft
y = fft(y);

// take conguate again
for (int i = 0; i < N; i++) {
y[i] = y[i].conjugate();
}

// divide by N
for (int i = 0; i < N; i++) {
y[i] = y[i].times(1.0/N);
}

return y;
}

// 快速一维卷积
public Complex[] convolve(Complex[] x, Complex[] y) {

if (x.length != y.length) {
throw new RuntimeException("Dimension don't agree");
}

int N = x.length;

// compute fft of each sequence;
Complex[] a = fft(x);
Complex[] b = fft(y);

// point-wise multiply
Complex[] c = new Complex
;
for (int i = 0; i < N; i++) {
c[i] = a[i].times(b[i]);
}

// compute inverse FFT
//		return ifft(c);
return c;
}
}

// complex inner class
public class Complex {
private final double r;
private final double i;

public Complex (double r, double i) {
this.r = r;
this.i = i;
}

public double abs() { // return sqrt(r^2 +i^2)
return Math.hypot(r, i);
}

public double phase() {
return Math.atan2(i, r);
}

public Complex plus (Complex c) {
return new Complex (this.r + c.r, this.i + c.i);
}

public Complex minus (Complex c) {
return new Complex (this.r - c.r, this.i - c.i);
}

public Complex times (Complex c) {
return new Complex (this.r * c.r - this.i * c.i,
this.r * c.i + this.i * c.r);
}

public Complex times (double d) {
return new Complex (this.r * d, this.i * d);
}

public Complex conjugate() {
return new Complex (r, -i);
}

public double getR () {
return r;
}

public double getI () {
return i;
}

public Complex exp() {
return new Complex(Math.exp(r) * Math.cos(i),
Math.exp(r) * Math.sin(i));
}
}