Purpose
Interstitial brachytherapy (iBT) is a radiotherapy option, where radioactive source is guided through implanted catheters and needles. It mostly uses iridium-192 (192Ir) source that produces a steep dose gradient. This feature allows for improved sparing of surrounding tissue, resulting in at least equivalent or even superior treatment outcomes in breast cancer patients compared with analogous external beam treatment options [1-3]. Interstitial brachytherapy is a well-established treatment option, and follows a less extensive patient-specific quality assurance than external beam radiotherapy. The introduction of iBT precedes most quality assurance (QA), resulting in a lack of widespread consensus of technical solution for QA, such as cone-beam CTs or means for dose verification [3]. Instead, QA in iBT heavily relies on check lists and dual control principle, such as double-verification by two independent persons. Still, brachytherapy treatments can be plagued by treatment unpredictability, including dose delivery uncertainties as well as inter- and intra-fractional variations of implant positioning. The most prominent example of dose delivery uncertainties are source positioning uncertainties [1, 4], while inter- and intra-fractional variations can result in implant uncertainties of up to 1.0 mm [5]. However, there is no common approach for resolving the above-mentioned uncertainties. To bridge this gap, an ongoing study in Erlangen presents electromagnetic tracking (EMT) as a novel technique for quality assurance [5-15]. Based on the principle of induction, EMT allows for a 3D spatial localization of the implant by tracing each catheter or hollow needle with a millimeter-sized sensor, yielding a 3D representation of the implant. Its radiation-free application and independence of direct line of sight (in contrast to optical tracking methods, i.e., implanted catheters can still be tracked even if not visible through the skin) are some of its key features [16-19]. It is suitable for several areas of quality assurance. EMT can be implemented as pre-treatment verification, allowing for the detection of potential planning errors, therefore addressing positional uncertainties [7, 10, 12, 20-22]. Especially reconstruction errors of catheters with wrong indexer lengths and offsets as well as indexer misidentifications, such as partial swaps of neighboring catheters, can be detected through the comparison of EMT data with the treatment plan. Similarly, potential delivery errors, including shifts or swaps can be detected through a comparison of EMT measurements at different times of the treatment, or EMT data with the treatment plan [12, 20]. These comparisons can even be performed online, paving the way for an adaptive patient-specific treatment workflow. Finally, EMT data might also be employed in an automated implant reconstruction, thereby reducing the planning time [23]. Clinically relevant information can only be derived from the implant geometry acquired via EMT by co-registration to the manually reconstructed and CT-based implant geometry that serves as the basis for treatment planning. Therefore, registration is crucial for the detection of potential uncertainties, and needs to be closely examined. A sub-optimal registration may imply non-detectable treatment uncertainties or cover up noteworthy changes in implant localization. Due to steep dose gradient essential in iBT treatments, the precision of registration holds even more significance due to the fact that even slight positional deviations might lead to clinically significant dose non-conformities. Since one of the most commonly used cases of registration in iBT setting is the comparison of follow-up image data, the Therapy Physics Committee of the American Association of Physicists in Medicine commissioned Task Group 132 to review current approaches for image registration in radiotherapy [24]. Since then, only a few studies have addressed the impact of a registration method on clinical outcomes, with none of them using data obtained from iBT of the breast. Additionally, these studies focused on the registration of images, while the present study examined the registration of pointset data derived from an iBT of breast implant. The structure of the most commonly used brachytherapy applications, such as iBT of the breast or prostate needles is unique, since it consists of a cluster of elongated tubes. Therefore, it differs noticeably from typical cases of registration studies, specifically 3D imaging data including organs. These shapes typically consist of coherent surfaces, while the individual tubes of implants used in the current study were disconnected. In addition to the precision of registration, the time needed to perform a sufficiently accurate registration was analyzed. This aspect is especially relevant for an online integration of EMT, because the clinical routine often follows a very tight schedule. Therefore, integrating an additional QA measure must not prolong the treatment substantially. To conclude, the registration performance on such clinical data as well as potential shortcomings, especially in regard to registration precision and timing as well as their related workarounds were explored.
In this study, four registration algorithms were compared: Procrustes [25], AbsOr [26], rICP [27] (a point cloud-based rigid iterative closest point method), and rCPD [28] (a stochastic approach). The quality of registration in an extensive patient cohort was investigated.
Material and methods
Patient data acquisition
A comparison of registration performance with patient data was based on two sets of point clouds, which represented catheter traces of the implants of breast cancer patients receiving high-dose-rate (HDR)-iBT treatment at the Uniklinikum Erlangen (UKER), Germany. Figure 1 shows the representation of an exemplary implant, highlighting a tube-like structure of the studied implants. The first point cloud was derived from treatment plan that was based on CT imaging. It served for treatment planning, as is generated in daily clinical practice by manual reconstruction of the path for each of catheters in an implant. The second point cloud was measured using EMT. Therein, data of 69 patients who participated in an ongoing EMT patient study was evaluated (NDI-HDR-1 355-14B, 2014, renewed 2020; informed consent was obtained before inclusion). The summary of the cohort is presented in Table 1.
Fig. 1
3D representation of an implant consisting of 17 individual catheters (implanted lengths between 60 and 150 mm) in a left-sided breast, as shown by treatment planning system. A) With enabled view of the skin of patient showing left-sided breast and upper arm (cyan); B) With disabled skin but with enabled ribs (blue). Yellow tips represent the so-called “tip-end side” of implanted catheter, pink tips the connector-end of catheters, yellow dots represent the landmarks for dose normalization, and red stripes along the catheter indicate the dwell positions of radiation source
Table 1
Patients’ specification. Accelerated partial breast irradiation (APBI) treatment with up to 9 fractions; Boost – treatment with up to 3 fractions; Unique – the sum of all catheters in all patients; Total – the sum of all catheters in all patients and fractions
| Number of patients | Total | 69 |
| APBI | 57 | |
| Boost | 12 | |
| Recruitment time | Oct 2016 – Feb 2023 | |
| Number of catheters | Range | 11-27 |
| Mean | 17 | |
| Unique | 1,191 | |
| Total | 10,939 | |
Computed tomography images for treatment planning followed institutional standards, and were acquired from either SOMATOM Sensation Open (until August 2018) or SOMATOM go. Open Pro (both, Siemens Healthineers, Erlangen, Germany), with a slice thickness of 2.0 mm. Catheter traces were manually reconstructed by a medical physicist in OncentraBrachy (Elekta, Veenendaal, The Netherlands) according to the GEC-ESTRO guidelines [1, 29]. Catheter traces were exported as DICOM files of the treatment plan and analyzed using an in-house MATLAB (version R2019b; MathWorks Inc., Massachusetts, USA) system.
Electromagnetic tracking data were acquired with a previously introduced hybrid treatment delivery system, consisting of a modified prototype Flexitron (Elekta Brachytherapy, Veenendaal, The Netherlands) and electromagnetic tracking system (Aurora, NDI, Waterloo, Canada). The prototype’s check cable drive features an EMT sensor with 5 degrees of freedom, while the primary source drive only features a dummy source. Measurements were performed right after each treatment delivery fraction with clinical afterloader and CT image acquisition. Therefore, EMT sensor was successively driven through the catheters while stopping every 10 mm for 1 second. Three additional EM sensors (Aurora, NDI, Waterloo, Canada) with six degrees of freedom were attached to the skin around the patient’s breast to allow for compensation of small motions (e.g., breathing or twitching due to discomfort) [6]. During data collection, the state of check cable (moving outwards, stopping, and moving inwards) as well as the current channel number were automatically associated with 3D data points. For this work, only data associated with a stop position were analyzed, allowing for re-sampling of the treatment plan data to ensure theoretical point-to-point correspondence between EMT and treatment plan data. A more detailed introduction to this setup and the subsequent data processing can be found in Kallis et al.’s paper [8].
Registration methods
While an ideal EMT dataset would have a direct point-to-point correspondence, this direct correspondence could not be achieved for an erroneous treatment application. The most severe deviations from an idealized registration can occur in treatment delivery errors called “swaps”, i.e., a misconnection of at least two transfer tubes, and shifts, i.e., wrong indexer lengths and movement along the implantation direction. For swaps, the mismatch between the true catheter position and the misconnected position might introduce phantom rotations within the direct point-to-point registration, while shifts might lead to an additional phantom translation. To mitigate these occurrences, assuming that neither acquisition method introduces distortions, the two sets of point clouds should be rigidly registered. A widely used method for rigid registration of point clouds is the classical iterative closest point algorithm, where each point within a set is directly matched with its closest corresponding point in other set [27-29]. The definition of closest corresponding point is thereby hard-coded into the algorithm in a form of a threshold. The algorithm can focus on different measures as the threshold of correspondence, therefore giving rise to different sub-algorithms (SAs). For this analysis, four SA variations were investigated.
SA 1 (AbsOr [26]) optimizes the distance of least squares, and calculates the distance using quaternion. This optimization can be formulated as an unweighted registration problem (see Equation (1)), where the sum over all indices (Σi) of the root mean square between the known matrix B and the transformed matrix A is minimized. The applied transformation consists of a rotation matrix R and a translation vector t.
SA 2 (Procrustes [25]) opts for the Procrustes distance, also known as orthogonal Procrustes problem, and can be written as a matrix approximation problem, as shown in Equation (2). Here, the orthogonal matrix Ω, which most closely maps the two known matrices A and B based on Frobenius norm ||•||F is to be received.
SA 3 (rICP [27]) determines the Euclidean distance between the translational vectors and the angular difference of consecutive iteration steps, such as SA 1 (see Equa- tion (1)), but does not work on quaternion. Instead, a point-to-point, point-to-plane, or plane-to-plane metric can be selected. For this experiment, only point-to-point metric was considered.
SA 4 (rCPD [28]) employs the complete negative log-likelihood function Λ in a stochastic expectation-maximization algorithm as its threshold. The points of the target set B are considered as a Gaussian mixture model centroid, while the points from other set A indicate data points generated by the Gaussian mixture model, allowing for a formulation of Λ in dependence of equal isotropic covariance σ2 and transformation parameters θ.
SA 2 and SA 3 are integrated with recent MATLAB releases, therefore further information on these algorithms can be found in the MathWorks library. SA 1 follows Horn’s method, and has been implemented through a MATLAB extension [26]. To ensure comparability between all four sub-algorithms, the same thresholds were implemented throughout the study. On the one hand, the tolerance, i.e., the maximum tolerable distance between corresponding points, was set to 10-5 scale units, meaning that the tolerance was only dependent on the overall size of the objects, which have to be registered, but not on the registration method. On the other hand, the total number of registration iterations was limited to 200 iterations to prevent infinite registration loops. Since the behavior of rigid registrations was supposed to be tested, the scaling factor was set to 1 where applicable (i.e., for SA 2 and SA 4). SA 4, as the only method based on stochastic properties rather than geometric features, could be more prone to influences from extreme outliers; therefore, an additional refinement of this method was introduced (SA 4+), excluding extremely outlying catheter traces from general determination of transformation matrix, and manually adjust the outliers.
Analysis
For the first part of the study, a 10 random rotations and random translations to the data obtained from the treatment plan were performed, and the shifted data set was registered to the original, resulting in a total of 690 registrations. The translational vector varied with random lengths between 0 cm and 1 cm along principal axes, and Cardan angles for the rotations were chosen randomly between 0° and 180°. To determine how closely the registered point cloud matched with the original, the closest point in the original point cloud to each point of the registered point cloud was identified, and the Euclidean distance between these points was determined. The smaller the sum over all Euclidean distances, the better the underlying SA. The same procedure was performed using data from one EMT measurement per patient, such as the measurement performed during CT, generally regarded as the reference measurement when comparing two EMT measurements. Additionally, these registrations with additional Gaussian noise were repeated to simulate random errors due to strong interferences in the measurement. Noise was randomly added to each data point, resulting in dislocations of up to 0.1 cm.
In the second part of the study, the catheter traces obtained from EMT data were registered along with the corresponding re-sampled treatment plan. To quantify the quality of this cross-modality registration, Generalization [30] was introduced to evaluate the ability to reproduce shapes not included in registration cohort. Specifically, the registration accuracy of SAs 1-4 for each patient through a four-step algorithm based on a leave-one-out approach was compared. First, one catheter of a patient-specific implant as the target structure was declared, and its data were excluded from the subsequent registration, which was considered as the second step. In the third step, the obtained transformation to the target structure was applied before finally determining Euclidean distances between both sets of images along the target structure. This analysis was iterated over each catheter within the patient’s implant. To increase comparability, the mean Euclidean distance along each target structure was used as a measure of registration accuracy, with lower values indicating higher accuracy.
Additionally, the run-time of each registration as an additional measure of performance was measured. A lower run-time indicated better performance, assuming that the precision of registration was comparable.
To conclude this investigation, according to Kallis et al. [7], the same patient cohort was re-evaluated using a different sub-algorithm. The reference EMT measurement (typically the EMT measurement in CT suite, called EMT-CT) to the plan data of patient (plan check) was registered as well as any remaining EMT measurement to the reference EMT measurement, whilst normalizing any detected difference to the first registration. Thresholds defined by Masitho et al. [12] or a classification of potential errors were applied.
To confirm validity of the experiment, paired Wilcoxon signed-rank tests were employed to all intermediate results testing the null hypothesis that the vector x-y (where x comes from one dataset and y comes from another dataset) derives from the distribution with zero median at a 5% significance level.
Results
The results of the registration of random rotations based on the patient plan data are shown in Figure 2A. The mean distances for SA 1, SA 2, and SA 4 were lower than for SA 3, with properly registered 294 out of 690 iterations. The remaining iterations are visible as clear outliers in Figure 2A.
Fig. 2
Summary of the rotated patient data registrations (69 patients). Distances resulting from (A) plan-based data, and (B) EMT data obtained from CT measurement. The presented algorithms are (from left to right): AbsOr (SA 1, blue), Procrustes (SA 2, orange), rICP (SA 3, green), and rCDP (SA 4, violet)
The registration of rotated EMT data yielded comparable results with the same registration of plan-based data. These results are shown in Figure 2B. While the general shape of distributions was mostly the same as in a previous set of registrations, the severe outlier of SA 3 could no longer be observed, accumulating an increased success rate of 96.8% (668 out of 690 iterations) for this sub-algorithm. In Wilcoxon signed-rank tests between data of SA 1 and SA 2, and data of SA 3 and SA 4, full rejection of the null hypothesis was observed, demonstrating a minimal corresponding p-values (p = 0).
The results of the Generalization showed minor differences between SA 1, 2, and 3, while SA 4 performed worse both regarding the overall sum of Euclidian difference with the target structure (Figure 3A) as well as the calculation time (Figure 3B). A noteworthy improvement could be observed for SA 3, which performed similarly to SA 1 and SA 2, a behavior that was not observed in the previous part of the study. However, in terms of run-time performance, SA 2 surpassed all other SAs showing a tightly spread distribution around a mean registration time of t~2 = 0.002 ±0.001 s. All other SAs performed worse, with SA 4 resulting by far in the widest spread and longest registration time, with a mean run-time of t~4 = 13.22 ±4.56 s. The modification of SA 4 into SA 4+ did yield a slight improvement in terms of a decrease in deviations, but with a slight increase in run-time. The rejection of the null hypothesis for any Wilcoxon signed-rank test can be established.
Fig. 3
Results of the Generalization. Box plots from (A) the overall deviations, and (B) the associated run-time by sub-algorithms. The presented algorithms are (from left to right): AbsOr (SA 1, blue), Procrustes (SA 2, orange), rICP (SA 3, green), and rCDP (SA 4, violet)
In order to complete this research, the results of Kallis et al. [7] study re-evaluating the same patient cohort (41 patients, 355 EMT measurements) were employed, with the best-suited algorithm according to the previous results, i.e., SA 2, and these results were compared with the originally used SA 4+. Re-analyzed patient cohort of Kallis et al. [7] with SA 2 algorithm resulted in a rejection of the null hypothesis in Wilcoxon signed-rank test and a median catheter trace error of 2.22 mm, with an overview of distribution of EMT measurements, as presented in Figure 4. The overall tendencies, including lowest spread and lowest mean values of the individual boxes observed at measurements in CT suite (CT and CT2), and a slight increase in the individual mean values of boxes throughout the entire treatment cycle, are visible for either SA. The originally used SA 4+ more prominently featured larger deviations, as indicated by longer whiskers, especially for the EMT measurements of BED and F9.
Fig. 4
Box plots of the Euclidean distances of dwell positions throughout the treatment cycle for 41 patients and two different registration algorithms. Violet – SA 4+ (rCPD), as presented by Kallis et al. [7]; Orange – SA 2 (Procrustes). Dotted vertical lines indicate the separation by treatment day, while the black dashed horizontal line marks the overall median of the second data set (median of 2.22 mm)
Discussion
The Therapy Physics Committee of the American Association of Physicists in Medicine Task Group 132 focused on recommendation of validation techniques for image registration [24]. However, the landscape of studies concerning registration performance, especially in brachytherapy settings as well as in a broader radiotherapy setting, is sparse. Former research by Masitho et al. explored the influence of three different automatic registration tools on rigid registrations of MRI and CT images in brain stereotactic radiotherapy, and found a more significant difference in registration methods than in acquisition setup [31]. Therefore, their findings underline the necessity for an extended study of registration performance. A direct comparison of similar registration methods was presented in this study, but the stereotactic neuro-surgical application investigated by Cutter et al. registered ideal point clouds with rICP and rCPD algorithms, where rICP outperformed rCPD in two out of three aspects [32]. These results are in line with the current study findings, where the latter algorithm lacks particularly in terms of the run-time. In the field of prostate brachytherapy, Grayales et al. reported remaining image registration errors in the registration of TRUS and EMT, even with previously conducted rotational alignments [33]. This finding emphasizes the necessity of a reliable registration algorithm to minimize registration errors.
The current study aimed to identify the efficient rigid registration method with reproducible performances. It was noted that the previous standard method (rCPD) lacked in terms of registration time, taking up to several minutes until obtaining a reliable result. Internally, the registration performance also varied, further increasing the need for a more reliable and faster alternative. In addition, the necessity for fast registration performance is directly linked to the strive for possible online analysis just before each treatment fraction, a desirable feature of the clinical application of EMT.
In this study, we choose to compare the registration performance of four registration methods, all of which follow the iterative closest-point approach. This decision was made because the iterative closest point approach is one of the most common algorithms for rigid registration, while still providing a large pool of sub-algorithms evident for their threshold of correspondence. This choice of algorithm is more favorable than direct point-to-point registrations, since it does not discriminate against potential treatment errors, such as shifts and swaps, which might otherwise result in incorrect registrations.
The first part of the present study comprises implant structures extracted from the treatment plan and later from EMT data. The implants were translated and randomly rotated before registering the transformed data to the original implant. The third sub-algorithm still struggled to perform as compared with the other sub-algorithms; for plan-based input data, SA 4 failed a small percentage of registrations.
Marzola et al. compared two non-rigid variants of ICP and CPD algorithm, and introduced three parameters that identified ideal features of a deformable statistical shape model (SSM), which is nowadays considered the gold standard for a comparison of different SSMs [30]. The authors found that the studied ICP algorithm performed better than the CPD algorithm in two out of three parameters. In the analysis presented in this paper, the catheter traces were meant to be non-deformable; therefore, only one of the three parameters introduced by Marzola et al., i.e., Generalization of the registration, is expected to carry significant information about the performance of the registration method, which was implemented in the second part of the current study. Here, we were also able to demonstrate a superior performance of our ICP-like algorithms (SA 1 – SA 3) in comparison with the CPD-based SA 4.
The second part of the study differed from the first part by utilizing data from two different modalities as data sets for the registration. The afore-mentioned unchanged implant structures from the treatment plan served as a reference, while the EMT data was registered into it. This part of the patient study aimed at quantifying the ability to reproduce shapes not included in the registration cohort (Generalization). It holds more significance over the usability of a sub-algorithm than the previously studied general registration performance, because this measure consider the observance of the structure rather than the strict matching of the related points. In this part of the study, most SAs performed comparably, except for SA 4 (rCPD). Knowing that the general registration performance of SA 4 as observed in previous parts was comparable with other SAs, an explanation of the poor performance in this part of the study must be related to the difference in input data. We theorize that SA 4 is more strongly influenced by extreme outliers, and therefore propose a modification to this sub-algorithm, which ignores the outliers. This modification has been implemented in previous analyses, such as studies conducted by Kallis et al. [5, 7], but was omitted in the auto-registration part of the current study to allow for a direct comparison of the baseline performance of all sub-algorithms. In regard to the run-time, the performance of SA 4 was also much worse than any other sub-algorithm. This increase might be because SA 4 was the only stochastic SA, therefore requiring more calculation efforts than the geometric SAs. This observation resulted in our desire to switch to a faster and more reliable registration method for our future online analysis of patient data, thus initializing the final part of this study.
Comparing the results of the auto-registration with a technical precision of EMT verified within the sub-millimeter range [11], we found that the auto-registration was performed on the same scale, therefore showing that the registration did not yield major contributions to positional uncertainties. In the second and last part of this study, we investigated the precision of the registration whilst using data from two different sources (i.e., plan data as the target, and EMT data to be registered to the target). These registrations can be compared with a clinical application of registering of different image sets of the same patient.
Within the brachytherapy setting, anatomical variations of the catheter implant are typically less than 5 mm (median 2.3 mm, [5]). Even though dosimetric influences from variations in the breast have not yet been systematically examined, related studies on implant variations in the prostate were unable to identify a single, universally acceptable threshold, but acknowledged thresholds in the range of 1-6 mm, until clinically relevant dosimetric influences were observed [34-37]. Other technical uncertainties, such as source positioning uncertainty of the afterloader (typically, ±1 mm, [38]), also contribute to the overall number of uncertainties. In the view of this, the registration errors of 2.6 mm found in this study are negligible, and the sub-algorithms SA 1-3 can be considered sufficiently accurate for the purpose of QA in iBT.
Knowing that the registration performance of SA 4 is not ideal for online analysis, we wanted to verify previously published results from our EMT project, i.e., Kallis et al.’s paper regarding inter-fractional variations [7]. We re-analyzed the same patient cohort with the same sub-algorithm used in the original publication and different registration sub-algorithms. A slightly higher values were found for the median catheter trace error using the new sub-algorithm (2.22 mm in the current study compared with 2.19 mm in the original paper). However, in our analysis, only catheter trace errors of at least 0.8 mm were considered, hence distort the newly found median error towards larger values. The mean values for known larger shifts, represented within the whiskers of box plots in Figure 4, were overall lower with the new sub-algorithm Procrustes than for the established rCPD algorithm (4.02 ±0.31 mm for Procrustes compared with 4.18 ±0.30 mm for rCPD), indicating the potential for an improvement of registration results when optimizing the registration sub-algorithm.
The registration performance is an important feature for any sub-algorithm, and should therefore be intensely examined. Quantifying this feature through the implementation of Generalization measures and run-time measurements is reasonable, since these two features allow for deeper examination of performance rather than only focusing on spatial relation of artificially matched points. For the purposes of the study, SAs 1-3 resulted in sufficient Generalization, while SA 4 yielded worse performance and extremely long calculation times, thus discouraging from its continuous use. SA 3 showed severe uncertainties in its registration performance for pre-determined data sets in two of the three parts of this study. It is therefore not an ideal candidate for further registration tasks. The performance of SA 1 and SA 2 was indistinguishable. The re-investigation of previously published data shows minor improvements; however, the observed statements of the paper still hold.
Conclusions
The present study investigated the impact of registration method on the accessibility of EMT for integrated quality assurance in iBT. It was based on original patient data consisting of CT-based plan data and EMT measurements. It provides an extensive analysis of four different sub-algorithms. Overall, the results of the registration analysis show consistency throughout all parts of the study, and align within the expected uncertainty of clinical setting. The results highlight geometric functions (AbsOr and Procrustes) as superior to stochastic function (rCPD). Nevertheless, the re-analysis of patient data still shows consistency with previous publications, therefore highlighting the productivity of EMT as the quality assurance tool.
