- Email: [email protected]
An algorithm for planning stereotactic brain implants
Inl. J. Rodiaion Oncology Biol. Phys.. Vol. Printed in the U.S.A. All rights reserved.
13, pp. 631-637 Copyright
0360-3016/87 0 1987 Pergamon
$3X0 + .oO Journals Ltd.
??Computer Applications
AN ALGORITHM
FOR PLANNING
STEREOTACTIC BRAIN IMPLANTS
WENDEL D. RENNER, M.S.,’ NEWELL 0. PUGH, M.D.,2 DAVID B. Ross, RONALD E. BERG, PH.D.~ AND DAVID C. HALL, M.D.3
M.D.,2
‘Department of Radiation Therapy, Community Hospital of Indianapolis, Indianapolis, IN; 2Department of Radiation Therapy, 3Department of Neurosurgery, Methodist Hospital of Indiana, Indianapolis, IN A computer software package was developed for the planning and execution of brain biopsy and radioactive implant procedures with the BRW Stereotaxic System. With the application of computer graphics and a zero-one integer variable programming algorithm, an implant plan with accompanying isodose distributions and stereotactic coordinates can be easily accomplished at the time of the operation when the computer imaging terminal and a printer/plotter is placed in the operating room. Stereotaxic procedure, Computer treatment planning, Interstitual implant, Optimization, Brain neoplasms.
of rapidly planning the implant is necessary because the planning must be done from CT scans with the stereotactic frame in place, and the implant subsequently performed without moving the BRW head ring relative to the patient to ensure accurate placement of the seeds. This does not preclude any pre-planning, but in general a plan must be regenerated in terms of the BRW stereotactic coordinates. The zero-one integer programming algorithm described below is unique in its ability to generate an implant plan from a specification of dose and source stock.
INTRODUCTION
This report describes a mathematical technique and computer software developed for planning and executing the stereotaxic guided implanting of radioactive materials in primary malignant brain tumors. The radioactive implant has the advantage over conventional external beam radiotherapy of delivering a high ratio of dose to the tumor volume relative to surrounding brain tissue. Although an I125implant has a dose gradient of typically 5% per mm at the 50% line as compared to 7-l 5% per mm cited from a competing technology,3 the implant dose gradient usually peaks at much higher doses in the interior of the implant volume, that is, the dose continues to increase to several multiples of the target dose. The clinical significance of this and other radio-biological considerations of the different dose-time fractionations will not be addressed in this report. The use of the computer in our application is rather unique, in that it occupies part of the surgical operative procedure. The software package currently supports the Brown-RobertWells (BRW) Stereotaxic System,4 but can most likely be modified for other systems. The BRW System was originally developed with support by software running on a CT scanner, but because this involved proprietary information concerning the CT scanner, the BRW system is sold with a calculator.4’* The use of the calculator enables the selection of a few specific points for brain biopsy, but was totally inadequate for planning and performing a radioactive brain implant. Software capable
METHODS
AND
MATERIALS
Equipment Equipment consists of: (a) The Brown-Robert-Wells CT Stereotaxic System available from Radionics, Inc., P.O. Box 438, Burlington, MA 01803. (b) A VAX 750 computer from Digital Equipment Corp., Maynard, MA 01754. (c) The Micro-Delta imaging terminal from Computer Design and Applications, Inc. (CDA), 411 Waverly Oaks Rd, Waltham, MA 02154. The MicroDelta is capable of a 256 by 256 pixel display; 256 gray levels or 256 colors out of 16.7 million can be windowed out of 16 bits per pixel. Color plots can be superimposed upon gray level images. (d) An LA210 printer/plotter from Digital Equipment Corp. (e) A GE9800 CT scanner is used to scan patients with the BRW system. The software is written in Fortran 77 and required
Reprint requests to: W.D. Renner, Department of Radiation Therapy, Community Hospital of Indianapolis, 1500 North Ritter Ave., Indianapolis, IN 462 19.
Accepted for publication 4 November 1986. * Hewlett-Packard HP 4 1CV, Palto Alto, CA. 631
632
I. J. Radiation Oncology 0 Biology 0 Physics
about 8000 lines of code in total, not counting routines supplied by Dr. Daniel L. McShan for CDA to read in and translate GE9800 CT tapes. Procedure The use of stereotaxic for radioactive implants has been described in the literature; we are concerned only with the advancements we have made,1,4,5,7 although some introductory description of the software package will be necessary. The patient is scanned with the BRW stereotactic head ring and rod assembly. The patient goes to the operating room (OR) and the CT series goes by way of magnetic tape to the VAX 750 computer system, which supports the Micro-Delta imaging terminal. The Micro-Delta can be placed in the OR with a cable connecting it to the computer. After the tape has been read in, the computer automatically finds the locations of the calibrating rods for each CT scan and calculates the transformation of coordinates matrix between the CT images and the stereotactic unit. It also automatically finds the head contour in each CT scan. A software option is then available for supporting a brain biopsy, with the joystick on the Micro-Delta used for selection of target and entry points. The entry point may also be found by using the stereotactic system as a digitizer, that is, positioning the unit to touch the entry point and entering either the stereotactic or base spatial coordinates at the computer terminal. Needle tracks are plotted on the CT images. Performing a brain implant requires a higher level of sophistication than that required for a biopsy. We have broken the planning procedure down into several steps which will be each described below: (a) Outline tumor volume, (b) define needle direction, (c) specify target dose or dose rate and available source stock, (d) view isodose curves after computer algorithm plans the implant, (e) needle and sources may be alternately specified individually or edited, (f) print out stereotactic coordinates and needle loadings for the surgeon. The joystick is used to outline the tumor in each successive co-planar CT scan. These outlines are used to define the treatment volume and the implant is designed so that all areas within the treatment volume receive at least the prescribed dose, as will be explained in more detail below. The algorithm for constructing a volume model from successive contours differs from previous reports2,5 by which points between successive CT planes are treated. In each CT plane, points inside the tumor contour are inside the tumor volume model, and points outside the contour are outside the tumor volume. For points between two successive co-planar CT planes with tumor contours, the following considerations are made. A particular point’s position relative to a projection on a common co-planar plane of the two successive tumor contours is determined (see Fig. 1a). Any point between two planes whose projection falls inside both contours is considered inside the tumor volume, as is point A in Fig-
April 1987, Volume 13, Number 4
Y
*X
Fig. 1a. Two hypothetical co-planar tumor contours projected Into a common co-planar plane. Point A is inside both con.ours, hence inside the tumor volume. Point B is outside the tumor volume. Points C and E, which lie on top of each other, are outside one contour and inside the other. D is a point on the contour C and E are outside of, giving the shortest distance from the contour to C and E.
ure 1a. Any point outside both contours mor volume, as is point B, Figure 1a.
is not in the tu-
The only difficulty occurs when a point is inside one contour and outside the other, as points C and E in Figure la. To determine whether such a point is inside or
Y
‘2
Fig. 1b. View of contours edgewise. Point C is seen to be outside the tumor volume, whereas point E is considered inside the tumor volume.
Stereotacticbrainimplants0 W. D. RENNERet al.
Fig. lc. Three dimensional view of the two co-planar contours. outside the tumor volume, one finds the shortest path along the common projected plane from the point in question to the contour it is outside of. Then this line is considered in three-space (lines DC or DE, Figs. 1b and lc), connecting the point in question and the closest point on the inside contour, and running through the plane of the other contour. If the line intersects the other plane inside the contour (line DE, Figs. lb and c), then the point in question is considered inside the volume. If the intersection point is outside the contour (line DC, Figs. 1b and c), the point in question is outside the tumor volume. Shown in Figure 2 are two contours with three ;“t‘ XV..,Pm;m” -P&n,4 The IULbI “~UlLl~ ~fi~+n..ro ~“IIC”U‘J Fn..n,4 I”UIIU &th **I&Uth& LUlDIIItiCII”U. method provides for a generalized technique for constructing volume models from successive contours and for interpolating intervening contours given any two contours in separate co-planar planes.
Fig. 2. Three intervening contours, clashed lines, in evenly spaced planes, interpolated in between two co-planar contours, solid lines, using the volume model algorithm.
633
As all implant needles are constrained to be parallel in the algorithm to be described below, it is first necessary to define a vector pointing in the desired direction. This is most likely determined by using the BRW stereotactic unit as a digitizer to determine the coordinates of the entry point. One enters either the stereotactic coordinates (alpha, beta, gamma, delta and depth), or the spatial coordinates (lateral, A-P, vertical) found by transferring the stereotactic assembly to the base frame.4 A target point mav, then iovstick. ------ be located on the CT image with the a-,-~ Alternately, a three-dimensional wire frame plot of tumor and head contours may be viewed until an acceptable angle of approach is found, or the entry point may be found on a CT image with the joystick. Lastly, the user now specifies the desired dose rate for temporary implants or the total dose for permanent implants. Temporary implants would include the use of high activity 1125,Ir 192, or other suitable nuclides. Au 198 or low activity I125might be used for a permanent implant. The software supports the use of any radioactive nuclide whose specifics are coded in a generalized source file. Angular anisotropy, which is pronounced for 1’25,is accounted for in the calculations, with the needle direction defining the orientation of the source. After specifying the target dose (or dose rate as the case may be), the program asks the user to specify the available source stock. The program then proceeds to locate needles and sources using the zero-one integer variable program formulism described in detail below. Alternately, the user may specify the needles and restrict the _I__-Zr~- L^ _l__Z-- suurws _^___^__on _- LL_ _-__1~_> -___3*__ mt: sptxineu ntxxws. iugomnm LOplacing Zero-one integer programming The use of zero-one integer programming for planning radioactive implants was described in detail in an earlier report.6 Our goal in the use of this formulism is two-fold. First, we want to put the least number of seeds in the tumor volume that will irradiate the entire tumor volume to at least the target dose. By minimizing the number of seeds inside the target volume, we minimize the dose to the brain outside of the target volume. Second, we want a computer algorithm that can accomplish this task quickly for us. Two difficulties were encountered in the prior report: solution times became excessive for anything other than a small volume and planned implants could not be realistically carried out by the surgeon. The latter difficulty has been removed by the BRW Stereotaxic System, which allows for the generation of stereotactic coordinates for guided placement of a needle with an accuracy of 1 to 2 mm. The formal difficulty is removed here by a rethinking of the application of the algorithm. The algorithm proceeds by first generating a mathematical lattice of closely spaced possible seed locations. The lattice is rotated so that one axis correlates to the implant needle direction. By using a finely spaced lattice, which is scaled to the tumor volume size, 25 to 125 possi-
634
I. J. Radiation Oncology 0 Biology 0 Physics
ble seed locations are typically considered. Only those points which fall within the tumor volume are viable points for consideration. The seed lattice is modified in that along each needle trajectory the seed locations were shifted so that one seed location will lie midway between the entry and exit points of the tumor volume, ensuring the possibility of centrally located seeds. A second lattice array of “constraint” points is generated, the spacing of which was chosen to be one-half that used for the seed lattice, resulting in roughly eight times more constraint points than seed locations. Again only those points within the tumor volume are considered. Additional points on the tumor volume surface were added to ensure that a small corner of the tumor volume would not be missed by a constraint point. Each seed location is numbered, designated here by the variable j, and is represented by a zero-one integer variable Xj, that is, a variable which may have only the value of zero or one. For a total of N possible seed locations, the dose, dosei, at each of M constraint points may be represented by the linear equation: $ a$jXj = Dosei for all i = 1 to M
(1)
j=l
where aij is the contribution of source location j to constraint location i accounting for attenuation, distance, angular dependency, and gamma factor. Sj is the source strength assigned to variable Xj from the available source stock as will be described below, and is not absorbed into &j to allow for differing source strengths. Because only constraint points within the tumor volume are considered, we require each constraint point to have a dose equal to or greater than the target dose (TD) or dose rate, as the case may be: s a$jXj 2 TD for all i = 1 to M
(2)
j=l
A zero-one integer variable programming problem would require the assignment of a cost to each variable Xj and the generation of an objective function whose value is to be minimized while satisfying all constraints. This was accomplished previously by letting the objective function be the sum ofthe number of seeds.6 By minimizing the number of seeds while maintaining a floor on the dose to the tumor volume, the tumor volume dose is achieved with the minimum number of seeds whose placement is consequently somewhat optimal, and the dose to the volume surrounding the tumor is minimized since the total number of seeds is minimized.6 We have noted that finding the mathematical solution to such a formulated problem posed difficulties,‘j and so we have pursued here an alternate solution method adapted from the earlier work.6,8 The available source stock is assumed entered by the user in the order in which use is preferred, that is, the
April 1987, Volume 13, Number 4
algorithm will select sources in the order entered. For lack of any other criterion, the hottest sources might be specified first, as they will generally end up in the center of the tumor volume. The algorithm will always fix a variable S at the next available source strength. Initially all seed location variables Xj are set to zero. We will define infeasibility B by the sum!
B = z Min [O.,(g %jSjXj] - TD] i=l
(3)
j=l
where we are summing up, for non-satisfied constraints, the constraint dose less the target dose. Initially this sum is just -TD M. When all constraints are satisfied, B will equal zero. Now an iteration procedure begins that continues until either B = 0 or all seed locations are filled. Consider the set of all unfilled seed locations (initially all seed locations j = 1 to N), Xk. For each open seed location k, temporarily set Xk = 1 and Sk = S. Compute the infeasibility, now Bk, given by Eq. (3) above. Then set Xk back to zero, increment k and repeat. Pick k so that B1,is maximum and set Xk = 1 and Sk = S permanently (thus removing from the set Xk of all unfilled seed locations). In effect, find the seed location which does the most to satisfy all constraints. Initially, this will generally be a centrally located seed. Note that in Eq. (3) satisfied constraints are no longer considered as they no longer contribute to the sum of infeasibility. For example, after picking the first centrally located seed, there will now be a hole in the volume of unsatisfied constraint points, and the next iteration has to consider those unsatisfied constraint points that are left. Set S to the next available source strength from the source stock and repeat the above. When Bk = 0, or all source locations are picked, stop. If the algorithm stops because all source locations are picked, then some constraint points are not satisfied and an appropriate error message must be written out. Hotter sources may be required or a more closely spaced seed lattice must be considered. Otherwise a plan is generated which guarantees the target dose to all of the tumor volume. The solution time typically requires about 30 seconds on the VAX 750. The prescribed dose covers the tumor volume and can be easily verified by reviewing the dose distribution in each CT scan. Figures 3-6 are shown for a computer generated plan assuming 5.0 mCi I 125 sources and a target dose rate of 4 1.7 rad/hr. Figure 7 shows a plot in one CT plane for another patient case planned with the zero-one integer programming algorithm and 40 mCi 1125 sources. The isodose lines are shown with tick marks on the low dose side of the isodose line. The tumor contour has no tick marks. The projection to this plane of the needles and sources are also shown along with the head contour. Note that the target dose rate of 4 1.7 rad/hour does cover the tumor area.
Stereotactic brain implants 0 W. D. RENNER et al.
635
Fig. 5. Wire frame plot (again shown in black and white but otherwise in color) of tumor contours, implant needles, and Fig. 3. Wire frame plot (shown here in black and white but otherwise in color) of tumor contours and needles looking
sources with a view line perpendicular to the implant needles. The view line may be rotated about the needles. Similar plots can be made for any arbitrary view line.
down the implant needles. Plus signs mark existing needles. Similar plots can be made for any arbitrary view line. Although this algorithm does not necessarily find the mathematical minimum number of seeds and their corresponding optimal distribution, all constraints are strictly enforced and the solution found above will, in general, be close enough to an optimal solution for clinical purposes. This can be seen because the algorithm proceeds in a fashion similar to current practice. In general, the surgeon performing a seed implant attempts to distribute seeds throughout the tumor volume in a uniform
Fig. 4. Plot of tumor in intersecting plane (shown here in black and white but otherwise in color) perpendicular to the implant needles. Tumor is solid area (colored red), round circles indicate where tumor contours intersect this plane. Plus signs mark implant needles. A 1 cm grid is superimposed. Similar planes can be shown for any angle and depth.
distribution fashion until the combined contribution of dose from all seeds irradiates the tumor volume to at least the prescribed dose. The advantage of the above zero-one integer algorithm is that the infeasibility B provides a rational way to successively pick the seed locations, the entire task can be performed by the computer, the target dose is strictly enforced throughout the tumor volume, and the integrated dose outside the tumor volume is minimized. Note, however, that nothing is done to minimize the number of needles as sources which lay on the same lat-
Fig. 6. CT scan of patient case, with the dose area of 4 1.7 rad/ hour or more blotted out (shown here in solid white). CT image under the blotted area may be alternately viewed by pushing a key. A complete set of isodose curves may be plotted on the LA2 10 printer/plotter along with the head and tumor contours.
636
I. J. Radiation Oncology 0 Biology 0 Physics
Fig. 7. Isodose plot printed on the LA2 10 printer/plotter for a rather large tumor planned with the algorithm. Shown are 10, 20, 41.7, and 100 rad/hour isodose lines. The target dose rate was 41.7 rad/hour (1000 rads/day) and seven 40 mCi I 125
sources were found in the solution.
tice line in the needle direction are simply reported to be on the same needle. Ideally, the number of needles should be minimized in a computer solution of an implant plan, since each needle track can result in brain damage. This could be accomplished by performing a mathematical OR operation, as in Boolian algebra, among all seed variables along each needle, and adding up the result for all needles for a total cost to be minimized. The resulting objective function would not be a linear function, however, and that approach was not further explored. If the needle tracks found by the algorithm are not acceptable, the computer program will allow restricting possible seed location to specific needle tracks either specified by the user as described below, or kept from a first pass solution. The algorithm can then be rerun for just the specified needle tracks, with seed locations only picked from those on the available needle tracks. In this manner, the computer algorithm allows for an opportunity to reduce the number of needles.
Other software supportfeatures As indicated above, the number of needle tracks is not minimized in the computer solution, and so means must be provided for the user to edit needles and sources and to view resultant isodose distributions. With the joystick, the user may locate needles on any of three types of images: CT image, wire frame plot, or intersecting plane plot. A color wire frame plot of the tumor contours as seen looking down the needles can be produced on the imaging terminal (shown in Fig. 3 in black and white only). This view allows for locating a needle manually so that it does intersect the tumor volume, and the distribu-
April 1987, Volume 13, Number 4
tion of the needle tracks in a plane perpendicular to the needles is readily appreciated and may be modified. Needles may be located with the joystick and the target point (end point of the needle) is initially located by the computer program on the exit surface of the tumor volume. The needles may also be located with the joystick on a plot of a plane intersecting with the tumor volume (Fig. 4). This is a plane perpendicular to the needles running through the tumor volume at selected depths, showing the shape of the tumor volume model in that plane. The use of this type of graphic representation is described in more detail by Kelly, et a1.5 In our implementation here, Fig. 4, we have constructed our tumor volume in terms of square cubes, voxels, scaled to the size of the CT image pixel. Hence each of three sides of the voxel is 1.5 mm in length for this case. Planes intersecting these cubes will leave triangular patterns along the edge of the tumor volume. As described below, one can view the coverage of the tumor area shown by overlaying a representation of computed dose. The user may also specify sources on particular needles, with a wire frame plot of a view perpendicular to the needles (Fig. 5). The user indicates the distance along the needles from the target end for the source placement. A negative distance will extend the needles deeper. At printout time, all needle target points are placed one centimeter beyond the end source. The printout of the stereotactic coordinates of all needles and the loadings along each needle is accomplished from a single command from the imaging terminal. A printer/plotter must be present for presenting this information to the surgeon. For any CT scan or intersecting plane, the dose distribution may be viewed for the planned implant (Fig. 6). A solid area for which a selected dose is equaled or exceeded is colored in on the overlay plane of the imaging terminal, and may be alternately switched off by pushing a button on the imaging terminal leaving the underlying CT image (or intersecting plane plot). Different doses are viewed separately by sequentially specifying the dose or dose rate. This method of displaying was chosen due to the limited graphics resolution capabilities of the imaging terminal, graphic pixels being the same size as the image pixels, and the need to easily appreciate the region covered by the dose. Traditional looking plots of line isodose contours may be produced rapidly on the printer/ plotter from a single command at the imaging terminal, and includes the tumor and head contours found from the CT scan (these plots are only available for the CT scans). Fig. 7 shows such a plot.
Accuracy of the BR Wsystem The ability of the BRW system to accurately locate a point in space found on a CT image was tested by scanning a plastic phantom with small embedded air holes about 1 mm in diameter. Coordinates in the stereotaxic system were computed by one of three methods: (a) The
631
Stereotactic brain implants 0 W. D. RENNER et al.
HP calculator and program supplied with the BRW system was used with numbers generated with the track ball on a GE9800 CT scanner. This included calibrating rod coordinates and target points. The image was 5 12 by 5 12 pixels, one pixel 0.8 mm in size with a slice thickness of 1.5 mm. (b) These same numbers were inputted as data into the program described in this report. (c) The joystick on the Micro-Delta terminal was used to locate the same points on the image. The image on the Micro-Delta had been reduced to 256 by 256 pixels, pixel size now 1.6 mm. A needle was then adjusted to the center of each hole and the spatial coordinates measured by transferring to the BRW base assembly, and then compared to that calculated above. The results for three target points are shown in Table 1. From the data, our program did better than the HP calculator with the same numbers. This is most likely because we solved for the transformation of coordinates matrix between the stereotactic system and image system with least squares using all of the available data, whereas the HP calculator uses a geometric solution. The Micro-Delta, however, did worst, most likely because the image pixel size was increased to 1.6 mm in reducing the 512 by 512 image down to a 256 by 256 image, thereby increasing the uncertainty introduced by the pixel size. Hence pixel size will have an effect on the
Table 1. Distance in mm between target point computed from the CT image and that measured with the BRW Stereotaxic System HP Target point 1 2 3
calculator 1.0 1.3 1.4
VAX with same numbers
Micro-Delta with joystick
0.7 0.5 1.2
1.5 2.1 2.0
accuracy. This was further demonstrated when a sliced thickness of 1.Ocm was tried. Errors on the order of 1 to 2 mm were found in the plane of the image, but increased in the direction perpendicular to the image to as large as 6 mm. A small slice thickness is essential for good accuracy in all three dimensions. placement
CONCLUSION The software algorithm described above for supporting stereotactic brain biopsy and implants allows for the rapid planning of the implant after the patient is scanned with the BRW system. It is possible to plan the implant as the patient is prepared for the surgical procedure, with the computer generation of needle track coordinates for the BRW system, along with isodose curves in selected CT scan planes.
REFERENCES Bouzaglou, A., Dyck, P., Solt-Bohman, L.G., Gruskin, P.: Stereotactic interstitial implantation of brain tumors. Endocurie. Hypertherm. Oncol. 1: 99-l 12, 1985. Christiansen, H.N., Sederberg, T.W.: Conversion of complex contour line definitions into polygonal element mosaits. Comp. Graph. 12: 187-192, 1978. Hartmann, G., Schlegel, W., Sturm, V., Kober, B., Pastyr, O., Lorenz, W.: Cerebral radiation surgery using moving field irradiation at a linear accelerator facility. Int. J Radiat. Oncol. Biol. Phys. 11: 1185-l 192, 1985. Heilbrun, M.P., Roberts, T.S., Apuzzo, M.L.J., Wells, T.H., Sabshin, J.K.: Preliminary experience with BrownRobert-Wells (BRW) computerized tomography stereotaxic guidance system. J. Neurosurg. 59: 2 17-222, 1983.
Kelly, P.J., Kall, B. A., Goerss, S.: Computer simulation for the stereotactic placement of interstitial radionuclide sources into computed tomography defined tumor volumes. Neurosurgery 14: 442-448, 1984. Renner, W.D., O’Connor, T.P., Paulauskas, N.M.: Computer assistance in planning radiotherapy seed implants. Int. J. Radiat. Oncol. Biol. Phys. 5: 427-432, 1979. Rossman, K.J., Shetter, A.G., Speiser, B.L., Nehls, D.: Stereotactic afterloading Iridium implants in treatment of high-grade astrocytomas. Endocurie Hypertherm. Oncol. 1: 49-57, 1985. Zionts, S.: Linear and Integer Programming, Englewood Cliffs, NJ, Prentice-Hall. 1974, pp. 320-475.
13, pp. 631-637 Copyright
0360-3016/87 0 1987 Pergamon
$3X0 + .oO Journals Ltd.
??Computer Applications
AN ALGORITHM
FOR PLANNING
STEREOTACTIC BRAIN IMPLANTS
WENDEL D. RENNER, M.S.,’ NEWELL 0. PUGH, M.D.,2 DAVID B. Ross, RONALD E. BERG, PH.D.~ AND DAVID C. HALL, M.D.3
M.D.,2
‘Department of Radiation Therapy, Community Hospital of Indianapolis, Indianapolis, IN; 2Department of Radiation Therapy, 3Department of Neurosurgery, Methodist Hospital of Indiana, Indianapolis, IN A computer software package was developed for the planning and execution of brain biopsy and radioactive implant procedures with the BRW Stereotaxic System. With the application of computer graphics and a zero-one integer variable programming algorithm, an implant plan with accompanying isodose distributions and stereotactic coordinates can be easily accomplished at the time of the operation when the computer imaging terminal and a printer/plotter is placed in the operating room. Stereotaxic procedure, Computer treatment planning, Interstitual implant, Optimization, Brain neoplasms.
of rapidly planning the implant is necessary because the planning must be done from CT scans with the stereotactic frame in place, and the implant subsequently performed without moving the BRW head ring relative to the patient to ensure accurate placement of the seeds. This does not preclude any pre-planning, but in general a plan must be regenerated in terms of the BRW stereotactic coordinates. The zero-one integer programming algorithm described below is unique in its ability to generate an implant plan from a specification of dose and source stock.
INTRODUCTION
This report describes a mathematical technique and computer software developed for planning and executing the stereotaxic guided implanting of radioactive materials in primary malignant brain tumors. The radioactive implant has the advantage over conventional external beam radiotherapy of delivering a high ratio of dose to the tumor volume relative to surrounding brain tissue. Although an I125implant has a dose gradient of typically 5% per mm at the 50% line as compared to 7-l 5% per mm cited from a competing technology,3 the implant dose gradient usually peaks at much higher doses in the interior of the implant volume, that is, the dose continues to increase to several multiples of the target dose. The clinical significance of this and other radio-biological considerations of the different dose-time fractionations will not be addressed in this report. The use of the computer in our application is rather unique, in that it occupies part of the surgical operative procedure. The software package currently supports the Brown-RobertWells (BRW) Stereotaxic System,4 but can most likely be modified for other systems. The BRW System was originally developed with support by software running on a CT scanner, but because this involved proprietary information concerning the CT scanner, the BRW system is sold with a calculator.4’* The use of the calculator enables the selection of a few specific points for brain biopsy, but was totally inadequate for planning and performing a radioactive brain implant. Software capable
METHODS
AND
MATERIALS
Equipment Equipment consists of: (a) The Brown-Robert-Wells CT Stereotaxic System available from Radionics, Inc., P.O. Box 438, Burlington, MA 01803. (b) A VAX 750 computer from Digital Equipment Corp., Maynard, MA 01754. (c) The Micro-Delta imaging terminal from Computer Design and Applications, Inc. (CDA), 411 Waverly Oaks Rd, Waltham, MA 02154. The MicroDelta is capable of a 256 by 256 pixel display; 256 gray levels or 256 colors out of 16.7 million can be windowed out of 16 bits per pixel. Color plots can be superimposed upon gray level images. (d) An LA210 printer/plotter from Digital Equipment Corp. (e) A GE9800 CT scanner is used to scan patients with the BRW system. The software is written in Fortran 77 and required
Reprint requests to: W.D. Renner, Department of Radiation Therapy, Community Hospital of Indianapolis, 1500 North Ritter Ave., Indianapolis, IN 462 19.
Accepted for publication 4 November 1986. * Hewlett-Packard HP 4 1CV, Palto Alto, CA. 631
632
I. J. Radiation Oncology 0 Biology 0 Physics
about 8000 lines of code in total, not counting routines supplied by Dr. Daniel L. McShan for CDA to read in and translate GE9800 CT tapes. Procedure The use of stereotaxic for radioactive implants has been described in the literature; we are concerned only with the advancements we have made,1,4,5,7 although some introductory description of the software package will be necessary. The patient is scanned with the BRW stereotactic head ring and rod assembly. The patient goes to the operating room (OR) and the CT series goes by way of magnetic tape to the VAX 750 computer system, which supports the Micro-Delta imaging terminal. The Micro-Delta can be placed in the OR with a cable connecting it to the computer. After the tape has been read in, the computer automatically finds the locations of the calibrating rods for each CT scan and calculates the transformation of coordinates matrix between the CT images and the stereotactic unit. It also automatically finds the head contour in each CT scan. A software option is then available for supporting a brain biopsy, with the joystick on the Micro-Delta used for selection of target and entry points. The entry point may also be found by using the stereotactic system as a digitizer, that is, positioning the unit to touch the entry point and entering either the stereotactic or base spatial coordinates at the computer terminal. Needle tracks are plotted on the CT images. Performing a brain implant requires a higher level of sophistication than that required for a biopsy. We have broken the planning procedure down into several steps which will be each described below: (a) Outline tumor volume, (b) define needle direction, (c) specify target dose or dose rate and available source stock, (d) view isodose curves after computer algorithm plans the implant, (e) needle and sources may be alternately specified individually or edited, (f) print out stereotactic coordinates and needle loadings for the surgeon. The joystick is used to outline the tumor in each successive co-planar CT scan. These outlines are used to define the treatment volume and the implant is designed so that all areas within the treatment volume receive at least the prescribed dose, as will be explained in more detail below. The algorithm for constructing a volume model from successive contours differs from previous reports2,5 by which points between successive CT planes are treated. In each CT plane, points inside the tumor contour are inside the tumor volume model, and points outside the contour are outside the tumor volume. For points between two successive co-planar CT planes with tumor contours, the following considerations are made. A particular point’s position relative to a projection on a common co-planar plane of the two successive tumor contours is determined (see Fig. 1a). Any point between two planes whose projection falls inside both contours is considered inside the tumor volume, as is point A in Fig-
April 1987, Volume 13, Number 4
Y
*X
Fig. 1a. Two hypothetical co-planar tumor contours projected Into a common co-planar plane. Point A is inside both con.ours, hence inside the tumor volume. Point B is outside the tumor volume. Points C and E, which lie on top of each other, are outside one contour and inside the other. D is a point on the contour C and E are outside of, giving the shortest distance from the contour to C and E.
ure 1a. Any point outside both contours mor volume, as is point B, Figure 1a.
is not in the tu-
The only difficulty occurs when a point is inside one contour and outside the other, as points C and E in Figure la. To determine whether such a point is inside or
Y
‘2
Fig. 1b. View of contours edgewise. Point C is seen to be outside the tumor volume, whereas point E is considered inside the tumor volume.
Stereotacticbrainimplants0 W. D. RENNERet al.
Fig. lc. Three dimensional view of the two co-planar contours. outside the tumor volume, one finds the shortest path along the common projected plane from the point in question to the contour it is outside of. Then this line is considered in three-space (lines DC or DE, Figs. 1b and lc), connecting the point in question and the closest point on the inside contour, and running through the plane of the other contour. If the line intersects the other plane inside the contour (line DE, Figs. lb and c), then the point in question is considered inside the volume. If the intersection point is outside the contour (line DC, Figs. 1b and c), the point in question is outside the tumor volume. Shown in Figure 2 are two contours with three ;“t‘ XV..,Pm;m” -P&n,4 The IULbI “~UlLl~ ~fi~+n..ro ~“IIC”U‘J Fn..n,4 I”UIIU &th **I&Uth& LUlDIIItiCII”U. method provides for a generalized technique for constructing volume models from successive contours and for interpolating intervening contours given any two contours in separate co-planar planes.
Fig. 2. Three intervening contours, clashed lines, in evenly spaced planes, interpolated in between two co-planar contours, solid lines, using the volume model algorithm.
633
As all implant needles are constrained to be parallel in the algorithm to be described below, it is first necessary to define a vector pointing in the desired direction. This is most likely determined by using the BRW stereotactic unit as a digitizer to determine the coordinates of the entry point. One enters either the stereotactic coordinates (alpha, beta, gamma, delta and depth), or the spatial coordinates (lateral, A-P, vertical) found by transferring the stereotactic assembly to the base frame.4 A target point mav, then iovstick. ------ be located on the CT image with the a-,-~ Alternately, a three-dimensional wire frame plot of tumor and head contours may be viewed until an acceptable angle of approach is found, or the entry point may be found on a CT image with the joystick. Lastly, the user now specifies the desired dose rate for temporary implants or the total dose for permanent implants. Temporary implants would include the use of high activity 1125,Ir 192, or other suitable nuclides. Au 198 or low activity I125might be used for a permanent implant. The software supports the use of any radioactive nuclide whose specifics are coded in a generalized source file. Angular anisotropy, which is pronounced for 1’25,is accounted for in the calculations, with the needle direction defining the orientation of the source. After specifying the target dose (or dose rate as the case may be), the program asks the user to specify the available source stock. The program then proceeds to locate needles and sources using the zero-one integer variable program formulism described in detail below. Alternately, the user may specify the needles and restrict the _I__-Zr~- L^ _l__Z-- suurws _^___^__on _- LL_ _-__1~_> -___3*__ mt: sptxineu ntxxws. iugomnm LOplacing Zero-one integer programming The use of zero-one integer programming for planning radioactive implants was described in detail in an earlier report.6 Our goal in the use of this formulism is two-fold. First, we want to put the least number of seeds in the tumor volume that will irradiate the entire tumor volume to at least the target dose. By minimizing the number of seeds inside the target volume, we minimize the dose to the brain outside of the target volume. Second, we want a computer algorithm that can accomplish this task quickly for us. Two difficulties were encountered in the prior report: solution times became excessive for anything other than a small volume and planned implants could not be realistically carried out by the surgeon. The latter difficulty has been removed by the BRW Stereotaxic System, which allows for the generation of stereotactic coordinates for guided placement of a needle with an accuracy of 1 to 2 mm. The formal difficulty is removed here by a rethinking of the application of the algorithm. The algorithm proceeds by first generating a mathematical lattice of closely spaced possible seed locations. The lattice is rotated so that one axis correlates to the implant needle direction. By using a finely spaced lattice, which is scaled to the tumor volume size, 25 to 125 possi-
634
I. J. Radiation Oncology 0 Biology 0 Physics
ble seed locations are typically considered. Only those points which fall within the tumor volume are viable points for consideration. The seed lattice is modified in that along each needle trajectory the seed locations were shifted so that one seed location will lie midway between the entry and exit points of the tumor volume, ensuring the possibility of centrally located seeds. A second lattice array of “constraint” points is generated, the spacing of which was chosen to be one-half that used for the seed lattice, resulting in roughly eight times more constraint points than seed locations. Again only those points within the tumor volume are considered. Additional points on the tumor volume surface were added to ensure that a small corner of the tumor volume would not be missed by a constraint point. Each seed location is numbered, designated here by the variable j, and is represented by a zero-one integer variable Xj, that is, a variable which may have only the value of zero or one. For a total of N possible seed locations, the dose, dosei, at each of M constraint points may be represented by the linear equation: $ a$jXj = Dosei for all i = 1 to M
(1)
j=l
where aij is the contribution of source location j to constraint location i accounting for attenuation, distance, angular dependency, and gamma factor. Sj is the source strength assigned to variable Xj from the available source stock as will be described below, and is not absorbed into &j to allow for differing source strengths. Because only constraint points within the tumor volume are considered, we require each constraint point to have a dose equal to or greater than the target dose (TD) or dose rate, as the case may be: s a$jXj 2 TD for all i = 1 to M
(2)
j=l
A zero-one integer variable programming problem would require the assignment of a cost to each variable Xj and the generation of an objective function whose value is to be minimized while satisfying all constraints. This was accomplished previously by letting the objective function be the sum ofthe number of seeds.6 By minimizing the number of seeds while maintaining a floor on the dose to the tumor volume, the tumor volume dose is achieved with the minimum number of seeds whose placement is consequently somewhat optimal, and the dose to the volume surrounding the tumor is minimized since the total number of seeds is minimized.6 We have noted that finding the mathematical solution to such a formulated problem posed difficulties,‘j and so we have pursued here an alternate solution method adapted from the earlier work.6,8 The available source stock is assumed entered by the user in the order in which use is preferred, that is, the
April 1987, Volume 13, Number 4
algorithm will select sources in the order entered. For lack of any other criterion, the hottest sources might be specified first, as they will generally end up in the center of the tumor volume. The algorithm will always fix a variable S at the next available source strength. Initially all seed location variables Xj are set to zero. We will define infeasibility B by the sum!
B = z Min [O.,(g %jSjXj] - TD] i=l
(3)
j=l
where we are summing up, for non-satisfied constraints, the constraint dose less the target dose. Initially this sum is just -TD M. When all constraints are satisfied, B will equal zero. Now an iteration procedure begins that continues until either B = 0 or all seed locations are filled. Consider the set of all unfilled seed locations (initially all seed locations j = 1 to N), Xk. For each open seed location k, temporarily set Xk = 1 and Sk = S. Compute the infeasibility, now Bk, given by Eq. (3) above. Then set Xk back to zero, increment k and repeat. Pick k so that B1,is maximum and set Xk = 1 and Sk = S permanently (thus removing from the set Xk of all unfilled seed locations). In effect, find the seed location which does the most to satisfy all constraints. Initially, this will generally be a centrally located seed. Note that in Eq. (3) satisfied constraints are no longer considered as they no longer contribute to the sum of infeasibility. For example, after picking the first centrally located seed, there will now be a hole in the volume of unsatisfied constraint points, and the next iteration has to consider those unsatisfied constraint points that are left. Set S to the next available source strength from the source stock and repeat the above. When Bk = 0, or all source locations are picked, stop. If the algorithm stops because all source locations are picked, then some constraint points are not satisfied and an appropriate error message must be written out. Hotter sources may be required or a more closely spaced seed lattice must be considered. Otherwise a plan is generated which guarantees the target dose to all of the tumor volume. The solution time typically requires about 30 seconds on the VAX 750. The prescribed dose covers the tumor volume and can be easily verified by reviewing the dose distribution in each CT scan. Figures 3-6 are shown for a computer generated plan assuming 5.0 mCi I 125 sources and a target dose rate of 4 1.7 rad/hr. Figure 7 shows a plot in one CT plane for another patient case planned with the zero-one integer programming algorithm and 40 mCi 1125 sources. The isodose lines are shown with tick marks on the low dose side of the isodose line. The tumor contour has no tick marks. The projection to this plane of the needles and sources are also shown along with the head contour. Note that the target dose rate of 4 1.7 rad/hour does cover the tumor area.
Stereotactic brain implants 0 W. D. RENNER et al.
635
Fig. 5. Wire frame plot (again shown in black and white but otherwise in color) of tumor contours, implant needles, and Fig. 3. Wire frame plot (shown here in black and white but otherwise in color) of tumor contours and needles looking
sources with a view line perpendicular to the implant needles. The view line may be rotated about the needles. Similar plots can be made for any arbitrary view line.
down the implant needles. Plus signs mark existing needles. Similar plots can be made for any arbitrary view line. Although this algorithm does not necessarily find the mathematical minimum number of seeds and their corresponding optimal distribution, all constraints are strictly enforced and the solution found above will, in general, be close enough to an optimal solution for clinical purposes. This can be seen because the algorithm proceeds in a fashion similar to current practice. In general, the surgeon performing a seed implant attempts to distribute seeds throughout the tumor volume in a uniform
Fig. 4. Plot of tumor in intersecting plane (shown here in black and white but otherwise in color) perpendicular to the implant needles. Tumor is solid area (colored red), round circles indicate where tumor contours intersect this plane. Plus signs mark implant needles. A 1 cm grid is superimposed. Similar planes can be shown for any angle and depth.
distribution fashion until the combined contribution of dose from all seeds irradiates the tumor volume to at least the prescribed dose. The advantage of the above zero-one integer algorithm is that the infeasibility B provides a rational way to successively pick the seed locations, the entire task can be performed by the computer, the target dose is strictly enforced throughout the tumor volume, and the integrated dose outside the tumor volume is minimized. Note, however, that nothing is done to minimize the number of needles as sources which lay on the same lat-
Fig. 6. CT scan of patient case, with the dose area of 4 1.7 rad/ hour or more blotted out (shown here in solid white). CT image under the blotted area may be alternately viewed by pushing a key. A complete set of isodose curves may be plotted on the LA2 10 printer/plotter along with the head and tumor contours.
636
I. J. Radiation Oncology 0 Biology 0 Physics
Fig. 7. Isodose plot printed on the LA2 10 printer/plotter for a rather large tumor planned with the algorithm. Shown are 10, 20, 41.7, and 100 rad/hour isodose lines. The target dose rate was 41.7 rad/hour (1000 rads/day) and seven 40 mCi I 125
sources were found in the solution.
tice line in the needle direction are simply reported to be on the same needle. Ideally, the number of needles should be minimized in a computer solution of an implant plan, since each needle track can result in brain damage. This could be accomplished by performing a mathematical OR operation, as in Boolian algebra, among all seed variables along each needle, and adding up the result for all needles for a total cost to be minimized. The resulting objective function would not be a linear function, however, and that approach was not further explored. If the needle tracks found by the algorithm are not acceptable, the computer program will allow restricting possible seed location to specific needle tracks either specified by the user as described below, or kept from a first pass solution. The algorithm can then be rerun for just the specified needle tracks, with seed locations only picked from those on the available needle tracks. In this manner, the computer algorithm allows for an opportunity to reduce the number of needles.
Other software supportfeatures As indicated above, the number of needle tracks is not minimized in the computer solution, and so means must be provided for the user to edit needles and sources and to view resultant isodose distributions. With the joystick, the user may locate needles on any of three types of images: CT image, wire frame plot, or intersecting plane plot. A color wire frame plot of the tumor contours as seen looking down the needles can be produced on the imaging terminal (shown in Fig. 3 in black and white only). This view allows for locating a needle manually so that it does intersect the tumor volume, and the distribu-
April 1987, Volume 13, Number 4
tion of the needle tracks in a plane perpendicular to the needles is readily appreciated and may be modified. Needles may be located with the joystick and the target point (end point of the needle) is initially located by the computer program on the exit surface of the tumor volume. The needles may also be located with the joystick on a plot of a plane intersecting with the tumor volume (Fig. 4). This is a plane perpendicular to the needles running through the tumor volume at selected depths, showing the shape of the tumor volume model in that plane. The use of this type of graphic representation is described in more detail by Kelly, et a1.5 In our implementation here, Fig. 4, we have constructed our tumor volume in terms of square cubes, voxels, scaled to the size of the CT image pixel. Hence each of three sides of the voxel is 1.5 mm in length for this case. Planes intersecting these cubes will leave triangular patterns along the edge of the tumor volume. As described below, one can view the coverage of the tumor area shown by overlaying a representation of computed dose. The user may also specify sources on particular needles, with a wire frame plot of a view perpendicular to the needles (Fig. 5). The user indicates the distance along the needles from the target end for the source placement. A negative distance will extend the needles deeper. At printout time, all needle target points are placed one centimeter beyond the end source. The printout of the stereotactic coordinates of all needles and the loadings along each needle is accomplished from a single command from the imaging terminal. A printer/plotter must be present for presenting this information to the surgeon. For any CT scan or intersecting plane, the dose distribution may be viewed for the planned implant (Fig. 6). A solid area for which a selected dose is equaled or exceeded is colored in on the overlay plane of the imaging terminal, and may be alternately switched off by pushing a button on the imaging terminal leaving the underlying CT image (or intersecting plane plot). Different doses are viewed separately by sequentially specifying the dose or dose rate. This method of displaying was chosen due to the limited graphics resolution capabilities of the imaging terminal, graphic pixels being the same size as the image pixels, and the need to easily appreciate the region covered by the dose. Traditional looking plots of line isodose contours may be produced rapidly on the printer/ plotter from a single command at the imaging terminal, and includes the tumor and head contours found from the CT scan (these plots are only available for the CT scans). Fig. 7 shows such a plot.
Accuracy of the BR Wsystem The ability of the BRW system to accurately locate a point in space found on a CT image was tested by scanning a plastic phantom with small embedded air holes about 1 mm in diameter. Coordinates in the stereotaxic system were computed by one of three methods: (a) The
631
Stereotactic brain implants 0 W. D. RENNER et al.
HP calculator and program supplied with the BRW system was used with numbers generated with the track ball on a GE9800 CT scanner. This included calibrating rod coordinates and target points. The image was 5 12 by 5 12 pixels, one pixel 0.8 mm in size with a slice thickness of 1.5 mm. (b) These same numbers were inputted as data into the program described in this report. (c) The joystick on the Micro-Delta terminal was used to locate the same points on the image. The image on the Micro-Delta had been reduced to 256 by 256 pixels, pixel size now 1.6 mm. A needle was then adjusted to the center of each hole and the spatial coordinates measured by transferring to the BRW base assembly, and then compared to that calculated above. The results for three target points are shown in Table 1. From the data, our program did better than the HP calculator with the same numbers. This is most likely because we solved for the transformation of coordinates matrix between the stereotactic system and image system with least squares using all of the available data, whereas the HP calculator uses a geometric solution. The Micro-Delta, however, did worst, most likely because the image pixel size was increased to 1.6 mm in reducing the 512 by 512 image down to a 256 by 256 image, thereby increasing the uncertainty introduced by the pixel size. Hence pixel size will have an effect on the
Table 1. Distance in mm between target point computed from the CT image and that measured with the BRW Stereotaxic System HP Target point 1 2 3
calculator 1.0 1.3 1.4
VAX with same numbers
Micro-Delta with joystick
0.7 0.5 1.2
1.5 2.1 2.0
accuracy. This was further demonstrated when a sliced thickness of 1.Ocm was tried. Errors on the order of 1 to 2 mm were found in the plane of the image, but increased in the direction perpendicular to the image to as large as 6 mm. A small slice thickness is essential for good accuracy in all three dimensions. placement
CONCLUSION The software algorithm described above for supporting stereotactic brain biopsy and implants allows for the rapid planning of the implant after the patient is scanned with the BRW system. It is possible to plan the implant as the patient is prepared for the surgical procedure, with the computer generation of needle track coordinates for the BRW system, along with isodose curves in selected CT scan planes.
REFERENCES Bouzaglou, A., Dyck, P., Solt-Bohman, L.G., Gruskin, P.: Stereotactic interstitial implantation of brain tumors. Endocurie. Hypertherm. Oncol. 1: 99-l 12, 1985. Christiansen, H.N., Sederberg, T.W.: Conversion of complex contour line definitions into polygonal element mosaits. Comp. Graph. 12: 187-192, 1978. Hartmann, G., Schlegel, W., Sturm, V., Kober, B., Pastyr, O., Lorenz, W.: Cerebral radiation surgery using moving field irradiation at a linear accelerator facility. Int. J Radiat. Oncol. Biol. Phys. 11: 1185-l 192, 1985. Heilbrun, M.P., Roberts, T.S., Apuzzo, M.L.J., Wells, T.H., Sabshin, J.K.: Preliminary experience with BrownRobert-Wells (BRW) computerized tomography stereotaxic guidance system. J. Neurosurg. 59: 2 17-222, 1983.
Kelly, P.J., Kall, B. A., Goerss, S.: Computer simulation for the stereotactic placement of interstitial radionuclide sources into computed tomography defined tumor volumes. Neurosurgery 14: 442-448, 1984. Renner, W.D., O’Connor, T.P., Paulauskas, N.M.: Computer assistance in planning radiotherapy seed implants. Int. J. Radiat. Oncol. Biol. Phys. 5: 427-432, 1979. Rossman, K.J., Shetter, A.G., Speiser, B.L., Nehls, D.: Stereotactic afterloading Iridium implants in treatment of high-grade astrocytomas. Endocurie Hypertherm. Oncol. 1: 49-57, 1985. Zionts, S.: Linear and Integer Programming, Englewood Cliffs, NJ, Prentice-Hall. 1974, pp. 320-475.