平滑化処理プログラム smoothing2.c において、main関数の中で smoothing ( img , img2 ); と書き、関数側は、smoothing ( int x[ ][ 260 ], int y[ ][ 260 ] )と書
書いて、メモリ上では
img、img2 と全く同じ配列
全くを
x y という配列名として宣言している。
x、yという配列名として宣言している
このとき、smoothing関数側は、配列x、yのとき
関数側は 配列の 行の長さはの１行の長さは260 だと知ることが出来るが、列が何列あるか知ることはだと知ることが出来るが 列が何列あるか知ることはできない。
しかし、この関数は x、ｙ の列数を知る必要がない。
その理由を記述して下さい。

正解と考えられる解答例 ：ポインタの示すアドレスから１ピクセル４バイトずつ
２６０ピクセルのデ タを順番に切り出していくだけなので
２６０ピクセルのデータを順番に切り出していくだけなので関数は列数を知らなくてもよい。
与えられたデータをだだ忠実に並べていくだけなので与えられたデタをだだ忠実に並 ていくだけなので関数自体は、いつ終わるのかは判らなくても良い。
配列 の１行の長さを決定することができれば、それを基にして次の列を作ることができるので、が列 何列ある 知 なく も良 。
列が何列あるか知らなくても良い。
１行の長ささえわかれば、その長さでどんどん折り返す
ことができ 値がある分だけ繰り返し続ければ良いのでことができ、値がある分だけ繰り返し続ければ良いので、何列あるかというのは重要ではない。

メイン関数側で、配列 img に入れた画像データは260行 260列 １画素４バイトなので260行、260列、１画素４バイトなので、img の内容を記録したメモリ上の先頭アドレス（ポインタ）から 4 x 260 x 260 バイト目までのデータを 配列 x に託している（エイリアス）。
に託している（エイリアス）
ところが、関数ところが
関数 smoothing は、どこまでがは どこまでが img のデータなのかは、全く知らない。
配列 ｘ、ｙｘ ｙ は、関数
は 関数 smoothig の中では、読み出せとの中では 読み出せと命令された列の配列データは何列でも読む。
例えば 261列目の配列データを使えとプログラム文に書けば 素直に先頭アドレスから順番に数えて 261列目
書けば、素直に先頭アドレスから順番に数えてに相当するメモリ内の数字を読み出す。

void smoothing ( int x[ ][ 260 ], int y[ ][ 260 ] ){for ( j=2; j <= 261; j++ ) { for ( i=2; i<=255; i++ ) {y[ i ][ j ] = x[ i ][ ｊ ] ;}}}例えば、このようなプログラムを書いてしまっても文法ミスは無いので C言語コンパイラはこの文を
文法ミスは無いので、C言語コンパイラはこの文をエラーとは指摘しないで実行する。
滅茶苦茶 意味不明な計算結果が出るだけである
滅茶苦茶、意味不明な計算結果が出るだけである。
つまり、配列の列数の管理は、プログラマに全責任を託されている。
プログラマに全責任を託されている
意味不明な列数を配列に入れないように注意する。

しかし、 C言語の配列の宣言法に慣れていないと、
このプログラムエラー
プ グ（バグ）を犯しやすい。
バグ
iimg [256][256] と 配列を宣言すると、配列を宣言すると
img[0][0] から img[255][255] までの要素ができる。
配列の添字に 256 を記述すると、エラーが起きたり、不明なデータが計算に使われてしまう。
解決策は、画像などのデータがぴったり収まるような配列宣言を避け、すこし余裕のある配列を宣言しよう。
画素デ タの最初の要素を img[0][0] から使えば画素データの最初の要素を問題が生じない。しかし、これは慣れないと扱い難い。

解答 ５

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

//filter.c//----------------------- Display original image ---maxcount = get_int_max(img); disp_img( img, maxcount, 0, 0);//----------------------- filtering --------------printf("¥n¥n 3x3 filter "); scanf("%c",&yn);
filter1( img, img2); disp_img( img2, maxcount, 260*0, 260*1 );filter2( img, img2); disp_img(
disp img( img2, maxcount, 260
260*1,
1, 260*1
260 1 );
filter3( img, img2); disp_img( img2, maxcount, 260*2, 260*1 );filter4( img, img2); disp_img( img2, maxcount, 260*0, 260*2 );filter5( img,
img img2); disp_img(
disp img( img2,
img2 maxcount,maxcount 260*1,
260*1 260*2 );画像 x の 画素最大値 を出力する関数
return関数 ： 関数の出力値を設定する
get_int_max関数は、整数の出力値をもつので、inti関数は 整数の出力値をも ので i で定義。
で定義
int get_int_max( int x[ ][260] ){inti, j, max = 0 ;f (j 1 j< 256 j ){ for(i=1;i<=256;i++){for(j=1;j<=256;j++){f (i 1 i< 256 i ){if(x[i][j] > max) max = x[i][j];}}printf("¥n¥n max count = %d ", max);
return(( max );)}画像 x を 最大値 max、座標（xp、yp）に描画する関数
void disp_img( int x[ ][260], int max, int xp, int yp ){int i, j, count ;for(j=1;j<=256;j++){ for(i=1;i<=256;i++){count = x[i][j]; count *= 255 / max;if(count<0)() count=0;; if(count>255)() count=255;;
SetColor( RGB(count,
RGB(count count,
count count) );SetPixel( i+xp, j+yp );}}}filter.c 実行結果

SPECT（ Single Photon Emission CT ）PET（ Positron Emission CT ） の原理
断層画像を得る方法
フィルタ重畳逆投影法
FBP （ Filtered Back Projectionj）
逐次近似再構成法 IterativeIt
ti R
Reconstructiontti
MLEM （ Maximun
M i
Likelihood
Lik lih d ExpectationEt ti Maximization
M i i ti ）OSEM （ Ordered Subsets Expectation Maximization ）Convolution 重畳積分

(*)

FBP
Filtered
B k
Back
Projection

Si l B kP j ti によるSPECT画像の作成
SimpleBackProjectionによるSPECT画像の作成

サイノグラムの各スライスの１次元配列は、
グ
各
３２方向の各々の角度から収集されたデ タ。
３２方向の各々の角度から収集されたデータ。
このデ タから、収集された各々の角度に傾いた
このデータから収集された各々の角度に傾いた
２次元透視画像（ Pθ ）を作成する。

単純重ね合わせ再構成法 Simple Back Projection
収集された各々の角度に傾いた２次元透視画像（
収集された各角度 傾 た 次 透視 像（ Pθθ ）を
全部単純に重ねると再構成画像ができる。
(回転中心近傍の値が盛り上がった不正確な画像。)
スライスｊにおけるサイノグラムを求める
スライスｊにおけるサイノグラムを求める。
サイノグラムの各スライスの１次元配列は、３２方向の各々の角度から収集されたデータ。
方向 各々 角度から収集されたデ タ
サイノグラムの各スライスの１次元配列から、収集された各々の角度に傾いた２次元透視画像Pθを作成する。
Pθを単純に重ね合わせた画像をｌとすると
ｌ＝∫ Pθ dθ(Simple back projection)ｌ は、回転中心部ほど重ね合せ回数が多くなり、は 回転中心部ほど重ね合せ回数が多くなり中心から距離が遠いほどカウントの低い像になる。

123I-IMP Brain SPECT 単純重ね合わせ像は回転中心近傍のカウントが持ち上がる。画像が不鮮明。

つまり、回転中心からの距離ｒに反比例した濃度に補正するフィルタ １／ｒ
１／を正確な断層像 に畳み込んだ像が l である。
を正確な断層像ｇに畳み込んだ像がである
式で表現すると ｌ＝ｇ＊(１／ｒ) となる。
ｌ、ｇ、１／ｒ のフーリエ変換を Ｌ、Ｇ、F(１／ｒ) と表現すると、畳み込みの定理より
Ｌ＝Ｇ・F(１／ｒ)
Ｌ
Ｇ F(１／ｒ) となる。
２次元フ リエ変換の公式の極座標表現を用いると、２次元フーリエ変換の公式の極座標表現を用いると（ｆｒは周波数空間上の原点からの距離 ）F(１／ｒ) ＝∫∫（１／ ）∫∫（１／ｒ）exp（－ｊ(2πｒｆr（ ｊ(2ｆ ））ｒｄｒｄθ
）） ｄ ｄθ ＝１／ｆｒ
１／ｆ
これより
れより

L ＝ Ｇ ／ｆｒ／

なのでＧ＝Ｌ・ｆｒＧ ＝Ｌ・ｆｒの 意味は、 ２次元周波数空間上で、単純重ね合わせ画像をフーリエ変換した２次元データ Ｌ に、フィルタ関数 ｆｒ（ ｆｒは周波数空間上の原点からの距離）をかけると、正しい再構成画像をフーリエ変換したデータ Ｇになる。

Ｇ＝Ｌ・ｆｒに、畳み込みの定理、
定を用いると、以下のような実空間での計算に変換できる。
この式を 逆フーリエ変換すると、逆フーリエ変換すると
ｇ＝ｌ＊ｈ（ ｈ は フィルタｆｒ の 逆フ
逆フーリエ変換
リエ変換 ）

この式に、ｌ＝∫ Pθ dθ を代入すると、
ｇ＝ ∫ Pθ dθ ＊ｈｇ＝ ∫( Pθ ＊ｈ）dθ （ｈはθと独立した値なので交換可）
ｇ＝ ∫ Pθ dθ

(Pθ ＝ Pθ ＊ｈ ）FBPの式Pθ に 実空間フィルタ h ( = frの逆フーリエ変換 ) を畳み込めば、重ね合せると正確な断層像
重ね合ると 確な断層像 ｇ になる２次元透視画像
なる 次 透視 像 Pθ を算出できる。 これを Filtered Back Projectionj（FBP） という。
周波数空間での実際の計算においては、フィルタ Ｈ（ =ｆｒ）はｆｒ）は常に正の値であり（絶対値）、
さらにサンプリング定理より、ナイキスト周波数以上の成分をさらにサンプリング定理よりナイキスト周波数以上の成分を削除する必要があるので、

フィルタ逆重畳画像再構成法
Filtered Back Projection (FBP)
サイノグラムの２次元透視画像 Pθ に 、
実空間フィルタ h ( = frの逆フ
frの逆フーリエ変換
リエ変換 ) を畳み込めば、を畳み込めば重ね合せると正確な断層像 ｇ になる２次元透視画像 Pθ を算出 きる
算出できる。

周波数空間での再構成フィルタ Ｈは、
Ｈ ＝｜ｆｒ｜（ｆｒがナイキスト周波数未満の場合）
Ｈ ＝０（ｆｒがナイキスト周波数以上の場合）となる。 これをＲａｍｐフィルタという。
Ｒａｍｐフィルタを逆フーリエ変換して実空間Ｒａｍｐフィルタｈにしてから、実空間でPθにｈを畳み込む。

Pθ = Pθ * ｈ（ * は畳み込み演算 ）ｈを畳み込んだ２次元透視画像 Pθ を単純に重ね合わせた画像 g は、
g =∫ Pθ ｄθ

（ Filtered back projection ）Filtered back projection は、
回転中心部ほど重ね合わせ回数が多くなる誤差がフィルタｈによって補正され、正しい再構成画像となる。
タ によ
補 され
し 再構成画像となる
ｈを畳み込んだ２次元透視画像 Pθ を単純に重ね合わせた画像 g =∫ Pθ ｄθ
（ Filtered back projection ）は、回転中心部ほど重ね合わせ回数が多くなる誤差がフィルタｈによって補正され、正しい再構成画像となる。

Ｒａｍｐフィルタは理論的には正確な再構成フィルタだが、実際の臨床データに適用すると実際の臨床デタに適用すると ナイキスト周波数以上の高周波成分を不連続に遮断するために生じる高周波ノイズや 放射状
ア チファクトが目立つ場合が多いため、アーチファクトが目立つ場合が多いため、臨床では高周波成分を抑制する工夫を施した再構成フィルタが用いられている。
用いられている
ＳＰＥＣＴで用いられる一般的な再構成フィルタとして、Shepp＆Loganフィルタがある。
Shepp＆Loganフィルタがある
周波数空間上で、ナイキスト周波数近傍の高周波成分を連続的に減衰させるように設計されており再構成画像に生じる高周波ノイズや放射状アーチファクトを抑制する効果をもつ。
抑制する効果をもつ

実空間フィルタを畳込んだ２次元透視画像（Pθ）を重ね合わせると フィルタ逆投影再構成像ができる。
フィルタ逆投影再構成像ができる
フィルタの形状で、再構成画像の高周波成分が変る。

123I-IMP Brain SPECTFBP with Ramp filter

123I-IMP Brain SPECT FBP with Shepp&Logan filter

吸収補正後の 18F-FDG 脳PET サイノグラム

脳PETのサイノグラムデ タを並べ替えると脳PETのサイノグラムデータを並べ替えるとプロジェクション像

脳PETのサイノグラムに Rampフィルタを重畳して重ね合わせ

ＦＢＰ画像

畳込みの定理 Convolution

畳み込みの定理
デ タ g をフーリエ変換して、データをフ リエ変換してその周波数空間成分 Ｇ に
周波数空間 Rampフィルタ
R
フィルタ H をかけて逆フ リ 変換すると逆フーリエ変換すると、実空間で 実空間 Rampフィルタ
実空間で、実空間
ルタ ｈ をｇ に畳み込みしたデータと同じになる。
に畳み込みしたデ タと同じになる( G x H と ｇ * h は等価演算
等価演算 )

Ramp フィルタの作成
Ramp.c 実行結果

//------ Generate frequency space Ramp filter -------f0 = 128; // Nyquist
for ( i=-200; i<=200; i++){if ( abs( i ) <= f0 ) R[ i+200 ] = ( double) abs( i ) ;else
R[ ii+200
200 ] = 0.0 ;}if ( A ) B ; else C ;Aが正しければBを実行
Aが正しければBを実行、
そうでなければCを実行。
abs ( ) ;絶対値を計算する関数

//----------- Transform Ramp filter ------------ffor(( i=1;i 1 i<=128;
i< 128 ii++ ){ x[i][i]= R[ 200
200+ii ] ; }for( i=1; i<=128; i++ ){ x[i+128] = R[ 329-i ] ; }１次元 FFT を実行する１次元データは、両端に直流成分（周波数 0 ）、中央が高周波（ナイキスト）に） 中央が高周波（ナイキスト）になるように並べ替える。

//--------- Save Frequency space Ramp filter --------printf("¥n¥n
i f("¥ ¥ SaveS
FrequencyFspace Ramp
R
filterfil");") scanf("%c",&yn);
f("% " & )GetFileName( f, 1);

fp = fopen( f, "wt");
wt );for( i=1; i<=256; i++){ fprintf ( fp, “ %lf ¥n”, x[ i ] ) ; }fclose( fp );
周波数空間の Ramp filter を ファイルに保存する。
GetFileName(( f, 1);) で、保存するファイル名を指定する
保存する
イ 名を指定する
ダイアログが開く。
ダイアログが開く
fp=fopen( f, “wt”); で、テキスト形式ファイルとして書き込む。
fprintfp( fp,
p, “ %lf ¥n”,, x[[ i ] ) ;file - pprintf

fprintf () 関数は、ファイルに文字や数字を書き込む関数。

グラフ描画関数 Graph ( ) の作成
座標（ 1 y1）
座標（x1,
1） から (x2,( y2）
） の範囲に、配列の範囲に 配列 x[[ ] のグラフをのグラ を最大値 max、 color 色 で描画。
void Graphp ( double x[ ], int x1, int yy1, int x2, int y2,
y double max,
COLORREF color, char *s ){int i,i ip,
ip jp,
jp ip2,
ip2 jp2 ;PenWidth(3);S C l (WHITE) Li
SetColor(WHITE);
Line(x1,( 1 y1,
1 x2,
2 y1);
1) Line(x1,Li ( 1 y1,
1 x1,
1 y2);
2)
SetColor(color); FontSize(20); DrawString(x1+10, y2, s );for(i=1; i<=255; i++){ipp = x1 + i ; ip2
p = x1 + ((i+1)) ;jp = y1 - (int)(x[i]*100./max); jp2 = y1 - (int)(x[i+1]*100./max);Line( ip
ip, jpjp, ip2
ip2, jp2 );}}グラフ描画関数 Graph ( ) の 最後の引数 s は、グラフに添付する文字列の、メモリ上の先頭アドレス（ポインタ）。
したが て この引数の変数型指定は ｃharしたがって
h * となる。
となるDrawString(( x, y, s )グラフィック画面の座標 ( x,x y ) に、FontoSize関数で指定した文字サイズで、文字列 s を書く。

色を指定する変数 color の型は、COLORREF型 になる。

fprintf ( fp, “ %lf ¥n”, x[ i ] ) ;fp は、保存するファイルのポインタ。
¥n” , x[i] は、printf関数と同じ様式で、“ %lf ¥n実数変数 x[i] の内容を、実数形式（%lf）でファイルに書く。
作成されたファイルを メモ帳などで開くとテキスト形式で x[i] が保存されることを確認できる。（¥nで各数字が改行される。）

ここでは、freq_ramp.txt と命名（どんなファイル名でもよい。）

//------ Transform Ramp filter into Real space -------printf("¥n¥nTransform
f("
f
Ramp
R
filterf
into RealR
space");scanf("%c",&yn);")
f("% " & )for( ii=1;
1; i<=256;i< 256; i++ ){ y[ i ] = 0.0
00 ; }
FFT1D( x , y , 256 , -1 ) ;１次元
次元 高速 逆逆フーリエ変換
リ 変換 IFFT
１次元の周波数成分（x に実数（sin）成分、に実数（sin）成分 ｙ に虚数（cos）成分 = 0）を逆フ
0）を逆フーリエ変換して、配列リエ変換して、配列 x に実空間実数成分、配列
実
実数、
y に実空間虚数成分を出力。
実数
３個目の引数は、配列のデータ数。
４個目の引数が -１の場合は IFFT を実行。

実空間 Rampフィルタ の保存
printf("¥n¥n
i f("¥ ¥ SaveS
RealR l space Ramp
R
filterfil");") scanf("%c",&yn);
f("% " & )GetFileName( f,f 1 );

fp = fopen( f,f "wt");wt );for(i=1; i<=128; i++){ fprintf( fp, "%.16lf ¥n", x[ i ] ) ; }fclose(fp);//-------- Integration of Real space Ramp filter ------for(( sum=0.0 , i=1; i<=256; i++ ){ sum += RR[[ i ];] }printf("¥n¥n
printf(¥n¥n Integration of Real space Ramp ¥n = %%.4lf
4lf ",sumsum ) ;

実空間フィルタの積分値は １ になる必要がある。
になる必要がある（処理後にカウント総和が変化しないことが必要。）

printf("¥n¥n Transform Ramp filter into frequency space ");scanf( %c ,&yn);
scanf("%c",&yn);

for(i=1; i<=256; i++){ y[ i ] = 0.0 ; }FFT1D( x, y, 256, 1) ;１次元高速フーリエ変換 FFT
配列 x に実空間Rampフィルタ（ y には ０） を入力し、FFT を行うと、元のRampフィルタになることを確認。
を行うと 元のRampフィルタになることを確認convolution.c の実行結果
半円形のサンプルデ タ g に
半円形のサンプルデータ
実空間で Rampフィルタ
R
フ ルタ を畳み込む。
畳み込む//------Generate sample data g[ ] : Half Circle -------f ( i=fori 1 ; i<=
i< 77;
77 i++){i ){ g[i][i] = 0.0;00 }for ( ii=78
78 ; i<=178;
i< 178; i++){}g[i] = sqrt(abs(50.0*50.0-((double)i-128)*((double)i-128)));for ( ii= 179; i<=i< 256; i++){i ){ g[i] = 0.0; }サンプル実空間データ
サンプル実空間デタg
として、
半径 50 ピクセルの半円
g[i] を作成。

//---------- Load Real space Ramp filter h[ ] -------printf(
i tf("¥¥n¥n¥ LoadL dR
Reall space R
Ramp filt
filter ")); scanf(f("%%c",&yn);& )GetFileName( f, 0 ) ; fp = fopen( f, "rt") ;for(i=0; i<=127; i++){ fscanf( fp, "%lf ¥n", h + i ); }fclose(fp) ;実数１次元配列として宣言した h[ ] に、Ramp.c で作成した
実空間 Rampフィルタを読み込む。
Rampフィルタを読み込む
ffscanff の ３番目の引数は、読み込んだファイルデータを３番目の引数は 読み込んだフ イルデ タを書き込むアドレス（ポインタ）を記述するので、h+ii と書くと、書き込むアドレス（ポインタ）を記述するので、h配列 h[0] から h[127] に データが書き込まれる。
（ 単純に fscanf( fp, “%lf ¥n”, &h[ i ] ) ; と書いても良い。 ）畳み込み Convolution

l=g*h

サンプルデータ g[ ] に 実空間Rampフィルタ
h[ ] を畳み込んで、配列を畳 込ん
列 l [ ] に書き込む。
書き込む//--------- Convolution h[ ] into g[ ] --------------printf("¥n¥n Convolution Real space Ramp filter "); scanf("%c",&yn);

for( i=1; i<=256 ; i++ ) {gh = 0.0 ;for(j=0; j<=127; j++){ if( i+j <= 256) gh += g[ i+j ] * h[ j ]; }for(j=1; j<=127; j++){ if( i-j >= 1 ) gh += g[ i-j ] * h[ j ]; }}l[i] = ghh;実空間 Rampフィルタ
ph[[ ] の積分値は １ なので、畳み込む前後の
ｇ[ ] と l [ ] の積分値は、ほぼ同じ。

周波数空間で Rampフィルタ を かける

convolution_FFT.c
convolution
FFT c
の実行結果

実空間 サンプルデータ ｇ をフーリエ変換して、周波数空間 実数成分 Gr、
周波数空間の虚数成分 Gi を表示。
を表示
実空間
実
間 サンプルデータ
タｇの積分値 が、周波数空間の実数成分 Gr の最初の成分（原点、直流成分）と同じ値になることを確認する。

周波数空間の 実数成分 Gr、
虚数成分 Gi の二乗和の平方根（パワースペクトル
方根（パ
ペク
G）
を表示。
を表示
実空間
実
間 サンプルデータ
タｇの積分値 が、周波数空間のパワースペクトル G の最初の成分（原点、直流成分）と同じ値になることを確認する。

//------ Convolution H[ ] into G[ ] ------for (i= 1; i<= 256; i++ ){ GHr[i] = Gr[i] * H[i] ; }for (i= 1; i<= 256; i++ ){ GHi[i] = Gi[i] * H[i] ; }//- Transform convolved GH into real space FFT1D( GHr, GHi, 256, -1 ) ;データ g の周波数成分 Gr、Gi に、
周波数空間 Rampフィルタ H をかけて、逆フーリエ変換する。

データ g の周波数成分 Gr、Gi に、
周波数空間 Rampフィルタ H をかけて、逆フーリエ変換すると、逆変換す実空間で 実空間 Rampフィルタ
実空間で、実空間
R
フィルタ ｈを ｇ に畳み込みしたデータと同じになる。 ( G x H と ｇ * h は等価演算 )

積分値が 元データ ｇ と同じことを確認する。
（カウント総和は保たれる。）

課題５ プログラム Ramp.c を改良して
実空間 Shepp&Loganフィルタを計算するプログラム
SheppLogan.cgを作って下さい。
周波数空間Shepp&Logan の式

絶対値を計算する関数は abs( )y=|x| →

y = abs(x);階乗を計算する関数は pow( , )y = x2 → y = pow(x,2.0);周波数空間Shepp&Logan の式の f0 は高周波遮断周波数。
f0 をプログラムで設定できるように
printf("¥n¥n Cut off freq. (0～0.5) = "); scanf("%lf",&CutOff );を入れてください。
ここで CutOff の単位は cycles/pixel 。 Shepp&Loganの式に使われる f0 の単位は cycles
l / 画像の１辺長 = cycles
l / 256 pixeli lなので f0 = 256.0 * CutOff ; と変換してから式に代入してください。

作成したプログラムとフィルタのグラフ画像をメールに添付して提出してください。
メールに添付して提出してください