補正アルゴリズムの提案②

概要

補正アルゴリズムの提案①で行った、目印をX線管に貼り付けてブレを検出する方法では、台座上の変動と、焦点の変動の仕方が違うことから(たぶん…)ほとんど画像の精度を向上することができませんでした。補正アルゴリズム②では台座の変動も含めた機器全体の変動を検出するために台座の上に目印を置き、ブレを補正することを考えます。

補正②

予備実験

予備実験としてボールベアリングのボールを被写体として通常撮影(約10分程度の撮影)と長時間撮影(約70分の撮影)を行いその軌跡を調べました。積算回数はどちらも16回となっています。

補正②予備実験

結果を比較すると通常撮影を行った場合、最終的に元の場所に目印が戻ってきたのを確認できたのに対して、長時間撮影を行った場合は軌道が大きくずれていることが分かりました。これが長時間撮影による画質低下の原因になっているブレであると考えることができると思います。

理想的な軌道を求める

幾何学的な計算を行うことで、球体の理想的な軌跡を求めることができるので、その算出を試みました。目印(金属球)の回転中心までの距離rは実際の軌跡のデータから算出して、0.96mmとしています。回転角θが0の時のPの座標を算出すると(-0.96, 1.37, -0.04)となりました。この座標をθ度ずつ回転させた座標をそれぞれ投影していきます。

補正②計算①被写体である金属球を通り、焦点と検出器をつなぐ直線の方程式を求める。

$$\frac{x-y_1}{x_1-b}=\frac{y-y_1}{y_1}=\frac{z-z_1}{z_1}$$

②投影面と交わる点の座標を代入して整理して、投影面の座標を求める。

$$y”_1=y_1(1+\frac{l}{a-x_1})$$

$$z”_1=z_1(1+\frac{l}{a-x_1})$$

③点Pの最初の座標は()でこの座標をY軸を固定してθ度ずつ回転させていきます。そうするとその座標は・・・

$$(x_1,y_2,z_1)=(x_1,y_1,z_1)\left(\begin{array}{ccc}cos\theta\;\;\;0-sin\theta\\0\;\;\;\;\;1\;\;\;\;\;0\\sin\theta\;\;\;0\;\;\;cos\theta\end{array} \right)$$

$$=(x_1cos\theta-z_1sin\theta, y_1, x_1sin\theta+z_1cos\theta)$$

④検出器に投影される座標は

$$(a,y_1(1+\frac{l}{a-x_1}),z_1(1+\frac{l}{a-x_1}))$$

⑤上の式に回転後の座標を代入すると・・・

$$(a,y_1(1+\frac{l}{b-x_1cos\theta+z_1sin\theta}),(x_1sin\theta+z_1cos\theta)(1+\frac{l}{b-x_1cos\theta+z_1sin\theta}))$$

この式が回転後の座標になります。台座上では回転中心からの距離はrで、高さは透視画像から算出しました。rについてはとりあえず1㎜として計算してみました。その結果を示します。

補正②理想

今回求めた結果楕円状の軌跡を算出することができました。rは適当ですが結構実験値と近い結果になりました。ですが、金属球は遠くなったり近くなったりするわけで、ここまできれいになることは考えにくいと思います。ですので計算方法をもう一度考え直してみたいと思います。

実験

被写体には補正アルゴリズム①で使用したものと同じ直径0.4mmのドリルの刃を使用します。目印にはボールペンのボール(直径0.5mm)を使用します。ボールペンのボールは大変苦労しましたが、ラジオペンチと万力を使用してペン先を押しつぶして取り出しました。当初はキャップを作って被写体と目印を同時撮影するつもりでしたが、時間もないので消しゴムの中にボールを埋め込んでドリルに突き刺して被写体を作りました。

IMAG0610

透視撮影を行うと以下のような透視画像が得られます。この透視画像800方向分のデータを使って、ボールの軌跡を調べドリルの刃の撮影中に生じたボケを補正していきたいと思います。

jjj

 金属球の重心の変動を調べてみると下図のような結果になりました。始点と終点で金属球の位置が大きく異なっていることがわかりました。これだけ変動していれば画像のボケも大きくなりそうです。

軌跡

サイノグラムを作成すると以下のようになります。消しゴムの領域が灰色になって表示されています。ファンパラ変換をすると長時間撮影の影響でサイノグラムが上と下でつながっていないのがわかります。ファンパラ変換とはパラレルビームで撮影を行っていると考えて作っていたサイノグラムのデータをファンビームでとったものに変換することだと思います。bugi

この画像を使ってフィルタ補正逆投影法で画像の再構成を行うと以下のようになります。拡大率22倍での撮影ではアーチファクトは出ていますが、画像のボケはあまりはっきり見えませんでした。もう少しぶれている気がしましたが…とりあえずこれらの画像データを使用することにします。bbbbg

ブレの補正

X線点光源ブレのを含む、X線CT装置全体の撮影中に生じる歪みによる画像のボケを補正する為に、まずはサイノグラムをブレの大きさに応じてY方向、Z方向に線形に移動させてみました。実験結果より最終的なブレはY方向に5.3px、Z方向に6.1pxです。X線点光源のブレの測定実験でも台座の上にのせた被写体は、ほぼ線形に移動していることがわかっているので、Y方向、Z方向にサイノグラムを動かした結果下図右側のようなサイノグラムと再構成画像を得ることができました。

        hoseimae hoseigo

補正前と補正後で比較すると、補正前ではサイノグラムは始めと終わりでつながっていませんでしたが、処理を行うことでつながりました。再構成画像を見ても、若干画像がぼけたようにも見えますが、ドリルの刃から出ていたアーチファクトが少なくなったことを確認することができました。画像評価の一つの使用と考えている抽出できたサイノグラムの本数は、478本から486本と微増しました。

課題

抽出できたサイノグラムの本数だけでは画像を評価することは困難なので、画像を評価する方法を何か考えたい。MSE(Mean Square Error)やPSNR(Peak Signal to Noise Ratio)など?

また、金属球の理想的な軌道を求めるためには、金属球の座標を正確に(0.1ミリスケールくらい)知る必要があり、座標をできる限り正確に調べる必要がありそうです。それと金属球の軌跡は投影すると楕円にはならなそうなので計算方法も考え直す必要がありそうです。

X線CTの焦点の変動について

2次元的な焦点の移動

検出器に検出される目印の位置から、図のように相似の関係を使って焦点の変動を調べることができます。本研究で使用しているCTの仕様では検出器上の1ピクセルは0.2353mmに相当し、SFD(焦点から検出器までの距離)は463mmとなっています。

2dfocal

まずはピクセルをmmに変換してみます。検出器上の変位をdis[px]として、検出器上で0.6pixel移動した場合を考えます。

$$d = dis \times 0.2353 \times \frac{SOD}{SFD}$$

この式よりmmに変換ができます。変動が0.6px、SOD(被写体から焦点までの距離)が15mmだった場合、計算すると0.14㎜検出器上で変動していたことがわかります。ちなみにSFDをSODで割ることで拡大率を求めることができます。

次に相似の関係から焦点の変動を求めます。以下の式で求めることができいます。

$$\Delta=\frac{d \times SOD}{SFD – SOD}$$

計算すると焦点の変動は0.0047mmであったことがわかります。この変換式を使うことで、2次元的な焦点の変動を求めることができます。

2d移動

3次元的な焦点移動

これまでに2次元的な焦点の移動を調べましたが、図のように被写体の2つの点(P1,P2)を使って3次元的な変動も調べることができると思います。

3dfocal

焦点の初期位置を原点として、2次元的に焦点が移動した点D1,D2と被写体のP1,P2を通る直線の方程式を求めます。

$$\frac{x-a}{-a}=\frac{y-b}{d_y-b}=\frac{z-c}{d_z-c}$$ ・・・①

$$\frac{x-a”}{-a”}=\frac{y-b”}{d”_y-b”}=\frac{z-c”}{d”_z-c”}$$ ・・・②

3次元的な焦点の位置は①、②式の交点であると考えられるので、①、②式を媒介変数s,tを使って解きます。媒介変数を使って式を変換することで以下の6つの式を求めることができます。

$$x = -as +a$$ ・・・③

$$y = (d_y – b)s + b$$ ・・・④

$$z = (d_z – c)s + c$$ ・・・⑤

$$x = -a”s +a”$$ ・・・⑥

$$y = (d”_y – b”)s + b”$$ ・・・⑦

$$z = (d”_z – c”)s + c”$$ ・・・⑧

③~⑧の6つの式を同時に満たすs,tを求めます。とりあえずtの値を求めてみると…

$$t = \frac{a”d_y – ad_y – a”b + ab”}{-ad”_y + a”b + a”d_y -ab}$$

求めたtを⑥⑦⑧式に代入すると三次元的な焦点の座標を求めることができます。

$$(x,y,z) = (-a”t+a”,(d”_y-b”)t+b”,(d”_z-c”)t+c”)$$

導出した式を使って、3次元的な焦点の変動を算出してみると図のような変動をしていることがわかりました。

3d移動

本研究で使用している焦点の大きさは5μmであることから、焦点の変動自体は非常に小さいことがわかりました。奥行き方向(X方向)のぶれもほとんどありませんでした。このことから画像のボケの原因は、やはり焦点の変動以外の要因が大きいことが考えられます。

X線CTの精度調査のための撮影実験について

実験目的

X線CT装置では高精度な撮影を行うために積算回数や投影方向数を増やすことで、長時間のCT撮影を行うことがあります。長時間のX線CT撮影による機器内の温度上昇などから、X線発生装置の点光源にぶれが生じ、画質の精度低下が起きているのかを確かめるために、CTの撮影実験を行いました。今回は通常の撮影方法で撮影した画像と長時間かけて撮影した画像でどのような違いが生じるのかを調べています。

実験手順

  1. 被写体として直径0.4mmのドリルの刃を使います。被写体のドリルの刃を台座に固定し、拡大率を25倍程度になるように台座を移動させて拡大率を調整します。
  2. 今回は以下の3つの撮影条件で撮影を行いました。
    ①積算回数16回、撮影時間10分
    ②積算回数128回、撮影時間70分
    ③積算回数16回、撮影時間70分
  3. それぞれの条件で撮影した透視画像のデータをConeCtExpressという画像再構成ソフトを使って被写体の断層画像を作成します。
  4. Molcer Plusというソフトを使って3次元モデルを作成します。

実験結果

dril

3つの条件で撮影を行った結果、上図のような結果を得ました。どれもドリルの先端部分を撮影した結果ですが、通常の撮影方法では画像にボケが生じていなかったのに対して、長時間撮影を行った右の二つはボケが生じていたことを確認できました。特に積算回数を増やさずにただ単に長時間撮影を行った真ん中の画像では、大きなボケが生じていることを確認できました。3次元的にモデリングを行った結果を以下に示します。

三次元モデル

3次元モデルを見ると、ボケが生じていたのにも関わらず積算回数を128回に増やして撮影したモデルが一番きれいなモデルになっていることがわかりました。

画像の定量的評価

3つの再構成画像の定量的な評価を行いました。ドリルの複雑な形を評価するのは困難であるので、ひとまず、円形のドリルの柄を使って定量的な評価を試みました。具体的には拡大率のわかっている画像に対して、ドリルの領域を2値化し理想的な再構成画像との差分をとっています。

public void evaImage(){
  int cnt=0, centerX=0, centerY=0;
  int width = 640;

  //被写体の重心を求める
  for(int i = 0; i < width; i++){
    for(int j = 0; j < width; j++){
      //既に2値化した再構成画像に対して、被写体の座標の合計値を求める
      if(image[i*width+j]!=0){
        centerX+=i;
        centerY+=j;
        cnt++;
      }
    }
  }
  //重心の座標を求める
  centerX/=cnt;
  crntery/=cnt;

  //理想的な再構成画像を作成する
  for(int i = 0; i < width; i++){
    for(int j = 0; j < width; j++){
      if((i-centerX)*(i-centerX)+(j-centerY)*(j-centerY) < 16*16){

        //半径が16pixel未満だったら黒
        difImage[i*width+j] = 255;
      }else{
        //それ以外の場合は白
        difImage[i*width+j] = 0;
      }
    }
  }

  //画像の差分を求める
for(int i = 0; i < width; i++){
    for(int j = 0; j < width; j++){
      //再構成画像と数値ファントムの差分が0でない場合
      if(difImage[i*640+j]-PanelDraw.image2[i*640+j]!=0){
        checkDiff[i * width + j] = 255;
        cnt2++;
      }
    }
  }

  //画像評価
  num = cnt2 / (width*height) * 100;
  System.out.println(num);

画像の評価を数値化すると、以下のような結果となりました。定性的な評価を行ったときと同様に積算回数を増やさずに長時間撮影を行った場合が一番低く(ボケが多い)、積算回数を増やして長時間撮影を行ったとき一番高い(理想的な形に近い)結果になりました。

画像評価

以上の結果から、定性的にも定量的にも長時間撮影を行うことによって画像にボケが生じることを確認することができました。しかしこの手法では評価できる形状が限定されるため、ほかの汎用性の高い方法を考える必要があります。本研究ではこのボケを画像データ上で補正することで画質の向上を行おうと考えています。

透視撮影実験について

実験目的

X線CTの精度評価のための撮影実験において、定性的にも定量的にも長時間のCT撮影においてボケが生じ画質が低下していることを確認することができました。本実験ではCT撮影は行わず、60分間1分間隔で目印の透視撮影を行うことで、どの程度焦点が変動しているかを定量的に調べることを目的としています。

透視撮影:台座を動かさずに1方向からX線を照射し撮影を行うもので、レントゲン撮影のような画像が得られます。

実験方法

本実験では以下の2つの目印を使って焦点の変動を調べます。ターゲットAは真鍮製の板に銅線を交差するように配置したもので、ターゲットBはアクリル製の樹脂に同様に銅線を配置したものです。

ターゲットA

ターゲットB

ターゲットAはX線管のベリリウム窓に直接貼り付け、X線管内で起きている変動を計測します。一方ターゲットBはステージの上に置くことで焦点のブレ以外にも台座やCT装置自体の歪みによる変動も計測します。それぞれを以下のように配置します。また、本実験では温度変化が問題になっていると考えられるので、ターゲットAの表面温度とターゲットB付近の雰囲気温度を計測します。

実験風景

また、焦点のブレはエージングを行った直後に大きくなることが考えられるので、エージング直後とエージング数時間後の2つの条件で撮影を行いました。撮影条件はそれぞれ管電圧130Kvで管電流は20mAです。幾何学的拡大率は27.8倍で撮影を行いました。それぞれを1分間隔で60分撮影を行った結果を画像解析し変動を調べます。

エージング:CT装置の電源を付けた直後に行われるもので、X線管球のフィラメントの寿命を延ばすための暖機運転のようなものである。具体的には管電圧、管電流を徐々に上げて最大の線量になるまでX線を発生させています。通常10~15分かかります。

幾何学的拡大率:CT撮影を行う時の倍率のこと。被写体から焦点までの距離と検出器から焦点までの距離の関係で決まります。被写体が焦点に近いほど大きくなります。

解析方法

2つのターゲットを同時撮影すると下図のように2つのターゲットを同時撮影した透視画像を得ることができます。60枚の画像に対して同じ処理を行い、銅線で囲まれた部分の重心を求めていきます。

変動解析

  1. 探索する領域を決定する
  2. 領域内で閾値以下の画素を探索する
  3. 閾値以下の画素とそれ以外で2値化する
  4. 2値化した領域の重心の座標を求める
public static void searchGravity(){

  double sumGravityX;
  double sumGravityY;
  double gravityX; //重心のX座標
  double gravityY; //重心のY座標
  int cnt = 0;

  //探索する領域を順番に調べていく(場合に応じて調節)
  for(int i = 0; i < 640; i++){
    for(int j = 0; j < 480; j++){
      //注目した画素が閾値以下の場合
      if(rawData[i+j*640] < 30000){
        sumGravityX += j;
        sumGravityY += i;
        cnt++;
      }
    }
  }

  //重心の座標を求める
  gravityX = sumGravityX / cnt;
  gravityY = sunGravityY / cnt;
}

実験結果

まずは測定した温度変化についてです。2地点の温度変化を時系列に示しています。

エージング直後エージング数時間後

温度変化はエージング直後の実験を行った場合は、ターゲットAの表面温度も台座付近の雰囲気温度も変化が大きくなりました。どちらも温度は上がり続けることがわかりました。

次にそれぞれのターゲットの変動についてです。それぞれの条件でX線管に貼り付けたターゲットAと台座の上に乗せたターゲットBの変動について調べています。

エージング数時間後2エージング数時間後1

X線管に貼り付けたターゲットAはエージング直後でも数時間後でも同程度の変動をしていたことが分かりました。その一方で台座に乗せたターゲットBではエージング直後に大きく変動していたことが分かりました。このことから温度変化による影響を受けやすいのは台座にターゲットを乗せた場合でした。以上のことから焦点のブレよりも、撮影中の台座の変動や、CT装置自体の歪みなどが長時間撮影による画質低下の主因になっていることが考えられます。

焦点の移動による画質低下が起きている可能性は低いことから今後は台座の変動も含めた撮影中の機器の歪みによる画質低下を補正する方法を考えていきたいと思います。

X線CT装置の構成要素について

CT装置の種類について

X線CT装置には主に医療用のCT産業用のCTがあります。大きな違いは医療用のCTでは少ない被爆量で撮影を行うために、ヘリカルスキャンという螺旋状に撮影を行う技術が用いられていて撮影も数十秒で行うことができます。一方産業用のCTでは精密な検査を行うことが求められるので、強いX線を発生することができて数十秒から数十分と長時間の撮影も行うことができます。おもにコーンビームのX線で撮影が行われます。

ctsouti

医療用のCTでは患者さんはベッドに寝たままで、装置自体が回転して撮影を行っています。その一方で産業用CTは通常被写体を回転させて撮影を行っています。

CT装置の構成要素について

X線CT装置は、おおまかにX線を発生させるX線管装置とX線量を検出する検出器からできています。X線管装置と検出器の間に被写体を起き、被写体もしくはX線管と検出器を回転させながら撮影を行います。

X線管装置

X線管装置には主に密閉型開放型の二種類があります。開放型のX線管では焦点がX線の出る窓に近くなるのでその分高拡大率で撮影することが可能になります。対陰極としては銅やモリブデン、タングステンなどの金属に加速した電子ビームを当て、原子の1s軌道の電子を弾き飛ばします。空になった1s軌道に外の軌道の電子が遷移してきます。この遷移によって放出される電磁波をX線特に特性X線と呼んでいます。本研究で使用しているX線管は以下のような仕様となっています。

Xray tube pic

検出器

X線検出器にはパルス型積分型の検出器があります。パルス型の検出器にはシンチレーションカウンタのように、X線光子を1個1個計測するタイプであり、検出器が持つノイズはほとんど無視することができ、微弱なX線強度を精度よく検出することができます。一方積分型はX線エネルギーを可視光強度や電流などの物理量に変換して検出しています。X線フィルムやX線CCDがこのタイプです。本研究で使われている検出器はCCD型で二次元的にフォトダイオードが画素として並んだ撮像素子です。

ditector

CTスキャンの仕組み

X線管で発生させたX線を被写体に照射し、検出器で検出器に到達したX線強度を測定します。被写体を通過したX線は吸収され検出されるX線強度は小さくなっています。このデータ1つ1つを透視画像と呼びます。この透視画像を被写体の周囲方向から通常800方向程度のデータを取得し、得られたデータを使い画像を再構成することで物体の断面画像を得ています。

—————————————–参考文献———————————————-

「Wikipedia X線」

「二次元X線検出器の原理と性能,  雨宮,伊藤,  日本結晶学会誌, 45, p163−170, 2003」

X線CTの単純逆投影について

CTの勉強会第2回目は単純逆投影を行いCTの再構成画像をつくりました。CT撮影によって得られた800方向分の透視画像を縦に並べたサイノグラムと呼ばれるデータを使って、再構成画像を作っていきます。単純逆投影は最も簡単な投影方法で誤差が大きく生じる方法です。

1.順投影と逆投影について

まず順投影は下図のように被写体の周囲方向からX線を当てて、被写体のデータを取得する方法です。通常800方向分のデータを取得します。逆投影は得られた周囲方向のデータを使って逆にデータを収めていく方法です。例えば左の{3,12,3}という並びの投影データを収めると「3」が入っていた行はすべて「3」入ります。「12」が入っていた行はすべて「12」となります。これを取得したデータすべてで行い、そのつど加算していきます。最後に撮影方向数でそれぞれを割ると、再構成画像ができます。単純逆投影だとかなり誤差が大きくなり、通常ですとフィルタをかけることでより鮮明な再構成画像を作っています。

project

2.sng形式のサイノグラムを読み込む

前回の勉強会の時にrawデータを読み込んだのと同様にサイノグラムを読み込みます。サイノグラムの大きさは640×800pixelとなっているのでそこだけ気をつけます。画像を読み込むプログラムは第一回勉強会のページを参考にしてください。今回はモーターを撮影したデータを扱います。最終的に透視画像の中心部分の断面画像を作成します。

スクリーンショット 2013-10-14 13.54.24

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 />

スクリーンショット 2013-10-14 13.55.05

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[]の中に再構成画像のデータを格納することができたので、投影データを表示した時と同じ方法を使って、再構成画像を表示します。これで一番簡単な逆投影ができるようになりました。全体的にボケが多い再構成画像になっています。よりはっきりした画像を得るためにはフィルタをかけるなど難しい処理がまだ必要になるようです。

スクリーンショット 2013-10-23 21.20.13

X線CTの撮影実験

X線CTのX線点光源のブレが被写体に与える影響を調べるために,微小な被写体を拡大率25倍で撮影を行った.本実験では直径0.4㎜のドリルの刃を,マイクロフォーカスX線CT(ScanXmate-A130ss940, コムスキャンテクノ株式会社)を使用しで撮影を行った。

CT
X線CT装置

drildril
被写体のドリルの刃

長時間のCT撮影を行うことで点光源のブレが大きくなると考えられるので,以下の3つの撮影条件で撮影を行った.各撮影方向で取得するデータ数である積算回数を増やすことによって撮影時間を増やす.

既存研究について

X線CT装置の焦点が撮影中に変動する現象はいくつかの著書で紹介されています[1]。ですが、焦点ブレを補正する研究がおこなわれている例はほとんどなく、私が知る限りでは日立メディコが行っている、機械的に焦点移動を起こしにくくする研究のみです。その研究の概要を紹介します。

 

その一方で画像データ上で焦点移動を補正する研究は私の知る限りでは存在せず、その有用性は大きいと考えています。

————————————–参考文献————————————–

「X線管装置の焦点移動量低減技術の開発」, MEDIX, VOL.36, pp32-36

クイックソートとその計算量について

【編集中】

数列を昇順または降順に並べるアルゴリズムの一つにクイックソートがあります。直感的に分かりやすいバブルソートと比較してピボットによる分割があったりでめんどくさいイメージがあったのですが、有名なアルゴリズムなので勉強しておきたいと思います。以下計算手順です。テストケースとして次の数列を扱います。

$$A={4,8,6,5,2,1,3,9,7}$$

  1. 基準値を決めます。たとえば、最初と最後の要素の平均にします。
  2. 前方から基準値より大きな数を、後方から基準より小さな数を探します。見つかったら入れ替えます。
  3. 得られた列はぶつかったところを境に、基準値よりも小さいグループと大きいグループに分かれています。そして、この得られた列それぞれに同じ手順を繰り返していきます。

qick

この操作を繰り返すことで数列を昇順に並び替えることができます。それではこの操作をC言語で記述してみます。


#include <stdio.h>
int main(void){

return 0;
}

CT画像の読み込み方書き出し方

第一回の勉強会では大学で使用しているマイクロフォーカスX線CT(ScanXmate-A130ss940, コムスキャンテクノ株式会社)で採用されているRAWデータの画像を読み書きする方法を勉強しました。今回は最終的にrawデータを読み込んで、リアルタイムで閾値処理をするプログラムを作成しました。

rawデータの読み込み方

今回扱うrawデータは640×480(307200)pixelですが、データの大きさは614400byteとピクセルの情報の2倍の大きさとなっています。これは2byteで1pixel分の情報を表していることを意味していてこの画像が16bitのrawデータであることがわかります(1byteが8bitです)。

data16bitの画像を読み込むためには、2byte分の情報から1pixelの情報を知るために特別な処理が必要となります。画像にはビッグエンディアンリトルエンディアンという記録方法の方式があります。例えば{1234ABCD}という4バイトのデータをデータの上位バイトからメモリに{12 34 AB CD}と並べる方式をビッグエンディアン、{CD AB 34 12}とデータの下位バイトから並べる方式をリトルエンディアンと言います。それぞれデータの読み込み方が異なってくるので注意が必要です。

$$d_1=b_1+b_2\times256$$

もしくは

$$d_1=b_1\times256+b_2$$

となります。

public static void readRawData() {
  File f = new File("○○○.raw");
  try {
    DataInputStream dis = new DataInputStream(new FileInputStream(f));
    //データの大きさだけバイト型の配列を用意
    byte[] b = new byte[dis.available()];
    //データを読み込む
    dis.read(b);
    //データの形式に合わせてint型の配列に画素地を入れる(今回は16bit)
    for (int i = 0; i < b.length / 2; i++) {
      rawdata[i] = (0xff & b[i * 2]) + (0xff & b[i * 2 + 1]) * 256;
    }
  } catch (Exception e) {
    e.printStackTrace();
  }
}

上記のように記述することによって16bitのrawデータを配列に格納することができました。配列に数値を格納できれば、閾値処理を行ったり、平滑化フィルタをかけたりなど様々な画像処理を行うことができるようになります。画像は画像の形式によって配列に格納する方法が変わってきます。例えばTIFF形式の画像ファイルの場合、ヘッダやフッタ(データの最初と最後)に画像に関する情報が入っている場合があり、このような画像では単純に配列に格納するだけでは読み込むことができません。

rawデータの表示方法

続けて表示方法です。描画を行うクラスを作成し、メインクラスでインスタンスを生成し、描画を行うます。NetBeansの場合はパッケージを右クリック→新規→JPanelフォームをクリックし、描画を行うパネルを作ります。そうすることでJPanelクラスをスーパークラスとするサブクラスが自動生成されます。自動生成したクラスの中に以下のコードを記述します。

public void paintComponent(Graphics g){
  //背景を白く塗りつぶす
  g.setColor(Color.white);
  g.fillRect(0,0,width,heigth);
  //画像データを描画する
  for(int i = 0; i < height; i++){
    for(int j = 0; j < width; j++){
      //65536諧調になっているので256諧調に変換
      int c = rawdata[i * width + j]/256;
      //色をセット(グレースケールなのでRGBすべてに同じ値を入れる)
      g.setColor(new Color(c,c,c));
      //1pixelずつ色を塗る
      g.drawLine(j, i, j, i);
    }
  }
}

上記のコードを記述し呼び出すことで、rawdata[]に格納していた画像のデータを描画することができます。次に呼び出し方についてです。

initComponentsメソッドが記述されているメソッドに以下のコードを書き足します。ちなみにinitComponentsメソッドはコンポーネントの作成に関するコードが記述されていてプログラムの中で一番最初に呼び出される部分になっています。

public StudyCt() {
  initComponents();
  //インスタンスの生成
  PanelDraw pDraw = new PanelDraw();
  lbl1.setLayout(new BorderLayout());
  //生成したパネルをラベルの上に張り付ける
  lbl1.add(pDraw);
  //外観を成形する
  pack();
}

これでおおまかな部分は完成です。後は閾値を変更した後、ボタンを押した後など随所でrepaintメソッドを呼び出すことで、再描画を行うことができます。

rawデータの書き出し方

16bitのrawデータとして画像を書き出します。読み込み方と方法は似ていて、おまじないの部分が多くなってしまいますが、以下のように記述します。理由は不明ですが書き出すときはsng形式で書き出すようです。

public void writeRaw(){
  File f = new File("○○○.sng");
  int bit = 16;
  try{
    byte b1[] = new byte[width * height * bit / 8];
    for(int i = 0; i < height; i++){
      for(int j = 0; j < width; j++){
        int buf = 0;
        if(rawdata[j][i]>65535)buf = 65535;
        else if(rawdata[j][i] < 0)buf = 0;
        else buf = (int)rawdata[j][i];
        //1pixelのデータを2byte分に分けて配列に格納する
        b1[(i * 640+j)*2] = (byte)(buf%256);
        b1[(i * 640+j)*2+1] = (byte)(buf/256);
      }
    }
    FileOutputStream fos0 = new FileOutputStream(f);
    fos0.write(b1,0,widht*height*2);
    fos0.close();
  }catch(IOException e){
  }
}

上記のように記述することで○○○.sngという画像ファイルを作成することができます。だいたいの部分が読み込むときと逆の操作を行っているのがわかると思います。配列に格納するときは条件分岐させて0以下65536以上の情報を持っているpixelがないようにしています。配列に格納できたら書き出してファイルを閉じます。これで完成です。

今回の勉強会で勉強したことを使って簡単なソフトを作ってみました。CTの画像を読み込んで、閾値処理を行い形状と拡大率がわかっている被写体に対して理想的な再構成画像とどれだけ違いがあるかを調べることができます。勉強会大変勉強になりました。

スクリーンショット 2013-10-12 11.40.51