Transcript 放射線治療計画シミュレーション

PHITS
Multi-Purpose Particle and Heavy Ion Transport code System
PHITS 講習会
放射線治療計画シミュレーション
2014年6月改訂
title
1
目的
治療計画シミュレーションの検証をPHITSを用いて行う
後方照射155MU
前方照射86MU
Introduction
2
講習の流れ
1. DICOM2PHITSの使い方
2. 放射線治療シミュレーション
a.
b.
c.
d.
e.
DICOMデータの変換
ビームの設定
タリーの設定
校正定数の設定
治療計画シミュレーション
Table of Contents
3
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）
指定ディレクトリ内に含まれるDICOMファイルを自動判別
指定した番号範囲のスライスファイルを読み込む
"data/HumanVoxelTable.data"
"DICOM/"
1 20
70 430 90 460
4 4 1
0
1
0
"PHITSinputs/"
!
!
!
!
!
!
!
!
!
File for conversion of human voxel data
DICOM DICOMデータの一部を切り出してボクセル化できる
file directory
Minimum画像を粗くする（分解能を下げる）ことも可能
slice number, Maximum slice number
Clipping: Nxmin, Nxmax, Nymin, Nymax
Coarse原点情報をDICOMヘッダーから抽出できる
graining: Nxc, Nyc, Nzc
Origin 0:Center 1:Lab
Z direction +1:forward or -1:backward
PHITSオプションの最適化を選択
PHITS parameter: 0:Minimal 1:PhotonTherapy 2:ParticleTherapy
Filename for PHITS input
PHITSインプットを格納するディレクトリを指定
② (Windows) dicom2phits_win.batにdicom2phits.inpをドラッグ&ドロップ
(Mac) dicom2phits_mac.commnadをダブルクリックし、現れる窓にdicom2phits.inpと入力
DICOM2PHITS HowTo
5
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
PHITSインプットファイル
file=phits.inp
[Parameters ]
…
phits.inp
必ず１行目！
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のヘッ
ダーから自動
的に設定
[Cell]
$ 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
ボクセルファントム Include fileコマンド
の中心がxyz座標
ファントムを回転・
の原点にくるように
平行移動が可能
配置されている
Origin=1でDICOMヘッダーから原点
情報取得し、その位置に平行移動
PHITS形式詳細は phits/lecture/voxel/phits-lec-voxel-jp.ppt 「Voxelファントムの作り方」を参照
DICOM2PHITS HowTo
7
CT値⇔物質密度・組成変換表
data/HumanVoxelTable.data サンプル表 [W. Schneider, Phys. Med. Biol. 45(2000)459を参照]
2行目：分割数
1行目：実行時に画面に表示するコメント
Table for human voxel based on W. Schneider, Phys.…
24
! Number of universe di…
3行目：物質１の定義
-1000.0d0 -1.21e-3 3
! Lowest
CT value, Dens
-950.0d0 -0.26
10 ! Universe 2
最小CT値
-120.0d0 -0.927
8
! Universe 3
-1000  物質１ < -950 -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/…
最後の物質
…
の最大CT値
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
組成割合： >0 = 粒子密度比、< 0 ＝質量比
14N -3.1
…
物質１の最小CT値よりも小さい ⇒ ワーニングを出して物質１で代用
最後の物質の最大CT値よりも大きい ⇒ ワーニングを出して最後の物質で代用
DICOM2PHITS HowTo
21
ジオメトリの確認
icntl = 11
CT3D.eps
icntl = 8
deposit-xy.eps
DICOM2PHITS HowTo
9
ボクセル化する範囲の変更
dicom2phits.inp
"data/HumanVoxelTable.data"
"DICOM/"
1 10
70 270 90 270
4 4 1
0
1
0
"PHITSinput1/"
! File for conversion of human voxel data
! 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 …
! Filename for voxel include file
変更後(phits.inp)
CT3D.eps
DICOM2PHITS HowTo
deposit-xy.eps
10
線量計算結果
icntl = 0
unit = 2
⇒ MeV/source 単位
unit = 0
⇒ Gy/source 単位
deposit-xy.eps
xyzメッシュを用いた[t-deposit]の結果
DICOM2PHITS HowTo
11
講習の流れ
1. DICOM2PHITSの使い方
2. 放射線治療シミュレーション
a.
b.
c.
d.
e.
DICOMデータの変換
ビームの設定
タリーの設定
校正定数の設定
治療計画シミュレーション
Conversion of DICOM
12
DICOMデータの変換
① dicom2phits.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}
② libpath.inpを環境に合わせて変更
file(7) = c:/phits/data/xsdir.jnd
file(14) = c:/phits/data/trxcrd.dat
Windowsデフォルトインストール
の場合変更不要
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
講習の流れ
1. DICOM2PHITSの使い方
2. 放射線治療シミュレーション
a.
b.
c.
d.
e.
DICOMデータの変換
ビームの設定
タリーの設定
校正定数の設定
治療計画シミュレーション
Beam setting
16
ビームの形成
コリメーター
domパラメータで
円錐ビームを生成
8cm
9cm
z
x
照射野
y
スペクトル(spectrum.inp)
Isocenter
100cm
Beam setting
(x0,y0,z0)=(-0.5, 6.5, -9)
17
ビームの形成
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
[ Source ]
s-type = 5
x0 = c91
x1 = c91
y0 = c92-c94
y1 = c92-c94
ｙ軸正方向
z0 = c93
z1 = c93
dir = 0.0
phi = 90
dom = atan(7.0/c94)*180.0/pi
proj = photon
e-type = 1
ne = 28
infl:{spectrum.inp}
Beam setting
domパラメータで
円錐ビームを生成
z
x
7cm
y
dom
スペクトルの設定
18
コリメーターの設定
beam.inp
c91[
c92[
c93[
c94[
c95[
c96[
-0.5
6.5
-9.0
100.0
9.0
8.0
]
]
]
]
]
]
$
$
$
$
$
$
isocenter x
isocenter y
isocenter z
beam distance
irradiation field x length
irradiation field y length
9cm
set:
set:
set:
set:
set:
set:
8cm
100cm
[ 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
$ 90 -1
90 -91
90 10 -11.34 90 -91
99 -1
99
$ Void
$ Void colimator
$ Pb colimator
$ Outer region
Beam setting
19
ビームの確認
beam.inp
[ Cell ]
$ Main space
98 0
-99 #90
$ 90
90-1-1
90-91
-91
90
$ 90901010-11.34
-11.34 90
90-91
-91
99 -1
99
$ Void
Voidcolimator
colimator
$$Void
$$Pb
Pbcolimator
colimator
$ Outer region
track-yz.eps
マテリアルを外部void (-1)
にすることで計算時間短縮
照射野@isocenter
track-xz.eps
Beam setting
20
多門照射の設定
前方照射86MU
後方照射155MU
z
x
y
Isocenter: (x0,y0,z0)=(-0.5, 6.5, -9)
Isocenterを中心に
Z軸回りにq回転
cosq -sinq 0
x’
y’ = sinq cosq 0
z’
0
0 1
cosq -sinq 0
x’
y’ = sinq cosq 0
z’
0
0 1
x-x0
x0
y-y0 + y0
z-z0
z0
x0(1-cosq)+y0sinq
x
y + y0(1-cosq)-x0sinq
z
0
Beam setting
21
多門照射の設定
beam2.inp
[ Transform ]
set: c91[
-0.5
] $
set: c92[
6.5
] $
set: c93[
-9.0
] $
set: c94[
100.0 ] $
set: c95[
9.0
] $
set: c96[
8.0
] $
set: c97[
0.0
] $
set: c98[cos(c97/180*pi)
set: c99[sin(c97/180*pi)
tr600
ソースのtransformを指定
isocenter x
isocenter y
isocenter z
beam distance
irradiation field x length
irradiation field z length
gantory angle
] $ cos
] $ sin
①
[ Source ]
trcl = 600
s-type = 5
…
② ③
c91*(1-c98)+c92*c99 -c91*c99+c92*(1-c98) 0
c98 ④-c99⑤ 0 ⑥
c99 ⑦ c98 ⑧ 0 ⑨
0 ⑩ 0 ⑪ 1⑫
1
x’
y’
=
z’
④
⑤ ⑥
cosq -sinq 0
sinq ⑦ cosq⑧ 0⑨
⑪
⑫
⑩
0
0 1
x
y
z
+
x0(1-cosq)+y0sinq ①
y0(1-cosq)-x0sinq ②
0③
Beam setting
22
多門照射の確認
track-yz-1.eps
track-yz-2.eps
beam2.inp
[ Transform ]
…
set: c97[ 180.0
0.0
…
] $ gantory angle
[ T-track ]
…
file
= track-xz-2.out
track-xz-1.out
…
[ T-track ]
…
file
= track-yz-2.out
track-yz-1.out
…
track-xz-1.eps
track-xz-2.eps
Beam setting
23
照射結果の合算
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
講習の流れ
1. DICOM2PHITSの使い方
2. 放射線治療シミュレーション
a.
b.
c.
d.
e.
DICOMデータの変換
ビームの設定
タリーの設定
校正定数の設定
治療計画シミュレーション
Tally setting
25
タリーの設定
phits1.inp
$ 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
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
Isocenter面
[ 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
26
deposit-xz-1.eps
試行結果
deposit-xy-1.eps
統計量を増やす
絶対線量に直す
deposit-yz-1.eps
Tally setting
27
講習の流れ
1. DICOM2PHITSの使い方
2. 放射線治療シミュレーション
a.
b.
c.
d.
e.
DICOMデータの変換
ビームの設定
タリーの設定
校正定数の設定
治療計画シミュレーション
Conversion factor
28
校正定数の設定
water.inp
水ファントムへの照射
100cm
10cm
X=2.6635*10-14 (Gy/source)
10cm
6MVの光子でのPDD(10cm)=66.1%
CF=0.661*10-2 / X (source/MU)
Dabs=DMC*CF*MU
MU(241)*CF(0.661e-2/2.6635e-14)
Conversion factor
sumtally.inp
"deposit-xy-3.out"
59.81e12
2
86.0
"deposit-xy-1.out”
155.0
"deposit-xy-2.out”
29
講習の流れ
1. DICOM2PHITSの使い方
2. 放射線治療シミュレーション
a.
b.
c.
d.
e.
DICOMデータの変換
ビームの設定
タリーの設定
校正定数の設定
治療計画シミュレーション
Results
30
治療計画シミュレーション
統計量増(phits2.inp) ⇒ 実行結果LongCalcディレクトリ
maxcas
maxbch
=
=
10000000
10
deposit-xy-1.out （前方照射の結果）
deposit-xy-2.out （後方照射の結果）
ANGELを使用してEPS化
deposit-xy-3.eps
deposit-xy-3x.eps
"deposit-xy-3x.out"
"deposit-xy-3.out"
59.81e12
2
"deposit-xy-1x.out”
86.0 "deposit-xy-1.out”
"deposit-xy-2x.out”
155.0 "deposit-xy-2.out”
sumtally.inp
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
115
…
118
119
#MODIFIED The line below is parameters for contours
p: ZLIN IPD2 ICUT(8) CUTS(2.1,2.0,1.9,1.8,1.7,1.6,1.5,1.4) COLS
#MODIFIED The line below is parameters for frames
p: NOLG NOFR NOMS NOTL ANGL(-90)
127 #MODIFIED Color map (hc) is changed to countor (h2)
128 #
hc: カラーマップ h2: 等高線
129
130 h2: y = 32.18229
to -32.30957
by 0.5078100
…
1785
1786
1787
1788
等高線用
パラメータ
; …
#MODIFIED Unit has been changed from Gy/source to Gy
y: Dose [Gy]
カラーバーの横の表示
#MODIFIED Liner color scale
p: YLIN
カラーバーの線形表示
詳しくはANGELのマニュアル(phits/manual/manual-angel4.pdf)参照
Results
32
おまけ：ANGELの上手な使い方
deposit-xy-1x.out
Icntl=8 で出力される幾何形状表示をファイルに追加
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
…
#MODIFIED CT-data color map added from result with icntl=8
#
#-------------------------------------------------------------# gshow
#-------------------------------------------------------------p: legs[c5*0.875]
h: x ny1,0 y2=[y1](void),ic[lightgray]
-1.000000E+00 -1.000000E+00
hb: line clip z cb(lightgray)
frame:
-3.2563480E+01 -3.2563480E+01
3.2436200E+01 -3.2563480E+01
⇒ 幾何形状表示（CTデータ）と線量分布を同時に表示
Results
33