300 likes | 728 Views
P HI T S. Multi-Purpose P article and H eavy I on T ransport code S ystem. PHITS Tutorial for making Voxel Phantom. Last revised 2013/10. title. 1. What is Voxel Phantom?. Reproduce a complex structure such as human body based on repeated rectangles filled with a certain material
E N D
PHITS Multi-Purpose Particle and Heavy Ion Transport code System PHITS Tutorial for making Voxel Phantom Last revised 2013/10 title 1
What is Voxel Phantom? Reproduce a complex structure such as human body based on repeated rectangles filled with a certain material (See Manual 5.7.5) You can make voxel phantom in PHITS virtual space using Universeand Latticefunctions (See Manual 5.7.3 and 5.7.4) Low resolution High resolution Introduction 2
Examples of PHITS calculations using voxel phantom Biological dose estimation for charged-particle therapy T. Sato et al. Radiat. Res. (2009) Treatment planning for BNCT H. Kumada et al. J. Phys.: Conf. Ser. (2007) Introduction 3
Table of Contents • Universe • Lattice • Simple voxel phantom • Conversion from DICOM format • Summary • Appendix Table of Contents 4
What is Universe? → Virtual space in PHITS You can define many universes in PHITS virtual space Universe1 But only 1 universe (main space) is the stage of particle transport simulation Other universes are used for replacing some parts of the main space using “fill” command Main space Some parts of the main space (inside the boxes) are filled with universe 1 Universe 5
Example of Universe always 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 Filled with Universe1 Declare universe 1 PX -3 PX 9 Figure 5.13 (a) Two rectangular solids. (b) Cylinder filled with water. (c) Iron cylinder in water Universe 6
Table of Contents • Universe • Lattice • Simple voxel phantom • Conversion from DICOM format • Summary • Appendix Table of Contents 7
What is Lattice? → Repeated structure used in PHITS virtual space It is troublesome to define all surfaces and cells used in repeated structure Define only surfaces and cells used in fundamental structure Examples of Lattice in PHITS Express the repeated structure using “lat” command Lattice 8
How to define lattice? You cannot directly define the contents inside lattice Only repeated structure can be defined in lattice universe It is better to define lattice not in main space but in a universe You have to fill lattice with other universe Define repeated structure using more than 2 universes fill fill Universe1 (Lattice structure) Universe2 (fundamental structure) Main space Lattice 9
PHITS input 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) -5Y5 Basic lattice(0,0,0) (-2,-2,0) -5X5 Region 101 Declare lattice type 1 (Rectangle) Define the region of basic lattice Define the number of repeated structure Universe number to be filled with(5×5×1 matrix) Location should be adjusted to that of the basic lattice Lattice 10
Change the contents of lattice lattice1.inp Change 1st box from golden ball to 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 After (lattice1.inp) Before (lattice.inp) Lattice 11
Table of Contents • Universe • Lattice • Simple voxel phantom • Conversion from DICOM format • Summary • Appendix Table of Contents 12
How to define voxel phantom? ① Make universes filled with an unique material such as bone and soft tissue Universe1 (void) Universe2 (water) Universe3 (Aluminum) ② Make voxel phantom by repeating those universes ③ Fill some part of the main space with the voxel phantom Universe10 (Voxel Phantom) Main Space Simple Voxel Phantom 13
PHITS input file robot.inp Surfaces for the basic lattice [ S u r f a c e ] $ fundamental voxel 1 px -5 2 px -3 3 py 3 4 py 5 5 pz 3 6 pz 5 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 … repeat 4 times $ 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 Any large region is OK Lattice order: X+, Y+, Z+ (start with left&lower voxel) z y x Simple Voxel Phantom 14
Change materials robot1.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 ... last one 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 … repeat 4 times $ 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 Change the material of the phantom head from water to copper Before (robot.inp) After (robot1.inp) Simple Voxel Phantom 15
Example of dose calculation 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 [t-heat] tallyusing mesh = reg Calculate dose for each region (Head, torso, and arm&leg) [t-heat] tally using mesh = xyz Useful for calculating dose inside tumor region Visualize the dose distribution Simple Voxel Phantom 16
Table of Contents • Universe • Lattice • Simple voxel phantom • Conversion from DICOM format • Summary • Appendix Table of Contents 17
DICOMformat (Binary) Data for 1 slice (sample001.dcm) ① Header (Information on time, voxel size etc.) ② CT values(1,1→2,1→3,1→…→nx-1, ny → nx, ny) Several files are contained in one folder to represent an object 3D view cross sectional view It is necessary to convert from DICOM to PHITS-input format (CT value, binary) (Universe number, text) Dicom to PHITS 18
DICOM2PHITS • Convert from Dicom data to PHITS input format (voxel phantom) • Relation between CT value and material density is defined in HumanVoxelTable.data* 1.Make an input file for DICOM2PHITS (dicom2phits.inp) 2.Windows: Drag dicom2phits.inp and drop into dicom2phits.bat Mac:Copy & double click dicom2phits.command and input dicom2phits.inp 3.PHITS input file “phits.inp” and voxel data “voxel.inp” are automatically generated "HumanVoxelTable.data" ! File for conversin of human voxel data "DICOM/" ! DICOM file directory "sample010.DCM" ! Sample dicom filename 7 9 ! Start and End position of numbering in the filename 10 29 ! Minimum slice number, Maximum slice number 70 430 90 460 ! Clipping: Nxmin, Nxmax, Nymin, Nymax 4 4 1 ! Coarse graining: Nxc, Nyc, Nzc 1 ! Z direction +1:forward or -1:backward "phits.inp" ! Filename for PHITS input "voxel.inp" ! Filename for include file You can pick up only important region Coarse graining (reduce the resolution) is feasible 19 *Data are based on W. Schneider, Phys. Med. Biol. 45, 459 (2000)
Generated Input File phits.inp file=phits.inp [Parameters ] … [ C e l l ] $ Material universe 8001 8001 -0.00121 -99 u=8001 … (24 materials are specified) $ Voxel universe 8100 0 -8100 lat=1 u=8100 fill= 0: 89 0: 91 0: 19 infl: { voxel.inp } $ Main space 97 0 -97 trcl=800 fill=8100 98 0 -99 #97 $ Void 99 -1 99 $ Outer region Always 1st line. Necessary for using “infl” command [ S u r f a c e ] 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 … $ Unit voxel 8100 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 Size and number of voxels are automatically adjusted Transform is feasible Center of voxel phantom is located at the origin of the XYZ coordinate Include file command Voxel data file generated by DICOM2PHITS is specified Dicom to PHITS 20
Geometry check icntl = 11 icntl = 8 CT3D.eps heat-xy.eps Dicom to PHITS 21
Change voxelized region dicom2phits.inp "HumanVoxelTable.data" ! File for conversin of human voxel data "DICOM/" ! DICOM file directory "sample010.DCM" ! Sample dicom filename 7 9 ! Start and End position of numbering in the filename 10 19 ! Minimum slice number, Maximum slice number 70 270 90 270 ! Clipping: Nxmin, Nxmax, Nymin, Nymax 4 4 1 ! Coarse graining: Nxc, Nyc, Nzc 1 ! Z direction +1:forward or -1:backward "phits1.inp" ! Filename for PHITS input "voxel1.inp" ! Filename for include file After change (phits1.inp) heat-xy.eps CT3D.eps Dicom to PHITS 22
Reduce computational time Purpose Every time PHITS runs… It converts its input file to binary, and re-reads the binary file It is better to… Make binary file of voxel phantom prior to the PHITS execution Procedure ① Insert the following 2 lines in the [Parameters] section ivoxel = 2 # Convert the “fill” part of lattice to binary and output to file(18) file(18) = voxel.bin # Output file name for binary voxel phantom ② Execute PHITS →Binary file was successfully generated!! ③ Change “ivoxel = 1”, and comment out “infl” command ivoxel = 1# Read the “fill” part of lattice from file(18) Speed up! $ infl:{voxel1.inp} Dicom to PHITS 23
Example of dose calculation icntl = 0 To calculate the absorbed dose in Gy, you have to estimate the mass of the region i.e. number of voxels assigned to the region (See Manual Chapter 7) heat-xy.eps [t-heat] tally using mesh = xyz Dicom to PHITS 24
Table of Contents • Universe • Lattice • Simple voxel phantom • Conversion from DICOM format • Summary • Appendix Table of Contents 25
Summary • Voxel phantom can be implemented in PHITS using Universe and Lattice concepts • DICOM format must be converted into PHITS input format using DICOM2PHITS • Computational time can be reduced by using “ivoxel” parameter Summary 26
Table of Contents • Universe • Lattice • Simple voxel phantom • Conversion from DICOM format • Summary • Appendix Table of Contents 27
How to Deal with High-Resolution Phantom? High resolution voxel phantom requires numerous memory Default setting of PHITS is allowed to use memory only less than 160 MByte e.g.) Whole body voxel phantom (180cm×30cm×50cm) with 1mm3 resolution consists of 270,000,000 voxels, and costs 5.4GByte memory, since PHITS uses memory approximately 20 Byte / voxel How to deal with the situation? Change “param.inc” included in “src” folder • increase mdas: Maximum memory allowed to be used by PHITS (Byte) / 8 • increase latmax:Maximum number of lattice in a cell • declare integer*8 for several parameters (see next page in detail) Delete all object files (*.o) and re-compile PHITS* Memory is insufficient? Divide voxel phantom into several regions to reduce the area to be voxelized Combine several CT pixels into one voxel to decrease the resolution of phantom *For Windows PC, gfortran is recommended to be used for this purpose, because PHITS executable file compiled by Intel Fortran may cause “stack overflow” for large voxel phantom 28
How to Deal with High-Resolution Phantom? If #voxels is greater than 50 millions, many changes are necessary e.g. total #voxel = 150 millions, max #voxel per cell = 40 millions Change include files in “src” folder param.inc integer*8 mdas,mcmx,mci,mmdas,mmmax,nbnds,mct ! avoid overflow (integer*4 =< 2147483647) parameter ( mdas = 500000000 )! Maximum memory allowed to be used by PHITS (Byte) / 8 parameter ( latmax = 47000000 ) ! Maximum number of lattice in a cell angel00.inc integer*8 mdas,mmdas,mmmax ! avoid overflow (integer*4 =< 2147483647) parameter ( mdas = 350000000 )! Maximum memory allowed to be used by ANGEL (Byte) / 8 Add compiler options (e.g. for Intel Fortran in Linux) makefile F77 = ifort FCFLAGS = -noautomatic -mcmodel=large -i-dynamic -i-dynamic: Dynamic link to libraries -mcmodel=large: no limitation in memory use (this option is only valid for Linux) Appendix 29
Change the order of lattice lattice2.inp [ S u r f a c e ](pick up partially) 101 px -1 102 px 1 103 py -1 104 py 1 105 pz -1 106 pz 1 [ C e l l ](pick up partially) $ 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 Same as RPP, BOX RPP is divided into each surface -102101 -104 103 -106 105 X:+, Y:+, Z:+ 101-102 -104 103 -106 105 X: –, Y:+, Z: + Order is important! Prior surface faces to the forward direction -102101 103 -104 -106 105 X: +, Y: –, Z:+ 101-102 103 -104 -106 105 X: –, Y: –, Z: + Appendix 30