Download
1 / 25
250 likes | 558 Views
P HI T S. Multi-Purpose P article and H eavy I on T ransport code S ystem. PHITS 講習会 Voxel ファントムの作り方. 2014 年 8 月改訂. title. 1. Voxel ファントムとは？. 直方体を積み重ねて生物など複雑な体系を再現したもの （マニュアル 4.7.5.3 ）. PHITS の Universe と Lattice 構造を使って構築 （マニュアル 4.7.3 ＆ 4.7.4 ）. 低分解能. 高分解能. Introduction. 2.
E N D
PHITS Multi-Purpose Particle and Heavy Ion Transport code System PHITS 講習会 Voxelファントムの作り方 2014年8月改訂 title 1
Voxelファントムとは？ 直方体を積み重ねて生物など複雑な体系を再現したもの （マニュアル4.7.5.3） PHITSのUniverseとLattice構造を使って構築 （マニュアル4.7.3＆4.7.4） 低分解能 高分解能 Introduction 2
Voxelファントムを用いた計算例 粒子線治療時の２次放射線被ばく線量評価 T. Sato et al. Radiat. Res. (2009) BNCTの治療効果推定 H. Kumada et al. J. Phys.: Conf. Ser. (2007) Introduction 3
講習の流れ • Universe • Lattice • 簡易Voxelファントム • Dicom形式からの変換 • まとめ Table of Contents 4
Universeとは？ → PHITSの中の仮想空間 PHITSでは，いくつもの仮想空間（universe）を定義できる Universe1 実際に粒子輸送の舞台となるのはメイン空間のみ メイン空間の一部を別のuniverseで満たすことにより入れ子構造を再現する メイン空間 メイン空間にある箱の中を，球が並んでいる別の宇宙(universe1)に置換 Universe 5
Universeの例 満たす領域は常にvoid universe.inp [ C e l l ] $ Main space 1 0 11 -12 13 -14 15 -17 FILL=1 2 0 11 -12 13 -14 17 -16 FILL=2 9 -1 #1 #2 $ Universe 1 101 1 -1.00 -10 13 -14 U=1 102 0 #101 U=1 $ Universe 2 201 2 -7.86 -10 13 -14 U=2 202 1 -1.00 #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
講習の流れ • Universe • Lattice • 簡易Voxelファントム • Dicom形式からの変換 • まとめ Table of Contents 7
Latticeとは？ → PHITSの中で使う 格子(lattice，繰り返し）構造 同じ形状のものが繰り返し並んでいるときに，個々の形状を全て定義するのは大変！ 基本構造のみ定義して，あとはその繰り返しで表現する 格子には直方体と六角柱があるが，基本的に使うのは直方体のみ Latticeの例 Lattice 8
格子の定義方法 格子の中身は直接定義できない 格子空間には他の構造物は定義できない メイン空間ではなく，universeの１つとして定義するのが便利 別のuniverseで満たす(fillする）必要がある Universeを２つ以上使った２重入れ子構造にする fill fill Universe2 Universe1 メイン空間 Lattice 9
格子の定義 lattice.inp [ S u r f a c e ] 1 rpp -5 5 -5 5 -1 1 2 rpp -6 6 -6 6 -2 1 99 so 100 101 rpp -1 1 -1 1 -1 1 201 sph 0 0 0 1 [ C e l l ] $ Main space 1 3 -8.96 1 -2 2 0 -1 fill=1 98 0 -99 2 99 -1 99 $ Universe 1 101 0 -101 lat=1 u=1 fill=-2:2 -2:2 0:0 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 $ Universe 2 201 1 -19.32 -201 u=2 202 0 201 u=2 （2,2,0） -5Y軸 5 基本格子（0,0,0） （-2,-2,0） -5X軸 5 領域101 直方体タイプのLatticeを宣言 基本格子の領域を定義 格子に物質を満たす領域を指定（0は基本格子） 各格子に満たすuniverse番号（5×5×1の行列） 基本格子と座標を合わせる必要有 メイン空間で使う領域は，必ず定義しておく必要有 Lattice 10
格子の中身の変更 lattice.inp １つ目の格子のみVoidに変更 [ C e l l ] $ Main space 1 3 -8.96 1 -2 2 0 -1 fill=1 98 0 -99 2 99 -1 99 $ Universe 1 101 0 -101 lat=1 u=1 fill=-2:2 -2:2 0:0 3 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 $ Universe 2 201 1 -19.32 -201 u=2 202 0 201 u=2 $ Universe 3 302 0 -99 u=3 変更後 (lattice1.inp) 変更前 (lattice.inp) Lattice 11
講習の流れ • Universe • Lattice • 簡易Voxelファントム • Dicom形式からの変換 • まとめ Table of Contents 12
簡易Voxelファントムの作り方 ① 単一の物質（骨・軟組織など）で満たされたuniverseを作る Universe1 (void) Universe2 (water) Universe3 (Aluminum) ② ①で作ったuniverseをLatticeを使って組み合わせ，ボクセルファントムが入ったuniverseを作る ③ ②で作ったuniverseを メイン空間の一部にFillする Universe10 (Voxel Phantom) メイン空間 Simple Voxel Phantom 13
インプットファイル robot.inp 基本格子（この場合は１つ目に定義する領域）の面 [ S u r f a c e ] $ fundamental voxel 101 rpp -5 -3 -5 -3 -5 -3 99 so 100 $ Main space 201 rpp -5 5 -5 5 -5 5 202 rcc 0 0 -5 0 0 4 8 203 rcc 0 0 -6 0 0 5 9 [ C e l l ] $ Material universe 1 0 -99 u=1 2 1 -1.00 -99 u=2 3 2 -2.70 -99 u=3 $ Voxel universe 101 0 -101 lat=1 u=10 fill=0:4 0:4 0:4 1 1 1 1 1 1 2 1 2 1 1 2 1 2 1 1 1 1 1 1 1 1 1 1 1 … 以下，４回繰り返し $ Main space 201 0 -201 fill=10 202 0 201 -202 203 3 -8.96 202 -203 204 0 -99 201 203 205 -1 99 十分に大きい領域であればどこでもよい X正方向，Y正方向，Z正方向 （図の左下から始まる） z y x Simple Voxel Phantom 14
物質の追加・変更方法 robot.inp ファントムの頭を水から銅に変更する [ C e l l ] $ Material universe 1 0 -99 u=1 2 1 -1.00 -99 u=2 3 2 -2.70 -99 u=3 4 3 -8.96 -99 u=4 $ Voxel universe 101 0 -2 1 3 -4 5 -6 lat=1 u=10 fill=0:4 0:4 0:4 … （始めの4ブロック） 1 1 1 1 1 1 1 1 1 1 1 1 4 1 1 1 1 1 1 1 1 1 1 1 1 $ Main space 201 0 -201 fill=10 202 0 201 -202 203 3 -8.96 202 -203 204 0 -99 201 203 205 -1 99 変更前 (robot.inp) 変更後 (robot1.inp) Simple Voxel Phantom 15
線量計算結果 robot-heat-xz.eps robot-heat-reg.out x: Serial Num. of Region y: Heat [MeV/source] h: x n n y(total),l3 n # num reg volume heat r.err 1 2 1.0000E+00 9.5978E-01 0.1277 2 3 1.0000E+00 3.4847E+01 0.0000 3 4 1.0000E+00 5.1924E+01 0.0000 regメッシュを用いた [t-heat]計算 各領域別（手足・胴・頭）の 線量を計算 xyzメッシュを用いた [t-heat]計算 腫瘍部分や重要臓器などに ROIをかけた線量計算が可能 ファントム内の線量分布を可視化 Simple Voxel Phantom 16
講習の流れ • Universe • Lattice • 簡易Voxelファントム • Dicom形式からの変換 • まとめ Table of Contents 17
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
変換プログラム(dicom2phits） Dicom形式のデータをPHITS形式のボクセルデータに変換 DICOM2PHITSの使い方の詳細は「DICOM2PHITSの使い方」を参照 phits/utility/dicom2phits/phits-lec-dicom2phits-jp.ppt ①入力ファイルを作成（dicom2phits.inp） "data/HumanVoxelTable.data" 変換テーブルの指定 "DICOM/" 指定ディレクトリ内に含まれるDICOMファイルを自動判別 "PHITSinputs/" PHITSインプットを格納するディレクトリを指定 1 20 指定した番号範囲のスライスファイルを読み込む (1<=z<=20) 70 430 90 460 DICOMデータの一部を切り出してボクセル化 (70<=x<=430, 90<=y<=460) 4 4 1 画像を粗くする（分解能を下げる）ことが可能 (x方向4個、y方向4個で平均) 0 原点オプション 0:ボクセル中心 1:DICOMヘッダーから抽出 0 PHITSパラメータ最適化オプション 0:最少 1:x線治療 2:粒子線治療 1 スライス読み込む順指定 +1:昇べき -1:降べき 原点オプション1では無視 ②実行 (Windows) dicom2phits_win.batにdicom2phits.inpをドラッグ&ドロップ (Mac) dicom2phits_mac.commnadをダブルクリック、現れる窓にdicom2phits.inpと入力 ⇒ PHITSinputs/ディレクトリにサンプルインプット生成 DICOM2PHITS HowTo 19
読込の高速化 目的 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} Dicom to PHITS 20
講習の流れ • Universe • Lattice • 簡易Voxelファントム • Dicom形式からの変換 • まとめ Table of Contents 21
まとめ • Voxelファントムは，UniverseとLatticeを組み合わせて作成 • Dicom形式は，別プログラム(DICOM2PHITS）を用いてＰＨＩＴＳ入力形式（テキスト形式）に変換 • ivoxelパラメータを設定することにより読込を高速化 Summary 22
高分解能ファントム利用時の注意点 全身を高分解能でボクセル化すると使用メモリが膨大になり実行できなくなる 例）全身（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
補足資料① Voxel数が約5000千万個を超える場合は，様々な対策が必要 例）JMファントム（ボクセル数：約1.5億個）を6分割（最大セルのボクセル数：約4千万個）にした場合 Fortranソースパラメータの変更 param.inc integer*8 mdas,mcmx,mci,mmdas,mmmax,nbnds,mct ! 4バイト整数の最大値（2147483647）超え対策 parameter ( mdas = 500000000 )! PHITSで使用する総メモリ数（Byte)/8, 64bit機の導入必須 parameter ( latmax = 47000000 ) ! 1セルあたりの最大ボクセル数 angel00.inc integer*8 mdas,mmdas,mmmax ! 4バイト整数の最大値（2147483647）超え対策 parameter ( mdas = 350000000 )! ANGELで使用する総メモリ数（Byte)/8 コンパイラオプションの追加（Intel Fortranの場合） -i-dynamic : Libraryを動的リンク -mcmodel=large : メモリの使用制限なし ただしLinux版のみ対応なので，巨大ファントムを利用する場合はLinuxシステムが必須 makefile F77 = ifort FCFLAGS = -noautomatic -mcmodel=large -i-dynamic Appendix 24
補足資料② 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 101 0 -102 101 -104 103 -106 105 lat=1 u=1 fill=-2:2 -2:2 0:0 3 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 RPPを各面 に分解 -102101 -104 103 -106 105 X正方向，Y正方向，Z正方向 101-102 -104 103 -106 105 X負方向，Y正方向，Z正方向 順番が 重要！ 先に書く面 が正方向 -102101 103 -104 -106 105 X正方向，Y負方向，Z正方向 101-102 103 -104 -106 105 X負方向，Y負方向，Z正方向 Lattice 25
