核医学検査技術学実習
実習３ フィルタ処理、畳込み画像処理において用いられる実空間、空間周波数フィルタ処理とその演算方法である畳込みについて学ぶ。
Practice3.zipをダウンロードして解凍（解凍ソフトLhacaをインストールして、Lhacaのアイコン上へ、解凍したいzipファイルをドロップダウンする。）

解凍ソフトLhacaがうまく作動しない
場合は、７zipをインストールして、zipファイルを解凍して下さい。

FFTFast Fourier Transform

FFT.exe を起動し、open fileから同フォルダ内にあるboneを選択する。
※bone：99mTc-MDPのSPECT画像。骨シンチ。
FFTをクリックして2次元フーリエ変換を行い、replaceをクリックして、パワースペクトル画像を再配置する。
discoverをクリックすると、逆フーリエ変換された画像が表示される。
coverをクリックすると、右側の画像が消える。

右中の枠をクリック・ドラッグすると、パワースペクトルが部分的に逆フーリエ変換される。

coverをクリックすると初期化される。
色々試してみてください。

実空間フィルタ
フォルダFilter内にあるFilter.exeを起動し、 「LoadFile」ボタンを押してBitmapまたはJPEG画像を選択する。（プログラムフォルダ内にあるsample.jpgを使用しても良い）「Kernel」の数値（９個のフィルタ行列の要素）を変更し、任意の3×3の実空間フィルタを作成する。
「Filter」ボタンを押し、作成したフィルタによりフィルタ処理を行う。

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

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

平滑化フィルタ

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

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

解答 ５

void filter1( int x[ ][260], int y[ ][260] ){// Weighted smoothing filter
int
i, j ;
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] ;}}for(j=2;j<=255;j++){ for(i=2;i<=255;i++){ y[i][j] /= 16; }}}輪郭の鮮明さを維持するために中心部に重みを付けたスムージングフィルタ。

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++){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] /= 9; }}}最も単純なスムージングフィルタ。
中央画素の重み付けがない。
輪郭の鮮明さが損なわれる。

void filter3( int x[ ][260], int y[ ][260] ){// Vertical Prewitt filter
int i, j ;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] ;}}}Prewitt フィルタ （垂直方向）
１次微分（差分）エッジ検出フィルタ
上下方向に画素値が大きく変化する部位（輪郭）の抽出フィルタ

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

void filter5( int x[ ][260], int y[ ][260] ){// Unsharp masking filter
int i, j ;for(j=2;j<=255;j++){ for(i=2;i<=255;i++){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]*４+ x[i+1][j-1]*-1+ x[i-1][j+1]*0 + x[i][j+1]*-1 + x[i+1][j+1]*0 ;}}}４近傍 先鋭化フィルタ
アンシャープマスキング。
ぼやけた部位をマスクする。
上下および左右方向に画素値が大きく変化する部位（輪郭）の強調

空間周波数フィルタ

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

「display sin cos component」ボタンを押し、原画像のフーリエ変換を行い、のsin成分とcos成分の画像を表示させる。
「FFT spectrum」ボタンを押し、パワースペクトルを表示する。

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

median.exe 実行結果
テキストボックス内にカーソルをおいて実行してください。

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++ ;}}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;y[ i ][ j ] = med[ 5 ];

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] に代入している。
配列 med は、 int med[10]; と宣言している。
int med[9]; でも良さそうに思えるが、C言語では、配列数は ０番目から数えている のでmed[9]; と宣言すると、med[0] から med[8] の９個の配列要素だけが宣言され、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;med[1] から med[9] の値を並べ替えている （ソート sort）。
この for ループが終了すると、 med[1] から med[9] の中の最小値が med[1] に、最大値が med[9] に代入されている。
バブルソートのアルゴリズム

Bubble sort algorithm

コップ底の泡が、軽いものから先に順番に昇っていくような様子をプログラムで表現している。ソート法のアルゴリズムとしては最も非効率的な方法だが単純で理解しやすい。

はじめに ni が 1 の 値で nj の ループに入る。
nj の 値は 8 から１ までの値に変化しながら ループを回る。
最初は nj が 8 なので、if 文の中身を具体的に書くとif ( med[ 9 ] < med[ 8 ] ){min=med[ 9 ]; med[ 9 ]=med[ 8 ]; med[ 8 ]=min;}もし med[9] が med[8] より小さければ med[9] と med[8]の値を入れ替える。 変数の入れ替え作業を swap という。
swap作業には 第三の変数（ここでは min）が必要になる。
ni が 1 の 状態で、nj が 8 から １ まで変化し終わると、med[1]の値は、med[1]からmed[9]の最小値になっている。
小さい値が泡のように上に昇っていくので バブルソート法。

次に ni が ２ の 値で nj の ループに入る。
nj の 値は 8 から ２ までの値に変化しながら ループを回る。
（ med[1] は最小値が確定なので ループに入る必要がない。）

ni が 2 の 状態で、nj が 8 から 2 まで変化し終わると、med[2] の値は、med[2] から med[９] の中の最小値になっている。
次に ni が ３ の 値で nj の ループに入る。
nj の 値は 8 から ３ までの値に変化しながら ループを回る。
（ med[1]、med[2] は １番目、２番目に小さい値が確定なのでループに入る必要がない。）

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

smoothing filter は、輪郭が不明瞭化し画像がぼやける。
median filter は smoothing filter より輪郭が保たれる。
median filter は、微細なノイズを除去する効果がある。

bone6noise は、10ピクセルおきに画素値 100 を入れた像。
median では、このノイズを完全に消すことができる。
smoothing ではノイズを完全には消せない。

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

FFT2

Fast Fourier Transform

FFT2.exeを起動し、open fileから同フォルダ内にあるboneを選択する。
※bone：99mTc-MDPのSPECT画像。骨シンチ。
FFTをクリックして2次元フーリエ変換を行い、replaceをクリックして、パワースペクトル画像を再配置する。
discoverをクリックすると、逆フーリエ変換された画像が表示される。
coverをクリックすると、右側の画像が消える。

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

FFT.exe２次元 FFT の出力である ２次元周波数分布
データは、４隅が原点 ( 0,0 )（= 直流成分）になっているのが面倒な点である。

中心が原点になるように各々 1/4 領域を反転する。

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

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

出力画像

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

バターワースフィルタ

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

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

Butterworth filter

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

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

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

Original Image

Frequency spectrum
Frequency
spacefiltering

FFT

Filtered Image

IFFT

FilteredFrequency
spectrum