Group Title: Inverse radiosurgery treatment planning through deconvolution and constrained optimization /
Title: Inverse radiosurgery treatment planning through deconvolution and constrained optimization
Full Citation
Permanent Link:
 Material Information
Title: Inverse radiosurgery treatment planning through deconvolution and constrained optimization
Physical Description: vi, 233 leaves : ill. ; 29 cm.
Language: English
Creator: Harmon, Joseph F., 1960-
Publication Date: 1994
Copyright Date: 1994
Subject: Medical physics   ( lcsh )
Nuclear Engineering Sciences thesis Ph.D
Dissertations, Academic -- Nuclear Engineering Sciences -- UF
Genre: bibliography   ( marcgt )
non-fiction   ( marcgt )
Thesis: Thesis (Ph. D.)--University of Florida, 1994.
Bibliography: Includes bibliographical references (leaves 226-232).
Additional Physical Form: Also available on World Wide Web
Statement of Responsibility: by Joseph F. Harmon.
General Note: Typescript.
General Note: Vita.
 Record Information
Bibliographic ID: UF00097377
Volume ID: VID00001
Source Institution: University of Florida
Holding Location: University of Florida
Rights Management: All rights reserved by the source institution and holding location.
Resource Identifier: alephbibnum - 001981548
oclc - 31920344
notis - AKF8479


This item has the following downloads:

PDF ( 9 MBs ) ( PDF )

Full Text








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.



ABSTRACT................................................ v


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


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



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 3-D 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.



The goal of this research was development of an

automated conformal 3-D 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

200-300 kV x-rays, 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


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 Brown-Roberts-Wells (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 (anterior-posterior,

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 (5-40 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 1-1.

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 1-1. 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


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


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 3-D 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 back-projected to obtain a 2-D

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.


Dosimetry Algorithms


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


D(d,x,y) = %dd(d)f(d,x)g(d,y)

This model is commonly referred to as the Milan-Bentley

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 off-axis-ratio

(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


C = field size at SAD (source to axis distance)

STD = source to target distance

d = depth at point p

r = off-axis 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 off-axis 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 off-center ratio. The off-axis 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


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


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)


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


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

x-ray source

phantom surface


kernel axis

beam axis

Figure 2-1 Relationship between the interaction site (T')and
dose deposition site (?) using energy deposition kernels
(Sharpe & Battista 1993).


Primary First
charged scatter
particle photon
kinetic energy

Primary First Second
energy scatter scatter
deposited charged photon
particle energy

First Second Multiple
scatter scatter scatter
energy charged photon
deposited particle energy

Second Multiple
scatter scatter
energy charged
deposited particle
plus annihilation
photon energy Multiple
Bremsstrahlung deposited
plus annihilation
charged particle
kinetic energy

plus annihilation
energy deposited

Figure 2-2. Energy deposition partitioning used in kernel
development (Mackie et al. 1988).

voxel at point r' (see figure 2-1). 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 Klein-Nishina

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 2-2. An example of such a kernel

is shown in figure 2-3. 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 2-4. If several pencil beams interact at a voxel

(point), a dose distribution referred to as a convergent


1200 (00

160 0 0 10

Radius (g tm-)

Figure 2-3. Isoline format of a total energy deposition
kernel for 1.25 MeV photons in water. The units are cGy MeV~
1 photon-1 (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 three-dimensions, a very smooth kernel is

formed representing the steepest dose gradient possible

about the voxel. An example of a kernel for the 2-D rotation

case is shown in figure 2-5. 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

o 10-
N ;
15 i0-


-10 -5 0 5 10

Figure 2-4. Isoline format for a 1.25 MeV pencil beam
kernel. Units are fractional energy deposited per unit
volume (cm3) (Eklof, Ahnesjo & Brahme 1990).





-10 -5 0 5 10

Figure 2-5. 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

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(r-r') 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 in-air

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 off-axis 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. plane-parallel 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 2-D 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 2-D circularly symmetric dose

distribution inside an arbitrarily shaped phantom with

rotation about any axis (Barth 1990). Having accomplished

this, one can consider 2-D nonsymmetric dose distributions

by considering them as composed of adjacent,

nonintersecting, radially symmetric dose distributions as

shown in figure 2-6. 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 2-6. 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

secondary-particle 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 2-D 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 1-D modulation

function for all incident beam directions. An example of

this method, extended to include 3-D cases (Boyer, Desobry &

Wells 1991), is shown in figure 2-7. As can be seen, this

". Nh -1 1 1 K ri

Figure 2-7. 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


Iterative Methods

One implicit method of obtaining an IMF is through

algorithms similar to those used in computed tomography.

Recall, that in CT, a 1-D transmission function is measured

and represents the projection of a 2-D density distribution

of that slice of the patient for a given projection angle.

Analogously, a 1-D IMF for a given projection angle can be

determined by first projecting the ideal 2-D dose

distribution. This comparison is illustrated in figure 2-8.

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

Ra(dion Rad

Figure 2-8. 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 2-D cases. For example,

allowing ten incident beams with each beam being composed of

ten (2-D case) or 100 (3-D 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 3-D

geometries although the inclusion of scatter is currently

limited to the 2-D 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 3-D 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.



This chapter describes the 3-D 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


The dosimetry algorithm used was chosen for its sound

physical basis, accuracy, and ease of incorporation into a

3-D 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

3-1 through 3-3. 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 -

" 0.8 -
8 0.5-

6 0.4


" 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 3-1. Differential pencil beam components, 0.1--1 MeV.









2 3 4 5 6 10
Actual Photon Energy (MeV)

I Short Long

Figure 3-2. Differential pencil beam components, 2--10 MeV.

20 25 30 40
Actuol Photon Energy (MeV)

SShort M Long

Figure 3-3. Differential pencil beam components, 15--50 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 3-D

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 x-y plane for every depth.

The radial and depth characteristics of the 2 MeV DPB in

Cartesian format is illustrated in figure 3-4. 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 3-5. 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 3-4. Uniformly scaled (2 MeV, 2x2x2 mm3 voxels)
differential pencil beam data for selected depths.

c 0.1



. 0.001


o 1E-05



.E 2

S 0.1 i



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 3-5. 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.


ce -h : x x .


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


r ,MeV iu c I f
TERMA= P- cm I e -
L 'cm'JJL g J r


Yr = reference incident energy fluence for Monoenergetic

photons at a specified distance from the source.

- = mass attenuation coefficient of phantom for
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


(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.


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 3-D 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 FFT-1. The array of data is actually

accessed as a l-D array with alternating real and imaginary

parts. Thus the dimension vector is used to interpret the

arrangement of data within the l-D 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 3-6. 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 3-7 for 2-D situations. Here it is shown

that only when the dimensions of the impulse array are small

Figure 3-6. General flowchart of dosimetry algorithm.



ii I I I I

0 200 400 600 800

figure 3-7 Comparison of direct and
processing for finite area convolution
(Pratt 1978).


Fourier domain
in two dimensions

* ~i

Figure 3-8. Example of wrap--around format for DPB and zero
padded TERMA in one dimension.




Depth in Water (mm)

-- 1 MeV 2 MeV 3 MeV
- 4 MeV 0 6 MV (Measured)

Figure 3-9. 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 3-D 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 wrap-around

artifacts will occur. This requirement may be met by zero

padding the edges of the signal function and putting the

response function in wrap-around order. The amount of zero

padding required depends on the non-zero extent of the

response function. An example of the zero padded TERMA and

wrap-around format for the response function in a single

dimension is shown in figure 3-8. Percent depth dose data

produced by the dosimetry algorithm for 1, 2, 3, and 4 MeV

photons are displayed in figure 3-9 and compared to measured

6 MV data. More accurate modeling requires a polyenergetic

approach as described in the next chapter.



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 x-ray

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 off-axis 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

4-1. 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 4-1. 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 4-2 and 4-3.

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 4-2. Absorbed dose per energy fluence for four
monoenergetic photon energies.




,j 2.5-



1 -


- 0.5 MeV -- 1 MeV --2 MeV --3 MeV


1.6 -

- 1.2-

E 0.8-




70 80 90

-4 MeV 5 MeV 6 MeV

Figure 4-3. Absorbed dose per energy fluence for three
monoenergetic photon energies.


10 20 30 40 50 60
Depth in Water (mm)




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


E = photon energy of given spectral bin.

d, = dose per reference energy fluence data (Gy MeV-1 cm2

MU-1) as a function of photon energy, E, and depth in the

phantom, z.

Y,(E) = reference energy fluence (MeV cm-2) for each photon

energy, E.

Dz = measured dose per monitor unit data (Gy MU-1) 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


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:


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


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 multi-dimensional

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 4-4). 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

Reject try

ry Reject old
U-1 Accept try
/, ,, ;ll\i

Reject old
Acrer'in+ +ri

ry (downhill)

Figure 4-4. 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







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 4-5). As the

temperature approaches zero, the AMEBSA algorithm reduces

exactly to the downhill simplex method and converges to a

local minimum.


Figure 4-6 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 4-6. It is obvious that

the spectrum does not reflect the properties of a physically

real Bremsstrahlung spectrum (see figure 4-7). 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 4-8 and 4-9 where depth

dose data is reconstructed using spectra more and less

realistic than that shown in figure 4-6 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



reflection and expansion



Figure 4-5. 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.


E 600

2 500




100 -

0.5 1 2 3 4 5 6
Photon energy bin (MeV)






O 0.005




0 5 10 15 20 25 30
Depth (cm)

o Best Fit -Meosured

Figure 4-6. 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)


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. 10-1
0 I 3 4 5 I 7 1 o10 0 S 1o s 20

Figure 4-7. 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

o 0.8 -





0.007 -
0 -

C --


S0.004 -

0.0 -


0.0 -

5 10 15
Depth (cm)

20 25 30

o Best Fit -- Meosured

Figure 4-8. Depth dose data (below) resulting from using a
more realistic Bremsstrahlung spectrum (above).

Photon energy bin (MeV)


Pholon energy bin (MeV)

SBesFiFit Msed

Figure 4-9. Depth dose data (below) resulting from using a
very unrealistic Bremsstrahlung spectrum (above).


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




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


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 5-1). Finally, and most

important of all, one must have access to a suitable

spectrophotometer or densitometer. As shown in figure 5-2,

the absorbance of the exposed radiochromic film depends on

the wavelength used during reading. The series of spectra

shown in figure 5-2 were made using a Beckman DU-64

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)

BB (Macbeth)
650 nm (Beckman)

BB (X-rite)

650+-20 nm (X-rite)

Figure 5-1. 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.




a I


5400 Gy.

Figure 5-2. 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 WD-102 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 XV-2 (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 5-3. 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 5-3. Comparison of 6 MeV e- dose profiles from Kodak
XV-2 film (top) and GafChromicTM film (model 041) using a
broadband CID densitometer system (middle: no filter,
bottom: 650 nm filter).

* .- I




E 3.
3 2.5-
c 2-

: 1.5-


0.5 -

40 50 60 70
Percent of Maximum Dose

80 90

A XV-2 + 041 m 041 w/filer

Figure 5-4. 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 *

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 5-4, 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 5-5. 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 5-5. 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


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 5-6), a

tremendous amount of work is needed for large arrays. Also,

the size of each TLD chip or rod prohibits high resolution


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

TLD Sample

Power Supply

Figure 5-6. 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


In the arrangement shown in figure 5-7, 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



oa ur.Il UNULL o ----- '/ DOSIMETER ARRAY

Figure 5-7. 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 5-8. 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 5-7 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 5-8.

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.



This chapter presents the method used to determine the

optimal 2-D 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


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


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(Y-Y') = 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 3-D 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 6-1. Isodose surface rendering of bowl shaped dose
prescription containing a uniform value of 10 Gy.

Figure 6-2. 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

6-1 and 6-2.

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 6-3 and 6-4. The bowl shaped dose

prescription shown in figure 6-1 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 6-3 and 6-4

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 6-3 shows only those voxels

Figure 6-3. Isoterma surface displaying 10 J/kg values.

Figure 6-4 Isoterma surface displaying 33 J/kg values.

containing 10 J/kg (released) and 6-4 shows only those

voxels containing 33 J/kg (released). Figure 6-5 provides

more detail by presenting a contour plot of the TERMA near

the central axis.

Figure 6-5. Isoterma contour plot parallel to central axis
of beam. The curves shown represent 10, 20, and 30 J/kg


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 6-6.

Included in the ray trace are inverse square, water phantom

(W) attenuation, and lead (L) attenuation:

( r 1 /

I TERMA 11.3
txy= In | 11.3
I [ 2- 1.0(r-s,) IP L
Tr s 1.6E-10e


txy = thickness of lead voxel, in cm, at grid position x, y.

Deconvolution \

Ter na

Figure 6-6. 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


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 6-7

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


Figure 6-7. 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, 3-D 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 6-8 illustrates

relative isodose contour plots of the absorbed dose

distribution using a polyenergetic beam and the intensity

modulation function illustrated in figure 6-7. 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 6-9. In the

next chapter, the use of such modeling will be extended to

University of Florida Home Page
© 2004 - 2010 University of Florida George A. Smathers Libraries.
All rights reserved.

Acceptable Use, Copyright, and Disclaimer Statement
Last updated October 10, 2010 - - mvs