MINIMUM SUPPORT SOLUTIONS FOR
RADIOTHERAPY TREATMENT PLANNING
by
Janine Marie Kennedy
B.A., Fairfield University, 1993
A thesis submitted to the
University of Colorado at Denver
in partial fulfillment
of the requirements for the degree of
Master of Science
Applied Mathematics
2000
Â£
\
This thesis for the Master of Science
degree by
Janine Marie Kennedy
has been approved
by
Weldon Lodwick
Date
Kennedy, Janine Marie (M.S.?Applied Mathematics)
Minimum Support Solutions for Radiotherapy Treatment Planning
Thesis directed by Professor Stephen C. Billups
ABSTRACT
A radiotherapy treatment session typically lasts for approximately
fifteen minutes. This time constraint makes the implementation of a radiation
treatment plan with more then ten beams impractical for the average clinic.
Therefore it is appropriate and beneficial to find solutions to the treatment
planning problem with as few radiation beams as possible. For small test
problems a mixed integer program can be solved to find a minimum support
solution, which delivers the prescribed dose with less then ten beams. The
combinatorial nature of the mixed integer programming formulation makes this
method impossible for real sized problems; hence we use a fast approximate
method developed by Mangasarian to solve a series of linear programming
subproblems which produce a solution with the minimum support property.
This abstract accurately represents the content of the candidates thesis. I
recommend its publication.
Signed*
Stephen C. Billups
in
DEDICATION
I dedicate this thesis to my husband John for his unwaivering patience,
understanding and support, and to my parents for their ever present
encouragement.
ACKNOWLEDGMENTS
My thanks to Karen Kafadar, Weldon Lodwick and Francis Newman for their
help and support with this thesis and to my advisor Stephen Billups for his
patience and encouragement.
CONTENTS
Figures ........................................................... viii
Tables................................................................ x
Chapter
1. Introduction....................................................... 1
1.1 Notation......................................................... 4
2. Background......................................................... 5
2.1 Radiotherapy Treatment Planning.................................. 5
2.2 Linear Programming Formulation................................... 8
3. Finding Minimum Support Solutions................................. 11
3.1 Penalty Problem................................................. 11
3.2 The Exponential Approximation................................... 14
3.2.1 The Successive Linearization Algorithm........................ 19
4. Implementation.................................................... 21
4.1 Creating the Pixel Sets........................................ 21
4.2 Calculating the Dose Deposition Operator........................ 22
4.2.1 Finding pencil/pixel intersections............................ 23
4.2.2 Calculating the dose.......................................... 26
4.3 The fj, File.................................................... 26
4.4 Choosing a and ................................................. 27
4.4.1 Choosing a.................................................... 27
vi
4.4.2 Choosing /?................................................... 30
4.4.3 The Augmented Successive Linearization
Algorithm.................................................... 31
5. Results.......................................................... 33
5.1 Solutions for Small Problems .................................. 34
5.2 Solutions for Large Problems .................................. 39
5.3 Conclusions.................................................... 40
Appendix
A. MATLAB Subroutines............................................... 55
A.l Create Set and // Files from Patient Image..................... 55
A.2 Writes /j, Text File to Disk................................... 57
A.3 Creates AMPL Sets.............................................. 58
A. 3.1 Writes AMPL Set File to Disk ............................ 60
B. AMPL Code....................................................... 61
B.l Linear Programming Model....................................... 61
B.2 Mixed Integer Programming Model................................ 63
B.3 ASLA Script Files ............................................. 65
B.3.1 ASLA Model.................................................... 65
B.3.2 SLA Script File............................................... 68
B. 3.3 The Parameter Selection Script.............................. 70
C. Dose Depostion Operator C++ Code.............................. 73
References.......................................................... 90
Vll
FIGURES
Figure
2.1 Radiation Beams.............................................. 5
2.2 Radiation Beams with Pencils................................. 7
3.1 Comparison of Exponential Approximation 1 s~az to the Star
Norm \z\*.................................................... 14
4.1 Divergence Angle............................................ 22
4.2 Dose Deposition Model....................................... 24
4.3 g(a, fl) := a/3e~aZi, for fixed and a > 0................... 28
5.1 Comparison of Solution Methods.............................. 35
5.2 RTPP1....................................................... 41
5.3 RTPP2....................................................... 42
5.4 RTPP3....................................................... 43
5.5 RTPP3 LP solution......................................... 44
5.6 RTPP3 ASLA solution....................................... 45
5.7 RTPP4....................................................... 46
5.8 RTPP4 LP solution......................................... 47
5.9 RTPP4 ASLA solution....................................... 48
5.10 RTPP5..................................................... 49
5.11 RTPP5 LP solution......................................... 50
5.12 RTPP5 ASLA solution....................................... 51
viii
5.13 RTPP6......................................................... 52
5.14 RTPP6 LP solution........................................... 53
5.15 RTPP6 ASLA solution......................................... 54
IX
TABLES
Table
5.1 Results for RTPP1 and RTPP2................................... 38
5.2 Results for Large Problems: 64 x 64 pixels.................... 39
x
1. Introduction
External beam radiotherapy is one method currently used to treat
various types of cancer. Ionizing radiation, delivered via a machine called a
linear accelerator or betatron, damages the genetic makeup of the targeted
cells and prevents diseased cells from spreading throughout the body. Typical
treatment sessions consist of depositing varying levels of radiation to the body
along a variety of different paths called beams. Before a patient undergoes
external beam radiotherapy, a radiotherapy treatment plan (RTP) is created.
The RTP specifies the positions of the linear accelerator, which determine the
paths that radiation travels through the body, and also specifies the radiation
levels for each path.
Definition 1.1 (RTP Conditions) If radiotherapy treatment is to be suc
cessful then the RTP must be carefully and consistently developed. In this
paper a good RTP will satisfy the following basic set of conditions:
a therapeutic dose must be delivered to all cancerous tissue,
the dose distributed to healthy tissue must be kept below some fixed
maximum level,
critical structures, e.g. the kidneys or the spinal cord, must be spared to
avoid potentially devastating side effects like organ failure or paralysis,
and finally,
hospital staff must be able to administer the plan during a short treatment
session, typically 15 minutes.
The RTP conditions can be met by formulating the radiotherapy treatment
1
planning problem (RTPP) as a linear program. Newman [7] describes several
related linear programming formulations that produce solutions to the RTPP,
for instance
minimize maximum dose to critical structures
subject to tumor dose > required tumor dose
normal tissue < dose bound for normal tissue
dose > 0,
or
maximize minimum dose to tumor
subject to critical dose < dose bound for critical tissue
normal tissue < dose bound for normal tissue
dose > 0.
Either of these examples can be formulated as a linear program. Unfortu
nately, the resulting treatment plans to these formulations can involve too
many beams. In order to meet the hospitals time requirement, each RTP
must consist of as few beams as possible, with ten beams considered an up
per bound since the linear accelerator can usually be repositioned at most ten
times in a fifteen minute session.
This thesis is concerned with finding solutions to the RTPP that have
the smallest number of beams turned on possible and that satisfy all radiation
levels specified by the oncologist. Solutions to the RTPP that have a minimum
number of beams turned on will be referred to as 0 minimum support so
lutions. Such plans can be created by adding a penalty term to the objective
2
function of a linear programming formulation. The resulting math program
has a discontinuous objective function that can be modeled with a mixed inte
ger program. Unfortunately such formulations are computationally expensive;
this limits the size of problem that can be solved by the mixed integer pro
gram, making this approach impractical in terms of generating clinically viable
solutions.
This difficulty, inherent in solving even moderately sized mixed integer
programming problems, suggests that alternative methods of finding minimum
support solutions for the RTPP should be considered. In this thesis, we develop
such an alternative method. The method is based on a successive linearization
algorithm (SLA) proposed by Mangasarian [6]. This iterative method solves a
series of linear programs and produces an RTP with fewer beams than the linear
programs developed by Newman, i.e an approximate (5 minimum support
solution.
This thesis first describes the RTPP model and then presents a de
tailed explanation of a specific linear programming formulation of the RTPP.
We then discuss how Mangasarians work can be applied to find /? minimum
support solutions to the linear programming formulation of the RTPP. Sec
tion 4 describes how we modeled the RTPP process in C++ and MATLAB,
and addresses our implementation of an augmented successive linearization al
gorithm (ALSA), which we expressed in AMPL. Chapter 5 presents the results
of this implementation and conclusions regarding its effectiveness.
3
1.1 Notation
This notation will be used throughout the thesis. Let R" denote the
positive ndimensional real space. All vectors will be column vectors unless
transposed to a row vector by a prime superscript '. For a vector x in i?",
the star norm rc* is a vector of zeros and ones, with ones corresponding to
positive components i = 1,..., n of x. The scalar product of two vectors x
and y will be denoted by x'y. For a linear program min^g^ c'x the notation
{arg vertex ramxeXdx\ will denote the set of vertex solutions of the linear
program. The vector e will denote a column vector of ones, and e will denote
the base of the natural logarithm.
4
2.
Background
2.1 Radiotherapy Treatment Planning
Figure 2.1. Radiation Beams
We begin our discussion by describing how the radiation energy and
the patients body are modeled. Information about the patients body is con
tained in a computed tomography (CT) image, from which technicians can
identify locations of various organs and any tumorous growth. Using image
processing software, a technician will examine this CT image and identify or
gan and tumor boundaries. Figure 2.1 illustrates that the selected organs and
tumors are comprised of smaller pieces of image data or pixels. The technicians
selection process classifies individual pixels into one of the following subsets:
critical, the set of healthy tissue pixels for sensitive structures, such as the
spinal column or kidneys; body, the set of less sensitive healthy pixels such as
bones or muscles; tumor, the set of cancerous pixels.
In order to determine the exact amount of radiation deposited to a
specific pixel in the model, we must also have some quantitative method for
describing the rays of radiation that are emitted from the linear accelerator.
Figure 2.1 contains the origins or source points for two radiation beams. Once
the radiation leaves the linear accelerator, it diverges into a beam much like
light emitted from a flashlight. The beam passes through the image and de
posits radiation to all the pixels in its path. In Figure 2.1 the two beams
irradiate the area of the grey oval that lies between the beam boundaries. The
position of the linear accelerator is described by the angle the source point
makes with the positive x axis; for example, Source 1 is at 90 degrees and
Source 2 is at 180 degrees. In this model the accelerator has the capability of
moving about the large circle in Figure 2.1 in one degree increments.
We would like to have the capability of turning on smaller portions
of the beam, since this can be accomplished at the treatment clinic through
the use of lead or alloy blocks. To incorporate this refinement into the model
we break the beam up into evenly spaced strips called pencils, see Figure 2.2.
These pencils allow us to focus the radiation on the tumor and avoid irradiating
unnecessary sections of the body. For example in Figure 2.2, pencils PI and
P2 distribute radiation to the tumor, whereas P3 and P4 do not. The RTP
specifies which pencils will be used as well as their respective intensity levels.
6
Figure 2.2. Radiation Beams with Pencils
We can represent an RTP by a two dimensional array x, whose entries
specify radiation intensity levels for each pencil of every possible beam. Let 6
index the beam set B and p index the pencil set P. If x(p, b) = 0 then the pth
pencil of the 6th beam is off. If x(p,b) d > 0 then the pth pencil of the 6th
beam has been selected as part of the radiation treatment plan with intensity
d. If an entire column x(,b) is zero, then the 6th beam is not used. Generally,
it takes time to set up the gantry for each beam, but it does not matter how
many pencils are selected within each beam. Thus, a /? minimum support
solution is one which uses a minimum number of beams, whereas the number
of pencils is irrelevant. A dose deposition operator D is calculated for each CT
t
image, where D(i,p, b) specifies the fraction of x(p, b) that will be delivered to
pixel i. Thus if x(p, b) d and D(i,p, b) = 6, then the radiation dose delivered
to pixel i by the pth pencil in beam b will be D(i,p, b)x(p, b) = 5d.
2.2 Linear Programming Formulation
We now describe how the RTP Conditions specified in the introduc
tion can be modeled mathematically. To begin, consider the tumor require
ments. Clearly the tumor must receive a therapeutic dose of radiation. Let Tl
be the lowest radiation dose needed for the therapeutic process. If too high
a dose is delivered to the tumor necrosis or tissue death may occur. When
necrosis happens quickly the patient may have serious difficulty excreting the
dead tissue. Therefore we include an upper radiation bound for each tumor
pixel. Let Tu be the maximum tumor radiation specified by the oncologist.
Then for each pixel t tumor we have
T'SEE &)*(P,) < T*. (2.1)
pePbeB
Equation (2.1) states that the total sum of radiation distributed to a tumor
pixel by each pencil in the plan must be more than the minimum bound Tl and
less than the maximum bound Tu prescribed by the oncologist. We use the
same method to put an upper bound on the body pixels. Let Ov be the upper
radiation bound for noncritical healthy cells, then for each pixel o E body we
have
X X D(' 6)*(Pi b)
pePbeB
In Section 1 we stated that critical structures like the spinal column need to
be spared as much as possible. We therefore seek to minimize the maximum
8
quantity of radiation that is distributed to a critical pixel. Then for each pixel
c 6 critical, let
7(c) = EE Â£(c>p> 6)x(p> 6) (23)
pePbeB
which represents the total radiation delivered to pixel c.
One approach to creating treatment plans is to solve the following
mathematical program first proposed by Newman [7, p.10]
mina; maxc 7 (c)
subject to Tl < EPepEbeBD(t,P,b)x{p,b)
Ep(o,p, 6)x(p, 6) < o <= body.
This math programming model produces uniform treatment plans that satisfy
all radiation bounds. The objective function minimizes the maximum dose
deposited to any pixel in any critical structure. This is superior to minimizing
the total radiation dose for all critical structures because it eliminates hot
spots, which are areas that receive very high doses of radiation. We avoid hot
spots for two reasons: first, hot spots can be necrotic and therefore difficult for
the patient to excrete; second, if the tumor is very small and the patient is not
set up perfectly, the hot spot could be shifted to normal tissue, which can be
a very serious problem.
Using a standard technique for piecewise linear and convex functions
from [1, p.16] one can change (2.4) into a linear program. Since 7(01), 7(02),...
are convex functions
7 = max 7(c)
is a piecewise linear convex function [1, Theorem 1.1]. By adding a constraint
9
for each pixel c G critical we can replace the piecewise linear convex objec
tive function maxc7(c) with the linear objective function 7 and produce the
following equivalent linear programming formulation:
min 7
subject to 7  Zppl2bÂ£BD(c,P,b)x(p,b) >0
Tl < EppEbeBD(t,p,b)x(p,b)
Epp Ebefl D(o,p, b)x(p, b) < Ou
For simplicity, let the convex set S be defined as the feasible region of linear
program (2.5). In this formulation, 7 represents the maximum radiation deliv
ered to any critical pixels; other related linear programming formulations are
discussed in detail in [7, p.711]. The only apparent difficulty with the above
linear programming formulation is that its solutions do not necessarily meet
the hospitals time constraints, In particular, the optimal solution may specify
too many different beam angles to allow for a timely treatment session. In the
following chapter, we will describe a method for satisfying this last criterion
by applying the SLA developed in [6] to the linear programming formulation
of the RTPP.
c G critical
(2.5)
t G tumor
0 G body.
10
3. Finding Minimum Support Solutions
3.1 Penalty Problem
Recall that in order to satisfy a hospitals time requirement, we need
to find RTPs that have the fewest number of beams possible. This can be
accomplished by modifying the objective function in (2.5) to penalize each
beam that is turned on. The penalty problem is therefore
7 + /3e'..,, /5 > 0,
(3.1)
where
T= {(7. *.2)1(7.*) e S,z = 53a:(p,)},
p
S is the feasible region of (2.5), /3 is a positive scalar, which regulates the
importance of beam reduction, and z* represents the starnorm of z, which is
a vector with components
(M*)i
0 if Zi = 0
<
(3.2)
1 otherwise.
In this model, the components of z represent the total radiation level of each
beam. The penalty parameter /? is chosen to balance the conflicting goals of
minimizing critical structure dose, while simultaneously keeping the number of
beams as small as possible. If (7^,x^,z^) is a solution to (3.1), we shall call
(7^,**) a P minimum support solution to the RTPP.
We now examine how solutions of (3.1) compare to solutions of (2.5).
When we solve a feasible instance of (2.5) we get an optimal solution (7*, a:*);
11
when we solve (3.1) for a particular /3 we get a solution (7^,0^, z@) for which
7^ can be no better than 7*, i.e. 7^ > 7*. The next result gives a bound on
the difference between 7^ and 7*.
Theorem 3.1 (Error Bounds for the Objective Value of the Penalized
Linear Program) If the RTPP (2.5) has solution (7*, re*) and the Penalized
Linear Program (3.1) has solution (7then
 7* < /%*l* < Pn,
where
**:=5>*0v).
p
and n is the total number of beams in the model.
Proof: By construction, (7*,x*,z*) is feasible for (3.1); hence in the model
since (7^, x&, z@) is optimal for (3.1),
7^ + P\z^\* < 7* + P\z*\*.
Rearranging terms yields the inequality
7^ 7* < P(\z*\* ~ 1^1*)
Furthermore, in order to satisfy the tumor dosage bounds set forth in the
RTPP, \z% 0. Hence 7^ 7* < P{\z% \z%) < p\z% < fin . a
Theorem (3.1) suggests that the choice of ft should satisfy two criteria.
First, P must be small enough to bound 7^ 7*, and second, it should be large
enough to produce an adequate reduction in the number of nonzero beams.
12
These concepts can be applied to create a practical algorithm for choosing f3.
We first choose ft to be large enough to force a substantial reduction in the
number of beams and then we drive /? towards zero to produce a comfortable
bound on 7^ 7*. Section 4.4 presents the parameter selection algorithm.
The penalty problem (3.1) is difficult to solve because z* is discon
tinuous, but it can be formulated as the following mixed integer program.
min(7)a:,j,) 7 + jfle'y, p > 0
subject to (7, x) Â£ S,
EP x{p, b) r]y(b) <0, b Â£ beams
y{b) e {0,1}, b Â£ beams
where 77 is an upper bound for J2PX{P, b). This mixed integer program can be
solved with a branch and bound algorithm [1, p.486]. For a given /? and suf
ficiently large 77, (3.3) produces an exact /3 minimum support solution. To
see this, observe that the only way the constraint (3.3) will be satisfied is if
22px(p,b) > 0 ==> y(b) = 1. But if EPx{p, b) = 0, then y(b) will be chosen to
be 0 in order to minimize the objective function. Thus any solution (7, x, y) to
(3.3) is equivalent to (7,x, z*) where (7,x,z) is a solution to (3.1). Unfortu
nately the computational complexity of the branch and bound algorithm limits
the amount of data or detail that may be used to describe a specific instance
of the RTPP. An alternative to the mixed integer programming formulation is
to attempt to find an approximate rather than an exact solution to (3.1). Such
an approach is discussed in the following section.
13
3.2 The Exponential Approximation
In [6], Mangasarian presents a Successive Linearization Algorithm
(SLA) for solving concave minimization problems. We adapt this algorithm
to the radiotherapy treatment planning problem. This requires a modification
to the theory described in [6]. This work is concerned with finding minimum
support solutions to problems of the form
min
(x,y)eT
f{x,y)
(3.4)
where / : Rn R is a concave function on Rn that is bounded below on the
polyhedral set T. The RTPP satisfies these conditions: In particular, 7 is a
linear function and therefore concave, T is bounded below because one can not
deliver negative radiation, and T is bounded above because there is an upper
bound on the total radiation.
 ALPHA a .5
ALPHA 1
 ALPHA5
 ALPHA 20
 STAR NORM
F{X) = 1 EXP(ALPHA X). ALPHA > 0
Figure 3.1. Comparison of Exponential Approximation 1 e az to the Star
Norm \z\*.
Lemma 3.3 presents a method for approximating \z\r with a smooth
14
exponential function
e Â£~a* (3.5)
where a is a positive scalar, e denotes the base of the natural logarithm, and
ez = [eZl, , Â£Zn]. Figure 3.1 demonstrates how the the exponential function
approximates \z\* for z e R as a > oo. Notice how increases in a produce a
better approximation of \z\M. Some theoretical justification for approximating
(3.1) by the following continuous concave minimization problem:
min(w)eT 7 +/3e'(e Â£ az), a,/d>0, (3.6)
is given by the next lemma, which is based on [6, Lemma 2.2].
Definition 3.2 (Vertex) [1, Def 2.7] Let T be a convex set. A vector w G T
is a vertex of T if there exists some vector c such that c'w < c'y for all y
satisfying y e T and y ^ w.
Lemma 3.3 (Vertex Solution of Penalty Problem) Let / : Rm xRP > R
be a concave function that is bounded below on the polyhedral set
T C R x RP+. Then for each ft > 0, there exists a0(/3) > 0 such that for all
a > a0(j3),
min f(x,y) + /3e'(eÂ£~ay), (3.7)
(x,y)7
has a vertex solution, and every vertex solution of (3.7) is also a vertex solution
of
mm f(x,y) + f3e'\y\*.
(x,y)Â£T
Proof: For a > 0,
(3.8)
(es ^ < y since y > 0
(3.9)
15
so it follows that
min f(x,y)+/3e' (ee ay) < inf f(x, y) + /3e'\y\*. (3.10)
(x,y)eT v ' {x,y)Â£T
For any a > 0, (3.7) has a vertex solution. This follows from the fact that a
concave function bounded below on a polyhedral set not containing lines going
to infinity in both directions attains a minimum at a vertex of the polyhedral
set [8, Corollaries 32.3.3 k 32.3.4]. Let V be the set or vertices (x,y) of T
that solve (3.7) for arbitrarily large values of a. Since T has a finite number of
vertices, V is nonempty. It follows that there exists some finite number a0(/3)
such that for a > ao(P) the only vertex solutions of (3.7) are in V.
Let a > ocq (/?), and let (x, y) be a vertex solution of (3.7). Since (x, y)
is in V, it repeatedly solves
,min /(x, y) + /3e'(e Â£~aiV), (3.11)
(x,y)eT
for some sequence {a0, aj, } t oo. Hence for each c*j in the sequence,
f(x, y) + /0e' (e e~aiy) = min f(x, y) + fie? (e e~Qiy)
< inf f{x,y) +/3e'\y\* by (3.10).
Cx,y)ET
Letting i > oo gives
f(x,y) +Pe'\y\* < inf +Pe'\y\: (3.12)
(x,y)eT
Since (x, y) is a.vertex of T, it follows that (x, y) is a vertex solution of (3.8).
16
The following theorem, which is based on [6, Theorem 2.3], relates
the solution of the original problem (3.4) to the solution of the exponential
approximation (3.7).
Theorem 3.4 Let / : Rm x Rp R be a concave function that is bounded
below on the polyhedral set T C i? x Rp+. Then there exists Po > 0 such that
for each 0 < P < Po there exists a(P) such that for all a > a(P), (3.8) has a
vertex solution, and each such solution is also a vertex of (3.4).
Proof: For each P > 0, (3.8) has a vertex solution. This follows from the
fact that a concave function bounded below on a polyhedral set not containing
lines going to infinity in both directions attains a minimum at a vertex of the
polyhedral set [8, Corollary 32.3.3 h 32.3.4].
Let W be the set of vertices of T that solve (3.8) for arbitrarily small
P > 0. Since T has a finite number of vertices, W is nonempty. It follows that
there exists Po > 0 such that for 0 < P < Po, the only vertex solutions of (3.8)
are in W. Let 0 < P < Po, then (3.8) has a solution (x, y) that repeatedly solves
min /(s,y) + fte'y., (3.13)
(x,y)eT
for some sequence {Pi, P2, } i 0. Hence for each Pi in the sequence,
{x,y) 6 {arg vertexmin(liy)6r /(^, y) + Pie'\y\*} ^
= {arg vertex min(lij,)eT f(x, y) f(x, y) + fte'y,},
where (x, y) is defined as a solution of (3.4), that is
(x,y) e T and f(x,y) = min f{x,y). (3.15)
(i,y)eT
First we show by contradiction that
f(x, y) f{x, y) < 0 (or equivalently f(x, y) f(x, y) = 0). (3.16)
17
If not, then for sufficiently small $,
Pi
> max
e'y* ~Av\* M
/(,y)/(>y)A>J
Thus, we get the contradiction
e'y* = e'\y\* + j(f{x,y)f(x,y))
> J\vU + 7r(f{x,y) ~f(x,y))
Pi
> e'\y\, + e'\y\, e'\y\
(3.17)
where the first inequality follows from (3.14) and the last inequality from (3.17).
Hence (3.16) holds. We now show that e'\y\* is a minimum over the set of
solutions of (3.4), that is (x, y) has a minimal number of zeros besides being
optimal for problem (3.4). We have for /% < Po and all optimal (x,y), that is
(x,y) T and f(x,y) < f(x,y) :
e'li/l* = e'y* + 4(/(, y) f{x, y))
Pi
< e'y, + ^(/(x,y)/(,y))
Pi
< e'y*,
where the first inequality follows from (3.14) and the second inequality from
f(x,y) f{x,y) < 0. Hence (x,y) is an optimal vertex solution. Finally, by
Lemma 3.3, for each ft < f30, there exists a(/3) such that for all a > a(/3), (3.7)
has a vertex solution, and each such solution is also a vertex solution of (3.8),
which since ft < /?o is a vertex solution of (3.4).
18
This theorem states that for sufficiently small ft and sufficiently large
a, a P minimum support solution to (3.6) will also be a solution to (2.5). This
thesis deals primarily with Lemma 3.3 rather than Theorem 3.4. We search
for solutions to (3.1) rather than a solution to (2.5) with fewer beams because
we can not guarantee that the beam reduction for the latter will be sufficient.
Therefore we seek to exert some control over the number of beams in a solution
to (3.1) by controlling f3, where ft does not necessarily satisfy /? < Pq.
3.2.1 The Successive Linearization Algorithm
The exponential approximation (3.5) is a much nicer function to work
with than z*, but (3.7) is still difficult to solve because the objective function
is not convex. Mangasarians successive linearization algorithm takes advan
tage of the concavity and differentiability of (3.7) by solving a series of linear
programs with the first order Taylor approximation of (3.7) as the objective
function.
Algorithm 3.5 Successive Linearization Algorithm (SLA)
[6, Algorithm 3.1] Choose a > 0 and P > 0 and starting point
(70, x, z) G i?+ x J?" x R+ and a stopping tolerance e. For i = 0,1,... choose
(7i+1, xi+1, zi+1) G {arg vertex min(w)Â£T (7 7*) + aP(Â£az')'{z z1)}.
Stop if
(7l+1, xl+1, zl+l) G T and (7*+1 7*) + aP{e~az')'{z%+l zl) < e.
The algorithm terminates when (7*+1, x'+1, ^+1) satisfies first order optimality
conditions; that is at a stationary point of the original linear program (2.5). To
19
ensure that (7*+1, rE*+1, zl+1) is a vertex solution, we use a simplex algorithm,
rather than an interior point algorithm.
20
4. Implementation
Before we apply the SLA to the Radiotherapy Treatment Planning
Problem we need to collect and format the data that will be used to define
a particular problem instance. In the previous sections we defined the dose
deposition operator D and the three pixel subsets: critical, body, and tumor.
In this section we begin by describing the processes used to calculate D and to
create the critical, tumor and body sets. We conclude with our algorithm for
choosing a and /3.
4.1 Creating the Pixel Sets
The inequality constraints of problem (2.5) require that we be able
to classify pixels from the CT image into one of three different pixel subsets:
critical, body, and tumor. We created a pixel selection routine with the library
functions from MATLABs Image Processing Toolbox. The user begins with a
datafile that contains the patient image, usually a CT image with 512 x 512
pixels. The MATLAB routine reads the image file and stores the image. Our
subroutine repeatedly prompts the user to select the appropriate areas of the
image that contain either critical, tumor, or body pixels with a mouse. For
example, if prompted to select a tumor, one uses the mouse to outline the
tumors that appear in the CT image on the computer monitor. Once the user
has finished selecting pixels, MATLAB creates a file in AMPL format that lists
each pixel in either the body, critical, or tumor sets. This file is part of the
21
data that is needed to solve a particular instance of the RTPP with the AMPL
model.
4.2 Calculating the Dose Deposition Operator
3B Body
Beam Boundary
Figure 4.1. Divergence Angle
In Section 2.1 we briefly discussed the dose deposition operator D.
This section gives a detailed description of the method used for calculating
D. The center of the patient image, or isocenter, is 100 cm from the source
point in the gantry. The gantrys position or angle describes the placement of
the source point in relation to the fixed isocenter of the patient image. This
ensures that the center of each beam passes through the center of the pa
tient. In this model of the dose deposition operator D, we take into account
the natural divergence of the radiation waves as they move away from the
22
source point in the gantry. The outside boundaries of the beam make an acute
angle of .39479112 radians that is bisected by the line connecting the radiation
source in the gantry to the center of the patient image, see Figure 4.1.
In Order to model the dose deposited by a beam to a pixel, we create
a discrete model of the beam by subdividing each beam into a fixed number
of rays, called pencils. Each pencil has a width associated with it because the
rays spread as the radiation moves from the source point through the patient.
However, in our model, we ignore this width and instead assume that all radi
ation is concentrated along the center line of the pencil. As radiation travels
along this line, it deposits energy into each pixel that it intersects. The amount
deposited is a function of the intensity of the radiation as it enters the pixel
as well as the distance traveled within the pixel. The majority of the follow
ing discussion is concerned with determining the pixel and pencil intersections,
since this is the most difficult aspect of calculating D. Once these intersections
are determined a straight forward calculation supplies the fractional value of
the dose deposited by a pencil x(p, b) to pixel i.
4.2.1 Finding pencil/pixel intersections
The patient image and the source point can be placed on a Cartesian
plane with the isocenter of the image at the origin. Since we know the length
of the ray extending between the source point and the isocenter as well as the
angle 0 this ray makes with the xaxis, it is easy to convert the source point to
a rectangular position on the plane. As the central vector of a pencil passes
through the medical image, it will intersect the fixed vertical and horizontal
23
pixel boundaries. This vector can be parameterized according to equations
x(r) = $x r cos(9 + Oii) (4.1)
l/(0 = Sy rsin(# + a*) (4.2)
were {sx,sy) is the (x,y) coordinates of the source point, and Oii is the angle
the ith central vector makes with the center of the beam. In Figure 4.2 the
central vector of P2 passes through the tumor, but the central vector of PI
does not. In this model, radiation is deposited to a pixel if and only if the
central vector intersects the pixel.
Figure 4.2. Dose Deposition Model.
In order to determine how much radiation is deposited in each pixel,
we construct a path through the image in stages. This path consists of a
24
sequence of points {(xk, yk)} lying on the central vector, which represent where
the central vector enters a new pixel in the image. For the first point in the
path (xi,yi) we determine the value of r at which the pencil first enters the
medical image. This is accomplished by choosing the minimum r that satisfies
the conditions:
xio < x(r) < xhi (4.3)
Via < y{r) < Vm (4.4)
Where Xi0, x^, yi0, and yhi represent the boundaries of the image. Once a
point (Xk,yk) on the path has been determined, the next step is to identify
which pixel is being entered at that point. The vertical and horizontal pixel
boundaries are stored in a one dimensional array q of length m + 1. Where m
is the number of pixels along each side of the image. We search q for elements
l and j that satisfy qi < y < qi+i and qj < x < qj+i For the first point
in the path, only one pair of indices (l,j) will satisfy this condition, thereby
identifying the pixel in row l, column j, of the image. Thereafter, the process
will identify two pairs of indices: the pixel being exited and the pixel being
entered. Since we already know the pixel being exited we choose the other as
the pixel being entered. Once we have determined which pixel is being entered,
we then calculate the next point on the path (xk+i, yk+i), which is the point
where the central vector exits the pixel by finding the maximum r satisfying
qi
We continue this process until the pencil exits the image.
25
4.2.2 Calculating the dose
Once the path {(xk,Vk)} has been calculated, we can then determine
the dose deposited to each pixel along this path. The radiation beam attenuates
at different rates, which depend on the distance traveled, the material passed
through and the starting energy level. We take these factors into consideration
when calculating the dose deposited to a particular pixel. The intensity at the
entrance to pixel k is given by the iterative formula:
/* = (4.5)
where 70 is the energy level at the source, //0 is the attenuation factor for air,
lik1 is the attenuation factor associated with the tissue or material in pixel
(fc 1), for k > 2, d0 is the distance from the source point to (xi, yi), and
dki = (Â£fc,Uk) ~ (zfci,i/fci) (4.6)
is the distance the center vector of the pencil travels through the (k l)st pixel
in the path. The dose, A*, for the 7th pixel is the intensity level Ik as we enter
the 7th pixel minus the intensity as we exit the pixel Ik+i
A k = Ik h+1 = (1 ekdk)Ik. (4.7)
Thus, the fraction of I0 that is deposited in pixel k is That is we assign
D(i,p,b) = where % is the index into the dose deposition array of the 7th
pixel in the path.
4.3 The // File
Up until this point we have classified pixels into the body, critical
or tumor sets. In this section we are concerned with describing the pixels of
26
the image in terms of bone, muscle, air and water. To create the necessary
information for a problem instance, we need to categorize the different types
of pixels in a CT image into one of these four categories, because radiation
attenuates at different rates as it passes through these pixels.
Using a method similar to that described in Section 4.1, the user
selects different structures in the body so that the algorithm for the dose depo
sition operator can calculate the fractional dose for each pixel. The MATLAB
routine allows the user to look at a CT image and to classify pixels. For exam
ple, if the user is prompted to select bone, he uses a mouse to outline all bones
that appear on the computer monitor. The routine continues in this fashion
until all pixels have been categorized. MATLAB writes this information to a
text file that can be read by the dose deposition operator routine.
4.4 Choosing a and /?
So far, this chapter has dealt with the collection and storage of the
data needed to solve a particular instance of the RTPP. This section describes
an augmented successive linearization algorithm (ASLA); which automates the
process of choosing a and /?. Sections 4.4.1 and 4.4.2 describe the a and (3
selection heuristics utilized by the ASLA and Section 4.4.3 gives the pseudo
code for the ASLA.
4.4.1 Choosing a
Practical experience has shown that the smoothing parameter a can
affect the quality of the solution found by the successive linearization algorithm.
Producing changes in solutions to (2.5) with many beams can be difficult;
27
therefore it is important to understand how a affects the penalty coefficients
associated with each beam. The ability of the SLA to penalize all beams with
positive levels of radiation during iteration i + 1 depends on zl and a. Therefore
let us consider the coefficients associated with each beam.
Figure 4.3. g{a,/3) := afte aZi, for fixed f3 and a > 0.
The first order Taylor series expansion of (3.6) about the point z1 is
/(7, x, z) := 7 + aj3Â£~azhi Hh a8e~az'nzm (4.8)
which is the linear objective function for the (i + l)st iteration of the SLA.
We are interested in examining how the coefficients of the linear function (4.8)
behave with respect to a for fixed values of 2] for j 1,..., n. To do this,
define
g{<*, P) :=
28
which corresponds to the jth coefficient in (4.8). As long as /? is bounded
we have lim^oo g(a, f3) = 0. Figure 4.3 illustrates how g(a,(3) depends on
different choices of a and zj, for a fixed value of /3. The rate at which g(a, /3)
converges to zero as a oo is substantially impacted by the relative size of
zj. This feature plays an important part in the way that we choose a. In
particular, we will choose an initial a based on the size of the largest beam in
a solution to linear program (2.5). The following discussion explains the logic
behind our choice.
for a fixed value of Zj. Then from Figure 4.3 we can see that d* > 0 as jo
iner eases. As it turns out the maximum dose distributed by a single beam in
an RTP can be very large. This poses a difficulty when choosing a, because
if a is too large, the penalty associated with larger beams will be zero since
limQxx, g(a, (3) = 0. Let z* represents the beam intensities for a solution to
the unpenalized RTPP (2.5), and let z^^ be the largest component of z*. The
largest penalty one can hope to place for a fixed z* occurs when a := d^. To
find dj^ we can solve (4.10) analytically using the first order condition:
yielding a^ax = l/T^ax The analysis shows that a^ax is independent of /?. In
order to ensure that all the beams are penalized in the first iteration of the
Let Oj be a solution to
maxQ g(a,/3)
subject to a > 0.
(4.10)
(4.11)
29
ASLA, we initialized a according to the formula
= JT <4'12)
^max
However, this doesnt always succeed in sufficiently reducing the number of
beams, in this case we need to increment a.
4.4.2 Choosing /?
Theorem 3.1 gives an error bound for the difference between 7* and
7^. We can use this result to make an educated guess about /3 that is dependent
upon the type of RTP the physician wants to produce. In particular, we choose
/? by the formula:
r 7*
*:= i?ir (413)
where T is the highest dose to critical structures allowed by the radiation
treatment staff. By Theorem 3.1, we see that this is the largest possible choice
of /3 that guarantees the maximum critical structure dose in a solution to (3.1),
will not exceed T. However, this choice of /3 is often too conservative, producing
solutions with an inadequate reduction in the number of beams and very little
change in the maximum critical structure dose. In such cases, a large value of
/? is called for. It therefore makes sense to increment (5 and try again.
The ASLA, which we present in the next section, has two phases. The
first phase ora phase is executed with a fixed large value for /3, and does not
directly affect the selection of /?. The second phase or ft phase begins with ft
chosen according to (4.13). If the beam reduction is adequate and maintains a
sufficiently small 7^ 7*, then we terminate the algorithm and is the best
approximate f3 minimum support solution. Otherwise, we increase j3 and
30
solve for x via the SLA, continuing in this manner until the staffs criteria are
met or the maximum iteration is reached.
4.4.3 The Augmented Successive Linearization
Algorithm
The ASLA is an iterative method that determines acceptable parame
ter values a and ft and utilizes them to solve the RTPP with a slightly modified
version of Mangasarians SLA. The solutions resulting from this process are not
exact /3 minimum support solutions, but they do satisfy the RTP conditions
set forth in the introduction.
Algorithm 4.1 Augmented Successive Linearization Algorithm
(1) Solve the unpenalized LP (2.5) giving the solution (7*,a:*). Set 70 = 7*,
x0 = x*, z* = EPx*(p, )> zo = and choose ainc, f3inc > 0.
(2) (a) Choose the smoothing parameter a according formula (4.12).
(b) Set the penalty parameter /3 to some large positive value.
(c) Set 7 = 00.
(3) For 1 = 1,..., max.
(a) Use SLA to generate 7*, xl and zl.
(b) If 7i < 7 and \zl\* < Z
(1) Set (7,5) = {ll,xl) and a = o.
(2) If 7 7* < clinical error tolerance, then goto step 4.
(c) Set ot ot H cxjjiQ.
(4) Set penalty parameter /? according formula (4.13) and set smoothing pa
rameter a = a.
(5) For 1 = 1,..., max.
31
(a) Use SLA to generate j1, xl and zl.
(b) If 7* < 7 and \zl\* < Z
(1) (7,x) = ('yl,xl)
(2) 0 = 0
(3) Break with the approximate minimum support solution
(7, x) and parameters & and /3.
(c) Set /? = /? + Pinc.
32
5. Results
In order to determine if the ASLA produces good solutions to the
RTPP we ask the following questions:
(1) Does the ASLA produce an adequate reduction in the number of beams
used to satisfy the radiation dose requirements?
(2) In most cases we expect that the maximum dose of radiation dis
tributed to a critical pixel will be greater than that produced by the
linear program solution. Is the increased dose produced by the ASLA
acceptable? How does the ASLA solution compare with the mixed
integer programming solution in this respect?
(3) Finally, how does the ASLAs computational time compare with the
linear programming and mixed integer programming computational
times?
In order to answer these questions we examine several RTPPs.
In the following sections, the RTPPs are defined by several patient
images, which are shown in Figures 5.2, 5.3, 5.4, 5.7, 5.10, and 5.13. In these
images, the tumors are outlined in red and the critical structures are outlined
in blue. Thus far we have presented three solution methods for solving the
RTPPs. The first method solves the linear programming formulation (2.5).
With this method, an optimal plan can be found that satisfies all but the last
RTP condition set forth in Section 1; that is the solutions may have too many
33
beams to be implemented in a treatment session. The second method uses
the mixed integer programming formulation (3.3). Optimal solutions to this
problem are dependent upon the penalty term /3, and can be difficult to solve
in a reasonable time frame. The final method is the ASLA, see Section 4.4.3.
In order to determine the effectiveness of the ASLA, we needed to examine the
types of solutions produced by each method. Hence we modeled (2.5), (3.3),
and the ASLA with AMPL, A Mathematical Programming Language [3], and
solved each model with the CPLEX optimization software package [4].
CPLEX is a very efficient solver, but large RTPPs can cause difficul
ties for all three methods; it is important to understand how we controlled the
problem size. We have described the RTPP as a structured beam selection
process. Therefore the number of beams available during the selection process
determines both the quality and complexity of the solution; a large number
of beams translates to a more complicated dose deposition operator D. Ad
ditionally we need to consider the accuracy with which the patient image is
constrained. The patient images we used have 512 x 512 pixels, which provides
a substantial amount of detail. Unfortunately, this also translates to a massive
dose deposition operator D. To simplify our storage requirements we sacrifice
some detail by reducing the patient image. For example we can represent the
image with 64 x 64 pixels in the AMPL model.
5.1 Solutions for Small Problems
We solved RTPP1 with the linear programming model, the mixed
integer programming model, and the ASLA several times varying the maximum
34
number of beams available during the selection process. Figure 5.1 is a graph
of the results. One can see from Figure 5.1 that the number of nonzero beams
for an linear programming solution tends to be much larger then 10 beams.
The ASLA reduces the number of nonzero beams, but does have more beams
then the mixed integer programming solutions.
COMPARISON OF MODELS
Figure 5.1. Comparison of Solution Methods.
Table 5.1 contains the results obtained by solving the RTPP1 and
RTPP2 with each method. The results in Table 5.1 are for a 16 x 16 pixel
representation of Figure 5.2 and a 32 x 32 pixel representation of Figure 5.3,
both with a variable number of beams.
The results in Table 5.1 help us in determining the efficacy of the
ASLA. The approximate solutions produced by the ASLA appear to be a nice
compromise between an optimal value for 7, the maximum critical dose, and
35
the number of beams in a particular solution. That is the number of beams in
an ASLA solution is less than 10, and the increase in maximum dose delivered
to a critical pixel is less then 3 units of radiation. It is important to note that
both these criteria can be controlled to some degree by the ASLA. Furthermore
if 7* is the optimal value for the linear programming formulation of the RTPP
(2.5) and 7 is the optimal value produced by solving (3.6) with the ASLA, then
the results presented here satisfy 7 7* < 3 < fin, the bounds established by
Theorem 3.1 for the mixed integer program.
In Section 4.4.3 we discussed how the ASLA determines good choices
for a and fi on the fly. This is efficient and effective only because the ASLA
developed by Mangasarian is very fast. The same simple iterative heuristic is
not effective for the mixed integer programming algorithm. The ASLA com
putation times listed in Table 5.1 include the time needed to find levels for
a and fi that produce solutions with an acceptable number of beams and an
acceptable 7. This is not the case with the mixed integer programming times.
The mixed integer programming times are incredibly large, in fact a problem
with a 32 x 32 pixel image with 50 beams took almost 19 hours to solve for
a fixed value of fi = 5. Unfortunately this wasnt a particularly good choice
for fi so 7 was much too large, and the solution was useless. This compares
unfavorably to the ASLA that returned a very good solution in about 2.5 min
utes. It is interesting to note that the 7 values for both the mixed integer
program and ASLA are usually quite close if we are fortunate enough to pick
fi correctly. In these cases it should be noted that the mixed integer program
uses far fewer beams to achieve similar results. Figure 5.1 demonstrates that
36
the ASLA reduces the number of nonzero beams in a consistent manner. As
the number of possible beams in the model increases, very little change occurs
in the number of nonzero beams produced by the ASLA that solve the RTPP
presented in Figure 5.2.
37
algorithm 7 P a time beams used figure pixels
LP 3.42847   0:00:00.76 8/10 5.2 16 x 16
ASLA 4.87314 5 .14 0:00:00.87 3/10 5.2 16 x 16
MIP 7.12787 5  0:00:01.91 2/10 5.2 16 x 16
LP 2.4244   0:00:02.94 15/20 5.2 16 x 16
ASLA 4.39832 5 .02 0:00:03.00 5/20 5.2 16 x 16
MIP 4.28042 5  0:00:06.19 2/20 5.2 16 x 16
LP 2.40245   0:00:03.16 19/30 5.2 16 x 16
ASLA 3.992 5 .01 0:00:05.93 6/30 5.2 16 x 16
MIP 4.50354 5  0:00:16.11 2/30 5.2 16 x 16
LP 2.33075   0:00:05.67 20/40 5.2 16 x 16
ASLA 4.04568 5 .01 0:00:06.95 6/40 5.2 16 x 16
MIP 4.28042 5  0:00:28.63 2/40 5.2 16 x 16
LP 2.31355   0:00:06.30 22/50 5.2 16 x 16
ASLA 4.15199 5 .01 0:00:08.67 6/50 5.2 16 x 16
MIP 4.32769 5  0:01:08.78 2/50 5.2 16 x 16
LP 2.29754   0:00:08.06 22/60 5.2 16 x 16
ASLA 4.08298 5 .01 0:00:09.80 6/60 5.2 16 x 16
MIP 4.20276 5  0:02:00.50 2/60 5.2 16 x 16
LP 2.2893   0:00:07.41 20/70 5.2 16 x 16
ASLA 3.87975 5 .01 0:00:06.83 6/70 5.2 16 x 16
MIP 3.85661 5  0:03:54.51 2/70 5.2 16 x 16
LP 2.28781   0:00:09.56 25/100 5.2 16 x 16
ASLA 4.05298 5 .01 0:00:10.03 6/100 5.2 16 x 16
MIP 4.18351 5  0:16:24.04 2/100 5.2 16 x 16
LP 2.28033   0:00:15.75 21/180 5.2 16 x 16
ASLA 3.96762 5 .01 0:00:17.99 5/180 5.2 16 x 16
MIP 3.98676 5  2:03:59.83 2/180 5.2 16 x 16
LP 3.73577   0:00:11.18 26/30 5.3 32 x 32
ASLA 4.70809 4 .49 0:02:31.57 9/30 5.3 32 x 32
MIP 8.89508 5  2:41:46.93 3/30' 5.3 32 x 32
LP 3.54444   0:00:19.69 38/50 5.3 32 x 32
ASLA 4.49706 1 .37 0:03:39.18 9/50 5.3 32 x 32
MIP 8.03349 5  18:49:00.00 3/50 5.3 32 x 32
Table 5.1. Results for RTPP1 and RTPP2
38
5.2 Solutions for Large Problems
In the last section we controlled the sizes of RTPP1 and RTPP2 by
reducing the number of beams in the model and by storing the image with
a smaller number of pixels. The main reason for this was so that we could
solve RTPP1 and RTPP2 with the mixed integer program and compare the
various results. We mentioned earlier that the solutions for RTPPl and RTPP2
detailed in Table 5.1 dont have enough clinical detail. In this section all
problems have been modeled with the largest representation possible, that is
on a 64 x 64 pixel image with 180 beams. The results in this section are also
promising but the solutions times for the ASLA have increased at a greater
rate then the times for the linear program. Table 5.2 shows that the ASLA
produced satisfactory 7 values and reduced the number of beams substantially.
Graphs of the linear programming and ASLA dose results for each problem are
presented at the end of the chapter.
algorithm 7 P a time beams used figure
LP 2.48557   7:50.65 100/180 5.4
ASLA 4.82768 1.0333 0.06009 36:20.86 8/180 5.4
LP 2.8565   7:14.02 114/180 5.7
ASLA 6.79427 2.02885 0.06102 25:52.37 8/180 5.7
LP 2.36607   6:31.38 69/180 5.10
ASLA 5.35838 2.05085 0.02072 17:54.62 9/180 5.10
LP 1.62386   6:44.57 92/180 5.13
ASLA 4.24193 1.03659 0.05109 17:22.48 7/180 5.13
Table 5.2. Results for Large Problems: 64 x 64 pixels
39
5.3 Conclusions
This thesis was motivated by the need to find solutions to the Ra
diotherapy Treatment Planning Problem that satisfy the RTP conditions, i.e
be simple enough to administer during a treatment session. We have applied
Mangasarians successive linearization algorithm to the Radiotherapy Treat
ment Planning Problem, and added a simple heuristic that determines good
values for the parameters a and /?. We then demonstrated that this solution
method gives good approximate results when compared with the exact solu
tions for small problem instances. This approach seems to be slower, when
the problem size increases but is still much faster and more accurate than the
exact method.
40
Figure 5.2. RTPP1
41
Figure 5.3. RTPP2
42
Figure 5.4. RTPP3
43
1
Figure 5.5. RTPP3 LP solution
44
T "11 *'" *V
Vl ____ I ^ __________________J _. jj
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Figure 5.6. RTPP3 ASLA solution
45
Figure 5.7. RTPP4
46
Figure 5.8. RTPP4 LP solution
47
0.9 1
Figure 5.9. RTPP4 ASLA solution
48
Figure 5.10. RTPP5
49
Figure 5.11. RTPP5 LP solution
50
Figure 5.12. RTPP5 ASLA solution
51
Figure 5.13. RTPP6
52
i
 ^
0.9 1
Figure 5.14. RTPP6 LP solution
53
7T"
0.4
0.5
T
1
0.6
I
0.7
X
0.8
T*
0.9
1
Figure 5.15. RTPP6 ASLA solution
54
Appendix A. MATLAB Subroutines
A.l Create Set and fi Files from Patient Image
function [x_b,y_b,x_c,y_c,x_t,y_t] = image_script(I, ...
set_file, mu_file, dim)
/.Image_script .m creates data files needed to run dose2.cpp and
/.also creates set labels for ampl
/.I =imread(image_name) this is a square matlab image
'/.matrix (input file) or I can be loaded from a *.mat file
/.set_file is a text filename with the ampl set data (output file)
'/0mu_file is a text filename for the mu matrix (output file)
/.dim is the dimension of the pixel reduced image.
y = *y;
n = n;
'/,Create sets
ans = input(Do you want to create ampl sets? (y or n));
if (ans == y )
fprintf(Create Body, Critical Tissue and Tumor SetsAn);
[x_b,y_b,x_c,y_c,x_t,y_t] = make_sets(dim,I,set_file);
/.Save set figure now
fprintf(Use the figure window file options to save as *.fig.\n);
fprintf(Press return to continue.\n);
pause
end;
/.Create table with tissue types
ans = input(Do you want to create a mu matrix? (y or n));
if (ans == y )
fprintf(Select Water, Air, Bone and Muscle Tissues for
mu value.\n);
select_tissue(I,dim,mu_file);
end;
55
fprintf(IMPORTANT: Save workspace before continueing this
session.\n);
56
A.2 Writes /i Text File to Disk
function mu_table = make_mu_table(T,mu_air,mu_bone,...
mu_muscle,mu_water,dim)
/.This function takes the tissue matrix T and assigns the
/.appropriate mu values. The data is written to file so
/.that it can be used to calculate the dose deposition
/.operator. For more info on T see select_tissue.m
/.Initialize data
air = 1;
bone = 2;
muscle = 3;
water = 4;
for i=l:dim
for j=l:dim
*/,Assign mu values for pixels that represent air
if T(i, j) == air
mu_table(i, j) = mu_air;
/.Assign mu values for pixels that represent bone
elseif T(i,j) == bone
mu_table(i,j) = mu_bone;
/.Assign mu values for pixels that represent muscle
elseif T(i,j) == muscle
mu_table(i,j) = mu.muscle;
/Assign mu values for pixels that represent water
else
mu_table(i,j) = mu_water;
end;
end;
end;
/.Write data to file
fid = fopen(filename,a);
max_row = 5;
57
A.3 Creates AMPL Sets
function [x_b,y_b,x_c,y_c,x_t,y_t] = make_sets(dim,I.filename)
%x_b, y_b are the vertices of the polygon that outlines the body
'/,x_c, y_c are the vertices of the polygons that outline
'/critical struc.
/,x_t, x_t are the vertices of the ploygons that outline
/the tumors.
/Initialize data
fid = fopen(filename,a);
elf;
[n,m] = size(double(I)+l);
x_b = 1;
y_b = 1;
x_c = 1;
y_c = 1;
x_t = 1;
y.t = 1;
AIR = zeros(dim);
CRITICAL = zeros(dim);
TUMOR = zeros(dim);
BODY = zeros(dim);
y = y;
n = n;
max_row = 5;
/Select body region
fprintf(Use cross hair to mark off the body boundary.\n)
[I2,x_b,y_b] = roipoly(I);
AIR = double(imresize("12,[dim,dim]));
hold on;
/Select critical pixels
fprintf(JUse cross hair to mark off critical areas.\n)
ans = y;
while(ans == y)
[I3,xi,yi] = roipoly;
58
x_c = [x_c;xi;l];
y_c = [y_c;yi;l];
line(xi,yi,color',b,linestyle,,linewidth,2);
CRITICAL = double(imresize(13,[dim,dim] )) + CRITICAL;
ans = input(Do you want to select more critical
tissue? (y or n));
end
/Create critical set
fprintf(fid,set critical := );
make_labels(CRITICAL,fid,max_row);
'/Select tumor pixels
ans = y;
fprintf(Use cross hair to mark off tumor(s)An)
while(ans == y)
[I4,xi,yi] = roipoly;
x_t = [x_t;xi;l];
y_t = [y_t;yi;l];
line(xi,yi,color,r,linestyle,,linewidth,2);
TUMOR = double(imresize(14,[dim,dim])) + TUMOR;
ans = input(Do you want to select a tumor? (y or n));
end
legend(critical tissue,tumor);
'/Create tumor set
fprintf(fid,set tumor := );
make_labels(TUMOR,fid,max_row);
/Create body set
fprintf(fid,set bod:= );
make_labels("(AIR + CRITICAL + TUMOR),fid,max_row);
fclose(fid);
59
A.3.1 Writes AMPL Set File to Disk
function make_labels(A,fid,max_row)
[n,m] = size(A);
count 1 = 0;
count2 = 0;
max = length(find(A));
for i=l:n
for j=l:m
if A(i,j)
count2 = count2 + 1;
if(count2 < max)
if (countl < max_row )
fprintf (fid,3 /.g_/.g, , i, j);
countl = countl + 1;
else
f printf (f id, /,g_/.g, \n , i, j ) ;
countl = 0;
end
else
fprintf (f id, %g_/,g; \nJ, i, j);
end
end
end
end
60
Appendix B. AMPL Code
B.l Linear Programming Model
#Model Name: rtp_lp.mod
#Programmmer: Janine Kennedy
#Date of last alteration: May 9, 2000
#
#Indexes for the patient image and beams, this data to be
#initialized by files created with MATLAB image processing
#software.
set pixels ordered; #element format row_column of image
set pencils ordered; #element format pi .. plO ...
set beams ordered; #element format bl .. blO ...
set critical within pixels; #c element of critical
set tumor within pixels; #t element of tumor
set bod within pixels; #o element of tumor
#The d = dose_deposition gives the fractional amount of radiation
#capable of being deposited by a particular pencil x
#into a particular pixel. This data is initialized
#with a file created by C program for calculating the
#dose deposition operator.
param d{pixels,pencils,beams} >= 0, default 0;
#Parameter intialization
param tumor_bounds_low{tumor} >= 0, default 9;
param tumor_bounds_up{tumor} >= 0, default 10;
param body.bounds{bod} >= 0, default 8;
param dim, default sqrt(card(pixels));
param curr.card >= 0, default 0;
#
#Decision variables
var gamma; #max dose delivered to critical pix
#radiation strength of each pencil
61
var x{p in pencils, b in beams}>= 0;
#
#0bjective function: a penalty of Beta is added for
#each nonzero pencil x
minimize optimal_val: gamma
#
#Defines the radiation limits for each tissue type
subject to critical_tiss_bounds {c in critical}:
gamma (sum{p in pencils,b in beams}
d[c,p,b]*x[p,b]) >= 0;
subject to tumor_tiss_bounds_low {t in tumor}:
sum{p in pencils,b in beams}
d [t,p,b]*x[p,b] >= tumor_bounds_low[t];
subject to tumor_tiss_bounds_up {t in tumor}:
sumCp in pencils,b in beams}
d [t,p,b]*x[p,b] <= tumor_bounds_up[t];
subject to body_tiss_bounds {o in bod}:
sum{p in pencils,b in beams}
d[o,p,b]*x[p,b] <= body.bounds[o];
#
62
B.2 Mixed Integer Programming Model
#Model Name: rtp_mip.mod
#Programmmer: Janine Kennedy
#Date of last alteration: May 9,2000
S
#Indexes for the patient image and beams, this data
#to be initialized by files created with MATLAB image
#processing software.
set beams ordered; #b element of beams
set pencils ordered; #p element of pencils
set pixels ordered; #i element of pixels
set critical within pixels; #c element of critical
set tumor within pixels; #t element of tumor
set bod within pixels; #o element of body
#The d = dose_deposition gives the fractional amount of
#radiation capable of being deposited by a particular pencil x
#into a particular pixel. This d is initialized with a file
#created by a C program.
param d{pixels,pencils,beams} >= 0, default 0;
#Parameter intialization
param BETA default 5; #penalty for nonzero pencils
param EPS default le07; Slower bound for x
param M default 100000; Supper bound for x
param curr_card default 0; Snumber of beams in solution
param doseCpixels} default 0; Sdose delivered by solution x
param tumor_bounds_low{tumor} >= 0, default 9;
param tumor.bounds_up{tumor} >= 0, default 10;
param body_bounds{bod} >= 0, default 8;
S
SDecision variables
var gamma; Smax critical dose delivered to pix
Sradiation strength of each pencil
var x{p in pencils,b in beams}>= 0;
var y{b in beams} binary; Sbinary penalty variables
#
63
#Objective function: a penalty of beta is added for
#each nonzero pencil x
minimize optimal_val: gamma + BETA*sum{b in beams} y[b];
#
#Defines the radiation limits for each tissue type
subject to critical_tiss_bounds {c in critical}:
gamma (sum{p in pencils,b in beams}
d[c,p,b]*x[p,b]) >= 0;
subject to tumor_tiss_bounds_low {t in tumor}:
sum{p in pencils,b in beams}
d[t,p,b]*x[p,b] >= tumor_bounds_low[t] ;
subject to tumor_tiss_bounds_up {t in tumor}:
sumCp in pencils,b in beams}
d[t ,p,b] *x[p,b] <= tumor_bounds_up[t];
subject to body_tiss_bounds {o in bod}:
sum{p in pencils,b in beams}
d [o, p, b] *x [p, b] <= body ..bounds [o] ;
#
#The variable y = 0 when x = 0, and 1 otherwise, it
#replaces the step function in the objective function.
#(sum(x) = 0) => (y = 0), and (sum(x) > 0) => (y in {0,1})
subject to beam_off{b in beams}:
(sum{p in pencils}x[p,b]) (EPS y[b])>= 0;
#(sum(x) > 0) => (y = 1), and (sum(x) = 0) => (y in {0,1})
subject to beam_on{b in beams}:
(M y[b]) (sum{p in pencils}x[p,b])>= 0;
#
64
B.3 ASLA Script Files
B.3.1 ASLA Model
#Model Name: rtp_sla.mod
SProgrammmer: Janine Kennedy
#Date of last alteration: May 8, 2000
#
#Indexes for the patient image and beams, this data to be
#initialized by files created with MATLAB image processing
#software.
set pixels ordered; #element format row_column of image
set pencils ordered; #element format pi . plO ...
set beams ordered; #element format bl .. blO . .
set temp_beam ordered by beams; #used to calc. card, of sol. x
set critical within pixels; #c element of critical
set tumor within pixels; #t element of tumor
set bod within pixels; #o element of tumor
#
#The d = dose_deposition gives the fractional amount of radiation
#capable of being deposited by a particular pencil x
#into a particular pixel. This d to be initialized by file
#created with a C program.
param d{pixels,pencils,beams} >= 0, default 0;
#Parameter initialization
param tumor_bounds_low{tumor} >= 0, default 9;
param tumor_bounds_up{tumor} >= 0, default 10;
param body_bounds{bod} >= 0, default 8;
param BETA, default 0; #penalty for nonzero pencils
param BETA_min, default 0; #BETA lower bound
param ALPHA >= 0, default .01; #expon. smoothing parameter
param EPS >= 0, default le07; #tolerance level
param M >= 0, default 200; Supper bound for number for loops
param MAX_BEAMS >= 0, default 10; Supper bound on beams
param GAMMA_UP >= 0, default 5; Supper bound on gamma
param GAMMA_ERR0R, default 3;
param X, default 0; Sfixed param for exp. objective
param BETA_opt; Sbest value found by ASLA
65
param ALPHA.opt; #best value found by ALSA
param sla.gamma >= 0, default 0;
param sla_sol{p in pencils, b in beams} >= 0, default 0;
param alpha_sol{p in pencils, b in beams} >= 0, default 0;
param beta.soHp in pencils, b in beams} >= 0, default 0;
param curr_sol{p in pencils, b in beams} >= 0, default 0;
param gamma.alpha, default Infinity;
param gamma.beta, default Infinity;
param gamma_diff, default Infinity;
param optimal.val.old, default Infinity;
param optimal.val.alpha, default Infinity;
param optimal.val.beta, default Infinity;
param stop.criteria, default Infinity;
param count, default 0;
param count2 default 0;
param dim, default sqrt(card(pixels));
param dose{pixels} >= 0, default 0;
param ALPHA.INC >= 0, default .01;
param BETA.INC >= 0, default 1;
param curr.card >= 0, default 0;
param opt.card >= 0, default Infinity;
#
#Decision variables
var gamma; #max dose delivered to critical pix
#radiation strength of each pencil
var x{p in pencils, b in beams} >= 0;
var alpha >= 0;
#
#0bjective function: a penalty of Beta is added for each
#nonzero pencil x
minimize linear.approx: gamma + ALPHA*BETA/card(pencils)*
(sum{p in pencils, b in beams}(exp(ALPHA/card(pencils)*
(sum{i in pencils} curr_sol[i,b]))*x[p,b]));
#
#Defines the radiation limits for each tissue type
subject to critical.tiss.bounds {c in critical}:
66
gamma (sum{p in pencils,b in beams}
d[c,p,b]*x[p,b]) >= 0;
subject to tumor_tiss_bounds_low {t in tumor}:
sum{p in pencils,b in beams}
d [t,p,b]*x[p,b] >= tumor_bounds_low[t];
subject to tumor_tiss_bounds_up {t in tumor}:
sum{p in pencils,b in beams}
d[t,p,b]*x[p,b] <= tumor_bounds_up[t] ;
subject to body_tiss_bounds [o in bod}:
sum{p in pencils,b in beams}
d[o,p,b]*x[p,b] <= body_bounds[o];
#
67
B.3.2 SLA Script File
#Command file: rtp_sla.run
#Programmer: Janine Kennedy
#Date of last alteration: May 8, 2000
#
#SLA Loop
repeat
{
solve;
#assign parameters for next iteration
let stop_criteria := abs(linear_approx optimal_val_old);
let optimal_val_old := linear_approx;
let count := count + 1;
let {p in pencils, b in beams} curr_sol[p,b] :=
max(0,x[p,b] .val);
} until ((count >= M) or (stop_criteria <= EPS));
#Calculate the cardinality of positive beams in the mss
#for the current ALPHA and BETA
let curr_card := 0;
for {b in beams}
C
if card{p in pencils : x[p,b] > 0} > 0 then
{
let curr_card := curr_card + 1;
};
};
#write results to file
printf 11____________________________\n" > slaresults.txt;
display sqrt(card(pixels)),card(beams) slaresults.txt;
display count slaresults.txt;
display ALPHA, BETA slaresults.txt;
display curr_card slaresults.txt;
display gamma, optimal_val_old slaresults.txt;
68
display x slaresults.txt;
#reset param for next problem
let optimal_val_old := Infinity;
let {p in pencils, b in beams} sla_sol[p,b] := x[p,b].val;
let sla_gamma := gamma;
let {p in pencils, b in beams} curr_sol[p,b] := 0;
let {p in pencils, b in beams} x[p,b] := 0;
let count2 := count2 + count;
let count := 0;
let gamma := 0;
#
69
B.3.3 The Parameter Selection Script
#Command file: rtp_sla_param.run
tfProgrammer: Janine Kennedy
#Date of last alteration: May 8, 2000
#
option display_lcol Infinity;
option omit_zero_rows 1;
option omit_zero_cols 1;
#Solve the original lp problem to determine bounds
#on ALPHA and BETA
solve;
#
#Calculate the cardinality of positive beams in the
#unpenalized solution
for {b in beams}
{
if card{p in pencils : x[p,b] > 0} > 0 then
C
let curr_card := curr.card + 1;
};
#Assign BETA based on error tolerance for gamma and
#nonzero beams.
let BETA_min := GAMMA_ERROR/(curr_card MAX.BEAMS);
let curr.caxd := 0;
#
#Assign param values for calculating ALPHA start point
for {b in beams}
{
if X < sum{p in pencils} x[p,b] then
{
let X:= sum{p in pencils} x[p,b];
};
};
let X:= X/caxd(pencils);
70
let ALPHA:= 1/X;
#
#Initialize parameters for ASLA
let BETA := 100;
let BETA_opt := BETA;
let ALPHA_opt := alpha;
#Increase ALPHA loop
for {1 .. M>
{
#SLA Loop
include rtp_sla2.run
#Test to continue ALPHA loop
if curr_card <= MAX_BEAMS and sla_gamma <= gamma_alpha then
{
let gamma_diff := gamma,alpha sla_gamma;
let opt_card := curr_card;
let ALPHA.opt := ALPHA;
let ALPHA := ALPHA + ALPHA.INC;
let {p in pencils, b in beams} alpha_sol[p,b] :=
sla_sol[p,b];
let optimal_val_alpha := linear_approx;
let gamma_alpha:= sla_gamma;
}
else let ALPHA := ALPHA + ALPHA,INC;
#Test to determine if gamma produced by SLA with
#current iteration meets clinical requirements
if curr_card <= MAX_BEAMS and sla_gamma <= GAMMA_UP then break;
>;
#Assign parameter values for BETA loop based on ALPHA
#loop results
let ALPHA := ALPHA.opt;
let BETA := BETA.min;
let optimal_val_beta := optimal_val_alpha;
let gamma_beta := gamma_alpha;
71
let {p in pencils, b in beams} beta_sol[p,b] := alpha_sol[p,b];
#Decrease BETA loop
for{l .. M}
{
#SLA loop
include rtp_sla2.run;
#Test to continue BETA loop
if curr_card <= MAX_BEAMS and sla_gamma <= gamma_beta then
{
let opt_card := curr.card;
let BETA_opt := BETA;
let {p in pencils, b in beams} beta_sol[p,b] :=
sla_sol[p,b];
let gamma_beta := sla_gamma;
break;
}
else let BETA := BETA + BETA.INC;
};
#write results to file
printf 11______________________________\n" > results.txt;
display sqrt(card(pixels)).card(beams) results.txt;
display count2 results.txt;
display ALPHA_opt, BETA_opt results.txt;
display opt_card results.txt;
display gamma_beta results.txt;
display beta_sol results.txt;
#calculate total dose distributed to each pixel and write to file
let {i in tumor union bod union critical}
dose[i] := sum{p in pencils, b in beams} d[i,p,b]*beta_sol[p,b];
display dose > dose.txt;
#
72
Appendix C. Dose Depostion Operator C++ Code
//PROGRAM: D0SE2.CPP
//DATE OF LAST MODIFICATION: MARCH 28,2000
//JANINE KENNEDY
//HEADER FILES
#include
#include
#include
#include
#include
#include
#include
#include
//CONST DECLARATIONS
const double PI = 3.14159265358979323846;
//#ifndef PI
//#define PI = 3.14159265358979323846
//#endif
const long unsigned int MAX_STRING_LEN = 100;
const long unsigned int SLICES = 1;
const double EPS = pow(10,(6));
const double DIVERGE.ANGLE = .19739556;
const double SOURCE_LEN = 100;
const double X_DIST = 50;
const double Y_DIST = 50;
#include "proto.h" //define function prototypes
//
//MAIN FUNCTION
int main()
{
//VARIABLE DECLARATIONS FOR MAIN
73
unsigned long int dim = 0;
unsigned long int slices = SLICES;
unsigned long int num_pencils = 0;
unsigned long int num_beams = 0;
double alpha_penalty = 0;//smoothing parameter
int max_row = 10;
int beta_penalty = 0;//penalty parameter
int int_up_bound = 0;//IP ineq. constraint bounds
double int_low_bound = 0;//IP ineq. constraint bounds
double x_dist = X_DIST;//horizontal image size
double y_dist = Y_DIST;//vertical image size
double z_dist = 0.0;
double phi = PI/2;
double theta = 0.0; //gantry angle
double intensity = 1; //energy level
double mu_air = .0278;// mu depends on intensity
double mu_water = .0309;
double mu_bone = .0295;
double mu_muscle = .0306;
int tissue = 0;
double I =0.0; //pencil intensity
double I_old = 0.0; //previos pencil intensity
double alpha = 0.0; //pencil angle
double alpha_o = 0; //first pencil angle
double distance = 0; //distance from source
double s = 0; //distance between intersections
double x_hi = 0; //vertical pixel boundary
double x_lo = 0; //vertical pixel boundary
double y_hi = 0; //horizontal pixel boundary
double y_lo = 0; //horizontal pixel boundary
double mu = 0; //mu value for current pixel
Coord rad_hit = {0.0,0.0,0.0};// intersection pt.
Coord prev_rad_hit = {0.0,0.0,0.0};// intersection pt.
Coord mid_point = {0.0,0.0,0.0};
Coord isocntr = {0.0,0.0,0.0};
Coord source_pt = {0.0,0.0,0.0};
Pixel pixel = {0,0};
ofstream fp; //pointer to output file
ifstream ip; //pointer to input file
74
long ip_inc; //ptr inc for rand, access file
chax filename[MAX_STRING_LEN] ; //output filename
char mu_file[MAX_STRING_LEN]; //name of mu file from matlab
//READ USER MODEL PARAMETERS
user_input(dim, num_pencils,num_beams, beta_penalty,alpha_penalty,
int_up_bound,int_low_bound,filename,mu_file);
//INITIALIZE USER DEPENDENT DATA
double inc = DIVERGE_ANGLE/(num_pencils/2);
double theta_inc = 2*PI/num_beams;
double x_scale[dim + 1]; //image vertical grid lines
double y_scale[dim + 1]; //image horizontal grid lines
//START TIMER
clock_t start_time, end_time;
float comp_time;
start.time = clockO;
//DEFINES THE X AND Y BOUNDARIES FOR EACH PIXEL AND STORES
//IN X.SCALE AND Y.SCALE. USED IN DETERMINING PENCIL INTERSECTIONS
//WITH PIXELS IN MEDICAL IMAGE.
def_pixel_borders(x_scale,y_scale,dim,x_dist,y_dist);
//FOR LOOP CALCULATES DOSES TO EACH PIXEL A BEAM AT A TIME
//AND STORES TO THE DOSE.DEPOSTION ARRAY
//
for(int beam = 0; beam < num_beams; beam++)
{
//WRITE SET SIMPLE PARAMETER DATA TO DISK IN AMPL FORMAT
if(!beam)
C
if ( (ampl_data(f ilename,dim,beam,num_beams ,max_row,
num_pencils,beta_penalty,alpha_penalty,
int_up_bound, int_low_bound)))
exit(O) ;
>
75
//CALCULATE THETA AND SOURCE POINT FOR CURRENT BEAM
theta = beam*theta_inc; //angle for the current beam
alpha_o = theta + PI + DIVERGE_ANGLE + inc/2;
source_pt = error(calc_src(SOURCE_LEN, theta, phi));
//LOOP CALCULATES PENCIL/PIXEL INTERSECTIONS FOR CURRENT BEAM
for(int pen = 1; pen <= num_pencils; pen++)
i
alpha = alpha_o pen*inc; //adjust angle for next pencil
for( int pix =1; pix <=2*dim; pix++)
if(pix == 1)
{
assign_coord_info(rad_hit, x_dist/2,x_dist/2,
y_dist/2,y_dist/2, source_pt, alpha, distance);
//CALCULATE INTENSITY AT FIRST INTERSECTION WITH IMAGE
s = distance;
I = intensity*exp(mu_air*s);
>
else
{
assign_coord_info(rad_hit, x_hi, x_lo,
y_hi, y_lo, source_pt, alpha, distance);
mid.point = scal_mult(.5,vec_add(prev_rad_hit,rad_hit));
find_pixel_location(pixel,x_scale,y_scale,
mid_point,dim);
//CALCULATE THE CURRENT RADIATION INTENSITY
s = norm(vec_sub(prev_rad_hit,rad_hit));
if(!s) break; //pencil exits image
I = I_old*exp(mu*s);
//WRITE CURRENT PIXEL DOSE TO DISK IN AMPL FORMAT
fp.open(filename,ios::app);
76
if (!fp)
{
cerr"ERROR: Could not open output file.\n"
exit(O);
>
fp.setf(ios::showpoint);
fp setprecision(O);
fp "let d[\"" pixel.row
fp pixel.col "V'.V'p" pen "\",\"b";
fp beam "\"] := I_old I ";\n";
fp.closeO ;
>
I_old = I; //reset I_old for next iteration
//ASSIGN MU FOR NEXT ITERATION
if (pix != 1)
{
ip.open(mu_file,ios::in); //open mu file as read only
if (! ip)
{
cout "ERROR: Could not open input file.\n";
exit(0) ;
>
ip_inc = 2*((pixel.row l)*dim + pixel.col 1);
ip.seekg(ip_inc,ios::beg); //moves ptr to next value
ip tissue; //assign new tissue type
ip.closeO; //close file
switch(tissue)
{
case (1): {mu
case (4): {mu
case (2): {mu
case (3): {mu
}
mu_air; break;}
mu_water; break;}
mu.bone; break;}
mu_muscle; break;}
77
>
//ASSIGN PIXEL BOUNDARIES FOR NEXT INTERSECTION
get_bounds(x_scale, y.scale, dim, x_hi,
x_lo, y_hi, y_lo, rad_hit);
prev_rad_hit = rad_hit; //store rad_hit for next iter.
} //END PIXEL FOR LOOP
distance = 0;
pixel.row = 0;
pixel.col = 0;
mid_point.x_coord = 0;
mid_point.y_coord = 0;
> //END PENCIL FOR LOOP
> //END OF BEAM FOR LOOP
//REPORT TIME RESULTS
end_time = clockO;
comp_time = (double)(end_time start_time)/ CLOCKS_PER_SEC;
cout "Calculation time = comp_time
seconds\n" endl;
return 0;
}
//
//FUNCTIONS
//DISPLAY_C00RD: Uses cout to display the elements of
//the struct Coord
void display_coord(Coord u)
{
cout u.x_coord "
cout u.y_coord "
78
cout u.z_coord endl;
>
//DISPLAY_PIXEL: Uses cout to display the elements of
//the struct Coord
void display.pixel(Pixel u)
{
cout u.row 11
cout u.col endl;
>
//CALC.SOURCE: Converts the radiation source pt from
//Spherical to Cartesian coordinates
//P = distance from the isocenter
//theta = gantry angle in radians
//phi = bed angle in radians
Coord calc.src(double P, double theta, double phi)
{
Coord my.srcpt = {0,0,0};
my.srcpt.x.coord = P*sin(phi)*cos(theta) ;
my.srcpt.y.coord = P*sin(phi)*sin(theta) ;
//my.srcpt.z.coord = P*cos(phi);
error(my.srcpt);
return my.srcpt;
}
//VECTOR ADDITION
Coord vec.add(Coord a, Coord b)
{
Coord v;
v.x.coord = a.x.coord + b.x.coord;
v.y.coord = a.y.coord + b.y.coord;
v.z.coord = a.z.coord + b.z.coord;
return v;
}
79
//VECTOR SUBTRACTION
Coord vec_sub(Coord a, Coord b)
C
Coord v;
v.x_coord = a.x_coord  b.x_coord;
v.y_coord = a.y_coord  b.y_coord;
v.z_coord = a.z_coord  b.z_coord;
return v;
}
//DOT PRODUCT
double dot_prod(Coord a, Coord b)
{
double d_prod = (a.x_coord b.x_coord) +
(a.y_coord b.y_coord) +
(a.z_coord b.z_coord);
if (fabs (d_prod) <= EPS)
d_prod = 0;
return d_prod;
}
//VECTOR NORM
double norm(Coord a)
{
double norm_a = 0.0;
norm_a = sqrt(pow(a.x_coord,2) +
pow(a.y_coord,2) +
pow(a.z_coord,2));
return norm_a;
>
//SCALAR VECTOR
Coord scal_mult(double scale, Coord v)
C
Coord new_v;
new_v.x_coord = scale v.x_coord;
new_v.y_coord = scale v.y.coord;
new_v.z_coord = scale v.z.coord;
80
return new_v;
>
//Error: Rounds components to zero when vector components
//are less then EPS
Coord error(Coord a)
{
if((fabs(a.x_coord)) <= EPS)
a.x.coord = 0.0;
if((fabs(a.y_coord)) <= EPS)
a.y_coord = 0.0;
if((fabs(a.z_coord)) <= EPS)
a.z_coord = 0.0;
return a;
>
//ANGLE_BETWEEN_VECT0RS:Calculates the angle between two vectors.
double angle_between_vectors(Coord a, Coord b)
{
double angle = 0;
angle = acos(dot_prod(a,b)/(norm(a)*norm(b)));
return angle;
>
//ASSIGN_C00RD_INF0:Calculates points of intersection for each
//pencil with the x and y boundaries of the pixels.
void assign_coord_info(Coord ferad.hit, double x_hi, double x_lo,
double y_hi, double y_lo, Coord source.pt, double angle,
double ftdistance)
{
double r = 0;
Coord rad_hit.fixed = {0.0,0.0,0.0};
double old_dist_to_source = distance;
double temp_dist_to_source = 0;
81
double current_dist_to_source = 400;
//Fix x and find parameter r_x for the equation of the line
//from the prev_rad_hit to the new rad_hit
r = (x_lo source_pt.x_coord)/cos(angle);
rad_hit_fixed.x_coord = x_lo;
rad_hit_fixed.y_coord = source_pt.y_coord + r*sin(angle);
rad_hit_fixed = error( rad_hit_fixed);
temp_dist_to_source = norm(vec_sub(rad_hit_fixed,source_pt));
compare_rad_hit(rad_hit_fixed, rad_hit,temp_dist_to_source,
current_dist_to_source, old_dist_to_source, X_DIST, Y_DIST);
r = (x_hi source_pt.x_coord)/cos(angle);
rad_hit_fixed.x_coord = x_hi;
rad_hit_fixed.y_coord = source_pt.y_coord + r*sin(angle);
rad_hit_fixed = error( rad_hit_fixed);
temp_dist_to_source = norm(vec_sub(rad_hit_fixed,source_pt));
compare_rad_hit(rad_hit_fixed, rad_hit,temp_dist_to_source,
current_dist_to_source, old_dist_to_source, X_DIST, Y_DIST);
//Fix y and find parameter r_y for the equation of the line
//from the prev_rad_hit to the new rad_hit
r = (y_lo source_pt.y_coord)/sin(angle);
rad_hit_fixed.y_coord = y_lo;
rad_hit_fixed.x_coord = source_pt.x_coord + r*cos(angle);
rad_hit_fixed = error(rad_hit_fixed);
temp_dist_to_source = norm(vec_sub(rad_hit_fixed,source_pt));
compare_rad_hit(rad_hit_fixed, rad_hit,temp_dist_to_source,
current_dist_to_source, old_dist_to_source, X_DIST, Y_DIST);
r = (y_hi source_pt.y_coord)/sin(angle);
rad_hit_fixed.y_coord = y_hi;
rad_hit_fixed.x_coord = source_pt.x_coord + r*cos(angle);
rad_hit_fixed = error(rad_hit_fixed);
temp_dist_to_source = norm(vec_sub(rad_hit_fixed,source_pt));
compare_rad_hit(rad_hit_fixed, rad_hit,temp_dist_to_source,
current_dist_to_source, old_dist_to_source, X_DIST, Y_DIST);
82
distance = current_dist_to_source;
>
//COMPARE_RAD_HIT: Chooses the pixel intersection that gives the
//shortest path through a pixel
void compare_rad_hit(Coord &rad_hit_fixed, Coord &rad_hit,
double &temp_dist, double ¤t_dist, double old_dist,
double x_dist, double y.dist)
{
if (fabs(rad_hit_fixed.x_coord) <= x_dist/2)
{
if (fabs(rad_hit_fixed.y.coord) <= y_dist/2)
{
if(temp_dist > old_dist)
{
if(temp_dist <= current_dist)
{
rad_hit.x_coord = rad_hit_fixed.x_coord;
rad_hit.y_coord = rad_hit_fixed.y_coord;
current_dist = temp_dist;
rad_hit_fixed.x_coord = 0;
rad_hit_fixed.y_coord = 0;
>
>
>
>
//Reset temp variables to zero
rad_hit_fixed.x_coord = 0;
rad_hit_fixed.y_coord = 0;
temp_dist = 0;
>
//Defines the x and y boundary lines for each pixel in the
//medical image and stores in x_scale and y.scale arrays from
// main func.
void def_pixel_borders(double *x_scale, double *y_scale, int dim,
double x_dist, double y_dist)
83
{
double x_inc =
double y_inc =
double x_start
double y_start
fabs(x_dist/dim);
fabs(y_dist/dim);
= (x_dist/2);
= (y_dist/2);
for(int i = 0; i <= dim; i++)
{
*(y_scale + i) = (y_start + i*y_inc);
*(x_scale + i) = (x_start + i*x_inc);
}
}
//FIND_PIXEL_L0CATI0N: Takes the mid.point of the line
//passing through prev_rad_hit and rad_hit in main
//and places it within a pixel of the medical image.
//The pixel is determined uniquely by assigning a matrix
//row and column number. Where pixel(i,j) gives
//the ith row and the jth column in the standarad way.
void find_pixel_location(Pixel &pixel_location,
double *x_scale, double *y_scale,
Coord mid_point, int dim)
{
pixel_location.row = 0;
pixel_location.col = 0;
for(int i = 0; i < dim; i++)
{
//Use the x element of mid.point to determine
//matrix column
if (*(x_scale + i) <= mid_point.x_coord )
{
if (mid.point.x_coord< *(x_scale+i+l))
pixel.location.col = i+1;
>
if (mid_point.x_coord == double(*(x_scale+dim)))
pixel_location.col = dim;
//Use the y element of mid_point to determine
//matrix row
84
if (*(y_scale + dim i 1) < mid_point.y_coord)
{
if (mid_point.y_coord < *(y_scale+dimi))
pixel_location.row = i+1;
}
if (mid_point.y_coord == double(*(y_scale + dim)))
pixel_location.row = 1;
}
}
//GET_BOUNDS: There cure four sides to each pixel. As a pencil
//passes through adjacent pixels it intersects the boundaries.
//This function determines which pixel boundaries may be
//intersected by next by a pencil.
void get_bounds(double *x_scale, double *y_scale, int dim,
double &x_hi, double &x_lo, double &y_hi,
double &y_lo, Coord rad_hit)
C
//ASSIGN VERTICAL BOUNDARIES
if (*(x_scale ) == rad_hit.x_coord)
C
x_lo = rad_hit.x_coord;
x_hi = *(x_scale + 1);
>
else if (*(x_scale + dim) == rad.hit.x_coord)
{
x_lo = *(x_scale + dim 1);
x_hi = rad_hit.x_coord;
>
else
{
for(int i=0;i < dim; i++)
{
if(*(x_scale + i) < rad_hit.x.coord)
C
if(*(x_scale + i+1) > rad.hit.x_coord)
{
85
x_lo
= *(x_scale + i);
x_hi = *(x_scale + i + 1);
>
else if (*(x_scale + i+1) == rad_hit.x_coord)
C
x_lo = *(x_scale + i);
x_hi = *(x_scale + i + 2);
}
}
>
}
//ASSIGN HORIZONTAL BOUNDARIES
if (*(y_scale ) == rad_hit.y_coord)
{
y_lo = rad_hit.y_coord;
y_hi = *(y_scale + 1);
>
else if (*(y_scale + dim) == rad_hit.y_coord)
{
y_lo = *(y_scale + dim 1);
y_hi = rad_hit.y_coord;
}
else
{
for(int i=0;i < dim ; i++)
{
if(*(y_scale + i) < rad_hit.y_coord)
C
if(*(y_scale + i+1) > rad_hit.y_coord)
C
y_lo = *(y_scale + i);
y_hi = *(y_scale + i + 1);
>
else if (*(y_scale + i+1) == rad_hit.y_coord)
C
86
y_lo = *(y_scale + i);
y_hi = *(y_scale + i + 2);
>
>
}
>
}
//AMPL_DATA: Write and format the set and parameter data
//needed in the rtp mip model to an AMPL readable file.
//This function doesnt manipulate the dose data.
int ampl_data(char filename[MAX_STRING_LEN], int dim, int beam,
int num_beams, int max_row_item, int num_pencils,
int beta_penalty, double alpha_penalty,
int int_up_bound, double int_low_bound)
{
ofstream fp; //pointer to output file
fp.open(filename,ios:rout);
if (!fp)
{
cerr"ERROR: Could not open file.Xn";
return 0;
>
fp.setf(ios::showpoint);
fp setprecision(O);
fp "set pencils :=
for(int i=l; i<=num_pencils; i = i + max_row_item)
{
for(int k=0; k
{
if(i+k < num_pencils)
fp "p" << i+k
else if (i+k == num_pencils)
fp "p" << i+k
else
break;
87
>
fp endl;
}
fp "set beams := ";
for(int i=0; i < num_beams; i = i + max_row_item)
C
for(int k=0; k
if(i+k < num_beams 1)
fp "b" i+k
else if (i+k == mun_beams 1)
fp "b" i+k
else
break;
}
fp endl;
>
fp "set pixels :=";
int count = 0;
for(int i=l; i <= dim; i++)
{
for(int j=l; j <= dim; j++)
C
if(i*j == dim*dim)
fp i j ";\n";
else if (count < max_row_item)
fp i j
else
C
fp i j ",\n";
count = 0;
>
count++;
}
}
fp.closeO ;
return 1;
88
>
//USER_INPUT: Uses references to store user input in
//main variables
void user_input(unsigned long int &dim,
unsigned long int &num_pencils,
unsigned long int &num_beams,int &beta_penalty, double
&alpha_penalty, int &int_up_bound, double &int_low_bound,
char filename[MAX_STRING_LEN] char mufile[MAX_STRING_LEN])
C
cout "DIM: "; cin dim;
cout "PENCILS (per beam): "; cin num_pencils;
cout "BEAMS (total in model): cin num_beams;
cout "input MU FILENAME (from MATLAB): ";
cin muf ile;
cout "output FILENAME (an AMPLA file): ";
cin filename;
>
89
References
[1] Dimitris Bertsimas and John N. Tsitsiklis. Introduction to Linear Opti
mization. Athena Scientific, Belmont, Massachusetts, 1997.
[2] Michael C. Ferris and O.L. Magasarian Linear programming with matlab.
Unpublished, May 1997.
[3] Robert Fourer, David M. Gay, and Brian W. Kernighan. AMPL: A Model
ing Language For Mathematical Programming. Boyd and Fraser publishing
company, Danvers, Massachusetts, 1993.
[4] ILOG, France. ILOG CPLEX 6.5, 1999.
[5] O.L. Magasarian. Machine learning via polyhedral concave minimization.
In H. Fischer, B. Riedmueller, and S. Schaeffler, editors, Applied Mathemat
ics and Parallel Computing Festschrift for Klaus Ritter, pages 175188.
PhysicaVerlag A SpringerVerlag Company, Heidelberg, 1996.
[6] O.L. Magasarian. Minimumsupport solutions of polyhedral concave pro
grams. Optimization., 45:149162, 1999.
[7] Francis Newman. A linear program for beam angle and intensity optimiza
tion in the radiotherapy treatment of tumors. Unpublished, May 1997.
[8] R.T. Rockafeller. Convex Analysis. Princeton University Press, Princeton,
New Jersey, 1970.
90
