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 Voxelファントムの作り方 PHITS講習会　入門実習 2014年8月改訂 title 1

2 Voxelファントムとは？ 直方体を積み重ねて生物など複雑な体系を再現したもの （マニュアル4.7.5.3）
PHITSのUniverseとLattice構造を使って構築 （マニュアル4.7.3＆4.7.4） 低分解能 高分解能 Introduction 2

3 Voxelファントムを用いた計算例 粒子線治療時の２次放射線被ばく線量評価 BNCTの治療効果推定 Introduction 3
T. Sato et al. Radiat. Res. (2009) BNCTの治療効果推定 H. Kumada et al. J. Phys.: Conf. Ser. (2007) Introduction 3

4 講習の流れ Universe Lattice 簡易Voxelファントム Dicom形式からの変換 まとめ Table of Contents
PHITS講習会　入門実習 Table of Contents 4

5 Universeとは？ → PHITSの中の仮想空間 PHITSでは，いくつもの仮想空間（universe）を定義できる Universe1
実際に粒子輸送の舞台となるのはメイン空間のみ メイン空間の一部を別のuniverseで満たすことにより入れ子構造を再現する メイン空間 メイン空間にある箱の中を，球が並んでいる別の宇宙(universe1)に置換 Universe 5

6 Universeの例 universe.inp Universe 6 満たす領域は常にvoid Universe1で満たす （Fillする）
[ C e l l ] $ Main space FILL=1 FILL=2 9 -1 #1 #2 $ Universe 1 U=1 102 0 #101 U=1 $ Universe 2 U=2 #201 U=2 [ S u r f a c e ] 10 CY 5 11 PX -6 12 PX 6 13 PY -6 14 PY 6 15 PZ -6 16 PZ 6 17 PZ 0 Universe1で満たす （Fillする） Universe1を宣言 PX -3 PX 9 Universe 6

7 講習の流れ Universe Lattice 簡易Voxelファントム Dicom形式からの変換 まとめ Table of Contents
PHITS講習会　入門実習 Table of Contents 7

8 Latticeとは？ → PHITSの中で使う 格子(lattice，繰り返し）構造
同じ形状のものが繰り返し並んでいるときに，個々の形状を全て定義するのは大変！ 基本構造のみ定義して，あとはその繰り返しで表現する 格子には直方体と六角柱があるが，基本的に使うのは直方体のみ Latticeの例 Lattice 8

9 Universeを２つ以上使った２重入れ子構造にする
格子の定義方法 格子空間には他の構造物は定義できない 格子の中身は直接定義できない メイン空間ではなく，universeの１つとして定義するのが便利 別のuniverseで満たす(fillする）必要がある Universeを２つ以上使った２重入れ子構造にする fill fill メイン空間 Universe1 Universe2 Lattice 9

10 格子の定義 lattice.inp 領域101 Lattice 10 （2,2,0） 基本格子 （0,0,0） （-2,-2,0）
[ S u r f a c e ] 1 rpp 2 rpp 99 so 100 101 rpp 201 sph [ C e l l ] $ Main space fill=1 $ Universe 1 lat=1 u=1 fill=-2:2 -2:2 0:0 $ Universe 2 u=2 u=2 （2,2,0） -5　 Y軸 基本格子 （0,0,0） （-2,-2,0） 　 X軸 　 5 領域101 直方体タイプのLatticeを宣言 基本格子の領域を定義 格子に物質を満たす領域を指定（0は基本格子） 各格子に満たすuniverse番号（5×5×1の行列） 基本格子と座標を合わせる必要有 メイン空間で使う領域は，必ず定義しておく必要有 Lattice 10

11 格子の中身の変更 lattice.inp １つ目の格子のみVoidに変更 変更前 (lattice.inp) 変更後
[ C e l l ] $ Main space fill=1 $ Universe 1 lat=1 u=1 fill=-2:2 -2:2 0:0 $ Universe 2 u=2 u=2 $ Universe 3 u=3 変更前 (lattice.inp) 変更後 (lattice1.inp) Lattice 11

12 講習の流れ Universe Lattice 簡易Voxelファントム Dicom形式からの変換 まとめ Table of Contents
PHITS講習会　入門実習 Table of Contents 12

13 簡易Voxelファントムの作り方 ① 単一の物質（骨・軟組織など）で満たされたuniverseを作る
(void) Universe2 (water) Universe3 (Aluminum) ②　①で作ったuniverseをLatticeを使って組み合わせ，ボクセルファントムが入ったuniverseを作る ③　②で作ったuniverseを メイン空間の一部にFillする Universe10 (Voxel Phantom) メイン空間 Simple Voxel Phantom 13

14 インプットファイル robot.inp x y z Simple Voxel Phantom 14
基本格子（この場合は１つ目に定義する領域）の面 [ C e l l ] $ Material universe u=1 u=2 u=3 $ Voxel universe lat=1 u=10 fill=0:4 0:4 0:4 … 以下，４回繰り返し $ Main space fill=10 [ S u r f a c e ] $ fundamental voxel 101 rpp 99 so 100 $ Main space 201 rpp 202 rcc 203 rcc 十分に大きい領域であればどこでもよい X正方向，Y正方向，Z正方向 （図の左下から始まる） x y z Simple Voxel Phantom 14

15 物質の追加・変更方法 robot.inp ファントムの頭を水から銅に変更する 変更前 (robot.inp) 変更後
[ C e l l ] $ Material universe u=1 u=2 u=3 u=4 $ Voxel universe lat=1 u=10 fill=0:4 0:4 0:4 … （始めの4ブロック） $ Main space fill=10 ファントムの頭を水から銅に変更する 変更前 (robot.inp) 変更後 (robot1.inp) Simple Voxel Phantom 15

16 線量計算結果 robot-heat-xz.eps robot-heat-reg.out regメッシュを用いた [t-heat]計算
x: Serial Num. of Region y: Heat [MeV/source] h: x n n y(total),l3 n # num reg volume heat r.err E E E E E E regメッシュを用いた [t-heat]計算 各領域別（手足・胴・頭）の 線量を計算 xyzメッシュを用いた [t-heat]計算 腫瘍部分や重要臓器などに ROIをかけた線量計算が可能 ファントム内の線量分布を可視化 Simple Voxel Phantom 16

17 講習の流れ Universe Lattice 簡易Voxelファントム Dicom形式からの変換 まとめ Table of Contents
PHITS講習会　入門実習 Table of Contents 17

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

19 DICOM2PHITSの使い方の詳細は「DICOM2PHITSの使い方」を参照
phits/utility/dicom2phits/phits-lec-dicom2phits-jp.ppt ①入力ファイルを作成（dicom2phits.inp） "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と入力 ⇒　PHITSinputs/ディレクトリにサンプルインプット生成 DICOM2PHITS HowTo 19

20 巨大なボクセルデータをあらかじめバイナリー化して読込時間短縮！
読込の高速化 目的 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} Dicom to PHITS 20

21 講習の流れ Universe Lattice 簡易Voxelファントム Dicom形式からの変換 まとめ Table of Contents
PHITS講習会　入門実習 Table of Contents 21

22 まとめ Voxelファントムは，UniverseとLatticeを組み合 わせて作成
Dicom形式は，別プログラム(DICOM2PHITS）を 用いてＰＨＩＴＳ入力形式（テキスト形式）に変換 ivoxelパラメータを設定することにより読込を高 速化 Summary 22

23 高分解能ファントム利用時の注意点 対策方法 メモリが（物理的に）足りない場合 23
全身を高分解能でボクセル化すると使用メモリが膨大になり実行できなくなる 例）全身（180cm×30cm×50cm）を1mm3でボクセル化すると，2億7000万個の 　　ボクセルが必要→　1CPUあたり最低5.4GBのメモリが必要（約20Byte/ボクセル） 初期設定のPHITSは，約2GBしかメモリを使うことができない 対策方法 “src”フォルダにある”param.inc”を変更し，全てのファイルを再コンパイル*する mdasの変更： PHITSで使用する総メモリ数（Byte)/8　（32bit機は特に上限値に注意） latmaxの変更： 最も大きなlattice数+1（複数Latticeを使う場合は，その最大値） 必要に応じて型宣言及びコンパイラオプションを追加： 詳細は次ページ参照 メモリが（物理的に）足りない場合 ボクセル化する領域を頭・胴・足などに分割し，無駄な領域を小さくする 複数のピクセルをまとめてボクセル化し，分解能を下げる *Windowsの場合，Stack memoryの関係からIntel Fortranではエラーが起きる場合があるので，gfortranの方がよい。Linux&MacはIntel Fortranでも動作します 23

24 補足資料① Voxel数が約5000千万個を超える場合は，様々な対策が必要 Appendix 24 Fortranソースパラメータの変更
例）JMファントム（ボクセル数：約1.5億個）を6分割（最大セルのボクセル数：約4千万個）にした場合 Fortranソースパラメータの変更 param.inc integer*8 mdas,mcmx,mci,mmdas,mmmax,nbnds,mct ! 4バイト整数の最大値（ ）超え対策 parameter ( mdas = ) ! PHITSで使用する総メモリ数（Byte)/8, 64bit機の導入必須 parameter ( latmax = ) ! 1セルあたりの最大ボクセル数 angel00.inc integer*8 mdas,mmdas,mmmax ! 4バイト整数の最大値（ ）超え対策 parameter ( mdas = ) ! ANGELで使用する総メモリ数（Byte)/8 コンパイラオプションの追加（Intel Fortranの場合） -i-dynamic : Libraryを動的リンク -mcmodel=large : メモリの使用制限なし ただしLinux版のみ対応なので，巨大ファントムを利用する場合はLinuxシステムが必須 makefile F = ifort FCFLAGS = -noautomatic -mcmodel=large -i-dynamic Appendix 24

25 補足資料② lattice.inp Lattice 25
Voxelファントムの向きを変えるには，基本単位となるvoxel を定義する際の順番を変える必要がある lattice.inp RPP, BOX と同じ [ S u r f a c e ]　（一部抜粋） 101 px -1 102 px 1 103 py -1 104 py 1 105 pz -1 106 pz 1 [ C e l l ]　（一部抜粋） $ Universe 1 lat=1 u=1 fill=-2:2 -2:2 0:0 RPPを各面 に分解 X正方向，Y正方向，Z正方向 X負方向，Y正方向，Z正方向 順番が 重要！ 先に書く面 が正方向 X正方向，Y負方向，Z正方向 X負方向，Y負方向，Z正方向 Lattice 25