Development and performance assessment of an advanced Lucas‐Kanade algorithm for dose mapping of cervical cancer external radiotherapy and brachytherapy plans

Abstract Purpose The aim of this study was to verify the possibility of summing the dose distributions of combined radiotherapeutic treatment of cervical cancer using the extended Lucas‐Kanade algorithm for deformable image registration. Materials and methods First, a deformable registration of planning computed tomography images for the external radiotherapy and brachytherapy treatment of 10 patients with different parameter settings of the Lucas‐Kanade algorithm was performed. By evaluating the registered data using landmarks distance, root mean square error of Hounsfield units and 2D gamma analysis, the optimal parameter values were found. Next, with another group of 10 patients, the accuracy of the dose mapping of the optimized Lucas‐Kanade algorithm was assessed and compared with Horn‐Schunck and modified Demons algorithms using dose differences at landmarks. Results The best results of the Lucas‐Kanade deformable registration were achieved for two pyramid levels in combination with a window size of 3 voxels. With this registration setting, the average landmarks distance was 2.35 mm, the RMSE was the smallest and the average gamma score reached a value of 86.7%. The mean dose difference at the landmarks after mapping the external radiotherapy and brachytherapy dose distributions was 1.33 Gy. A statistically significant difference was observed on comparing the Lucas‐Kanade method with the Horn‐Schunck and Demons algorithms, where after the deformable registration, the average difference in dose was 1.60 Gy (P‐value: 0.0055) and 1.69 Gy (P‐value: 0.0012), respectively. Conclusion Lucas‐Kanade deformable registration can lead to a more accurate model of dose accumulation and provide a more realistic idea of the dose distribution.


| INTRODUCTION
Locally advanced cervical cancer is standardly treated using a combination of concomitant chemotherapy, external radiotherapy (EBRT), and brachytherapy (BRT) boost to the cervical region. 1 The current standard for dose accumulation of the combined radiotherapeutic treatment according to the International Commission on Radiation Units and Measurements (ICRU) report 89 is based on the simple DVH parameter addition without employing an adequate registration model. 2 Doses for absolute organs at risk (OAR) volumes are estimated by adding dose-volume histogram (DVH) parameters from each fraction, assuming that the location of a given hotspot volume (e.g., 2 cm 3 ) is identical in each fraction. 3 So the worst alternative should always be considered. For critical organs, this indicates that the maximum dose calculated in each fraction is always realized in the same volume of tissue. On the contrary, in the case of tumors, the minimum doses always meet in the same volume.
Changes in the patient's irradiation position between EBRT and BRT and the introduction of applicators close to the tumor before each fraction of BRT causes organ movement and soft tissue deformation. As a result, the specific tissue structures generally occupy a completely new position compared to previously applied radiation fields. The tissue structures bring previously accumulated doses to these new positions, which are often significantly different from those that would correspond to their current position. To correctly express the accumulated dose, the absorbed doses must be summed up not in the same spatial coordinates, but always in the same anatomical volume of tissue. 4 Therefore, a prerequisite for proper adaptive dose accumulation planning is the registration of shifts that occur within the tissues between individual fractions. These displacements are formally described by the deformation field. Appropriate anatomy imaging is the starting point for determining shifts between fractions. The general registration scheme involves the selection of a reference fraction, to which so-called floating images, corresponding to the remaining fractions of treatment, are subsequently related. Using a reference and a floating image, an appropriate registration algorithm can then be employed to determine the deformation fields. Knowledge of the deformation field enables reconstruction of the actual position of the tissues with respect to the applied radiation fields, thus, determining the appropriate contribution to the accumulated dose in each fraction. This procedure is referred to as dose mapping.
The use of the concept of dose mapping by image registration thus, enables the addition of doses in corresponding anatomical areas, even when applied in different fractions of combined radiotherapy.
The concept described above was used in this study to develop a software (SW) tool designed to optimize dose prescription in combined cervical radiotherapy. The application of such a tool to real tissue changes is merely an approximation of reality. Therefore, an integral part of the design of a model tool for dose mapping using image registration is careful testing of its statements.
Thus, this work aimed to create an SW tool containing an algorithm for deformable image registration (DIR) and to verify its possibilities for adaptive planning of combined EBRT and BRT treatment of patients with cervical cancer.

| MATERIALS AND METHODS
The study was divided into two parts. The setting of the parameters of the algorithm for the registration of computed tomography (CT) images and dose mapping of the combined EBRT and BRT was specified in retrospective study #1. The optimization phase of study #1, in which the optimal parameters of the registration algorithm were assessed, was performed on a group of 10 patients with cervical cancer. The second part of this study (study #2) assessed the performance of the algorithm by evaluating the accuracy of dose mapping on 10 other patients and compared the results with the current standard of dose accumulation in combined radiotherapy.
The fractionation parameters applied in the treatment of patients whose CTs were used in this study were EBRT 45 Gy in 25 fractions on planning target volume including nodes and then boost to the cervical region 6 Gy in three fractions using the volumetric modulated arc therapy (VMAT). After completing an external radiotherapy course, patients were treated with 28 Gy intracavitary BRT on highrisk clinical target volume in four fractions using a Fletcher Utrecht CT compatible applicator.
Compared to conventional radiotherapeutic procedures, combined radiotherapy is characterized by significantly different fractionation parameters during external and brachytherapeutic parts. For this case, ICRU report 89 recommends recalculating the absorbed dose distribution to the equieffective dose. The equieffective dose related to normofractionation is denoted by the symbol EQD2. 2,5 The following relationship was used to calculate EQD2 based on the applied doses: in which the ratio α=β was set to 10 Gy for tumor response and 3 Gy for late effect in critical organs.

2.A | Image preprocessing
The presence of an applicator in the reference image represents a major problem when registering floating images, taken as a part of EBRT planning, to a reference image corresponding to the BRT. In general, it is not possible to register a part of the data that is present in a single image. As per the procedure suggested by Moulton et al, 6 before the application of the deformable registration, preparatory steps were made to modify each BRT study. 2. The replaced area was blurred using a Gaussian filter.

2.B | Rigid and affine preregistration
Application of the algorithm for deformable image registration to significantly different CT studies of individual fractions of combined radiotherapy demonstrated that for its effective functioning, it is advantageous to first preregister all CT studies manually using rigid transformation to match the corresponding bone structures.
Improved image data matching was further achieved by using a relatively fast automatic affine preregistration of floating images to the reference image. The result of these steps was then followed by deformable registration to fine-tune the image similarity.

2.C | Lucas-Kanade deformable registration
Considering the results of studies 7-10 comparing the accuracy of different approaches for deformable CT data registration, the Lucas-Kanade (LK) algorithm was used. This method works with the concept of deformation as a continuous optical flow between CTs. The results of the above studies indicate that when modeling "physiologi-

2.D | Patient study #1registration parameters setting
The function of the selected registration algorithm can be further optimized for a specific application. The accuracy of image data registration is essential for the subsequent mapping of dose distributions of individual patients, as the same deformation vector field is used.
In study #1, the values and combinations of the number of image pyramid levels and size of the neighborhood around a given point were optimized.
Various metrics have been proposed to calculate the similarity of images and to assess registration accuracy. However, when evaluating registration errors under real conditions, it appears that none of the available criteria can be applied universally, and the application of a single criterion often leads to incorrect conclusions. An example might be the failure of image similarity evaluation procedures based on image intensity in areas of homogeneous intensity. 15 Regarding the recommendations of publications 6 and, 16 a combination of the following approaches was used to assess the accuracy of registration when optimizing the parameters.

2.D.1 | Landmarks distance
In the landmark-tracking technique, the degree of similarity is determined by the distance between the corresponding points marked in both the floating and reference images. When testing the algorithm, at least 100 well-recognizable matching points were manually identified in both CT studies by an experienced physician. After registration, the average 3D distance for each patient was calculated.
Since the area to which the BRT boost dose distribution was to be delivered was smaller than the area of interest for the EBRT, the selection of landmarks was not uniform across the CT study. A denser network of points (a minimum of five points in a given CT slice) was marked in slices affected by the BRT dose distribution. In other parts of the CT covered only by the EBRT dose, a smaller number of points were marked (always a minimum of two points in the slice).

2.D.2 | Mean squared error (MSE)
This approach is based on the comparison of the distribution of HU in the corresponding images using a certain mathematical or statisti- therefore, proper deformation in these areas is not guaranteed. Thus, a more robust method of comparing two intensity distributions was used, which, to some extent, combines both approaches and considers not only differences in intensity (or HU) at a given point but also spatial differences such as gamma analysis. 17,18 Gamma analysis was performed only on the patient's contour. For practical reasons, Verisoft v.7.1 software (PTW, Freiburg, Germany) was used to evaluate the agreement between individual CT slices by two-dimensional gamma analysis. The evaluation criteria were set to 2%/2 mm with reference to local "dose." An average gamma score specifying the percentage of image pixels that fulfilled the selected criteria, was determined for the given registration parameters.  of published studies, 6,7,9,10 indicating that both methods were sufficiently accurate in comparison with other freely or commercially available algorithms, the results were used as a reference to evaluate the accuracy of the tested registration model.

2.E.2 | Comparison with the current standard of dose accumulation
After mapping the dose distributions, the results can be presented as a summary DVH of EBRT and BRT treatment expressed in EQD2.

2.E.3 | Statistical methods
Welch's two-sample t test was used to compare the average dose differences at landmarks between the LK algorithm and the HS

3.A.3 | Gamma analysis
The evaluation of the image similarity by gamma analysis is presented in Table 3 and confirms the results of the previous methods.
When setting the gamma analysis parameters to 2mm/2%, the highest mean gamma score of registered CT studies of all patients was 86.7%.

3.A.4 | Evaluation summary
The results of the evaluation using all three selected criteria correspond to each other. Based on the average landmark distance, RMSE, and gamma scores, the best agreement was achieved for two pyramid levels and a neighborhood radius of 3 voxels. The parameters optimized in this manner were introduced into the LK algorithm and further used for deformable registration in patient study #2.
Regarding patient no. 5 (in Tables 1-3), the average distance of the corresponding points was doubled compared to the others, and the average gamma score after affine preregistration was <50%. Probably due to the significantly different distribution of the patient's tissues during the acquisition of the planning CT for BRT and EBRT, the presented SW model failed. Moreover, its application generated visually unrealistic deformations (Fig. 4). It should also be noted that no form of post-optimization regularization was used in the method.  The statistical analysis revealed a significant difference in the average dose differences at landmarks between the extended LK algorithm and HS (P = 0.0055) and modified Demons (P = 0.0012) algorithm. Therefore, the accuracy of the LK method appears to be slightly higher than that of others.

| DISCUSSION
The presented results indicate that the proposed SW tool, using the LK algorithm extended by the implementation of weighted windows and pyramidal representation of the image, can be used to estimate the cumulative dose of combined EBRT and BRT with higher accuracy than similar approaches. 7 F I G . 2. Difference between CT studies of patient no. 1 before registration (left) and after affine registration (right).
Difference between CT studies of patient no. 1 after deformable registration by optimized pyramid LK method with weighted window (pyramid level 2, neighborhood radius of 3 voxels). In patient study #2, the accuracy of dose deformation and mapping of dose distributions using the LK algorithm was evaluated.
According to statistical analysis, the results summarized in In addition to the possible advantages of the presented algorithm, the results of this work also showed the limits of its use.
Although the results indicate the potential use of the method to generate an accumulated dose estimate with the stated accuracy, its practical use significantly limits the need for manual removal of the applicator introduced into the patient's body in BRT applications.
This procedure could also introduce an additional error into the evaluation process.
Furthermore, it appears that the algorithm fails in cases of extreme changes in the position of the tissues between the floating and reference CT images. Unfortunately, these cases are not exceptional when EBRT is combined with BRT (results for patient no. 5 in Tables 1-3). In this context, it should be noted that the higher robustness of the presented tool against extreme deformations could be provided mainly by the addition of the LK algorithm with a pyramidal representation of the CT image. However, the pyramid approach is not limited in this regard.
Thus, the use of manual alignment of bone structures and affine preregistration of CT images is of similar importance to the inclusion of the pyramid method. Without using these auxiliary steps, the results for the entire patient data tested were erroneous, as shown in Fig. 4.
Registration inaccuracies were more pronounced, particularly in the case of soft tissues. Bone structures were usually paired correctly owing to the different densities. In Fig. 6, which shows the results of the comparison by gamma analysis, it can be observed that F I G . 5. EQD2 differences between bladder DVH parameters obtained from simple DVH parameter addition and DIRbased dose accumulation.
F I G . 6. Comparison of the agreement of individual CT slice of patient no. 1 using gamma analysis after affine registration (left) and after deformable registration by LK with optimized parameters (right) (pyramid level 2, neighborhood radius of 3 voxels).
F I G . 7. Rectum (green) and bladder (blue) matching differences for patient no. 1 of study #1.
| 77 in the case of soft tissues in areas with homogeneous image intensity, the detected registration errors were higher. This is a relatively serious finding as the registration errors are subsequently transmitted to the deformation of the contours of critical structures, such as the bladder and rectum (Fig. 7).

| CONCLUSION
This study shows that the extended LK method of DIR could provide competitive accuracy of dose accumulation. The properly deformed and summed dose distribution after each BRT fraction could then be used for adaptive planning of the following fraction.
Based on the comparison with the current standard of dose accumulation, it can be concluded that the use of deformable registration should allow radiation oncologists to gain a more realistic view of the dose distribution within the patient.
The outcomes of this study also indicate some possibilities for further development of procedures for the preliminary detection of CT studies unsuitable for this type of registration; that is, finding a criterion for the magnitude of the anatomical change between fractions in which DIR fails.

CONF LICT OF I NTEREST
No conflict of interest.

DATA AVAILABILITY STATEMENT
Data available from the corresponding author upon request.