ハイパス/ローパスフィルタ
前回の続き。
今回はフーリエ変換で得られた周波数スペクトルにハイパスフィルタ(HPF)/ローパスフィルタ(LPF)をかけます。
Flashだとサウンド関連でよく耳にするフィルタかもしれませんが、画像でも基本は同じです。
前回作ったクラスにフィルタ処理等を追加します。
・1D/2D-FFT + HPF/LPF
(追記 2008 12/19) バンドパスフィルタ(BPF)の処理を追加。
(追記 2009 2/20) 空間周波数フィルタのメソッドを少し修正。処理内容は変わってません。
package {
import __AS3__.vec.Vector;
public class FFT{
private var n:int; // データ数
private var bitrev:Vector.<int>; // ビット反転テーブル
private var cstb:Vector.<Number>; // 三角関数テーブル
public static const HPF:String = "high";
public static const LPF:String = "low";
public static const BPF:String = "band"
public function FFT(n:int) {
if(n != 0 && (n & (n-1)) == 0) {
this.n = n;
this.cstb = new Vector.<Number>(n + (n>>2), true);
this.bitrev = new Vector.<int>(n, true);
makeCstb();
makeBitrev();
}
}
// 1D-FFT
public function fft(re:Vector.<Number>, im:Vector.<Number>, inv:Boolean=false):Boolean {
if(!inv) {
return fftCore(re, im, 1);
} else {
if(fftCore(re, im, -1)) {
for(var i:int=0; i<n; i++) {
re[i] /= n;
im[i] /= n;
}
}
}
return true;
}
// 2D-FFT
public function fft2d(re:Vector.<Number>, im:Vector.<Number>, inv:Boolean=false):Boolean {
var tre:Vector.<Number> = new Vector.<Number>(re.length, true);
var tim:Vector.<Number> = new Vector.<Number>(im.length, true);
if(inv) cvtSpectrum(re, im);
// 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];
}
if(!inv) fft(tre, tim);
else fft(tre, tim, true);
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];
}
if(!inv) fft(tre, tim);
else fft(tre, tim, true);
for(var y:int=0; y<n; y++) {
re[x + y*n] = tre[y];
im[x + y*n] = tim[y];
}
}
if(!inv) cvtSpectrum(re, im);
return true;
}
// 空間周波数フィルタ
public function applyFilter(re:Vector.<Number>, im:Vector.<Number>, rad:uint, type:String, bandWidth:uint = 0):void {
var r:int = 0; // 半径
var n2:int = n>>1;
for(var y:int=-n2; y<n2; y++) {
for(var x:int=-n2; x<n2; x++) {
r = Math.sqrt(x*x + y*y);
if((type == HPF && r<rad) || (type == LPF && r>rad) || (type == BPF && (r<rad || r>(rad + bandWidth)))) {
re[x + n2 + (y + n2)*n] = im[x + n2 + (y + n2)*n] = 0;
}
}
}
}
private function fftCore(re:Vector.<Number>, im:Vector.<Number>, sign:int):Boolean {
var h:int, d:int, wr:Number, wi:Number, ik:int, xr:Number, xi:Number, m:int, tmp:Number;
for(var l:int=0; l<n; l++) { // ビット反転
m = bitrev[l];
if(l<m) {
tmp = re[l];
re[l] = re[m];
re[m] = tmp;
tmp = im[l];
im[l] = im[m];
im[m] = tmp;
}
}
for(var k:int=1; k<n; k<<=1) { // バタフライ演算
h = 0;
d = n/(k<<1);
for(var j:int=0; j<k; j++) {
wr = cstb[h + (n>>2)];
wi = sign*cstb[h];
for(var i:int=j; i<n; i+=(k<<1)) {
ik = i+k;
xr = wr*re[ik] + wi*im[ik]
xi = wr*im[ik] - wi*re[ik];
re[ik] = re[i] - xr;
re[i] += xr;
im[ik] = im[i] - xi;
im[i] += xi;
}
h += d;
}
}
return true;
}
// スペクトル位置並び替え
private function cvtSpectrum(re:Vector.<Number>, im:Vector.<Number>):void {
var tmp:Number = 0.0, xn:int = 0, yn:int = 0;
for(var y:int=0; y<(n>>1); y++) {
for(var x:int=0; x<(n>>1); x++) {
xn = x + (n>>1);
yn = y + (n>>1);
tmp = re[x + y*n];
re[x + y*n] = re[xn + yn*n];
re[xn + yn*n] = tmp;
tmp = re[x + yn*n];
re[x + yn*n] = re[xn + y*n];
re[xn + y*n] = tmp;
tmp = im[x + y*n];
im[x + y*n] = im[xn + yn*n];
im[xn + yn*n] = tmp;
tmp = im[x + yn*n];
im[x + yn*n] = im[xn + y*n];
im[xn + y*n] = tmp;
}
}
}
// 三角関数テーブル作成
private function makeCstb():void {
var n2:int = n>>1, n4:int = n>>2, n8:int = n>>3;
var t:Number = Math.sin(Math.PI/n);
var dc:Number = 2*t*t;
var ds:Number = Math.sqrt(dc*(2 - dc));
var c:Number = cstb[n4] = 1;
var s:Number = cstb[0] = 0;
t = 2*dc;
for(var i:int=1; i<n8; i++) {
c -= dc;
dc += t*c;
s += ds;
ds -= t*s;
cstb[i] = s;
cstb[n4 - i] = c;
}
if(n8 != 0) cstb[n8] = Math.sqrt(0.5);
for(var j:int=0; j<n4; j++) cstb[n2 - j] = cstb[j];
for(var k:int=0; k<(n2 + n4); k++) cstb[k + n2] = -cstb[k];
}
// ビット反転テーブル作成
private function makeBitrev():void {
var i:int = 0, j:int = 0, k:int = 0;
bitrev[0] = 0;
while(++i<n) {
k= n >> 1;
while(k<=j) {
j -= k;
k >>= 1 ;
}
j += k;
bitrev[i] = j;
}
}
}
}
・ハイパス/ローパスフィルタのテスト
package {
import __AS3__.vec.Vector;
import flash.display.*;
import flash.text.*;
import flash.utils.getTimer;
[SWF(width = "768", height = "300", backgroundColor = "#000000", frameRate = "12")]
public class FFTImg extends Sprite {
[Embed(source = './assets/lenagray.jpg')]
private var EmbedImage:Class;
private var srcBmp:Bitmap;
private var fftBmp:Bitmap;
private var dstBmp:Bitmap;
private var imgData:BitmapData;
private var fftTest:FFT;
private var re:Vector.<Number>;
private var im:Vector.<Number>;
private var retmp:Vector.<Number>;
private var max:Number = 0.0;
private var data:Number = 0.0;
private var w:int = 0;
private var h:int = 0;
private var tf:TypingTextField;
private var fmt:TextFormat;
public function FFTImg() {
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();
retmp = 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();
// 2D-FFT
fftTest.fft2d(re, im);
var te:int = getTimer();
trace("time: "+(te-ts)+" ms");
// 空間周波数フィルタ(ハイパス)
fftTest.applyFilter(re, im, w/16, FFT.HPF);
// スペクトルの振幅を計算
for(var i:int=0; i<w*h; i++) {
retmp[i] = Math.log(Math.sqrt(re[i]*re[i] + im[i]*im[i]));
if(retmp[i] > max) max = retmp[i];
}
// 8ビットに正規化 (表示用)
for(var i:int=0; i<w*h; i++) {
retmp[i] = retmp[i]/max*255;
}
// 振幅画像の生成
for(var y:int=0; y<h; y++) {
for(var x:int=0; x<w; x++) {
data = retmp[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);
}
}
fftBmp = new Bitmap(imgData.clone());
// 逆2D-FFT
fftTest.fft2d(re, im, true);
// 8ビットに正規化 (表示用)
max = 0.0;
for(var i:int=0; i<w*h; i++) {
if(re[i] > max) max = re[i];
}
for(var i:int=0; i<w*h; i++) {
re[i] = re[i]/max*255;
}
// 出力画像の生成
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;
color = data<<16 | data<<8 | data;
imgData.setPixel(x, y, color);
}
}
dstBmp = new Bitmap(imgData);
addChild(srcBmp);
addChild(fftBmp);
addChild(dstBmp);
fftBmp.x += w;
dstBmp.x += w*2;
}
}
}
・出力結果 (クリックで拡大)
ハイパスフィルタは高域の周波数成分のみ通過させます。
これによって緩やかな変動成分は除去され、エッジ部分だけが残る画像が得られます。
一方、ローパスフィルタは低域の周波数成分のみ通過させます。
こちらは緩やかな変動成分だけが残り、全体的にボケた画像が得られます。
フィルタ適用後のスペクトル画像を見れば視覚的にもわかりやすいかと。
画像の周波数成分が扱えるとできることがさらに広がります。
・関連記事
二次元離散フーリエ変換 – AS3.0
離散フーリエ変換 – AS3.0
Tags: actionscript, cv/im
Comments
No comments so far.