mcm优秀论文

Forward and Inverse Treatment Planning Analysis of the Problem
Gamma Knife is an advanced treatment unit that supplies stereotactic approach to the treatment of lesions within the head such as tumors and malformations by delivering a single dose of ionizing radiation, which is referred to as a "shot". Every shot emanates from 201 cobalt-60 unit sources in a shielded treatment unit. The 201 beams simultaneously intersect at the same point in space, resulting in a spherical region of high dose. Four kinds of shots with beam widths of 4,8,14,and 18 mm are available for irradiating different size volumes. By combining several shots, relatively larger target volume can be treated with the Gamma Knife. Conventionally, the treatment plan is figured out by experienced neurosurgeons and radiation oncologists, but this empirical plan may be feasible at the price of some damage of critical or normal tissues. So our goal is to automate the treatment planning process for lesions of various sizes and shapes. Namely, given a tumor of some shape and size, we have to determine: the number of shots, the location of each shot, and the diameter of each shot under some constraints to produce a dose distribution that conforms closely to the target volume and confines the dose to the normal tissues to a maximum extent.
Assumptions
●Gamma Beam is perfectly cylindrical.
●Each individual gamma beam is of relatively weak intensity, and therefore
has no affect on any tissues.
●Tumor is treated with Gamma Knife in our models.
Definitions of Constants and Terms
●target volume: tumor;
●normal tissues: fine tissues but not dose-sensitive;
●critical areas: dose-sensitive tissues.
●X i, Y i ,Z i : the coordinates of isocenter;
●R i: the diameter of shot i
Set Parameters Model
In this model, we try to produce treatment plans simply under the constraints:
no shot protrudes outside the target,
no overlapping between any two shots, and
cover at least 90% of the target volume by all shots.
Based on the discussion in the analysis of the problem, we define the treatment plan as a function of the coordinates of isocenter of each shot and associated diameters:
S (x1 , y1 , z1 , R1,……x i , y i , z i , R i,……),
where x i , y i , z i is the coordinates of isocenter and R i is the diameter of shot i . To formulate optimal treatment plan under the constraints mentioned above is to decide the parameters included in S . Thus, the whole problem is reduced to a sphere-packing problem, and given a tumor, we employ simulated annealing algorithm to produce the treatment plan S (see Simulated Annealing Algorithm). Although there are several algorithms to deal with the sphere-packing problem, the simulated annealing algorithm is relatively fast, which is appreciated by the users of Gamma Knife. Afterwards, with the fixed parameters in S, we sketch the dose distribution to see whether or not it meets clinical demands. If so, this treatment plan is fit for the present tumor. Clinical Demands
●Constraints
To protect the normal tissues from lethal dose delivered by all shots, there is
one constraint:
keep the dose to normal tissue below tolerance doses.
For each treatment plan S , we get the different isodose curves. Given a critical tissue, since the tolerance dose is fixed, we use the corresponding isodose curve to evaluate the plan. Assuming the tolerance dose of some tissue is 50%, if the area within the 50% isodose curve includes any part of the tissue, this plan is unsatisfactory. (Both satisfactory-plan and unsatisfactory-plan are shown in Figure 1.)
● Optimization Demands
To guarantee safe and accurate treatment, several optimization demands are as follows:
● homogeneity: the doses to points in the target show no significant
difference ,
● minimal dose to the entire normal tissues, and
●
minimal dose delivered to points in critical areas. Since the points of some distance from the isocenter are hardly affected by the shot, we encompass the tumor with a cube, whose size can be determined empirically, and only consider the doses to points within the cube. Assuming that target volume, normal tissues and critical areas are all included in the cube, uniformly we pick N a points, N b points and N c points respectively in the three parts, respectively . We measure the total dose by:
∑=Na la la la la z y x
U 1),,( in the target volume,
),,(∑b b N l lb lb lb z y x
U in the normal tissues, and
),,(lc lc Nc lc lc z y x
U ∑ in the critical areas,
where U (x l , y l , z l ) is the total dose delivered to a point by all shots in some treatment plan.
∑==N
i i z y x D z y x U 1),,(),,(,
where N is the number of shots;
∑=⨯⨯⨯⨯=201
12),,(),()(),,(m m m m OAR m m TMR m OF m F d R r U d R U R U C z y x D
In the target volume, the maximal dose of one point is
),,(max 1la la la N l z y x U a a ≤≤ and the minimal is ),,(min 1la la la N l z y x U a a ≤≤. We measure the dose gradient by the difference:
),,(min ),,(max )(11la la la N l la la la N l i z y x U z y x U s U a a a a ≤≤≤≤-=∆
The reason why we use ∆U :
Since the points is representative of the dose distribution, it is reasonable to believe that the maximum dose of the whole target volume is close to ),,(max 1la la la N l z y x U a a ≤≤, and so is it with the minimum dose. To decrease the dose gradient across the target volume is to make the dose to each point is close as much as possible, so minimizing the difference between the maximum and the minimum dose decrease the gradient to the largest extent.
As for points in critical areas every dose of point is strictly constrained so that we have to minimize the maximal dose of all these points. The maximal dose of any point in the critical areas is:
),,(max 1lc lc lc N l z y x U c
a ≤≤
Therefore, for N satisfactory-plans:
S i , i =1…N
we use the following optimization form:
⎪⎪⎩
⎪⎪⎨⎧-=∆≤≤=≤≤≤≤∑)),,)((max min()
,,)((max ),,(min ),,(max )(min 1111c c c c c a a a a l l l i N l Na la la la la i la la la N l la la la N l i z y x S U z y x S U z y x U z y x U s U s.t. }{N i S S S ..1∈
where ∆U (S i ) is the ∆U of the i th treatment among the N plans
Generally speaking, no treatment plan S i would satisfy the above three objective functions simultaneously. So we cannot but produce a relative optimal treatment plan from the N plans by assigning three weights α,β,γ to the objective functions. Thus, the single objective function:
),,)((max ),,)(()(min 11c
c c c c a a a a a l l l i N l N l l l l i i z y x S U z y x S U S U ≤≤=⨯+⨯-∆⨯∑γβα
is used to substitute for the multiple objective functions. α,β,γ,differing in importance, are determined with relatively weights on the three objective function, and, obviously, α+β+γ=1.
Simulated Annealing Algorithm
To produce a treatment plan is to determine the parameters of S. Here, we treat the process of setting parameters of plan as a sphere-packing problem. Given a tumor within a brain of certain shape and size, on the surface of the tumor we uniformly pick evenly N sample points from the digital image that consists of millions of points:
)z , y , (x i i i 1N
i P ≤≤
From the shape and size of the tumor, the number N s of shots that will be delivered, along with the coordinates of isocenter and the diameter of the shot:
S (1111,,,R z y x …Ns Ns Ns Ns R z y x ,,,),
Where x j , y j , z j are the coordinates of isocenter and R j is the diameter of ith shot in the treatment plan S. We build up a corresponding relationship between sample points P and shots S by the principle that the distance between the point P i and the sphere that represents shot j is the shortest among all the distances from P i to the sphere of other shots. Put simply in mathematics, our objective function in the simulated annealing algorithm is:
F (S )= (∑=N
i j i d 12),() / N
j j i j i j i R z z y y x x sqrt j i d --+-+-=])()()[(),(222,
where ),(j i d is the distance from point P i to the sphere of its corresponding shot.
1. With estimated number of shots, select one possible treatment plan S 0 at
random, and let the ―temperature‖ T , a parameter in the cooling schedule of the algorithm, high enough.
2. Alter S 0 to get another plan S 1 (see How to Get S 1); calculate the objective
function with respect to S 0 and S 1, and get the difference: )()(01S F S F F -=∆
3. If F ∆<0, then substitute S 1 to S 0; otherwise when the probability Exp
{-F ∆/T } is greater than a random value between 0 and 1, the substitution is also necessary.
4. Repeat the step 2 and 3 until the process produces a relatively optimal
treatment plan. In practice, we get times L for repetition beforehand.
5. Each time reduce the temperature T by ()T λ-1, where λ is a constant
that we assume is 0.9. If the updated temperature T is greater than some certain value, go back to step 2; otherwise terminate the whole process and
thus we get the global optimal treatment plan.
6. If the percentage of the coverage of the target is less than 90%, we can
increase the number of shots by one and go back to step 1.
Besides, with the same process, we introduce another objection function:
∑∑==--+-+-=s
N j N
i j j i j i j i R z z y y x x sqrt S G 11
222)))()()((()(
How to Get S 1
To get a new treatment plan S 1 based on S 0, we need to alter some or all parameters of S 0. The parameters include the coordinates of isocenter, diameter of each shot. Each time we either increase or decrease coordinates by 0.1mm, whereas change the diameter by 4mm. However, to make sure that any two shots would not overlap, we calculate the distance ),(j i D between any two isocenters. If ),(j i D >R i +R j , the S 1 is obtained.
Set the Desired Dose Distribution Model
When the constraints are strict, the process of ―trial and err‖ of the ―Set Parameter Model‖ would not work efficiently. For example, when the critical areas are close to the target volume, the process is time consuming. So, in this case, we cannot set the parameters first and measure the treatment plan. Instead, we may as well turn back on the dose distribution to determine the parameters. According to concrete case, we may consider one of the functions below as objective function.
Dose to each point :
2,,11
)(),()(),,(F d s r U d s U s U C z y x D j j j OAK Ns i j TMK N
j OFJ j p t
⨯⨯⨯⨯=∑∑==
N s is the number of shots , N t =201
])sin )(sin cos )(cos cos )(()()()[(2222θβθβϕj j j j j j p y y z z x x z z y y x x sqrt r -+-+---+-+-=θβθβθsin )(cos cos )(cos cos )(j j j p y y z z x x SSD A d -+-+---=
Integral dose to tumor: ∑=a
N i i
i i z y x f 1),,( Maximum dose to tumor: )),,((max 1i i i N
i z y x f ≤≤ Minimum dose to tumor: )),,((min 1i i i N
i z y x f ≤≤ Maximum dose to sensitive tissue or organs: ),,((max 1i i i N i z y x f c
≤≤ Maximum dose to critical volumes: ),,((max 1i i i N i z y x f d
≤≤ Integral dose to the entire volume of normal tissue or organs:
∑=b N i i i i z y x f 1),,(
Maximum dose to specified normal tissue points:
),,((max 1i i i Ne i z y x f ≤≤
⎪⎪⎩
⎪⎪⎨⎧-≤≤=≤≤≤≤∑)),,((max min )
,,(min ))],,((min ),,((max min[1111i i i N i N i i i i i i i N i i i i N i z y x f z y x f z y x f z y x f d b s.t. :
C z y x f i i i Ne i ≤≤≤),,((max 1
In fact, this is a multi objective function optimize problem. We give each objective function weights γβα,,, respectively. We also employ simulated
annealing algorithm to find the optimal solution.
Test
We write a C++ program to simulate the implementation of our algorithm. For simple manipulation, we suppose that the tumor is egg-shaped in the test, and
the tolerance dose of normal tissues is 50%. We consider two different situations in terms of critical areas, shown in Figure 2. The Figure 2 (a) is the situation where critical areas are almost out of the cube encompassing the tumor, whereas Figure 2 (b) is the situation where critical areas are close to
the target volume.
For situation shown in Figure 2 (a), by ―Set Parameters Model‖, we get a set
of treatment plans shown in Table 1:
Treatment
Coordinates of isocenter and diamater
plan
1 (59.58,64.01,50.10,14) (59.73,46.70,61.70,14) (41.24,55.26,58.06,14) (52.93,52.93,40.53,18)
2 (48.32,50.32,43.32,14) (58.14,60.07,31.84,14) (56.46,58.44,51.44,14) (58.96,60.96,58.96, 8)
3 (53.16,61.78,58.16,14) (39.87,37.87,40.87,18) (65.84,43.83,49.33,14) (38.75,55.75,49.76,18)
4 (47.60,49.65,42.50,14) (61.36,63.40,60.20,14) (55.70,57.70,50.70,14) (62.20,44.20,57.50,18)
5 (52.32,54.29,47.32,14) (46.82,50.05,50.00,14) (42.70,44.80,39.70, 8) (52.32,54.29,47.32,14)
6 (47.56,49.56,42.53,14) (61.46,63.40,37.22,14) (55.74,57.71,50.72,14) (62.20,64.23,57.20, 8)
7 (59.40,49.40,61.40,14) (51.20,54.72,50.62,18) (41.62,40.62,61.42,18) (45.95,45.95,33.55,14)
Although the all ten-treatment plan is obtained through a process of trail and err, we need to determine which plan is the relatively optimal plan.
By the single objective function:
),,)((max ),,)(()(min 11c c c c c a a a a a l l l i N l N
l l l l i i z y x S U z y x S U S U ≤≤=⨯+⨯-∆⨯∑γβα
s.t.
{}10...S Si S i ∈
First time, let 4.0,4.0,2.0===γβα, and we get the relatively optimal treatment plan S 10. By calculating the dose of points within the cube, we sketch the isodose curves shown in Figure (3). And in the same way, we get the dose-volume histogram of S 10 shown in Figure (4).
Second time, let 5.0,3.0,2.0===γβα, and we get the relatively optimal treatment plan S 4. So situations with different weights generally lead to different optimal plan.
As the situation shown in Figure 2( b), the program of ―Set Parameters Model‖ becomes time consumin g, and the treatment generated take little care of the critical areas. So, we suggest that when the critical areas are close to the tumor, it may be much safer and fast to adopt the ―Set the Desired Dose Parameters Model‖
Strengths and Weaknesses
●Strengths
The ―Set Parameters Model‖ can give an optimal treatment plan if the tumor is of some distance from the critical areas with high coverage; while the ―Set the Desired Dose Distribution Model‖ guarantees the safety of the critical areas regardless of high coverage.
●Weaknesses
We only test our models on the relatively well-defined volume, so the
algorithm can hardly apply to volumes of other shapes.
References
[1] Y. Sato, N. Shiraga, S. Nakajima, S. Tamura, R. Kikinis, Local maximum intensity projection (LMIP): a new rendering method for vascular visualization, J. Comput. Assisted Tomogr. 22 (6) (1998)
912–917.
[2] Robert J. Sterootactic Radiosurgery Using the 201 Cobalt—60 Source Gamma knift. Neurosurg Clin North Am,1990,1: 933.
[3] D. Vandermeulen, R. Verbreeck, L. Berden, P. Suetents, G. Marchal, Continuous voxel classification by stochastic relaxation: theory and application to MR imaging and MR angiography, Inf. Process. Med. [4] G. Gerig, T. Koller, G. Szekely, C. Brechbuhler, O. Kubler, Symbolic description of 3D structures applied to cerebral vessel tree obtained
from MR angiography volume data, Inf. Process. Med. Imaging Lect. Notes Comput. Sci. 687 (1993) 94–111.
[5] Shu HZ, Yan YL, Bao XD,et al.Three dimensional optimization of
treatment planning for gamma unit treatment system[J] . Med Phy , 1998,25:2352--2357
[6] P.P. Maeder, R.A. Meuli, N. der Tribolet, Three-dimensional volume rendering for magnetic resonance angiography in the screening and preoperative workup of intracranial aneurysms, J. Neurosurg. 85 (6) (1996) 1050–1055.
[7] L.O. Hall, A.M. Bensaid, L.P. Clarke, R.P. Velthuizen, M.S. Silboger, J.C. Bezdek, A comparison of neural network and fuzzy clustering techniques in segmenting magnetic resonance images of the brain,
IEEE Trans. Neural Network 3 (5) (1992) 672–682.
[8] D.G. Health, P.A. Soyer, B.S. Kuszyk, D.F. Bliss, P.S. Calhoun, D.A. Bluemke, M.A. Choti, E.K. Fishman, Three-dimensional spiral CT
during arterial portography: comparison of three rendering technique, Radiographics 15 (4) (1995) 1001–1011.
[9] Lunsford LD, Flickinger J, Lindner G, et al . Stereotactic radiosurgery of the brain using the first United States 201 Cobalt-60 source gamma knife . Neurosurgery , 1989,24:151—158 .
[10] S. Kobashi, N. Kamiura, Y. Hata, Fuzzy information granulation on segmentation of human brain MR images, J. Jpn Soc. Fuzzy Theory . [11] User’s manual of Leksell gamma unit, Elekta Edition,April 1992 .。