P HI T S

1. PHITS Multi-Purpose Particle and Heavy Ion Transport code System PHITS 講習会 放射線治療計画シミュレーション 2014年5月改訂 title 1

2. 目的 治療計画シミュレーションの検証をPHITSを用いて行う 後方照射155MU 前方照射86MU Introduction 2

3. 講習の流れ • DICOM2PHITSの使い方 • 放射線治療シミュレーション • DICOMデータの変換 • ビームの設定 • タリーの設定 • 校正定数の設定 • 治療計画シミュレーション Table of Contents 3

4. 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） • Dicom形式のデータをPHITS形式のボクセルデータに変換 • CT値と物質密度・組成の関係はdata/HumanVoxelTable.data[W. Schneider, Phys. Med. Biol. 45(2000)459-478]を参照 ①　基本情報ファイルを作成（dicom2phits.inp） ②　(Windows) dicom2phits_win.batにdicom2phits.inpをドラッグ&ドロップ (Mac) dicom2phits_mac.commnadをダブルクリックし、現れる窓にdicom2phits.inpと入力 指定ディレクトリ内に含まれるDICOMファイルを自動判別 指定した番号範囲のスライスファイルを読み込む "data/HumanVoxelTable.data" ! File for conversion of human voxel data "DICOM/" ! DICOM file directory 1 20 ! Minimum slice number, Maximum slice number 70 430 90 460 ! Clipping: Nxmin, Nxmax, Nymin, Nymax 4 4 1 ! Coarse graining: Nxc, Nyc, Nzc 0 ! Origin 0:Center 1:Lab 1 ! Z direction +1:forward or -1:backward 0 ! PHITS parameter: 0:Minimal 1:PhotonTherapy2:ParticleTherapy "PHITSinputs/" ! Filename for PHITS input DICOMデータの一部を切り出してボクセル化できる PHITSインプットを格納するディレクトリを指定 原点情報をDICOMヘッダーから抽出できる PHITSオプションの最適化を選択 画像を粗くする（分解能を下げる）ことも可能 DICOM2PHITS HowTo 5

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

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

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

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

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

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

12. 講習の流れ • DICOM2PHITSの使い方 • 放射線治療シミュレーション • DICOMデータの変換 • ビームの設定 • タリーの設定 • 校正定数の設定 • 治療計画シミュレーション Conversion of DICOM 12

13. DICOMデータの変換 ①　dicom2phits.inpを以下のように変更して実行 ②　libpath.inpを環境に合わせて変更 "data/HumanVoxelTable-KumamotoUniv.data" ! File for conversion of human voxel data "DICOM-KumamotoUniv/" ! DICOM file directory 1 148 ! Minimum slice number, Maximum slice number 1 512 1 512 ! Clipping: Nxmin, Nxmax, Nymin, Nymax 4 4 2 ! Coarse graining: Nxc, Nyc, Nzc 1 ! Origin 0:Center 1:Lab 1 ! Z direction +1:forward or -1:backward 1 ! PHITS parameter: 0:Minimal 1:PhotonTherapy2:ParticleTherapy "PHITSinputs2/" ! Filename for PHITS input ⇒　phits.inp [ Transform ] $ Transform system according to DICOM header tr500 -32.50000 -32.50000 -29.75100 1.00000 0.00000 0.00000 0.00000 1.00000 0.00000 0.00000 0.00000 1.00000 1 [ Parameters ] … $ Optimized for Photon Therapy emin(2) = 1.000000000E-10 … $ Library data paths for file(7) & file(14) infl:{libpath.inp} Windowsデフォルトインストールの場合変更不要 file(7) = c:/phits/data/xsdir.jnd file(14) = c:/phits/data/trxcrd.dat Conversion of DICOM 13

14. 読込の高速化 目的 PHITSでは一度インプットファイルを全てバイナリー化してから再読込 巨大なボクセルデータをあらかじめバイナリー化して読込時間短縮！ 手順 ①　[Parameters]セクションのivoxelを有効にする（cを消す） ivoxel = 2 # 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 15

16. 講習の流れ • DICOM2PHITSの使い方 • 放射線治療シミュレーション • DICOMデータの変換 • ビームの設定 • タリーの設定 • 校正定数の設定 • 治療計画シミュレーション Beam setting 16

17. ビームの形成 照射野 コリメーター 8cm domパラメータで 円錐ビームを生成 9cm z y x Isocenter スペクトル(spectrum.inp) 100cm (x0,y0,z0)=(-0.5, 6.5, -9) Beam setting 17

18. ビームの形成 beam.inp domパラメータで 円錐ビームを生成 set: c91[ -0.5 ] $ isocenter x set: c92[ 6.5 ] $ isocenter y set: c93[ -9.0 ] $ isocenter z set: c94[ 100.0 ] $ beam distance set: c95[ 9.0 ] $ irradiation field x length set: c96[ 8.0 ] $ irradiation field y length [ Source ] s-type = 5 x0 = c91 x1 = c91 y0 = c92-c94 y1 = c92-c94 z0 = c93 z1 = c93 dir = 0.0 phi = 90 dom = asin(7.0/c94)*180.0/pi proj = photon e-type = 1 ne = 28 infl:{spectrum.inp} 5cm dom z y x ｙ軸正方向 スペクトルの設定 Beam setting 18

19. コリメーターの設定 beam.inp set: c91[ -0.5 ] $ isocenter x set: c92[ 6.5 ] $ isocenter y set: c93[ -9.0 ] $ isocenter z set: c94[ 100.0 ] $ beam distance set: c95[ 9.0 ] $ irradiation field x length set: c96[ 8.0 ] $ irradiation field y length 8cm [ Surface ] $ Outer region 99 so 500 $ Colimator inside 90 rpp c91-0.1*c95 c91+0.1*c95 c92-0.9*c94 c92-0.8*c94 c93-0.1*c96 c93+0.1*c96 $ Colimator outside 91 rpp c91-0.5*c95 c91+0.5*c95 c92-0.9*c94 c92-0.8*c94 c93-0.5*c96 c93+0.5*c96 [ Cell ] $ Main space 98 0 -99 #90 $ Void $ 90 -1 90 -91 $ Void colimator 90 10 -11.34 90 -91 $ Pb colimator 99 -1 99 $ Outer region 9cm 100cm Beam setting 19

20. ビームの確認 beam.inp [ Cell ] $ Main space 98 0 -99 #90 $ Void $ 90 -1 90 -91 $ Void colimator 90 10 -11.34 90 -91 $ Pb colimator 99 -1 99 $ Outer region [ Cell ] $ Main space 98 0 -99 #90 $ Void 90 -1 90 -91 $ Void colimator $90 10 -11.34 90 -91 $ Pb colimator 99 -1 99 $ Outer region マテリアルを外部void (-1)にすることで計算時間短縮 照射野@isocenter track-yz.eps track-xz.eps Beam setting 20

21. 多門照射の設定 前方照射86MU 後方照射155MU Isocenter: (x0,y0,z0)=(-0.5, 6.5, -9) z x-x0 y-y0 z-z0 x0 y0 z0 x’ y’ z’ cosq sinq 0 -sinq cosq 0 0 0 1 y Isocenterを中心にZ軸回りにq回転 + = x x y z x0(1-cosq)+y0sinq y0(1-cosq)-x0sinq 0 x’ y’ z’ cosq sinq 0 -sinq cosq 0 0 0 1 + = Beam setting 21

22. 多門照射の設定 beam2.inp ソースのtransformを指定 [ Transform ] set: c91[ -0.5 ] $ isocenter x set: c92[ 6.5 ] $ isocenter y set: c93[ -9.0 ] $ isocenter z set: c94[ 100.0 ] $ beam distance set: c95[ 9.0 ] $ irradiation field x length set: c96[ 8.0 ] $ irradiation field z length set: c97[ 0.0 ] $ gantory angle set: c98[cos(c97/180*pi) ] $ cos set: c99[sin(c97/180*pi) ] $ sin tr600 c91*(1-c98)+c92*c99 -c91*c99+c92*(1-c98) 0 c98 -c99 0 c99 c98 0 0 0 1 1 [ Source ] trcl = 600 s-type = 5 … ① ② ③ ⑤ ⑥ ④ ⑦ ⑧ ⑨ ⑩ ⑪ ⑫ ⑥ ⑤ ④ x y z ① x0(1-cosq)+y0sinq y0(1-cosq)-x0sinq 0 x’ y’ z’ cosq sinq 0 -sinq cosq 0 0 0 1 ② ⑦ ⑨ ⑧ + = ⑪ ⑫ ⑩ ③ Beam setting 22

23. 多門照射の確認 track-yz-1.eps track-yz-2.eps beam2.inp [ Transform ] … set: c97[ 0.0 ] $ gantory angle … 180.0 [ T-track ] … file = track-xz-1.out … track-xz-2.out [ T-track ] … file = track-yz-1.out … track-yz-2.out track-xz-1.eps track-xz-2.eps Beam setting 23

24. 照射結果の合算 sumtally.inp "track-yz-3.out" 1.0 2 1.0 "track-yz-1.out” 1.0 "track-yz-2.out” track-yz-3.eps Beam setting 24

25. 講習の流れ • DICOM2PHITSの使い方 • 放射線治療シミュレーション • DICOMデータの変換 • ビームの設定 • タリーの設定 • 校正定数の設定 • 治療計画シミュレーション Tally setting 25

26. タリーの設定 phits1.inp Isocenter面 $ parameters set: c81[ 128] $ number of x pixel set: c82[ 128] $ number of y pixel set: c83[ 74] $ number of z pixel set: c84[ 0.50781] $ unit voxel x set: c85[ 0.50781] $ unit voxel y set: c86[ 0.50000] $ unit voxel z set: c87[ -0.06348] $ smallest x set: c88[ -0.06348] $ smallest y set: c89[ -0.12500] $ smallest z set: c90[ 0.00001] $ small quota [ T-Deposit ] title = [t-deposit: icntl=0] file = deposit-xy-1.out mesh = xyz x-type = 2 nx = 128 xmin = c87 -32.50000 xmax = c87+c81*c84-32.50000 y-type = 2 ny = 128 ymin = c88 -32.50000 ymax = c88+c82*c85-32.50000 z-type = 2 nz = 1 zmin = c93 -1.0 zmax = c93 +1.0 unit = 0 axis = xy output = dose epsout = 1 set: c91[ -0.5 ] $ isocenter x set: c92[ 6.5 ] $ isocenter y set: c93[ -9.0 ] $ isocenter z [ Transform ] $ Transform system according to DICOM header tr500 -32.50000 -32.50000 -29.75100 1.00000 0.00000 0.00000 0.00000 1.00000 0.00000 0.00000 0.00000 1.00000 1 Tally setting 26

27. 試行結果 deposit-xz-1.eps deposit-xy-1.eps 統計量を増やす 絶対線量に直す deposit-yz-1.eps 27 Tally setting

28. 講習の流れ • DICOM2PHITSの使い方 • 放射線治療シミュレーション • DICOMデータの変換 • ビームの設定 • タリーの設定 • 校正定数の設定 • 治療計画シミュレーション Conversion factor 28

29. 校正定数の設定 水ファントムへの照射 water.inp 100cm 10cm X=2.6635*10-14 (Gy/source) 10cm 6MVの光子でのPDD(10cm)=66.1% sumtally.inp CF=0.661*10-2 / X (source/MU) "deposit-xy-3.out" 59.81e12 2 86.0 "deposit-xy-1.out” 155.0 "deposit-xy-2.out” Dabs=DMC*CF*MU MU(241)*CF(0.661e-2/2.6635e-14) Conversion factor 29

30. 講習の流れ • DICOM2PHITSの使い方 • 放射線治療シミュレーション • DICOMデータの変換 • ビームの設定 • タリーの設定 • 校正定数の設定 • 治療計画シミュレーション Results 30

31. 治療計画シミュレーション "deposit-xy-3x.out" 59.81e12 2 86.0 "deposit-xy-1x.out” 155.0 "deposit-xy-2x.out” "deposit-xy-3.out" 59.81e12 2 86.0 "deposit-xy-1.out” 155.0 "deposit-xy-2.out” 統計量増(phits2.inp) ⇒ 実行結果LongCalcディレクトリ deposit-xy-1.out　（前方照射の結果） deposit-xy-2.out　（後方照射の結果） maxcas = 10000000 maxbch = 10 sumtally.inp ANGELを使用してEPS化 deposit-xy-3.eps deposit-xy-3x.eps deposit-xy-3_err.eps Results 31

32. おまけ：ANGELの上手な使い方 deposit-xy-1.out きれいな出力のために変更済み カラーマップの上限 下限 104 #MODIFIED The line below is setting minimum and maxmum of the color 105 set: c3[ 1.400000E-00] c4[ 2.100000E-00] 線形表示 114 #MODIFIED The line below is parameters for contours 115 p: ZLIN IPD2 ICUT(8) CUTS(2.1,2.0,1.9,1.8,1.7,1.6,1.5,1.4) COLS … 118 #MODIFIED The line below is parameters for frames 119 p: NOLG NOFR NOMS NOTL ANGL(-90) グラフのフレームに関するパラメータ 等高線用パラメータ 127 #MODIFIED Color map (hc) is changed to countor (h2) 128 # 129 130 h2: y = 32.18229 to -32.30957 by 0.5078100 ;… … hc: カラーマップ　h2: 等高線 1785 #MODIFIED Unit has been changed from Gy/source to Gy 1786 y: Dose [Gy] 1787 #MODIFIED Liner color scale 1788 p: YLIN カラーバーの横の表示 カラーバーの線形表示 詳しくはANGELのマニュアル(phits/manual/manual-angel4.pdf)参照 Results 32

33. おまけ：ANGELの上手な使い方 deposit-xy-1x.out Icntl=8 で出力される幾何形状表示をファイルに追加 123 #MODIFIED CT-data color map added from result with icntl=8 124 # 125 126 #-------------------------------------------------------------- 127 # gshow 128 #-------------------------------------------------------------- 129 130 p: legs[c5*0.875] 131 132 h: x ny1,0 y2=[y1](void),ic[lightgray] 133 -1.000000E+00 -1.000000E+00 134 135 136 hb: line clip z cb(lightgray) 137 frame: 138 -3.2563480E+01 -3.2563480E+01 139 3.2436200E+01 -3.2563480E+01 … ⇒　幾何形状表示（CTデータ）と線量分布を同時に表示 Results 33