チーム20A3
課題名
画像の周波数成分をFFTによる可視化&アベレージフィルタリングを行うプログラム
研究者名
3年14組14番 Makoto Kageyama
3年14組16番 Tastuki Kawamoto
概要
1,高速フーリエ変換を用いて画像の周波数成分を可視化する。
以下のように処理をする。
・画像の白黒処理
・高速フーリエ変換
・変換した値の絶対値に色を対応させて強弱がわかるようにして出力
元の画像
白黒処理後の画像
高速フーリエ変換後の値は、複素数であるため大小の判断のために絶対値を取りその値を比較し、
値に対して色付けし周波数成分の分布を表現したのが以下の画像
(中央->隅:高周波数->低周波数、赤->青:周波数成分が多い->少ない)
2,アベレージフィルタの適用
以下のように処理する。
・3タップのアベレージフィルタを100回ごとにかける
3タップのアベレージフィルタを100回適用
3タップのアベレージフィルタを200回適用
3タップのアベレージフィルタを300回適用
3タップのアベレージフィルタを400回適用
3タップのアベレージフィルタを500回適用
3タップのアベレージフィルタを600回適用
3タップのアベレージフィルタを700回適用
3タップのアベレージフィルタを800回適用
3タップのアベレージフィルタを900回適用
3タップのアベレージフィルタを1000回適用
アベレージフィルタを1000回程度かけることでほぼ完全に単色になったことがわかる。
プログラム
15C2,17C3を参考に実装した。
・複素数のクラス。
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); } }
}
・高速フーリエ変換を2次元で計算
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"));
// アベレージフィルタをかける for(int i=0;i<1000;i++){ src = toAvarageFilter(src, 5); if(i%100==0){ ImageIO.write(toAvarageFilter(src,10),"jpg",new File("./output2/aveOut"+((i/100)+1)+"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; }
/** *アベレージフィルタをかける. **/ private static BufferedImage toAvarageFilter(final BufferedImage src,int n){ final int width=src.getWidth(); final int height=src.getHeight(); final BufferedImage dest = new BufferedImage( width,height,BufferedImage.TYPE_INT_BGR);
if(n%2==1)n++; for(int i=0;i<height;i++){ for(int j=0;j<width;j++){ int aveR=0,aveG=0,aveB=0; for(int k=-n/2;k<n/2;k++){ int rgb; if(j+k>=width| |j+k<0){ rgb=src.getRGB(j,i); }else{ rgb=src.getRGB(j+k,i); } int r=(rgb&0x000000FF)>>0; int g=(rgb&0x0000FF00)>>8; int b=(rgb&0x00FF0000)>>16; aveR+=r; aveG+=g; aveB+=b; } for(int k=0-n/2;k<n/2;k++){ int rgb; if(i+k>=height| |i+k<0){ rgb=src.getRGB(j,i); }else{ rgb=src.getRGB(j,i+k); } int r=(rgb&0x000000FF)>>0; int g=(rgb&0x0000FF00)>>8; int b=(rgb&0x00FF0000)>>16; aveR+=r; aveG+=g; aveB+=b; } aveR/=2*n; aveG/=2*n; aveB/=2*n; int rgb=aveR|(aveG<<8)|(aveB<<16); dest.setRGB(j,i,rgb); } } return dest; }
}
*ImageFFT.javaにアベレージフィルタを追加で実装した。
考察
・周波数成分の画像より、大部分は低周波数成分であることがわかった。
よって、今回の画像のデータ圧縮の場合にはローパスフィルターを用いて高周波数帯を切ることでデータを圧縮することができると考察される。
・アベレージフィルタをかけた画像より、3タップのアベレージフィルタであれば100回程度までは画像を識別できるので、ぼかし処理として使えるであろうと考察できる。
- 最終更新:2020-06-01 16:12:00