A unified approach to the optimization of brachytherapy and external beam dosimetry

Copyright

??Technical Innovations

A UNIFIED

TIMOTHY

and Notes

APPROACH TO THE OPTIMIZATION OF BRACHYTHERAPY AND EXTERNAL BEAM DOSIMETRY HOLMES,

MS.,

T. ROCK AND

Departments

0360.3016191 $3.00 + .OO CC:1991 Pergamon Press plc

PAUL

MACKIE,

PH.D.,

RECKWERDT,

DOUGLAS

SIMPKIN,

M.S.

B.S.

of Medical Physics and Human Oncology, University of Wisconsin 1300 University Avenue, Madison, WI 53706

School of Medicine,

Semi-automated optimization of dose distributions is possible using techniques borrowed from imaging science. The ideal distribution of dose is first deconvolved by a convolution kernel yielding an ideal weighting distribution in the patient. The weighting distribution describes the total energy released per unit mass of the irradiated medium. For internal and external radiation sources, this is directly related to the amount and distribution of radioactivity and energy fluence in the medium, respectively. For external sources, the exponential Radon transform is used to obtain ideal fluence projections incident on the patient. In both instances negative values are produced, which when set to zero result in perturbed dose distributions. This may necessitate iterative techniques to reduce the ‘residual dose’ produced by the zeroing process. Application of the approach is presented for the optimization of 2-dimensional dose distributions in external beam therapy, radioimmunotherapy, and brachytherapy sources. Radiation therapy, Dosimetry,

Brachytherapy,

Optimization,

the broadest support in the research community, is the convolution model ( I, 4, 10, 13, 15. 35, 40, 42, 43, 53). The attractive features of this model are its accuracy (39), its potential efficiency ( lo), its applicability to the calculation of dose distributions produced by either internal radioactive source distribution ( 14), or external X rays and y-rays ( 12, 39, 42), and the possibility of mathematically inverting the model for dose optimization (18, 34). The treatment planning process for either brachytherapy or external beam radiotherapy is primarily concerned with determining the number, placement, relative strengths (or weights), and exposure times of radiation sources that produce a dose distribution that most closely approximates the radiotherapist’s prescribed dose distribution. This can be a time-consuming process and many investigators have proposed solutions to this problem, typically using linear programming or least squares methods to automate the iterative procedure otherwise per-

INTRODUCTION Recent advances in computer hardware (3-dimensional graphics and image display technologies), architectures ( reduced instruction set computers, multi-processors ), and networking (distributed processing using multiple workstations) are providing the computational speed required for rapid calculation of 3-dimensional dose distributions in radiotherapy. Historically, 3-dimensional dose calculations of large radioactive source distributions have been successfully implemented in the clinical setting because of their relative simplicity compared to external beam models. The latter have remained essentially 2-dimensional approximations largely because of restrictions on the computation time. The past decade has seen much activity in the development of 3-dimensional dose algorithms, especially for external beam treatment planning (24, 37). One of the most promising models, and the one currently receiving

Reprint

Radon transform, Filtered backprojection.

requests

to: T. R. Mackie. authors would like to thank Jim Holden and Herb Attix for many fruitful discussions of imaging science and radiation dosimetry which contributed to this work, and Anders Brahme for his thoughtful comments. We also express our appreciation to Cynthia Thomason for providing data used for the Ir- 192 kernel. The ‘Donner Algorithms for Reconstruction Tomography’ proved invaluable for the external source optimization method presented in this work. We gratefully acknowledge its authors, R. H. Huesman, G. T. Gullberg, W. L.

Greenberg, and T. F. Budinger, for the use of this software. We also acknowledge P. Bennett and C. J. Kost, and J. L. Chuma, of TRIUMF, authors respectively of the software applications OPDATA (data analysis) and PLOTDATA (graphics), which were both used extensively for this work. Finally, we thank Orlando Canto for his assistance in preparing the figures of this work for publication. This work was supported by NC1 grant number CA48902. Accepted for publication 10 October 1990.

Acknowledgments-The

859

860

I. J. Radiation

Oncology

0 Biology

0 Physics

formed manually by a dosimetrist or physicist (5, 6, 27. 28, 29, 32, 41, 44, 4.5, 47, 49. 52). A more attractive approach would determine the treatment parameters directly from the treatment prescription by mathematical inversion of the dose distribution. Much of the initial work on dose optimization by mathematical inversion has centered on obtaining analytical solutions for external beam radiotherapy optimization. especially rotational or dynamic therapy (9. 18. 19. 20, 22, 23, 34). The initial work of Brahme (11N/. ( 20) produced an analytical expression describing the shape of the lateral absorbed dose profile of a rotational X ray beam which would produce a uniform dose volume in a cylindrical phantom when the X ray beam was rotated about (or near to) the center of the cylinder. This result was applied clinically by Lax and Brahme (33) for the treatment of head and neck tumors. Cormack (22) and Cormack and Cormack (23) extended the analytical results which are not of Brahme et (I/. (20) to dose distributions circularly symmetric. In arriving at their solutions. all of these authors made various simplifying assumptions to make the problem tractable: (a) ignore scatter and attcnuation. (b) non-divergent geometry. (c) ignore buildup, (d) assume dose is proportional to intensity. (e) assume the patient is cylindrically symmetric, and ( f) assume the patient is treated by a full rotation of parallel-opposed beams. A general solution to the equations of Cormack and Cormack (23) has not yet been found. Barth (9) has further extended the above work to dose distributions which are non-radially symmetric in convex phantoms of arbitrary cross-section. His results were limited to a simple class of dose distributions which posess zero dose outside of the distributions. All of these works (9, 20, 22, 23) provide analytical solutions to highly idealized treatment situations. Recent work by Lind and Brahme (24) and Brahme ( 18, 19) presents a way in which powerful numerical techniques used in image reconstruction can be used to optimize dose distributions from external beams. In these investigations, dose distributions are considered to be the superposition of fluence-weighted, convergent point irradiation kernels. A convergent point irradiation kernel is produced by rotating a pencil beam dose distribution about the center of a 2-dimensional circular, or 3-dimensional spherical medium of uniform density. Because of the symmetry of the medium about the point of rotation. each beam contributes equal fluence at that point. The authors proposed using Fourier deconvolution to derive the required weighting distribution of convergent point irradiations that produces the ideal dose distribution for external radiation sources (and internal sources as Well. given the appropriate kernel). The deconvolution can be performed rapidly in frequency-space as a simple division if the frequency-space form of the kernel does not have zeros, or very small values. Since this was not true for their kernel, they used an iterative approach instead to

April

IYYI.

Volume

20. Number

4

arrive at the desired weighting distribution (34). This distribution. when convolved with the convergent point irradiation kernels. produced a physically realizable dose distribution. Attenuation compensated inverse-backprojection was performed resulting in projections of the dose distribution which they called ‘transversal dose distributions.’ The projections approximated the photon fluence profiles required to produce the desired dose distribution. A resultant dose distribution was calculated in a single step by backprojecting and attenuating the ‘transversal dose distributions.’ This dose distribution closely resembled the projected dose distribution, although no filtering or resealing of the projections prior to backprojection was mentioned. The methods outlined above are very similar to image filtering (deconvolution) and image reconstruction (projection/backprojection) techniques used in other fields such as radioastronomy. electron microscopy, nuclear medicine. and diagnostic radiology. These techniques are well understood and have been implemented numerically in well-proven software algorithms in all of these tields. Unlike the analytic solutions mentioned above. numerical methods readily lend themselves to the problem oftreatment planning with the non-idealized geometries found in the clinical setting. This paper will present an algorithm similar to Brahme’s ( 19) that uses deconvolution and Radon transform mathematics to optimize the delivery of absorbed dose by internal and external radiation sources. Our approach differs from Brahme’s in several important ways: 1. Fundamental definitions of the weighting distribution and convolution kernel are used so that energy release into a medium by either internal or external sources is treated in the same manner. 2. We perform deconvolution of our kernel from the dose distribution by direct Fourier inversion instead of the iterative approach of Lind and Brahme (34). 3. We generalize the deconvolution step above to the optimization of both external and internal sources. 4. For external beam optimization, our kernel describes the result of photons interacting at a point-of-convergence from a number of directions rather than one describing the spread of dose about a point from pencil beams that interact before and after that point. It is similar to Brahme’s ‘point spread function’ h( 7) ( 19). This definition allows us to draw the analogy between emission computed tomography (ECT) and external beam radiation therapy. 5. We obtain true energy fluence projection estimates by projecting the TERMA distribution (described in the following section) instead of the dose distribution. 6. We filter the projection estimates prior to backprojection. 7. We calculate dose in a two-step process of projecting energy fluence into the patient to obtain a TERMA distribution followed by convolution with the kernel to obtain dose. Dose can be determined accurately by

A unified approach to optimization

calculating single beams one at a time, or it can be determined approximately by projecting all beams to obtain an integral TERMA distribution, followed by convolution with a multiple beam kernel that is assumed to be spatially invariant. a dose residual, or excess, and iterate 8. We determine to reduce it to an acceptable level. The approach presented here is consistent with fundamental concepts of radiation dosimetry (7), and is general in scope in its description of dose deposition by internal and external radiation sources. METHODS

AND

MATERIALS

Optimization using Deconwlution Reconstruction Muthematics

and Image

The modelling of absorbed dose distributions used in clinical radiotherapy has been shown by several authors ( 1, 1 I, 12, 35, 38. 39) to mathematically resemble a convolution process. For the case of a homogeneous medium. the dose at a voxel at position 7. D(i), from energy released at voxels ?’ in the medium can be described by the convolution:

D(i)

=

s

T(?)A(i

~ i’)d’i’

(1)

where T( 7’) is a 3-dimensional weighting distribution describing the total energy released per unit mass. A(? - it’) is the monodirectional convolution kernel. The integration is carried out over three dimensions. The distribution T(i’) has been called ‘terma’ (total energy released per unit mass) by Ahnesjo et a/. (2). It is a spatial distribution of scalar values describing the energy released per unit mass of medium at each voxel i’ in the medium. T( 7”‘) may represent either a continuous distribution of radioactive decay sites as in radioimmunotherapy, a distribution of discrete brachytherapy sources, or sites of energy transfer in external X ray therapy. Radiation sources may be located at the sites of energy release i’ (internal sources) or at distant locations (external sources). When located at i’, energy is transferred to the medium through radioactive decay, producing photons and charged particles. Otherwise, energy is transported by photons from the radiation source through the medium to 7’ where interaction occurs. The mechanisms of energy transport away from the sites of energy release it’ are similar for both types of sources. They primarily differ in the way energy is made available at the locations 7’. For a monoenergetic, internal source of radiation, the terma distribution is described by

T(?‘)

E = ~ A (7’)dt P(T’) s V

0 T. HOLMES

861

PI cd.

radiation in Joules. A/V( i’) is the average activity per unit volume given in units of Bq - m 3. The integration is over the residence time of the radioactive source in the medium to give the number of disintegrations per unit volume. ~(7’) is the mass density of the medium given in kg/m”. For a monoenergetic, external source of photons, the terma distribution is described by

T(i’)

= f (7’)E

(3)

where E is the energy of the photons in Joules. &(T’) is the distribution of fluence rate (flux density), and has the is over the time of units mm2 - s-‘. Now the integration irradiation to give the photon fluence at i’. p/p( T’) is the effective mass attenuation coefficient at the interaction site and has units of m’/kg. The monodirectional convolution kernel A (7 ~ i’) is a 3-dimensional function describing non-stochastic energy transport into a uniform medium and can be generated using the Monte Carlo method (2, 36, 48). It is defined in terms of (a) a single particle type and energy, and (b) a single incoming or outgoing direction for either a photon from an external source, or an emitted particle from an internal point source, respectively. It describes the fraction of energy absorbed at i. per unit energy released per unit volume in the medium by a particle at i’. Its units are m-j. Energy conservation requires that (setting i = 0) : ‘* s Q

A( i’)d3i’

= 1.O.

(4)

For an internal source having a spectrum of decay particles, or a polychromatic external photon source, the monodirectional kernel can be approximated by summing a set of weighted monodirectional kernels:

A(i

- ?‘) = ~ fjA(i

- i’)j,

~ f, = 1.0,

j=l

(5)

j=l

where f, represents the decay fraction of the jth decay particle for a radioactive source, or the spectral fraction of the jth spectral component for an external photon source. The multidirectional convolution kernel B(i’ - 7’) is described as a weighted superposition of rotated monodirectional kernels A( i - ?‘)i: B(i

- 7’) = i W(l’)iA(I i=l

- T’);,

5 W(?‘)i = 1.0,

(2)

where T (7’) has units of J / kg. E is the energy of the decay

&(i’)dt

s

(6)

i-l

where the subscript

i designates

the ith rotation

direction,

1. J. Radiation

862

Oncology

and the weight factor w( ?‘)i represents total energy released at the interaction ith direction. The general form of w ( i-‘)l is:

0 Biology 0 Physics

the fraction of the voxel i from the

April 199

I, Volume 20. Number 4

and the term inside the brackets is recognized as the isotropic multidirectional kernel. Combining Eq. 11 and Eq. 8, dose is approximately determined by a single convolution instead of the ‘n’ convolutions of Eq. 9:

W,e-lp/lr(;‘).p(;‘)d(,~-~‘.~‘)di

w(i’)

=

’

”

c

(7)

D(i)

W,,e-lriP(;').Pi;')~(,~~~'.~')dil J

=

s

T(?‘)B(?

- i’)d’i’.

(13)

J-1

where Wi is the source strength or the beam weight for the ith direction of n total directions, for internal and external sources, respectively. The exponential factors describe the attenuation between the source positions and the interaction voxel at i along a rayline described by S(lj _ i’ . i’)di. For internal sources these factors are unity. and for external sources, they are less than unity. If the fractional contributions of terma at the interaction voxel are equivalent for each direction i then w( ?‘)i = 1 /n, and Eq. 6 reduces to what is termed the ‘isotropic’ form of the multidirectional convolution kernel:

B(i

- i’) = $ f: A(7 I_ I

- ?‘)i.

(8)

For an internal source Eq. 8 can be used to generate a 47r isotropic kernel by summing over all directions. Alternatively, this kernel could be directly calculated by Monte Carlo methods (46, 48) and could include anisotropy effects due to source shielding, if present (21, 5 1). B( 7 - i’) can also represent the absorption of energy in the medium due to primary photons from multiple external radiation beams interacting at a point. For multiple external beams, the total dose at i is the sum of dose at that point from each constituent beam (Eq. I ):

D(i)

= i [J I=1

T(?‘)iA(Y

- i-)id’(Y’)]

(9)

where the summation is over n specific beam directions. If T(i’), is assumed to be equivalent for each beam. the distributive property of convolution can be used to recast Eq. 9 as:

D(?)

An additional

i=l

=

i

modification

where the total terma

A(?

- i’)i

1

d’i’.

results in a more useful form:

= nT(i’)i.

(12)

convolution above is implemented as the discrete convolution:

D(i)

= 2 T(F’)B,(i

- 7’)

on a

(14)

where Br (i - 7’) describes the fraction of energy absorbed at i per energy released at i’ and is unitless. It is given by the relation: Br(f’ - 7’) = B(7 ~ i’)d’?‘.

(15)

If the multidirectional kernel is assumed spatially invariant, the Convolution Theorem ( 16) can be used to cast the problem into a frequency space representation and exploit the calculational efficiency of the fast Fourier transform ( IO): D(T)=F-‘{F{T(i’)}.F{Br(?-T’)})

(16)

where F and F-r symbolize Fourier and inverse Fourier transforms. Spatial invariance requires the kernel shape to be the same throughout the irradiated medium. This is fulfilled if (a) the physical properties of the medium are uniform, and (b) terma at any location in the medium is independent of terma at any other location. The latter requirement is fulfilled for the case of internal sources because the energy released into the medium at any point is uncoupled to that released by any other point. A familiar example is found in brachytherapy where the radiation emitted by a given source does not depend on the radiation emitted by any other source in the medium. For external radiation sources, the shape of Br(i ~ it’) is a complex function of: (a) ‘beam hardening’ of polyenergetic X ray beams. (b) the number and orientation of primary rays crossing i’, and (c) exponential attenuation. As a first order approximation, beam hardening is neglected since it is an effect smaller than the attenuation problem, and the photon energy spectrum in the patient is assumed to be the same as the incident spectrum. This effect could be more accurately accounted for by representing the photon energy spectrum in the irradiated medium by a finite number of energy intervals (Eq. 5) or by an average energy, and using a convolution kernel for that spectrum.

(10)

at i’ is given by

T(i’)

The continuous digital computer

A unified approach

to optimization

If a single convolution kernel is to be used throughout the medium, then by definition the same number and orientations of primary rays will have to cross each point in the medium. For a broad beam geometry, this occurs in the region of intersection (such as target region) for either a combination of intersecting parallel-ray X ray beams, or a full rotation of a divergent-ray X ray beam. For a finite number of divergent X ray beams, each point will see a different set of ray orientations complicating the situation. It has been shown for a single beam that dose errors on the order of 1.5% may result when divergence is neglected and parallel rays are assumed in calculating central axis dose (3, 39). It is reasonable to expect a similar contribution to the overall error in dose if rayline divergence at off-axis points is neglected, and a single multidirectional kernel ( BF( i - ?‘), defined at the axis of rotation) is used for all points in the medium when calculating dose from multiple external X ray sources. Exponential attenuation has the greatest effect on altering the shape of Bt(i - it’) through the weights w( ?‘)i, and strictly speaking, the kernel is not spatially invariant for external sources. Spatial invariance implies that the fractional contribution of total terma for a given beam direction be constant for all points (e.g., ~(7’); N constant). This situation is reasonably approximated only within a small target region when centrally located in a symmetrical homogeneous medium, and irradiated by uniform intensity X ray beams. The dimensions of the region are limited to be much less than the photon mean free path of the X ray energy used (d @ 16 - 25 cm). These assumptions are typically satisfied in the very limited case of stereotactic irradiation of small volumes of several centimeters dimension, centrally located in the brain. For the general situation, BE( i - 7’) would have to be recomputed at each location i’ to properly account for distortion due to contour shape (inhomogeneities are neglected here) and beam modifiers (blocks, wedges, compensators) resulting in prohibitively long calculation times. In a recursive optimization procedure, such as described here, this limitation can be temporarily overlooked by noting that the use of a single kernel in the deconvolution step is merely a convenience that leads to an efficient estimate of incident fluence profiles (26). These profile estimates are then used in calculating dose either a) approximately, using the single convolution approximation of Eq. 13, or b) accurately, by calculating each beam separately using Eq. 9. Typically spatial invariance of B,( i - 7’) is assumed for intermediate steps to improve efficiency at the cost of accuracy, with the requirement that dose for each beam is calculated separately and summed (Eq. 9) to accurately determine the total dose distribution in the final step of the optimization algorithm. Given a ‘spatially compact’ kernel, Eq. 16 can be inverted producing what is equivalently an optimization operation. Ideally, the kernel will approximate a deltafunction in appearance in real space. This form has a

0 T. HOLMES eta/.

863

frequency-space amplitude distribution which does not have significantly small values at high frequencies (or zeroes) that would produce large oscillations in the deconvolved result. Measurement or calculational noise may be present in the kernel data resulting in zeros or small high frequency values. It may be possible to low pass filter the data or approximate it with analytical expressions amenable to Fourier deconvolution. Ahnesjij (4) has fit Monte Carlo generated kernel data to a radial function of two exponentials for the range of energies used in external beam radiotherapy. His expressions may be useful in this regard, although they have not been investigated by the authors. The definition of the convolution kernel B(7 - 7’) (and BE( i - 7’)) results in sufficiently ‘compact’, that is, steeply varying, kernels to allow Fourier deconvolution for both internal and external radiation sources. Radial profiles of the kernels used for this work are shown in Figure 1. On the other hand, the convergent point irradiation dose distribution defined by Brahme ( 19) for external source optimization did not lend itself well to Fourier deconvolution because of its broad shape. More importantly, the kernel definition forms the basis for a unified approach to dose optimization of internal and external sources by treating the energy released by these sources in a similar fashion. In the following section we will discuss the application of the model to three radiotherapeutic modalities: radioimmunotherapy, brachytherapy, and external beam

10’

ETTrrIl

c

loo 10-l

b r

10-e 9 1om3

r

1o-4

b

10-s

r ? 0

1o-6

h

lo-’

a

1om8

5 o

1o-Q

F

Fig. I. Log-log representation of planar convolution kernels used in this work. 2-D planar kernels: internal point sources of 32P (beta emitter, At = 8.0 X 10m4 cm), and ‘?r (sheathed in platinum, At = 2.5 X 10m3 cm), 3-D planar kernel: external, rotationally symmetric, line source of 60Co. Geometrical inversesquare falloff (rM2) is included for comparison.

864

I. J. Radiation

Oncology 0 Biology 0 Physics

April

199I ~Volume

radiotherapy. They will be presented in the order shown since it is the order of increasing complexity as well as clinical utility.

RESULTS

Optimization of‘internal

Source Distributions

Optimization, or the deconvolution of the kernel from the ideal dose distribution to give a desired ideal weighting distribution of terma T( i)idcal, is given by: T(

7 )ideal

=

F’{

F{

D(?)ldea,}/F{

&(i

-

i’)]

).

(17)

An ideal dose distribution (a 2-dimensional example is used here for the sake of discussion) is initially specified, typically having the ‘step-like’ form shown in Figure 2a. although there is no restriction placed on its shape. In this example, the dose within the treatment volume is assigned a value of 1000 cGy, and 0 cGy outside. This merely reflects the fact that, lacking a priori knowledge of the final form of the distribution edges, which are not ‘step-like’, one would simply assume discontinous transitions from one region to another. For the 2-dimensional internal source examples presented here, a 2-D planar kernel was obtained by taking a thin slice through the center of the three-dimensional kernel BE( i - 7’). The slice thickness, At, was chosen to be much less than ( 10%) of the minimum resolution of the 3-D data in order to minimize the partial volume error due to sampling at a coarse resolution. This kernel represents the fraction of energy absorbed per unit energy released at the decay site, per unit volume in the plane. Fourier deconvolution using an ideal dose distribution with discontinuous shape will produce negative terma (Fig. 2b). Although mathematically correct, negative terma is physically unrealizable and must be dealt with. The simplest approach is to remove the negative values by addition of positive terma at these locations in the distribution T( ?‘)idcal. The effect of zeroing negative terma on the dose distribution must be established. This is accomplished by a ‘forward’ convolution ofthe perturbed terma distribution with the convolution kernel (Eq. 16). A blurring of the ideal dose distribution results which is a direct consequence of the zeroing operation. For a beta emitter this smearing is small (Fig. 2~). The dose optimization method described above could be applied in radioimmunotherapy for the determination of the optimal distribution of source activity. Unfortunately, this modality does not allow direct creation of the optimal activity distribution in viva even though it is ideally suited, from a mathematical point of view, for this optimization technique. For brachytherapy optimization, the situation is different. Equation 14 can be modified to explicitly account for the discrete positions of the sources ( 14).

D(i)

20, Number

4

= C [T(T’).II(T’)]B,(i ;’

~ ?‘)d3T’

(18)

where II( i’) E (0, 1 ), also called the shah function ( 16 ). describes the spatial distribution of sources (Fig. 3a). 18 is a generalization of Eq. 14 because for a Equation ‘continuous’ terma distribution, such as occurs in radioimmunotherapy and external beam therapy, II is unity everywhere, implying all points in the medium can participate in the distribution ofenergy into the medium. Again. the Convolution Theorem can be used to cast the problem into frequency space where it becomes a simple product of a discretely sampled terma distribution and a convolution kernel: D(i)

= F-l{

F( [T(i’).II(i’)]; XF{B,(ip?‘))).

Likewise.

the deconvolution

[T(i’).II(i’)]

is described

(19) by:

= F-‘{F{D(i)J/ F{ Br(i

~ 7’)) I.

(20)

The discrete sampling places a severe restriction on the deconvolution process. First, the terma distribution and the sampling function are inseparable. This is a consequence of zeros in the sampling function which restrict its use as a divisor on the right side of Eq. 20. This also implies that for deconvolution to take place. the ideal dose distribution should itself be somewhat discrete in form. This is further exacerbated if the kernel is highly ‘delta-like.’ Proper specification of the ideal dose distribution would then require a priori knowledge of the source placement, II( 7’). The problem now becomes a difficult one, that of specifying a highly irregular shape for the ideal dose distribution. A way around this difficulty is to ignore the discrete sampling for the time being and perform the deconvolution assuming II = 1 for all i’. The ‘continuous’ terma distribution (Fig. 2d) can then be sampled by the actual II to give a first-order approximation to which, when convolved with the con[T( 7’) - II( i’)], volution kernel, gives the actual dose distribution (Fig. 3b). The sampled activity distribution can be scaled so that this result approximates the desired clinical dose distribution.

ldernal

source distributions

The terma T( 7) at a point T along a primary ray i in a medium irradiated by a monoenergetic source of X rays or y-rays is given by:

T(i)

= ; (?).‘I?(?)

(21)

A unified approach to optimization 0 T. HOLMES etal.

865

Fig. 2. Fourier deconvolution example for internal point sources. a.) A 2-dimensional distribution representing the ideal dose distribution. The grid size is 64 pixels X 64 pixels. b.) The ideal terma distribution resulting from the deconvolution of the P-32 2-D planar convolution kernel from the ideal dose distribution. Due to local absorption of emitted beta particles the convolution kernel is ‘spatially compact.’ The negative terma is physically not realizable and must be set to zero. The pixel size is 0.5 mm X 0.5 mm. c.) The dose distribution resulting from the convolution of the P-32 2-D planar convolution kernel with the positive component of the ideal terma distribution. Removal of negative terma perturbs the ideal solution, which results in a slight increase in dose near the edges of the dose distribution. d.) The ideal terma distribution resulting from the deconvolution of the Ir- 192 2-D planar convolution kernel from the ideal dose distribution. The gradual falloff of terma away from the distribution edges is due to a

broader kernel describing the transport of energy to greater distances away from the decay site by emitted photons. The pixel size is 2 mm X 2 mm.

and the energy fluence q(7) is related energy fluence \k,( 4, 19)by the relation

to the incident

(22)

where Qo( 4, 0) is the energy fluence along a rayline described by 6 (p - i - ? ) . C#J is the angle between the rayline and the beam axis. 8 is the gantry rotation angle. The exponential term represents the attenuation of the fluence from the source to the point at i. p/p(T) is the effective mass attenuation coefficient distribution. p( 7) is the den-

866

I. J. Radiation

Oncology

0 Biology 0 Physics

April I99

I 1Volume 20, Number 4

6.4

3.2

-32

-6.4 6.4

-3.2

0.0 CM

3.2

6.4

Fig. 3. (a) The shah function optimization.

II(i) used to describe the spatial distribution of point sources in brachytherapy The grid size is 64 pixels X 64 pixels. The pixel size is 2 mm X 2 mm. (b) The isodose distribution

(cGy) resulting from the convolution of the sampled ideal terma distribution with the Ir- 192 2-D planar convolution kernel. The grid size is 64 pixels X 64 pixels. The pixel size is 2 mm X 2 mm.

sity distribution. SSD( 4, (I) is the source-to-surface distance. i is the position vector relative to the rotation axis, i is the unit vector along 7, and p is the perpendicular distance from the axis of rotation to the rayline. These parameters are illustrated in Figure 4. The examples given in this paper assume parallel-ray geometry (SSD*( 4, fI)/ (iI’= l.O),aunitdensitymedium(p(?)= l.Og/cm”), and ignore beam hardening effects ( p/p (T ) = constant ) . Equation 22 functionally describes an attenuated backprojection along a single rayline of a divergent-ray projection. A similar situation is encountered in Emission Computed Tomography (ECT), albeit with a sign change

on the exponent to reflect the attenuation of emitted photons occurring in ECT. Figure 5 illustrates the analogy of external beam radiotherapy with image reconstruction for ECT. In nuclear medicine, an image of gamma-rays projecting through the patient is obtained by a gamma camera. Here the attenuated Radon transform (50) is physically obtained. The projection image is exponentially backprojected to obtain an image of the distribution, but before this is done a filter is applied to the projections to remove the blurring inherent in the backprojection operation. This is called exponential filtered backprojection.

A unified approach

to optimization

861

0 T. HOLMES c'f nl.

attenuated Radon transform with a sign change on the exponent to reflect exponentiation) . By Eq. 2 1, this relation is also equivalent to

a($,

0) = P

T(ije+p~a(p-i.i)di

s

S(p - P

G)di.

(26)

For a unit density medium, the energy fluence projections can be obtained directly from the result of the deconvolution operation, that is, the terma distribution. Prior to backprojection a filtering operation can be performed to remove the blurring caused by this operation (8, 30, 31). Fig. 4. The geometry of the projection / backprojection given by Eq. 22 and Eq. 24.

In radiotherapy, energy fluence is ‘backprojected’ and attenuated into the patient to a point of interaction. By summing multiple backprojections. a distribution of energy fluence in the patient is built up. Conversely, given a distribution of energy fluence within the patient. one can obtain incident energy fluence projections using the exponential Radon transform. Since a projection is a line integral across a distribution, the energy fluence distribution (Eq. 22) must first be differentiated with respect to depth along projection rays prior to integrating:

In this way dimensionality is preserved between the projected distribution and the backprojected distribution. The line integral s (d*( 7)/dt) dt, corrected for attenuation and inverse-square fall off, represents a projection operation:

J[f(T)p(?)S(p X *(7)e

- i -7) + f 1

+lalp(T)p(;)a(p~i.;)di

Iii2

SSD2(4,

x 6(p - i.i)di

0)

(24)

where the negative sign has been dropped by reversing the limits of integration to be from the source to the point T. For a parallel ray geometry, assuming a uniform medium and no beam hardening, Eq. 24 reduces to

This is a scaled exponential

Radon

transform

(that is, the

0) = F-‘{F{*(&,

e)ilwl)+

(27)

where \k,( &tI) is an estimate of the ideal incident energy fluence, )w 1 is the frequency space filter function (a linear ramp), and the ‘+’ subscript indicates the positive component of the filtering result. The projections can be converted to photon fluence, a,,($, B), by dividing by the photon energy. This is not necessary in practice since the profile data can be normalized to provide relative transmission, and the absolute fluence accounted for in the timer setting calculation for each beam. The optimization process for external beam radiotherapy consists of the steps outlined below for the example of a 360” rotation of a 60Co radiation source and a ‘W’ shaped dose distribution in an elliptical unit density medium (Fig. 6a). The rotation is approximated by 72 equally spaced, isocentric, coplanar, parallel-ray radiation beams. An isotropic kernel was used under the assumption that the fractional contribution of terma is the same for each beam direction for all points in the irradiated volume. The feasibility of the algorithm is demonstrated by optimizing the dose in a plane perpendicular to the rotation axis and centered in the irradiated volume. The dose in this plane is determined by transporting the primary photons into the irradiated volume followed by a three-dimensional convolution over the volume. An approximation to the above is given by a 2-D convolution over the plane using a 2-D planar kernel that is obtained by integrating the 3-D volume kernel into a plane perpendicular to the rotation axis. The planar kernel implicitly includes out-of-plane scatter to a distance of 60 cm on either side of the plane of calculation, so it is assumed that the radiation beams (including the in-plane primary profiles), external contour, and target volume extend semi-infinitely. We have used the latter approach in this work. Deconvolution mathematics is used to arrive at the ideal terma distribution in the patient. This may contain negative values. The ideal terma distribution is projected (Eq. 26) at desired angles and scaled by the mass density of the medium to give a set of ideal, unfiltered energy fluence projections (Fig. 7a). The projections must then be filtered

(23)

d*(i) _=_ dt

WA 0) =

*,(d,

operations

868

I. .I. Radiation

Oncology

ECT lmaae

0 Biology 0 Physics

Reconstructlon

199

I, Volume 20, Number

External

A ECT Image aqulsltlon Is physlcally obtalneci projectlons eouwalent to the attenuated Radon transform

C ECT Image reconstructlon eaulvalent to the mathematical operation of exponential backproject Ion of f 1I tered projectlons

April

Is

Source

4

Optlmlzatlon

B External beam optlmizatlon mathematlcally obtalned pro]ectlons by the exponential Radon transform.

is

0 Optlmlzed external beam therapy 1s physrcally eoulvalent to the posItwe values of attenuated flltered backprojections

Fig. 5. Diagram illustrating the analogy of external beam therapy (XT) with emission computed tomography (ECT). The physical acquisition of ECT image profiles (a) by the attenuated Radon transform of the activity distribution is analogous to energy fluence profiles obtained mathematically by the exponential Radon transform of the terma distribution in XT (b). Frequency filtering in both situations results in negative values in the projection data. These values are not physically realizable m XT and must be set to zero before backprojecting. The mathematical operation of image reconstruction in ECT (c) by backprojection of filtered projections has rotational XT (d) as its physical analog.

to remove the blurring artifact that will be caused by the backprojection operation (Fig. 7b). Since the projections represent a physical quantity, any negative values must be set to zero (Fig. 7~). The positive component of the filtered projections are then attenuation-backprojected into the patient (Eq. 22) to get the perturbed energy fluence distribution and this result converted to terma (Eq. 21). The zeroing operation no longer guarantees the initial assumption that each beam contributes equal terma to a given point. Because of this perturbation, dose distributions for each beam should be calculated separately and summed to give exact results for the total dose distribution

(Eq. 9). As this algorithm can be very time-consuming for a large number of beams, Eq. 13 can be used to obtain an approximation to the actual dose distribution, assuming that each beam still contributes equal amounts of terma at points T’ in the patient (Fig. 6b). The error in dose resulting from this assumption may be clinically significant and requires further investigation. The zeroing operation causes an increase of dose within the boundary of the dose distribution and a blurring of its edges. The dose distribution can be optimized further by a method which iteratively reduces the excess dose to the optimal level initially specified. A ‘residual dose’ distribution can be obtained (Fig. 8a) which is the dose that

16

IG

-a

- 16 0

A

16

CM 16

8

If

-8

0

8

16

I6

0

-a

CM

a

IO

CM

Fig. 6. External beam dose optimization example. (a) The ideal dose distribution consisting of a target volume of 1000 cGy and 0 cGy outside of the target. The optimum dose distribution (in units of cGy) after (b) one iteration, (c) two iterations, and (d) ten iterations. A ‘(‘Co 3-D planar convolution kernel was used in the deconvolution and 72 equally spaced projections were used in the projection/backprojection steps. The variation of dose in the target volume in Figure d is 1-2 percent about 1000 cGy. In each figure the bold isodose line represents the 980 cGy isocontour, and the set of isodoses shown is 400, 600, 800,980. 1100, 1200, 1300, 1400, and 1500 cGy. The pixel size is 5 mm X 5 mm.

0 CM

12

23

12

23

5

4

-3

t

/

:“j 0 *

2

1

-12

0 CM

12

23

0 -23

-12 c”M

Fig. 7. Projections at 45” clockwise from vertical illustrating the various steps in the optimization method described in the text. (a) The initial projection of energy fluence. (b) The projction after filtering with a frequency ramp filter. (c) The projection after zeroing negative energy fluence. (d) The projection after removal of the first residual projection. (e) The projection after removal of the second residual projection. ( f) The projection after removal of the ninth residual projection. Note that most of the excess energy fluence is removed after the first iteration (Fig. d). The small amount of energy fluence outside of the beam boundary in Figures d-f is a result of not zeroing negative fluence in the residual projections prior to subtracting them from the main projection data. The pixel size is 5 mm X 5 mm.

870

I. J. Radiation

Oncology 0 Biology 0 Physics

April

199I.

Volume

20. Number

4

Fig. 8. The reduction in the residual dose for increasing number of iterations. The pixel size is 5 mm X 5 mm. (a) The maximum residual dose after the first iteration is 586 cGy. (b) the second iteration is 103 cGy, (c) the third iteration is 39 cGy, and (d) the tenth iteration is 19 cGy.

exceeds the upper and lower dose constraints. The idea1 dose distribution (Fig. 6a) specifies the lower dose constraint in this example. The upper dose constraint was set at a constant value of 1000 cGy. Frequency filtered projections can be found for the residual dose distribution using the same techniques mentioned previously, although negative energy fluence need not be zeroed for these projections. These data are subtracted from the previous projection data and negative values set to zero (Fig. 7d). The result is attenuation-backprojected to give the corrected dose distribution (Fig. 6~). In this way the ‘residual dose’ can be reduced rapidly after only several iterations (Fig. 8d). In this example, the first iteration is very effective at removing most of the residual dose, as demonstrated by the change in the projection data from Figures 7c to 7d. Subsequent iterations have a much smaller effect on a given projection, as seen in the small changes between Figures 7d, 7e, and 7f. In Figure 9, a set of four pairs of diametrically opposed projections are shown. These were drawn from the set of 72 projections produced after three iterations of the algorithm. The opposed projections are slightly dissimilar because of differences in attenuation corrections for the non-symmetric target dose distribution, and to rounding errors introduced in the calculation. The

magnitude of the variation of energy fluence within the projections is roughly a factor of 10. For fixed field radiation therapy, the variations present in these projections can be achieved using 2-3 HVL (2-3 cm) of heavy metal. Given the large number of beams used in this example, modulating each beam using a fixed metal filter (dose compensator) is impractical in the clinical setting. Either a multileaf collimator or a scanned photon beam would be necessary to achieve treatment times consistent with conventional radiotherapy techniques ( 17). Since these technologies are now becoming commercially available, it is reasonable to expect that within the next few years a large number of radiotherapy departments will have them in place, and the optimization algorithm presented here may find practical application.

DISCUSSION In this work we present a numerical approach to optimizing radiation dose distributions that uses genera1 image filtering and reconstruction techniques. A key feature of our approach is our definition of the convolution kernel so that Fourier deconvolution can be used instead of iterative methods ( 34). This operation is used in a unified

A unified

3000

approach

to optimization

0 T. HOLMES

871

et ul.

A

2000

1000 0 -23

-12

C"M

12

j

loooi__s{

23

0 -23

4000

_i -L2

II CM

12

23

-12

0 CM

12

23

0

12

4000 y

0 -23

-12

0 CM

L2

23

0 -23

-12

0 ,-l‘

12

23

0 -23

0

Fig. 9. Selected projections from the set of 72 projections (e) 135, and ( f) 3 15 degrees clockwise from vertical.

way to determine an ideal terma distribution for either internal or external sources. For external sources the terma distribution is transformed into projection estimates by exponential Radon transformation and projection filtering. The method is demonstrated for several 2-dimensional examples using a number of simplifying assumptions. The optimization of continuously distributed internal point sources by Fourier deconvolution is reasonably straightforward. In a clinical application such as radioimmunotherapy, the in-vivo distribution of sources is governed by biochemical factors not under the control of the physician, and the method is not readily applicable. The optimization of discretely distributed sources is more difficult. The method presented only provides a first order estimate of the distribution of source activities in brachytherapy, and it requires further refinement. Fourier deconvolution (coupled with filtered backprojection) is most useful for external beam radiotherapy optimization given that the approriate assumptions are fulfilled. A critical assumption is that the convolution kernel is spatially invariant. This is not strictly true for multidirectional kernels. The error resulting from dose computation using these kernels and Eq. 13 needs further clarification. Our initial investigation of this error, using

-23

-12

of the third iteration:

ru

23

(a) zero. (b) 180, (c) 90, (d) 270,

the example in this work, indicates it is approximately I5%, where Eq. 13 underestimates the dose by this amount (26). On the other hand. the use of multidirectional kernels in providing estimates of incident energy fluence profiles is certainly acceptable if the dose is subsequently determined by calculating each beam separately using Eq. 9 and a single beam kernel. Work is underway to extend the model to 3-dimensions by using 3-D Fourier deconvolution, and for external sources, 3-dimensional image reconstruction techniques ( 8,25 ). Further effort is required to include optimization of the number and location of radiation sources since they are currently specified beforehand. Additionally, the beam hardening effect needs to be accounted for in the dose calculation for external sources. Another important issue is the choice of dose constraints. In the simple example presented for external beam optimization we chose only to optimize the dose uniformity in the target region. The addition of other constraints such as minimal dose to critical structures and normal tissue volumes is necessary and will be investigated. The methods presented in this paper are general enough in their features that algorithms written originally for 2dimensional image filtering and for the ECT image re-

872

I. J. Radiation Oncology 0 Biology 0 Physics

construction problem (30) were used. The coupling of convolution kernels generated by first-principle Monte Carlo methods36 with optimization using convolution and image reconstruction mathematics provides a very powerful calculation platform for 3-dimensional treatment

April 199

I 1Volume 20. Number 4

planning. These algorithms will place a high demand on computer hardware. Fortunately, computer hardware price/performance ratios continue to improve, reducing calculation times and enabling 3-dimensional treatment optimization to become clinically feasible ( 37 ).

REFERENCES 1. Ahnesjo, A. Application of transform algorithms for calculation of absorbed dose in photon beams. In: Proceedings of the 8th Int. Conf. in Radiotherapy, Toronto, Canada. Silver Springs. MD: IEEE Computer Society; 1984: 17-20. 2. Ahnesjii, A.; Andreo, P.: Brahme, A. Calculations and applications of point spread functions for treatment planning with high energy photon beams. Acta Oncol. 26:49-56; 1987. 3. Ahnesjii, A. Invariance of convolution kernels applied to photon beams. In: Proceedings of the 9th Int. Conf. on Computers in Radiotherapy. Den Haag, Amsterdam. Amsterdam: Elsevier Science Publishers; 1987:99. 4. Ahnesjo, A. Collapsed cone convolution of radiant energy for photon dose calculation in heterogeneous media. Med. Phys. 16(4):577-592; 1989. 5. Altschuler, M. D.: Censor, Y. Feasibility solutions in radiation therapy treatment planning. In: Proceedings of the 8th Int. Conf. in Radiotherapy, Toronto, Canada. Silver Springs, MD: IEEE Computer Society; 1984:220-224. 6. Anderson, L. L. Physical optimization of afterloading techniques. Strahlentherapie I6 I :264: 1985. 7. Attix, F. H. Introduction to radiological physics and radiation dosimetry. New York: John Wiley and Sons: 1986. 8. Barrett. H. H.; Swindell. W. Radiological imaging, Vol. 2. New York: Academic Press; 198 I. 9. Barth, N. H. An inverse problem in radiation therapy. Int. J. Radiat. Oncol. Biol. Phys. 18:425-431; 1990. IO. Boyer, A. L. Shortening the calculation time of photon dose distributions in an inhomogeneous medium. Med. Phys. I I (4):552-554: 1984. 1I Boyer. A. L.; Mok, E. C. Photon beam modeling using Fourier transform techniques. In: Proceedings of the 8th Int. Conf. in Radiotherapy, Toronto, Canada. Silver Springs. MD: IEEE Computer Society; 1984: l4- 16. 12. Boyer, A. L.; Mok. E. C. A photon dose distribution model employing convolution calculations. Med. Phys. 12( 2): I69177; 1985. 13. Boyer. A. L.: Mok, E. C. Calculation of photon dose distributions in an inhomogeneous medium using convolutions. Med. Phys. I3 (4):503-509; 1986. 14. Boyer, A. L.; Mok, E. C. Brachytherapy seed dose distribution calculation employing the fast Fourier transform. Med. Phys. I3( 4):525-529: 1986. 15. Boyer, A. L.; Zhu, Y.; Wang, L.; Francois, P. Fast Fourier transform convolutions calculations of x-ray isodose distributions in homogeneous media. Med. Phys. 16(2):248253; 1989. 16. Bracewell. R. N. The Fourier transform and its applications, 2nd edition. New York: McGraw-Hill; 1978. 17. Brahme, A. Design principles and clinical possibilities with a new generation of radiation therapy equipment. Acta Oncol. 26(6):403-412; 1987. 18. Brahme, A. Optimization of conformation and moving beam radiation therapy techniques. In: Proceedings of the 9th Int. Conf. on Computers in Radiotherapy, Den Haag, Amsterdam. Amsterdam: Elsevier Science Publishers: 1987: 227-230. 19. Brahme, A. Optimization of stationary and moving beam

20.

2 I.

22. 23.

24.

25. 26.

27.

28. 29.

30.

3 1. 32.

33. 34.

35. 36.

37.

38.

radiation therapy techniques. Radiother. Oncol. 12: I29140: 1988. Brahme. A.: Roos, J. E.; Lax, I. Solution on an integral equation encountered in rotation therapy. Phys. Med. Biol. 27:1221-1229; 1982. Burns. G. S.; Raeside, D. E. Two-dimensional dose distribution around a commercial ‘25I seed. Med. Phys. 15 ( 1): 56-60; 1988. Cormack, A. A problem in rotation therapy with xrays. Int. J. Radiat. Oncol. Biol. Phys. 13:623-630; 1987. Cormack, A.; Cormack, R. A problem in rotation therapy with xrays: dose distributions with an axis of symmetry. Int. J. Radiat. Oncol. Biol. Phys. 13:1921-1925; 1987. Cunningham, J. R. Keynote address: development of computer algorithms for radiation treatment planning. In. J. Radiat. Oncol. Biol. Phys.. 16(6):1367-1376; 1989. Deans, S. R. The radon transform and its applications. New York: John Wiley and Sons: 1983. Holmes, T. W.: Mackie, T. R. A study of the validity of modelling multiple x ray beams by a single convolution (Abstr.). Med. Phys. 17(3):533: 1990. Hodes. L. Semiautomatic optimization of external beam radiation treatment planning. Radiology I IO: 19 l-196: 1974. Hope. C. S.; Orr, J. S. Computer optimization of 4 MeV treatment planning. Phys. Med. Biol. lO( 3):365-373; 1965. Hope, C. S.: Laurie, J.; Orr. J. S.; Halnan. K. E. Optimization of X-ray treatment planning by computer judgement. Phys. Med. BIoI. 12(4):531-542; 1967. Huesman. R. H.: Gullberg, G. T.; Greenberg, W. L.; Budinger. T. F. Users manual: Donner algorithms for reconstruction tomography, LBL pub-214. Berkeley, CA: Lawrence Berkeley Laboratory, October 1977. Kak, A. C.; Slaney, M. Principles of computerized tomographic imaging. New York: IEEE Press; 1988. Langer, M.; Leong, J. Optimization of beam weights under dose-volume restrictions. Int. J. Radiat. Oncol. Biol. Phys., 13:1255-1260; 1987. Lax. I.; Brahme, A. Rotation therapy using a novel highgradient hlter. Radiology 145:473-478; Nov. 1982. Lind, B. K.; Brahme. A. Optimization of radiation therapy dose distributions using scanned electron and photon beams and multileaf collimators. In: Proceedings of the 9th Int. Conf. on Computers in Radiotherapy, Den Haag, Amsterdam. Amsterdam: Elsevier Science Publishers; 1987:235239. Mackie, T. R. Ph.D. Thesis. Univ. of Alberta, Edmonton, Canada: 1984. Mackie, T. R.; Bielajew, A. F.; Rogers, D. W. 0.; Battista, J. J. Generation of photon energy deposition kernels using the EGS Monte Carlo code. Phys. Med. Biol. 33( I ): I-20; 1988. Mackie, T. R.; Kubsad, S.; Holmes, T.; Sohn, W. New developments in radiotherapy dose planning. In: Proceedings of the Int. Congress of Radiology, Paris, France, 1989. Mackie, T. R.; Scrimger, J. W. Computing radiation dose for high energy x rays using a convolution method. In: Proceedings of the 8th Int. Conf. in Radiotherapy, Toronto,

A unified

39.

40.

41.

42. 43.

44.

45.

approach

to optimization

Canada. Silver Springs, MD: IEEE Computer Society; 1984: 36-40. Mackie, T. R.; Scrimger, J. W.; Battista, J. J. A convolution method of calculating dose for 15-MV x rays. Med. Phys. 12(2):188-196; 1985. Mackie, T. R.; Reckwerdt, P. J.; Holmes, T. W.; Kubsad, S. S. Review of convolution/superposition methods for photon beam dose computation. In: Proceedings of the 10th Int. Conf. on the use of Computers in Radiation Therapy, Lucknow, India, Nov. I 1-14, 1990. McDonald, S. C.; Rubin, P. Optimization of external beam radiation therapy. Int. J. Radiat. Oncol. Biol. Phys., 2:307317; 1977. Mohan, R.; Chui, C. Differential pencil beam dose computation model for photons. Med. Phys. 13( 1):64-73; 1986. Mohan, R.; Chui. C. Use of fast Fourier transforms in calculating dose distributions for irregularly shaped fields for three-dimensional treatment planning. Med. Phys. 14( 1): 70-77; 1987. Redpath, A. T.; Vickery, B. L.; Wright, D. H. A new technique for radiotherapy planning using quadratic programming. Phys. Med. Biol. 21(5):781-791; 1976. Rosenstein, L. M. A simple computer program for optimization of source loading in cervical applicators. Brit. J. Radiol. 50: 119; 1977.

0 T. HOLMES el al.

873

46. Simpkin, D.; Mackie, T. R. EGS4 Monte Carlo determination of the beta dose kernel in water. Med. Phys. 17( 2): 179-186; 1990. 47. Starkshall, G. A constrained least-squares optimization method for external beam radiation therapy treatment planning. Med. Phys. 11(5):659-665; 1984. 48. Thomason, C.; Higgins, P. Radial dose distribution of 19*Ir and 13’Cs seed sources. Med. Phys. 16(2):254-256; 1989. 49. Tolli, H.; Ragnhult, I. An analytical method to solve an optimization problem in brachytherapy. In: Proceedings of the 9th Int. Conf. on Computers in Radiotherapy, Den Haag, Amsterdam. Amsterdam: Elsevier Science Publishers; 1987: 123-126. 50. Tretiak, 0.; Metz, C. The exponential radon transform. SIAM J. Appl. Math. 39(2):341-354; 1980. 51. Weaver, K.; Smith, V.; Huang, D.; Barnett, C.; Schell, M.; Ling, C. Dose parameters of “‘1 and ‘921rseed sources. Med. Phys. 16(4):636-643; 1989. 52. Webb, S. Optimization tributions by simulated 134991370; 1989.

of conformal radiotherapy dose disannealing. Phys. Med. Biol. 34( 10):

53. Zhu, Y.; Boyer, A. X-ray dose computations in heterogeneous media using 3-dimensional FIT convolution. Phys. Med. Biol. 35(3):351-368; 1990.