Multi-Purpose Particle and Heavy Ion Transport code System PHITS Multi-Purpose Particle and Heavy Ion Transport code System 放射線治療計画シミュレーション PHITS講習会 入門実習 2014年5月改訂 title 1
目的 治療計画シミュレーションの検証をPHITSを用いて行う 後方照射155MU 前方照射86MU Introduction 2
講習の流れ DICOM2PHITSの使い方 放射線治療シミュレーション DICOMデータの変換 ビームの設定 タリーの設定 校正定数の設定 治療計画シミュレーション PHITS講習会 入門実習 Table of Contents 3
１つのスライスに対するデータ（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
変換プログラム(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:PhotonTherapy 2:ParticleTherapy "PHITSinputs/" ! Filename for PHITS input DICOMデータの一部を切り出してボクセル化できる 画像を粗くする（分解能を下げる）ことも可能 原点情報をDICOMヘッダーから抽出できる PHITSオプションの最適化を選択 PHITSインプットを格納するディレクトリを指定 DICOM2PHITS HowTo 5
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
Origin=1でDICOMヘッダーから原点情報取得し、その位置に平行移動 PHITSインプットファイル file=phits.inp [Parameters ] … phits.inp [ 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
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.… 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
ジオメトリの確認 icntl = 11 icntl = 8 CT3D.eps deposit-xy.eps DICOM2PHITS HowTo 9
ボクセル化する範囲の変更 変更後(phits.inp) DICOM2PHITS HowTo 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) CT3D.eps deposit-xy.eps DICOM2PHITS HowTo 10
xyzメッシュを用いた[t-deposit]の結果 線量計算結果 icntl = 0 unit = 2 ⇒ MeV/source 単位 unit = 0 ⇒ Gy/source 単位 deposit-xy.eps xyzメッシュを用いた[t-deposit]の結果 DICOM2PHITS HowTo 11
講習の流れ DICOM2PHITSの使い方 放射線治療シミュレーション DICOMデータの変換 ビームの設定 タリーの設定 校正定数の設定 治療計画シミュレーション PHITS講習会 入門実習 Conversion of DICOM 12
Windowsデフォルトインストールの場合変更不要 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:PhotonTherapy 2: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
巨大なボクセルデータをあらかじめバイナリー化して読込時間短縮！ 読込の高速化 目的 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
ジオメトリの確認 icntl = 8 deposit-xz.eps deposit-xy.eps Conversion of DICOM 15
講習の流れ DICOM2PHITSの使い方 放射線治療シミュレーション DICOMデータの変換 ビームの設定 タリーの設定 校正定数の設定 治療計画シミュレーション PHITS講習会 入門実習 Beam setting 16
ビームの形成 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
ビームの形成 beam.inp x y z Beam setting 18 domパラメータで 円錐ビームを生成 5cm 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} x y z dom 5cm ｙ軸正方向 スペクトルの設定 Beam setting 18
コリメーターの設定 beam.inp Beam setting 19 8cm 9cm 100cm 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 [ 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 Beam setting 19
ビームの確認 beam.inp 照射野@isocenter track-yz.eps track-xz.eps Beam setting [ 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
多門照射の設定 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 0 0 1 x-x0 y-y0 z-z0 x0 y0 z0 = + Isocenterを中心にZ軸回りにq回転 x’ y’ z’ cosq sinq 0 -sinq cosq 0 0 0 1 x y z x0(1-cosq)+y0sinq y0(1-cosq)-x0sinq = + Beam setting 21
多門照射の設定 beam2.inp Beam setting 22 ソースのtransformを指定 x’ y’ z’ 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’ cosq sinq 0 -sinq cosq 0 0 0 1 x y z x0(1-cosq)+y0sinq y0(1-cosq)-x0sinq + = Beam setting 22
多門照射の確認 track-xz-2.eps track-yz-2.eps track-yz-1.eps beam2.inp [ Transform ] … set: c97[ 0.0 ] $ 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
照射結果の合算 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
講習の流れ DICOM2PHITSの使い方 放射線治療シミュレーション DICOMデータの変換 ビームの設定 タリーの設定 校正定数の設定 治療計画シミュレーション PHITS講習会 入門実習 Tally setting 25
タリーの設定 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[ 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
試行結果 deposit-xz-1.eps deposit-xy-1.eps 統計量を増やす 絶対線量に直す deposit-yz-1.eps Tally setting 27
講習の流れ DICOM2PHITSの使い方 放射線治療シミュレーション DICOMデータの変換 ビームの設定 タリーの設定 校正定数の設定 治療計画シミュレーション PHITS講習会 入門実習 Conversion factor 28
校正定数の設定 水ファントムへの照射 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 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
講習の流れ DICOM2PHITSの使い方 放射線治療シミュレーション DICOMデータの変換 ビームの設定 タリーの設定 校正定数の設定 治療計画シミュレーション PHITS講習会 入門実習 Results 30
治療計画シミュレーション 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 = 10000000 maxbch = 10 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
おまけ：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
おまけ：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] 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