Lateral Penumbra Modelling Based Leaf End Shape Optimization for Multileaf Collimator in Radiotherapy.

Lateral penumbra of multileaf collimator plays an important role in radiotherapy treatment planning. Growing evidence has revealed that, for a single-focused multileaf collimator, lateral penumbra width is leaf position dependent and largely attributed to the leaf end shape. In our study, an analytical method for leaf end induced lateral penumbra modelling is formulated using Tangent Secant Theory. Compared with Monte Carlo simulation and ray tracing algorithm, our model serves well the purpose of cost-efficient penumbra evaluation. Leaf ends represented in parametric forms of circular arc, elliptical arc, Bézier curve, and B-spline are implemented. With biobjective function of penumbra mean and variance introduced, genetic algorithm is carried out for approximating the Pareto frontier. Results show that for circular arc leaf end objective function is convex and convergence to optimal solution is guaranteed using gradient based iterative method. It is found that optimal leaf end in the shape of Bézier curve achieves minimal standard deviation, while using B-spline minimum of penumbra mean is obtained. For treatment modalities in clinical application, optimized leaf ends are in close agreement with actual shapes. Taken together, the method that we propose can provide insight into leaf end shape design of multileaf collimator.

1. Introduction

Lateral penumbra of single-focused multileaf collimator has been recognized as one of the key dosimetric characteristics and has a significant impact on dose delivery accuracy in radiation therapy treatment planning [1, 2]. Penumbra characteristics strongly influence the amount of healthy tissue involvement, especially where sharp dose gradient is required for stereotactic body radiation therapy.

Penumbra is typically defined as the distance over which the dose profile falls from 80% to 20% [3]. For clarification, penumbra is classified into two categories, namely, fluence penumbra and dosimetric penumbra. Fluence penumbra is acquired on the in-air scoring plane, and it is composed of geometric penumbra and transmission penumbra. Fluence penumbra is also known as in-air penumbra. In contrast, dosimetric penumbra is measured within phantom, and it is the combined effect of fluence penumbra and phantom scatter factor. Dosimetric penumbra is also referred to as physical penumbra or clinical penumbra. In what follows we refer exclusively to fluence penumbra of single-focused multileaf collimator.

Intensive studies have been carried out to explore the experimental penumbral properties for radiotherapy modalities [4, 5]. Results have shown that lateral penumbra is significantly correlated with source model, leaf position, and leaf end shape. On the other hand, the impact of rounded leaf end effect on radiation field offset has been well understood and quantified. Growing evidence has also revealed that penumbra region is field size dependent and largely attributed to the leaf end shape [6, 7].

Researches on leaf end design have been conducted in the past decades [8]. Based on chord intersection assumption, empirical method for leaf end shape optimization has been utilized for obtaining optimal radius of leaf end in the shape of symmetric circular arc. Analytical expressions of geometric penumbra and transmission penumbra have been derived separately. With focal spot size incorporated, computer simulation has shown that leaf end with slightly smaller radius than the empirical result would be suggested. However, problems remain in terms of modelling accuracy of leaf end induced total penumbra and leaf end shape optimization if geometric parameterization is more complex than symmetric circular arc. This begs the question, “How to get physically the lateral fluence penumbra performance by means of leaf end shape optimization?,” which is the motivation for the study.

In terms of total penumbra modelling [9], full Monte Carlo simulation has been proved to be the most accurate among currently available algorithms. However, detailed radiotherapy geometries should be provided in practice and a large number of particle histories have to be recorded in order to achieve desired outcomes. Alternatively, leaf end correction based ray tracing algorithm was found to be easy to implement and could potentially achieve the desired accuracy, which deploys simplified geometries and virtual source models [10]. Since shape optimization process typically involves a large number of penumbra function evaluations in order to search for global optimum, the difficulty of ray tracing algorithm lies in the fact that computation time and penumbra accuracy are largely dependent on discretization error of radiation field and virtual source.

The framework of our work is primarily made up of four parts. In the first part, for the purpose of fast penumbra evaluation, a novel method for leaf end induced penumbra modelling is presented based on Tangent Secant Theory (TST). Compared with Monte Carlo simulation and ray tracing algorithm, TST penumbra modelling is proved to be cost-efficient. In the second part, a variety of leaf end parameterization techniques are investigated. In the third part, leaf end shape design is formulated as a problem of biobjective optimization with mean-variance cost function. Genetic algorithm based global optimization and gradient based local optimization are deployed. In the last part, radiotherapy geometries of Agility 160-leaf MLC (Elekta AB, Stockholm, Sweden) and Millennium 120 MLC (Varian Medical System, Palo Alto, CA, USA) are utilized to testify the feasibility of TST modelling and leaf end shape optimization method.

2. Materials and Methods

2.1. Tangent Secant Theory for Penumbra Modelling

As illustrated in Figure 1, given an arbitrary point T on the scoring plane, ray tracing algorithm for calculation of beam intensity at T is utilized by summation of weighted beam intensity within the angle of θ. Considering beam penetrating leaf end with path length l and density ρ, emerging beam intensity I is related to incident beam intensity I 0, which is given by the exponential attenuation law, usually referred to as Beer-Lambert law:where μ is the attenuation coefficient and μ /ρ denotes X-ray mass attenuation coefficient.

Figure 1

TST penumbra modelling.

In contrast, the main idea of TST penumbra modelling is that, for arbitrary leaf position, penumbra width is obtained by directly searching for two points, the 80% intensity point and the 20% intensity point, where the 100% intensity point is measured at the centre of radiation field C . Therefore, computational efficiency can be realized, without calculating full-field dose profile.

The modelling method is depicted as follows. Given an arbitrary physical position of leaf end L , the projection of L onto the scoring plane is point T . Penumbra width W is determined by the 80% intensity point P 80 and the 20% intensity point P 20. On the one hand, P 80 is approximated by visible source area integration method, which is utilized by projecting multileaf collimator back to the source plane and summing up the visible area of source profile. Firstly, cumulative distribution function of source profile is expressed as the integral of its probability density function, and the 80% cumulative intensity on source profile is obtained by one-dimension search technique, which denotes point E in Figure 1. Secondly, by drawing tangent line of leaf end from point E, P 80 on the scoring plane is obtained. On the other side, point P 20 is obtained by drawing a secant line from source point C to the scoring plane, with path length defined by chord AB on leaf end curve. C is defined by the equivalent source size length e and AB is defined by the effective path length l. As a consequence, TST penumbra modelling employs only two variables, namely, e and l.

Based on mathematical optimization, iterative approaches for deriving the tangent line and secant line are introduced. For details, refer to Appendices A and B. Intercept theorem is used for penumbra evaluation, and the relations are written with the notation:where the origin of the coordinate system is placed at the centre of circular arc. D denotes point of tangency. Penumbra W is the function of system related vector s and leaf end related vector v. Vector s is composed of the equivalent source size and the effective path length, and vector v is composed of leaf end curve C and leaf end position T. Leaf end curve C is determined by design variables p. For circular arc leaf end, design variables include arc radius and centre offset.

The equivalent source size and the effective path length are obtained by parameter identification, which involves two steps. Firstly, with curves grouped according to design variables, reference data of leaf position related penumbra are obtained by Monte Carlo simulation or ray tracing algorithm. Secondly, nonlinear least squares based curve fitting is introduced to determine the equivalent source size and the effective path length. Mathematical term is as follows:where W and denote penumbra matrices of TST model results and reference data, respectively. W denotes the element of TST penumbra matrix. denotes the element of reference penumbra matrix, which is obtained by Monte Carlo simulation or ray tracing algorithm. C denotes the leaf end curve i, and T denotes the leaf position j.

2.2. Leaf End Shape Parameterization

Topolnjak and van der Heide [8] suggested leaf end designed in the shape of elliptical arc could be beneficial for penumbral properties. Intuitively, polynomial curves are superior in terms of the ability to handle local shape changes and the availability to obtain sensitivity derivatives. Therefore, apart from circular arc leaf end, the potentials of elliptical arc, Bézier curve, and B-spline are explored and rigorously investigated in our work.

Note that the path length of beam penetration through leaf entity is related to the distal and proximal leaf edge. A piecewise parametric leaf curve is established, which consists of three segments, the leaf end curve, the distal leaf edge, and the proximal leaf edge. Origin of the coordinate system is placed on the leaf average height, and the longitudinal axis y is in accordance with leaf motion. Let lh denote leaf height; leaf end shape parameterization is depicted in Figure 2.

Figure 2

Leaf end shape parameterization: (a) circular arc, (b) elliptical arc, (c) Bézier curve, and (d) B-spline.

Circular arc and elliptical arc are regularized as parametric curves. The uniform representation of leaf curves is as follows:where C is the set of points that satisfy leaf curve parametric equations. Point on leaf curve with design variables p is a function of independent parameter t. For t < 0 and t > 1, parametric curve point is on the distal and proximal leaf edge. Mathematical terms are given as follows:

For t ∈ [0,1], four kinds of leaf end curve parameterization techniques and their design variables are listed in Table 1.

Table 1

Leaf end parameterization and design variables.

Type	Parametric curves	Design variables
Circular arc	x=-Rcos⁡π2-φ+tϕ+φ+dy=Rsin⁡π2-φ+tϕ+φ ϕ=arcsinlh/2-dRφ=arcsinlh/2+dR	p = {R, d}  R: radius of circular arc  d: arc centre offset
Elliptical arc	x=acos⁡(1-t)πy=bsin⁡(1-t)π	p = {a, b}   a: semimajor axis  b: semiminor axis
Bézier curvea	Bt=∑i=0nb_i,ntP_i,  b i,n(t) = C n i t i(1 − t)n−i	p = {P 0, P 1,…, P n}  P i = (x i, y i), i = 0,1,…, n
B-splineb	St=∑i=0nP_iN_i,kt  N_i,0t=1,if  t∈u_i,u_i+10,otherwise  N_i,kt=t-u_iu_i+k-u_iN_i,k-1t+u_i+k+1-tu_i+k+1-u_i+1N_i+1,k-1t  U = [u 0, u 1,…, u n+k+1],  u 0 = 0, u n+k+1 = 1	p = {P 0, P 1,…, P n}  P i = (x i, y i),  i = 0,1,…, n

aBernstein basis polynomial of Bézier curve is denoted by b .

bB-spline is a piecewise polynomial function of degree k, which is defined by n + 1 control points and n + k + 2 knots U. Coefficient N is obtained by recurrence relation.

2.3. Verification and Validation of TST Penumbra Model

Verification and validation of TST penumbra model involved three steps. Firstly, leaf position-penumbra curves grouped by radius of circular arc leaf end are obtained by Monte Carlo simulation. Secondly, the equivalent source size and the effective path length are obtained using parameter identification, with reference penumbra data obtained from Monte Carlo simulation. Thirdly, with Gaussian source approximation, results of ray tracing algorithm are presented for comparison.

Numerical simulation is based on EGSnrc/BEAMnrc Monte Carlo codes [11]. In our study, monoenergetic source of 1.5 MeV is adopted, which is the average energy for 6 MeV source. Source size with Gaussian distribution of 0.2 cm full width at half maximum (FWHM) is used. Leaf is made of W700ICRU, with density ρ of 19.3 g/cm3. X-ray mass attenuation coefficient μ /ρ of tungsten leaf is 0.05 cm2/g for photon energy of 1.5 MeV [12]. Beam angle α about collimator rotation axis is 15.8°, which is determined by field size (FS) and source to axis distance (SAD).

In order to estimate the contribution of leaf end shape to penumbral properties, we develop a specific geometric model for verification and validation of TST model. Parameters of geometric model are listed in Table 2. In order to obtain leaf position-penumbra curve, field size shaped by leaf ends is designated as 10 × 10 cm2, which remains unchanged while the centre of the radiation field shifts.

Table 2

Configurations for TST model verification and validation.

Geometry configurations								Photon source configurations
	SAD	SCD	SDD	lh	dh	FS	ρ	Energy	μ a/ρ	FWHM	α B
Unit	cm	cm	cm	cm	cm	cm	g/cm3	MeV	cm2/g	cm	—
Value	100	46	33.9	8	7.8	40	19.3	1.5	0.05	0.2	15.8°

Virtual energy fluence source models, including single source model, dual source model, three-source model, and hybrid source model, have been proved to be feasible for photon source modelling [13]. By substitution of Gaussian source with virtual source models, various source distributions could be made possible in our program.

2.4. Leaf End Shape Optimization

2.4.1. Penumbra Mean-Variance Optimization Objectives

Idealized collimator system for clinicians should meet the requirements of minimal penumbra width and consistent field variance [8]. To this end, leaf end shape design is formulated as a biobjective problem of penumbra mean-variance optimization, which attempts to minimize penumbra mean for a given penumbra variance or equivalently minimize penumbra variance for a given level of penumbra mean, by carefully choosing leaf end shape. It can be written in the following notation:where J is the biobjective function, J 1 denotes penumbra mean μ, and J 2 denotes standard deviation σ. Constraint functions include bound constraints b , b , linear constraints b , and nonlinear constraints Q(p). For details, refer to Table 3.

Table 3

Leaf end curve constraints.

Type	Bound constraints	Linear constraints
Circular arca	lh2≤R<+∞	-R+d+lh2≤0
	-lh2≤d≤lh2	-R-d+lh2≤0
Elliptical arcb	a=lh2	0<1-b^2a^2^1/2<1
	0≤b≤lh2
Bézier curve	x_0,y_0=-lh2,0,  x_n,y_n=lh2,0	y i + y i+2 − 2y i+1 ≤ 0,  i = 0,1,…, n − 2
	-12+i-1n-1lh≤x_i≤-12+in-1lh  0≤y_i≤lh2, i=1,…,n-1
B-spline	x_i=-12+inlh,  i = 0,1,…, n	y i + y i+2 − 2y i+1 ≤ 0,  i = 0,1,…, n − 2
	0≤y_j≤lh2,  j = 0,1,…, n − 1,y n = 0

  aLeaf end is straight when R tends to infinity.

bThe special case of elliptical arc is circular arc when b equals lh/2.

In our study, leaf position-penumbra curve derivation is composed of two steps, namely, field discretization and sample point evaluation. Let T be leaf position on the scoring plane. Penumbra mean and standard deviation are expressed mathematically as follows:where T is the leaf end point on scoring plane when leaf protrudes fully across the central axis. Point T is denoted by the position where leaf is fully withdrawn. Distance between T and T denotes leaf stroke. Uniform sampling method is adopted to obtain leaf position T . On account of mechanical constraints, leaf stroke can be asymmetric about the collimator rotation axis.

2.4.2. Convex Hull Assumption of Leaf End and Optimization Constraints

It is stated that for a given concave leaf end the corresponding convex hull leaf end can achieve better penumbra performance. A brief proof is presented for illustration. Based on computational geometry theory, convex hull leaf end is designated as the intersection of all sets of convex leaf ends containing the concave leaf end. Note that the relative complement of the convex hull leaf end with respect to the concave leaf end is filled with leaf material for the convex hull leaf end. This observation implies that radiation beams require extra path length for penetration through the convex hull leaf end. Consequently, radiation beams attenuate more quickly and sharper radiation field edge should be achieved. Under the convex hull assumption and geometry boundaries, leaf curve constraints are listed in Table 3. For simplicity, cubic Bézier curve is used with fixed starting point and ending point. Cubic B-spline with 9 control points is utilized with ending point fixed and degree of freedom along y-axis constrained.

2.4.3. Multiobjective Optimization and Pareto Frontier Approximation

Considering that objective space is convex for circular arc and elliptical arc, gradient based iterative algorithm is robust and efficient for shape optimization. Since shape optimization typically involves multiple control variables for Bézier curve and B-spline, it is difficult to distinguish the potential global optimum from local optimum without function convexity information provided in advance. In view of multiple local optima which exist in the objective space, genetic algorithm based global optimization is implemented for approximating Pareto frontier without being trapped near local optima.

A weighted sum approach is introduced by creating a scalar function for mean-variance biobjective shape optimization. A tuning coefficient is introduced to modulate the weight between penumbra mean and standard deviation. In mathematical terms, it can be formulated aswhere J tot denotes the composite objective function. By systematically changing the weight among objectives, the Pareto frontier is obtained. Since function evaluation is most time-consuming, parallel computing technique is used to speed up the process.

2.5. Implementation of Leaf End Shape Optimization

Matlab codes (Mathworks, Natick, MA, USA) have been developed to compute the leaf position-penumbra curve and solve the problem of leaf end shape optimization. Besides the proposed geometric model, the method that we propose is utilized with treatment modalities of Elekta Agility 160-leaf MLC and Varian Millennium 120 MLC. Comparison is conducted with empirical method [8]. The optimal radius equation using empirical method is written aswherewhere optimal radius is a function of μ , SAD, FS, and lh. It is noted that source energy distribution and source to collimator distance are not included in the equation.

3. Results

3.1. Algorithmic Efficiency

As listed in Table 4, efficiency of three algorithms is investigated for the geometric model, including Monte Carlo simulation, ray tracing algorithm, and TST penumbra model. Leaf position-penumbra curve is obtained using 17 sample points, which are evenly spaced from −20 to 20 cm, with interval length of 2.5 cm. In terms of ray tracing algorithm, three-sigma rule is used to determine the range of source energy distribution. Results have shown that the subsource number of 100 results in source energy error of 1%. Method of reduced space searching is implemented for field discretization, which is implemented by searching only the neighbourhood of leaf end projection point on scoring plane for the 80% and 20% intensity points. It takes 100 segments of truncated penumbra region to achieve calculation error less than 1%.

Table 4

Algorithmic efficiency comparison of three algorithms.

Algorithmic efficiency	Monte Carlo	Ray tracing	TST model
Number of leaf positions	17	17	17
Number of histories/line intersectionsa	109	104	2
Computation timeb	33.15 h	8.45 min	0.33 s

aNumber of histories or line intersections is calculated for a given leaf position.

bTiming is recorded on the same workstation without using parallel computing.

3.2. Verification and Validation of TST Penumbra Model

As illustrated in Figure 3, symmetric circular arcs with radius values of 4, 6, 8, 10, 15, 20, and 25 cm are used, which correspond to curves A to G. Error bars indicate the absolute discrepancy between TST model and Monte Carlo simulation. Results of source approximation have shown that penumbra curves of ray tracing algorithm using Gaussian distribution with FWHM of 0.217 cm agree well with Monte Carlo simulation. Parameter identification for TST model reveals that the equivalent source size of TST model is 0.152 cm, and the effective path length is 1.174 cm. Results demonstrate that leaf position-penumbra curves of TST model have a good approximation to Monte Carlo simulation, with maximum value of absolute error 11.9% on curve G. This observation is probably related to the leaf position offset caused by rounded leaf end effect. The observation that curve F and curve G are bowl-shaped is related to beam penetration through the distal and proximal leaf edge.

Figure 3

Verification and validation of TST model.

3.3. Circular Arc Leaf End Shape Optimization

Figure 4 shows that penumbra mean contour profile is convex in objective space. Therefore, gradient based optimization algorithm is robust and efficient to search for the global optimum of penumbra mean, and solution is not sensitive to initial guess. For initial guess at radius of 4 cm without centre offset, it takes 53 iterations and 316 function evaluations to converge to global optimum of penumbra mean at radius of 16.213 cm and centre offset of −0.732 cm. On the other hand, the global optimum of standard deviation is at radius of 4 cm without centre offset, which is the minimal radius for leaf height of 8 cm. Note that marginal regions with saw-tooth-shaped edges in Figure 4 are related to constraints on radius and centre offset. The saw-tooth-shaped edge can be eliminated by decreasing interval length of radius and centre offset. In our study, interval lengths are 1 cm and 0.1 cm for radius and centre offset, respectively.

Figure 4

Penumbra mean and standard deviation of circular arc leaf end. Pareto frontier is derived by weighted sum method, with λ ranging from 0 to 1. The optimum denotes global minimum of penumbra mean.

3.4. Elliptical Arc Leaf End Shape Optimization

Figure 5(a) demonstrates penumbral properties of elliptical arc leaf end. Notably, global optimum of penumbra mean is reached with semiminor axis of 0.728 cm. Global minimum of standard deviation is reached with semiminor axis of 4 cm, which means leaf end is in the shape of circular arc. Pareto frontier of elliptical leaf end is depicted in Figure 5(b). For λ greater than 0.5, penumbra mean is almost invariable, while standard deviation changes rapidly.

Figure 5

Elliptical arc leaf end penumbra mean and standard deviation: (a) semiminor axis dependent penumbral properties, (b) Pareto frontier of elliptical arc leaf end.

3.5. Pareto Frontier of Bézier and B-Spline Curves and A Priori Method for Optimal Point Selection

Points on Pareto frontier are superior in penumbra mean and standard deviation. As illustrated in Figures 6 and 7, Pareto frontiers of Bézier curve and B-spline are obtained using multiobjective genetic algorithm. Local optimum of penumbra mean is obtained by multistart global optimization algorithm. The data suggest that local optima are spread out over the entire objective space. Therefore, gradient based optimization algorithm converges slowly and gets easily trapped in local optimum. Since penumbra width is of interest for clinical application, it is suggested that penumbra mean is preferred over standard deviation for optimal point selection, which is used as a priori knowledge. Consequently, the search for potential of penumbra mean comes first. After that, searching for point on Pareto frontier with relatively small standard deviation is conducted. Note that the Pareto frontier is V-shaped curve; points located at the near-vertex region are selected in our study.

Figure 6

Pareto frontier of cubic Bézier curve. Local optima are obtained by multistart algorithm.

Figure 7

Pareto frontier of B-spline leaf end. Local optima are obtained by multistart algorithm.

3.6. Optimal Leaf End Shapes and Leaf Position-Penumbra Curves

Figure 8 compares the shapes of optimal leaf ends using four parameterization techniques. It is noted that optimal leaf ends of circular arc and B-spline resemble each other closely. Leaf position-penumbra curves are depicted in Figure 9. The smooth curve of Bézier leaf end suggests that constant penumbra width can be achieved while maintaining penumbra mean at an acceptable level. B-spline is advantageous in local shape control; thus penumbra mean and standard deviation can be obtained simultaneously. Notably, leaf position-penumbra curve of circular arc is not monotonically decreasing. The observation is related to the fact of beam penetration through the distal leaf edge.

Figure 8

Shape comparison of optimal leaf ends with four parameterization techniques.

Figure 9

Leaf position-penumbra curves of optimal leaf ends with four parameterization techniques.

Table 5 summarizes all the optimized leaf ends. Results have shown that penumbra mean values are very close, while standard deviation values vary. Minimum of standard deviation is achieved using Bézier curve, and minimum of penumbra mean is obtained using B-spline.

Table 5

Summary of optimal leaf ends and penumbral properties.

Leaf end curve	Design variables (cm)	Penumbra mean (cm)	Standard deviation (cm)
Circular arc	R = 16.213, d = −0.732	0.198	0.0130
Elliptical arc	b = 0.728	0.209	0.0207
Cubic Bézier	P 1 = (−1.594, 1.077), P 2 = (0.814, 0.748)	0.210	0.0005
B-splinea	Y = (0.401, 0.604, 0.699, 0.718, …, 0.671, 0.572, 0.422, 0.240)	0.197	0.0066

a Y denotes the set of y-axis values of control points P 0 to P 7.

3.7. Combined Effect of SCD and Source Size on Optimal Circular Arc Leaf End

As illustrated in Figure 10, combined effect of SCD and source size on the values of the equivalent source size and the effective path length is shown based on the geometric model. Note that a large amount of reference data should be obtained. For the sake of computation time, parameter identification is utilized with the reference data derived from ray tracing algorithm. In Figure 11, results show that optimal circular arc is a function of SCD and source size. Observe that the offset values are uniformly negative. This is mainly related to the effect of geometric penumbra.

Figure 10

Parameter identification of the equivalent source size and the effective path length based on the geometric model with leaf height of 8 cm.

Figure 11

SCD and source size dependent radius and centre offset contour of optimal circular arc leaf ends with leaf height of 8 cm.

3.8. Leaf End Shape Optimization for Elekta Agility 160-Leaf MLC and Varian Millennium 120 MLC

In our study, Gaussian distribution with FWHM of 0.2 cm is assigned to virtual source. Leaf material is typically heavy tungsten alloy, with density of 18 g/cm3. As illustrated in Figure 12, estimation of the equivalent source size and the effective path length is implemented by parameter identification method, with results from ray tracing algorithm for reference. Figure 13 shows the combined effect of SCD and leaf height on optimal radius and centre offset. Gradient based optimization algorithm is utilized to search for optimal circular arc. Results have shown that the optima are located at point A for Elekta Agility 160-leaf MLC and point B for Varian Millennium 120 MLC, respectively.

Figure 12

Parameter identification of the equivalent source size and the effective path length.

Figure 13

Optimal radius and centre offset of circular arc leaf end.

Table 6 summarized the results of optimal leaf ends in the shape of circular arc. Compared with empirical method, TST model can achieve good agreement with actual leaf ends.

Table 6

Results of TST model based optimization and empirical method compared with actual leaf ends.

Type	Actual shape (cm)	Empirical method (cm)	TST model (cm)
Elekta Agility 160-leaf MLCa	R = 17	R = 20.9	R = 16.805, d = −1.153
Varian Millennium 120 MLCb	R = 8, α l = 11.3°	R = 8	R = 12.354, d = −0.125

aSCD = 35.1 cm, leaf height = 9 cm [14].

bSCD = 51.02 cm, leaf height = 5.65 cm [15]. Angle α denotes the angle between line segment of piecewise leaf end curve and collimator rotation axis.

Figure 14 illustrates leaf ends of TST result and empirical result, compared with actual leaf end of Varian Millennium 120 MLC. Notably, empirical method results in optimal circular arc radius of 8 cm, while optimal leaf end using TST model is in close approximation to the actual piecewise leaf end, with maximum lateral deviation of 0.06 cm.

Figure 14

Comparison of optimal circular arc leaf ends using empirical method and TST model with actual shape of Varian Millennium 120 MLC.

4. Discussions

Leaf end shape optimization in our study is composed of four steps, geometric model initialization, penumbra evaluation, parameter identification, and Pareto optimization. The results demonstrate that TST penumbra model based mean-variance optimization approach serves well the purpose of leaf end shape design for multileaf collimator in radiotherapy. The results developed herein may be used for leaf end shape design, as well as analytical penumbra calculation. In this study, we employ a specific geometric model, whose parameters are reasonable for multileaf collimator based radiotherapy. Modification can be introduced to adapt to different treatment modalities. These findings will be testified by experimental measurements, in progress in our research group.

Results of optimal elliptical arc manifest the fact that decreasing the semiminor axis value to be less than global minimum of penumbra mean leads to penumbra mean and standard deviation increase sharply. This is due to the edge effect, which is caused by penetration of radiation beams through the proximal and distal leaf edge. Results of optimal Bézier curve indicate that the goal of consistent penumbra variance across the radiation field can be attained; while using B-spline, minimal penumbra mean is made possible. Due to the flexibility of local shape control, B-spline curve with more design variables would be applied to explore the full potential of leaf end shape induced lateral penumbral properties. Besides, piecewise leaf end composed of line segments and curve segments can be introduced. In addition, the stroke length of leaf travel is designated to be symmetric about the collimator rotation axis. However, asymmetric leaf travel range could be implemented.

In principle, choice of the optimal point on Pareto frontier varies according to preference for penumbra mean and standard deviation. In our study, the selection criterion of optimal point on Pareto frontier is based on a priori rule of penumbra mean first. Nevertheless, adjacent points of optimal point could be used with comparable penumbral characteristics and similar leaf end shape.

Objectives of leaf end shape optimization include, but are not limited to, penumbra mean and standard deviation. For better control of leaf position-penumbra curve, we propose penumbra off axis ratio as one of the objective functions, which is expressed arithmetically as the ratio of peripheral penumbra to axis penumbra. Taking clinical practice into account, penumbra located in the neighbouring region of collimator rotation axis is more frequently used. Therefore, a discrepancy of dosimetric effect exists between axis penumbra and peripheral penumbra. In order to account for leaf position dependent dosimetric effect, variant leaf position sampling or weighted penumbra would be introduced for future studies.

Ray tracing penumbra calculation can yield quantitative and accurate results only when discretization error is sufficiently small. By virtue of the development of computation capacity and algorithm improvement, both ray tracing algorithm and Monte Carlo simulation can be accelerated on GPU. For further investigation, ray tracing based leaf end shape optimization or Monte Carlo simulation based optimization would be made possible.

Although kinetic energy of photon beam is assigned to be monoenergetic, functional form of photon spectra model could be adopted for future study. We would like to point out that in our study leaf end shape optimization is largely dependent on fluence penumbra modelling. In view of the fact that dosimetric penumbra in percent depth dose profile can be estimated by kernel based dose calculation method, such as pencil beam kernel based superposition and convolution algorithm, leaf end shape optimization based on direct dosimetric penumbra modelling will be the subject of future publications.

5. Conclusions

The purpose of this paper is to introduce an analytical method for leaf end shape optimization that can be easily implemented in leaf end design of multileaf collimator. A cost-efficient modelling method of leaf end induced lateral penumbra is proposed based on Tangent Secant Theory, which is verified by Monte Carlo simulation and ray tracing algorithm. Penumbra mean and variance are introduced for biobjective optimization. Leaf end curve parameterization techniques are introduced, including circular arc, elliptical arc, Bézier curve, and B-spline. Observing that objective space is convex for circular arc and elliptical arc, gradient based iterative method is used for local search, while for leaf ends in the shape of Bézier curve and B-spline, objective space is concave and genetic algorithm is applied to explore the full potential of leaf end shape. Results have shown that our method serves well for efficiently deriving optimal radius and centre offset of circular arc. It is found that leaf position-penumbra curve is flat and the goal of consistent penumbra width can be reached using leaf end in the shape of Bézier curve. Results of optimal B-spline leaf end manifest minimum of penumbra mean due to its flexibility of shape representation. Geometries of treatment modalities manufactured by Varian and Elekta are incorporated for TST model based leaf end shape optimization, and it is shown that a relatively close match is found between optimal leaf end and actual shape.

The method that we propose is feasible to estimate leaf position-penumbra curve and suggest optimal leaf end design for a particular treatment modality. Although in this paper geometries of specific treatment modalities are utilized, the conclusions we reach could provide insight into leaf end shape design of multileaf collimator in radiotherapy.