CTの勉強会第2回目は単純逆投影を行いCTの再構成画像をつくりました。CT撮影によって得られた800方向分の透視画像を縦に並べたサイノグラムと呼ばれるデータを使って、再構成画像を作っていきます。単純逆投影は最も簡単な投影方法で誤差が大きく生じる方法です。
1.順投影と逆投影について
まず順投影は下図のように被写体の周囲方向からX線を当てて、被写体のデータを取得する方法です。通常800方向分のデータを取得します。逆投影は得られた周囲方向のデータを使って逆にデータを収めていく方法です。例えば左の｛3,12,3｝という並びの投影データを収めると「3」が入っていた行はすべて「3」入ります。「12」が入っていた行はすべて「12」となります。これを取得したデータすべてで行い、そのつど加算していきます。最後に撮影方向数でそれぞれを割ると、再構成画像ができます。単純逆投影だとかなり誤差が大きくなり、通常ですとフィルタをかけることでより鮮明な再構成画像を作っています。
2．sng形式のサイノグラムを読み込む
前回の勉強会の時にrawデータを読み込んだのと同様にサイノグラムを読み込みます。サイノグラムの大きさは640×800pixelとなっているのでそこだけ気をつけます。画像を読み込むプログラムは第一回勉強会のページを参考にしてください。今回はモーターを撮影したデータを扱います。最終的に透視画像の中心部分の断面画像を作成します。
3．投影データを作成する
sino[]という配列にデータを格納できたら、rawデータのままでは再構成画像を作ることができないので、対数を取って投影データを作成します。投影データは以下の数式で表せれます。但し、$$I_o$$が吸収体がないときの入射強度(65536)で、$$I_y$$が透過強度を表しています。
$$I = ln(\frac{I_0}{I_y})$$
これはX線に吸収体の厚み(L)が大きくなると、指数関数的にX線が透過しにくくなる性質があり、再構成する際には線形に戻す必要があるからです。詳しくは理解できたら追記します。
$$I = I_0e^{\sum_{m=1}^{n}\mu_ml_m}$$
それと投影データを作る際に、スケーリングを行うために同時に投影データの最大値と最小値を求めておきます。
<br />public static void change(){<br /><%%KEEPWHITESPACE%%> for(int i = 0; i < height; i++){<br /><%%KEEPWHITESPACE%%> for(int j = 0; j < width; j++){<br /><%%KEEPWHITESPACE%%> //対数をとる<br /><%%KEEPWHITESPACE%%> sinoP[i * width + j]= - Math.log(sino[i * width + j]/65536.0);<br /><br /><%%KEEPWHITESPACE%%> //最大最小を求める<br /><%%KEEPWHITESPACE%%> sinoPmax=Math.max(sinoP[i * width + j], sinoPmax);<br /><%%KEEPWHITESPACE%%> sinoPmin=Math.min(sinoP[i * width + j], sinoPmin);<br /><%%KEEPWHITESPACE%%> }<br /><%%KEEPWHITESPACE%%> }<br /><%%KEEPWHITESPACE%%> repaint();<br />}<br />投影データの描画方法です。前回同様に描画を行うクラスを作成し、その中に記述していきます。理由はよくわかりませんがrawデータを単純に表示するときと違いスケーリングを行う必要があります。rawデータでは65536諧調のデータだったので256で割っていましたが、今回は画素値の最大値から最小値を引いたものに256をかけてからそれぞれの画素値を割っています。こうすることで投影データのサイノグラムを表示することができます。
<br />public void paintComponent(Graphics g) {<br /><%%KEEPWHITESPACE%%> int width = 640; //幅<br /><br /><%%KEEPWHITESPACE%%> //背景を白く塗りつぶす<br /><%%KEEPWHITESPACE%%> g.setColor(Color.white);<br /><%%KEEPWHITESPACE%%> g.fillRect(0, 0, width, width);<br /><br /><%%KEEPWHITESPACE%%> for (int i = 0; i < width; i++) {<br /><%%KEEPWHITESPACE%%> for (int j = 0; j < width; j++) {<br /><%%KEEPWHITESPACE%%> //スケーリングを行う<br /><%%KEEPWHITESPACE%%> int c = (int)(sbp[j+i*width]/(sbpmax - sbpmin)*256);<br /><br /><%%KEEPWHITESPACE%%> g.setColor(new Color(c, c, c));<br /><%%KEEPWHITESPACE%%> g.drawLine(j, i, j, i);<br /><%%KEEPWHITESPACE%%> }<br /><%%KEEPWHITESPACE%%> }<br />}<br />
4．単純逆投影を行う
単純逆投影は順投影とは逆にwidth(640px)×height(640px)の領域に800方向からのデータを投影していくという説明がよくされますが、通常はwidth×heightの領域を回転させ、その領域に投影データを入れていくという手法を取る方が実装が簡単なようです。投影データを入れていく際にピクセルをまたいでしまう場合があり、それはバイリニア法という画像処理方法を使って線形に補完して数値を入れていきます。最後にまた再構成画像の最大値と最小値を求めます。プログラムは以下のように記述することができます。
<br />public static void backProjection(){<br /><%%KEEPWHITESPACE%%> int proj = 800; // プロジェクション数<br /><%%KEEPWHITESPACE%%> int center = width/2;<br /><%%KEEPWHITESPACE%%> double step = 360.0/(double)proj; //ステップ角度（度）<br /><%%KEEPWHITESPACE%%> double sin, cos;<br /><%%KEEPWHITESPACE%%> int x,y;<br /><%%KEEPWHITESPACE%%> double bx,by;<br /><br /><%%KEEPWHITESPACE%%> for(int p = 0; p < height; p++){<br /><%%KEEPWHITESPACE%%> //現在の角度を計算<br /><%%KEEPWHITESPACE%%> theta = p * step; //+=stepより安全<br /><br /><%%KEEPWHITESPACE%%> //あらかじめ算出（角度はラジアン表記）<br /><%%KEEPWHITESPACE%%> sin = Math.sin((theta+180)*Math.PI/180);<br /><%%KEEPWHITESPACE%%> cos = Math.cos(theta * Math.PI /180);<br /><br /><%%KEEPWHITESPACE%%> for(int i = 0; i < width; i++){<br /><%%KEEPWHITESPACE%%> for(int j = 0; j < width; j++){<br /><br /><%%KEEPWHITESPACE%%> //再構成画像上のピクセルを回転 各方向から投影していく方法は難しい<br /><%%KEEPWHITESPACE%%> bx = (int)((i - center) * cos - (j - center) * sin );<br /><%%KEEPWHITESPACE%%> by = (int)((i - center) * sin + (j - center) * cos );<br /><br /><%%KEEPWHITESPACE%%> //回転後再構成領域からはみ出していたら次の点に進む<br /><%%KEEPWHITESPACE%%> if(bx*bx + by*by >= center*center){<br /><%%KEEPWHITESPACE%%> continue;<br /><%%KEEPWHITESPACE%%> }<br /><br /><%%KEEPWHITESPACE%%> bx += center;<br /><%%KEEPWHITESPACE%%> by += center;<br /><br /><%%KEEPWHITESPACE%%> x = (int)bx;<br /><%%KEEPWHITESPACE%%> y = (int)by;<br /><br /><%%KEEPWHITESPACE%%> double bbx0,bbx1,bby0,bby1,s00,s10,s01,s11;<br /><%%KEEPWHITESPACE%%> bbx0 = bx - x;<br /><%%KEEPWHITESPACE%%> bbx1 = 1 - bbx0;<br /><%%KEEPWHITESPACE%%> bby0 = by - y;<br /><%%KEEPWHITESPACE%%> bby1 = 1 - bby0;<br /><br /><%%KEEPWHITESPACE%%> s00 = bbx0 * bby0;<br /><%%KEEPWHITESPACE%%> s10 = bbx1 * bby0;<br /><%%KEEPWHITESPACE%%> s01 = bbx0 * bby1;<br /><%%KEEPWHITESPACE%%> s11 = bbx1 * bby1;<br /><br /><%%KEEPWHITESPACE%%> if(x >= 0 && x <= width -1 && y >= 0&& y <= width -1){<br /><%%KEEPWHITESPACE%%> sbp[x + y * width] += sino[i + width * (int)(theta/step)]*s11;<br /><%%KEEPWHITESPACE%%> }<br /><%%KEEPWHITESPACE%%> if(x+1 >= 0 && x+1 <= width -1 && y >= 0&& y <= width -1){<br /><%%KEEPWHITESPACE%%> sbp[x + y * width] += sino[i + width * (int)(theta/step)]*s01;<br /><%%KEEPWHITESPACE%%> }<br /><%%KEEPWHITESPACE%%> if(x >= 0 && x <= width -1 && y+1 >= 0&& y+1 <= width -1){<br /><%%KEEPWHITESPACE%%> sbp[x + y * width] += sino[i + width * (int)(theta/step)]*s10;<br /><%%KEEPWHITESPACE%%> }<br /><%%KEEPWHITESPACE%%> if(x +1>= 0 && x+1 <= width -1 && y +1>= 0&& y +1<= width -1){<br /><%%KEEPWHITESPACE%%> sbp[x + y * width] += sino[i + width * (int)(theta/step)]*s00;<br /><%%KEEPWHITESPACE%%> }<br /><%%KEEPWHITESPACE%%> }<br /><%%KEEPWHITESPACE%%> }<br /><%%KEEPWHITESPACE%%> }<br /><%%KEEPWHITESPACE%%> //再構成画像の最大値最小値を求める<br /><%%KEEPWHITESPACE%%> for(int i = 0; i < width; i++){<br /><%%KEEPWHITESPACE%%> for(int j = 0; j < width; j++){<br /><%%KEEPWHITESPACE%%> spbmax = Math.max(sbp[i*width+j], spbmax);<br /><%%KEEPWHITESPACE%%> spbmin = Math.min(sbp[i*width+j], spbmin);<br /><%%KEEPWHITESPACE%%> }<br /><%%KEEPWHITESPACE%%> }<br />}<br />これでsbp[]の中に再構成画像のデータを格納することができたので、投影データを表示した時と同じ方法を使って、再構成画像を表示します。これで一番簡単な逆投影ができるようになりました。全体的にボケが多い再構成画像になっています。よりはっきりした画像を得るためにはフィルタをかけるなど難しい処理がまだ必要になるようです。
