核医学機器工学概論３

実空間フィルタ

フィルタ処理、畳込み

フォルダFilter内にあるFilter.exeを起動し、 「LoadFile」ボタンを押してBitmapまたはJPEG画像を選択。
（プログラムフォルダ内にあるsample.jpgを使用）

画像処理において用いられる実空間、空間周波数フィルタ処理とその演算方法である畳込みについて学ぶ「Kernel」の数値（９個のフィルタ行列の要素）を変更し、任意の3×3の実空間フィルタを作成する。
「Filter」ボタンを押し、作成したフィルタによりフィルタ処理を行う。

1

2

画像に平滑化フィルタ処理（smoothing）を行う

処理を行うと以下の様な画像が得られる

34

平滑化フィルタ

実空間での雑音除去フィルタ

フィルタ処理前後で画素値の総和が変わらないようにフィルタ内成分の和が１になっていることに注目。
移動平均フィルタ （成分値が全部同じ）
最も単純な 3x3 スムージングフィルタ。
中央画素の重み付けがない。
輪郭の鮮明さが損なわれる。
荷重平均フィルタ
（成分値が中心部で大きい）
輪郭の鮮明さを維持するために中心部に重みを付けた
3x3, 5x5 スムージングフィルタ。

5解答 ５61void filter1( int x[ ][260], int y[ ][260] ){// Weighted smoothing filter
int
i, j ;

void filter2( int x[ ][260], int y[ ][260] ){// Simple smoothing filter
int
i, j ;

for(j=2;j<=255;j++){ for(i=2;i<=255;i++){for(j=2;j<=255;j++){ for(i=2;i<=255;i++){y[i][j] = x[i-1][j-1] + x[i][j-1]*2 + x[i+1][j-1]+ x[i-1][j]*2 + x[i][j]*4 + x[i+1][j-1]*2+ x[i-1][j+1] + x[i][j+1]*2 + x[i+1][j+1] ;y[i][j] = x[i-1][j-1] + x[i][j-1] + x[i+1][j-1]+ x[i-1][j]+ x[i][j]+ x[i+1][j-1]+ x[i-1][j+1] + x[i][j+1] + x[i+1][j+1] ;}}for(j=2;j<=255;j++){ for(i=2;i<=255;i++){ y[i][j] /= 16; }}}}
for(j=2;j<=255;j++){ for(i=2;i<=255;i++){ y[i][j] /= 9; }}}}輪郭の鮮明さを維持するために中心部に重みを付けたスムージングフィルタ。

最も単純なスムージングフィルタ。
中央画素の重み付けがない。
輪郭の鮮明さが損なわれる。

78void filter3( int x[ ][260], int y[ ][260] ){// Vertical Prewitt filter
int i, j ;void filter4( int x[ ][260], int y[ ][260] ){// Horizontal Sobel filter
int i, j ;for(j=2;j<=255;j++){ for(i=2;i<=255;i++){for(j=2;j<=255;j++){ for(i=2;i<=255;i++){y[i][j] = x[i-1][j-1]*-1 + x[i][j-1]*-1 + x[i+1][j-1]*-1+ x[i-1][j]*0+ x[i][j]*0+ x[i+1][j-1]*0+ x[i-1][j+1]+ x[i][j+1]+ x[i+1][j+1] ;y[i][j] = x[i-1][j-1]*-1 + x[i][j-1]*0 + x[i+1][j-1]+ x[i-1][j]*-2 + x[i][j]*0+ x[i+1][j-1]*2+ x[i-1][j+1]*-1 + x[i][j+1]*0 + x[i+1][j+1] ;

}}

}}

}

}

１次微分（差分）エッジ検出フィルタ
上下方向に画素値が大きく変化する部位（輪郭）の抽出フィルタ

Sobel フィルタ （水平方向）
１次微分（差分）エッジ検出フィルタ
左右方向に画素値が大きく変化する部位（輪郭）の抽出フィルタに中央部の重み付けを加えている。

910Prewitt フィルタ （垂直方向）

void filter5( int x[ ][260], int y[ ][260] ){// Unsharp masking filter
int i, j ;Smoothing.exe
実行結果

for(j=2;j<=255;j++){ for(i=2;i<=255;i++){Bone6noise画像を選択。

y[i][j] = x[i-1][j-1]*0 + x[i][j-1]*-1 + x[i+1][j-1]*0+ x[i-1][j]*-1 + x[i][j]*5+ x[i+1][j-1]*-1+ x[i-1][j+1]*0 + x[i][j+1]*-1 + x[i+1][j+1]*0 ;}}}４近傍 先鋭化フィルタ
アンシャープマスキング。
ぼやけた部位をマスクする。
上下および左右方向に画素値が大きく変化する部位（輪郭）の強調
11122

画像に 中央値フィルタ処理（median）を行う

median.c 実行結果

1314

void median( int x[ ][260], int y[ ][260] ){int i, j, k, mi, mj, ni, nj, yn, min, med[10];for( j=2; j<=255; j++){ for( i=2; i<=255; i++){median.c
メディアン
フィルタ

k = 1;
for( mi=-1; mi<=1; mi++){ for( mj=-1; mj<=1; mj++){med[k] = x[i+mi][j+mj] ; k++ ;}}k = 1；
for( mi=-1; mi<=1; mi++){ for( mj=-1; mj<=1; mj++) {}}med[ k ] = x[ i+mi ][ j+mj ] ; k++ ;座標 [ i ][ j ] と、その近傍の画素カウント値を
配列 med[1] ～ med[9] に代入している。

for( ni=1; ni<=8; ni++){ for( nj=8; nj >= ni; nj - -){if( med[nj+1] < med[nj] ){}}}}}}

min=med[nj+1]; med[nj+1]=med[nj]; med[nj]=min;

C言語では、配列数は ０番目から数えている のでmed[9]; と宣言すると、med[0] から med[8] の９個の配列要素だけが宣言され、med[9] は宣言されない。

y[ i ][ j ] = med[ 5 ];15

16

for ( ni=1; ni<=8; ni++) { for ( nj=8; nj >= ni; nj - -) {if ( med[nj+1] < med[nj] ){}}

配列 med は、 int med[10]; と宣言している。
int med[9]; でも良さそうに思えるが、}はじめに ni が 1 の 値で nj の ループに入る。
nj の 値は 8 から１ までの値に変化しながら ループを回る。

min=med[nj+1]; med[nj+1]=med[nj]; med[nj]=min;

最初は nj が 8 なので、if 文の中身を具体的に書くとif ( med[ 9 ] < med[ 8 ] ){med[1] から med[9] の値を並べ替えている （ソート sort）。

min=med[ 9 ]; med[ 9 ]=med[ 8 ]; med[ 8 ]=min;

この for ループが終了すると、 med[1] から med[9] の中の}最小値が med[1] に、最大値が med[9] に代入されている。

もし med[9] が med[8] より小さければ med[9] と med[8]の値を入れ替える。 変数の入れ替え作業を swap という。
swap作業には 第三の変数（ここでは min）が必要になる。

バブルソートのアルゴリズム

Bubble sort algorithm

コップ底の泡が、軽いものから先に順番に昇っていくような様子をプログラムで表現している。ソート法のアルゴリズム

ni が 1 の 状態で、nj が 8 から １ まで変化し終わると、med[1]の値は、med[1]からmed[9]の最小値になっている。

としては最も非効率的な方法だが単純で理解しやすい。

小さい値が泡のように上に昇っていくので バブルソート法。

17183次に ni が ２ の 値で nj の ループに入る。

smoothing filter は、輪郭が不明瞭化し画像がぼやける。

nj の 値は 8 から ２ までの値に変化しながら ループを回る。

median filter は smoothing filter より輪郭が保たれる。

（ med[1] は最小値が確定なので ループに入る必要がない。）

median filter は、微細なノイズを除去する効果がある。

ni が 2 の 状態で、nj が 8 から 2 まで変化し終わると、med[2] の値は、med[2] から med[９] の中の最小値になっている。

bone6noise は、10ピクセルおきに画素値 100 を入れた像。

次に ni が ３ の 値で nj の ループに入る。
nj の 値は 8 から ３ までの値に変化しながら ループを回る。

median では、このノイズを完全に消すことができる。
smoothing ではノイズを完全には消せない。

（ med[1]、med[2] は １番目、２番目に小さい値が確定なのでループに入る必要がない。）

同様の操作を ni が 8 になるまで繰り返すと、med[1] に最小値、med[2] に２番目に小さい値、med[3]に３番目に小さい値、 ・ ・、med[9]に9番目に小さい値が入る。
中央値は med[5] なので、 y[ i ][ j ] = med[5] ; となる。

1920

フーリエ変換 （Fourier Transform）は、波形データなどが含む周波数成分（スペクトル）を分析する方法。

フーリエ変換 （Fourier Transform）の基本式は、波形データが 連続データであることに基づいている。

交流雑音（ハム）を含む心電図をフーリエ変換すると、５０Hzの周波数成分が多く含まれていることが判る。
心電図の １次元フーリエ変換
心電図波形成分５０Hzの商用交流雑音成分
（Hum）

2122

高速フーリエ変換 （ FFT : Fast Fourier Transform ）フーリエ変換を高速に計算するアルゴリズム。
1942年に Danielson と Lanczos が発明。
プログラムによるフーリエ変換は、波形データが離散的なので、基本式に現れる e（２πｆ t ）i の項を、Ｗ ｎｋ と変形して巧みに解いている。
高速逆フーリエ変換（ IFFT : Inverse FFT ）も ほとんど同じアルゴリズムなので、両方可能な関数が作れる。
FFT 、 IFFT におけるデータの制限としては、データ数が ２の階乗でなければならない。

23244

FFT

Fast Fourier Transform

FFT.exeを起動し、open fileから同フォルダ内にあるboneを選択する。
※bone：99mTc‐MDPのSPECT画像。骨シンチ。

空間周波数フィルタ

フォルダButterworth内にあるButterworth.exeを起動し、 「Load File」ボタンを押して、プログラムフォルダ内にあるboneファイルを読み込む。

FFTをクリックして2次元フーリエ変換を行い、replaceをクリックして、パワースペクトル画像を再配置する。

「display sin cos component」ボタンを押し、原画像のフーリエ変換を行い、のsin成分とcos成分の画像を表示させる。

discoverをクリックすると、逆フーリエ変換された画像が表示される。

「FFT spectrum」ボタンを押し、パワースペクトルを表示する。

coverをクリックすると、右側の画像が消える。
25

26画像を フーリエ変換する。

FFT.c

２次元 FFT の出力である ２次元周波数分布データは、４隅が原点 ( 0,0 )（= 直流成分）になっているのが面倒な点である。
中心が原点になるように各々 1/4 領域を反転する。

2728

原画像 img を x[ ] に入れて２次元 FFT を行うと、2次元周波数分布のsin成分が x[ ] に、cos成分がy[ ] に出力される。
x[ ] と y[ ] を、４隅が原点 ( 0,0 ) になっている状態で ２次元逆フーリエ変換 ( IFFT ) の関数に入力すると、x[ ] に 元に戻った像が出力される。

Butterworth.c

2930

周波数空間でフィルタ処理を行う

バターワース
フィルタ

5結果として、以下の様な画像が得られる。

テキストボックス内に0.5～0.01の任意の遮断周波数を設定する。
「Butterworth filter」ボタンを押し、設定した遮断周波数からバターワースフィルタを作成し、空間周波数フィルタ処理を行う。
「IFFT」ボタンを押し、逆フーリエ変換を行い、フィルタ処理後の画像を表示する。

出力画像

3132

画像 を 周波数空間でフィルタ処理する。

１次元の Butterworth filter を ２次元にしている。

バターワースフィルタ Butterworth filter

2次元周波数分布のsin成分 x[ ] 、cos成分 y[ ] は４隅が原点 ( 0,0 ) （ = 直流成分）になっている。
４隅近傍が低周波成分の領域。
したがって、２次元の Butterworth filter B2[ ][ ]を、x[ ]、y[ ] の４隅を残すように作成している。

ノイズ成分が相対的に多い高周波成分を抑制。

3334

Original Image

Frequency spectrum
Frequency
spacefiltering

FFT

Filtered Image

FilteredFrequency
spectrum

IFFT

35366