PHITS

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
講習の流れ
1. Universe
2. Lattice
3. 簡易Voxelファントム
4. Dicom形式からの変換
5. まとめ
Table of Contents
4
Universeとは？
→ PHITSの中の仮想空間
Universe1
PHITSでは，いくつもの仮想空
間（universe）を定義できる
実際に粒子輸送の舞台とな
るのはメイン空間のみ
メイン空間の一部を別の
universeで満たすことにより
入れ子構造を再現する
メイン空間
メイン空間にある箱の中を，球が並ん
でいる別の宇宙(universe1)に置換
Universe
5
満たす領域は常にvoid
universe.inp
Universeの例
[Cell]
Universe1で満たす
$ Main space
1 0 11 -12 13 -14 15 -17 FILL=1 （Fillする）
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
Universe1を宣言
$ Universe 2
201 2 -7.86 -10 13 -14 U=2
202 1 -1.00 #201 U=2
[Surface]
10 CY 5
PX -3
11 PX -6
12 PX 6
PX 9
13 PY -6
14 PY 6
15 PZ -6
16 PZ 6
17 PZ 0
Universe
6
講習の流れ
1. Universe
2. Lattice
3. 簡易Voxelファントム
4. Dicom形式からの変換
5. まとめ
Table of Contents
7
Latticeとは？
→ PHITSの中で使う
格子(lattice，繰り返し）構造
同じ形状のものが繰り返し並ん
でいるときに，個々の形状を全
て定義するのは大変！
基本構造のみ定義して，あと
はその繰り返しで表現する
Latticeの例
格子には直方体と六角柱が
あるが，基本的に使うのは直
方体のみ
Lattice
8
格子の定義方法
格子空間には他の構造物は定
義できない
格子の中身は直接定義できない
メイン空間ではなく，universeの
１つとして定義するのが便利
別のuniverseで満たす(fillする）
必要がある
Universeを２つ以上使った２重入れ子構造にする
fill
メイン空間
fill
Universe1
Lattice
Universe2
9
5
（2,2,0）
Y軸
[Surface]
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
[Cell]
$ 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
22222
22222
22222
22222
22222
$ Universe 2
201 1 -19.32 -201 u=2
202 0 201
u=2
格子の定義
基本格子（0,0,0）
（-2,-2,0）
-5
lattice.inp
-5
直方体タイプのLatticeを宣言
X軸
5
領域101
基本格子の領域を定義
格子に物質を満たす領域を指定（0は基本格子）
各格子に満たすuniverse番号（5×5×1の行列）
基本格子と座標を合わせる必要有
メイン空間で使う領域は，必ず定義しておく必要有
Lattice
10
格子の中身の変更
lattice.inp
[Cell]
$ 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
32222
22222
22222
22222
22222
$ Universe 2
201 1 -19.32 -201 u=2
202 0 201
u=2
$ Universe 3
302 0 -99 u=3
１つ目の格子のみVoidに変更
変更前
(lattice.inp)
Lattice
変更後
(lattice1.inp)
11
講習の流れ
1. Universe
2. Lattice
3. 簡易Voxelファントム
4. Dicom形式からの変換
5. まとめ
Table of Contents
12
簡易Voxelファントムの作り方
① 単一の物質（骨・軟組織など）
Universe1 Universe2
で満たされたuniverseを作る
(void)
② ①で作ったuniverseを
Latticeを使って組み合わせ，
ボクセルファントムが入った
universeを作る
(water)
Universe3
(Aluminum)
③ ②で作ったuniverseを
メイン空間の一部にFillする
Universe10
(Voxel Phantom)
メイン空間
Simple Voxel Phantom
13
インプットファイル
robot.inp
[Cell]
$ 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
11111
12121
12121
11111
11111
… 以下，４回繰り返し
$ 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正方向
（図の左下から始まる）
[Surface]
$ 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
z y
x
Simple Voxel Phantom
14
物質の追加・変更方法
robot.inp
[Cell]
$ 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ブロック）
11111
11111
11411
11111
11111
$ 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)
Simple Voxel Phantom
変更後
(robot1.inp)
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
講習の流れ
1. Universe
2. Lattice
3. 簡易Voxelファントム
4. Dicom形式からの変換
5. まとめ
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/"
"PHITSinputs/"
1 20
70 430 90 460
4 4 1
0
0
1
変換テーブルの指定
指定ディレクトリ内に含まれるDICOMファイルを自動判別
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
読込の高速化
目的
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
講習の流れ
1. Universe
2. Lattice
3. 簡易Voxelファントム
4. Dicom形式からの変換
5. まとめ
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の場合）
makefile
F77
= ifort
FCFLAGS = -noautomatic -mcmodel=large -i-dynamic
Appendix
-i-dynamic : Libraryを動的リンク
-mcmodel=large : メモリの使用制限なし
ただしLinux版のみ対応なので，巨大ファン
トムを利用する場合はLinuxシステムが必須
24
補足資料②
Voxelファントムの向きを変えるには，基本単位となるvoxel
を定義する際の順番を変える必要がある
lattice.inp
[ S u r f a c e ] （一部抜粋）
101 px -1
102 px 1
RPPを各面
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
重要！
32222
22222
22222
先に書く面
22222
22222
が正方向
RPP, BOX
と同じ
-102 101 -104 103 -106 105
X正方向，Y正方向，Z正方向
101 -102 -104 103 -106 105
X負方向，Y正方向，Z正方向
-102 101 103 -104 -106 105
X正方向，Y負方向，Z正方向
101 -102 103 -104 -106 105
X負方向，Y負方向，Z正方向
Lattice
25

Download Report