前回の続き。
今回はフーリエ変換で得られた周波数スペクトルにハイパスフィルタ(HPF)/ローパスフィルタ(LPF)をかけます。
Flashだとサウンド関連でよく耳にするフィルタかもしれませんが、画像でも基本は同じです。
前回作ったクラスにフィルタ処理等を追加します。
・1D/2D-FFT + HPF/LPF
(追記 2008 12/19) バンドパスフィルタ(BPF)の処理を追加。
(追記 2009 2/20) 空間周波数フィルタのメソッドを少し修正。処理内容は変わってません。
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 |
package { import __AS3__.vec.Vector; /* * Fast Fourier Transform implementation in pure AS3 */ 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):void { if(!inv) { fftCore(re, im, 1); }else { fftCore(re, im, -1); for(var i:int=0; i<n; i++) { re[i] /= n; im[i] /= n; } } } // 2D-FFT public function fft2d(re:Vector.<Number>, im:Vector.<Number>, inv:Boolean=false):void { var tre:Vector.<Number> = new Vector.<Number>(re.length, true); var tim:Vector.<Number> = new Vector.<Number>(im.length, true); var i:uint; if(inv) swapQuadrant(re, im); // x-axis for(var y:int=0; y<n; y++) { i = y*n; for(var x1:int=0; x1<n; x1++) { tre[x1] = re[x1 + i]; tim[x1] = im[x1 + i]; } if(!inv) fft(tre, tim); else fft(tre, tim, true); for(var x2:int=0; x2<n; x2++) { re[x2 + i] = tre[x2]; im[x2 + i] = tim[x2]; } } // y-axis for(var x:int=0; x<n; x++) { for(var y1:int=0; y1<n; y1++) { i = x + y1*n; tre[y1] = re[i]; tim[y1] = im[i]; } if(!inv) fft(tre, tim); else fft(tre, tim, true); for(var y2:int=0; y2<n; y2++) { i = x + y2*n; re[i] = tre[y2]; im[i] = tim[y2]; } } if(!inv) swapQuadrant(re, im); } // spatial frequency filtering. public function applyFilter(re:Vector.<Number>, im:Vector.<Number>, rad:uint, type:String, bandWidth:uint = 0):void { var r:int = 0; // radius var n2:int = n>>1; var i:int, ptr:int; for(var y:int=-n2; y<n2; y++) { i = n2 + (y + n2)*n; for(var x:int=-n2; x<n2; x++) { r = Math.sqrt(x*x + y*y); ptr = x + i; if((type == HPF && r < rad) || (type == LPF && r > rad) || (type == BPF && (r < rad || r > (rad + bandWidth)))) { re[ptr] = im[ptr] = 0; } } } } // Fast Fourier Transform core operation. private function fftCore(re:Vector.<Number>, im:Vector.<Number>, sign:int):void { var h:int, d:int, wr:Number, wi:Number, ik:int, xr:Number, xi:Number, m:int, tmp:Number; // bit reversal 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; } } // butterfly operation 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; } } } // swap quadrant. private function swapQuadrant(re:Vector.<Number>, im:Vector.<Number>):void { var tmp:Number, xn:int, yn:int, i:int, j:int, k:int, l:int; var len:int = n>>1; for(var y:int=0; y<len; y++) { yn = y + len; for(var x:int=0; x<len; x++) { xn = x + len; i = x + y*n; j = xn + yn*n; k = x + yn*n; l = xn + y*n; tmp = re[i]; re[i] = re[j]; re[j] = tmp; tmp = re[k]; re[k] = re[l]; re[l] = tmp; tmp = im[i]; im[i] = im[j]; im[j] = tmp; tmp = im[k]; im[k] = im[l]; im[l] = tmp; } } } // make table of trigonometric function. 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]; } // make table of bit reversal. 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; } } } } |
・ハイパス/ローパスフィルタのテスト
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 |
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