チーム15C2
課題名
FFTを用いた画像の周波数成分の可視化
研究者名
3年14組42番 Kyu Matuura
3年14組45番 Satoru Murakami
概要
高速フーリエ変換(FFT)を用いて、画像の周波数成分を可視化するプログラムを作成した。
高速フーリエ変換というのはデジタル信号(いわゆる離散値)に対して行う離散フーリエ変換を 処理速度を速めたものである。フーリエ変換を行うと、周波数成分に変換をすることができる。
変換の手順を示すと、
変換したい画像を入力
→入力画像をグレースケール処理
→処理した画像をFFTで変換
→変換した値(絶対値)を色で区別したものを出力
実際に画像を処理しながら見ていく。
まず、入力画像は以下のものとする。
この画像をグレースケール処理すると以下のようになる。つまり、白黒の画像に変換する。
そして、この処理をしてから高速フーリエ変換を行う。高速フーリエ変換を行うことで周波数成分
に変換することができると言ったが値は複素数で表現される。複素数のままでは大きさの判断というのが し辛いので絶対値(=パワースペクトラム)に直してから値の比較を行っている。そして、値ごとに色を 付けることで成分の分布を表現した画像が出力される。その画像が以下のものである。
画像の四隅は低周波成分、中央は高周波成分を表している。
赤い部分は、その周波数成分が多く画像に含まれていることを表し、青い部分は含まれていないことを表す。
以上の処理を行うプログラムを今回実装した。
FFTプログラムの内容
プログラムは以下の4つに分かれている(記述はjava)。簡単な説明も付随する。
・Complex:複素数を表すクラス。
・FFT:高速フーリエ変換の計算を行うクラス
・FFT2D:高速フーリエ変換を2次元の高速フーリエ変換に拡張したクラス
・ImageFFT:周波数成分ごとに色を区別して表示するクラス
それぞれソースコードを載せる。
(見たいプログラムの名前をクリックしてください)
Complex.java
package myfft;
import java.lang.Math;
/**
* 複素数クラス */
public final class Complex{
/** 実部 */ private final double re; /** 虚部 */ private final double im;
/** 実部虚部を初期化するコンストラクタ */ public Complex(final double re, final double im){ this.re = re; this.im = im; }
/** 実部のみを初期化するコンストラクタ */ public Complex(final double re){ this(re, 0); }
/** 実部reのゲッター*/ public double getRe(){ return re; }
/** 虚部imのゲッター*/ public double getIm(){ return im; }
/** 和 */ public static Complex add(final Complex c1, final Complex c2){ return new Complex(c1.re + c2.re, c1.im + c2.im); }
/** 差 */ public static Complex sub(final Complex c1, final Complex c2){ return new Complex(c1.re - c2.re, c1.im - c2.im); }
/** 積 */ public static Complex mult(final Complex c1, final Complex c2){ return new Complex(c1.re * c2.re - c1.im * c2.im , c1.re * c2.im + c1.im * c2.re); }
/** スカラー積 */ public static Complex mult(final Complex c, final double k){ return new Complex(k * c.re, k * c.im); }
/** スカラー商 */ public static Complex div(final Complex c, final double k){ return new Complex(c.re / k, c.im / k); }
/** 共役 */ public Complex conjugate(){ return new Complex(re, -im); }
/** 長さ */ public double norm(){ return Math.sqrt(re * re + im * im); }
/** 比較 */ @Override public boolean equals(final Object obj){ if(obj instanceof Complex){ Complex c = (Complex)obj; return c.re == re && c.im == im; } return false; }
/** ハッシュ */ @Override public int hashCode(){ return Double.valueOf(re).hashCode() ^ Double.valueOf(im).hashCode(); }
/** 文字列表現 */ @Override public String toString(){ return String.format("(% 3.2f,% 3.2f)", re, im); }
}
FFT.java
package myfft;
import java.lang.Math;
/**
* FFTクラス */
public final class FFT{
/** FFTの点数を表すn */ private final int n; /** 2を底とするNの対数 */ private final int logN; /** 回転因子 */ private final Complex w[];
/** N点FFT用オブジェクトのコンストラクタ */ public FFT(final int n){ if(n < 0 | | (n & (n - 1)) != 0){ throw new IllegalArgumentException("n must be power of 2"); } this.n = n;
// calc log_2 n int count = 0; for(int i = n; i > 1; i /= 2){ count++; } logN = count;
// prepare twiddle factors w = new Complex[n / 2]; for(int i = 0; i < n / 2; i++){ double theta = 2 * i * Math.PI / n; w[i] = new Complex(Math.cos(theta), -Math.sin(theta)); } }
/** * 複素数配列xを離散フーリエ変換する */ public void fft(final Complex x[]){ if(x == null){ throw new NullPointerException(); } if(x.length != n){ throw new IllegalArgumentException("x length must be n"); }
// N point DFT for(int m = n; m >= 2; m /= 2){ int hm = m / 2; // half of m // do the DFT begins at i for(int i = 0; i < n; i += m){ // butterfly tx[i + j] and tx[i + hn + j]; for(int j = 0; j < hm; j++){ final int k = i + j; final int l = i + hm + j; final Complex xk = x[k]; final Complex xl = x[l]; x[k] = Complex.add(xk, xl); x[l] = Complex.mult(Complex.sub(xk, xl), w[j*n/m]); } } }
for(int i = 0; i < n; i++){ x[i] = Complex.div(x[i], n); }
bitReverseSort(x); }
/** * 複素数配列xを離散フーリエ逆変換する */ public void ifft(final Complex x[]){ if(x == null){ throw new NullPointerException(); } if(x.length != n){ throw new IllegalArgumentException("x length must be n"); }
for(int i = 0; i < n; i++){ x[i] = x[i].conjugate(); }
fft(x);
for(int i = 0; i < n; i++){ x[i] = Complex.mult(x[i].conjugate(), n); } }
/** * 配列をビットリバース順に並べ替えるメソッド */ private void bitReverseSort(final Complex x[]){ for(int i = 0; i < n; i++){ int j = bitReverse(i, logN); if(i < j){ Complex tmp = x[i]; x[i] = x[j]; x[j] = tmp; } } }
/** * 変数nのビットリバースを計算するメソッド */ private static int bitReverse(final int n){ int w = n; w = (w & 0xAAAAAAAA) >>> 1 | (w & 0x55555555) << 1; w = (w & 0xCCCCCCCC) >>> 2 | (w & 0x33333333) << 2; w = (w & 0xF0F0F0F0) >>> 4 | (w & 0x0F0F0F0F) << 4; w = (w & 0xFF00FF00) >>> 8 | (w & 0x00FF00FF) << 8; w = (w & 0xFFFF0000) >>> 16 | (w & 0x0000FFFF) << 16; return w; }
/** * 下位sizeビットを変数nのビットリバースを計算するメソッド */ private static int bitReverse(final int n, final int size){ final int SIZE_OF_INT = 32; assert 0 <= size && size < SIZE_OF_INT;
final int rev = bitReverse(n); return rev >>> (SIZE_OF_INT - size); }
/** * テスト */ public static void main(String args[]){ int n = 32; FFT f = new FFT(32); // 矩形波 Complex[] x = new Complex[n]; for(int i = 0; i < n; i += 16){ for(int j = 0; j < 8; j++){ x[i+j] = new Complex(0); x[i+j+8] = new Complex(1); } } System.out.println("Original"); for(Complex c : x){ System.out.println(c); } f.fft(x); System.out.println("FFT"); for(Complex c : x){ System.out.println(c); } f.ifft(x); System.out.println("IFFT"); for(Complex c : x){ System.out.println(c); } }
}
FFT2D
package myfft;
import java.lang.Math;
/**
* 2次元FFTクラス */
public final class FFT2D{
/** FFTの点数を表すn */ private final int n; /** FFTの点数を表すm */ private final int m; /** n点FFT */ private final FFT nPointFFT; /** m点FFT */ private final FFT mPointFFT;
/** nxm点FFT用オブジェクトのコンストラクタ */ public FFT2D(final int n, final int m){ nPointFFT = new FFT(n); mPointFFT = new FFT(m); this.n = n; this.m = m; }
/** * 2次元複素数配列xを離散フーリエ変換する */ public void fft(final Complex x[][]){ if(x == null){ throw new NullPointerException(); } if(x.length != n){ throw new IllegalArgumentException("x length must be n"); } if(x[0].length != m){ throw new IllegalArgumentException("x length must be m"); }
// 行ごとにFFT Complex tmp[]; tmp = new Complex[m]; for(int i = 0; i < n; i++){ for(int j = 0; j < m; j++){ tmp[j] = x[i][j]; }
mPointFFT.fft(tmp);
for(int j = 0; j < m; j++){ x[i][j] = tmp[j]; } }
// 列ごとにFFT tmp = new Complex[n]; for(int i = 0; i < m; i++){ for(int j = 0; j < n; j++){ tmp[j] = x[j][i]; }
nPointFFT.fft(tmp);
for(int j = 0; j < n; j++){ x[j][i] = tmp[j]; } } }
/** * 2次元複素数配列xを離散フーリエ逆変換する */ public void ifft(final Complex x[][]){ if(x == null){ throw new NullPointerException(); } if(x.length != n){ throw new IllegalArgumentException("x length must be n"); } if(x[0].length != m){ throw new IllegalArgumentException("x length must be m"); }
// 列ごとにIFFT Complex tmp[]; tmp = new Complex[m]; for(int i = 0; i < n; i++){ for(int j = 0; j < m; j++){ tmp[j] = x[i][j]; }
mPointFFT.ifft(tmp);
for(int j = 0; j < m; j++){ x[i][j] = tmp[j]; } }
// 列ごとにIFFT tmp = new Complex[n]; for(int i = 0; i < m; i++){ for(int j = 0; j < n; j++){ tmp[j] = x[j][i]; }
nPointFFT.ifft(tmp);
for(int j = 0; j < n; j++){ x[j][i] = tmp[j]; } } }
}
ImageFFT.java
package img;
import myfft.*;
import java.awt.*;
import java.awt.image.*;
import javax.imageio.*;
import java.io.*;
/**
* 2次元FFTを用いて、グレースケール画像のパワースペクトルを表示する */
public class ImageFFT{
public static void main(String args[]){ try{ // 引数チェック if(args.length == 0){ System.err.println("JPGファイルを引数に与えてください."); return; }
// 画像の読み込み BufferedImage src = ImageIO.read(new File(args[0]));
// 画像をグレースケールに変換 BufferedImage gray = toGrayScale(src); ImageIO.write(gray, "jpg", new File("gray.jpg"));
// 画像を2の累乗の解像度にリサイズする BufferedImage expandedSrc = resizePowerOfTwo(gray);
// 画像を複素数配列に変換して,離散フーリエ変換する Complex c[][] = imageToComplexArray(expandedSrc); FFT2D ft = new FFT2D(c.length, c[0].length); ft.fft(c);
// 得られた周波数領域の複素数二次配列を,画像に変換する BufferedImage freq = toFreq(c); ImageIO.write(freq, "jpg", new File("freq.jpg"));
// 逆離散フーリエ変換する ft.ifft(c);
// 複素数配列から画像に変換する BufferedImage origin = complexArrayToImage(c); ImageIO.write(origin, "jpg", new File("origin.jpg")); }catch(Exception e){ e.printStackTrace(); } }
/** * nより大きい最少の2の累乗の値を返す. */ private static int smallestPowerOfTwolargerThan(final int n){ int i; for(i = 1; i < n; i *= 2); return i; }
/** * 画像の解像度が2の累乗になるようリサイズする. */ private static BufferedImage resizePowerOfTwo(final BufferedImage src){ assert src != null;
final int srcWidth = src.getWidth(); final int srcHeight = src.getHeight(); final int dstWidth = smallestPowerOfTwolargerThan(srcWidth); final int dstHeight = smallestPowerOfTwolargerThan(srcHeight);
final BufferedImage dst = new BufferedImage( dstWidth, dstHeight, BufferedImage.TYPE_INT_BGR);
for(int i = 0; i < srcHeight; i++){ for(int j = 0; j < srcWidth; j++){ dst.setRGB(j, i, src.getRGB(j, i)); } }
return dst; }
/** * 二次元複素数配列を画像に変換する. * 画像はグレースケール画像になる. */ private static BufferedImage complexArrayToImage(final Complex c[][]){ assert c != null;
final int width = c[0].length; final int height = c.length;
final BufferedImage img = new BufferedImage( width, height, BufferedImage.TYPE_INT_BGR); for(int i = 0; i < height; i++){ for(int j = 0; j < width; j++){ int sig = (int)c[i][j].getRe(); if(sig > 255){ sig = 255; }else if(sig < 0){ sig = 0; } int rgb = sig | sig << 8 | sig << 16; img.setRGB(j, i, rgb); } } return img; }
/** * 画像を二次元複素数配列に変換する. * RGBのうち,B成分のみに着目して変換する. */ private static Complex[][] imageToComplexArray(final BufferedImage img){ assert img != null;
final int width = img.getWidth(); final int height = img.getHeight();
BufferedImage gray = toGrayScale(img);
Complex c[][] = new Complex[height][width]; for(int i = 0; i < height; i++){ for(int j = 0; j < width; j++){ int sig = gray.getRGB(j, i) & 0xFF; c[i][j] = new Complex(sig); } } return c; }
/** * 画像をグレイスケールに変換する. */ private static BufferedImage toGrayScale(final BufferedImage src){ assert src != null;
final int width = src.getWidth(); final int height = src.getHeight(); final BufferedImage dst = new BufferedImage( width, height, BufferedImage.TYPE_INT_BGR);
for(int i = 0; i < height; i++){ for(int j = 0; j < width; j++){ int rgb = src.getRGB(j, i); int r = (rgb & 0x000000FF) >> 0; int g = (rgb & 0x0000FF00) >> 8; int b = (rgb & 0x00FF0000) >> 16; int y = (int)(0.299 * r + 0.587 * g + 0.114 * b); int gray = y | (y << 8) | (y << 16); dst.setRGB(j, i, gray); } } return dst; }
/** * 二次元複素数配列のうち,ノルムの最大値を返す. */ private static double maxNormOf(Complex x[][]){ assert x != null;
double max = x[0][0].norm(); for(int i = 0; i < x.length; i++){ for(int j = 0; j < x[i].length; j++){ double tmp = x[i][j].norm(); if(max < tmp){ max = tmp; } } } return max; }
/** * 複素数配列を画像に変換する. * 複素数のノルムをもとに画像を生成する. * ノルムの大きさをHSB色空間の色相で表現する. * 大きい値は赤,小さい値は青で示す. */ private static BufferedImage toFreq(final Complex[][] x){ assert x != null;
final int width = x[0].length; final int height = x.length;
final BufferedImage dst = new BufferedImage( width, height, BufferedImage.TYPE_INT_BGR);
double max = maxNormOf(x); double k = 1.0 / Math.sqrt(max); for(int i = 0; i < height; i++){ for(int j = 0; j < width; j++){ double rate = k * Math.sqrt(x[i][j].norm() * 128); if(rate > 1){ rate = 1; } int rgb = Color.HSBtoRGB((float)(0.66 * (1 - rate)), 1.0f, 1.0f); dst.setRGB(j, i, rgb); } } return dst; }
}
考察
画像処理を行った結果、低周波成分が大半を占めていることが分かった。
つまり、高周波成分を除去しても元の画像を保つことができる。 逆に低周波成分を除去しすぎると元の画像が大きく崩れる。 この処理を行うことによってどの成分が重要で、どの成分がいらないかの区別を行うことができる。 そして、この技術を応用するとJPEGなどの画像圧縮処理に利用することができる。
先の画像の高周波数成分をゼロにし、逆離散フーリエ変換した。
高周波数成分を大きく除去しても、元の画像が概ね保たれる。
以下に、高周波成分を除去したパワースペクトラムとそれを逆離散フーリエ変換した画像を示す。
- 最終更新:2023-05-15 14:26:55