Multi-Purpose Particle and Heavy Ion Transport code System

Presentation on theme: "Multi-Purpose Particle and Heavy Ion Transport code System"— Presentation transcript:

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

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

3 講習の流れ DICOM2PHITSの使い方 放射線治療シミュレーション DICOMデータの変換 ビームの設定 タリーの設定 校正定数の設定
治療計画シミュレーション PHITS講習会　入門実習 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） Dicom形式のデータをPHITS形式のボクセルデータに変換
CT値と物質密度・組成の関係はdata/HumanVoxelTable.data　[W. Schneider, Phys. Med. Biol. 45(2000) ]を参照 ①　基本情報ファイルを作成（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 ! Minimum slice number, Maximum slice number ! Clipping: Nxmin, Nxmax, Nymin, Nymax ! Coarse graining: Nxc, Nyc, Nzc ! Origin 0:Center 1:Lab ! Z direction +1:forward or -1:backward ! PHITS parameter: 0:Minimal 1:PhotonTherapy 2:ParticleTherapy "PHITSinputs/" ! Filename for PHITS input DICOMデータの一部を切り出してボクセル化できる 画像を粗くする（分解能を下げる）ことも可能 原点情報をDICOMヘッダーから抽出できる PHITSオプションの最適化を選択 PHITSインプットを格納するディレクトリを指定 DICOM2PHITS HowTo 5

6 DICOM2PHITS実行結果 DICOM2PHITS HowTo 6
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 Origin=1でDICOMヘッダーから原点情報取得し、その位置に平行移動
PHITSインプットファイル file=phits.inp [Parameters ] … phits.inp [ 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 必ず１行目！ 　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[ ] $ 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 ボクセルファントムのピクセル数や大きさはDICOMのヘッダーから自動的に設定 ボクセルファントムの中心がxyz座標の原点にくるように配置されている Include fileコマンド ファントムを回転・ 平行移動が可能 Origin=1でDICOMヘッダーから原点情報取得し、その位置に平行移動 PHITS形式詳細は phits/lecture/voxel/phits-lec-voxel-jp.ppt 「Voxelファントムの作り方」を参照 DICOM2PHITS HowTo 7

8 CT値⇔物質密度・組成変換表 DICOM2PHITS HowTo 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値よりも大きい　⇒　ワーニングを出して最後の物質で代用 DICOM2PHITS HowTo 21

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

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

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

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

13 Windowsデフォルトインストールの場合変更不要
DICOMデータの変換 ①　dicom2phits.inpを以下のように変更して実行 ②　libpath.inpを環境に合わせて変更 "data/HumanVoxelTable-KumamotoUniv.data" ! File for conversion of human voxel data "DICOM-KumamotoUniv/" ! DICOM file directory ! Minimum slice number, Maximum slice number ! Clipping: Nxmin, Nxmax, Nymin, Nymax ! Coarse graining: Nxc, Nyc, Nzc ! Origin 0:Center 1:Lab ! Z direction +1:forward or -1:backward ! PHITS parameter: 0:Minimal 1:PhotonTherapy 2:ParticleTherapy "PHITSinputs2/" ! Filename for PHITS input ⇒　phits.inp [ Transform ] $ Transform system according to DICOM header tr 1 [ Parameters ] … $ Optimized for Photon Therapy emin(2) = E-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 = # 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データの変換 ビームの設定 タリーの設定 校正定数の設定
治療計画シミュレーション PHITS講習会　入門実習 Beam setting 16

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

18 ビームの形成 beam.inp x y z Beam setting 18 domパラメータで 円錐ビームを生成 5cm dom ｙ軸正方向
set: c91[ ] $ isocenter x set: c92[ ] $ isocenter y set: c93[ ] $ isocenter z set: c94[ ] $ beam distance set: c95[ ] $ irradiation field x length set: c96[ ] $ 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} x y z dom 5cm ｙ軸正方向 スペクトルの設定 Beam setting 18

19 コリメーターの設定 beam.inp Beam setting 19 8cm 9cm 100cm
set: c91[ ] $ isocenter x set: c92[ ] $ isocenter y set: c93[ ] $ isocenter z set: c94[ ] $ beam distance set: c95[ ] $ irradiation field x length set: c96[ ] $ irradiation field y length [ 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 # $ Void $ $ Void colimator $ Pb colimator $ Outer region Beam setting 19

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

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

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

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

24 照射結果の合算 sumtally.inp track-yz-3.eps Beam setting 24 "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データの変換 ビームの設定 タリーの設定 校正定数の設定
治療計画シミュレーション PHITS講習会　入門実習 Tally setting 25

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

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

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

29 校正定数の設定 水ファントムへの照射 water.inp X=2.6635*10-14 (Gy/source)
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 "deposit-xy-1.out” "deposit-xy-2.out” Dabs=DMC*CF*MU MU(241)*CF(0.661e-2/2.6635e-14) Conversion factor 29

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

31 治療計画シミュレーション sumtally.inp ANGELを使用してEPS化 deposit-xy-3.eps
"deposit-xy-3x.out" 59.81e12 2 86.0 "deposit-xy-1x.out” 155.0 "deposit-xy-2x.out” 統計量増(phits2.inp) ⇒ 実行結果LongCalcディレクトリ "deposit-xy-3.out" 59.81e12 2 86.0 "deposit-xy-1.out” 155.0 "deposit-xy-2.out” maxcas = maxbch = deposit-xy-1.out　（前方照射の結果） deposit-xy-2.out　（後方照射の結果） 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[ E-00] c4[ E-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 = to by ; … … 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 Results 33
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] E E+00 134 135 136 hb: line clip z cb(lightgray) 137 frame: E E+01 E E+01 … ⇒　幾何形状表示（CTデータ）と線量分布を同時に表示 Results 33