Home > Tags > cv/im
cv/im
位相画像
- 2008-12-20 (Sat)
- tech/study
引き続き、空間周波数領域での画像処理について。
前回、前々回と扱っていたスペクトル画像は周波数成分の振幅特性を表したものでした。
今回は画像の位相特性に注目し、位相画像を生成します。
(位相画像は位相限定画像とも呼ばれる)
FFT演算用クラスはこれまで利用してきたものと同じものです。
(ここではスペクトル位置の並び替え部分は不要なのでその部分の処理は省略できます)
・位相画像の生成
package {
import __AS3__.vec.Vector;
import flash.display.*;
public class FFTImg extends Sprite {
[Embed(source = './assets/gray.jpg')]
private var EmbedImage:Class;
private var srcBmp: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 tmp: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 bias:int = 0;
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();
tmp = 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;
}
}
// 2D-FFT
fftTest.fft2d(re, im);
// スペクトルの振幅で正規化
for(var i:int=0; i<w*h; i++) {
tmp[i] = Math.sqrt(re[i]*re[i] + im[i]*im[i]);
if(tmp[i] > max) max = tmp[i];
re[i] /= tmp[i];
im[i] /= tmp[i];
}
// 逆2D-FFT
fftTest.fft2d(re, im, true);
// 8ビットに正規化 (表示用)
max = 0.0;
bias = 32;
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 + bias;
}
// 出力画像の生成
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(dstBmp);
dstBmp.x += w;
}
}
}
・出力結果
スペクトルを振幅で正規化することにより画像の位相成分のみを取り出しています。
また、一般的にこの位相画像はエッジ画像と似た画像になることが知られています。
これは高周波成分に含まれているエッジ情報が、正規化することによって強調されるためです。
このように位相成分には照明変化にロバストな構造的特徴量が多く含まれるため、
位相成分のみを用いてパターン認識を行う手法もいくつか提案されています。
- Comments : 0
- Trackbacks : 0
ハイパス/ローパスフィルタ
- 2008-12-18 (Thu)
- tech/study
前回の続き。
今回はフーリエ変換で得られた周波数スペクトルにハイパスフィルタ(HPF)/ローパスフィルタ(LPF)をかけます。
Flashだとサウンド関連でよく耳にするフィルタかもしれませんが、画像でも基本は同じです。
前回作ったクラスにフィルタ処理等を追加します。
・1D/2D-FFT + HPF/LPF
(追記 12/19) バンドパスフィルタ(BPF)の処理を追加
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;
}
}
}
・出力結果 (クリックで拡大)
ハイパスフィルタは高域の周波数成分のみ通過させます。
これによって緩やかな変動成分は除去され、エッジ部分だけが残る画像が得られます。
一方、ローパスフィルタは低域の周波数成分のみ通過させます。
こちらは緩やかな変動成分だけが残り、全体的にボケた画像が得られます。
フィルタ適用後のスペクトル画像を見れば視覚的にもわかりやすいかと。
画像の周波数成分が扱えるとできることがさらに広がります。
- Comments : 0
- Trackbacks : 0
二次元離散フーリエ変換 - AS3.0
- 2008-12-15 (Mon)
- tech/study
前回の続き。今回は画像に対するフーリエ変換です。
前回は一次元の話でしたが、画像の場合は平面なので水平/垂直方向の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;
}
入力画像を二次元フーリエ変換して振幅特性を輝度とするスペクトル画像を生成します。
・スペクトル画像の生成
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
- Comments : 0
- Trackbacks : 0
GPU
- 2008-12-14 (Sun)
- diary | tech/study
表題に”実時間”という文字を入れたいので高速化に励む日々。
そんな時 ”それGPUで” のお言葉。
良いアルゴリズムをひねり出す頭がないので機械に頼ります。
ただ、今度は別のことで頭をひねることになりそうですが；
ついこのあいだGeForce 9600 GTに変えたので、デバイスクエリとバンド幅のテスト。
・環境
Core2Quad Q9650 3.00GHz, DDR2-800 2GB, GeForce 9600 GT, PCI Express 2.0
// Device Query
There is 1 device supporting CUDA
Device 0: "GeForce 9600 GT"
Major revision number: 1
Minor revision number: 1
Total amount of global memory: 536543232 bytes
Number of multiprocessors: 8
Number of cores: 64
Total amount of constant memory: 65536 bytes
Total amount of shared memory per block: 16384 bytes
Total number of registers available per block: 8192
Warp size: 32
Maximum number of threads per block: 512
Maximum sizes of each dimension of a block: 512 x 512 x 64
Maximum sizes of each dimension of a grid: 65535 x 65535 x 1
Maximum memory pitch: 262144 bytes
Texture alignment: 256 bytes
Clock rate: 1.60 GHz
Concurrent copy and execution: Yes
Test PASSED
// Bandwidth Test
Running on......
device 0:GeForce 9600 GT
Quick Mode
Host to Device Bandwidth for Pageable memory
.
Transfer Size (Bytes) Bandwidth(MB/s)
33554432 2314.8
Quick Mode
Device to Host Bandwidth for Pageable memory
.
Transfer Size (Bytes) Bandwidth(MB/s)
33554432 2148.7
Quick Mode
Device to Device Bandwidth
.
Transfer Size (Bytes) Bandwidth(MB/s)
33554432 38671.6
&&&& Test PASSED
CUDA2.1ならVS2008でビルドできるらしいんですが、まだベータ版のようです。
とりあえずCUDA2.0+VS2005でしばらくやってみることにします。
- Comments : 0
- Trackbacks : 0
離散フーリエ変換 - AS3.0
- 2008-12-13 (Sat)
- tech/study
“画像” に対する離散フーリエ変換について。 ”音” は扱わないので注意。
また、一度で書くには分量が多いので何回かに分けて書こうと思います。
・離散フーリエ変換 (DFT)
フーリエ変換についてはFlasherの間でも音関連でよく話題に出ているようなので理論的な説明は省きます。
画像の場合も音の周波数解析と基本は同じ。
ただ、音と違って画像は二次元なので空間周波数領域への変換を行うことになります。
(音の信号 x(t)は一次元、画像の信号 f(x,y) は二次元になる)
実際の計算ではX軸方向にフーリエ変換した後、データを転置してY軸方向にもう一度フーリエ変換します。
また、高速に計算を行うために実装上では高速フーリエ変換(FFT)を用いるのが一般的です。
とりあえずは一次元のFFT(1D-FFT)から始めます。
ActionScriptでの例は以下のようになります。
・1D-FFT
package {
import __AS3__.vec.Vector;
public class FFT{
private var n:int; // データ数
private var bitrev:Vector.<int>; // ビット反転テーブル
private var cstb:Vector.<Number>; // 三角関数テーブル
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();
}
}
// FFT
public function fft(re:Vector.<Number>, im:Vector.<Number>):Boolean {
return fftCore(re, im, 1);
}
//逆FFT
public function ifft(re:Vector.<Number>, im:Vector.<Number>):Boolean {
if(fftCore(re, im, -1)){
for(var i:int=0; i<n; i++){
re[i] /= n;
im[i] /= n;
}
}
return true;
}
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;
for(var i:int=0; i<n; i++){ // ビット反転
var j:int = bitrev[i];
if(i<j){
var tmp:Number = re[i];
re[i] = re[j];
re[j] = tmp;
tmp = im[i];
im[i] = im[j];
im[j] = 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 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));
t = 2*dc;
var c:Number = cstb[n4] = 1;
var s:Number = cstb[0] = 0;
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.Sprite;
import flash.utils.getTimer;
public class Test extends Sprite {
private const N:int = 32;
private var fftTest:FFT;
public function Test() {
fftTest = new FFT(N);
var x1:Vector.<Number> = new Vector.<Number>(N, true);
var y1:Vector.<Number> = x1.concat();
var x2:Vector.<Number> = x1.concat();
var y2:Vector.<Number> = x1.concat();
var x3:Vector.<Number> = x1.concat();
var y3:Vector.<Number> = x1.concat();
for(var i:int=0; i<N; i++){
x1[i] = x2[i] = 6*Math.cos(6*Math.PI*i/N) + 4*Math.sin(18*Math.PI*i/N);
}
var ts:int = getTimer();
fftTest.fft(x2, y2); // FFT
var te:int = getTimer();
trace("time FFT: "+(te-ts)+" ms");
for(var i:int=0; i<N; i++){
x3[i] = x2[i];
y3[i] = y2[i];
}
fftTest.ifft(x3, y3); // 逆FFT
trace("データ数: (元データ) (FFT) (逆FFT)");
for(var i:int=0; i<N; i++){
trace(i+": ("+x1[i]+", "+y1[i]+") ("+x2[i]+", "+y2[i]+") ("+x3[i]+", " +y3[i]+")");
}
}
}
}
・出力
誤差が出るので逆FFTしたときに全く同じ値には戻りません。
(誤差：10のマイナス14~15乗ほど)
time FFT: 1 ms データ数: (元データ) (FFT) (逆FFT) 0: (6, 0) (1.3100631690576847e-14, 0) (6, 0) 1: (8.911958795428193, 0) (-3.148595663620103e-15, -1.634119351805764e-15) (8.911958795428195, 2.220446049250315e-16) 2: (0.7653668647301803, 0) (-3.582018523007779e-14, 1.0273603512156518e-14) (0.76536686473018, 5.551115123125783e-17) 3: (-4.496420381306951, 0) (96.00000000000003, -4.440892098500626e-14) (-4.496420381306951, -4.440892098500624e-16) 4: (-1.4142135623730954, 0) (2.074753706171886e-14, 3.110202491371294e-14) (-1.4142135623730956, 0) 5: (-3.6624307503409734, 0) (-1.1638332535152583e-14, 4.013150866834269e-16) (-3.6624307503409748, -2.2204460492503091e-16) 6: (-9.238795325112866, 0) (2.0928516568038727e-14, -1.362553389586755e-14) (-9.238795325112866, 0) 7: (-4.113782686182125, 0) (2.6645352591003757e-15, 2.1316282072803006e-14) (-4.113782686182125, 4.440892098500626e-16) 8: (3.999999999999999, 0) (3.3306690738754696e-15, 3.419486915845482e-14) (3.999999999999999, 0) 9: (2.5530601100531034, 0) (-6.838973831690964e-14, -63.99999999999999) (2.553060110053105, -5.551115123125785e-16) 10: (1.8477590650225717, 0) (-9.893236143066244e-15, -2.623249277195817e-14) (1.8477590650225717, 0) 11: (8.106992614497795, 0) (-9.677949537650423e-15, -1.0256825949718076e-14) (8.106992614497795, -4.440892098500628e-16) 12: (7.071067811865481, 0) (-1.6306644963218234e-14, -2.796183999634539e-14) (7.07106781186548, 0) 13: (-2.155336517113407, 0) (-2.1316282072803006e-14, -2.6645352591003757e-14) (-2.1553365171134082, -7.7715611723761e-16) 14: (-3.826834323650897, 0) (3.3666689002106555e-14, -5.5624060545941617e-14) (-3.826834323650898, 6.661338147750939e-16) 15: (-1.0656765522023512, 0) (5.4662944006227367e-14, -1.584497406700777e-14) (-1.065676552202351, -4.718447854656915e-16) 16: (-5.999999999999996, 0) (2.4646951146678475e-14, 0) (-5.999999999999995, 0) 17: (-8.911958795428195, 0) (5.4662944006227367e-14, 1.584497406700777e-14) (-8.911958795428198, -2.2204460492503114e-16) 18: (-0.7653668647301939, 0) (3.3666689002106555e-14, 5.5624060545941617e-14) (-0.7653668647301936, -5.551115123125783e-17) 19: (4.496420381306947, 0) (-2.1316282072803006e-14, 2.6645352591003757e-14) (4.496420381306947, 4.440892098500628e-16) 20: (1.414213562373091, 0) (-1.6306644963218234e-14, 2.796183999634539e-14) (1.4142135623730911, 0) 21: (3.6624307503409756, 0) (-9.677949537650423e-15, 1.0256825949718076e-14) (3.6624307503409765, 2.220446049250317e-16) 22: (9.23879532511287, 0) (-9.893236143066244e-15, 2.623249277195817e-14) (9.23879532511287, 0) 23: (4.113782686182145, 0) (-6.838973831690964e-14, 63.99999999999999) (4.113782686182146, -4.440892098500626e-16) 24: (-3.999999999999997, 0) (3.3306690738754696e-15, -3.419486915845482e-14) (-3.9999999999999973, 0) 25: (-2.5530601100531056, 0) (2.6645352591003757e-15, -2.1316282072803006e-14) (-2.553060110053107, 5.551115123125781e-16) 26: (-1.847759065022569, 0) (2.0928516568038727e-14, 1.362553389586755e-14) (-1.847759065022569, 0) 27: (-8.106992614497791, 0) (-1.1638332535152583e-14, -4.013150866834269e-16) (-8.106992614497791, 4.440892098500624e-16) 28: (-7.071067811865472, 0) (2.074753706171886e-14, -3.110202491371294e-14) (-7.071067811865471, 0) 29: (2.1553365171133887, 0) (96.00000000000003, 4.440892098500626e-14) (2.1553365171133905, 7.771561172376092e-16) 30: (3.826834323650912, 0) (-3.582018523007779e-14, -1.0273603512156518e-14) (3.826834323650913, -6.661338147750939e-16) 31: (1.065676552202345, 0) (-3.148595663620103e-15, 1.634119351805764e-15) (1.0656765522023448, 4.718447854656915e-16)
以前研究で医用画像を扱っていたことがあって、その時に周波数解析の基礎を学びました。
でもその後、画像データの数が確保できずにその研究を断念することに；；
(医用画像はそのまま患者の個人情報になるので取り扱いに注意しなければならない。
そのため医療機関のサポート等がないと、研究で使えるデータを十分確保できないなどの理由)
Flashで医用画像を扱うつもりはもちろんありませんが、基礎だけでも形に残しておこうかと。
画像の周波数解析自体は、医用画像に限らず幅広い分野で用いられている基礎的な処理になります。
今回は一次元のFFTでしたが、次はこれを画像に適用するために二次元のFFTを考えたいと思います。
- Comments : 0
- Trackbacks : 0
Home > Tags > cv/im
