Medical Printed
Donmefry. Vol. 16, pp. 147-I 51 in the U.S.A. All nghts reserved.
Copyright
0
1991 American
0739-021 l/91 Associatron of Medical
$3.00 + .M) Dosimetnsts
GENERATION OF DOSE CALCULATION DATA TABLES USING CUBIC SPLINE INTERPOLATION and SANDRA E. BURCH, M.M.Sc. Georgia Radiation Therapy Center, Department of Radiology, Medical College of Georgia, Augusta, GA 309 12, U.S.A.
KENT H. LARSEN, PH.D.
Abstract--In order to calculate treatment machine settings for a teletherapy machine (e.g., time or monitor units), tables are usually used for variables such as output factor, TMR, percent depth dose. The tables are often generated from data collected at a few points. A linear interpolation is usually used to generate values between the measured points. This can introduce errors as great as 2% between the calculated and actual data points. Using a mathematical software package a computer can generate smooth, accurate curves that agree with measured values to within a few tenths of a percent. This method is not an averaging type of procedure by which a certain function is chosen and parameters are adjusted to force the function to fit the data as closely as possible, but rather is a procedure that fits curves exactly through the measured data points. Key Words: Treatment tables, Computers, Treatment calculations.
bic equation to each successive pair of data points.4 There are four adjustable parameters in each cubic equation. They are calculated in MathCAD using the “cspline” function. They are chosen to make adjacent curves exactly matching and smooth in slope (first derivative) and curvature (second derivative) at the interface. The constants for a curve can, for instance, be determined by requiring that the curve agree with the:
INTRODUCIION The actual dose delivered to a certain point in a phantom or patient from a teletherapy unit depends on several factors such as depth, field size, time, and the presence of beam modifiers (wedges and compensators). There are different methods for calculating the time or number of monitor units to deliver, but all involve the use of tables that are prepared from data collected during installation of a new unit or from an annual calibration. For clinical use, it may be necessary to have data available at half-centimeter intervals over the range of the variables. This results in tables containing hundreds or thousands of data points. As it is impractical to measure each one, selected points over the range of the variables are measured, and the other data points are generated by interpolation. Very often, the method used to generate the tables is a linear interpolation between measured data points. This introduces errors, since the actual function describing the data is not linear. A mathematical software package, MathCAD,’ is available for both IBM PC-type computers and Macintosh computers. This package can generate measured data points by fitting smooth curves instead of straight lines through the measured data points. Data can be entered directly into MathCAD or via spreadsheet programs such as Lotus 1-2-32 and Excel.3 After the final data tables are generated, they can be sent back into the spreadsheet program in order to arrange the data in any desired form for printout.
1. data point at the left end of the interval; 2. data point at the right end of the interval; 3. previous spline in slope at the left end of the interval; and 4. previous spline in curvature at the left end of the interval. The detailed derivation of the equations used to determine the constants, including common assumptions for the equation for the first pair of points, can be found in many numerical methods textbooks (e.g., Reference 4). After the parameters for the curves have been determined, intermediate points are calculated using the “interp” function. Since this method involves fitting a different curve to each pair of data points rather than a single curve to several points, care must be taken to measure each point accurately. In other words, any errors in the individual measurements will not be corrected using this technique. The data can, of course, be smoothed using a technique such as described by Nikesch.’
THEORY RESULTS
MathCAD is a versatile software package that can perform a cubic spline interpolation on a set of data points. The spline method involves fitting a cu-
In order to demonstrate this technique, monitor unit per Gy data from a Varian Clinac 4’j linear accel147
148
Medical Dosimetry
erator was calculated from measured TMR data. The measurements for this table were made at square field sizes of 5, 7, 10, 12, 15, 20, 25, and 30 cm2 and at depths of 1.4,2,5, 10, 15, and 20 cm. Using the cubic spline and linear interpolation methods, data points were calculated at .5 cm intervals for both field size and depth. A comparison of the data for a 12 cm x 12 cm field is shown in Fig. 1. Also shown in Fig. 1 is actual measured data for a 12 cm X 12 cm field size at 1 cm depth increments. It can be seen that the measured data agree much better with the cubic spline interpolated data than with the linear interpolated data. Figure 2 shows a similar plot except that in this case, the 12 cm X 12 cm measured data were not included in the matrix used for the interpolation. Even in this case, the cubic spline data agree very well with the measured data. When the data in Fig. 2 are compared point by point, the maximum difference between the measured and linear interpolated data points is 1.7%, whereas, the maximum difference between the measured and cubic spline interpolated data is 0.4%.
210-
Volume 16. Number 3. I99 I 210.
Measured
0
-
190.
-
Cubic interp -
-.
Linear interp
170
6 ln .Z
5 150 & F E
130
110
90 ~III,I/I,III,lll,lll,lll,lll,.ll,lll,lll,
0 Measured
0
-
190-
-
-
4
6
8
10 12 depth (cm)
14
16
18
20
Fig. 2. Monitor unit data for 12 cm X 12 cm field size showing actual data points and interpolated data at half-centimeter depths using cubic spline and linear interpolation methods. The measured 12 cm X 12 cm data were not included in the matrix used in the interpolation.
Cubic interp -
2
Linear interp
170-
DISCUSSION 6 .=v)
5 150& % E 130-
110-
90
~..,,.,.,,..,...,...,..~,...I..,,~..,~..,
0
2
4
6
8
10 12 depth (cm)
14
16
18
20
Fig. 1. Monitor unit data for 12 cm X 12 cm field size showing actual data points and interpolated data at half-centimeter depths using cubic spline and linear interpolation methods. The measured 12 cm x 12 cm data were included in the matrix used in the interpolation.
The listing of a typical Math’CAD program is given in Fig. 3. The data can be entered directly into the program or read from a file created by the spreadsheet program. In this program, the measured data are read from a file named MUDATA (this listing shows the data used for Fig. 1). This matrix contains field sizes along the top row and depths along the leftmost column. Note that the 0,O element must exist in the matrix even though it is not used, since each line in the file must contain the same number of elements. The program first determines the number of field sizes and depths that were measured. It then puts the field sizes and depths into separate vectors and calculates two more vectors containing the desired field sizes and depths. The generation of the final table is done in two parts. First, the interpolation is performed at each measured field size for the desired depths, then this matrix is used for the interpolation at the desired field sizes. The final table that is written to file MUPERGY can be read directly into EXCEL for the Macintosh
149
Dose calculation tables 0 K. H. LARSENand S. E. BURCH l
l
. * t t
*
’
l
Program
l
l
l
l
l
l
l
“MU”
l
Read in monitor unit table from spreadsheet including depths along the first column and field sizes along the top row (the 0,O element must exist)
l
TAR : : READPRN (AIDDATA) 0 5 7 1.4 101.5 101 2 102.9 102.2 TAR = 5 116.4 113.7 10 149.2 142.3
15 20 t l
l
182.6 237.2
12 99.3 100.5 109.1 131.4
15 98.4 99.3 107.2 126.8
30 97.9 98.7 105.4 122.6
169.9 217.4
163.2 206.8
155.7 194.9
147.6 181.7
25 97 97.7 103.7
30 96.2 96.9 102.6
119.3 141.5 171.9
116.8 137.3 165.1
.
* * Calculate
* *
194.5 253.8
10 100 101.1 110.5 134.8
the number
of rows and columns
l
NR := rows(TAR)
NR = 7
i
:= 0 ..NR
- 2
NC :=
NC = 9
j
:= 0 ..NC
- 2
cola
(TAR)
t * * t t Set up depth and field size vectors * * * l
DEPTH
*
l
.
t
t
t
* t .
i+l,
Set up initial
0
:= TAR O,j+l
j
monitor
unit
matrix
:= TAR
Mu
i+l,j+l
i,j
* * * * Determine * * * l
FLDSZ
:= TAR i
l
NDPTHS :=
number
of depths
floor
DEPTH
in final
matrix
- DEPTH0 + ‘l-2
-
.Ol]
NR-2
[[
floor NDPTHS = 39
m :=
0 . .NDPTHS
DEPTH 0I
[ Dl
:=
Dl
= 0.5
2 Fig. 3. Program “MU”-MathCAD@
program to calculate monitor unit/Gy cubic spline interpolation.
matrix from measured data using
150
Medical Dosimetry
*
l
Volume 16, Number 3, 1991
*
* t t l
Calculate vector of depths for final table
l l
m
CALCDPTH := - + Dl a m t
l
l
* ’ Calculate
monitor units interpolated
[
j #m
l
’
l
l
l
l
4,
,DEPTH,MU
,CALCDPTH m
Determine number of field sizes in final table
FLDSZ
-
FLDSZ
-2+ 1 1
0
NC-2
P := 0 ..NFLDS
NFLDS = 51
l
[
l
[
l
d* 1
l
NPLDS t=
t
for desired depths
:= interp cspline DEPTH,MU
MUD
l
CALCDPTH := DEPTH 1 0
*
*tt
t
CALCDPTH := DEPTH 0 0
l
* Calculate vector of field sizes for final table n
CALCFLD := - + PLDSZ - .5 0 n a .* t Calculate final monitor unit table l
t
l
l
l
l
r
r
xm>l 1 j,PLDSZ,MUD ,CALCFLD := interp csplineLFLDSZ,MUD MUDF I1 m,n J 1
l
l
l
l
’
l
l
l
l
Write table to file
:= CALCDPTH MUDF m m,O
MUDF O,n
:= CALCFLD n
WRITEPRN (MUPERGY) := MUDF Fig. 3 (Contd)
MUDF
:= 0 o,o
Dose calculation tables 0 K. H. LARSENand S. E. BURCH
version. However, LOTUS l-2-3 will not accept records longer than 256 characters. Therefore, the final matrix calculation must be broken up into sections with different field size ranges put into separate matrices. CONCLUSION It is important to make every effort to calibrate radiation therapy treatment machines and calculate treatment settings as accurately as possible. Due to the vast amount of data contained in dose calculation tables for a typical teletherapy unit, it is impractical to measure each data point individually. Automated scanners have alleviated this problem somewhat, but they introduce further uncertainties such as inaccuracies in the position sensors for depth dose/TMR measurements and generation of ripples by a continuously moving chamber.
151
Very accurate TMR measurements can be made in a water phantom by manually adjusting and carefully measuring the water level. These data can be measured at a few data points, and the remainder can be calculated on a digital computer using the cubic spline interpolation method. This can easily be done on an IBM PC-type or Macintosh computer by use of a mathematical software package such as MathCAD@. REFERENCES I. MathCAD, MathSoft, Inc., Cambridge, MA 02 139. 2. LOTUS l-2-3. Lotus Develooment Corn. Cambridae. MA 02142. 3. Microsoft Excel, Microsoft Corp., Redmond, WA 98073. 4. Pizer, S.M. Numerical computing and mathematical analysis. Chicago, IL: Science Research Associates, Inc.; 1975. 5. Nikesch, Walter. Use of spreadsheets to manipulate field scanning data. Medical Dosimetry 15:7 1-75: 1990. 6. Varian. Palo Alto, CA 94303.
Donmefry. Vol. 16, pp. 147-I 51 in the U.S.A. All nghts reserved.
Copyright
0
1991 American
0739-021 l/91 Associatron of Medical
$3.00 + .M) Dosimetnsts
GENERATION OF DOSE CALCULATION DATA TABLES USING CUBIC SPLINE INTERPOLATION and SANDRA E. BURCH, M.M.Sc. Georgia Radiation Therapy Center, Department of Radiology, Medical College of Georgia, Augusta, GA 309 12, U.S.A.
KENT H. LARSEN, PH.D.
Abstract--In order to calculate treatment machine settings for a teletherapy machine (e.g., time or monitor units), tables are usually used for variables such as output factor, TMR, percent depth dose. The tables are often generated from data collected at a few points. A linear interpolation is usually used to generate values between the measured points. This can introduce errors as great as 2% between the calculated and actual data points. Using a mathematical software package a computer can generate smooth, accurate curves that agree with measured values to within a few tenths of a percent. This method is not an averaging type of procedure by which a certain function is chosen and parameters are adjusted to force the function to fit the data as closely as possible, but rather is a procedure that fits curves exactly through the measured data points. Key Words: Treatment tables, Computers, Treatment calculations.
bic equation to each successive pair of data points.4 There are four adjustable parameters in each cubic equation. They are calculated in MathCAD using the “cspline” function. They are chosen to make adjacent curves exactly matching and smooth in slope (first derivative) and curvature (second derivative) at the interface. The constants for a curve can, for instance, be determined by requiring that the curve agree with the:
INTRODUCIION The actual dose delivered to a certain point in a phantom or patient from a teletherapy unit depends on several factors such as depth, field size, time, and the presence of beam modifiers (wedges and compensators). There are different methods for calculating the time or number of monitor units to deliver, but all involve the use of tables that are prepared from data collected during installation of a new unit or from an annual calibration. For clinical use, it may be necessary to have data available at half-centimeter intervals over the range of the variables. This results in tables containing hundreds or thousands of data points. As it is impractical to measure each one, selected points over the range of the variables are measured, and the other data points are generated by interpolation. Very often, the method used to generate the tables is a linear interpolation between measured data points. This introduces errors, since the actual function describing the data is not linear. A mathematical software package, MathCAD,’ is available for both IBM PC-type computers and Macintosh computers. This package can generate measured data points by fitting smooth curves instead of straight lines through the measured data points. Data can be entered directly into MathCAD or via spreadsheet programs such as Lotus 1-2-32 and Excel.3 After the final data tables are generated, they can be sent back into the spreadsheet program in order to arrange the data in any desired form for printout.
1. data point at the left end of the interval; 2. data point at the right end of the interval; 3. previous spline in slope at the left end of the interval; and 4. previous spline in curvature at the left end of the interval. The detailed derivation of the equations used to determine the constants, including common assumptions for the equation for the first pair of points, can be found in many numerical methods textbooks (e.g., Reference 4). After the parameters for the curves have been determined, intermediate points are calculated using the “interp” function. Since this method involves fitting a different curve to each pair of data points rather than a single curve to several points, care must be taken to measure each point accurately. In other words, any errors in the individual measurements will not be corrected using this technique. The data can, of course, be smoothed using a technique such as described by Nikesch.’
THEORY RESULTS
MathCAD is a versatile software package that can perform a cubic spline interpolation on a set of data points. The spline method involves fitting a cu-
In order to demonstrate this technique, monitor unit per Gy data from a Varian Clinac 4’j linear accel147
148
Medical Dosimetry
erator was calculated from measured TMR data. The measurements for this table were made at square field sizes of 5, 7, 10, 12, 15, 20, 25, and 30 cm2 and at depths of 1.4,2,5, 10, 15, and 20 cm. Using the cubic spline and linear interpolation methods, data points were calculated at .5 cm intervals for both field size and depth. A comparison of the data for a 12 cm x 12 cm field is shown in Fig. 1. Also shown in Fig. 1 is actual measured data for a 12 cm X 12 cm field size at 1 cm depth increments. It can be seen that the measured data agree much better with the cubic spline interpolated data than with the linear interpolated data. Figure 2 shows a similar plot except that in this case, the 12 cm X 12 cm measured data were not included in the matrix used for the interpolation. Even in this case, the cubic spline data agree very well with the measured data. When the data in Fig. 2 are compared point by point, the maximum difference between the measured and linear interpolated data points is 1.7%, whereas, the maximum difference between the measured and cubic spline interpolated data is 0.4%.
210-
Volume 16. Number 3. I99 I 210.
Measured
0
-
190.
-
Cubic interp -
-.
Linear interp
170
6 ln .Z
5 150 & F E
130
110
90 ~III,I/I,III,lll,lll,lll,lll,.ll,lll,lll,
0 Measured
0
-
190-
-
-
4
6
8
10 12 depth (cm)
14
16
18
20
Fig. 2. Monitor unit data for 12 cm X 12 cm field size showing actual data points and interpolated data at half-centimeter depths using cubic spline and linear interpolation methods. The measured 12 cm X 12 cm data were not included in the matrix used in the interpolation.
Cubic interp -
2
Linear interp
170-
DISCUSSION 6 .=v)
5 150& % E 130-
110-
90
~..,,.,.,,..,...,...,..~,...I..,,~..,~..,
0
2
4
6
8
10 12 depth (cm)
14
16
18
20
Fig. 1. Monitor unit data for 12 cm X 12 cm field size showing actual data points and interpolated data at half-centimeter depths using cubic spline and linear interpolation methods. The measured 12 cm x 12 cm data were included in the matrix used in the interpolation.
The listing of a typical Math’CAD program is given in Fig. 3. The data can be entered directly into the program or read from a file created by the spreadsheet program. In this program, the measured data are read from a file named MUDATA (this listing shows the data used for Fig. 1). This matrix contains field sizes along the top row and depths along the leftmost column. Note that the 0,O element must exist in the matrix even though it is not used, since each line in the file must contain the same number of elements. The program first determines the number of field sizes and depths that were measured. It then puts the field sizes and depths into separate vectors and calculates two more vectors containing the desired field sizes and depths. The generation of the final table is done in two parts. First, the interpolation is performed at each measured field size for the desired depths, then this matrix is used for the interpolation at the desired field sizes. The final table that is written to file MUPERGY can be read directly into EXCEL for the Macintosh
149
Dose calculation tables 0 K. H. LARSENand S. E. BURCH l
l
. * t t
*
’
l
Program
l
l
l
l
l
l
l
“MU”
l
Read in monitor unit table from spreadsheet including depths along the first column and field sizes along the top row (the 0,O element must exist)
l
TAR : : READPRN (AIDDATA) 0 5 7 1.4 101.5 101 2 102.9 102.2 TAR = 5 116.4 113.7 10 149.2 142.3
15 20 t l
l
182.6 237.2
12 99.3 100.5 109.1 131.4
15 98.4 99.3 107.2 126.8
30 97.9 98.7 105.4 122.6
169.9 217.4
163.2 206.8
155.7 194.9
147.6 181.7
25 97 97.7 103.7
30 96.2 96.9 102.6
119.3 141.5 171.9
116.8 137.3 165.1
.
* * Calculate
* *
194.5 253.8
10 100 101.1 110.5 134.8
the number
of rows and columns
l
NR := rows(TAR)
NR = 7
i
:= 0 ..NR
- 2
NC :=
NC = 9
j
:= 0 ..NC
- 2
cola
(TAR)
t * * t t Set up depth and field size vectors * * * l
DEPTH
*
l
.
t
t
t
* t .
i+l,
Set up initial
0
:= TAR O,j+l
j
monitor
unit
matrix
:= TAR
Mu
i+l,j+l
i,j
* * * * Determine * * * l
FLDSZ
:= TAR i
l
NDPTHS :=
number
of depths
floor
DEPTH
in final
matrix
- DEPTH0 + ‘l-2
-
.Ol]
NR-2
[[
floor NDPTHS = 39
m :=
0 . .NDPTHS
DEPTH 0I
[ Dl
:=
Dl
= 0.5
2 Fig. 3. Program “MU”-MathCAD@
program to calculate monitor unit/Gy cubic spline interpolation.
matrix from measured data using
150
Medical Dosimetry
*
l
Volume 16, Number 3, 1991
*
* t t l
Calculate vector of depths for final table
l l
m
CALCDPTH := - + Dl a m t
l
l
* ’ Calculate
monitor units interpolated
[
j #m
l
’
l
l
l
l
4,
,DEPTH,MU
,CALCDPTH m
Determine number of field sizes in final table
FLDSZ
-
FLDSZ
-2+ 1 1
0
NC-2
P := 0 ..NFLDS
NFLDS = 51
l
[
l
[
l
d* 1
l
NPLDS t=
t
for desired depths
:= interp cspline DEPTH,MU
MUD
l
CALCDPTH := DEPTH 1 0
*
*tt
t
CALCDPTH := DEPTH 0 0
l
* Calculate vector of field sizes for final table n
CALCFLD := - + PLDSZ - .5 0 n a .* t Calculate final monitor unit table l
t
l
l
l
l
r
r
xm>l 1 j,PLDSZ,MUD ,CALCFLD := interp csplineLFLDSZ,MUD MUDF I1 m,n J 1
l
l
l
l
’
l
l
l
l
Write table to file
:= CALCDPTH MUDF m m,O
MUDF O,n
:= CALCFLD n
WRITEPRN (MUPERGY) := MUDF Fig. 3 (Contd)
MUDF
:= 0 o,o
Dose calculation tables 0 K. H. LARSENand S. E. BURCH
version. However, LOTUS l-2-3 will not accept records longer than 256 characters. Therefore, the final matrix calculation must be broken up into sections with different field size ranges put into separate matrices. CONCLUSION It is important to make every effort to calibrate radiation therapy treatment machines and calculate treatment settings as accurately as possible. Due to the vast amount of data contained in dose calculation tables for a typical teletherapy unit, it is impractical to measure each data point individually. Automated scanners have alleviated this problem somewhat, but they introduce further uncertainties such as inaccuracies in the position sensors for depth dose/TMR measurements and generation of ripples by a continuously moving chamber.
151
Very accurate TMR measurements can be made in a water phantom by manually adjusting and carefully measuring the water level. These data can be measured at a few data points, and the remainder can be calculated on a digital computer using the cubic spline interpolation method. This can easily be done on an IBM PC-type or Macintosh computer by use of a mathematical software package such as MathCAD@. REFERENCES I. MathCAD, MathSoft, Inc., Cambridge, MA 02 139. 2. LOTUS l-2-3. Lotus Develooment Corn. Cambridae. MA 02142. 3. Microsoft Excel, Microsoft Corp., Redmond, WA 98073. 4. Pizer, S.M. Numerical computing and mathematical analysis. Chicago, IL: Science Research Associates, Inc.; 1975. 5. Nikesch, Walter. Use of spreadsheets to manipulate field scanning data. Medical Dosimetry 15:7 1-75: 1990. 6. Varian. Palo Alto, CA 94303.
