概要
ふと、画像をフーリエ変換するソフトをJavaで作成することになりました。
ただ、これがなかなか面倒だったので、試行錯誤の記録をメモしておきます。
FFTライブラリについて
フーリエ変換をプログラムで利用する際は通常、既存のFFTライブラリに頼ることになります。カリカリにチューニングされているものとして有名なのは、FFTWや大浦版FFTやOTFFTなどですね。今回はJavaで組むと決めていたので、このうちJava移植版があった大浦版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に掛けてみると、全く違う結果が出力されてしまいます。
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だけスケールが違うことと、そもそもデータの並ぶ順番が違うことから仕方ないと言えますね。後者の資料によると、大浦版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になります。カリカリにチューニングした2次元FFTだと更に計算量を減らせる(例1・例2)そうですが、そこまで変態的なコードを組む気はありませんので素直に1次元FFTを利用することにします。1次元FFTを行および列に掛けた際は、その都度正規化するようにしてください。
もっとも、配列が行オーダーでも列オーダーでも逆のオーダーでアクセスする場合においてメモリアクセスが大変なので、途中で転置を挟む方が良いでしょう。
また、画像化するには結果(複素数)からパワースペクトル(絶対値)を取り出して対数スケールにするのが常道なので、こんな感じに変換することにします。
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さんを例にすると、元画像
に対して、本来ならこんな感じの画像
を出したいのに、こんな感じの画像
になってしまうんですね。明るさ(スケーリング)を無視しても、対角線方向の具合が微妙に違うってのが致命的です。じゃあどうしたかって？
そりゃもう……cfftを移植するしか無いよね？
具体的には、元々の移植版に次の関数を付け足すだけです。
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を使えばこの通り。マトモな画像が出力できました。
完成したプログラム
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];
}
}
}
}