Transcript Beam source Model

11 of February, 2014
Principle of
convolution/superposition
algorithm for dose
calculation for Treatment
Planning System
Woong Cho
[email protected]
What is RTP System?
 Radiation Treatment Planning
System
 Eclipse, Pinnacle,
RayStation, Xio ….
 CorePlan (?)
 Determining treatment
parameters to deliver acceptable
radiation dose to the patient




Contouring
Treatment simulation (planning)
Dose calculation
Evaluation of plans
Architecture of general RTP System
Patient DB
• Study information
• Images(CT, MRI…)
• Planning Data
• Dose Data
RTP Main Module
• Patient Data I/O
• 3D Visual rendering
• Contouring
• Planning GUI
• Dose calculation
• Planning result analysis
ETC module
• DVH calculator
• Isodose generator
• NTCP/TCP calculator
Beam Data DB
• Linac geometry info
• Dosimetric Data
•PDD data
• Profile data
• Output factor
• Wedge/MLC info
• CT HU-density table
• ETC…
• Kernel table (hidden)
Commissioning module
• RTP commissioning
• Data I/O interface
• Spectrum optimization
• Source modeling optimization
• Determine output factor:
value/MU, cGy/MU
• Wedge shape…
IMRT Module
• Beam-let extraction
• Make Voxel lists per organ
• Constraint I/O
• Run optimization
• MLC sequencing
Dose calc Engine
• Photon Dose calculation
- ETAR , CCC, AAA…
• Electron Dsoe
- PBC Hogstrom algorithm
Optimization engine
• Steepest decent method
- for spectrum
- for beam source model
- for IMRT
Dose calculation engines
 Dose calculation engine
 Most important to predict dose distribution
 Consisted of two parts
 Dose calculation part
 Interpolation based method: ETAR
 Semi-analytic model method
 Superposition/Convolution algorithm
 FFT, CCC, AAA…
 LBTE algorithm
 Acuros XB
 Monte Carlo based method
 Beam modeling part
 Based on measured dosimetric data
 Directly linked to RTP commissioning module
Convolution/superposition method
 Two step of the process at depositing doses
 Photon and a media interaction making TERMA
 Electrons from a interaction point are traversed through the media
making excitation and ionization (kernel)
 TERMA
 Total energy released per media
 Temporal energy deposition before electron transport
 Considering attenuation and divergence effect
 Kernel
 Relative energy spread from unit TERMA
 Energy deposition from electron transport
 Function of energy and media
 Dose = TERMA ⓧ kernel
D(r )   Tp (r ')  ker nel (r  r ' ) d 3 r '
Convolution/superposition method
 TERMA deposition:
 Pencil beam transfer energy at each
voxels in tissue through its
traveling way.
 Assume that directly give potential
energy to each volume element
 Kernel spread:
 Deposited dose spread to 3D space
from TERMA at each voxel.
 Dose:
D(r )   Tp (r ') ker nel (r  r ' ) d 3 r '
 Superposition: Sun of all kerma
map at each voxel.
 Same as convolution of TERMA
and kernel
General Process of Dose Calculation
• Consider geometric infomations
Calculate
Photon fluence
• Photon beam source
•
• 2 source model
•
• 3 source model
• Beam aperture
•
• Collimators
• MLC
• Horn effects
• Beam Divergence
• Inverse square law
Calculate
3D TERMA
Do convolution
Poly-energetic kernels
Attenuation
Beam
Kernel(r , ,  ) 
hardening/softening
all Energy
Transmission through  Monoker nel(r , ,  , Energy)
MLC /block/collimators • C/S or CCC or AAA
• Differential kernels or
Accumulative kernels
Geometric transformation
Source
• Geometric definition
Beam coordinate
yb
xb
•
Freedom of beam directions
•
•
•
•
zb
0
Block
Block
1
•
Beam coordinate
Transformation between Patient
coordinate and Beam coordinate
xb
CT coordinate
yb
xg
zb
zg
•
Source
0
1
isocenter
yg
•
Gantry angles
Collimator angles
Couch angles
Translations of couch
Patient coordinate system
•
•
•
•
World (global) coordinators
Based on CT coordinators
Define 3D dose
Define contours
Beam coordinate system
•
•
•
•
Origin is beam source
Z axis is beam center
Fluence. Beam source, TERMA,
kernels
Defined per beams
Geometric transformation
• Simple Example
•
yb
zb
Source
Sb ( 0, 0, 0)
Beam coordinate
yb
•
xb
zb
•
yg
zg
yb
zb
yb
CT Origin
Pg (0,0,0)
Beam coordinate
5zb
xb
Isocenter
Beam coordinate
5
xb
Pb
xg
Pb at iso-center ( Xb, Yb, Zb)
= (0, 0, 100)
Assume Iso-center (Pg ) is (0, 5, 5)
Converting Pb to Pg?
• Rotate based on Xb axis by 90’ counter
clock wise
• Pb = (0,-100, 0)
• Translate Beam coordinate by (0, -100, 0)
• Beam coordinate origin shifted to
isocenter
• Pb = ( 0,0,0)
• Translate Beam coordinate by (0, -5, -5)
• Beam coordinate is same to patient
coordinate
• Pb = (0, 5, 5) = Pg
 x g  cos  cos  cos  sin   sin   sin   cos sin   sin   cos  sin   cos   xb  t x 
  
   
 y g    cos  sin  cos  cos  sin   sin   sin   sin   cos  cos  sin   sin    yb   t y 
 z    sin 
  z b  t z 
sin   cos
cos  cos
 g 
Process of TERMA Calculation
• Calculation of Fluence and TERMA distribution in media
• Beam source model
Calculate
Fluence at each voxel • Binary MLC plane
• Using Hit test
Calculate
• Considering
Beam divergence
partially block
Source
0
1
Calculate Horn effect
Calculate effective depth
for each voxels
Calculate
beam softening effect
Calculate Attenuation
Calculation
point
Beam source Model: 3 source model
Gantry head
Point source
Srcprimary plane
Zsp = 4cm

Point source: Cp

Srcsp plane
Srcprimary (r )  C p
Disk source
Zsf = 12.5 cm
Srcsf plane
Collimator
MLC
Primary photon source
Annulus source
Collimator
MLC
or

or
Scattered photon source from
primary collimator
Annulus shape

Srcsp (r )  Csp , ( R01  r  R02 )
0 ,
otherwise
SCD = 100 cm

Isocenter
Scattered photon source from
other structures

Point of calculation
(xb,yb,zb)

Disk shape
Intensity function with r
A0
Src sf (r )  exp( k  r )
r
Beam source model: 3 source model
Gantry head
Point source
Srcprimary plane
Annulus source
Srcsp plane
Disk source
Srcsf plane

Collimator
MLC
Collimator
MLC
or
or
Fluence at an arbitrary point
C
 C
Fluence

Isocenter
Point of calculation
(xb,yb,zb)

p
 dAsrc 1  f ISW _ src 1
sp
(rsp )  dAsrc 2  f ISW _ src 2
A0
exp(k  rsf )  dAsrc 3  f ISW _ src 3
r
Implementation of Beam source model
Beam Source center (0, 0, 0)
Binary Block Plane
1 1
Fluence
Voxel
(xb, yb, zb)
0
• Define Binary MLC Grid
from the position of MLC
leaves
– Considering partially
block
– Hit Test algorithm
Prepare Binary 2DGrid
For all voxels (xb,yb,zb)
{
Src1_fluence =
Calc_OpenRatio(Subboxels)
{
for (3x3x3 Subvoxels)
do HitTest(SubVoxels)
}
Src2_fluene
=Calc_ScatterSource2()
Src3_fluene
=Calc_ScatterSource3()
}
Implementation of Beam source model
1
0
Xs
Ys
0
Source plane
Z = Z_src
5mm resolution
0
Hit_corner_grid:
partially block
MLC plane
Z= 67 cm
0.1mm resolution
Sub voxel
Blocking ratio=
9/25 = 0.36
(Xb,Yb,Zb)
Calculation point:
101 x 101 x 101
Src2_fluene =Calc_ScatterSource():
For all SourceGrid
{
Hit-test 4 corner bloks at first
If not blocked
Calculate Radius (Xb,Yb,Zb)
Get source value from source
functions
Src_flue += source_value
}
Process of TERMA Calculation
• Calculation of Fluence and TERMA distribution in media
Calculate
Fluence at each voxel
Calculate
Beam divergence
Calculate Horn effect
Source
FluenceDiv
 FluenceInit
 Distref ( 100cm) 

 
SPD


0
2
1
SPD
Distref
Calculate effective depth
for each voxels
Calculate
beam softening effect
Calculate Attenuation
Calculation
point
Process of TERMA Calculation
• Calculation of Fluence and TERMA distribution in media
Source
Calculate
Fluence at each voxel
0
Calculate
Beam divergence
Calculate Horn effect
1
FluenceHorn
 HornF(OAD)

 FluenceDiv  
 OAD
100


Calculate effective depth
for each voxels
Calculate
beam softening effect
Calculate Attenuation
OAD
Calculation
point
Process of Dose Calculation
• Calculation of Fluence and TERMA distribution in media
Source
Calculate
Fluence at each voxel
0
Calculate
Beam divergence
1
Calculate Horn effect
Calculate effective depth
for each voxel
d Effective

Calculate
beam softening effect


r r
r  source
r  voxel
  (r )  r
ri  source
Calculate Attenuation
 (r )d r
i
Calculation
point
Process of Dose Calculation
• Calculation of Fluence and TERMA distribution in media
Source
Calculate
Fluence at each voxel
0
Calculate
Beam divergence
1
Calculate Horn effect
Calculate effective depth
for each voxels
Calculate
beam softening effect
Calculate Attenuation
d soften (OAD)

1
1  f SofteningRatio  OAD
 d Effective
Calculation
point
Process of Dose Calculation
• Calculation of Fluence and TERMA distribution in media
Source
Calculate
Fluence at each voxel
0
Calculate
Beam divergence
1
Calculate Horn effect
Calculate effective depth
for each voxels
Calculate
beam softening effect
•
Beam hardening
effect
TERMAmono ( E )  FluenceHorn
 exp(
Calculate Attenuation
TERMA 


( E )  d soften )  E  ( E )


E max
 TERMA
 E min
mono
(E)
Implementation of convolution
Photon
source
D(r )   TERMAp (r ') Kernel (r  r ' ) d 3 r '

N
 TERMA (r ' )  Kernel (r  r ' )  V
r ' 1
p
for all r’…(N voxels in volume)
{
D(r)
T(r’)
Get TERMA(r’)
r
r'
for all r….
{
Get Kernel (r-r’) from Table
Accumulate TERMA(r’) x Kernel (r-r’) to the r’ voxel
}
}
Why Collapse Cone Convolution?
 Limitation of convolution/superposition method
 Center of kernel has too steep
gradient.
 Discrete kernel data  Significant
error at the center voxels in dose
calculation.
 Too long calculation time
 FFT is a good method. But no inhomogeneity correction.
 Not invariant kernel at inhomogeneous medium
 Iterative calculation: N6 number of iteration at N x N x N voxels
From “Current Concepts in Dose Calculations”, Anders Ahnesjö.
Collapsed cone approximation
M Number of
Cones
N
N voxels
M rays
N number
of Voxels
N
No of Iterations: N3
No of Iterations : M x N
• Can reduce calculation time because of computing dose
to MxN points instead of NxNxN points.
• More accurate dose calculation in heterogeneous media
by considering effective pathway through cone lines
Scatter particle transport directions
 Near center voxels
 Spherical voxels are generally smaller than cubic voxels
 Far away voxels
 One spherical voxels covers several cubic voxels
 Only consider the voxels in axial lines
 Too much energy imparted.
 But small errors because of small fraction energy at far site.
Process of Convolution
• Sphere Convolution Process
• Total 288 cone rays
A Collapsed
Cone
– 24 divisions of theta
Y
– 12 divisions of phi angle
θ
• Extracting voxel lists traversed by
each ray vector r.
• Heterogeneity correction by effective
X
pathway through vector r
φ
A TERMA Voxel
• Convolution kernel table with
-Z
spherical coordinate system.
• Consideration of kernel tilting effect
• Adaption of accumulative kernel

D( x ) 
 
 D (x  r )
288
 1


sub



 
 



TERMA
(
r
)

k
(
x

r
)

k
(
x

r
)


(
r
)V

primary
scatter
288
 1
Process of convolution
Prepare Poyenrgetic_Spherical_Kernel (r, θ, φ)
Make Accumulative_Kernel (r, θ, φ)
For (all 3D voxels with (Xp,Yp,Zp) )
{
for (all r, θ, φ)
{
Calc_vector()
Get_transversed_voxel_Lists()
Calculate eff_pathlength(voxel_Lists)
for (all Listed voxels)
{
Calculate θtilt
Energy = Accumulative_Kernel(rinner)
- Accumulative_Kernel(router)
Get_TERMA(voxel)
Dose +=TERMA x Energy
}
}
}
(r, θ, φ)
(Xp, Yp, Zp)
Kernel tilting
 
rg  rsrc
TERMA
Voxel

rg
θbeam
Divergent
Beam
 tilt   beam   cone
θcone
Dose
Voxel


riso  rsrc
 


(rg  rsrc )  (riso  rsrc )
 arccos(  

 )   cone
(rg  rsrc )  (riso  rsrc )
Considering Beam hardening
• Poly-energetic kernel
– Photon spectrum is changed according to depths
• Solve:
– Get changed spectrums at every 10 cm depth
– Prepare each kernel tables from the changed
spectrum
– Calculate interpolate kernel values between two
tables using the depth of voxels
Discussion & Question