７月１日 期末試験 核医学機器工学概論
試験問題は、以下の 問題１ から問題１０ のうち４問、（画像サイズ等は変更して出題 = 丸暗記では答えられない）および応用問題１問を出題します。 （合計５０点）
プログラム文終端の セミコロン ； の記述忘れに注意。
変数の型に注意。 キャストの記述忘れに注意。
配列の添字に注意。宣言範囲外の変数を入れないように。
括弧 ( ) [ ] { } の使い分けに注意。混同しないように。
問題 １.
整数１からＮまでの 合計 sum、平均 mean、 分散Variation v（v＝｛ Σ{ (mean－i) * (mean－i) } ｝/ N ） および標準偏差 Deviation sd （ sd = sqrt(v); sqrt( )関数は平方根を計算するＣ言語関数）を求めるＣコードを記述して下さい。

void Main(void){intdouble// 問題 １;;

TextWindow(0,0,300,300); Title(“Calculate Mean、SD");
printf("¥n N = ") ; scanf(sum = 0 ;);;printf("¥n Sum= %d" , sum ) ;printf("¥n Mean= %lf“ , mean ) ;printf("¥n Variation = %lf“ , v);
printf("¥n Deviation = %lf“ , sd);}問題 ２.
１画素２バイトの 256x256画素の 画像ファイルを読込んで配列 img[260][260] に入力するプログラムを記述して下さい。
void Main(void){char a, b, f[100] ;int
i, j, img[260][260] ;GetFileName( f, 0 ) ; printf( “ file = %s ” , f ) ;}問題 ３.
配列 img[ ][ ] に入力された 256x256画素の画像（１画素の実長は 1.695mm、 画像サイズは 43.4 x 43.4cm ）の計数密度（count/cm2）を求めるプログラムを記述して下さい。
ただし、カウントが 0 の部位は計数に加えないようにする。
//---------- Calculate count density ---------int
doublei, j, total_count, pixel , img[260][260];
area, density ;printf("¥n¥n Count density = %.2lf (count/cm2)", density );問題 ４.
配列 x[260][260] に入力された 256x256画素の 画像に３ｘ３メディアンフィルタをかけて配列 y[260][260]に出力する
関数 median（） を記述して下さい。
void median( int x[ ][260], int y[ ][260] ){int i, j, k, mi, mj, ni, nj, min, med[10] ;}問題 ５.
配列 x[260][260] に入力された 256x256画素の画像に、右図に示す３ｘ３平滑化フィルタ（移動平均フィルタ）をかけて配列 y[ 260][260]に出力する関数 smoothing（） を記述して下さい。
void smoothing( int x[ ][260], int y[ ][260] ){}問題 ６.
配列 x[260][260] に入力された 256x256画素の画像の最大画素値を整数値で出力する関数 get_int_max（） を記述して下さい。
get_int_max( int x[ ][260] ){}問題 ７.
１次元データを g(x)、フィルタ関数を h(x)、フィルタ処理後のデータを l(x)、それぞれのフーリエ変換を G(f)、H(f)、L(f) として、（ x は実空間座標、 f は周波数 ）畳み込みの定理 を 文章と数式を用いて 説明して下さい。

問題 ８.
データ数６４の１次元データ ｇ[0]～ｇ[63] と、データ数３２の実空間フィルタ h[0]～h[31]がある （ h[0] が原点）。
データ ｇ に フィルタ h を畳み込んで配列 l [0]～ l [63] に書き込むプログラムを記述して下さい。
//--------- Convolution ｈ[ ] into ｇ[ ] ----------int
doublei, j;g[64] , h[32] , l[64] , gh ;

データ数６４の１次元データ ｇ[0]～ｇ[63] と、右図に示すようなデータ数３２の実空間フィルタ
h[0]～h[31] がある （ h[0] が原点）。
データ ｇ に フィルタ h を畳み込んで配列 l [0]～ l [63] に書き込むプログラムを記述して下さい。
//--------- Convolution ｈ[ ] into ｇ[ ] ----------int
doublei, j;g[64] , h[32] , l[64] , gh ;for( i = 0 ; i <=63 ; i++ ){// i <64 と書いても同じ

gh = 0.0 ;for( j = 0; j <=31; j++ ){ if( i+j <=63 ) gh += g[ i+j ] * h[ j ] ; }for( j = 1; j <=31; j++ ){ if( i-j >= 0 ) gh += g[ i-j ] * h[ j ] ; }}l[i] = gh ;ｇ に h を畳みこんだ結果を l とする。座標 i における ｌ の値 l[i] は、
ｌ[ i ] = g[ i ]*h[0] + g[i+1]*h[-1] + g[i+2]*h[-2] + g[i+3]*h[-3] + ・ ・ ・+ g[i-1]*h[ 1] + g[i-2]*h[ 2] + g[i-3]*h[ 3] + ・ ・ ・
ここで フィルタ h は偶関数（左右対称）なので、h[-j] = h[j] を代入し、
ｌ[ i ] = g[ i ]*h[0] + g[i+1]*h[1] + g[i+2]*h[2] + g[i+3]*h[3] + ・ ・ ・+ g[i-1]*h[1] + g[i-2]*h[2] + g[i-3]*h[3] + ・ ・ ・
これを Ｃ言語で表す。 ｈ の要素数 j は h[0] から h[31] まで。
g の 要素数 は g[0] から g[63] までなので、g の 添字を表す i+j は63 以下、 i-j は 0 以上の場合だけ加算する条件式を入れる。
gh = 0.0 ;for( j = 0; j <=31; j++ ){ if( i+j <=63 ) gh += g[ i+j ] * h[ j ] ; }for( j = 1; j <=31; j++ ){ if( i-j >= 0 ) gh += g[ i-j ] * h[ j ] ; }l [ i ] = gh ;これを i が 0 から 63 まで繰り返すと、 l [0] から l [63] が計算される。

問題 ９.
64x64画素の 画像 x[ ][ ] を 角度 T (度） 回転させて配列 y[ ][ ] に格納するプログラムを記述して下さい。
ただし座標が整数値であることに伴う、回転後に画素値が割り当てられない座標が生じる誤差は無視してよい。
//----------------- Image Rotation ------------------int
i , j , ir , jr ;double I , J , T , rot , x[64][64] , y[64][64] , PAI = 3.1415926 ;問題 1０.
32方向から撮像されたプロジェクション Proj[ i ][ j ][T]( i、j は画像の横、縦座標、Tは撮像方向 ) から、フィルタ逆投影再構成法で SPECT像を作成する手順を、以下のキーワードを用いて（不要なキーワードもある）、文章、式および図などで説明してください。
実空間フィルタ h、周波数空間フィルタ H、
１次元フーリエ変換、２次元フーリエ変換、１次元逆フーリエ変換、２次元逆フーリエ変換、畳み込み、実空間、周波数空間、実数成分、虚数成分、２次元透視画像 Pθ、 サイノグラム、Rampフィルタ、 Shepp＆Loganフィルタ、
ナイキスト周波数Convolution 重畳積分(*)

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

//------------- Disp Sinogram -------------------------for(frame=0; frame<MAXFRAME; frame++){ for(i=0;i<64;i++){Sino[ i ][ frame ] = Proj[ i ][ RJ ][ frame ] ;}}for(j=0;j<64;j++){ for(i=0;i<64;i++){c = Sino[i][j] * 800/max; if(c>255)c=255;
for(ii=0;ii<=3;ii++){ for(jj=0;jj<=3;jj++){
SetColor(RGB(c,c,c)); SetPixel(4*i+ii+260,4*j+jj+140);}}}}
Flush( ); // 画像描画の後に Flush( ); を記述すると描画がスムーズになる。
SetColor(WHITE); FontSize(30);DrawString(300, 300, "Sinogram" );サイノグラムの各スライスの１次元配列は、３２方向の各々の角度から収集されたデータ。
このデータから、収集された各々の角度に傾いた２次元透視画像（ Pθ ）を作成する。

//------------ Generate Pth[ ][ ][ ] (Pθ) data --------------PAI = 3.1415926/180.0;

dT = 180.0/(double)MAXFRAME;

for(T=0;T<MAXFRAME;T++){ for(j=0;j<64;j++){ for(i=0;i<64;i++){Pth[i][j][T] = -1;}}}for(T=0; T<MAXFRAME;T++){for(j=0;j<64;j++){for(i=0;i<64;i++){ Simple_Sino[i][j] = Sino[i][T]; }}rot = (double)T * dT * PAI ;// 画像の回転
for(j=0;j<64;j++){for(i=0;i<64;i++){ I = (double)i ; J = (double)j ;ir = (int)( ( I-31.5 )*cos(rot) - ( J-31.5 )*sin(rot) + 31.5 );jr = (int)( ( I-31.5 )*sin(rot) + ( J-31.5 )*cos(rot) + 31.5 );if(ir>=0 && ir<64 && jr>=0 && jr<64)Pth[ir][jr][T]=Simple_Sino[i][j];}}}画像の回転は、 画像の中心を原点に平行移動してから回転移動し、再び元の座標に平行移動。

回転後の画像には、画素値が入力されない画素が少数だが発生する。（画素の座標は整数のため）
そのため、配列 Pth[ ][ ] には初めに -１を代入して回転後にも -１であれば画素値が入力されなかったと判断し、その画素の前後の平均値を代入している。
for(T=0; T<MAXFRAME;T++){for(j=1;j<63;j++){for(i=1;i<63;i++){if(Pth[i][j][T]==-1 && Pth[i-1][j][T]!=-1 && Pth[i+1][j][T]!=-1)
Pth[i][j][T] = (Pth[i-1][j][T]+Pth[i+1][j][T])/2;
if(Pth[i][j][T]==-1 && Pth[i][j-1][T]!=-1 && Pth[i][j+1][T]!=-1)
Pth[i][j][T] = (Pth[i][j-1][T]+Pth[i][j+1][T])/2;}}
for(j=0;j<64;j++){for(i=0;i<64;i++){ if(Pth[i][j][T]<0)Pth[i][j][T]=0; }}}プログラム PETFBP.c
PET画像を Simple BackProjection、または FBP( Filtered BackProjection ) で再構成する。

PETFBP.c を実行して、Simple BackProjection とFBP の 再構成画像の違いを観察し、Ramp等のフィルタ関数の必要性を理解してください。

PETsinpgram ファイル
PETの収集データは各スライスのサイノグラムが並んでいる ３次元データ。（ CTと同じ ）。

各スライスのサイノグラムが並んでいる３次元データを並べ替えれば、SPECTのプロジェクション像と同じ並びになる。 Sinogram[ i ][θ][ j ] = Projection [ i ][ j ][θ]//------------- GET sinogram -------------------START: printf("¥n¥n Load PET Singram data ") ;GetFileName( f, 0) ; fp = fopen( f, "rb“ ) ;OFFSET = 24; fseek( fp, OFFSET, SEEK_SET ) ;for(slice=0;slice<MAXSLICE;slice++){ Event( ) ;for(j=0; j<128; j++){ for(i=0; i<256; i++){fread( &data, sizeof(float), 1, fp ) ;y[ i ][ j ][ slice ] = data ;}}} fclose(fp);データファイルのオフセット（ OFFSET )DICOMなどの医用画像データファイルは、画像の画素数や患者名などの数字や文字情報と画素値情報がつながって入っている。
用意したPETsinogramファイルは、先頭２４バイトには数字や文字情報が入っている。
画像を取り出したい場合は、２４バイト読み飛ばす必要がある。これを OFFSET という。

OFFSET = 24; fseek( fp, OFFSET, SEEK_SET ) ;OFFSET を 入れる変数は巨大整数を扱える必要があり、倍精度整数で宣言する（ long 型）。
fseek 関数は、何バイト目からファイルの読み書きを開始するかを指定する。
fseek( 読むファイルのポインタ、OFFSET、SEEK_SET ）；fread( &data, sizeof(float), 1, fp ) ;PETsinogramファイルの画素は、単精度実数（float型）で書き込まれている。fread関数は、各種の型で読み出せる。
fread（ 読んだデータを保存する場所（変数のアドレス）、読み込むデータの型 を sizeof(型) と記述、読み込むデータの数 （ここでは１個づつ読む）、読むファイルのポインタ）

fseek や fread関数 などの標準C言語関数の文法を調べるオプションが、Visual C++ に装備されている。
たとえば freadの文法を調べたい場合は、fread の５文字のどこかにマウスカーソルをおいて左クリックしてからファンクションキー１（F1） を押すと、Microsoft の 文法解説
サイトの fread の説明ページが開く。

Event 関数、Flush 関数 ( 大変重要な 裏関数 )ここでは３次元画像データを読むため for文が３重にループしている。長時間このプログラム実行にWindows OS が占領される可能性がある。 その間はマウスやキーボードなどを操作しても Windows に無視される（ソフトが固まった状態）。
この不都合を避けるため、多重ループの中に、時々マウスやキーボードなどの操作状態（ハンドル）を Windows に渡す
関数 Event( ) を入れる。 引数は無いので括弧内は空白。
類似した関数として Flush( ) がある。 描画に時間のかかる画像を描くときに、描画が固まったような状態になる場合にSetPixel 等の描画関数を書いた後に入れる。
for(slice=0;slice<MAXSLICE;slice++){ Event( ) ;for(j=0; j<128; j++){ for(i=0; i<256; i++){fread( &data, sizeof(float), 1, fp ) ;y[ i ][ j ][ slice ] = data ;}}} fclose(fp);画像処理プログラム作成に興味を持った人、重要性に気付いた人、さらに学習したい人は後期に選択講義 医用画像処理演習 月曜２講目を受講して下さい。
卒業研究に５名程度、画像処理プログラム作成技術の習得を目的とした研究テーマを用意するので是非志願してください。
保健学科大学院には画像処理を研究する講義
機能画像解析学特論が用意されています。
卒業後、画像処理プログラムを学習し直したいと思ったら、いつでもこの講義ホームページを参照してください。質問も随時受付けます。