【Java】画像をフーリエ変換したい！【FFT】

#概要
　ふと、**画像をフーリエ変換するソフトをJavaで作成**することになりました。
　ただ、これがなかなか面倒だったので、試行錯誤の記録をメモしておきます。

#FFTライブラリについて
　フーリエ変換をプログラムで利用する際は通常、**既存のFFTライブラリに頼る**ことになります。カリカリにチューニングされているものとして有名なのは、[FFTW](http://www.fftw.org/)や[大浦版FFT](http://www.kurims.kyoto-u.ac.jp/~ooura/fft-j.html)や[OTFFT](http://www.moon.sannet.ne.jp/okahisa/stockham/node4.html)などですね。今回は**Javaで組む**と決めていたので、このうち[Java移植版があった](http://hp.vector.co.jp/authors/VA046927/fft4gjava.html)大浦版FFTを使用することにしました。
　ところでこの移植版、**用意されている関数がrdftしか無い**んですね。rdftとは、「実離散フーリエ変換」という意味で、**入力が(2^n個の)実数だけ受け付けますよ**ということです。どうせ入力データは画像(実数データ)なのでその点は別にいいのですが、問題はこの関数の挙動がよく分からないことです。
　ということで、適当にテストしてみました。具体的には次のコードを走らせて結果を見ます。ライブラリは(package文を削って)別途コンパイルしておきます。

```sample1.java
public class fft_test{
	public static void main(String args[]){
		double[] data = {1, 2, 3, 4, 5, 6, 7, 8};
		int n = data.length;
		FFT4g fft = new FFT4g(n);
		fft.rdft(1, data);
		for(int k = 0; k < n; k++){
			System.out.println("" + k + "," + data[k]);
		}
	}
}
```
　この結果、次のような出力を得ました。

```result1.txt
0,54.0
1,-12.0
2,3.7781745930520234
3,-28.67766952966369
4,-9.0
5,-15.0
6,-11.778174593052023
7,-6.6776695296636905
```

　これがどういった意味を持つかですが、同じデータを[Mathematicaに掛けてみる](http://www.wolframalpha.com/input/?i=Fourier[{1%2C+1%2C+2%2C+3%2C+5%2C+8%2C+13%2C+21}])と、全く違う結果が出力されてしまいます。

```result2.txt
# Fourier[{1, 1, 2, 3, 5, 8, 13, 21}]
{19.0919+0. i, 1.33579-10.1391 i, -3.18198-5.3033 i, -4.16421-2.36091 i, -4.24264+0. i, -4.16421+2.36091 i, -3.18198+5.3033 i, 1.33579+10.1391 i}
```

　ただこれは、[大浦版FFTとMathematicaの結果では√Nだけスケールが違う](http://akita-nct.jp/yamamoto/comp/mathematica/material/fft/fourier_html/node2.html)ことと、[そもそもデータの並ぶ順番が違う](http://www.math.ryukoku.ac.jp/~tsutomu/csc/csc10/csc107.pdf)ことから**仕方ない**と言えますね。後者の資料によると、大浦版FFTでは、rdftした後のデータが次のように並んでいると書かれています。

* 出力データがd[i]で、M=N/2とする
* a[0]～a[M-1]はd[0],d[2],...,d[N-1]に書かれている
* b[1]～b[M-1]はd[1],d[3],...,d[N]に書かれている
* a[M]はd[1]に書かれている

　一方、Mathematicaは素直にa[0],a[1]+b[1]i,a[2]+b[2]i,...,a[N-1]+b[N-1]iといった形式で書かれています(a[0]とa[M]では虚数成分が無いことに注意)。**どうせ[標本化定理](標本化定理)のせいでN個点があっても半分しか結果として使えない**ので合理的とも言えます。

|index|input|output|kind|
|-----|-----|------|----|
|0|1|54.0|a_0|
|1|2|-12.0|a_4|
|2|3|3.78|a_1|
|3|4|-28.7|b_1|
|4|5|-9.0|a_2|
|5|6|-15.0|b_2|
|6|7|-11.8|a_3|
|7|8|-6.68|b_3|

　ただ、これだと後述する処理を行う際に多少都合が悪いので、「入力がN個の実数である際に出力をN個の複素数にする」ことにします。まあ要するに虚数部を0として水増しすればいいんですね。

```result3.txt
# コード変更
double[] data = {1, 0, 1, 0, 2, 0, 3, 0, 5, 0, 8, 0, 13, 0, 21, 0};
# 結果
0,54.0
1,54.0    // ここだけ無視すれば(0にすれば)Mathematicaの結果と一致
2,3.7781745930520234
3,-28.67766952966369
4,-9.0
5,-15.0
6,-11.778174593052023
7,-6.677669529663689
8,-12.0
9,0.0
10,-11.778174593052023
11,6.677669529663689
12,-9.0
13,15.0
14,3.7781745930520234
15,28.67766952966369
# 参考(Fourier[{1, 1, 2, 3, 5, 8, 13, 21}]*sqrt(8))
{54.+0. i, 3.77817-28.6777 i, -9.-15. i, -11.7782-6.67767 i, -12.+0. i, -11.7782+6.67767 i, -9.+15. i, 3.77817+28.6777 i}
```

# 2次元FFTについて
　端的に言えば、**[1次元FFTを行方向および列方向に掛ければ2次元FFTになります](http://hooktail.org/computer/index.php?2%BC%A1%B8%B5%A5%D5%A1%BC%A5%EA%A5%A8%CA%D1%B4%B9)**。カリカリにチューニングした2次元FFTだと更に計算量を減らせる([例1](http://www.kurims.kyoto-u.ac.jp/~ooura/fftman/ftmn3_2.html#sec3_2)・[例2](http://www.kurims.kyoto-u.ac.jp/~ooura/fftman/ftmn3_3.html#sec3_3))そうですが、そこまで変態的なコードを組む気はありませんので素直に1次元FFTを利用することにします。[1次元FFTを行および列に掛けた際は、その都度正規化する](http://mf-atelier.sakura.ne.jp/mf-atelier/modules/tips/program/algorithm/a35.html)ようにしてください。
　もっとも、配列が行オーダーでも列オーダーでも**逆のオーダーでアクセスする場合においてメモリアクセスが大変**なので、[途中で転置を挟む方が良い](http://pokosho.com/t/image/11/)でしょう。
　また、画像化するには結果(複素数)からパワースペクトル(絶対値)を取り出して対数スケールにするのが常道なので、こんな感じに変換することにします。

```sample2.java
// 配列は行オーダーだが、大浦版FFTの仕様で実部と虚部を横に並べざるを得ない
double re = a[(x + y * w) * 2];
double im = a[(x + y * w) * 2 + 1];
double norm = re * re + im * im;
if(norm != 0.0) norm = Math.log(norm) / 2;
```

#ここにcfftがあるじゃろ？
　……さて。この考えでごく普通にrfftを二方向に使って2次元FFTすると、**絵がおかしくなります。**
　これは、「一度rfftすると結果が複素数で出てくるのに、その複素数混じりの配列相手にrfftしたから」なのですが、これを言い出すと**じゃあ何故rfftしたのか**って話になりますよね……。
　みんな大好きlenaさんを例にすると、元画像
![lena.png](https://qiita-image-store.s3.amazonaws.com/0/29802/eb026b99-da26-2de5-c76c-65147c4dd32a.png)
に対して、本来ならこんな感じの画像
![lena_FFT2.png](https://qiita-image-store.s3.amazonaws.com/0/29802/549ee58e-d6a0-d655-6ae5-556fef2ae3af.png)
を出したいのに、こんな感じの画像
![lena_FFT0.png](https://qiita-image-store.s3.amazonaws.com/0/29802/13db5d60-0ebd-279d-db2c-fc3688c62166.png)
になってしまうんですね。明るさ(スケーリング)を無視しても、**対角線方向の具合が微妙に違う**ってのが致命的です。じゃあどうしたかって？
　そりゃもう……**cfftを移植するしか無いよね？**
　具体的には、[元々の移植版](http://hp.vector.co.jp/authors/VA046927/fft4gjava.html)に次の関数を付け足すだけです。

```FFT4g+.java
	public void cdft(int isgn, double[] a)
	{
	    if (n > (ip[0] << 2)) {
	        makewt(n >> 2);
	    }
	    if (n > 4) {
			if (isgn >= 0) {
	            bitrv2(n, a);
	            cftfsub(a);
	    	} else {
	            bitrv2conj(n, a);
	            cftbsub(a);
	        }
	    }else if (n == 4) {
	        cftfsub(a);
	    }
	}
	private void bitrv2conj(int n, double[] a)
	{
	    int j, j1, k, k1, l, m, m2;
	    double xr, xi, yr, yi;
	    
	    ip[2 + 0] = 0;
	    l = n;
	    m = 1;
	    while ((m << 3) < l) {
	        l >>= 1;
	        for (j = 0; j < m; j++) {
	            ip[2 + m + j] = ip[2 + j] + l;
	        }
	        m <<= 1;
	    }
	    m2 = 2 * m;
	    if ((m << 3) == l) {
	        for (k = 0; k < m; k++) {
	            for (j = 0; j < k; j++) {
	                j1 = 2 * j + ip[2 + k];
	                k1 = 2 * k + ip[2 + j];
	                xr = a[j1];
	                xi = -a[j1 + 1];
	                yr = a[k1];
	                yi = -a[k1 + 1];
	                a[j1] = yr;
	                a[j1 + 1] = yi;
	                a[k1] = xr;
	                a[k1 + 1] = xi;
	                j1 += m2;
	                k1 += 2 * m2;
	                xr = a[j1];
	                xi = -a[j1 + 1];
	                yr = a[k1];
	                yi = -a[k1 + 1];
	                a[j1] = yr;
	                a[j1 + 1] = yi;
	                a[k1] = xr;
	                a[k1 + 1] = xi;
	                j1 += m2;
	                k1 -= m2;
	                xr = a[j1];
	                xi = -a[j1 + 1];
	                yr = a[k1];
	                yi = -a[k1 + 1];
	                a[j1] = yr;
	                a[j1 + 1] = yi;
	                a[k1] = xr;
	                a[k1 + 1] = xi;
	                j1 += m2;
	                k1 += 2 * m2;
	                xr = a[j1];
	                xi = -a[j1 + 1];
	                yr = a[k1];
	                yi = -a[k1 + 1];
	                a[j1] = yr;
	                a[j1 + 1] = yi;
	                a[k1] = xr;
	                a[k1 + 1] = xi;
	            }
	            k1 = 2 * k + ip[k];
	            a[k1 + 1] = -a[k1 + 1];
	            j1 = k1 + m2;
	            k1 = j1 + m2;
	            xr = a[j1];
	            xi = -a[j1 + 1];
	            yr = a[k1];
	            yi = -a[k1 + 1];
	            a[j1] = yr;
	            a[j1 + 1] = yi;
	            a[k1] = xr;
	            a[k1 + 1] = xi;
	            k1 += m2;
	            a[k1 + 1] = -a[k1 + 1];
	        }
    } else {
        a[1] = -a[1];
        a[m2 + 1] = -a[m2 + 1];
        for (k = 1; k < m; k++) {
            for (j = 0; j < k; j++) {
                j1 = 2 * j + ip[k];
                k1 = 2 * k + ip[j];
                xr = a[j1];
                xi = -a[j1 + 1];
                yr = a[k1];
                yi = -a[k1 + 1];
                a[j1] = yr;
                a[j1 + 1] = yi;
                a[k1] = xr;
                a[k1 + 1] = xi;
                j1 += m2;
                k1 += m2;
                xr = a[j1];
                xi = -a[j1 + 1];
                yr = a[k1];
                yi = -a[k1 + 1];
                a[j1] = yr;
                a[j1 + 1] = yi;
                a[k1] = xr;
                a[k1 + 1] = xi;
            }
            k1 = 2 * k + ip[k];
            a[k1 + 1] = -a[k1 + 1];
            a[k1 + m2 + 1] = -a[k1 + m2 + 1];
        }
    }
}
```
　で、二次元FFTにこのcfftを使えばこの通り。マトモな画像が出力できました。
![lena_FFT1.png](https://qiita-image-store.s3.amazonaws.com/0/29802/ed66ce5a-aca7-e618-2464-2ca03bddc9a3.png)


#完成したプログラム
```2d_fft.java
/* フーリエ変換で見る画像解像度
 * 参考：
 * 「FFT (高速フーリエ・コサイン・サイン変換) の概略と設計法」
 * http://www.kurims.kyoto-u.ac.jp/~ooura/fftman/index.html
 * 「大浦版FFTのJava移植」
 * http://hp.vector.co.jp/authors/VA046927/fft4gjava.html
 * 「第11章 周波数で処理する」
 * http://pokosho.com/t/image/11/
 * 「備忘録: 画像処理での二次元フーリエ変換」
 * http://bebolog.blogspot.jp/2014/12/blog-post_15.html
 */

import java.awt.image.BufferedImage;
import java.io.File;
import javax.imageio.ImageIO;

public class show_resolution{
	public static void main(String args[]){
		if(args.length < 1) return;
		try{
			BufferedImage image = ImageIO.read(new File(args[0]));
			int w = image.getWidth(), h = image.getHeight();
			double[] image_data = new double[w * h * 2];
			// データを配列に代入する
			for(int y = 0; y < h; ++y){
				for(int x = 0; x < w; ++x){
					int color = image.getRGB(x, y);
					int r = (color & 0xff0000) >> 16;
					int g = (color & 0xff00) >> 8;
					int b = color & 0xff;
					image_data[(x + y * w) * 2] = 0.299 * r + 0.587 * g + 0.114 * b;
					image_data[(x + y * w) * 2 + 1] = 0.0;
				}
			}
			// 二次元FFT
			fft2d(image_data, w, h);
			// パワースペクトルに変換
			power_spectral(image_data, w, h);
			// 象限入れ替え
			swap_quadrants(image_data, w, h);
			// 正規化
			normalization(image_data, w, h);
			// 出力
			for(int y = 0; y < h; ++y){
				for(int x = 0; x < w; ++x){
					int Y = (int)(image_data[(x + y * w)] * 255);
					if(Y > 255) Y = 255;
					if(Y < 0) Y = 0;
					image.setRGB(x, y, 0x010101 * Y + 0xff000000);
				}
			}
			String file_name = args[0].substring(0, args[0].lastIndexOf(".")) + "_FFT1.png";
			ImageIO.write(image, "png", new File(file_name));
		}catch(Exception error){
			error.printStackTrace();
		}
	}
	// 正規化
	static void normalization(double[] a, int w, int h){
		double min = a[0], max = a[0];
		for(int k = 1; k < w * h; ++k){
			min = Math.min(min, a[k]);
			max = Math.max(max, a[k]);
		}
		double diff = max - min;
		for(int k = 0; k < w * h; ++k){
			a[k] = (a[k] - min) / diff;
		}
	}
	// 象限入れ替え
	static void swap_quadrants(double[] a, int w, int h){
		int hw = w / 2, hh = h / 2;
		double[] b = new double[w * h];
		for(int y = 0; y < hh; ++y){
			for(int x = 0; x < hw; ++x){
				b[(y + hh) * w + x] = a[y * w + (x + hw)];	//第1象限
				b[(y + hh) * w + (x + hw)] = a[y * w + x];	//第2象限
				b[y * w + (x + hw)] = a[(y + hh) * w + x];	//第3象限
				b[y * w + x] = a[(y + hh) * w + (x + hw)];	//第4象限
			}
		}
		for(int k = 0; k < w * h; ++k){
			a[k] = b[k];
		}
	}
	// パワースペクトルに変換
	static void power_spectral(double[] a, int w, int h){
		for(int y = 0; y < h; ++y){
			for(int x = 0; x < w; ++x){
				double re = a[(x + y * w) * 2];
				double im = a[(x + y * w) * 2 + 1];
				double norm = re * re + im * im;
				if(norm != 0.0) norm = Math.log(norm) / 2;
				a[x + y * w] = norm;
			}
		}
	}
	// 2次元FFT
	static void fft2d(double[] a, int w, int h){
		double[] b = new double[w * h * 2];
		// 水平方向のFFT
		for(int y = 0; y < h; ++y){
			fft1d(a, w * 2, y * w * 2);
		}
		// 転置操作
		transpose(a, b, w, h);
		// 垂直方向のFFT
		for(int x = 0; x < w; ++x){
			fft1d(b, h * 2, x * h * 2);
		}
		// 転置操作
		transpose(b, a, h, w);
	}
	// 1次元FFT
	static void fft1d(double[] a, int n, int p){
		double[] temp = new double[n];
		for(int k = 0; k < n; ++k){
			temp[k] = a[p + k];
		}
		FFT4g fft = new FFT4g(n);
		fft.cdft(1, temp);
		for(int k = 0; k < n; ++k){
			a[p + k] = temp[k] / n * 2;
		}
	}
	// 行列の転置
	static void transpose(double[] src, double[] dst, int w, int h){
		for(int y = 0; y < h; ++y){
			for(int x = 0; x < w; ++x){
				int p = x + y * w;
				int q = y + x * h;
				dst[q * 2    ] = src[p * 2];
				dst[q * 2 + 1] = src[p * 2 + 1];
			}
		}
	}
}
```