INVERSE RADIOSURGERY TREATMENT PLANNING THROUGH
DECONVOLUTION AND CONSTRAINED OPTIMIZATION
by
JOSEPH F. HARMON JR.
A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE
UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF THE
REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY
UNIVERSITY OF FLORIDA
1994
ACKNOWLEDGMENTS
Many individuals have contributed to the completion of
this research. First, I am grateful for the guidance and
helpful comments provided by all of my committee members.
Special thanks are due Dr. Frank Bova for his valuable
insight and support while serving as committee chairman. I
wish to also extend my thanks to Dr. Rodhe Mohan and Michael
Lovelock at Memorial Sloan Kettering Cancer Center for
providing me with the differential pencil beam data so
critical to this endeavor. Finally, and most of all, I wish
to thank my wife, Rose, and my children, Ashley, Sarah, and
Taylor, for their patience and never ending support.
TABLE OF CONTENTS
page
ABSTRACT................................................ v
CHAPTERS
1 INTRODUCTION...................................... 1
Overview......................................... 1
Stereotactic Radiosurgery......................... 2
Inverse Treatment Planning........................ 4
Treatment Implementation.......................... 6
2 REVIEW OF LITERATURE.............................. 10
Dosimetry Algorithms............................... 10
Overview............................................ 10
Matrix Models................................... 10
Separation of Primary and Scattered Components.. 12
Monte Carlo Techniques.......................... 13
Energy Deposition Kernels....................... 14
Inverse Radiotherapy Planning Approaches.......... 23
Analytical Methods................................ 24
Iterative Methods ............................... 28
3 DOSIMETRY CALCULATIONS............................ 32
Introduction...................................... 32
Energy Deposition Kernel.......................... 32
TERMA............................................ 41
Convolution....................................... 42
4 BREMSSTRAHLUNG SPECTRA............................ 48
Introduction...................................... 48
Spectrum Model.................................... 49
Optimization Algorithm............................ 56
Results............................................ 59
5 DOSIMETRY MEASUREMENTS............................ 66
Overview....................................... .... 66
Radiochromic Film ................................. 67
TLD Sheets......................................... 77
6 INTENSITY MODULATION .............................. 83
Overview....... ................................... 83
Optimal Intensity Modulation...................... 83
Deconvolution..................................... 83
Backprojection......................... ........ 88
Single Beam Intensity Modulated Dosimetry......... 93
7 RADIOSURGERY PLANNING............................. 97
Overview............... ........................... 97
Coordinate Transformation......................... 98
Multiple Beam Optimization........................ 102
Annealing Parameters.............................. 104
Example Results................................... 107
8 EXTENSION TO GENERAL RADIOTHERAPY................. 135
Overview........................................... 135
Small Field Approximations........................ 136
Large Field Dosimetry verification................ 137
Example Inverse Radiotherapy Case................. 149
9 SUMMARY AND CONCLUSIONS ........................... 153
APPENDIX................................................ 156
REFERENCES.............................................. 226
BIOGRAPHICAL SKETCH.................................... 233
Abstract of Dissertation Presented to the Graduate School of
the University of Florida in Partial Fulfillment of the
Requirements for the Degree of Doctor of Philosophy
INVERSE RADIOSURGERY TREATMENT PLANNING THROUGH
DECONVOLUTION AND CONSTRAINED OPTIMIZATION
By
Joseph F. Harmon Jr.
April 1994
Chairman: Frank J. Bova, Ph.D.
Major Department: Nuclear Engineering Sciences
An inverse treatment planning approach is presented for
radiosurgery which permits conformal dose distributions to
small, irregularly shaped targets. Traditionally, such
targets would require somewhat complex, multiple isocenter
plans, arrived at after much trial and error effort. In
contrast, the approach presented here requires only one
isocenter location and arrives at an optimal plan
automatically. The single isocenter is made possible by
allowing complete beam intensity modulation for multiple
beam's eye views (BEVs) of the target. The automation is due
to incorporation of a deconvolution procedure to determine
the intensity modulation of each BEV and a simulated
annealing algorithm to optimize their respective weighting.
Additionally, this method is fully 3D and accurately models
phantom scatter by incorporating Monte Carlo generated
energy deposition kernels into the dosimetry process. The
extension of this methodology to larger targets for general
radiotherapy cases is also addressed.
CHAPTER 1
INTRODUCTION
Overview
The goal of this research was development of an
automated conformal 3D treatment planning tool for
radiosurgery sized targets incorporating photon beam
intensity modulation and phantom scatter. Development of
such a tool is considered a necessary step to eventually
permit improved treatment of irregularly shaped radiosurgery
sized target lesions by sparing dose to adjacent normal and
critical structure tissue.
Automation of the planning procedure should reduce
excessive trial and error treatment planning times
associated with such radiosurgery cases. Also, due to the
single isocenter approach of this method, the treatment
delivery times should also be significantly reduced.
Finally, the nature of the approach taken will permit a very
flexible beam arrangement thus allowing a wide variety of
mechanical treatment delivery methods.
This chapter will present a brief introduction to
radiosurgery, inverse treatment planning, and potential
methods of conformal treatment delivery. Subsequent chapters
will present details of the specific approach taken and how
it differs from other inverse treatment planning approaches.
Stereotactic Radiosurgery
Stereotactic radiosurgery was originally introduced in
1951 as a method of carefully focusing several small beams
of ionizing radiation to destroy tumors in the brain
(Leksell 1951). Over the years, the idea has been
implemented with a wide range of radiation sources including
200300 kV xrays, multiple 60Co sources, higher energy
photons from linear accelerators, and heavy particles such
as protons (Friedman & Bova 1989a). Of these possibilities,
producing photons with a linac is the least expensive and
thus most widely available option. Recently, this option has
been made even more attractive for radiosurgery due to
development of a system which quickly couples to an existing
linac to improve positional accuracy. This system
incorporates precision machined bearings and a gimbal device
to ensure an isocenter accuracy of 0.2 + 0.1 mm (Friedman &
Bova 1989b).
The attractiveness of radiosurgery is due to the
ability to treat lesions in regions of the brain that are
otherwise inaccessible (Winston & Lutz 1988). Examples of
the clinical applications include treatment of arteriovenous
malformations (AVMs), acoustic schwannomas, malignant brain
neoplasms, and meningiomas (Friedman, Bova & Spiegelmann
1992). The majority of the lesions treated are less than 30
mm in size and are subjected to single fraction doses
ranging from 1500 to 3000 cGy (Friedman, Bova & Spiegelmann
1992).
In conventional radiation therapy, treatment is broken
into many fractions in an attempt to take advantage of
repair of normal tissues. This is accomplished through the
differential biological sensitivity of healthy and tumor
tissue to radiation (Hendrickson & Withers 1991). In
stereotactic radiosurgery however, the radiation is usually
delivered in a single fraction due to the difficulty in
accurately localizing the lesion and positioning the patient
over several treatments. In fact, the developments in
radiosurgery have been slowed partly due to the need for
more precise localization of the target (Leksell 1983).
Because of the high, single fraction doses involved, it is
obviously more important than ever to ensure lesion coverage
while sparing healthy tissue. Thus, a highly accurate method
of dose delivery is needed.
Depending on the specific case, target localization is
accomplished with angiography, CT, MRI or a combination of
these techniques. For example, due to the limitations of
angiography, all AVM radiosurgery planning is supplemented
with CT to ensure adequate target identification
(Spiegelmann, Friedman & Bova 1992). The reference frame for
localization is a precisely built BrownRobertsWells (BRW)
ring, mounted to the patient's head. An additional
localizing frame is then attached to the top of the head
ring. A different localizer is used for each of the possible
diagnostic procedures mentioned earlier. They each, however,
contain reference points which allow one to locate the
lesion in BRW stereotactic coordinates (anteriorposterior,
lateral, and vertical). During treatment planning, a trial
and error method is used to arrive at the number and
direction of multiple beam arcs needed for treatment. Unlike
conventional radiotherapy, radiosurgery treatment planning
is not accomplished with the aid of images from a treatment
simulator. Instead, the geometrical relationship between the
treatment fields and internal anatomy is accomplished solely
from images used from target localization. Also, as part of
treatment planning, an appropriate diameter opening of a
small bore (540 mm) circular collimator is chosen. This
will be used, down stream from the standard rectangular
field defining jaws, to obtain a circular shaped field.
Treatment is then carried out by delivering the radiation in
multiple arcs to the patient making use of the head ring
still attached to the patient's head as illustrated in
figure 11.
Inverse Treatment Planning
For radiation therapy one typically obtains an
approximation to the ideal dose distribution through an
iterative trial and error method. Various combinations of
modified (e.g. wedges, irregular field shapes, compensators,
etc.) and unmodified shaped beam combinations are simulated
on a treatment planning workstation until an acceptable
solution is found. Ideally, one wants as uniform coverage as
possible in the target volume with a rapid fall off in dose
outside the target volume. Standard empirical solutions
Figure 11. A schematic diagram indicating gantry rotation
for a given couch position using the stereotactic
radiosurgery attachment developed at the University of
Florida (Bova 1990a).
exist for the more common cases (e.g. parallel opposed, four
field box, etc.) but those involving multiple targets,
targets of irregular shape or targets near critical
radiosensitive structures, such as the spinal cord for
example, can be difficult to plan for.
An alternate way of approaching the problem is to
calculate the desired "optimal plan" directly from the ideal
dose distribution and is referred to as inverse radiotherapy
treatment planning. An optimal inverse solution will specify
not only the number and direction of beams required but also
the field shape and intensity modulation of each individual
beam comprising the plan. Also, one must limit the solution
to physically real solutions. For example, a solution
resulting in negative fluence values may provide an ideal
dose distribution in theory but is impossible to achieve in
practice. There are several approaches one may take to
arrive at a solution to the inverse planning problem and
these are discussed later in the review of literature.
Treatment Implementation
Regardless of which approach is used to obtain the
treatment planning solution, a method for implementation
must be chosen. Specifically, a method of shaping and
modulating the intensity of each incident field must be
selected. Although hardware concerns are beyond the scope of
this research, it is of interest to briefly present the
techniques currently used or proposed to implement such a
plan.
Other than hand made cerrobend blocks, the most
predominant technique for general radiotherapy field shaping
is the use of a multileaf collimator (Brahme 1988, Moss
1992). In this approach, the incident beam is shaped to
conform to the beam's eye view of the target (lesion)
contour at each orientation of the gantry (beam projection).
A similar approach uses multiple independent vanes (Leavitt
et al. 1991) to shape the borders of the field. However,
these methods are currently limited to treating convex
shaped targets for general radiotherapy since intensity
modulation across the beam is not included (Brahme 1988). It
also excludes the possibility of optimum treatment when an
organ or region at risk is included within the target
volume. Regardless of these limitations it is believed that
field shaping alone would offer an improvement over current
nonconformal techniques used in radiosurgery (Sixel,
Podgorsak & Souhami 1993). In fact, studies have shown that
approximately half of all stereotactic radiosurgery cases
would receive a smaller radiation dose to normal brain
tissue if conformal field shaping treatment options were
available (Leavitt et al. 1991). Recently, the potential
usefulness of multileaf collimators has been extended to
include treatment of concave shaped targets and avoidance of
enclosed critical structures (Galvin, Chen & Smith 1993). In
this approach, the multileaf collimator is used to define a
series of field shapes that are superimposed for a fixed
beam's eye view (BEV) to produce the desired modulated
fluence pattern instead of simply shaping the fieldss.
A second approach is to construct a customized
intensity modulation device compensatorr) to shape and
attenuate the field as needed for each incident beam.
Compensators are used in select general radiotherapy cases
but usually take considerable time to create. Thus, they are
only used when a small number of beams are present. This
technique, as well as the multileaf approach, is not used
for linac based stereotactic radiosurgery since the
treatments are currently achieved by multiple continuous
arcs. However, these techniques may prove useful for a
treatment plan based on a reasonable number of highly
modulated incident beams as proposed by this research.
Another possible approach is to scan a pencil beam of
photons across the target volume (Ngfstadius, Brahme &
Nordell 1984). In principle this seems attractive but is
limited to facilities equipped with such specialized
equipment. Also, the intrinsic bremsstrahlung process and
multiple scattering within the photon producing target
result in a photon beam too broad for radiosurgery
applications.
A fourth approach is irradiation of the target with a
slit field incorporating intensity modulation along the
length of the slit (Carol et al. 1993; Mackie et al. 1993).
The common theme in this approach is the irradiation of
multiple, contiguous slices with a temporally modulated
slit. In order to treat a 3D region, both parallel slices
(PeacockTM) and spiral slices (tomotherapy) have been
proposed. Thus, this technique is analogous to image
reconstruction in computed tomography where measured
transmission profiles are backprojected to obtain a 2D
density distribution. However, the difference here is that a
"measured" intensity modulation function (IMF) is "back
projected" to obtain a dose distribution. Possible problems
with this approach include the extreme geometrical accuracy
needed to avoid "hot and cold" spots in the dose
distribution between adjacent slices.
Regardless of the specific approach, implementing an
inverse solution involves complete beam compensation (field
shaping and intensity modulation across the field) to
achieve nonuniform fields related to a desired dose
distribution instead of achieving uniform fields at depth
through incorporation of patient contours (Goitein 1990).
The steep, often irregular dose profile encountered in
radiosurgery represents a difficult inverse problem to
solve. To date, no treatment planning methodology has been
applied specifically for such a solution for stereotactic
radiosurgery that includes complete beam compensation. Thus,
the following chapters will address this issue.
CHAPTER 2
REVIEW OF LITERATURE
Dosimetry Algorithms
Overview
Regardless of how one approaches the inverse
radiotherapy planning problem, an accurate algorithm for
calculating dose distributions must be incorporated. There
are several models to choose from although they do not all
work equally well for a particular situation. The following
sections will review the models available and point out
their weaknesses and strengths.
Matrix Models
A simple and straightforward method of dose calculation
is to measure and store isodose curves or beam profiles in
an array. For completeness, this would be accomplished over
the range of square field sizes and depths one expects to
use. For rectangular field sizes, an equivalent square field
can be used. Multidimensional interpolation is used to
obtain data at positions not measured. Thus, a significant
number of measurements may be required if one wants to
reduce interpolation errors. For common field sizes and
shapes, one can obtain off axis ratios, f(d,x) and g(d,y),
that can be multiplied by central axis depth dose data,
%dd(d), to obtain the dose, D(d,x,y) at any point (Horton
1987):
D(d,x,y) = %dd(d)f(d,x)g(d,y)
This model is commonly referred to as the MilanBentley
algorithm and has been applied to stereotactic radiosurgery
through measurements of small circular fields (Podgorsak et
al. 1988). In fact, the dose model for stereotactic
radiosurgery at the University of Florida uses a similar
approach (Suh 1990). In this case, the dose at position p,
Dp, for a circular field, is found from an offaxisratio
(OAR), tissue maximum ratio (TMR), and a relative output
factor (ROF) obtained from measured data such that:
Dp(C,STD,d,r)=Dref ROF(C) TMR(w,d) OAR(C,STD,d,r) (SAD/STD)2
where
C = field size at SAD (source to axis distance)
STD = source to target distance
d = depth at point p
r = offaxis distance
w = field size at point p
Dref = reference dose
A closely related technique is that of storing the
measured dose information in the form of decrement lines
instead of a grid of points or offaxis ratios for all
depths of interest. These lines, usually expressed as a
second order polynomial, represent the distance from the
central axis to a given offcenter ratio. The offaxis data
can then be represented by a small number of constants for
each field size and can be combined with central axis
percent depth dose data to calculate dose at any given
point.
These matrix type algorithms are fast and simple but
require a large number of measurements and are not suitable
for irregular field shapes or intensity modulated beams
since each variation would require a complete set of
measurements.
Separation of primary and scattered components
A commonly used technique for dose calculation is
separation of the beam into primary and scattered
components. In this model, the scattered photons are those
defined to have undergone at least one interaction in the
irradiated object before they contribute, in a subsequent
interaction, to the dose at the point of interest (Bjarngard
& Petti 1988). This technique, proposed by Clarkson and
elaborated by Cunningham (Horton 1987) defines the primary
dose as
Dprim = Da(d)f(x,y)TAR(d,0)
where
Da(d) = dose in air for a given field size at depth, d.
TAR(d,0) = tissue air ratio for zero area field size at
depth, d.
f(x,y) = a function dependent on penumbra collimatorr
transmission, source size) and type of filter.
The function f(x,y) is found by comparing measured and
calculated data in an iterative fashion.
The scattered component of dose, unlike the primary
component, depends on field size and shape. It can be
calculated by dividing the field into wedge shaped sectors
and summing the scatter air ratio (SAR) values for each
sector to obtain the average SAR for the field (Johns &
Cunningham 1983):
D.,= A2C .SARi(d,r)
Thus, the total dose is the sum of that from primary and
scatter:
Dtot = Dprim + Dscatt
Although this method is conceptually appealing, it does have
limitations. It has been pointed out (Mohan & Chui 1986)
that this method assumes absorbed dose is equal to KERMA.
This is a good approximation in many cases but not for high
photon energies or near beam edges where electron
equilibrium does not exist. Also, the concept of zero area
TAR is flawed since the primary dose actually becomes zero
as the field size becomes zero (Nizin & Chang 1991; Mohan &
Chui 1986).
Monte Carlo Techniques
Analytically, it is virtually impossible to accurately
include the effects of electron transport such that one can
compute absorbed dose from KERMA (Bjarngard & Cunningham
1986). However, one can simulate electron interactions with
Monte Carlo techniques. In fact, Monte Carlo techniques
provide the most complete description of dose distributions
in both homogeneous and nonhomogeneous materials (Horton
1987). Unfortunately, the computer time required to follow
the hundreds of thousands, or even millions of photons,
required to ensure high accuracy is currently prohibitive.
For this reason, Monte Carlo techniques currently find
little application in routine treatment planning.
Energy Deposition Kernels
Although Monte Carlo techniques are too time consuming
to use on a routine basis, they can be used to resolve the
dose vs KERMA dilemma. Photon beam algorithms (Mackie et al.
1988, Boyer & Mok 1985, Mohan, Chui & Lidofsky 1986,
Ahnesjo, Saxner & Trepp 1992) have recently been introduced
that calculate dose by combining the distribution of photons
originating from the linac head with a kernel representing
the interaction of these photons in the phantom.
Specifically, the kernel describes the complete history of a
given photon from the point of initial interaction to the
point of complete energy deposition of the resulting
electrons. These algorithms are attractive due to their
sound physical basis, accuracy, and versatility. Once a
kernel has been established, intensity modulated beam
dosimetry can be accomplished by simply determining the
effect on the primary photons.
The kernels can be divided into three types. The most
basic type is the single voxel interaction kernel. It
represents the fraction of energy deposited in a voxel at
point r that is released by primary photons interacting in a
xray source
phantom surface
xray
kernel axis
deposition
site
beam axis
Figure 21 Relationship between the interaction site (T')and
dose deposition site (?) using energy deposition kernels
(Sharpe & Battista 1993).
Primary
photon
energy
Primary First
charged scatter
particle photon
kinetic energy
energy
Primary First Second
energy scatter scatter
deposited charged photon
particle energy
kinetic
energy
First Second Multiple
scatter scatter scatter
energy charged photon
deposited particle energy
kinetic
energy
Second Multiple
scatter scatter
energy charged
deposited particle
kinetic
energy
Bremsstrahlung
plus annihilation
photon energy Multiple
scatter
energy
Bremsstrahlung deposited
plus annihilation
charged particle
kinetic energy
Bremsstrahlung
plus annihilation
energy deposited
Figure 22. Energy deposition partitioning used in kernel
development (Mackie et al. 1988).
voxel at point r' (see figure 21). Since it is
experimentally impossible to make photons interact in a
single voxel, the kernels are usually obtained through Monte
Carlo modeling (Mackie et al. 1988, Mohan, Chui & Lidofsky
1986) or analytical approximations (Boyer & Mok 1985, 1986).
However, the analytical models only follow scattered photons
and not electrons thus assuming KERMA is equal to absorbed
dose. Also, the analytical approach is limited to lower
energies (e.g. 60Co) since they incorporate KleinNishina
cross sections for the Compton interaction but neglect other
effects. Higher energy effects are more effectively modeled
with the Monte Carlo method. For example, Mackie et al.,
1988, has calculated energy deposition kernels using the
Electron Gamma Shower (EGS) Monte Carlo code and includes
the effects shown in figure 22. An example of such a kernel
is shown in figure 23. This type of kernel is different
from others mentioned in the literature. The pencil beam
kernel terminology often encountered represents a dose
distribution about a very narrow beam of radiation
interacting with the phantom. This type of kernel can be
obtained from Monte Carlo measurements (Ahnesjo, Saxner &
Trepp 1992), depth dose measurements (Chui & Mohan 1988), or
by convolving a single voxel interaction kernel with an
attenuation term along a ray line through the phantom
(Brahme 1988). An example of a pencil beam kernel is shown
in figure 24. If several pencil beams interact at a voxel
(point), a dose distribution referred to as a convergent
900
1200 (00
1013
160 0 0 10
Radius (g tm)
Figure 23. Isoline format of a total energy deposition
kernel for 1.25 MeV photons in water. The units are cGy MeV~
1 photon1 (Mackie et al. 1988).
point irradiation kernel is found. If the irradiation of the
voxel is from a 2t rotation in a plane or from all
directions (4i) in threedimensions, a very smooth kernel is
formed representing the steepest dose gradient possible
about the voxel. An example of a kernel for the 2D rotation
case is shown in figure 25. In the next section a method
for solving the inverse problem in radiation therapy using
these convergent point irradiation kernels will be
1.25 MeV
E
o 10
N ;
15 i0
20
10 5 0 5 10
r/cm
Figure 24. Isoline format for a 1.25 MeV pencil beam
kernel. Units are fractional energy deposited per unit
volume (cm3) (Eklof, Ahnesjo & Brahme 1990).
10
UO'
10
5
10
10 5 0 5 10
z/cm
Figure 25. Isoline format of a convergent point irradiation
kernel from a 3600 rotation of a pencil beam in a single
plane perpendicular to this page (Eklof, Ahnesjo & Brahme
1990).
presented. However, as will be shown, computational problems
arise when this approach is applied to the stereotactic
radiosurgery case.
The algorithm of interest for my research makes use of
the single voxel interaction kernel. The relationship for
the absorbed dose at r, D(r), can be expressed as
D(r) = fT(r')k(r r')d'r'
where T(r') represents the distribution of total energy
released to matter (TERMA) in each voxel of the phantom,
k(rr') is the kernel representing the deposition of energy
about a single interaction voxel, and d3r' represents
integration over all interaction sites in three dimensions.
T(r') is found from the primary energy fluence distribution,
'(r'), and the mass attenuation coefficient p/p. For a
homogeneous phantom, the TERMA can be expressed as
T(r') = p/p I(r')
The first step in dose calculation is determining the
primary energy fluence at each point in the phantom. One
method of accomplishing this is to first obtain the inair
photon fluence distribution for the open field. Exponential
attenuation is then applied along ray lines considering beam
modification and patient anatomy (Mohan, Chui & Lidofsky
1986). The TERMA is then found as described above and
convolved with the energy deposition kernel. If the kernel
is invariant, the convolution can be achieved much more
quickly by employing discrete Fourier transforms (DFT).
Situations resulting in noninvariant kernels include phantom
heterogenieties, significant beam divergence, and beam
hardening as a function of offaxis distance. However, if
one is modeling a small diameter beam in a homogeneous
media, the kernel can easily be treated as invariant. The
assumption of an invariant kernel for larger fields will be
discussed in chapter 8.
As mentioned, a convenient feature of this algorithm is
its ability to incorporate a wide variety of complicating
factors while retaining its simple physical foundation. Much
of its flexibility lies in the reliance on primary energy
fluence data which can easily be modified. In fact, even
first order corrections for heterogenieties in a phantom can
be accounted for by considering their attenuation effect on
the energy fluence and then multiplying this fluence by the
local attenuation coefficient to determine the number of
primary photon collisions in the scattering voxel (Mohan,
Chui & Lidofsky 1986). Further corrections may be included
from density scaling along the path between the scattering
voxel and the dose computation point. Another strength of
the algorithm is the ability to model the buildup region and
field edges where electron nonequilibrium exist. This is
especially critical when trying to accurately model the
steep dose gradient from multiple overlapping beams. It has
been shown that for a 15 MV, 10X10 cm field, the lateral
transport of secondary electrons alone contribute as much to
penumbra as all other factors (geometry, photons scattered
from beam defining system, and photons scattered in phantom)
combined (Mohan, Chui & Lidofsky 1986).
As a result of the features just mentioned, this model
has a promising future and has even been used for
traditional stereotactic treatment planning (Kusbad et al.
1990). However, it has been not been used for inverse
treatment planning in stereotactic radiosurgery.
Inverse Radiotherapy Planning Approaches
Inverse radiotherapy planning is a specific case of the
more general inverse problem in radiation transport which is
applicable to nuclear engineering, oil well logging,
atmospheric studies, medical imaging, etc. (McCormick 1992).
While direct transport problems are concerned with
estimating the particle distribution within and emerging
from a prescribed medium, inverse transport problems are
concerned with estimating the properties of an obscured
medium with which particles have interacted and then emerged
from. In the case of radiation therapy, for example, a
direct transport problem would be the calculation of an
absorbed dose distribution in a phantom from a photon beam
that has been modified by a compensating filter. An inverse
transport problem example would be estimation of the
compensating filter properties given an absorbed dose
distribution in a phantom.
Although there are several types of radiative transfer
inverse problems, there are two general approaches for
solving them: explicit and implicit (McCormick 1992).
Explicit approaches are analytical in nature and are
typically applied to problems involving significant
simplifications (e.g. planeparallel illumination, plane
geometry, one spatial variable, etc.) (McCormick 1992).
Regardless of these simplifications, the explicit solutions
are suitable for a wide variety of realistic problems and
may also serve as initial estimates for implicit methods.
Implicit approaches are iterative in nature. Each iteration
is simply a solution to the direct radiation transfer
problem where the number of iterations performed depends on
the desired level of accuracy.
In the following two sections, a summary of methods
used to approach the inverse radiation therapy problem is
presented. The first section will present explicit
(analytical) methods while the second will present implicit
iterativee) methods.
Analytical Methods
An analytical approach to inverse treatment planning is
possible if simplifications are made. Early techniques
assumed simple attenuation of a parallel incident photon
beam and calculated the intensity modulation function (IMF)
for 2D circularly (Brahme, Roos & Lax 1982) and
noncircularly (Cormack 1987; Cormack & Cormack 1987)
symmetric dose distributions inside a cylindrical phantom,
for arc therapy, with the axis of rotation coinciding with
the origin of the phantom. A similar approach was later used
to obtain the IMF for a 2D circularly symmetric dose
distribution inside an arbitrarily shaped phantom with
rotation about any axis (Barth 1990). Having accomplished
this, one can consider 2D nonsymmetric dose distributions
by considering them as composed of adjacent,
nonintersecting, radially symmetric dose distributions as
shown in figure 26. Analytical solutions of this type are
quickly obtained but do not account for photon scatter or
electron transport (Kooy & Barth 1990). The solutions also
assume rotational therapy thus excluding the option of a
finite number of projections.
Figure 26. A nonsymmetric dose distribution composed of
adjacent radially symmetric dose distributions located
inside an arbitrarily shaped phantom (Barth 1990).
Another, quite different approach to obtaining IMFs for
given beam projections was proposed by Mackie, Scrimger and
Battista (1985). It was proposed that by deconvolving a
Monte Carlo generated single voxel energy deposition kernel
from an ideal dose distribution in the phantom, the
resulting primary energy fluence distribution could be
obtained. This idea makes use of the kernels (discussed
previously in Dosimetry Algorithms) which include effects of
secondaryparticle transport.
A related idea was actually investigated (Brahme 1988)
with encouraging results. The problem was approached by
dividing the target into many small voxels. For each voxel,
the dose distribution resulting from a pencil beam kernel,
rotated 3600 while intersecting the voxel, was found. This
rotated pencil beam dose distribution represents the
convergent point irradiation kernel discussed earlier. As
long as the symmetric, 3600 irradiation by pencil beams
about each voxel in the target is present, a single,
invariant kernel can be used to calculate dose throughout
the target volume and represents a Fredholm integral
equation of the convolution type. The next step, thus,
consisted of deconvolving the rotated pencil beam kernel
from the desired dose distribution to obtain a 2D function
representing the number of convergent point irradiation
kernels needed in each target voxel. This density function
was then projected, as in CT, to obtain a 1D modulation
function for all incident beam directions. An example of
this method, extended to include 3D cases (Boyer, Desobry &
Wells 1991), is shown in figure 27. As can be seen, this
". Nh 1 1 1 K ri
Figure 27. Resultant dose distribution for a torus shaped
target (Boyer, Desobry & Wells 1991). The isodose lines
shown at the bottom of the page represent a horizontal slice
through the surface rendering above.
method produces good results which include higher order
effects such as electron transport etc. and has been
proposed as a treatment planning tool for tomotherapy
(Holmes & Mackie 1994). However, situations involving
nonsymmetric irradiation become difficult to solve since the
rotated pencil beam kernel is no longer invariant. More
realistic situations involving nonsymmetric
irradiation can be crudely approximated by deconvolving
symmetric irradiation kernels but projecting the resulting
density function for only those incident beam projections
chosen. Methods of choosing appropriate beam projections are
discussed in Soderstrom and Brahme (1992). Otherwise, more
cumbersome numerical techniques are required to solve the
Fredholm integral equation of the first kind resulting from
non invariant point irradiation kernels (Lind & Brahme
1992).
Iterative Methods
One implicit method of obtaining an IMF is through
algorithms similar to those used in computed tomography.
Recall, that in CT, a 1D transmission function is measured
and represents the projection of a 2D density distribution
of that slice of the patient for a given projection angle.
Analogously, a 1D IMF for a given projection angle can be
determined by first projecting the ideal 2D dose
distribution. This comparison is illustrated in figure 28.
Due to much similarity, it has been proposed (Bortfeld
et al. 1990) that CT algorithms can be transferred to this
new application. However, one significant difference between
CT and radiation therapy "backprojection" lies in the
distribution of values along each ray line. In CT the
cumulative transmission value obtained from attenuation
through the patient slice is evenly redistributed along the
same ray line during backprojection. Of course, this is not
possible when backprojecting a physical beam due to
radiation interactions (attenuation and scattering).
Computer Tomography p ~, econ Conlormatlon Therapy
Target
Ra(dion Rad
Figure 28. Sketch comparing computed tomography and
intensity modulated radiotherapy (Bortfeld et al. 1990).
Unlike the analytical approaches just mentioned, a
fixed number of projection angles can De determined thus
giving more flexibility than rotational therapy. Of the
algorithms used in CT, the filtered backprojection followed
by iterative optimization has been investigated (Bortfeld et
al. 1990). Filtered backprojection provides a good first
guess which is then used as input to an iterative algorithm.
The results seem to indicate that the expense of using more
than seven or nine beams was not justified with improved
results. Unfortunately, no consideration was given to
inhomogenieties and scattering was only approximately
included (Bortfeld et al. 1990).
Another type of implicit method relies on direct
transport and numerical optimization alone. Of the many
treatment plan optimization studies accomplished, only a
select few have truly addressed the inverse radiotherapy
problem by including intensity modulation of the incident
beam (Soderstrom & Brahme 1993; Webb 1989, 1991a, 1991b).
These studies approach the problem by dividing each incident
macroscopic beam into multiple elemental beam segments.
Multidimensional optimization routines are then used to vary
the intensity of each elemental beam while keeping score of
the cumulative absorbed dose distribution in the phantom.
The number of variables involved in such an optimization
problem can be quite large even for 2D cases. For example,
allowing ten incident beams with each beam being composed of
ten (2D case) or 100 (3D case) elemental beams would
result in 100 and 1000 variables respectively.
It has been shown that the solution space to such an
optimization problem is filled with local minima traps
(Soderstrom & Brahme 1993). The most effective approach to
such a problem is using a simulated annealing algorithm
which will be discussed in chapter 4. This approach, used by
Webb, is the most comprehensive and has been extended to 3D
geometries although the inclusion of scatter is currently
limited to the 2D case. An advantage of this type of
implicit approach is the degree of complexity allowed in the
planning process. However, the price paid is the extensive
computing resources needed.
A compromise between the two inverse problem approaches
would be an analytic 3D deconvolution approach to arrive at
any individual BEV intensity modulation function followed by
an iterative beam weight optimization. In this case, the
number of optimization variables would be limited to the
number of incident macroscopic beams. Such an approach is
the topic of this report and has significant advantages.
First of all, the deconvolution of single voxel interaction
kernels and subsequent backprojection permit solutions for
any desired number and direction of incident beams. Thus,
significant planning flexibility is possible. Also, by using
the energy deposition kernels, photon and electron transport
are included resulting in realistic dosimetry based on
fundamental principles. A third advantage results from
optimizing macroscopic beam weights instead of individual
beam elements. Due to the much fewer number of variables, a
significant increase in speed can be expected. The details
of this approach are discussed in the next several chapters
and will be followed be several example test cases.
CHAPTER 3
DOSIMETRY CALCULATIONS
Introduction
This chapter describes the 3D dosimetry algorithm used
for the inverse treatment planning approach of this
research. For clarity, the discussion will be limited to
unmodulated, monoenergetic photon beams and flat,
homogeneous water equivalent phantoms. Applications
involving intensity modulation, polyenergetic beams, and
irregular surfaces will be considered in subsequent
chapters.
The dosimetry algorithm used was chosen for its sound
physical basis, accuracy, and ease of incorporation into a
3D inverse radiotherapy protocol. As described earlier, the
absorbed dose is determined through convolution of an energy
deposition kernel and a TERMA distribution resulting from
the primary beam. I will first discuss the characteristics
of the energy deposition kernel and TERMA. I will then
illustrate how they are incorporated into a dosimetry
algorithm through convolution.
Energy Deposition Kernels
The kernels used for my research are differential
pencil beams (DPBs) developed at Memorial Sloan Kettering
Cancer Center in New York (Mohan, Chui & Lidofsky L, 1986).
These kernels were specifically created for photon based
dosimetry in a tissue equivalent medium and are available to
anyone interested in their applications. These kernels are
equivalent to those generated by Mackie et al. (1988) but
are provided in ASCII format for machine compatibility. The
DPBs, as received, have units of
Energy deposited in calculation voxel (MeV) per cm3
10" Monoenergetic photons at interaction point
associated with each and every annulus shaped voxel within a
cylindrical water phantom. The initial interaction for every
incident photon was made to occur at the center (interaction
point) of the cylinder and the resulting location of energy
deposition from secondary particles was determined via Monte
Carlo methods. The bins were scaled in a logarithmic fashion
such that the radial and height dimensions of each annulus
is 10% larger than the previous one. This arrangement allows
fine resolution near the interaction point and a course
resolution further away. The DPBs, as received, were also
divided into two components. A short range component
included energy deposition from electrons produced in the
initial photon interaction at the interaction point (due to
Compton electrons, pair production electrons, and
photoelectrons). A long range component included energy
deposition resulting from subsequent interactions of Compton
photons, Bremsstrahlung, and positron annihilation. This
arrangement arose from a desire to decrease the size of the
computation array when performing dose computations. With
this arrangement, a small, fine grid for the rapidly varying
short range component can be overlaid on a larger, coarse
grid for the slowly varying long range component.
Upon receiving the files, an algorithm was first
developed to determine the total and component energies for
each of the 22 incident photon energies included in the
library. This program, DPBPROC, is included in the appendix
and a summary of the entire DPB library is shown in figures
31 through 33. Note that, at higher energies, the summed
energy total of both components is slightly less than that
of the incident photon. This, of course, is due to the
finite size of the cylindrical phantom which was
approximately 2.5 meters in both height and diameter. Using
DPBPROC, the components were then summed to produce a
single, complete DPB for each photon energy in the library.
Next, a second algorithm, DPBVOX, was developed to
convert the summed DPBs into a linearly scaled Cartesian
coordinate format needed for the convolution process. This
was accomplished by first obtaining a volume weighted
average of all logarithmic scaled voxels lying within the
limits of a given linearly scaled voxel. This provided an
average energy deposition density (MeV/cm3 per 106 incident
photons) for each linearly scaled annulus shaped voxel. The
linearly scaled annulus voxel format was then converted to
that of equal sized cubes along the radial (x or y axis) and
longitudinal (z axis) directions by multiplying the energy
density at each position by the appropriate cube volume.
Finally, the DPBs were normalized by dividing the deposited
1 
S0.9

" 0.8 
C
S0.7
C

S0.6
E
8 0.5
6 0.4
z
0.3
" 0.2
5 0.1 i
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Actual Photon Energy (MeV)
l Shori Long I
Figure 31. Differential pencil beam components, 0.11 MeV.
36
10
9
8
7
6
E
05
4
3
2
2 3 4 5 6 10
Actual Photon Energy (MeV)
I Short Long
Figure 32. Differential pencil beam components, 210 MeV.
20 25 30 40
Actuol Photon Energy (MeV)
SShort M Long
Figure 33. Differential pencil beam components, 1550 MeV.
energy assigned to each voxel by the energy per incident
photon (in MeV) and number of incident photons (106). Thus,
the DPBs may now be expressed in units of
Energy deposited in calculation voxel
Energy released at interaction point
representing the fraction of incident photon energy that
ends up being deposited in a particular calculation voxel.
Having characterized the DPBs as needed in the x (or y)z
plane, a spreadsheet template was used to produce a 3D
version. One simply enters the linearly scaled values for
each radial (x or y axis) and depth (z axis). The
spreadsheet then makes use of radial symmetry and assigns
values to the remainder of the xy plane for every depth.
The radial and depth characteristics of the 2 MeV DPB in
Cartesian format is illustrated in figure 34. In this
example, the grid spacing is 2 mm. Other spacing can be
chosen as needed. However, the uniformity of spacing cannot
be relaxed due to restrictions imposed by the FFT
convolution procedure. A comparison of the linearly scaled
Cartesian coordinate format DPB data and the original
logrithmically scaled cylindrical coordinate format is shown
in figure 35. Note that the linearly scaled Cartesian
coordinate curves are much smoother than the raw Monte Carlo
data thus significantly reducing any high frequency noise
during convolution/deconvolution.
6 8
Radial Distance (mm)
 1 mm Depth 5 mm Depth 9 mm Depth
 13 mm Depth 17 mm Depth 21 mm Depth
Figure 34. Uniformly scaled (2 MeV, 2x2x2 mm3 voxels)
differential pencil beam data for selected depths.
x
0
c 0.1
0
o
0.01
0
a
. 0.001
S0.0001
C
o 1E05
LI
10
.E 2
S 0.1 i
a1
0
0.01
S0.001
0110.0001
C
1E05
1E05
0 5 10 15 20 25
Radial distance (mm)
30 35 40
o Row 5.23 mm + Avg 5 mm W Row 14.91 m x Avg 15 mm
Figure 35. A comparison of the linearly scaled Cartesian
coordinate (avg) and logrithmically scaled cylindrical
coordinate (raw) format for the 2 MeV, 2x2x2 mm3 case at
depths of 5 and 15 mm.
0
0
D
ce h : x x .
TERMA
The energy released at any location within a phantom
can be described in terms of TERMA. As mentioned earlier,
the TERMA at a given point represents the total energy
released to matter (water phantom) for a given number of
incident photons at that point. It can be found by applying
inverse square and attenuation to a calibrated source of
photons:
r ,MeV iu c I f
TERMA= P cm I e 
L 'cm'JJL g J r
where
Yr = reference incident energy fluence for Monoenergetic
photons at a specified distance from the source.
 = mass attenuation coefficient of phantom for
P
monoenergetic photon energy.
s = distance from source to reference distance.
r = distance from source to TERMA calculation point.
ssd = distance from source to surface of phantom.
The TERMA at each interaction point is expressed in
terms of energy released, Erel, per mass, m, and per
incident photon energy fluence (e.g. [pJ/kg]/[MeV/cm2]). Note
that the energy per mass is not referred to as absorbed dose
in this case since the energy is that released at the
interaction point but not all deposited, Edep, at that
location. By choosing units as described above, it is seen
that the product of TERMA and a DPB in a convolution
operation yields absorbed dose per reference energy fluence
(D/wr):
(E JmJ (EEmJ/ ,
Thus, the absorbed dose at any point in the phantom can be
found if the reference energy fluence is known for the
incident beam. However, determination of energy fluence for
an intense source of high energy photons is not easily
accomplished. I have devoted all of chapter five to this
topic. For the purpose of this chapter, I will simply
present all output in relative units.
Convolution
The convolution approach used for dose calculation is
very similar to that often used in image processing. Instead
of convolving an image with, for example, a smoothing or
edge enhancement kernel, a TERMA distribution is convolved
with a DPB. The FFT is used to speed up the convolution
process. As previously mentioned, this is especially useful
when working with large 3D arrays. Additionally, once the
procedure for FFT convolution is established, deconvolution
becomes an option. In fact, the deconvolution process is
used to determine the desired intensity modulation for a
given incident beam and will be described in chapter 6.
The FFT algorithm chosen for this work is a multi
dimensional routine (Press et al. 1992). As input, one
provides the number of dimensions involved, a dimension
vector identifying the length of the array in each
dimension, the array of data, and an integer (1) to
indicate FFT or FFT1. The array of data is actually
accessed as a lD array with alternating real and imaginary
parts. Thus the dimension vector is used to interpret the
arrangement of data within the lD array. Finally, the
length of the data array in each dimension must be a power
of two.
The FFT routine was incorporated into the dosimetry
algorithm as illustrated in figure 36. The DPB is first
read from a file and the FFT is computed. Next, the
TERMA is computed and its FFT found. The convolution is then
accomplished via the product of the two FFTs in frequency
space and then performing an inverse FFT (Starkschall 1988).
It has been shown (Boyer 1984) that the factor by which a
Fourier transform of two functions decreases the calculation
time over traditional convolution is equal to
N3/(3 Log2 N)
where N is the number of calculation points in each
dimension of a cube. Thus, over a 50x50x50 voxel region, a
time savings of 7,400 would be realized. In general, there
is a linear relationship between calculation time and Nlog2N
(Boyer, Wackwitz & Mok 1988). A comparison of the relative
usefulness of Fourier vs direct processing convolution is
shown in figure 37 for 2D situations. Here it is shown
that only when the dimensions of the impulse array are small
Figure 36. General flowchart of dosimetry algorithm.
FOURIER PROCESSING MORE EFFICIENT
DIRECT PROCESSING MORE EFFICIENT
ii I I I I
0 200 400 600 800
SIZE OF INPUT ARRAY, N
figure 37 Comparison of direct and
processing for finite area convolution
(Pratt 1978).
IUUU IzLu1
Fourier domain
in two dimensions
* ~i
Figure 38. Example of wraparound format for DPB and zero
padded TERMA in one dimension.
'^^^

rLoa
Depth in Water (mm)
 1 MeV 2 MeV 3 MeV
 4 MeV 0 6 MV (Measured)
Figure 39. Comparison of measured percent depth dose to
that calculated from monoenergetic differential pencil beam
convolution (FS = 5x5 cm2, SSD = 100 cm, 5x5x5 mm3 voxels).
relative to the input (signal) array is direct processing
more efficient. For this research, since the range of the
DPB arrays extend to the dimensions of the calculation
volumes of interest, Fourier processing is more efficient.
Extension of this trend to 3D situations would be even more
dramatic. Although the convolution process is relatively
straightforward, one must take great care in the arrangement
and format of the arrays to be convolved. The convolution
theorem presumes two conditions which must be met (Press et
al. 1992). First, the arrays must have the same dimensions.
This criteria can be accomplished by zero padding the
smaller array. A second, more difficult constraint, is the
assumption that the signal (TERMA) and response (DPB) are
periodic. If this is not the case, significant wraparound
artifacts will occur. This requirement may be met by zero
padding the edges of the signal function and putting the
response function in wraparound order. The amount of zero
padding required depends on the nonzero extent of the
response function. An example of the zero padded TERMA and
wraparound format for the response function in a single
dimension is shown in figure 38. Percent depth dose data
produced by the dosimetry algorithm for 1, 2, 3, and 4 MeV
photons are displayed in figure 39 and compared to measured
6 MV data. More accurate modeling requires a polyenergetic
approach as described in the next chapter.
CHAPTER 4
BREMSSTRAHLUNG SPECTRA
Introduction
This chapter describes the method used to determine the
effective beam spectrum for a given linac and collimating
system. This spectrum is used to model the clinical xray
beam in all treatment planning dosimetry as opposed to
assuming a monoenergetic beam of energy equal to one third
of the nominal accelerating potential.
Although an approximate estimation of a linac spectrum
is possible using analytical methods (Desobry & Boyer 1992;
Koch & Motz 1959), the complexity and variability of linac
physical parameters limit the usefulness of such an
approach. Thus, most estimates of linac spectra are based on
Monte Carlo simulations or experimental measurements.
The Monte Carlo simulations have proven to be fairly
accurate as long as sufficient detailed information of the
linac head is available (Mohan, Chui & Lidofsky 1985; Kusbad
et al. 1990; Lo 1992). However, such simulations require
computing resources not readily available in most clinical
settings. Also, regardless of their sophistication, Monte
Carlo simulations should always be experimentally verified.
A variety of experimental techniques have been
investigated and include offaxis Compton scattering
spectroscopy (Faddegon, Ross & Rogers 1991, Levy et al.
1974, 1976), photoactivation (Nath & Schulz 1976),
attenuation analysis (Dawson 1989; Huang, Kase & Bjarngard
1983, Archer, Almond & Wagner 1985), and numerical
reconstruction from depth dose data (Sauer & Neumann 1990;
Ahnesjo & Andreo 1989). Of these experimental techniques,
only those based on numerical reconstruction of depth dose
data rely on readily available measured clinical data. As
will be shown, by using such measured data and the energy
deposition kernels previously discussed, a link between the
dosimetry model in question and the calibration of specific
hardware (i.e. SL75/5 linac) can be achieved.
Spectrum Model
The algorithm developed for determination of the
effective beam spectrum, SPECTRUM, is illustrated in figure
41. As input, two different sets of data are required.
First, experimentally measured data are obtained for the
open field unmodulatedd) collimating system of interest. For
example, a 5x5 cm2 field size at SAD may consistently be
used, in conjunction with a beam intensity modulator, for
all treatments (further modulation of the beam by the
attenuation device is considered in the dosimetry
algorithm). The measured depth dose data requires units of
absorbed dose per monitor unit (Gy/MU) and was obtained by
applying calibrated output and field size factors (AAPM Task
Group 21, 1983) to percent depth dose measurements for the
field size of interest. For the small field sizes typical of
radiosurgery, the beam quality does not vary
Figure 41. General flowchart of SPECTRUM algorithm.
significantly across the field. Thus, central axis data
alone may be used to characterize the beam (The potential of
extending this method to larger fields will be discussed in
chapter 8).
The second set of data required is absorbed dose per
reference energy fluence and per monitor unit for each
photon energy to be included in the effective spectrum. The
"reference" term used here indicates a point in space
between the photon source and the surface of the phantom.
Although any distance may be chosen, it is convenient to
place the reference point at the location of the beam
modulating device. For the purposes of this research, a
value of 65 cm from source to center of modulating device
was chosen. This data is created by the monoenergetic
dosimetry model discussed in the previous chapter and is
shown in Figures 42 and 43.
Obviously, the number of energy bins included in the
spectrum is a compromise between speed and accuracy of
dosimetry calculation. For example, it has been shown
(Mackie, Scrimger & Battista 1985) that a monoenergetic
approximation of 5 MeV duplicates measured depth dose values
to within 5% up to depths of 20 cm for a 10 x 10 cm, 15 MV
beam. Modeling their spectra with five energy bins reduced
the error to within 0.5 %. Thus, based on these results, and
the availability of DPB energies, seven energy bins were
chosen for my research (0.5, 1, 2, 3, 4, 5 and 6 MeV). Any
more energy bins than this would increase the computation
10 20 30 40 50 60
Depth in Water (mm)
70 80 90
Figure 42. Absorbed dose per energy fluence for four
monoenergetic photon energies.
4
3.5
3
T
,j 2.5
0
2

1.5
1 
ii
i/
!
/
/
 0.5 MeV  1 MeV 2 MeV 3 MeV
1.8
1.6 
1.4
(N
 1.2
12
0
E 0.8
0.6
0.4
0.2
0
70 80 90
4 MeV 5 MeV 6 MeV
Figure 43. Absorbed dose per energy fluence for three
monoenergetic photon energies.
r
10 20 30 40 50 60
Depth in Water (mm)
i
i
i
i

time with no significant increase in dosimetry accuracy and
would require interpolation between the energies provided in
the kernel library.
Once the input data files are read, SPECTRUM then
proceeds to optimize the fit between measured and calculated
depth dose data. The objective function (a measure of the
fit) to minimize is similar to that used by Ahnesjo (Ahnesjo
& Andreo 1989):
objfunc= E> dE2 (E)D,
a E
where
E = photon energy of given spectral bin.
d, = dose per reference energy fluence data (Gy MeV1 cm2
MU1) as a function of photon energy, E, and depth in the
phantom, z.
Y,(E) = reference energy fluence (MeV cm2) for each photon
energy, E.
Dz = measured dose per monitor unit data (Gy MU1) for
depths, z.
z = Depth in phantom (cm) where difference between measured
and calculated data is found.
Thus, at each depth, the difference is found between
calculated and measured absorbed dose. The actual depth
values are based on the voxel size used in the dose
calculation. The total objective function, a sum of the
absolute value of these differences, is then normalized for
the total number of depths used. A good fit to measured data
is found by choosing the proper weighting of energy fluence
for each photon energy bin comprising the spectrum. An
optimum weighting is referred to as the "effective spectrum"
as opposed to a physically accurate spectrum. This is due to
both the discrete polyenergetic dosimetry modeling and the
weak energy dependence of depth dose on the spectral shape.
These conditions result in multiple "solutions" to the
problem.
Due to such a weak energy dependence of depth dose
data, two practical constraints were imposed to help guide
the search through solution space towards more realistic
solutions (Sauer & Neumann 1990). First of all, the energy
fluence values for all energy bins are required to be equal
to or greater than zero:
T.W!0 for all n energy bins
The second constraint limits the magnitude of energy fluence
bins relative to their neighbors. This constraint assumes
that a pure bremsstrahlung spectrum with no characteristic
peaks should increase up to a modal energy (peak of energy
fluence spectrum) and then steadily decrease. Thus, assuming
a modal energy bin, nmodal,
T',5,, for n < nmodal
Yfn',, for n > nmodal
Since the modal energy bin is not known ahead of time, the
optimization is automatically repeated n times letting each
energy bin assume the modal bin role and the best fit is
chosen from these solutions.
In addition to these constraints, the minimum depth,
zc, used in the objective function was limited to values
beyond the depth of electron contamination. Thus, the fitted
spectrum represents all primary and scattered photons which
penetrate beyond zc but excludes contaminant electrons and
lower energy photons. Based on previous purging magnet
experiments (Ahnesjo & Andreo 1989) an empirical
relationship was found between the maximum depth of
contamination, zc (cm), and the nominal accelerating
potential, E0:
z,:0.3Eo
Based on this relationship, the maximum depth of electron
contamination would be approximately 1.8 cm for a 6 MV
photon beam. For the SL75/5, 6 MV linac modeled in my
research, a conservative minimum depth value of 2.25 cm was
chosen.
Optimization Algorithm
The algorithm used within SPECTRUM for multi
dimensional (7 energy bins) optimization, AMEBSA (Press et
al. 1992) is based on simulated annealing (Metropolis et al.
1953; Kirkpatrick, Gelatt & Vecchi 1983; Bohachevsky,
Johnson & Stein 1986). This approach makes an analogy
between the statistical behavior of large numbers of atoms
in a liquid or solid system and multidimensional
optimization problems. In such a system, there exist a large
number of possible atomic configurations/energy states
(local minima) although only one ground state (global
minimum) exists. The probability, P(E), of a system being in
a particular energy state, E, is described by the Boltzmann
probability distribution:
P(E)c e
where k is Boltzmann's constant and T is the thermal
equilibrium temperature. This relationship indicates that,
even at low temperatures, a finite probability exists for
the system to be in a higher energy state. The probability
of a system changing from energy state E1 to energy state E2
is expressed as
P(AE) e k"
This relationship indicates a high probability for E2 < E1
(downhill) and small but finite probabilities for E2 > E1
(uphill) changes. Using these statistical trends, the
Metropolis based algorithm explores solution space, always
allowing downhill steps and sometimes allowing uphill steps,
as the temperature of the system is gradually decreased (see
figure 44). As with real physical systems, if the
temperature of the system is cooled too rapidly (quenched),
the final resulting energy state will not be the ground
state. However, if the temperature is cooled sufficiently
slowly and the various system configurations are thoroughly
sampled at each temperature step (by incorporating a random
number generator), the global minimum is more likely to be
found. The workhorse for AMEBSA is the downhill simplex
algorithm, AMOEBA (Press et al. 1992). This method replaces
Keep old
point
Reject try
ry Reject old
point
U1 Accept try
/, ,, ;ll\i
Reject old
point
Acrer'in+ +ri
ry (downhill)
Figure 44. Example of accept or reject criteria used in the
AMEBSA algorithm. One random number proportional to the
current temperature is added to the objective function value
associated with a point in solution space. A second random
number is subtracted from the value associated with the
potential replacement point. The resulting values are then
compared.
ptry
,
old
1*
old
:
old
a single point (solution array) description of the system
with a simplex of N+1 points. The exploration of solution
space is then accomplished by reflections, expansions and
contractions of the simplex (see figure 45). As the
temperature approaches zero, the AMEBSA algorithm reduces
exactly to the downhill simplex method and converges to a
local minimum.
Results
Figure 46 illustrates an "effective energy fluence"
spectrum consisting of seven energy bins optimized for a 6
MV, 5x5 cm2 at SSD, 100 cm SSD beam from the SL 75/5 linac
using Ixlxl cm3 calculation voxels. A comparison between the
depth dose resulting from this spectrum and the measured
depth dose is also shown in figure 46. It is obvious that
the spectrum does not reflect the properties of a physically
real Bremsstrahlung spectrum (see figure 47). The reason
for this results from the relative insensitivity of depth
dose data to spectral variations as previously mentioned.
This trend is indicated in figures 48 and 49 where depth
dose data is reconstructed using spectra more and less
realistic than that shown in figure 46 with reasonable
results. This trend is also supported by recent work (Zhu &
Van Dyk 1993) indicating most changes in percent depth dose
are caused by changes in the average energy of the beam
rather than the specific spectral shape. For this research,
the actual shape of the spectrum is not the critical issue.
What is important is having a good fit to measured data with
simplex at beginning of step
low
reflection
reflection and expansion
contraction
multiple
contraction
Figure 45. Possible outcomes for a step in the downhill
simplex method. Such steps allow the simplex to change shape
as needed and "crawl" through the solution space landscape.
800
E 600
2 500
400
300
200
100 
0.5 1 2 3 4 5 6
Photon energy bin (MeV)
0.01
0.009
0.008
0.007
0.006
O 0.005
0.004
0.003
0.002
0.001
0 5 10 15 20 25 30
Depth (cm)
o Best Fit Meosured
Figure 46. Above: Seven energy bin Bremsstrahlung spectrum
output from SPECTRUM algorithm. Below: Comparison of
calculated and measured depth dose data for 6 MV beam.
(I (hI
Ao 4 MV Phwon Soct'lW u
10 10> ig W ho*on Sm.arIn 
1 e0 Ra3e" 0 3 cm
Uem [nrg 1 44 M1 O
Rloj Range 0 3 rm ,
USn4 (mpy I 9t lUV
o 0 0 10'1 V
wy (NM V) ,(Ow Iivy)
r3
10 MV fwa SWOcuvm 24 ; V Phom n Soocl"m
10 4 I 1'
Rta RaIne 03 cm Rad'l Rgqe t2 cm
Men Inewy I IMV UMn tim I 123 IIWV
10 c. 101
0 I 3 4 5 I 7 1 o10 0 S 1o s 20
Figure 47. Example of realistic 6 MV spectra generated by
Monte Carlo modeling for several linacs. a) Clinac 4 (4 MV),
b) Clinac 6 (6 MV), c) Clinac 18 (10 MV), d) Clinac 2500 (24
MV) (Mohan, Chui & Lidofsky 1985).
1.2 
N 1
0.2
o 0.8 
0
0.01
0.009
0.008
0.007 
0 
S0.006
C 
S0.005
S0.004 

0.003
0.0 
0.002
0.001
0.0 
5 10 15
Depth (cm)
20 25 30
o Best Fit  Meosured
Figure 48. Depth dose data (below) resulting from using a
more realistic Bremsstrahlung spectrum (above).
Photon energy bin (MeV)
0
Pholon energy bin (MeV)
SBesFiFit Msed
Figure 49. Depth dose data (below) resulting from using a
very unrealistic Bremsstrahlung spectrum (above).
65
the total (summed over all energy bins) energy fluence
providing a link to the calibrated linac output.
Having characterized a polyenergetic open field beam,
the next task is to model the effect of partial attenuation
(intensity modulation) on the dose distribution. This topic,
and how it is incorporated into inverse radiotherapy
planning for a single, static beam is discussed in chapter
six. However, prior to this, chapter 5 will compare
calculated narrow beam polyenergetic dose profiles with
measurements.
CHAPTER 5
DOSIMETRY MEASUREMENTS
Overview
This chapter presents absorbed dose profiles resulting
from narrow photon beams and the experimental methods used
to achieve them. Such measurements are necessary to evaluate
the accuracy of dosimetric calculations discussed in
previous chapters. Of particular interest is the absorbed
dose profile resulting from the extreme case of a 6 MV
photon beam collimated by a 3 x 3 mm2 opening in a 10 cm
thick lead attenuator. As will be discussed in chapter 6,
beam intensity modulation is modeled as variable attenuation
in adjacent variable thickness lead voxels (0 10 cm
thick). The smallest feasible lead voxel dimension was
chosen to be 3 x 3 mm2 as a compromise between resolution
and realistic engineering constraints. Thus, it seems
appropriate to compare measurements of this "building block"
of full sized beams with theoretical calculations.
The convolution dosimetry method is ideal for producing
dose distributions in situations where longitudinal and/or
lateral equilibrium does not exist. However, measuring
absorbed dose in such a situation is quite difficult.
Requirements for such dosimetry measurements include high
resolution and a flat energy response. Unfortunately, the
most common methods of clinical measuring devices (i.e. ion
chambers, semiconductor diodes, and silver halide film) do
not meet one or both of these requirements. The ion chamber
has an ideal energy response but is typically not small
enough relative to the conditions of the narrow photon beam.
The silver halide film provides excellent resolution but
over responds to the low energy components of the beam. The
diode suffers from being both too large and having a non
linear energy response. Alternate methods of measurement
include radiochromic film and high resolution TLD sheets. As
will be shown, both of these materials better approximate
the resolution and energy response criteria.
Radiochromic film
Radiochromic film dosimetry is based on an ionizing
radiation induced polymerization process resulting in a
color change without the need of post exposure processing
(Muench et al. 1991; McLaughlin et al. 1991). The most
widely available form of radiochromic film, GafChromicTM
Dosimetry Media (Nuclear Associates), consists of a 7 micron
sensitive layer on a 4 mil polyester base. The film is
colorless and practically "grainless" with spatial
resolution of > 1200 lines/mm (McLaughlin et al. 1991).
Exposure to radiation at wavelengths below 300 nm result in
an absorbed dose dependent change of color from clear to
blue. The increase in absorbance or optical density is
determined using a spectrophotometer or appropriate
densitometer.
Radiochromic film has been used for a variety of
measurements including National Institute of Standards and
Technology source profiling (Walker et al. 1992),
calibration of ophthalmic applicators (Soares 1991; Sayeg &
Gregory 1991), brachytherapy dosimetry (Coursey et al. 1992;
Muench et al. 1991), and a variety of electron and photon
beam investigations (Galvin, Smith & Lally 1993; McLaughlin
et al. 1991). The results of these studies have shown that
the radiochromic materials have a response similar to that
of tissue over a wide energy range (0.1 20 MeV for photons
and 0.01 20 MeV for electrons) with reproducible response
characteristics of approximately 5% at a dose of 200 Gy
(McLaughlin et al. 1991). It has also been shown that the
material is dose rate independent and has nearly identical
response (A absorbance per unit absorbed dose = 5%) for
electrons and photons. Thus, this dosimetry material has the
resolution and energy response characteristics needed for
the narrow beam dosimetry application at hand.
There are, however, a few difficulties associated with
this type of material. Other than the relatively high cost
of the material, one must be careful to avoid temperature
variations between exposure and reading of more than a few
degrees Celsius. However, this is typically not a problem in
a controlled climate of a cancer clinic. If significant
temperature differences are noted, simple correction factors
can be applied (McLaughlin et al. 1991). Also, it has been
suggested that one wait at least 24 hours between exposure
and reading in order to allow post irradiation density
growth to plateau (McLaughlin et al. 1991; Muench et al.
1991). Data collected for my research support this
recommendation and indicate a 10% increase in optical
density during the first 12 hours with no further detectable
change noted at a 36 hour reading. A third matter of concern
is the low sensitivity of the material. For example, a
change in optical density of 0.5 requires an absorbed dose
of several thousand cGy (see figure 51). Finally, and most
important of all, one must have access to a suitable
spectrophotometer or densitometer. As shown in figure 52,
the absorbance of the exposed radiochromic film depends on
the wavelength used during reading. The series of spectra
shown in figure 52 were made using a Beckman DU64
spectrophotometer and agree well with published results
(McLaughlin et al. 1991). Thus, due to the spectral
sensitivity, one must be sure to evaluate all exposed films
(including calibration shots) under the same spectral
conditions. A variety of wavelengths (broadband, 400, 510,
540, 580, 600, 605, 633, 650, and 675 nm) have been
evaluated (McLaughlin et al. 1991; Muench et al. 1991)
although no specific recommendations are given as to the
ideal wavelength to use. It appears most investigators are
simply using whatever wavelength(s) are available to them.
If one is interested in dose profiles, the most likely
choice would be either a scanning laser (633 nm)
densitometer or a solid state camera film digitization
4 6 8
Absorbed dose (cGy)
(Thousands)
I
BB (Macbeth)
+
650 nm (Beckman)
BB (Xrite)
650+20 nm (Xrite)
Figure 51. The relationship between absorbed dose and
optical density for a variety of film reading devices. BB
indicates broadband as defined by the white light source of
the specific densitometer.
,1.5
C
0
o
0o
0.5
a I
S
5400 Gy.
Figure 52. Series of absorption spectra obtained from
GafChromicT film (model 041) with a Beckman
spectrophotometer. The absorbed dose delivered to each film,
in order of increasing magnitude, are: 0, 1800, 3600, and
5400 cGy.
system. The later, of course, may allow some wavelengt@
variability if narrow band optical filters are used. The
dose profile capable densitometers currently available in
the Shands Cancer Center include a Wellhofer WD102 scanning
diode system and a Charge Injection Device (CID) video
system. Unfortunately, the exposed radiochromic film is
transparent to the IR diode light source and detector system
(950 + 20 nm) of the Wellhofer system. Thus, all narrow beam
dose profile data was collected using the CID video system
which has a broadband light source. This system, with a
512x480 CID array, has a resolution based on the field of
view (FOV) used during image capture. For the FOV used (10
cm), a resolution of approximately 0.21 mm was achieved.
Since the CID video system was previously shown to work
well with the standard silver halide film (Dubois 1993), a
comparison was first made between Kodak XV2 (silver halide)
and GafChromicTM for the case of 6 MeV electrons. Although
there is a significant difference in the effective Z between
the water equivalent radiochromic film and silver halide
film, there should be little difference in response for 6
MeV electrons since their ratio of collision stopping power
varies slowly with electron energy (Khan 1984). The results
of each film type are shown in figure 53. Since most of the
published measurements with radiochromic film were achieved
with monochromatic light sources, these measurements were
repeated using a 650 + 20 nm interference filter in
conjunction with the CID camera. The 650 nm filter was
Figure 53. Comparison of 6 MeV e dose profiles from Kodak
XV2 film (top) and GafChromicTM film (model 041) using a
broadband CID densitometer system (middle: no filter,
bottom: 650 nm filter).
* . I
5
4.5
4
E3.5
E 3.
C
o
3 2.5
0
c 2
2
: 1.5
o05
01
0.5 
40 50 60 70
Percent of Maximum Dose
80 90
A XV2 + 041 m 041 w/filer
Figure 54. Comparison of 6 MeV e depth dose data (10x10
cm cone, 3x3 cm2 blocked field, and 100 cm SSD) obtained
using broadband and filtered (650 + 20 nm) light sources.
L *
A
chosen since it provided the largest signal amplification as
reported in the literature (McLaughlin et al. 1991) and
since the CID camera has a peaked response at 650 nm. As
shown in figure 54, no significant difference in response
was observed between the broadband and filtered electron
depth dose data. However, the filtered image was
significantly blurred compared to the unfiltered image.
Thus, the decision was made to proceed with narrow beam
measurements using the broadband light source.
The 6 MV absorbed dose profiles for the 3x3 mm2
collimator opening is shown in figure 55. The experimental
data was obtained by carefully aligning the radiochromic
film, placed inside of a solid water film phantom (Bova
1990b), parallel to the incident collimated photon beam. The
parameters of interest include a collimator setting of 5x5
cm2, source to surface distance of 100 cm, center of
attenuator to surface distance of 27 cm, and an opening of
3x3 mm2 within a 12.7 cm thick (10 cm lead equivalent)
cerrobend attenuator. A total of 4500 monitor units were
delivered resulting in a maximum absorbed dose reading of
2078 + 81 cGy at a depth of 0.83 + 0.05 cm compared to a
calculated maximum value of 2246 cGy at a depth of 0.85 cm.
These absorbed dose values differ by approximately 7.5%
which is reasonable considering the experimental conditions
(geometrical alignment) and degree of nonequilibrium
present. Also, the data appear consistent with Monte Carlo
Figure 55. Comparison of calculated (left) and measured
(right) 6 MV narrow beam (3x3 mm2 blocked field) photon
absorbed dose profiles. Measured data was obtained from
GafChromicTM film using the CID densitometer system. The
90,70,50 and 30% isodose lines are shown.
predictions for narrow photon beams (Bjarngard, Tsai & Rice
1990).
TLD sheets
TLDs have been in use for quite some time for obtaining
accurate integrated dosimetry data. The small size of
available TLD chips or rods make them ideal for obtaining
"point" measurements comprising a 2 or 3 dimension dose
distribution. However, as each TLD requires its own
calibration and is read individually (see figure 56), a
tremendous amount of work is needed for large arrays. Also,
the size of each TLD chip or rod prohibits high resolution
measurements.
Recently, an alternative approach to TL dosimetry has
become available (Jones 1993; Jones et al. 1992) whereby a
nearly continues layer of TLD material is exposed to
ionizing radiation and read via a precision controlled CO2
laser. The TLD material is currently available in two
different forms. For standard sized field dosimetry, 1.5 mm
diameter, 38 pm thick MgB407:Tm elements are arranged in a 3
mm x 3 mm grid on a 30 cm x 30 cm sized sheet totaling 10201
elements. For small field, high resolution work, a
continuous 4.65 cm x 4.65 cm sheet of CaSO4:Mn is read at
0.375 mm intervals for a total of 15,625 readings. These
materials, while offering good resolution and measurement
convenience, do have drawbacks. For example, the CaSO4:Mn
material used on the high resolution sheets exhibits
significant fading (50 60% in the first 24 hours) and some
Grounded
TLD Sample
Heater
Power Supply
Figure 56. Schematic diagram showing typical apparatus for
measuring thermoluminescence (Khan 1984).
sensitivity to ambient light (Khan 1984). Thus, the data
sheets must be kept in the dark as much as possible and a
standard exposure to reading delay should be chosen. Also,
while the effective Z of CaSO4:Mn is only 15.3, an enhanced
sensitivity to lower energy photons does exist. For example,
the material is approximately 10 times more sensitive at 30
keV than at 1.25 MeV (Khan 1984). In contrast, however,
silver halide film is approximately 25 times more sensitive
at 25 keV than at 1.25 MeV (Johns & Cunningham 1983). The
substrate for both types of material is a 0.125 mm thick
polymer.
In the arrangement shown in figure 57, the CO2 laser
rapidly heats the TLD spot producing the characteristic glow
curve which is detected by the photomultiplier tube. If
careful calibration and annealing schedules are followed,
accuracy to within 1% of the delivered dose is claimed
(Jones 1993, coarse grid arrays). Unfortunately, no
published data on the accuracy of the high resolution sheets
is currently available. These sheets are still in the
experimental phase of development for which the Shands
Cancer Center is a test site. Ongoing tests have shown an
accuracy of approximately 7%.
Prior to collecting narrow beam data for this research,
a high resolution TLD sheet calibration was accomplished
with 400 cGy from a Co60 unit. After the calibration
reading, the sheet was annealed for 30 minutes at 1700 C. As
with the radiochromic film, the TLD sheet was placed
MIRROR 3
MIRROR 4
oa ur.Il UNULL o  '/ DOSIMETER ARRAY
Figure 57. Schematic diagram of the laser/ optical system
of the TLD array reader. The 3 mm x 3 mm grid array type of
sheet is shown (Jones 1993).
Figure 58. Comparison of calculated (left) and measured
(right) 6 MV narrow beam (3x3 mm2 blocked field) photon
absorbed dose profiles. Measured data were obtained from two
CaSO4:Mn high resolution TLD sheets placed adjacent to each
other. The 90, 70, 50, and 30% isodose curves are shown for
calculated results on the left. The same isodose
percentages, as well as the 95% line, are shown for measured
data on the right.
parallel to the incident 6 MV beam using the film phantom
and same geometry as before but with a monitor unit value of
900. The TLD sheet was then read with the device shown in
figure 57 using the previously stored calibration data to
obtain absolute data. Experiments have shown that several
readings can be accomplished without recalibrating the TLD
sheets. However, for this research, the data were the first
obtained following the calibration. A comparison of the TLD
reading and calculated results are shown in figure 58.
These results, as with the radiochromic film data, compare
favorably with the calculated data with a measured maximum
dose of 503 + 35 cGy at a depth of 0.80 + 0.05 compared to
calculated values of 449 cGy at 0.85 cm.
The position of the calculated dmax matches the
measured values from both experimental methods. The
calculated absorbed dose value appears to be 8% larger than
that measured with radiochromic film and 11% smaller than
that shown with the high resolution TLD sheets. Considering
the experimental conditions and dosimetry materials
involved, the calculated data appear to be sound and quite
suitable for the current application.
CHAPTER 6
INTENSITY MODULATION
Overview
This chapter presents the method used to determine the
optimal 2D beam intensity modulation function (IMF) needed
for any given beam's eye view of the target tissue. The
optimal choice is chosen in such a way that each beam's eye
view IMF is independent of any others. Determining the IMF
is a two step process involving deconvolution followed by
reverse ray tracing. Polyenergetic absorbed dose
distributions resulting from the projection of photon beams
through intensity modulation devices will also be presented
Optimal Intensity Modulation
Deconvolution
As you recall, the absorbed dose can be calculated by
convolving an energy deposition kernel, referred to as a
differential pencil beam (DPB), with a TERMA distribution
resulting from the primary photon fluence. In a similar
manner, the TERMA can be found if an absorbed dose
distribution and DPB are provided. The convolution theorem,
which proved useful in performing the convolution, is
equally useful for the deconvolution process. The TERMA can
be found from
T(')= [k [(r )]J
where
3 = Fourier transform.
3 = Inverse Fourier transform.
T(f') = TERMA value associated with a tissue voxel centered
at I'.
D(r) = prescribed dose distribution value associated with a
tissue voxel centered at T.
k(YY') = DPB describing deposition of energy at i due to an
initial photon, representing the average energy of the beam,
interacting at ?'.
Thus, instead of multiplying two values in frequency space,
a ratio is found. The inverse FFT of this ratio represents
the TERMA distribution required to achieve the desired
absorbed dose prescription. This prescription contains a
value of absorbed dose, in Gy, for each cubical voxel of
tissue in 3D space. It is an ideal request that contains
all zeros outside of the target volume and a single, uniform
value inside the target volume. Critical structures may also
be included in the dose prescription by assigning an
arbitrary numerical value (chosen to be equal to 1) outside
the dose range of interest to the voxel(s) comprising the
region of interest. Voxels that have been "tagged" as
critical structures are treated the same as zero dose voxels
during the deconvolution process but are spared
Figure 61. Isodose surface rendering of bowl shaped dose
prescription containing a uniform value of 10 Gy.
Figure 62. Isodose wire surface rendering of a bowl shaped
dose prescription containing a uniform value of 10 Gy
surrounding a region defined as a critical structure
dose during the backprojection process to be discussed in
the next section. Examples of bowl shaped dose prescriptions
with and without a critical structure are shown in figures
61 and 62.
Obviously, the TERMA distribution represents an ideal
situation that can not strictly be achieved with a single
photon beam or even multiple beams. The physics of photon
interactions simply do not allow all values of the
deconvolved ideal TERMA, especially those that are negative.
However, the information resulting from the deconvolution
does allow computation of the "best approximation" to the
absorbed dose prescription for a given photon beam. Examples
of the TERMA resulting from deconvolution using a 2 MeV DPB,
a bowl shaped dose prescription, and 2x2x2 mm3 voxels is
shown in figures 63 and 64. The bowl shaped dose
prescription shown in figure 61 was chosen to ultimately
illustrate the ability of the full field intensity
modulation technique to handle concave features. The values
for the isoterma surface renderings in figures 63 and 64
were chosen to illustrate the difference between the TERMA
and the dose prescription. This difference is due, of
course, to the fact that energy released in one voxel is
actually deposited in many voxels. Thus, while the entire
bowl shaped absorbed dose prescription is assigned 10 Gy (10
J/kg deposited), a wide range of energy released values
arranged in a nonuniform manner are required to achieve this
prescription. Figure 63 shows only those voxels
Figure 63. Isoterma surface displaying 10 J/kg values.
Figure 64 Isoterma surface displaying 33 J/kg values.
containing 10 J/kg (released) and 64 shows only those
voxels containing 33 J/kg (released). Figure 65 provides
more detail by presenting a contour plot of the TERMA near
the central axis.
Figure 65. Isoterma contour plot parallel to central axis
of beam. The curves shown represent 10, 20, and 30 J/kg
respectively.
Backproiection
Once the TERMA is found, it can be used to calculate
the physical characteristics of the intensity modulating
device that the photon beam will pass through. Differential
spatial attenuation will then result in an absorbed dose
distribution that, within the physical limits of a single
beam, will conform to that in the prescription.
For this work, lead was chosen as the attenuating
material. The attenuating device was modeled as a uniform
40x40 voxel grid with each voxel having a variable dimension
square face and variable thickness ranging from 0 to 10 cm.
As discussed in chapter 5, a minimum value of 3x3 mm2 for
the voxel face was chosen as a convenient compromise between
high resolution and construction limitations and is totally
independent of the absorbed dose calculation voxel size in
the phantom. The distance from the linac target to the
center of the attenuator can be chosen to reflect any
desired geometry. For the TERMA assigned to each cubical
voxel within the dose calculation volume of the phantom, an
appropriate lead thickness value in the lead grid can be
determined by backprojection. The backprojection referred to
is a simple ray trace from the center of a TERMA voxel
towards the linac target as illustrated in figure 66.
Included in the ray trace are inverse square, water phantom
(W) attenuation, and lead (L) attenuation:
( r 1 /
I I
I TERMA 11.3
txy= In  11.3
I [ 2 1.0(rs,) IP L
Tr s 1.6E10e
where
txy = thickness of lead voxel, in cm, at grid position x, y.
Deconvolution \
Ter na
Figure 66. Illustration of two step process to determine
IMF from absorbed dose prescription.
Y, = reference incident energy fluence for monoenergetic
photons at a specified distance from the source.
All other variables in the above equation were previously
defined.
In order to exclude physically unrealistic solutions,
the thickness calculation is only performed for those
cubical voxels with TERMA values above a preset minimum. The
minimum value was chosen to be that resulting from
transmission through the thickest possible (10 cm) lead
attenuator. Thus, for each cubical voxel containing a
physically realistic TERMA value, a lead thickness value is
computed. For voxels tagged as critical structures, a lead
thickness value is assigned based on a weighting factor
(ranging from 0 to 1.0) relative to that of target tissue. A
default value of 1.0 results in a maximum lead thickness (10
cm) being assigned. This same tissue weighting factor will
also play a role during beam weight optimization to be
discussed in the next chapter. Next, the thickness values
computed for all TERMA voxels (and those for critical
structures) whose ray lines intersect a given lead
attenuator voxel are averaged. When completed, the 40x40
voxel attenuator grid contains an average thickness value
for each voxel which is representative of the beam's eye
view of the TERMA distribution below it. Figure 67
illustrates the intensity modulating grid design for the
case of a beam's eye view into the top of the bowl shaped
absorbed dose prescription. This design was computed using a
MI
Figure 67. Example of beam's eye view intensity modulation
device needed for a top view of the bowl shaped dose
prescription surrounding a critical structure (cutaway
provided for clarity).
single (2 MeV) energy DPB instead of a polyenergetic
spectrum thus resulting in reasonable calculation times. The
need for such short calculation times will be demonstrated
in the next chapter where multiple beam treatment plans will
be investigated.
Single Beam Intensity Modulated Dosimetry
A conformal, 3D absorbed dose distribution can now be
computed by modeling the effect of the lead grid on the
incident energy fluence. At this point the versatility of
the DPB convolution dosimetry algorithm becomes even more
apparent. The differential attenuation of each energy
component comprising the primary photon beam is accomplished
with a simple exponential term for each energy bin. The term
actually includes differential attenuation for the phantom
as well as the lead. However, the phantom effect is minor
compared to the lead. The attenuated incident fluence is
then used to compute the TERMA which is convolved with the
DPB(s) resulting in absorbed dose. Figure 68 illustrates
relative isodose contour plots of the absorbed dose
distribution using a polyenergetic beam and the intensity
modulation function illustrated in figure 67. The dose
sparing effect of averaging maximum lead thickness values
for the critical structure (critical tissue weight = 1.0)
regions is quite noticeable. If no critical structure is
present, the resulting intensity modulating grid design
projects the isodose contours shown in figure 69. In the
next chapter, the use of such modeling will be extended to
