Automatic detection of graticule isocenter and scale from kV and MV images

Abstract Purpose To automate the detection of isocenter and scale of the mechanical graticule on kilo‐voltage (kV) or mega‐voltage (MV) films or electronic portal imaging device (EPID) images. Methods We developed a robust image processing approach to automatically detect isocenter and scale of mechanical graticule from digitized kV or MV films and EPID images. After a series of preprocessing steps applied to the digital images, a combination of Hough transform and Radon transform was performed to detect the graticule axes and isocenter. The magnification of the graticule was automatically detected by solving an optimization problem using golden section search and parabolic interpolation algorithm. Tick marks of the graticule were then determined by extending from isocenter along the graticule axes with multiples of the magnification value. This approach was validated using 23 kV films, 26 MV films, and 91 EPID images in different anatomical sites (head‐and‐neck, thorax, and pelvis). Accuracy was measured by comparing computer detected results with manually selected results. Results The proposed approach was robust for kV and MV films of varying image quality. The isocenter was detected within 1 mm for 98% of the images. The exceptions were three kV films where the graticule was not actually visible. Of all images with correct isocenter detection, 99% had a magnification detection error less than 1% and tick mark detection error less than 1 mm, with the exception of 1 kV film (magnification error: 3.17%; tick mark error: 1.29 mm) and 1 MV film (magnification error: 0.45%; tick mark error: 1.11 mm). Conclusion We developed an approach to robustly and automatically detect graticule isocenter and scale from two‐dimensionla (2D) kV and MV films. This is a first step toward automated treatment planning based on 2D x‐ray images.


| INTRODUCTION
In low-and middle-income countries (LMICs), radiation therapy has been shown to be a cost-effective therapy for many cancer treatments. 1,2 However, according to the IAEA-DIRAC data, 3 availability of radiation therapy in LMICs is extremely limited, and the shortage of radiation therapy staffing is significant. [4][5][6] A recent report 7 showed that by 2020, LMICs will have a deficit of around 12 000 radiation oncologists, 10 000 medical physicists, and 29 000 radiation therapy technologists. In response to the staffing shortage, automation of radiation treatment planning could potentially alleviate the staffing burden without compromising the quality of treatment in LMICs. [8][9][10][11] However, fully automated treatment planning is nontrivial. In particular, considering the resource limited settings in LMICs, many advanced technologies and equipment commonly available in the high resource setting countries are not available in LMICs. For example, three-dimensional (3D) simulation based on CT images are common in many countries, but many clinics in LMICs do not have access to a CT scanner, or have to limit the number of patients for whom they perform CT imaging. 12 Instead, conventional two-dimensional (2D) simulation is used. 13,14 2D treatment simulation images are typically taken using radiographic film, which are then viewed on a light box. In some settings, this has been replaced with flat panel x-ray detectors, but many centers still use film. Furthermore, in situations where the x-ray tube is out of service, it may be necessary to use mega-voltage (MV) imaging [using the radiation beams from the radiotherapy treatment device, and film or Electronic Portal Imaging Device (EPID)].
To use the radiographic films for automatic treatment planning, the first step is to digitize them. Traditionally this is done using a specialized film digitizer. Other options are available in resource scarce environments, such as using an inexpensive commercial flatbed document scanner, 15 although care should be taken such the digitized image does not distort the film.
Important steps in treatment planning using portal images typically include determination of patient treatment position, beam geometry, treatment isocenter, field limits, contours, and dosage. 10,16 For isocentric treatment, determining the isocenter of treatment field is a key step. The intended treatment isocenter is the primary reference location for radiation treatments and typically needs to be determined at the beginning of treatment planning to facilitate the subsequent planning tasks such as defining beam geometry and dose calculation. 17 In 2D simulation, the intended treatment isocenter can be determined by using a mechanical graticule when taking the 2D simulation images. The mechanical graticule is visible in the images and its isocenter can surrogate the intended treatment isocenter. To automate the treatment planning, it is necessary to automatically identify the isocenter and scale of the graticule shown on the 2D simulation images, which are likely to be created using a kilo-voltage (kV) x-ray tube and film, but could be created using the mega-voltage (MV) x rays from a linac or cobalt unit also, as described above. Inspired by the linac quality assurance (QA) of localizing the isocenter of radiation fields described by Du et al. 18 , we proposed an automated process by combining template matching, 19 Hough transform, 20 and Radon transform 21 to detect the isocenter and scales of graticule from kV/MV films. We tested our proposed automated detection algorithm on both scanned kV and MV films and digital MV images obtained from EPID. In this study, we showed the feasibility to automate the process of localizing isocenters and determining scales from kV/MV films with mechanical graticule, which will allow the automation of 2D radiation treatment planning.

2.A | Patient data
Under the approval of our institutional review board, 23 kV films,

2.B | Image preprocessing
Appropriate preprocessing on the original images should be done before the automatic isocenter detection and graticule magnification determination. The preprocessing was separated into 3 steps: perform correlation with a template; apply a linear weighting to correlated image; and perform nonmaxima suppression. The entire preprocessing steps are shown in Fig. 2, and the details of each step are described as below.

2.B.1 | Template correlation
The first step was to perform a correlation with a predefined template for the scanned images. 19 The correlation step found pixels on the scanned image that might be part of tick marks. Let I denote the scanned image and T the template, which was designed to detect the tick marks in the graticule. Based on the shape of the tick marks (crosses for kV films and dots for MV images), binary templates were created (Fig. 3). In the case that a graticule has tick marks of other shape, a different template matching the tick mark can be created as well and the algorithm can be tweaked for the new template. Both templates had a size of 1 cm × 1 cm with a pixel resolution of 50 DPI. In the cross-shaped template, the foreground value was 1 and the background value was 0, and the width of the two axis of the cross-shaped template was set to 1.2 mm based on the actual thickness of graticule axes. Letting (x 0 , y 0 ) denote the center of the  template and w the width of template, the intensity of the template was modeled with For the circular template, the intensity of the template had a radial Gaussian weight with a standard deviation (SD) of 1 12 mm, which was determined using the radius of the circles of graticule tick marks so that two SDs of the Gaussian equals to the radius of the circles. It was noted that the best results occurred when the width and radii of the templates corresponded to the width and radii of the graticule tick marks. Letting r denote the SD, 1 12 in this case, the intensity of the template was modeled with The normalized cross-correlation (NCC) 22 values were computed for each pixel in the scanned image by correlating with the template as: where T was the mean intensity value of the template T and I u;v was the mean intensity value of the image portion overlaid with the template. In order to account for magnification, the templates were scaled where T m ðx À x 0 ; y À y 0 Þ ¼ T xÀx0 m ; yÀy 0 m À Á , with (x 0 , y 0 ) the center of the template.

2.B.2 | Linear weighting
The following linear weighting was applied to C(x, y): here, C max was the largest correlation value in C. The linear weighting process filtered out most of the pixels that did not belong to tick marks. The assumption of the linear weighting was that the tick marks had higher correlation values with the template than the rest of the scanned image. In addition, this process also set a threshold for a pixel to be considered as a tick mark as the half of the overall maximum correlation value C max , and linearly rescaled correlation values to be above this threshold.

2.B.3 | Non-maxima suppression
Non-maxima suppression was performed to remove lines on the film that were not graticule axis. 24 The assumption was that the graticule axes should have the largest value after the correlation process. Let C″ denote the image after non-maxima suppression process, D the physical distance in centimeter between two nearest tick marks, and dpcm (dots per centimeter) the resolution of the scanned image (dpcm = DPI/2.54). The first step of the non-maxima suppression was to remove points which did not likely represent tick marks: The second step was to restore pixel values at the locations neighboring the nonzero values in C″ to those in C′. The neighborhood was defined to be within a Euclidean distance of 3 pixels: This step ensured that enough nonzero pixels remained in C″ for the subsequent Hough transform and Radon transform to detect graticule axial lines.

2.C | Isocenter detection
Detecting the treatment isocenter is equivalent to detecting the graticule axes. The intended treatment isocenter is the intersection of the two graticule axes. In general, the two graticule axes are perpendicular to each other and are aligned analogously to the x and y axis on a typical graph. We refer to the axis aligned roughly parallel to the x-axis on a graph as the horizontal axis and axis aligned roughly parallel to the y-axis on a graph as the vertical axis.
A combination of the Hough transform 20  ρ ¼ x cosθ þ y sinθ: Every nonzero pixel (x, y) in C″ was transformed to the Hough space (ρ, θ) using the Hough transform. To detect the horizontal axis, θ were restricted to 85 ∘ ≤ θ ≤ 95 ∘ . Let N(ρ, θ) denote the number of nonzero pixels (x, y) in C″ that were parameterized with (ρ, θ). Let (ρ m , θ m ) denote the parameters having largest N(ρ, θ) and (ρ m2 , θ m2 ) the parameters corresponding to the second largest value chosen to parameterize the horizontal graticule axis. Otherwise, the parameters (ρ, θ) that maximized the line integral R(ρ, θ) of the line represented by x cosθ þ y sinθ in C″ was chosen to parameterize the horizontal axis: To improve the speed, we used only the parameter sets of (ρ, θ) corresponding to the largest four values of N(ρ, θ) for the above line integral. For 85°≤ θ ≤ 95°degrees, values of ρ h and θ h that maximize R(ρ, θ) were used to parameterize the horizontal axis of the graticule. This procedure was equivalent to the Radon transform 21 except that only a small set of possible (ρ, θ) values were considered in order to reduce the computational time.
The vertical axis of the graticule was detected using the similar procedure for horizontal axis except that the θ value was restricted to [−5°, 5°]. The intersection of the vertical and horizontal axes then defined the treatment isocenter, noted as (x iso , y iso ). The graticule axes and isocenter detection is exemplified in Fig. 4.

2.D | Magnification and tick mark detection
Once the isocenter was identified, we were able to detect the tick marks from the preprocessed image C″. It is known that the physical distance between two nearest graticule tick marks is 1 cm, and this distance showing on the film is normally greater than 1 cm because of magnification during imaging. The magnification is determined by the distance between the physical graticule and the imaging plane, varying for different films and normally with a value between 1.1 and 1.5. 23 Unlike the EPIDs, the magnification factor is normally unavailable for a film. However, appropriate processing on the scanned films is able to recover the magnification, as described below.
The The coordinates of the tick mark points can be represented as The tick mark points determined by Eqs. (10) and (11) were used to select pixels from the processed image C″. A correct magnification value m should give a maximum average intensity value of all selected pixels. Formally, the optimal magnification, M c , was determined using the following optimization equation: where |P m | is the number of points in P m . Initial guess of the optimal magnification was done by sampling the value of m between 1 and 1.5 with an incremental of 0.01, which maximized Eq. (12). This step ensured that the following maximization algorithm would not converge to a poor local maximum. The initial guess was then used to initialize the golden section search and parabolic interpolation algorithm 25 (fminbnb function in Matlab, Mathworks, Natick, MA) to determine the optimum magnification M c .
Once the optimal magnification was determined, we performed the following step to further optimize the isocenter location.
The assumption here was that the optimal isocenter location should give the best detection of tick marks. We assumed that small perturbations (x ′ , y ′ ) added to the isocenter x iso ; y iso ð Þcould bring it to an optimal location, (X iso , Y iso ), that is, ðX iso ; Y iso Þ ¼ ðx iso ; y iso Þ þ ðx 0 ; y 0 Þ.
Following a similar procedure in determining the optimal magnification, an optimal solution of (x ′ , y ′ ) was determined using the following optimization function: ðx s ; y s Þ ¼ max The initial guess of (x ′ , y ′ ) was set to (0, 0) when the fminbnb function in Matlab was used to find the optimal solution of Eq. (15).
Once the optimal isocenter location and magnification value were determined, the tick marks were automatically determined by extending from isocenter along the graticule axes multiples of the magnification value, as shown in Fig. 5. The actual magnification of the image, M m , was estimated from the manually selected points as follows. The distance between two nearest tick marks was 1 cm in reality; therefore, the distance between two nearest tick marks on the image was equivalent to the magnification of the image and this distance was estimated by the following equation:

2.E | Approach validation
where d(·,·) represented the Euclidean distance of two points The isocenter was manually selected three times. The mean distance of any two selected points was used to quantify the intra-observer variability in selecting the isocenter, and used to compare with the isocenter detection error, which was defined as the Euclidean distance between the automatically detected isocenter and the geometric mean of the manually selected isocenters. The tick mark detection error was found by calculating the Euclidean distance between the automatically detected tick mark and the corresponding manually selected tick mark. The median of all tick mark errors in one graticule axis was used as the representative value of that axis. For each image, two tick mark errors were reported, one for the horizontal axis and one for the vertical axis. Both the isocenter error and tick mark error were scaled by dividing the magnification factor M m to reflect the actual physical error. In our automatic detection algorithm, we assumed that the tick marks were evenly spaced along graticule axes. However, of the images under testing, 11 MV films had uneven physical spacing between tick marks so that our automatic magnification and tick mark detection algorithm does not work. For these 11 films, only the isocenter detection error was evaluated.

| RESULTS
The proposed approach was tested on the 23 scanned kV films, 26 scanned MV films, and 91 EPID images. Overall, the processing time was fast. The typical runtime, including preprocessing, isocenter detection, magnification estimation, and tick mark detection, was approximately 0.2 s for an image of 800 × 1000 pixels on a computer with 2 GHz Intel Core i7 CPU and 8 GB memory.  Table 1, with details described as below.

3.A | Isocenter detection
The proposed approach was robust to kV and MV films of varying image quality. The isocenter could be detected in most images despite varied complexities such as occlusion and glare. The histogram of isocenter detection error is shown in Fig. 8(a), compared  with the intra-observer variability shown in Fig. 8(b). The isocenter was detected with accuracy less than 1 mm for all but three kV films where the graticule was not actually visible. These three kV films were illustrated in Fig. 9. Because the isocenter could not be automatically detected, these three cases were not included in the subsequent analysis.
The average isocenter detection errors for kV films (excluding the aforementioned three kV films), MV films, and EPIDs were 0.3 ± 0.2 mm, 0.4 ± 0.2 mm, and 0.2 ± 0.1 mm, respectively. For comparison, the corresponding intra-observer variability in manually selecting the isocenter was 0.3 ± 0.2 mm, 0.3 ± 0.1 mm, and 0.2 ± 0.1 mm, respectively, as shown in Table 1. The computer detection error was comparable to the intra-observer variability, which showed that the automatic isocenter detection had a similar performance of manual isocenter placement. Of all cases with successful isocenter detection, the maximum error was 0.8 mm for KV films, 0.8 mm for MV films, and 0.5 mm for EPID images.

3.B | Magnification and tick mark detection
Magnification and tick mark detection were applied to all those images with successful isocenter detection. The histograms of magnification and tick mark detection error are shown in Fig. 10. Of all images under evaluation, 99% had a magnification detection error less than 1% with the exception of one kV film, which had an error of 3.17%. The mean magnification error for kV films, MV films, and EPID images were 0.29%, 0.22%, and 0.18%, respectively (Table 1).
This result showed that the magnification could be automatically detected from images very accurately, no matter what images were scanned from kV or MV films, or were obtained through EPID.
Tick mark detection error was small as well for all images under evaluation. Ninety-nine percent of the images had a tick mark detection error <1 mm with the exception of one kV film and one MV film. The kV film, which had a magnification detection error of 3.17%, had a tick mark error of 1.3 mm for the horizontal axis. The error was due to the poor image quality. Many of the graticule tick marks in the image had very low contrast to the background, resulting to ineffective preprocessing for the automatic detection (Fig. 11).
On the other hand, the MV film, which had a magnification detection error of 0.45%, had a tick mark error of 1.1 mm for the vertical axis.
The mean tick mark detection error for kV films, MV films, and EPID images were 0.4, 0.6, and 0.2 mm, respectively (Table 1). These results showed that the tick mark detection was very accurate. Upon examination of some kV scans, we noticed that a slight asymmetry of the distance of nearest tick marks to the left and right of the isocenter were often observed. This was the most significant factor contributing to the tick mark detection error. The isocenter T A B L E 1 Summary of quantitative results (mean ± SD) of the isocenter detection error, intra-observer variability, magnification detection error, and tick mark detection error.   Notice that all EPID images had very good detection accuracy and results were better than either kV or MV films. EPID images were acquired digitally without the need of scan, which had much less uncertainty and the image quality was much better than the kV and MV films. Our proposed approach worked very well with all 91 EPID images. While the digital x-ray image acquisition becomes more accessible, the use of films becomes less. This implies that our proposed approach will be more robust when the image quality is no longer a concern.

Isocenter
As aforementioned, the conventional 2D simulation based on x-ray images is still the choice of many LMICs due to the limited resources. [12][13][14] The shortage of qualified staff for radiation treatment planning also presents challenge in delivering high quality radiation treatment in LMICs. 7 Our study has shown that isocenter on 2D x-ray images can be automatically detected by computers and automated process was comparable to manual selection conducted by human beings. This automated process will facilitate the development of automated treatment planning based on 2D x-ray images, which potentially can address the issue of staff shortage in LMICs.
This signifies an important application of our study in LMICs. On the other hand, the proposed approach can also facilitate the linac QA 26,27 by reducing the workload in analyzing the QA images for LMICs.
Our study has some limitations. One major limitation is the assumption that the tick marks were evenly spaced along graticule axes. For some x-ray films, there was geometric distortion during imaging so that this assumption was not true. In our study, the tick marks on eleven MV films could not be correctly detected due to this reason. Nevertheless, this issue did not present on the digitally acquired EPID images, which will be more accessible in the future. In addition, though the automatic graticule detection works very well, it has not been integrated into an automated treatment planning workflow for verification. As such, this work is a pilot study to verify the feasibility. Its clinical usability needs further validation. This will be our future study.

| CONCLUSION
We have developed an image processing approach to automatically and robustly detect graticule isocenter and tick marks from 2D x-ray images. Our results showed that the automated process was comparable to manual selection conducted by human beings. This essentially allows the automation of treatment planning based on 2D x-ray images. Together with automated treatment planning, this technique will have important applications in LMICs.