APPLICATION OF IMAGE ANALYSIS METHODS FOR ISOCENTER QUALITY ASSURANCE IN RADIOTHERAPY

One of the major procedures for testing the geometrical accuracy of devices used in radiotherapy treatments is the test of the geometrical position of the radiation isocenter. The importance of the test reflects the fact that geometrical position of the radiation isocenter generally affects the tumor targeting. At present the geometric accuracy is assessed with the Winston-Lutz test which checks the position of an image of a ball marker with the respect to the center of the radiation field as projected on a detector plane. Obviously, determination of coordinates of a single marker is not sufficient to fully account for a complicated geometry of a therapeutic device. The purpose of the study was to design a new image analysis tool to better determine the isocenter. The proposed automated procedure for determining isocenter position uses projection data acquired for a special cube phantom. The projection images of a phantom are acquired for various angles of rotation of the gantry. A procedure is proposed to extract some geometric characteristics of a therapeutic device from the projection images.


INTRODUCTION
For all radiation dose delivery techniques in external beam therapy (EBT), the accuracy of the treatment is the ultimate goal.However, all assumptions of the perfect therapeutic procedure are in fact impossible to reach.Inevitable uncertainties always appear and could be identified at each step of the preparation and implementation of therapeutic process.Assessment of potential errors should include all aspects of overall performance of imaging devices used for patient positioning and target localization, accuracy of the treatment planning system, and geometrical precision of the treatment device with meticulous check of dose delivery.Advanced modes of high-precision radiotherapy (e.g., intensity modulated radiation therapy IMRT or volumetric modulated arc therapy VMAT) allow conformal delivery of radiation dose to a planned target volume while minimize exposure to organs at risk.On the other hand, the more advanced and precise the device or method is, the more attention should be paid to ensure its high accuracy and the lesser uncertainties are clinically accepted.Therefore, each radiotherapy department develops quality control tests to maintain high-level accuracy of the treatment.Various types and frequencies for providing such tests are recommended by national legislations.One of the major geometric tests is the test of the geometrical localization of the radiation isocenter.
The geometrical localization of the radiation isocenter generally affects the tumor targeting and overall accuracy of the treatment (Du et al., 2010).On the other hand the mechanical isocenter of the linac is the point in space where rotation axes of gantry, collimator, and treatment table intersect (Lutz et al., 1988).Ideally these three elements rotate around one point which additionally coincides with the radiation isocenter.Furthermore, this assumption implicates that the real position of radiation field during the treatment session is in accordance with a treatment plan.However, due to inaccuracies in the construction and real weights of device components (also the jaws of linac), the isocenter might slightly move in the space in a virtual spherical volume, which gives obvious mechanical limitation of the linac.Therefore, the test of radiation isocenter localization should be considered as very important.

NIEDŹWIECKI M ET AL: Quality assurance in radiotherapy
At present the discrepancy between mechanical and radiation isocenter is tested using the Winston-Lutz (W-L) test (Lutz et al., 1988).The W-L phantom is a small metallic ball representing the planned target.The phantom is positioned in 3D space using treatment room lasers, which should ideally point to the radiation isocenter.Then, projection images of the phantom are acquired and discrepancy between radiation and mechanical isocenter is reported whenever the center of the ball -as projected onto detector plane -deviates from the center of the radiation field.Measurements are repeated for a few angular positions of a gantry of a therapeutic device (Rowshanfarzad et al., 2011).
Obviously, comprehensive characterization of the geometry of a therapeutic device is not possible with a limited data provided by the W-L test.For example, the distance between the center of an image of a ball and the center of the radiation field depends on the distance between the source of ionizing radiation and the detector plane which is known only approximately.
However, the images of a W-L phantom certainly contain more data than explored in a simple W-L test.Thus, the purpose of the study was to design a new image analysis tool to support the isocenter quality assurance by extracting additional parameters from images of a W-L phantom.

METHODS
The problem of determination of the position of an isocenter has been defined in a study of Du et al. (Du et al., 2016).According to (Du et al., 2016), for a set of gantry orientations (specified by an angle  of gantry rotation around a horizontal axis of a therapeutic device) a line referred to as a central axis (CAX), which is normal to the detector plane and passes through a source of ionizing radiation must be determined (Fig. 1).Fig. 2. Images of a ball marker and a source in the detector plane.
Then the position P(x,y,z) of the isocenter is the solution to the following minimization problem (Du et al., 2016): where d(Q(x,y,z), CAX i ) is the distance from some point Q(x,y,z) in a 3D space to line CAX i and the sum on the RHS of Eq. ( 1) is minimized over all possible 3D points Q(x,y,z).
All coordinates are computed in a natural frame of reference with the origin in the center of the ball marker and following three axes: horizontal axis of the gantry rotation, vertical axis and the third axis normal to the previous two.Then, the position of the Image Anal Stereol 2017;36:35-41 isocenter is calculated with the respect to the position of the ball center.Because in clinical settings the phantom is positioned using the treatment room lasers, the computations provide data sufficient e.g. to recalibrate the lasers.
CAX is fully specified by its orientation angle  (i.e. the gantry has only one rotational degree of freedom) and an anchor point A -a projection of the phantom ball marker onto CAX.Both the anchor point A and angle α must be determined from phantom image data to find the position of the radiation isocenter, according to Eq. ( 1).
Denoting the source to detector panel distance by L, the ball to detector panel distance by L-l, and ball to anchor point A distance by d it follows from basic geometrical considerations that: where both D and  (Fig. 2) can be measured as in clinical settings the pixel size is known.
Consequently, the coordinates of the CAX anchor point A are: and the directional vector of CAX is fully specified by .Clearly, to find any CAX and to calculate the radiation isocenter from Eq. ( 1), both the ratio l/L and  must be determined from image data (Fig. 3).In the following we focus on solution to this problem.

Fig. 3. A typical projection image of a cubic phantom with a ball marker inside. A radiation center is by definition the center of the crosshair. Any shift of the ball with the respect to the radiation center indicates imperfect geometrical conditions.
The ratio l/L can be determined from the analysis of vertical profiles of projection phantom images (i.e.profiles parallel to the gantry rotation axis).Indeed, assuming geometry like in Fig. 4 (where s is the size of the cubic phantom and s 1 + s 2 = s, and the phantom is not necessarily ideally positioned i.e., some shift from a CAX to the phantom center is accepted) it can be shown that the distance r(x), the ionizing radiation travels through the volume of the phantom is: ) and for l equal approximately to L (which is the case of clinical settings adopted in this study, as explained in the next section) the dependence of r(x) upon x can be reasonably assumed to be linear.Indeed, the extent of this range of x values is approximately equal to s 2 /L which for the adopted clinical settings is only a small fraction of s.Replacing in Eq. (4) x by x min + and expanding r(x) around x min with a Taylor's series it can be shown that the nonlinear terms of the Taylor's series can be neglected.Moreover, due to Beer-Lambert law the logarithms of gray values are proportional to r(x).Then it follows from triangles similarity that if we NIEDŹWIECKI M ET AL: Quality assurance in radiotherapy measure the width Z of the vertical profile of the logarithm of gray values at the half of the profile height we have: (5) The rotation angle  can be determined from the analysis of horizontal profiles of projection phantom images (i.e., profiles normal to the gantry rotation axis).The geometry is however more complicated in this case (Fig. 5) as all possible scenarios of rays passing through the phantom must be considered.Although closed form expressions can be found also for this case we choose to approach the problem numerically by analyzing the limiting cases of X-rays sliding through the edges of the cube (dot-dashed lines in Fig. 5).

Fig. 5. The geometry used for the calculation of the attenuation of the ionizing radiation along profiles normal to the gantry axis.
It follows that the dependence of the distance r(x), the ionizing radiation travels through the volume of the phantom on the position along the horizontal profile can be approximated very well with a piecewise linear function (Fig. 6) and from the analysis of the asymmetry of the profile the rotation angle  can be extracted.

RESULTS
In the study images of a phantom were acquired using PaxScan 4030 system and Acuity Simulation System (Varian Medical Systems, Salt lake City, Utah).The phantom was a cube with 70 mm edge length and 2mm tungsten ball in the center.The system was operating at 60 kV, 80mA current, 5 ms exposition time and pixel size 0.194 x 0.194 mm.The spectrum of X-rays emitted by an X-ray tube fort such settings contains a strong characteristic line which justifies application of the Beer-Lambert low in the analysis of intensity profiles.The images were acquired for preset gantry rotation angles from 0 to 180 degrees with 15 degrees step.The values of rotation angles as displayed in the operator panel were different from the preset values by no more than 0.05 degree.The magnification of images was selected in the operator console of the simulator to achieve ratio l/L equal to 1.For such settings the imaging system automatically resizes pixels of the acquired images to sizes corresponding to a panel detector crossing the center of a phantom.Sample vertical and horizontal profiles through the phantom before and after taking natural logarithms of gray levels are shown in Fig. 7 and Fig. 8, respecttively.Because the transition region between foreground and background is very narrow for vertical profiles, in Fig. 7 we show for clarity only halves of the profiles.The shape of the profiles confirms qualitatively the analysis presented in the previous section.
In the automated analysis of the profiles we focused on profiles crossing the center of the ball as projected onto the detector plane.Phantom images were segmented with threshold set to the value at the half distance between the two modes corresponding to the background and to the phantom.Then the ball was extracted as the black connected component closest to the center of mass of the white region.Next, the center of mass C of the ball was computed and vertical and horizontal profiles crossing C were extracted from the original images and further analyzed after taking natural logarithms of gray levels.The data points inside black regions in segmented phantom images were discarded in the analysis.
The profile (either horizontal or vertical) was split into five regions: the left-most and the right-most modeled by a constant function and the middle three regions modeled by a linear function (Fig. 11).Prior to the split the phantom region was extracted from an original image by thresholding it with a threshold set in the middle between the histogram modes corresponding to the foreground and the phantom regions.Then the split points P 1 , P 2 , P 3 , and P 4 were found by a brute force search minimizing the total least squares fitting error under a constraint that the split points between constant and linear components of the model (P 1 and P 4 ) must be within the foreground region.Given the split points P 1 , P 2 , P 3 , and P 4 the height of the vertical profile was determined separately for the left and the right wing of the profile as the difference between ordinates of P 1 and P 2 and between P 3 and P 4 .Then the width of the profile was set equal to the difference between the position of points at the half of the left height and at the half of the right height.The calculated values of the ratio l/L are shown in Fig. 12 for eleven tested gantry orientations.On average the computed ratio l/L was equal to 0.999 with 0.0014 standard deviation compared to nominal value equal to one.
The calculation of the rotation angle of the gantry required some more attention.For every angle in the range ±2 degrees around the nominal gantry rotation angles with 0.01 degree step we numerically generated profiles like in Fig. 6 but with a minus sign (to account for a minus sign present in the Beer-Lambert low).This process is straightforward as it requires calculation of coordinates of cube corners after rotation and calculation of the intersection points between the rays of ionizing radiation and phantom faces.In the calculations we used nominal values of L = l = 100 cm which is the preset option in the PaxScan system for Winston-Lutz test and 70 mm phantom cube size.Then the generated profiles were compared with the profiles computed from the image.The angle corres-ponding to the generated profile fitted best to the computed profile was selected as the calculated gantry rotation angle.The calculated values of the rotation angle are shown in Fig. 12 for eleven tested gantry rotations.The coefficient of correlation between nominal and computed gantry rotation angle was equal to 0.09997.The slope of the best fit line was equal to 0.999 with 0.005 error.

DISCUSSION
In the study we addressed the problem whether it is possible to extract some useful information about the geometry of a therapeutic device from projection images of a cubic phantom.In a clinical practice the only information extracted currently is the position of a center of a phantom as indicated by a ball marker with the respect to the center of the radiation field.We have shown that using still such a basic phantom additional information can be extracted, namely the ratio l/L and the gantry rotation angle .We have shown that there is a very good agreement between the nominal and calculated values of these quantities which is an important conclusion from the viewpoint of quality assurance.The calculated values of the ratio l/L and the gantry rotation angle  can be in turn used for more precise determination of the position of the radiation isocenter from Eq. (1) as the uncertainties related to the discrepancy between preset and real values of these parameters may be eliminated this way.
The developed method has been used for the analysis of phantom images acquired with a classical simulator which uses kV imaging.Although the method can be directly used for MV imagers of medical linear accelerators, the quality of MV images is worse than the quality of kV images.We expect however that the developed method can be still useful for MV images for example if the exposition time is increased to suppress noise.
We have proposed an algorithm for automated analysis of the projection data to extract the aforementioned quantities of interest.The algorithm is based on image analysis techniques.The results of the analysis can be further used to define axes of radiation field and for the computation of the isocenter of radiation as defined in Eq. ( 1).
The results of the study prove that it is possible to relax assumptions about an ideal geometry of a therapeutic device.It follows from the results of this study that analyzing projections of more than one phantom point we can subsequently extract information about more geometric characteristics of the device.The amount of information about the geometry of a therapeutic device which can be extracted from images of commercially available phantoms used in routine clinical practice (e.g., the one used in the present study) is however limited to parameters evaluated in this study.

Fig. 1 .
Fig. 1.Axial cross-section of the geometry of a therapeutic device.

)Fig. 4 .
Fig. 4. The geometry used for the calculation of the attenuation of the ionizing radiation along profiles parallel to the gantry axis.For x in the range specified in Fig. 4 by the dashed lines (i.e., from a minimal value x min equal to ) 2 / ( 1 s L Ls 

Fig. 6 .
Fig. 6.Typical horizontal profile of the distance the ionizing radiation travels through the volume of the phantom.
The automated analysis started from extracting the position of the ball.The extraction of the ball was based on the histogram-based (Fig. 9) segmentation of the phantom images.In all cases the histogram contained three modes corresponding to the background (low gray Image Anal Stereol 2017;36:35-41 levels), foreground (high gray levels) and phantom (intermediate gray levels) regions.Vertical profiles of a phantom image before (a) and after (b) taking logarithms of gray values.
Horizontal profiles of a phantom image before (a) and after (b) taking logarithms of gray values.

Fig. 11 .
Fig. 11.Logarithm of gray level intensity profile modeled as a piece-wise linear function.Given the split points P 1 , P 2 , P 3 , and P 4 the height of the vertical profile was determined separately for the left and the right wing of the profile as the difference between ordinates of P 1 and P 2 and between P 3 and P 4 .Then the width of the profile was set equal to the difference between the position of points at the half of the left height and at the half of the right height.The calculated values of the ratio l/L are shown in Fig.12for eleven tested gantry orientations.On average the computed ratio l/L was equal to 0.999 with 0.0014 standard deviation compared to nominal value equal to one.
Ratio l/L (a) and rotation angle(b)as extracted from the phantom images.