二次元離散フーリエ変換 – AS3.0

前回の続き。今回は画像に対するフーリエ変換です。
前回は一次元の話でしたが、画像の場合は平面なので水平/垂直方向の2つの周波数を持つことになります。
実際に画像に対して二次元フーリエ変換を行うには、
　x軸方向に一次元フーリエ変換　→　y軸方向に一次元フーリエ変換
という手順で処理します。
ここでは、前回作ったクラスを少し拡張します。
・2D-FFT

public function fft2d(re:Vector.<Number>, im:Vector.<Number>):Boolean {
	var tre:Vector.<Number> = new Vector.<Number>(re.length, true);
	var tim:Vector.<Number> = new Vector.<Number>(im.length, true);
	// x軸方向のFFT
	for(var y:int=0; y<n; y++) {
		for(var x:int=0; x<n; x++) {
			tre[x] = re[x + y*n];
			tim[x] = im[x + y*n];
		}
		fft(tre, tim);
		for(var x:int=0; x<n; x++) {
			re[x + y*n] = tre[x];
			im[x + y*n] = tim[x];
		}
	}
	// y軸方向のFFT
	for(var x:int=0; x<n; x++) {
		for(var y:int=0; y<n; y++) {
			tre[y] = re[x + y*n];
			tim[y] = im[x + y*n];
		}
		fft(tre, tim);
		for(var y:int=0; y<n; y++) {
			re[x + y*n] = tre[y];
			im[x + y*n] = tim[y];
		}
	}
	return true;
}

public function fft2d(re:Vector.<Number>, im:Vector.<Number>):Boolean {

var tre:Vector.<Number> = new Vector.<Number>(re.length, true);

var tim:Vector.<Number> = new Vector.<Number>(im.length, true);

// x軸方向のFFT

for(var y:int=0; y<n; y++) {

for(var x:int=0; x<n; x++) {

tre[x] = re[x + y*n];

tim[x] = im[x + y*n];

}

fft(tre, tim);

for(var x:int=0; x<n; x++) {

re[x + y*n] = tre[x];

im[x + y*n] = tim[x];

}

// y軸方向のFFT

for(var x:int=0; x<n; x++) {

for(var y:int=0; y<n; y++) {

tre[y] = re[x + y*n];

tim[y] = im[x + y*n];

}

fft(tre, tim);

for(var y:int=0; y<n; y++) {

re[x + y*n] = tre[y];

im[x + y*n] = tim[y];

}

return true;

}

入力画像を二次元フーリエ変換して振幅特性を輝度とするスペクトル画像を生成します。
・スペクトル画像の生成

package {
	import __AS3__.vec.Vector;
	import flash.display.*;
	import flash.geom.Point;
	import flash.geom.Rectangle;
	import flash.utils.getTimer;

	public class Test extends Sprite {
		[Embed(source = './assets/gray.jpg')] // グレースケール画像
		private var EmbedImage:Class;

		private var srcBmp:Bitmap;
		private var imgData:BitmapData;
		private var buffData:BitmapData;
		private var fftTest:FFT;
		private var re:Vector.<Number>;
		private var im:Vector.<Number>;
		private var max:Number = 0.0;
		private var data:Number = 0.0;
		private var w:int = 0;
		private var h:int = 0;

		public function Test() {
			srcBmp = new EmbedImage() as Bitmap;
			imgData = srcBmp.bitmapData.clone();
			w = imgData.width;
			h = imgData.height;
			fftTest = new FFT(w);
			re = new Vector.<Number>(w*h, true);
			im = re.concat();
			// 画像データの直列化
			for(var y:int=0; y<h; y++) {
				for(var x:int=0; x<w; x++) {
					re[x + y*w] = imgData.getPixel(x, y) & 0xff;
				}
			}
			var ts:int = getTimer();
			// 二次元離散フーリエ変換
			fftTest.fft2d(re, im);
			var te:int = getTimer();
			trace("time FFT: "+(te-ts)+" ms");
			// スペクトルの振幅を計算
			for(var i:int=0; i<w*h; i++) {
				re[i] = Math.log(Math.sqrt(re[i]*re[i] + im[i]*im[i]));
				if(re[i] > max) max = re[i];
			}
			// 正規化
			for(var i:int=0; i<w*h; i++) {
				re[i] = re[i]*255/max;
			}
			// 振幅画像の生成
			for(var y:int=0; y<h; y++) {
				for(var x:int=0; x<w; x++) {
					data = re[x + y*w];
					if(data>255) data = 255;
					if(data<0) data = 0;
					var color:uint = data<<16 | data<<8 | data;
					imgData.setPixel(x, y, color);
				}
			}
			// 象限の入れ替え
			buffData = new BitmapData(w, h, false);
			buffData.copyPixels(imgData, new Rectangle(w/2, 0, w/2, h/2), new Point(0, h/2)); // 1 -> 3
			buffData.copyPixels(imgData, new Rectangle(0, 0, w/2, h/2), new Point(w/2, h/2)); // 2 -> 4
			buffData.copyPixels(imgData, new Rectangle(0, h/2, w/2, h/2), new Point(w/2, 0)); // 3 -> 1
			buffData.copyPixels(imgData, new Rectangle(w/2, h/2, w/2, h/2), new Point(0, 0)); // 4 -> 2
			imgData = buffData.clone();
			buffData.dispose();
			addChild(new Bitmap(imgData));
		}
	}
}

package {

import __AS3__.vec.Vector;

import flash.display.*;

import flash.geom.Point;

import flash.geom.Rectangle;

import flash.utils.getTimer;

public class Test extends Sprite {

[Embed(source = './assets/gray.jpg')] // グレースケール画像

private var EmbedImage:Class;

private var srcBmp:Bitmap;

private var imgData:BitmapData;

private var buffData:BitmapData;

private var fftTest:FFT;

private var re:Vector.<Number>;

private var im:Vector.<Number>;

private var max:Number = 0.0;

private var data:Number = 0.0;

private var w:int = 0;

private var h:int = 0;

public function Test() {

srcBmp = new EmbedImage() as Bitmap;

imgData = srcBmp.bitmapData.clone();

w = imgData.width;

h = imgData.height;

fftTest = new FFT(w);

re = new Vector.<Number>(w*h, true);

im = re.concat();

// 画像データの直列化

for(var y:int=0; y<h; y++) {

for(var x:int=0; x<w; x++) {

re[x + y*w] = imgData.getPixel(x, y) & 0xff;

}

var ts:int = getTimer();

// 二次元離散フーリエ変換

fftTest.fft2d(re, im);

var te:int = getTimer();

trace("time FFT: "+(te-ts)+" ms");

// スペクトルの振幅を計算

for(var i:int=0; i<w*h; i++) {

re[i] = Math.log(Math.sqrt(re[i]*re[i] + im[i]*im[i]));

if(re[i] > max) max = re[i];

}

// 正規化

for(var i:int=0; i<w*h; i++) {

re[i] = re[i]*255/max;

}

// 振幅画像の生成

for(var y:int=0; y<h; y++) {

for(var x:int=0; x<w; x++) {

data = re[x + y*w];

if(data>255) data = 255;

if(data<0) data = 0;

var color:uint = data<<16 | data<<8 | data;

imgData.setPixel(x, y, color);

}

// 象限の入れ替え

buffData = new BitmapData(w, h, false);

buffData.copyPixels(imgData, new Rectangle(w/2, 0, w/2, h/2), new Point(0, h/2)); // 1 -> 3

buffData.copyPixels(imgData, new Rectangle(0, 0, w/2, h/2), new Point(w/2, h/2)); // 2 -> 4

buffData.copyPixels(imgData, new Rectangle(0, h/2, w/2, h/2), new Point(w/2, 0)); // 3 -> 1

buffData.copyPixels(imgData, new Rectangle(w/2, h/2, w/2, h/2), new Point(0, 0)); // 4 -> 2

imgData = buffData.clone();

buffData.dispose();

addChild(new Bitmap(imgData));

}

一般的に画像では振幅を持つ周波数成分と比べて直流成分が非常に多く(ダイナミックレンジが広い)、
実際の振幅をそのまま表示しようとすると、直流成分以外のデータがほとんどない状態になります。
そこで、振幅の対数をとることで成分間の差を小さくします。
さらに、このまま画像を表示すると原点となる直流成分が画像の端になるので、
それが画像の中心にくるようにデータを象限単位で並び替えます。

これで空間周波数の世界に行けるようになりました。

・関連記事
離散フーリエ変換 – AS3.0

二次元離散フーリエ変換 – AS3.0

あわせて読む:

コメントを残す