- Email: [email protected]
Dose calculations for high energy electron and photon beams
Compulerized Medical Imaging and Graph,ics. Vol. 13, No. 3, pp. 241-250, Printed in the U.S.A. All rights reserved.
DOSE CALCULATIONS
1989 Copyright
0
0895-61 I l/89 $3.00 + .OO 1989 Pergamon Press plc
FOR HIGH ENERGY ELECTRON AND PHOTON BEAMS J. R. Cunningham
Physics Division, The Ontario Cancer Institute, Toronto, Canada (Received 3 I August 1988)
Abstract-The procedures of radiation therapy consist of a number of steps, one of which is the calculation of the radiation dosage pattern that will result when a particular arrangement of radiation beams or sources is applied to a patient. The purpose of this calculation is two-fold. One is to predict, as part of the dose planning process, what dose distribution can be achieved with a selected beam arrangement. The second is to record what treatment has been given so that post-treatment analysis can be carried out. Both of these are important. Key Words: Dose calculations, Electronbeams, Photon beams, Tissue inhomogeneities, Dose distributions,Beam models
may be obtained from measurements made in electron beams. The task of a 3-D dose calculation method is to produce this sort of data for beams of arbitrary crosssection, incident from any direction on a surface of any shape and to include allowance for internal inhomogeneities and to do this quickly and to an accuracy better than f5%. Figure 3 shows the problem for photons in a simplified form. A beam of finite size is irradiating an object. The dose is to be calculated at point P. Photons will come along the direct path a and interact at P, releasing energy there. Photons will also travel along path b, interact in the volume dV and send scattered photons along path c to point P, again releasing energy there. The total energy released at point P would then be the sum of the contribution due to primary photons and the total of the scattered photons that come from the whole of the irradiated volume. This can be calculated in a straightforward manner if we include photons that have been scattered no more than once.
INTRODUCTION The problem of calculating radiation dose distributions, as part of radiotherapy planning, is a very old one. The advent of the CT scanner, with its ability to gather anatomical detail, has however forced a reevaluation of the methods used. The task is a mathe-
matical one and is sufficiently complicated that it has no exact solution. The development of methods for treatment planning is therefore the search for ever better methods of making approximate calculations of absorbed dose. In the days before cotmputers, dose distributions were measured in water and isodose charts were produced from these measurements. The isodose charts were then applied to a diagram of the patient’s contour and approximate corrections were made to allow for the shape of the patient’s contour under the assumption that within the contour the patient was water equivalent. Only occasionally were corrections made for the presence of internal nonwater equivalent structures. DOSE CALCULATION
METHODS Kc&-J“=.
0 !f
.Eab
P
Figure 1 shows isodose surfaces for a rectangular 16 X 5 cm beam incide:nt on a flat homogeneous phantom. These data (1) were obtained from measurements made in a water tank. They also show some of the difficulties :in displaying three-dimensional objects. The figure shows isodose distributions in four different planes; I(a) and (b) are in principal planes, that is they contain the central ray, (c) is in a plane that is parallel to (b) but is 2 cm away, and (d) is in a cross-section at a depth of 10 cm. Figure 2 is a “perspective” view of thle same data. Similar data
Here, K is the kerma, or kinetic energy released to the medium and the first term describes the primary component; 4 is the photon fluence, p is the linear attenuation coefficient for these photons, a is the path length through the medium, (p/p) the mass interaction coefficient expresses the probability that the photons will interact at point P and &, is the 241
242
Computerized Medical Imaging and Graphics
May-June/ 1989, Volume 13, Number 3
a)
Fig. 1. Isodose curves in four planes for a beam from a cobalt unit. The planes are the principal planes a cross-section and a plane parallel to one of the principal planes. Note. From The Physics of Radiology, 4th dd. Charles C Thomas; Springfield, IL; 1983.
48K
Fig. 2. A “perspective
view” of isodose surfaces. Note. From The Physics of Radiology, 4th ed. Charles C Thomas; Springfield, IL; 1983.
Dose calculations 0 J. R.
243
CUNNINGHAM
motion and for a complete calculation of the energy absorbed, that is the absorbed dose, the energy transport by these electrons must also be considered. This makes the problem quite unsolvable by analytical means and quite basic assumptions are generally made, the most common being to assume that the electrons that are set into motion deposit their energy locally. This is equivalent to assuming absorbed dose is equal to kerma and Kin Eq. 1 could be interpreted as absorbed dose. Typically, the first term in Eq. 1, the primary, accounts for 60% to 70% of the total dose at point P, first scatter accounts for 20% to 30% of the total and multiple scatter accounts for 10% or more. These proportions are shown in Fig. 4 for cobalt-60 radiation for a range of beam sizes. Although scattered radiation constitutes the minor component of the dose at a point, it is large enough that it must be taken seriously and a three-dimensional calculation method must take its properties into account. EVOLUTION
OF CALCULATION
METHODS
A direct mathematical solution to the calculation of the dose at a point such as P of Fig. 3 would 14OOG Fig. 3. Diagram showing the path for primary and scattered radiation re,aching point P.
average energy deposite’d per interaction. The primary therefore, depends on the path between the source and the point of calculation and the material at the point of calculation. The second term is much more complicated. It describes the energy deposited at point P due to photons that are scattered from every part, dI’, of the irradiated volume. In this term, b is the path length of primary photons to the scattering volume, da(O)/& is the probability of scattering photons of energy hv through angle 0 and C;, (p’/p) and &, are the linear attenuation coefficient, mass interaction coefficient and the average energy deposited from the scattered photons. c is the path length of the scattered photons in this medium. The integration is over the whole of the irradiated volume. This term although already complicated, is solvable on a computer but it accounts for photons scattered only once. Furthermore, the whole equation applies for incident photons of only a single energy. A linear accelerator beam would., for example, consist of a spectrum of photon energies and it would be necessary to integrate both terms over the energy spectrum. In addition, it is necessary to remember that photons are indirectly ionizing particles and that when energy is lost by photons, electrons are set into
2
P 50 ,o
0
-
Total (Ist.+multiple)
2
4
6
6
scatter
IO 12 14 16 16 20 22 24 26 28 30 32 34 36 Beam radius [cm] soln
Fig. 4. Graphs showing the proportion of the dose at a point, 10 cm deep in a phantom irradiated by a cobalt unit, due to primary once scattered radiation and multiply scattered radiation.
244
Computerized Medical Imaging and Graphics
May-June/ 1989, Volume 13, Number 3
involve the Boltzmann transport equation. Transport of energy by both photons and electrons would have to be included and no solution to this equation has been found for beams of finite size even in homogeneous media. Therefore, approximate solutions must be found and the history of the development of calculation methods is a story of the search for better and faster methods of approximation. The situation for electrons is, in some ways easier since they are directly ionizing particles and they do not spread as far as do photons. The Boltzmann equation has been solved for electrons, albeit approximately, and the result is the Gaussian pencil beam model. An excellent review of these solutions has been given by Niisslin (2). The development of approximate solutions for the calculation of dose has, for both photons and electrons, taken place through the evolution of four different categories of model; matrix representation, empirical or semi-empirical beam generating functions, semi-empirical beam models and three-dimensional integral models. (a) Matrix methods The representation of the dose distribution in a radiation beam is shown in Fig. 5 where isodose curves for a beam from a cobalt unit are shown superimposed on a grid of points. The dose value at each of these points can be stored in a computer as a matrix, or table of numbers. The manipulation of such a matrix is indicated in Fig. 6 which shows 3 beams irradiating a patient. Again, a grid of points is shown and the location of these points must be determined in each beam table and the dose from each beam determined by interpolation. Adjustments for contour shape are also made at this stage usually by a procedure known as the “effective SSD method.” The matrix method is fast and is applicable to both photons and electron beams but it is limited to beam configurations for which dose measurements have been made previously and for which such data have been stored. (b) Analytical beam models The oldest of these is due to Sterling et al. (3) and the most fully developed for photon beams is due to Van de Geijn et al. (4). The format followed is that the dose from a beam at a point, such as one of those marked in Fig. 6, is given by a product of two functions. D(x,, Y,, 2,) = WB,
Zrcf)
- &~Xtr,
YB)
(2)
where (x,, , y,, , z,,) are the coordinates of a point within the patient. P(zu, z,,f) is the dose at depth zB along the
on a Cartesian grid of points. The dose at each point can be read off and stored to form a matrix representation of the dose distribution.
Fig. 5. Isodose curves superimposed
beam axis relative to its value at a reference depth zref and is the percentage depth dose. gZ,(xB, ys) is the “off centre ratio” at point xB, yE with respect to the beam axis and depth zB with respect to the contour surface. The latter function describes the dose distribution along a line perpendicular to the beam axis. Van de Geijn has further taken the second function to be itself a product of two functions. g&I?, YB) = gl.z&B)*gz.rB(YB).
(3)
Each of the manipulated quantities, P, g,, and gz, may be determined by choosing analytical functions and fitting the constants associated with them to measured data or alternatively, measured percent depth doses and beam dose profile data may be stored and manipulated by table lookup and interpolation procedures. Allowance for patient surface curvature is made by using ray tracing procedures, such as are indicated in Fig. 7, which shows a beam from a cobalt unit irradiating a curved surface. The dose is to be calculated at point P which is at a depth d’ below the patient surface as measured along a ray line and is at a depth d below the flat surface for which the beam data applies. The dose at P is determined by using d’, the radiological equivalent depth, zB in Eq. 2.
Dose calculations 0 J. R.
245
CUNNINGHAM
/
I \
\ \ .
.
Fig. 6. The dose distribution
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
from three radiation beams determined by calculating the dose at each point from each beam and adding them point by point.
t
F_ T
T h
d
t d’
1/1f
/ -
I20 J Fig. 7. Isodose curves superimposed on a ray line grid. The dose at point P (depth d) below the flat surface is about 55 and is about 65 below the curved surface (depth d’).
Computerized Medical Imaging and Graphics
246
May-June/ 1989, Volume 13, Number 3
This method can be applied to electrons as well as photons and a review of these applications is given by Sternick (5). For photon beams, the function, gZ,(xB, yB), is chosen empirically, for electrons its form has been chosen as a result of solving the transport equation which for a broad beam implies the use of an error function (6). It is interesting to note that Sterling et al. (3) chose this same form for photon beams although his choice was based purely on empirical grounds. The depth dose function, P(zB, zr,r) is taken as an empirical fit to depth dose data for both photon and electron beams. This calculation method is not three-dimensional, since the calculation depends only on the distance and material found along the ray line to the point of calculation. There is no account taken of the three-dimensional character of the irradiated volume.
y
(c) Semi-empirical methods The main example of this type of representation for photon beams is the manipulation of scatter-air ratios and for electron beams is the use of Gaussian shaped pencil beams. A problem that is at least twodimensional in nature is indicated in Fig. 8 which shows a radiation beam irradiating a contour tangentially. The dose at a point such as that labelled Y, is in part due to radiation that comes along the ray line to the point and in part due to radiation that has been scattered to the point from all of the irradiated volume. In this example there is scatter from the right of the point but not from the left. This can be taken into account by integrating a function describing the scattered radiation over the irradiated volume. The kind of differences in calculated doses are seen in Fig. 8 by comparing the solid lines shown with the dashed lines. The latter were produced using a procedure that does take the configuration of the scattering volume into account. The use of scatter-air ratios for this purpose is described at length elsewhere ( 1) and will not be repeated here. The essence of the procedure is indicated schematically in Fig. 9 which shows a beam of photons irradiating a phantom. It is possible to derive differential scatter-air ratios which describe the amount of scattered radiation arriving at points, such as P, or Q, in this diagram, from a rectangular volume element of phantom material as shown. If the scatter is given by A&,(xl, - x&Ax, the total dose at a point (x,,, y,,, z,,) in the phantom is given by
+
s
W&)~f(Xi3, Yd
XB
AWxh - XB). ax,
B
Ax,
B t4) 1
1 220K
Fig. 8. Isodose curves for tangential radiation of a cylindrical surface. Note. From The Physics of Radiology, 4th ed. Charles C Thomas; Springfield, IL; 1983.
where DA is the dose in air at the point of calculation W(x,) is a wedge attenuation factor, f(xB, yB) is a function describing the beam intensity profile in air, T,(zJ is the zero area tissue-air ratio for depth zB and the integral is over the range of the lateral beam coordinate XL. The differential scatter-air ratio is evaluated for depth z’s and the lateral distance xb - xB. As before, there must be a coordinate transformation between the patient coordinates, x1>etc., and the beam coordinate system xB etc. The integral has the appearance of a convolution but the kernel, the differential scatter air ratio ASJAx, is not independent of position and so cannot easily be evaluated by transformation theory as might be hoped. Another name for AS/Ax is the “line spread function” of scattered photons. This two-dimensional integration procedure can be applied equally well to electron beams in which case the primary term in Eq. 4 is deleted and the quantity replacing AS/Ax is derived, either from experiment or by integrating a pencil beam spread function to provide a line spread function.
Dose calculations 0
247
J. R. CUNNINGHAM
(d) Three-dimensional integral methods Mackie et al. (8) has explored the properties of the transport of energy away from the site of interaction of photons using Monte Carlo methods. He has produced “point spread functions” which can in principle be integrated over the three-dimensions of the irradiated volume. The integration would be of the form: D(xp, Yp, ZP)=
s,,s,., Jzp de439
Yi37zL4
x h(xk - x,,, y; - yn, zk - zt,)*dx~dyl,dzi,
Fig. 9. Diagram depicting ditferential scatter-air ratios. The cross-hatched region represents a volume element of an irradiated phantom. Radiation is scattered from it to points such as Ql , which is a distance X from it at a depth d below the surface of the phantom. The amount of scatter is also dependent on the intensity of the radiation incident on the volume element and is, in this case, modified by a wedge filter.
A simple extension of the scatter-air ratio method to the third-dimension may be made by using the sector integration method for a beam of irregular cross-sectional shape as indicated in Fig. 10 but this step carries with it the restriction that the beam intensity be uniform and the phantom surface be flat over the cross-section (7). In order to further extend this calculation procedure the integral in Eq. 4 must be a double integral and one way of doing this is to perform the integration also along the radii of each sector as shown in Fig. 11. The differential scatter-air ratio is in this case in a I&n that is equivalent to a spread function describing a pencil beam of photons. The integration procedure can be applied equally well to electron beams but in this case the pencil beam data may be obtained from electron beam experiments and/or the theory of electron scattering. Integration may either be in polar coordinates as shown in Fig. 11 or in a Cartesian coordinate arrangement.
(5)
where 4 is the photon or electron beam intensity at point x)s, yb, zh in the phantom and h describes the energy transport from there to xB , yB , zB , the point of calculation. Practical methods for evaluating this expression for clinical situations have not yet been worked out. One procedure that may possibly be of help for photons is to split the point spread function, h, into its two physical components; hl representing electron transport and h2 representing photon transport. This would produce two such integral expressions and although the forms would be similar, the domain of the integrations would be quite different and this might facilitate the numerical evaluation. The possibility of using Fourier Transform techniques is also being explored by various workers (9, 10). This method is truly three-dimensional in nature.
‘Cl
124J
Fig. 10. The cross-section of an irregular shaped beam. The scatter is being determined by integration over sectors of circular beams.
ComputerizedMedicalImagingand Graphics
248
May-June/l989, Volume 13,Number3 are examined. The agreement is good but the comparison shows that both the measurement and the calculation must be carried out with care and both must be normalized following the same conventions. Although the example used is for cobalt the test applies also to beams from linear accelerators. The importance of experimental testing, preferably by the creator of the calculation method, cannot be overstated.
h
CORRECTIONS FOR TISSUE INHOMOGENEITIES Methods (a), (b) and (c) above can be applied only with difficulty to the calculation of dose in an inhomogeneous patient and the procedure that has been most commonly adopted has been to calculate first for a shaped but homogeneous patient and subsequently to apply correction factors for the inhomogeneities. The correction factor is defined as
+ I -__dm ____. t I
d
C = 01141
I
Fig. 11. Depiction of scattering volumes that are like pencil beams. For photon beams they can be determined by differentiating scatter-air ratios with respect to radius. Note. From the Physics of Radiology, 3rd ed. Springtield, IL: Charles C Thomas; 1969. Another method which in form is like Eq. 4 has been suggested by Beaudoin (11) and explored further by Wong et al. (12, 13) who have called it the “Delta-volume method.” EXPERIMENTAL
TESTING
All calculation procedures must be tested by comparison with carefully carried out experiments. There are many such sources of data already available such as percent depth dose data and isodose charts. A particularly powerful test is the beam profile as shown in Fig. 12 and 13. Figure 12 shows the wellknown geometrical relations pertaining to beam penumbra and the shape of a simple beam profile (1). Figure 13 shows a direct comparison between measured and calculated profiles for a cobalt unit. Both an open beam and a beam with a wedge filter in place
(6)
where Du is the dose at a point as calculated for the homogeneous phantom and 0, is the dose at the same point but with allowance for inhomogeneities. This has had the merit of separating these procedures and allowing them to be developed and evaluated separately. A number of approximate methods for finding the correction factor for photon beams have been developed over the years. Several of them are described in ICRU Report 24 (14) and by Johns and Cunningham ( 1) and evaluated by Cunningham ( 15) and by Mackie et al. (16). Allowance for tissue-inhomogeneities in calculations for electron beams have most commonly been incorporated directly into the first calculation and again a number of approximate methods have been devised. The literature on this subject is very extensive and useful summaries are given by Sternick (5), Niisslin (2) and by Brahme ( 17). SUMMARY A number of approximation methods for calculating dose and for making corrections for inhomogeneities have been developed. If a target of &5% in accuracy is to be achieved, it is necessary to choose a procedure that will take account of the three-dimensional nature and shape of the irradiated body. Although several such procedures are available, being approximate, they must be expected to have limitations and it must be realized that there will be irradiation conditions that will fool them. Therefore, experimental testing will always be required when new irradiation conditions are employed. The conditions that are not yet well handled for
Dose
calculations
(a)
0
J. R.
249
CUNNINGHAM
(b)
6lU
Fig. 12. Diagram showing the geometrical parameters that shape the profile of the primary component of a radiation beam. Note. From The Physics of Radiology, 4th ed. Charles C Thomas; Springfield, IL, 1983. photon beams are those in which electronic equilibrium is not maintained, that is, where kerma and absorbed dose are not related to each other in a sim-
THERATRON
ple manner. It is likely that the approach using point spread functions and three-dimensional integrations will go a long way toward solving this problem but practical implementation of these ideas is not yet available. For electron beams it is clear that the paths taken by the electrons as they are scattered is very complicated and there is still much to be done to allow for tissue inhomogeneities accurately in a practical way. Acknowledgements-The author wishes to acknowledge discussions with very many colleagues; J. Van Dyk, Alan Rawlinson, J. Jezioranski and others of the Ontario Cancer Institute, Bengt Bjamgard of the Joint Radiotherapy Center in Boston, Garett Halt, of Memorial Hospital in New York, Ma&me And& Dutreix of the Gustav Roussy Institute in Paris and many others, far too many to mention. Much of the subject of this paper also appears as part of ICRU Report 42.
REFERENCES
Fig. 13. Comparison between calculation and measurement of the dose profile for a radiation beam from a cobalt unit. The sloping profile is fo’r the same beam with a wedge filter in place.
Johns, H.E.; Cunningham, J.R. The physics of radiology, 4th ed. Springfield, IL.: Charles C Thomas; 1983. Niisslin, F. Computerized treatment planning in therapy with fast electrons. Medicamundi 24: 112; 1979. Sterling T.D.; Perry, H.; Katz, L. Brit. Automation of radiation treatment planning. IV. Derivation of a mathematical expression for the percent depth dose of cobalt-60 beams and visualization of multiple field distributions. J. Radial. 37544; 1964. Van de Geijn, J.; Chien, I-C.; Pocheng, Ch.; Frederickson, H.A. A Unified 3-D Mode1 for External Beam Dose Distributions. Computers in Radiation Therapy. Proceedings ofthe 7th International Conference on the Use of Computers in Radiotherapy. Tokyo, 1980.
250
Computerized Medical Imaging and Graphics
5. Stemick, E.S. Algoritms for computerized treatment planning. Chapter 4 in practical aspects of electron beam treatment planning. AAPM Monograph No. 2; 1978. 6. Kawachi, K. Calculation of electron dose distribution for radiotherapy treatment planning. Phys. Med. Biol. 20:57 1; 1975. 7. Cunningham, J.R.; Shrivastava, P.N.; Wilkinson, J.M. IRREG-calculation of dose from irregularly shaped radiation beams. Comput. Progr. Biomed. 2: 192; 1972. 8. Mackie, T.R.; Scrimger, J.W.; Battista. J.J. A convolution method of calculating dose for 15 MV x-rays. Med. Phys. 12:188; 1985. 9. Boyer, A.L.; Mok., E.C. A photon dose distribution model employing convolution calculations. Med. Phys. 12: 169: 1985. 10. Mohan, R.; Chiu, C-S. Use of fast fourier transforms in calculating dose distributions for irregular shaped fields for threedimensional treatment planning. Med. Phys. 14:70; 1987. 11. Beaudoin, L. Analytical approach to the solution of dosimetry in heterogenous media. Masters Thesis, University of Toronto, 1968. 12. Wong, J.W.; Slessinger, E.D.; Rosenberger, F.U.; Krippner, K.; Purdy, J.A. The delta volume method for three-dimensional photon dose calculations. Proceedings of the 8th International Conference on the Use of Computers in Radiation Therapy. IEEE Computer Society. 1984, p. 26. 13. Wong. J.W.; Henkelman, R.M. A new approach to CT-pixelbased photon dose calculations in heterogenous media. Med. Phys. 10:199; 1983.
May-June/ 1989, Volume 13, Number 3 14. ICRU. Determination of Absorbed Dose in a Patient Irradiated by Beams of X or Gamma Rays in Radiotherapy Procedures. ICRU Report 24. International Commission on Radiation Units and Measurements. Washington, DC.; 1976. 15. Cunningham, J.R.; Tissue inhomogeneity corrections in photon based treatment planning. In: Orton, C., ed. Progress in Medical Physics, New York: Plenum Publishing; 1982. 16. Ma&e, T.R.; El-Kbatib, E.; Battista, J.J.; Scrimger, J.; Van Dyk, J.; Cunningham, J.R. Lung dose corrections for 6- and 15 MV x-rays. Med. Phys. 12:327; 1985. 17. Brahme, A. Current algorithms for computed electron beam dose planning. Radiotherapy and Oncology 3:347; 1985. About the Author-JOHN ROBERT CUNNINGHAM was born in Regina Saskatchewan and received a B.Eng. degree in 1950 and an M.Sc. in Physics in 195 1 from the University of Saskatchewan and the Ph.D. degree in Physics in 1955 from the University of Toronto. He has worked in the field of Polymer Physics as well as in Radiation Physics but joined the staff of the Ontario Cancer Institute in Toronto in 1958 where he has been in Clinical Radiation Physics to the present. He also holds a position of professor from the University of Toronto. His principle interest is in dosimetry for ionizing radiation and has contributed widely to the development of methods for calculating dose distributions as part of the procedures for using radiation for cancer therapy. He is the author of some 50 publications and is a co-author, with Dr. H. E. Johns ofthe Physics of Radiology.
DOSE CALCULATIONS
1989 Copyright
0
0895-61 I l/89 $3.00 + .OO 1989 Pergamon Press plc
FOR HIGH ENERGY ELECTRON AND PHOTON BEAMS J. R. Cunningham
Physics Division, The Ontario Cancer Institute, Toronto, Canada (Received 3 I August 1988)
Abstract-The procedures of radiation therapy consist of a number of steps, one of which is the calculation of the radiation dosage pattern that will result when a particular arrangement of radiation beams or sources is applied to a patient. The purpose of this calculation is two-fold. One is to predict, as part of the dose planning process, what dose distribution can be achieved with a selected beam arrangement. The second is to record what treatment has been given so that post-treatment analysis can be carried out. Both of these are important. Key Words: Dose calculations, Electronbeams, Photon beams, Tissue inhomogeneities, Dose distributions,Beam models
may be obtained from measurements made in electron beams. The task of a 3-D dose calculation method is to produce this sort of data for beams of arbitrary crosssection, incident from any direction on a surface of any shape and to include allowance for internal inhomogeneities and to do this quickly and to an accuracy better than f5%. Figure 3 shows the problem for photons in a simplified form. A beam of finite size is irradiating an object. The dose is to be calculated at point P. Photons will come along the direct path a and interact at P, releasing energy there. Photons will also travel along path b, interact in the volume dV and send scattered photons along path c to point P, again releasing energy there. The total energy released at point P would then be the sum of the contribution due to primary photons and the total of the scattered photons that come from the whole of the irradiated volume. This can be calculated in a straightforward manner if we include photons that have been scattered no more than once.
INTRODUCTION The problem of calculating radiation dose distributions, as part of radiotherapy planning, is a very old one. The advent of the CT scanner, with its ability to gather anatomical detail, has however forced a reevaluation of the methods used. The task is a mathe-
matical one and is sufficiently complicated that it has no exact solution. The development of methods for treatment planning is therefore the search for ever better methods of making approximate calculations of absorbed dose. In the days before cotmputers, dose distributions were measured in water and isodose charts were produced from these measurements. The isodose charts were then applied to a diagram of the patient’s contour and approximate corrections were made to allow for the shape of the patient’s contour under the assumption that within the contour the patient was water equivalent. Only occasionally were corrections made for the presence of internal nonwater equivalent structures. DOSE CALCULATION
METHODS Kc&-J“=.
0 !f
.Eab
P
Figure 1 shows isodose surfaces for a rectangular 16 X 5 cm beam incide:nt on a flat homogeneous phantom. These data (1) were obtained from measurements made in a water tank. They also show some of the difficulties :in displaying three-dimensional objects. The figure shows isodose distributions in four different planes; I(a) and (b) are in principal planes, that is they contain the central ray, (c) is in a plane that is parallel to (b) but is 2 cm away, and (d) is in a cross-section at a depth of 10 cm. Figure 2 is a “perspective” view of thle same data. Similar data
Here, K is the kerma, or kinetic energy released to the medium and the first term describes the primary component; 4 is the photon fluence, p is the linear attenuation coefficient for these photons, a is the path length through the medium, (p/p) the mass interaction coefficient expresses the probability that the photons will interact at point P and &, is the 241
242
Computerized Medical Imaging and Graphics
May-June/ 1989, Volume 13, Number 3
a)
Fig. 1. Isodose curves in four planes for a beam from a cobalt unit. The planes are the principal planes a cross-section and a plane parallel to one of the principal planes. Note. From The Physics of Radiology, 4th dd. Charles C Thomas; Springfield, IL; 1983.
48K
Fig. 2. A “perspective
view” of isodose surfaces. Note. From The Physics of Radiology, 4th ed. Charles C Thomas; Springfield, IL; 1983.
Dose calculations 0 J. R.
243
CUNNINGHAM
motion and for a complete calculation of the energy absorbed, that is the absorbed dose, the energy transport by these electrons must also be considered. This makes the problem quite unsolvable by analytical means and quite basic assumptions are generally made, the most common being to assume that the electrons that are set into motion deposit their energy locally. This is equivalent to assuming absorbed dose is equal to kerma and Kin Eq. 1 could be interpreted as absorbed dose. Typically, the first term in Eq. 1, the primary, accounts for 60% to 70% of the total dose at point P, first scatter accounts for 20% to 30% of the total and multiple scatter accounts for 10% or more. These proportions are shown in Fig. 4 for cobalt-60 radiation for a range of beam sizes. Although scattered radiation constitutes the minor component of the dose at a point, it is large enough that it must be taken seriously and a three-dimensional calculation method must take its properties into account. EVOLUTION
OF CALCULATION
METHODS
A direct mathematical solution to the calculation of the dose at a point such as P of Fig. 3 would 14OOG Fig. 3. Diagram showing the path for primary and scattered radiation re,aching point P.
average energy deposite’d per interaction. The primary therefore, depends on the path between the source and the point of calculation and the material at the point of calculation. The second term is much more complicated. It describes the energy deposited at point P due to photons that are scattered from every part, dI’, of the irradiated volume. In this term, b is the path length of primary photons to the scattering volume, da(O)/& is the probability of scattering photons of energy hv through angle 0 and C;, (p’/p) and &, are the linear attenuation coefficient, mass interaction coefficient and the average energy deposited from the scattered photons. c is the path length of the scattered photons in this medium. The integration is over the whole of the irradiated volume. This term although already complicated, is solvable on a computer but it accounts for photons scattered only once. Furthermore, the whole equation applies for incident photons of only a single energy. A linear accelerator beam would., for example, consist of a spectrum of photon energies and it would be necessary to integrate both terms over the energy spectrum. In addition, it is necessary to remember that photons are indirectly ionizing particles and that when energy is lost by photons, electrons are set into
2
P 50 ,o
0
-
Total (Ist.+multiple)
2
4
6
6
scatter
IO 12 14 16 16 20 22 24 26 28 30 32 34 36 Beam radius [cm] soln
Fig. 4. Graphs showing the proportion of the dose at a point, 10 cm deep in a phantom irradiated by a cobalt unit, due to primary once scattered radiation and multiply scattered radiation.
244
Computerized Medical Imaging and Graphics
May-June/ 1989, Volume 13, Number 3
involve the Boltzmann transport equation. Transport of energy by both photons and electrons would have to be included and no solution to this equation has been found for beams of finite size even in homogeneous media. Therefore, approximate solutions must be found and the history of the development of calculation methods is a story of the search for better and faster methods of approximation. The situation for electrons is, in some ways easier since they are directly ionizing particles and they do not spread as far as do photons. The Boltzmann equation has been solved for electrons, albeit approximately, and the result is the Gaussian pencil beam model. An excellent review of these solutions has been given by Niisslin (2). The development of approximate solutions for the calculation of dose has, for both photons and electrons, taken place through the evolution of four different categories of model; matrix representation, empirical or semi-empirical beam generating functions, semi-empirical beam models and three-dimensional integral models. (a) Matrix methods The representation of the dose distribution in a radiation beam is shown in Fig. 5 where isodose curves for a beam from a cobalt unit are shown superimposed on a grid of points. The dose value at each of these points can be stored in a computer as a matrix, or table of numbers. The manipulation of such a matrix is indicated in Fig. 6 which shows 3 beams irradiating a patient. Again, a grid of points is shown and the location of these points must be determined in each beam table and the dose from each beam determined by interpolation. Adjustments for contour shape are also made at this stage usually by a procedure known as the “effective SSD method.” The matrix method is fast and is applicable to both photons and electron beams but it is limited to beam configurations for which dose measurements have been made previously and for which such data have been stored. (b) Analytical beam models The oldest of these is due to Sterling et al. (3) and the most fully developed for photon beams is due to Van de Geijn et al. (4). The format followed is that the dose from a beam at a point, such as one of those marked in Fig. 6, is given by a product of two functions. D(x,, Y,, 2,) = WB,
Zrcf)
- &~Xtr,
YB)
(2)
where (x,, , y,, , z,,) are the coordinates of a point within the patient. P(zu, z,,f) is the dose at depth zB along the
on a Cartesian grid of points. The dose at each point can be read off and stored to form a matrix representation of the dose distribution.
Fig. 5. Isodose curves superimposed
beam axis relative to its value at a reference depth zref and is the percentage depth dose. gZ,(xB, ys) is the “off centre ratio” at point xB, yE with respect to the beam axis and depth zB with respect to the contour surface. The latter function describes the dose distribution along a line perpendicular to the beam axis. Van de Geijn has further taken the second function to be itself a product of two functions. g&I?, YB) = gl.z&B)*gz.rB(YB).
(3)
Each of the manipulated quantities, P, g,, and gz, may be determined by choosing analytical functions and fitting the constants associated with them to measured data or alternatively, measured percent depth doses and beam dose profile data may be stored and manipulated by table lookup and interpolation procedures. Allowance for patient surface curvature is made by using ray tracing procedures, such as are indicated in Fig. 7, which shows a beam from a cobalt unit irradiating a curved surface. The dose is to be calculated at point P which is at a depth d’ below the patient surface as measured along a ray line and is at a depth d below the flat surface for which the beam data applies. The dose at P is determined by using d’, the radiological equivalent depth, zB in Eq. 2.
Dose calculations 0 J. R.
245
CUNNINGHAM
/
I \
\ \ .
.
Fig. 6. The dose distribution
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
from three radiation beams determined by calculating the dose at each point from each beam and adding them point by point.
t
F_ T
T h
d
t d’
1/1f
/ -
I20 J Fig. 7. Isodose curves superimposed on a ray line grid. The dose at point P (depth d) below the flat surface is about 55 and is about 65 below the curved surface (depth d’).
Computerized Medical Imaging and Graphics
246
May-June/ 1989, Volume 13, Number 3
This method can be applied to electrons as well as photons and a review of these applications is given by Sternick (5). For photon beams, the function, gZ,(xB, yB), is chosen empirically, for electrons its form has been chosen as a result of solving the transport equation which for a broad beam implies the use of an error function (6). It is interesting to note that Sterling et al. (3) chose this same form for photon beams although his choice was based purely on empirical grounds. The depth dose function, P(zB, zr,r) is taken as an empirical fit to depth dose data for both photon and electron beams. This calculation method is not three-dimensional, since the calculation depends only on the distance and material found along the ray line to the point of calculation. There is no account taken of the three-dimensional character of the irradiated volume.
y
(c) Semi-empirical methods The main example of this type of representation for photon beams is the manipulation of scatter-air ratios and for electron beams is the use of Gaussian shaped pencil beams. A problem that is at least twodimensional in nature is indicated in Fig. 8 which shows a radiation beam irradiating a contour tangentially. The dose at a point such as that labelled Y, is in part due to radiation that comes along the ray line to the point and in part due to radiation that has been scattered to the point from all of the irradiated volume. In this example there is scatter from the right of the point but not from the left. This can be taken into account by integrating a function describing the scattered radiation over the irradiated volume. The kind of differences in calculated doses are seen in Fig. 8 by comparing the solid lines shown with the dashed lines. The latter were produced using a procedure that does take the configuration of the scattering volume into account. The use of scatter-air ratios for this purpose is described at length elsewhere ( 1) and will not be repeated here. The essence of the procedure is indicated schematically in Fig. 9 which shows a beam of photons irradiating a phantom. It is possible to derive differential scatter-air ratios which describe the amount of scattered radiation arriving at points, such as P, or Q, in this diagram, from a rectangular volume element of phantom material as shown. If the scatter is given by A&,(xl, - x&Ax, the total dose at a point (x,,, y,,, z,,) in the phantom is given by
+
s
W&)~f(Xi3, Yd
XB
AWxh - XB). ax,
B
Ax,
B t4) 1
1 220K
Fig. 8. Isodose curves for tangential radiation of a cylindrical surface. Note. From The Physics of Radiology, 4th ed. Charles C Thomas; Springfield, IL; 1983.
where DA is the dose in air at the point of calculation W(x,) is a wedge attenuation factor, f(xB, yB) is a function describing the beam intensity profile in air, T,(zJ is the zero area tissue-air ratio for depth zB and the integral is over the range of the lateral beam coordinate XL. The differential scatter-air ratio is evaluated for depth z’s and the lateral distance xb - xB. As before, there must be a coordinate transformation between the patient coordinates, x1>etc., and the beam coordinate system xB etc. The integral has the appearance of a convolution but the kernel, the differential scatter air ratio ASJAx, is not independent of position and so cannot easily be evaluated by transformation theory as might be hoped. Another name for AS/Ax is the “line spread function” of scattered photons. This two-dimensional integration procedure can be applied equally well to electron beams in which case the primary term in Eq. 4 is deleted and the quantity replacing AS/Ax is derived, either from experiment or by integrating a pencil beam spread function to provide a line spread function.
Dose calculations 0
247
J. R. CUNNINGHAM
(d) Three-dimensional integral methods Mackie et al. (8) has explored the properties of the transport of energy away from the site of interaction of photons using Monte Carlo methods. He has produced “point spread functions” which can in principle be integrated over the three-dimensions of the irradiated volume. The integration would be of the form: D(xp, Yp, ZP)=
s,,s,., Jzp de439
Yi37zL4
x h(xk - x,,, y; - yn, zk - zt,)*dx~dyl,dzi,
Fig. 9. Diagram depicting ditferential scatter-air ratios. The cross-hatched region represents a volume element of an irradiated phantom. Radiation is scattered from it to points such as Ql , which is a distance X from it at a depth d below the surface of the phantom. The amount of scatter is also dependent on the intensity of the radiation incident on the volume element and is, in this case, modified by a wedge filter.
A simple extension of the scatter-air ratio method to the third-dimension may be made by using the sector integration method for a beam of irregular cross-sectional shape as indicated in Fig. 10 but this step carries with it the restriction that the beam intensity be uniform and the phantom surface be flat over the cross-section (7). In order to further extend this calculation procedure the integral in Eq. 4 must be a double integral and one way of doing this is to perform the integration also along the radii of each sector as shown in Fig. 11. The differential scatter-air ratio is in this case in a I&n that is equivalent to a spread function describing a pencil beam of photons. The integration procedure can be applied equally well to electron beams but in this case the pencil beam data may be obtained from electron beam experiments and/or the theory of electron scattering. Integration may either be in polar coordinates as shown in Fig. 11 or in a Cartesian coordinate arrangement.
(5)
where 4 is the photon or electron beam intensity at point x)s, yb, zh in the phantom and h describes the energy transport from there to xB , yB , zB , the point of calculation. Practical methods for evaluating this expression for clinical situations have not yet been worked out. One procedure that may possibly be of help for photons is to split the point spread function, h, into its two physical components; hl representing electron transport and h2 representing photon transport. This would produce two such integral expressions and although the forms would be similar, the domain of the integrations would be quite different and this might facilitate the numerical evaluation. The possibility of using Fourier Transform techniques is also being explored by various workers (9, 10). This method is truly three-dimensional in nature.
‘Cl
124J
Fig. 10. The cross-section of an irregular shaped beam. The scatter is being determined by integration over sectors of circular beams.
ComputerizedMedicalImagingand Graphics
248
May-June/l989, Volume 13,Number3 are examined. The agreement is good but the comparison shows that both the measurement and the calculation must be carried out with care and both must be normalized following the same conventions. Although the example used is for cobalt the test applies also to beams from linear accelerators. The importance of experimental testing, preferably by the creator of the calculation method, cannot be overstated.
h
CORRECTIONS FOR TISSUE INHOMOGENEITIES Methods (a), (b) and (c) above can be applied only with difficulty to the calculation of dose in an inhomogeneous patient and the procedure that has been most commonly adopted has been to calculate first for a shaped but homogeneous patient and subsequently to apply correction factors for the inhomogeneities. The correction factor is defined as
+ I -__dm ____. t I
d
C = 01141
I
Fig. 11. Depiction of scattering volumes that are like pencil beams. For photon beams they can be determined by differentiating scatter-air ratios with respect to radius. Note. From the Physics of Radiology, 3rd ed. Springtield, IL: Charles C Thomas; 1969. Another method which in form is like Eq. 4 has been suggested by Beaudoin (11) and explored further by Wong et al. (12, 13) who have called it the “Delta-volume method.” EXPERIMENTAL
TESTING
All calculation procedures must be tested by comparison with carefully carried out experiments. There are many such sources of data already available such as percent depth dose data and isodose charts. A particularly powerful test is the beam profile as shown in Fig. 12 and 13. Figure 12 shows the wellknown geometrical relations pertaining to beam penumbra and the shape of a simple beam profile (1). Figure 13 shows a direct comparison between measured and calculated profiles for a cobalt unit. Both an open beam and a beam with a wedge filter in place
(6)
where Du is the dose at a point as calculated for the homogeneous phantom and 0, is the dose at the same point but with allowance for inhomogeneities. This has had the merit of separating these procedures and allowing them to be developed and evaluated separately. A number of approximate methods for finding the correction factor for photon beams have been developed over the years. Several of them are described in ICRU Report 24 (14) and by Johns and Cunningham ( 1) and evaluated by Cunningham ( 15) and by Mackie et al. (16). Allowance for tissue-inhomogeneities in calculations for electron beams have most commonly been incorporated directly into the first calculation and again a number of approximate methods have been devised. The literature on this subject is very extensive and useful summaries are given by Sternick (5), Niisslin (2) and by Brahme ( 17). SUMMARY A number of approximation methods for calculating dose and for making corrections for inhomogeneities have been developed. If a target of &5% in accuracy is to be achieved, it is necessary to choose a procedure that will take account of the three-dimensional nature and shape of the irradiated body. Although several such procedures are available, being approximate, they must be expected to have limitations and it must be realized that there will be irradiation conditions that will fool them. Therefore, experimental testing will always be required when new irradiation conditions are employed. The conditions that are not yet well handled for
Dose
calculations
(a)
0
J. R.
249
CUNNINGHAM
(b)
6lU
Fig. 12. Diagram showing the geometrical parameters that shape the profile of the primary component of a radiation beam. Note. From The Physics of Radiology, 4th ed. Charles C Thomas; Springfield, IL, 1983. photon beams are those in which electronic equilibrium is not maintained, that is, where kerma and absorbed dose are not related to each other in a sim-
THERATRON
ple manner. It is likely that the approach using point spread functions and three-dimensional integrations will go a long way toward solving this problem but practical implementation of these ideas is not yet available. For electron beams it is clear that the paths taken by the electrons as they are scattered is very complicated and there is still much to be done to allow for tissue inhomogeneities accurately in a practical way. Acknowledgements-The author wishes to acknowledge discussions with very many colleagues; J. Van Dyk, Alan Rawlinson, J. Jezioranski and others of the Ontario Cancer Institute, Bengt Bjamgard of the Joint Radiotherapy Center in Boston, Garett Halt, of Memorial Hospital in New York, Ma&me And& Dutreix of the Gustav Roussy Institute in Paris and many others, far too many to mention. Much of the subject of this paper also appears as part of ICRU Report 42.
REFERENCES
Fig. 13. Comparison between calculation and measurement of the dose profile for a radiation beam from a cobalt unit. The sloping profile is fo’r the same beam with a wedge filter in place.
Johns, H.E.; Cunningham, J.R. The physics of radiology, 4th ed. Springfield, IL.: Charles C Thomas; 1983. Niisslin, F. Computerized treatment planning in therapy with fast electrons. Medicamundi 24: 112; 1979. Sterling T.D.; Perry, H.; Katz, L. Brit. Automation of radiation treatment planning. IV. Derivation of a mathematical expression for the percent depth dose of cobalt-60 beams and visualization of multiple field distributions. J. Radial. 37544; 1964. Van de Geijn, J.; Chien, I-C.; Pocheng, Ch.; Frederickson, H.A. A Unified 3-D Mode1 for External Beam Dose Distributions. Computers in Radiation Therapy. Proceedings ofthe 7th International Conference on the Use of Computers in Radiotherapy. Tokyo, 1980.
250
Computerized Medical Imaging and Graphics
5. Stemick, E.S. Algoritms for computerized treatment planning. Chapter 4 in practical aspects of electron beam treatment planning. AAPM Monograph No. 2; 1978. 6. Kawachi, K. Calculation of electron dose distribution for radiotherapy treatment planning. Phys. Med. Biol. 20:57 1; 1975. 7. Cunningham, J.R.; Shrivastava, P.N.; Wilkinson, J.M. IRREG-calculation of dose from irregularly shaped radiation beams. Comput. Progr. Biomed. 2: 192; 1972. 8. Mackie, T.R.; Scrimger, J.W.; Battista. J.J. A convolution method of calculating dose for 15 MV x-rays. Med. Phys. 12:188; 1985. 9. Boyer, A.L.; Mok., E.C. A photon dose distribution model employing convolution calculations. Med. Phys. 12: 169: 1985. 10. Mohan, R.; Chiu, C-S. Use of fast fourier transforms in calculating dose distributions for irregular shaped fields for threedimensional treatment planning. Med. Phys. 14:70; 1987. 11. Beaudoin, L. Analytical approach to the solution of dosimetry in heterogenous media. Masters Thesis, University of Toronto, 1968. 12. Wong, J.W.; Slessinger, E.D.; Rosenberger, F.U.; Krippner, K.; Purdy, J.A. The delta volume method for three-dimensional photon dose calculations. Proceedings of the 8th International Conference on the Use of Computers in Radiation Therapy. IEEE Computer Society. 1984, p. 26. 13. Wong. J.W.; Henkelman, R.M. A new approach to CT-pixelbased photon dose calculations in heterogenous media. Med. Phys. 10:199; 1983.
May-June/ 1989, Volume 13, Number 3 14. ICRU. Determination of Absorbed Dose in a Patient Irradiated by Beams of X or Gamma Rays in Radiotherapy Procedures. ICRU Report 24. International Commission on Radiation Units and Measurements. Washington, DC.; 1976. 15. Cunningham, J.R.; Tissue inhomogeneity corrections in photon based treatment planning. In: Orton, C., ed. Progress in Medical Physics, New York: Plenum Publishing; 1982. 16. Ma&e, T.R.; El-Kbatib, E.; Battista, J.J.; Scrimger, J.; Van Dyk, J.; Cunningham, J.R. Lung dose corrections for 6- and 15 MV x-rays. Med. Phys. 12:327; 1985. 17. Brahme, A. Current algorithms for computed electron beam dose planning. Radiotherapy and Oncology 3:347; 1985. About the Author-JOHN ROBERT CUNNINGHAM was born in Regina Saskatchewan and received a B.Eng. degree in 1950 and an M.Sc. in Physics in 195 1 from the University of Saskatchewan and the Ph.D. degree in Physics in 1955 from the University of Toronto. He has worked in the field of Polymer Physics as well as in Radiation Physics but joined the staff of the Ontario Cancer Institute in Toronto in 1958 where he has been in Clinical Radiation Physics to the present. He also holds a position of professor from the University of Toronto. His principle interest is in dosimetry for ionizing radiation and has contributed widely to the development of methods for calculating dose distributions as part of the procedures for using radiation for cancer therapy. He is the author of some 50 publications and is a co-author, with Dr. H. E. Johns ofthe Physics of Radiology.