Radiol Oncol 1996; 30: 58-65. A simplified technique for aligning radiation fields in portal imaging Hui Wang and B. Gino Fallone Medical Physics Unit and Department of Physics, McGill University, Montreal General Hospital, Montreal, Quebec, Canada A simplification of the chamfer matching technique is proposed far aligning radiation fields in portal imaging. In this application, the shaped treatment field is first aligned to the prescribed field by using low order geometric moments. Then the alignment is fine tuned with a simplified chamfer matching technique in which the cost function is successively minimized along coordinate directions. The viability of this minimization approach far the detection of radiation beam shaping errors in portal imaging is demonstrated through comparison with the downhill simplex method. When used in combination with lower geometric moments, this minimization method appears to offer advantages over the standard downhill simplex method in terms of precision and computation speed. Effects of two types of cost functions and edge distance maps in this application are also discussed. Key words: portal imaging; tomography, x-ray computed-methods; radiotherapy Introduction With the ongoing improvement in treatment planning accuracy, there is greater demand on the accurate implementation of a treatment plan. A significant amount of research has been carried out on every aspect in portal image verification, from image acquisition systems1-7 to post processing algorithms^10 and registration procedures.11"16 Current progress in this Correspondence to: Prof. B. Gino Fallone, Ph. D., FCCPM, ABMP, Medical Physics Unit and Department of Physics, McGill University, Montreal General Hospital, 1650 Cedar Avenue, Montreal, Quebec, Canada H3G 1A4; Phone: + l 514 934 8052; Fax: 1 514 934 8229. UDC: 616.149.66-073.756.8 relatively new field and its impact on clinical practice of radiation therapy have been systematically reviewed.17, 18 To set up a radiation therapy treatment, one needs to register portal images to a reference image in order to visualise the coverage of the target by the radiation beam. Because of the radiation required for portal image acquisition, it is not feasible to perform beam shaping and localization at the same time. The common practice is to shape the radiation beam by following the prescription before setting the patient up, then to direct the shaped beam to the desired target by positioning the patient. Accurate beam shaping is not only a factor determining the precision of beam coverage but is also beneficial for the interactive beam localization in this two step procedure. Extracting A .simplified technique for aligning radiation fields in portal imaging 59 anatomical landmarks for image registration is time consuming. On the other hand, the radiation field is a prominent feature in a portal image. A correctly shaped field can be automatically extracted and used as a registration landmark. Computer programs that can automatically align the portal and reference image to the field and show the relative position of the target in the field can serve as a tool for beam localization with an on-line imager. In this paper, we propose a simplified method for aligning radiation fields. A natural way to detect shaping errors is to align the prescribed field with the treatment field, and to evaluate the visual match. Among ali the techniques that have been employed in this task, chamfer matching is probably the most successful one due to its insensitivity to noise and to discrepancies in shapes. 19 A technical issue in chamfer matching is the minimization of a cost function of the geometric transformation parameters (translation offsets, rotation angle and scaling factor). Since there is no analytical expression for a cost function, the minimization has to rely on iterative searching techniques among which the downhill simplex method20' 21 and the Powell's method21 have been applied. In a recent study on multileaf collimator configuration verification, Zhou and Verhey reported that the downhill simplex method is not very sensitive to rotation and requires starting points close to the global minimum.22 To overcome these shortcomings, they used Hough transform and geometric properties of two contours to determine the starting point of chamfer matching. The downhill simplex method is straightforward and easy to implement but not very efficient in terms of the number of function evaluations that it requires. Powell's method is almost surely faster in ali likely applications.21 If chamfer matching is used in combination with lower order geometric moments, the minimization process will start very close to the global minimum.23 Because of this, standard minimization techniques which are designed for general purposes could be simplified to suit specific applications. In this paper, we propose a simplified adaptation of the chamfer technique for matching treatment and prescription field pairs. Materials and method Thirty pairs of simulation and double-exposure portal films randomly selected from our patient archives were digitized to 512 x 512 digital images with 8 bit contrast. The prescribed field contours were drawn with a mouse by following the prescription on the simulation images. The treatment field contours were obtained by using a "contour" operation on the field masks automatically extracted from the portal images. 15 The "contour" operation peels the mask (an object in a binary image) at the depth of one pixel with an "erosion" operation and then subtract the eroded mask from the original one to acquire the field contour. Edge distance map generation An edge distance map image E(i,j) (Figure lc) was generated from a binary image containing a prescribed field contour (Figure la). The value of a pixel in the edge map image is the distance from that pixel to the closest feature (contour) pixel. The higher the value of a pixel, the farther away it is from the contour points. This edge map image, resembling a landscape model with a valley along the prescribed field contour, will be used as a mould onto which the treatment field contour (Figure lb) will be fit. By analogy, the treatment field contour will have different gravitational potential energy at different locations and with different orientations (Figure ld). Registration will be achieved when the treatment field contour slides down into the valley under the action of gravity. Ideally, the value of a pixel in an edge map image should be the Euclidean distance from that pixel to the closest contour point. However calculating Euclidean distances is computationally intensive and yet may not be worthwhile because of digitization effect on the contours. The common approach to generate an edge map is to use distance transformation which approximates global distances by propagating 60 Wang H and Fallone B G Figure l. Illustration of procedure for chamfcr matching: (a) Simulator image; (b) Portal image; (c) Prescribed field contour obtained from (a); (d) Treatment field contour extracted from (a); (e) Edge distance map generated from (c); (f) After (d) is transformed with a trial set of parameters, it is overlaid on top of (e). local distances.24 By analogy, local distance propagation is similar to using a ruler of unit length for measuring long distances. As the ruler is being passed in steps, the number of passes is counted and the total number of counts is used as the distance from the starting point. A distance transformation passes a kernel which carries distances between a pixel and its neighbors over an image and assigns the accumulated distances to the pixels it passes by. Our edge distance maps were generated with the "chessboard" distance transformation25 which employes the 3 x 3 kernel (1) for representing the local distances and propagates the local distances in two steps. In the first step, the upper-left half of the "chessboard" kernel 1 1 0 (2) was passed over the contour image (in which ali the contour pixels were assigned the value of O and ali the non-contour pixels were assigned the value of co) from left to right and from top to bottom. At each pixel, this half kernel was added to that pixel and its four neighbors and the minimum value among the five was assigned to the corresponding pixel in an intermediate image. In the second step, the lower-right half of the "chessboard" kernel o 1 1 (3) was passed over the intermediate image from right to left and from bottom to top. The minimum value among the five was assigned to the corresponding pixel in the distance image. Cost function minimization After the edge map (Figure le) was generated from the prescribed field contour (Figure la), a geometric transformation with a set of trial parameters, where (a, b) is the translation vector, 0 is the rotation angle around the center of mass of the prescribed field, (Xp, Yp), and m is the scaling factor, was applied to the treatment field contour {(x,, y); i= 1, 2, ...N} (Figure lb). The transformed contour {(x;, yi); i = 1, 2, ...N} where - Xp y„ mcosfl -msinO m sin O m cosO ■ A' (4) was overlaid on the edge map (Figure ld), and the sum of the values of the pixels in the edge map along the transformed treatment field contour was used as the cost function (5) F (a, b, 0, m) is a function of the transformation parameters and has its minimum value when I A .simplified technique for aligning radiation fields in portal imaging 61 the correct transformation parameters relating treatment field contour to the prescribed one are used. Minimization of this cost function is usually achieved with iterative searching algorithms which start with a set of the initial trial parameters (a0, b0, 0O, m0), navigate in the 4 dimensional space formed by ali the possible transformation parameters, and check some criterion at each step until convergence is reached. We used a method similar to that used by Zhou and Verhey to obtain the starting point. Except for the rotation angle (which is set as O), the initial trial transformation parameters are obtained from the center of mass and the area of the fields: {aa.hu.Oa.ma) Xr - Xt, Y„ -Y,.O. (6) where (Xp, Yp) and Ap are the center of mass and area of the prescribed field, (X,, Y() and A, are those of the treatment field. 0O is set to O because the angle between the two contours is always very small before matching in our case. After this preliminary matching, the starting point is very close to the global minimum of the cost function. We minimize the cost function with respect to one variable at a time, one after another and cycle after cycle. One dimensional minimization is achieved by bracketing (Appendix). The cost function is minimized with respect to 0 first. The initial bracket is set to ( 0Q, 0O + 0.1). The bracketing stops when F( < 10-'1 where n is the number of iteration. Similar procedures are then carried out to a, b and m sequentially with initial brackets (fl0, an + 5), (bo, b0 + 5) and (m0, m0 + 5), respectively. All the procedures are stopped at the same precision. Once a bracketing procedure stops, the corresponding variable is kept at the convergent value. Since our stopping criterion is set very low, one cycle is sufficient to reach convergence for all the variables. Results and discussion Chamfer matching is a technique of pattern recognition type. The goodness of a match is characterized by a similarity measure. The measure we used is the average edge distance which is the minimum cost. It is defined as the average pixel value along the registered treatment field contour in the edge map image, which can be expressed as k F (am, bm, 0m, m,„), where N is the number of pixels along the treatment field contour and (am, b„„ 0„„ mm) is the final transformation parameters reached by the minimization. If the two feature being matched are perfectly similar to each other, the minimum cost should be zero. Otherwise, it will take on a positive value. When matching simple features like closed field contours, unless the shaping errors are extremely large, a smaller minimum cost should correspond to a better match. Minimization approach To test the viability of successive minimizations along coordinate directions in chamfer matching when combined with geometric moments, we compared the average edge distance obtained with successive minimization along coordinate directions, ds.m.a.c.d. to that obtained with the downhill simplex method Dsimpiex. A flowchart of the comparison test is shown in Figure 2. Figure 2. Flowchart of thc comparison test bctwccn thc simplex method and succcssivc minimization along coordinate directions (S.M.A.C.D). 62 Wang H and Fallone B G Throughout the comparison tests, both minimization procedures are started after the first step match obtained with geometric moments in Eq. (6). The vertices of the initial simplex are set to (ao, bo, 0„, m0), (a0 + 5, b0, 60, m0), (ao, bo + 5, 0o, mo), (ao, bo, Bo + 0.1, mo) and (ao, bo, 0o, mo + 0.1), where a0, b0, 6o and m0 are given in Eq. (6). The absolute values of the minimum costs obtained with the two methods are very close. In order to see the difference clearly, we calculated the relative difference in percentage pi(I (Ps l()(). " ¡¡in (7) AVERAGE EDGE DISTANCE 6 * i 32 1 24 U 15 16 -64 -56 -48 -40 -32 -24 -16 -S O 8 16 24 32 40 48 56 64 (b) Arithmetic avg. = -2.271 % vs - - RMS - ........ « r -64 -56 -48 -40 -32 -24 -16 -8 O 8 16 24 32 40 48 56 64 .o I 4 z -64 -56 -48 -40 -32 -24 -16 -8 0 S 16 24 32 40 48 56 64 Percentage Difference Figure 3. Comparison of final average edge distances achieved with different approaches. (a) successive line minimization versus the simplex method with "chessboard" edge distance maps and arithmetic average edge distance as cost function. (b) arithmetic versus root of mean square average edge distance as cost function when "chessboard" edge maps and successive line minimization were used. (c) "chessboard" versus "5-7-11" edge maps when arithmetic average edge distance as cost function was minimized successive line minimization. a histogram of which obtained from the 30 cases is plotted in Figure 3a. It can be seen that, on the average, successive line minimization along coordinate directions can achieve a smaller residual value of the cost function therefore has better precision than the simplex algorithm in this matching scheme. Successive minimization along coordinate directions is actually the first step of the Powell's method.26 Starting from the coordinate direction set, Powell's method successively minimizes a function in each direction in the set and adjusts the direction set after each cycle of line minimization through ali the variables. The adjustment of searching directions is to handle functions which can not be well approximated by quadratic forms because Powell's method is based on Taylor expansion. For example, a function having a long narrow valley will make successive minimization along coordinate directions very inefficient. However, when simple geometric objects like the field contours, are being matched, it is very unlikely that the cost function will have such erratic behavior. Moreover, after preliminary matching has been done with geometric moments, the starting point is very close to the global minimum of the cost function. In this neighborhood, the cost function can be well approximated by a quadratic form which can be exactly minimized by one pass of line minimization through all the variables. Computation speed The two algorithms are also compared in term of computation speed. Table l lists the numbers of iterations to reach convergence. Convergence is defined in the following manner: for successive line minimization, either the precision criterion is satisfied or the values of the cost function at the two ends of the bracket no longer change; for simplex algorithm, either the precision criterion is satisfied or the values of the cost function at the vertices of the simplex no long change. It should be noticed S.\l..\.(. O x A simplified technique far aligning radiation fields in portal imaging 63 Table l. Numbers of iterations for reaching convergence in successive line minimization along coordinate directions (Ns.m.a.c.d.) and the downhill simplex method (N,,-„,p/l,r). The average difference (Ns.m.a.c.d. - Nsimplcx)~ -3-2 case number Ns.M.A.C.D. ^simplex Ns.M.A.C.D.-N..simplex 1 32 31 1 2 31 24 7 3 33 41 -8 4 32 40 -8 5 25 24 1 6 30 38 -8 7 35 47 -12 8 34 47 -13 9 24 33 -9 10 32 31 1 11 29 44 -15 12 29 19 10 13 34 34 o 14 32 35 -3 15 34 42 -8 16 28 25 3 17 36 44 -8 18 36 37 -1 19 32 34 -2 20 30 32 -2 21 33 33 o 22 23 33 -10 23 29 32 -3 24 26 23 3 25 31 37 -6 26 32 38 -6 27 30 31 -1 28 26 26 o 29 31 29 2 30 31 32 -1 that not only less number of iterations are required for sequential bracketing (3.2 on the average), but the number of calculations involved in a single iteration in bracketing is also much less than that in the simplex algorithm. The reason we used the numbers of iterations to measure computation speed is that the absolute computation time for either algorithm on our computer (Indigo, Silicon Graphics Inc., Mountain View, Calif. 94043) is actually very short since the search starts from close to the minimum cost. This comparison is to show that the better accuracy of successive minimization along coordinate directions as discussed previously is achieved without sacrificing computation speed. Effects of cost functions According to the investigation by Borgefors,27 the cost function based on the root of mean square average edge distance (summing up the • square of the pixel value along the contour that is being matched when it is fitted into the edge map) has fewer local minima than does the cost function based on the arithmetic average edge distance. It can, therefore, reduce the chance for minimization being trapped by false convergence. Local minima are not a problem within our approach since they are bypassed by the first approximation [Eq. (6)]. To investigate the effect of the two types of cost function in this application, we calculated the average edge distance, and observed a slight advantage of the arithmetic average type cost function over the one of root of mean square average type (Figure 2b). In chamfer matching every point on the treatment field contour contributes to the cost function. The root of mean square average type cost function increases faster around the minimum and therefore assigns more weight to the distorted parts of the treatment field contour making the minimum of the cost function not correspond as well to the actual match as does the arithmetic average type cost function. Effects of distance transformations The accuracy of chamfer matching also depends on the distance transformation used to generate the edge distance map. The difference in the edge distance map will be carried into the cost function. The "chessboard" distance transformation is one of the simplest approximations of the Euclidean distance transformation because it assigns the same distance from a pixel to ali its 8 immediate neighbors while the Euclidean distance from a pixel to its diagonal neighbor is Closer approximations can be achieved by (a) assigning different numbers to the orthogonal and the diagonal neighbors that better represent the 1 : ff ratio, e.g. the 3 x 3 kernel in the "3-4" distance transformation24 3 O 3 (8) 64 Wang H and Fallone B G 14 11 10 11 14 11 7 5 7 11 10 5 O 5 10 11 7 5 7 11 14 11 10 11 14 or by (b) using a larger kernel to take more neighbors into account, such as the "5-7-11" kernel24 (9) The deviation from the true Euclidean gauge may also result in local minima of the cost function when minimization starts far from the global minimum. Except for orientation, the treatment field contour in our case has been brought very close to the matched position by the first approximation [Eq. (6)]. To determine whether a closer approximation to the Euclidean distance transfarmation would have an advantage in our case, we compared the effect of the "5-7-11" kernel with that of the "chessboard" kernel. The average edge distance obtained with the "5-7-11" kernel has been divided by 5 to normalize the unit distance to one. Figure 3c shows that the "chessboard" kernel has a slightly better performance than the "5-711" one. This can be explained rather simply. An edge distance map generated by a distance transformation ("chessboard" or "5-7-11") is very accurate in specifying the orientation of a feature. When the minimization with respect to the rotation angle is completed, the treatment field contour is very close to the final registration. The advantage of the "5-7-11" kernel only occurs when the searching process starts at a point which is far from registration. This situation does not occur in our case because we start with a relatively close approximation. On the other hand, the slopes in the edge map generated with the "5-7-11" kernel are steeper so that shaping errors and noise on the treatment field contour will have more weight in the cost function weakening the correspondence of the minimum of the cost function to the actual match. In conclusion, our results show that incorporating additional geometric information into chamfer technique may improve performance without sacrificing computation speed when matching simple features. In this kind of tasks, successive line minimization of the cost function along coordinate directions is viable when chamfer matching follows a preliminary alignment. Appendix: Minimization by bracketing Given a continuous unimodal functionf(x) (Fig. 5) in the region (a, b), the minimum can be reached by simple bracketing with a pair of points a < < b and a < x\"1 < b, where f (x}I)) < f (x\")). Here the superscripts at l and h denote "low" and "high", respectively, and the subscript i is the iteration number. Starting from an initial pair of points xb') and x^'), the function is evaluated at the reflection point of xb"1 about xbl), which is x* = xb") + 2 (xb') -xb">), and the value of the function at this point f (x*) is compared with /(xbl)) and /(xb")). Then the following assignments are made: r<" = .r " and A"'' = 4", ¡f/U') < /(*'.") : r(" - r< ."and ,<*>=* •. if/( .ri") /( *• ) < /(*'»') : i" = ,<"»nd = ,<"> + -0'") .if /(A") < /< ,• 1 -A precision value e is preset at 104 for our calculations, and if 1 the process stops; otherwise the process continues with another iteration. The programing syntax is shown below: / Ih iteration if 1 /('I") 1 t. then .r,,M„ = .r, and terminate; else, i + I th iteraiion i f .A x- )) ) ('t" else, then W (<, W( C l | . , . it j---l <, then .r„,,,, = .r, and lemunate; I /( C.) I else. i + '2th iteration. A .simplified technique for aligning radiation fields in portal imaging 65 Figure 4. Illustration of the bracketing process for minimizing a one dimensional function. References 1. Wolfe L, Kalisher L, Considine B. Cobalt-60 treatment field verification by xeroradiography. Am J Roentgenology 1973; 118: 916-8. 2. Baily N, Horn R, Kampp T. Fluoroscopic visualization of megavoltage therapeutic x-ray beams. Int J Radiat Oncol Biol Phys 1980; 6: 935-9. 3. Meertens H, van Herk M, Weeda J. A liquid ionization detector for digital radiography of therapeutic megavoltage photon beams. Phys Med Biol 1985; 30: 313-21. 4. Lam K, Partowmah M, Lam W. An on-line electronic portal imaging system for external beam radiotherapy. Br J Radiol 1986; 59: 1007-13. 5. Leong J. Use of digital flouroscopy as an on-line verification device in radiation therapy. Phys Med Biol 1986; 31: 985-92. 6. Shalev S, Lee T, Lesczcynski K, Cosby S, Chu T. Video techniques for on-line portal imaging. Computerized Medical Imaging and Graphics 1989; 13: 217-26. 7. Munro P, Rawlinson JA and Fenster A. A digital fluoroscopic imaging device for radiotherapy localization. int J Radiat Oncol Biol Phys 1990; 18: 641-9. 8. Leszczynski KW, Shalev S, Cosby NS. The enhancement of radiotherapy verification images by an automated edge detection technique. Med Phys 1992; 19: 611-21. 9. Crooks I, Fallone BG. Contrast enhancement of portal images by selective histogram equalization. Med Phys 1993; 20: 199-204. 10. Moseley J, Munro P. Display equalization: A new display method for portal images. Med Phys 1993; 20: 99-102. 11. Bijhold J, Gilhuijs KGA, van Herk M, Meertens H. Radiation field edge detection in portal images. Phys Med Biol 1991; 36: 1705-10. 12. Lam WC, Herman KS, Lam KS, Lee DJ. On-line portal imaging: computer-assisted error measurement. Radiology 1991; 179: 871-3. 13. Jones SM, Boyer AL.Investigation of a FFl-based correlation technique for verification of radiation treatment setup. Med Phys 1991; 18: 1116-25. 14. Balter JM, Pellizari CA, Chen GTY. Correlation of projection radiographs in radiation therapy using open curve segments and points. Med Phys 1992; 19: 329-34. 15. Wang H, Fallone BG. A robust morphological algorithm for automatic radiation field-extraction and correlation or portal images. Med Phys 1994; 21: 237-44. 16. Moselcy J, Munro P. A semiautomatic method of registration of portal images. Med Phys 1994; 21: 551-8. 17. Boyer AL, Antonuk L, Fenster A, van Herk M, Meertens H, Munro P, Reinstein LE, Wong J. A review of electronic portal imaging devices (EPIDs). Med Phys 1992; 19: 1455-66. 18. van Herk M, Gilhuijs K, Williams P, Cinquin P, Cionini L, Swindel W. Role of electronic portal imaging in high dose/high precision radiotherapy. Radiother Oncol 1993; 29: 269-70. 19. Leszczynski K, Loose S, Dunscombe P. A chamfer matching technique for evaluation of radiation field placement in portal images. In: Proceedings of COMP/CMBES Conference. Ottawa, Canada: 1993. 20. Nelder J, Mead R. A simplex method for function minimization. Compt J 1964; 7: 308-13. 21. Press WH, Flannery BP, Teukolsky SA, Vettcr-ling WT. Numerical Recipes in C. 2nd edn. Cambridge: Cambridge University Press, 1992. 22. Zhou S, Verhey L.J. A robust method of multilcaf collimator (mlc) leaf-configuration verification. Phys Med Biol 1994; 39: 1929-47. 23. Wang H, Fallone B. A sequential minimization technique in chamfer matching. In: Pro. 3rd Int. Workshop On Electronic Portal Imaging. San Francisco, Calif.: 1994. 24. Borgefors G. Distance transformations in digital images Compt. Vision Graphics Jmage Process 1986; 34: 344-71. 25. Rosenfcld A, Kak AD. Digital Picture Processing. New York: Academic Press, 1982. 26. Press WH, Flannery BP, Teukolsky SA, Vetter-ling WT. Numerical Recipes in C. 2nd edn. Cambridge: Cambridge University Press, 1992. 27. Borgefors G. An improved version of the chamfer matching algorithm. In: Proceedings of the 7th Jnternational Conference on Pattern Recognition. Montreal, Canada: 1984.