ボクセルファントムを用いた X線治療のシミュレーション

Presentation on theme: "ボクセルファントムを用いた X線治療のシミュレーション"— Presentation transcript:

1 ボクセルファントムを用いた X線治療のシミュレーション
PHITS Multi-Purpose Particle and Heavy Ion Transport code System ボクセルファントムを用いた X線治療のシミュレーション 2015年10月改訂 title 1

2 実習内容 DICOM2PHITSにより出力されるインプットファイルを ベースに、X線治療シミュレーションを実行するための 線源、コリメーターの設定方法についてまとめる。 コーン状のエネルギー分布をもつ光子ビームをDICOM2PHITSにより変換したボクセルファントムに照射するシミュレーション。 シミュレーションの全体図（左図）とボクセルファントム（右上図）および吸収線量分布（右下図）。 Contents 2

3 講習の流れ DICOM2PHITSの使い方 放射線治療シミュレーション DICOMデータの変換 全体像の把握 線源の設定 コリメーターの設定
はじめは、 /phits/utility/DICOM2PHITS/フォルダに 移動してください。 Table of Contents 3

4 １つのスライスに対するデータ（sample001.dcm）
Dicom形式（バイナリー） １つのスライスに対するデータ（sample001.dcm） ①　Header (撮影日時，ピクセルサイズなどの情報） ②　CT値（1,1→2,1→3,1→ … → nx-1, ny → nx, ny）の順番 このファイルがスライス数入ったフォルダで１つの物体を表現 １つのファイルを表示 フォルダ全体のデータを３D表示 Dicom形式からPHITS入力形式に変換する必要有 （CT値・バイナリ） （Universe番号・テキスト） DICOM2PHITS HowTo 4

5 変換プログラム(dicom2phits） DICOM2PHITS HowTo 5
CT値と物質密度・組成の関係はdata/HumanVoxelTable.data　[W. Schneider, Phys. Med. Biol. 45(2000) ]を参照 ①入力ファイルを作成（dicom2phits.inp）⇒詳しくはREADME参照 "data/HumanVoxelTable.data" 　変換テーブルの指定 "DICOM/" 　指定ディレクトリ内に含まれるDICOMファイルを自動判別 "PHITSinputs/" 　 PHITSインプットを格納するディレクトリを指定 　 指定した番号範囲のスライスファイルを読み込む (1<=z<=20) 　 DICOMデータの一部を切り出してボクセル化　(70<=x<=430, 90<=y<=460) 　 画像を粗くする（分解能を下げる）ことが可能 (x方向4個、y方向4個で平均) 　 原点オプション 0:ボクセル中心 1:DICOMヘッダーから抽出 　 PHITSパラメータ最適化オプション　0:最少 1:x線治療 2:粒子線治療 　 スライス読み込む順指定 +1:昇べき -1:降べき 原点オプション1では無視 ②実行 (Windows) dicom2phits_win.batにdicom2phits.inpをドラッグ&ドロップ 　　 (Mac) dicom2phits_mac.commnadをダブルクリック、現れる窓にdicom2phits.inpと入力 DICOM2PHITS HowTo 5

6 DICOM2PHITS実行結果 DICOM2PHITS outputs 6
PHITSのサンプルインプットphits.inpを指定ディレクトリ(PHITSinputs/)に生成 　　併せて、以下のインクルードファイルを同時に生成 　　　・material.inp 材質データを格納するファイル 　　　・universe.inp 各材質のユニバースデータを指定するファイル 　　　・voxel.inp PHITS形式に変換したボクセルデータを格納するファイル 　　　・libpath.inp 核データ、光子データのパスを指定するファイル 但し、libpath.inpは必要時（PHITSオプションの最適化1or2）且つ既に存在しない場合にfile(7)とfile(20)のサンプルパスを作成する。ユーザーの環境に合わせて変更する必要があるが、指定ディレクトリに一度作成すれば毎回変更する必要は無い。 DICOM2PHITS outputs 6

7 DICOM2PHITS output (phits.inp)
file=phits.inp [ Parameters ] … set: c81[ 90] $ number of x pixel set: c82[ 92] $ number of y pixel set: c83[ 20] $ number of z pixel set: c84[ ] $ unit voxel x set: c85[ ] $ unit voxel y set: c86[ ] $ unit voxel z [ Surface ] $ Unit voxel 5000 rpp -c87 -c87+c84 -c88 -c88+c85 -c89 -c89+c86 $ Outer region 99 so 500 $ Main Space 97 rpp -c87+c90 c87-c90 -c88+c90 c88-c90 -c89+c90 c89-c90 必ず１行目！ 　inflコマンドを使う時のおまじない [ C e l l ] $ Material universe infl:{universe.inp} $ Voxel universe lat=1 u=5000 fill= 0: 89 0: 91 0: 19 infl: {voxel.inp} $ Main space trcl=500 fill=5000 #97 $ Void $ Outer region ボクセルファントムのピクセル数や大きさはDICOMのヘッダーから自動的に設定 ボクセルファントムの中心がxyz座標の原点にくるように配置されている Include fileコマンド ファントムを回転・ 平行移動が可能 PHITS形式詳細は phits/lecture/voxel/phits-lec-voxel-jp.ppt 「Voxelファントムの作り方」を参照 DICOM2PHITS output (phits.inp) 7

8 CT値⇔物質密度・組成変換表 Conversion table 21 1行目：実行時に画面に表示するコメント 2行目：分割数
data/HumanVoxelTable.data　サンプル表 [W. Schneider, Phys. Med. Biol. 45(2000)459を参照] 1行目：実行時に画面に表示するコメント 2行目：分割数 Table for human voxel based on W. Schneider, Phys.… ! Number of universe di… d e ! Lowest CT value, Dens -950.0d ! Universe 2 -120.0d ! Universe 3 -82.0d ! Universe 4 NOTE: Den … -52.0d ! Universe [10^… -22.0d ! Universe [g/… … 1500.0d ! Universe 24 1600.0d ! Highest CT value for… #1 Air :: Air density is used ! Composition 14N ! Element, Ele… 16O ! 40Ar ! #2 Lung :: Lung density is used ! Composition o… 1H 12C 14N 3行目：物質１の定義 最小CT値 物質密度 構成元素数 -1000  物質１ < -950 最後の物質の最大CT値 元素名 仕切り行：#で始まる必要あり 組成割合： >0 = 粒子密度比、< 0 ＝質量比 物質１の組成 仕切り行：#で始まる必要あり 物質１の最小CT値よりも小さい　⇒　ワーニングを出して物質１で代用 最後の物質の最大CT値よりも大きい　⇒　ワーニングを出して最後の物質で代用 Conversion table 21

9 ジオメトリの確認 icntl = 11 icntl = 8 CT3D.eps deposit-xy.eps Geometry check 9

10 ボクセル化する範囲の変更 変更後(phits.inp) Region specification 10 dicom2phits.inp
"data/HumanVoxelTable.data" ! File for conversion of human voxel data "DICOM/" ! DICOM file directory "PHITSinputs/" ! Filename for voxel include file ! Minimum slice number, Maximum slice number ! Clipping: Nxmin, Nxmax, Nymin, Nymax ! Coarse graining: Nxc, Nyc, Nzc ! Origin 0:Center 1:Lab ! PHITS parameter: 0:Minimal 1:PhotonTherapy … ! Z direction +1:forward or -1:backward 変更後(phits.inp) CT3D.eps deposit-xy.eps Region specification 10

11 xyzメッシュを用いた[t-deposit]の結果
線量計算結果 icntl = 0 unit = 2 ⇒ MeV/source 単位 unit = 0 ⇒ Gy/source 単位 deposit-xy.eps xyzメッシュを用いた[t-deposit]の結果 Dose calculation 11

12 講習の流れ DICOM2PHITSの使い方 放射線治療シミュレーション DICOMデータの変換 全体像の把握 線源の設定 コリメーターの設定
Table of Contents 12

13 Windowsデフォルトインストールの場合変更不要
DICOMデータの変換 ①　dicom2phits.inpを以下のように変更して実行 ②　libpath.inpを環境に合わせて変更 ③　phits.inpの[parameters]セクションにあるnegsを0とする （多数の物質番号を定義した場合、EGSは未対応） "data/HumanVoxelTable-KumamotoUniv.data" ! File for conversion of human voxel data "DICOM-KumamotoUniv/" ! DICOM file directory "PHITSinputs/" ! Filename for PHITS input ! Minimum slice number, Maximum slice number ! Clipping: Nxmin, Nxmax, Nymin, Nymax ! Coarse graining: Nxc, Nyc, Nzc ! Origin 0:Center 1:Lab ! PHITS parameter: 0:Minimal 1:PhotonTherapy 2:ParticleTherapy ! Z direction +1:forward or -1:backward phits.inp X線治療に最適な設定 計算時間短縮のために ボクセルを粗視化 [ Parameters ] … $ Optimized for Photon Therapy emin(2) = E-10 negs = 0 $ Library data paths for file(7) & file(14) infl:{libpath.inp} file(7) = c:/phits/data/xsdir.jnd file(20) = c:/phits/XS/egs Windowsデフォルトインストールの場合変更不要 Conversion of DICOM 13

14 巨大なボクセルデータをあらかじめバイナリー化して読込時間短縮！
読込の高速化 目的 PHITSでは一度インプットファイルを全てバイナリー化してから再読込 巨大なボクセルデータをあらかじめバイナリー化して読込時間短縮！ 手順 ①　[Parameters]セクションのivoxelを有効にする（cを消す） ivoxel = # LatticeのFill部分をバイナリー化としてfile(18)に出力させるオプション file(18) = voxel.bin # 出力するバイナリーファイルのファイル名 ②　PHITSを実行する　→ Binary file was successfully generated!!　 ③　ivoxel = 1に変更し，Latticeを定義するinflコマンドをコメントアウト ivoxel = 1　# LatticeのFill部分をfile(18)から読み込むオプション 高速化！ $ infl:{voxel1.inp} Conversion of DICOM 14

15 ジオメトリの確認 icntl = 8 deposit-xz.eps deposit-xy.eps Conversion of DICOM
どの領域（座標）にファントムが配置されているか確認しよう。 deposit-xy.eps Conversion of DICOM 15

16 講習の流れ DICOM2PHITSの使い方 放射線治療シミュレーション DICOMデータの変換 全体像の把握 線源の設定 コリメーターの設定
ここで、/PHITSinputs/フォルダを /phits/lecture/therapy/XrayTherapyCT/フォルダに コピーしてください。以後は、 /phits/lecture/therapy/XrayTherapyCT/PHITSinputs/ にあるインプットファイルを書き換えてください。 Table of Contents 16

17 課題1:[t-track]タリーの追加 [t-track]を用いて、放射線治療シミュレーションの全体を確認できるようにしましょう。
/XrayTherapyCT/フォルダにあるt-track.inpから、axis=xyとなっている方の[t-track] をコピー&ペーストする。 タリーする範囲をxに関して-50から50cm、yに関して-120から40cm、zに関して-20から20cmに拡大する。 nx=100, ny=160とする。 t-track.inp [ T - T r a c k ] ・ ・ ・ ・ ・ ・ x-type = 2 nx = 25 xmin = -25. xmax = 25. y-type = 2 ny = 25 ymin = -25. ymax = 25. z-type = 1 nz = 1 axis = xy xy断面を出力する方をコピーする。 100cm離れた場所に線源を設定したい。 [t-track] 17

18 課題1の答え合わせ [t-track]を用いて、放射線治療シミュレーションの全体を確認できるようにしましょう。 track_xy.eps
phits.inp track_xy.eps [ T - T r a c k ] ・ ・ ・ ・ ・ ・ x-type = 2 nx = 100 xmin = -50. xmax = 50. y-type = 2 ny = 160 ymin = ymax = 40. z-type = 1 nz = 1 axis = xy 100cm離れた場所に線源を設定したい。 [t-track] 18

19 講習の流れ DICOM2PHITSの使い方 放射線治療シミュレーション DICOMデータの変換 全体像の把握 線源の設定 コリメーターの設定
Table of Contents 19

20 課題2:線源の設定 ボクセルファントムの正面に100cm離れた位置から光子を照射するように[source]セクションを設定してみましょう。ただし、光子ビームはコーン状に拡がるようにしましょう。 線源の発生位置をx=0cm, y=-105cm, z=0cmとする。 dir=0.0, phi=90とする（phiは追加する）。 r0=0.0, dom=20とする（domは追加する）。 phits.inp 照射方向を+yにするにはどのようにdirとphiを設定すれば良いだろうか？ また、コーン状に拡がったX線ビームを模擬するにはどうすれば良いだろうか？ [ Source ] s-type = 1 x0 = y0 = z0 = z1 = r0 = dir = 1.0 proj = photon e0 = 20.0 icntl=0として確認してみましょう。 [source] 20

21 課題2の答え合わせ ボクセルファントムの正面に100cm離れた位置から光子を照射するように[source]セクションを設定してみましょう。ただし、光子ビームはコーン状に拡がるようにしましょう。 track_xy.eps phits.inp [ Source ] s-type = 1 x0 = y0 = z0 = z1 = r0 = dir = phi = dom = proj = photon e0 = 20.0 +yの方向にコーン状のビームが照射されている。 [source] 21

22 課題3:エネルギー分布をもつ線源 線源光子がエネルギー分布をもつように設定しましょう。 [source] 22
s-type=4とする（円柱形状のs-type=1にエネルギー分布を与える場合）。 e0=の行はコメントアウト e-type=1, ne=28とし、inflコマンドをもちいてspectrum.inpを読み込む。 /XrayTherapyCT/フォルダにあるt-product-s.inpから、[t-product] をコピー&ペーストする。 phits.inp [ Source ] s-type = 1 x0 = y0 = z0 = z1 = r0 = dir = phi = dom = proj = photon e0 = 20.0 e-type = 1 ne = 28 infl:{******.inp} スペクトル(spectrum.inp) icntl=0としてphitsを実行し、[t-product]の出力結果であるproduct_source.epsを確認してみましょう。 [source] 22

23 課題3の答え合わせ 線源光子がエネルギー分布をもつように設定しましょう。 product_source.eps [source] 23
phits.inp product_source.eps [ Source ] s-type = 4 x0 = y0 = z0 = z1 = r0 = dir = phi = dom = proj = photon $ e0 = 20.0 e-type = 1 ne = 28 infl:{spectrum.inp} ・ ・ ・ ・ ・ ・ [ T-Product ] file=product_source.out エネルギー分布をもつ光子線源が発生していることが確認できた！ [source] 23

24 講習の流れ DICOM2PHITSの使い方 放射線治療シミュレーション DICOMデータの変換 全体像の把握 線源の設定 コリメーターの設定
Table of Contents 24

25 課題4:コリメーターの配置 コリメーターを配置して約10cm×10cmの正方形の照射野を設定しましょう。 Collimator 25
x=0cm, y=-50.5cm, z=0cmの位置に、1辺5cmの正方形の穴が空いた1辺40cm、厚み1cmのコリメーターの板を配置する（2つのrppを組み合わせる）。 コリメーターの領域番号は100とし、物質はouter region(物質番号-1)とする（計算時間短縮のため）。 /XrayTherapyCT/フォルダにあるt-track.inpから、axis=xzとなっている方の[t-track] をコピー&ペーストする。 タリーする範囲をxに関して-35から35cm(nx=70)、yに関して-11から-10cm(ny=1)、zに関して-20から20cm(nz=40)に変更する。 40 cm 5 cm 40 cm 5 cm y= -10cmの辺りをタリーする phits.inp [ Surface ] ・ ・ ・ ・ ・ ・ rpp ・ ・ ・ ・ ・ ・ 10 rpp ** ** ** ** ** ** 11 rpp ** ** ** ** ** ** y= -50.5cmの位置 [ Cell ] ・ ・ ・ ・ ・ ・ #** ** ** Collimator 25

26 課題4の答え合わせ コリメーターを配置して約10cm×10cmの正方形の照射野を設定しましょう。 track_xy.eps
phits.inp track_xy.eps track_xz.eps [ Surface ] ・ ・ ・ ・ ・ ・ rpp ・ ・ ・ ・ ・ ・ 10 rpp 11 rpp [ Cell ] #100 y= -10cmにおける光子のフルエンス。 およそ10cm×10cmの照射野が作られている。 コリメーターによりビームの周辺部分はカットされている。 Collimator 26

27 課題5:コリメーターの調整 ボクセルファントムの肺の部分に約5cm×5cmの大きさで照射されるようコリメーターに空いている穴の位置と大きさを調整しましょう。 コリメーターにある正方形の繰り抜きをxに関して-7.5から-4cm、zに関して1.5から4.5cmにする。 出力する単位を[Gy/source]にするため、2つの[t-deposit]にあるunitを0に変更する。 phits.inp この領域に照射したい [ Surface ] ・ ・ ・ ・ ・ ・ rpp ・ ・ ・ ・ ・ ・ 10 rpp 11 rpp x= -10cm, z=5cmを中心とする位置に照射野を設定するには？ Collimator 27

28 課題5の答え合わせ ボクセルファントムの肺の部分に約5cm×5cmの大きさで照射されるようコリメーターに空いている穴の位置と大きさを調整しましょう。 phits.inp track_xz.eps [ Surface ] ・ ・ ・ ・ ・ ・ rpp ・ ・ ・ ・ ・ ・ 10 rpp 11 rpp [ t-deposit] unit = 0 axis = xy axis = xz 照射野が変わった Collimator 28

29 吸収線量の空間分布 ただし、ボクセルファントム中の吸収線量の分布を見ると、ほとんどの領域で値をもっておらず、統計量を増やす必要がある。
deposit-xz.eps 統計量を増やすと Collimator 29

30 まとめ DICOM2PHITSにより出力されるインプットファイ ルをベースに、X線治療シミュレーションを実行す るためのセッティングを行った。 エネルギー分布をもつコーン状の光子線源を設 定できるようになった。 コリメーターを調整し任意の照射領域を設定でき るようになった。 《休憩はさむ》 まとめ Summary 30

31 宿題 右図の赤枠の領域（およそx=10cm, z=-5cm）に約5cm×5cmの大きさでX線が照射されるように、コリメーター中の穴の位置や大きさを変更してみましょう。 きちんとした吸収線量の分布が得られるようmaxcas, maxbchを増やしてみましょう。その際、照射部分の統計誤差が10%以下（deposit-xz_err.epsを見て、該当領域が黄色）になるまで増やしましょう。 この領域に照射するには？ 宿題（3次元体系） Homework 31

32 宿題（解答例） deposit-xz.eps deposit-xz_err.eps Homework 32