2023/4/14

核医学検査技術学実習
実習３ フィルタ処理、畳込み画像処理において用いられる実空間、空間周波数フィルタ処理とその演算方法である畳込みについて学ぶ。

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

Practice3.zipをダウンロードして解凍（解凍ソフトLhacaをインストールして、Lhacaのアイコン上へ、解凍したいzipファイルをドロップダウンする。）
1

2

FFT

Fast Fourier Transform

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

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

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

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

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

2023/4/14

平滑化フィルタ

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

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

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

78void 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; }}}輪郭の鮮明さを維持するために中心部に重みを付けたスムージングフィルタ。

解答 ５
910

void filter2( int x[ ][260], int y[ ][260] ){// Simple smoothing filter
inti, j ;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++){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; }}}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 フィルタ （垂直方向）
最も単純なスムージングフィルタ。
中央画素の重み付けがない。
輪郭の鮮明さが損なわれる。

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

11122

2023/4/14

void filter4( int x[ ][260], int y[ ][260] ){// Horizontal Sobel filter
int i, j ;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++){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] ;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 ;}}

}}
}

}

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

４近傍 先鋭化フィルタ
アンシャープマスキング。
ぼやけた部位をマスクする。
上下および左右方向に画素値が大きく変化する部位（輪郭）の強調

1314

空間周波数フィルタ

テキストボックス内に0.5～0.01の任意の遮断周波数を設定する。

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

「Butterworth filter」ボタンを押し、設定した遮断周波数からバターワースフィルタを作成し、空間周波数フィルタ処理を行う。

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

「IFFT」ボタンを押し、逆フーリエ変換を行い、フィルタ処理後の画像を表示する。

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

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] ){}}

}

17}}

}

min=med[nj+1]; med[nj+1]=med[nj]; med[nj]=min;y[ i ][ j ] = med[ 5 ];183

2023/4/14

k = 1；
for( mi=-1; mi<=1; mi++){ for( mj=-1; mj<=1; mj++) {}}for ( ni=1; ni<=8; ni++) { for ( nj=8; nj >= ni; nj - -) {if ( med[nj+1] < med[nj] ){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] は宣言されない。
19}

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

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

20

はじめに 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 が ２ の 値で 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] ; となる。

2122

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

高速フーリエ変換 （ FFT : Fast Fourier Transform ）ni が 1 の 状態で、nj が 8 から １ まで変化し終わると、med[1]の値は、med[1]からmed[9]の最小値になっている。

median filter は smoothing filter より輪郭が保たれる。
median filter は、微細なノイズを除去する効果がある。
bone6noise は、10ピクセルおきに画素値 100 を入れた像。

フーリエ変換を高速に計算するアルゴリズム。
1942年に Danielson と Lanczos が発明。

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

プログラムによるフーリエ変換は、波形データが離散的smoothing ではノイズを完全には消せない。

なので、基本式に現れる e（２πｆ t ）i の項を、Ｗ ｎｋ と変形して巧みに解いている。
高速逆フーリエ変換（ IFFT : Inverse FFT ）も ほとんど同じアルゴリズムなので、両方可能な関数が作れる。
FFT 、 IFFT におけるデータの制限としては、データ数が ２の階乗でなければならない。

23244

2023/4/14

FFT2

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

Fast Fourier Transform

FFT.exe

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

２次元 FFT の出力である ２次元周波数分布

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

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

2728

Butterworth.exe

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

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

バターワース
フィルタ

出力画像

29305

2023/4/14

画像 を 周波数空間でフィルタ処理する。
テキストボックス内に0.5～0.01の任意の遮断周波数を設定する。

バターワースフィルタ Butterworth filter
ノイズ成分が相対的に多い高周波成分を抑制。

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

31

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

32

Original Image

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

Frequency spectrum
Frequency
spacefiltering

FFT

Filtered Image

FilteredFrequency
spectrum

IFFT

3334356