An automated, fast and accurate registration method to link stranded seeds in permanent prostate implants

The geometry of a permanent prostate implant varies over time. Seeds can migrate and edema of the prostate affects the position of seeds. Seed movements directly influence dosimetry which relates to treatment quality. We present a method that tracks all individual seeds over time allowing quantification of seed movements. This linking procedure was tested on transrectal ultrasound (TRUS) and cone-beam CT (CBCT) datasets of 699 patients. These datasets were acquired intraoperatively during a dynamic implantation procedure, that combines both imaging modalities. The procedure was subdivided in four automatic linking steps. (I) The Hungarian Algorithm was applied to initially link seeds in CBCT and the corresponding TRUS datasets. (II) Strands were identified and optimized based on curvature and linefits: non optimal links were removed. (III) The positions of unlinked seeds were reviewed and were linked to incomplete strands if within curvature- and distance-thresholds. (IV) Finally, seeds close to strands were linked, also if the curvature-threshold was violated. After linking the seeds an affine transformation was applied. The procedure was repeated until the results were stable or the 6th iteration ended. All results were visually reviewed for mismatches and uncertainties. Eleven implants showed a mismatch and in 12 cases an uncertainty was identified. On average the linking procedure took 42 ms per case. This accurate and fast method has the potential to be used for other time spans, like Day 30, and other imaging modalities. It can potentially be used during a dynamic implantation procedure to faster and better evaluate the quality of the permanent prostate implant.


Introduction
The dosimetry of permanent prostate implants depends critically on the localization of 125 I seeds (Kuo et al 2014). Difficult visualization of seeds on transrectal ultrasound (TRUS) hampers the localization of seeds and consequently accurate dosimetry (Polo et al 2010). C-arm cone-beam computed tomography (CBCT) allows, like CT (De Brabandere et al 2012), for accurate localization of seeds, but visualizes the prostate boundary poorly (Smith et al 2007). Combining TRUS and CBCT improves the accuracy of the dosimetry of implants (Polo et al 2010). We developed a dynamic dosimetry procedure in which we assess the dosimetry during implantation by combining TRUS and CBCT (Westendorp et al 2007). This live method allows for immediate adaptation of the implant: if an implant shows a cold spot remedial seeds are placed in the same intraoperative procedure.
The geometry of an implant is not stable over time (McLaughlin et al 2006, Moerland et al 2009, Steggerda et al 2015. Edema of the prostate pushes seeds away from the center (Westendorp et al 2012) and seeds migrate after implantation. Better understanding of geometry changes may help to improve the implantation practice and could assist in a modified strand and seed design.
Migration of seeds can only be quantified if each seed can be tracked over time. In each image-dataset all seed-positions have to be identified and related to positions in a reference dataset. It is impossible to perform this task manually. An automated model is needed to connect the seed positions for each image dataset, after which the result can be visually inspected.
We developed a linking procedure to track movements of individual seeds and validated it on 699 intraoperative TRUS and CBCT datasets. With adjustments, the linking procedure may also be used to study implant geometry changes for longer time spans, like the time from implantation to Day 14 or Day 30.

TRUS and CBCT datasets
Intraoperatively obtained TRUS and CBCT datasets of 740 consecutive patients, treated in the period from October 2007 to May 2012, were reviewed. Forty one cases lacked an image dataset due to technical reasons, leaving 699 cases eligible for analysis. All patients were treated according to our protocol that makes use of combined TRUS-CBCT imaging (Falcon 2101 EX and Flex Focus 400, BK Medical; Herlev, Denmark and Siemens Arcadis Orbic 3D; Siemens Medical Systems, Erlangen, Germany) (Westendorp et al 2007). Stranded seeds of ∅ × 1 5 mm were used (IBt 1251L; Seneffe, Belgium, IBt-Bebig I25.S06; Berlin, Germany, Bard STM1251; Murray Hill, NJ USA). Spacing was equidistant in most strands, but also strands with varying spacing were used. Seed positions were initially localized on the TRUS datasets directly after releasing the strands from the needles. The seed finder of the treatment planning system (Variseed 7.2-8.02; Varian Medical Systems, Inc., Palo Alto, CA USA) was used to identify seeds on CBCT, after which the seed positions were visually inspected. These CBCT seed positions were used for clinical dosimetry.

Linking procedure
DICOM RT objects were exported from the treatment planning system. The DICOM RT Plan objects provided the seed positions as identified on the treatment planning system for both TRUS and CBCT.
TRUS seed positions were acquired first and served as the reference-distribution (X). The CBCT distribution was the compared distribution (Y, figure 1). The center of mass (COM) of both seed distributions was set to 0, 0, 0 ( ). After that, four linking steps (I, II, II and IV) were performed. The linking procedure is visualized in figure 2 and depicted as a flowchart in figure 3.
Step I. The initial distance between both distributions was defined as the minimum total squared distance, similar to the approach described by Chng et al (2011) and Jain et al (2005).
First, a cost matrix was generated with all squared distances from the reference, TRUS, distribution X (with seeds Occasionally, due to human error, the number of seeds in both distributions was not equal. In some cases a seed was not identified on one of the image datasets. In other cases a misidentification led to too many seeds in a distribution. In case ≠ n m, the distribution with the lowest number of seeds determined the number of links. After the cost matrix (c ij ) was defined, the total cost σ was determined: where 1 if th seed from is assigned to th seed from 0 if not, under the restriction that each seed is used once.
The Hungarian algorithm (Kuhn 1955) was used to solve the linear assignment problem, resulting in the minimum value for σ. The resulting links represented the lowest cost but Note Phys. Med. Biol. 60 (2015) N391 N394 could be physically incorrect, a mismatch (figure 2(a)). In the next steps strand properties were used to correct the mismatches originating from step I.
Step II. In this step strand properties of Y were determined and, if links were physically unlikely, i.e. thresholds are exceeded, links were removed. This is visualized in detail in the right part of figure 3 and in figure 2(b). Based on the TRUS DICOM RT Plan object (seed coordinates) all strands were identified in X. The reference strands contained the seeds (x i ) and provided links to y j . The links to y j provided an initial estimate of strands in Y. Following, the curvature (κ) (Stewart 2015, p 877-82) where r is the radius of curvature. The coordinates of a strand can be parameterized as γ.
Under the assumption that γ is a function of z (γ = f z ( )), where z is the needle insertion direction, the curvature is defined as: In case the angle of the strand with respect to the z-axis is small, nator is approximately equal to unity and κ can be approximated by Note that the second derivative is a measure of the position-variation of seeds around a straight line. The total curvature (K) of a strand (with N seeds) in the discrete situation (with = y x y z , , Only strands with 3 ⩾ seeds could be taken into account, as the discrete second derivative can only be computed if N 3 ⩾ . Seeds close together in the z-direction can increase the curvature excessively. Therefore, a lower limit of 10 mm was applied for ∆z. This limit equals the minimum physical distance between the upper and lower neighbors of a seed in a strand. Linking procedure for two strands. In this example the compared distribution (•, Y ) contains one seed more than the reference ( X , ) distribution. Dashed lines: connections of seeds in strand. Thin solid lines: links between reference and compared seeds.
Step I: seeds are initially linked.
Step II: a linefit of the (compared) strand is added ( lf, thick gray lines) and links causing violation of strand constraints (dotted lines) are removed.
Step III: new links, not violating strand constraints, are generated.
Step II.a If the curvature of a strand in Y exceeded a curvature-threshold (see section 2.3), the link that accounted for most curvature was removed from the links of that strand. After removal of a link, the remaining links of the strand were re-evaluated using equations (1) and (2) and solving this system with the Hungarian algorithm, for the involved x i and y j (figure 3). II.b A line was fitted through all y j of a strand by least square regression using singular value decomposition. The linefit (lf) is visualized in figure 1. The linefit extended 10 mm beyond the first and last seed in the strand. II.c Also shown in figure 1 is the distance of a seed in Y to the linefit (dlf, for calculation see appendix A). Furthermore the link-length (ll) was calculated (figure 1). All links with ll and dlf exceeding the maximum distance-threshold (see section 2.3) were excluded from the strand. Next, all loose seeds were unlinked. Consequently, at the end of step II less links between X and Y were present than at the end of step I.
Step III Seeds that were unlinked in step II were, if possible, linked anew, taking into account the restrictions of the thresholds. First, a cost matrix (equation (1)) was generated with dlf 2 of non-linked seeds to strands missing a link. Seeds without dlf's within the distance-threshold were discarded. The Hungarian Algorithm was used to solve the cost matrix. If the curvature as well as ll or dlf were within the curvature-and distancethresholds, y j was linked to the corresponding x i . Finally, free y j closest to loose x i were linked based on proximity if ll was within the distance-threshold. At the end of step III, only a small fraction of seeds was not linked.
Step IV The linefits that originated from step II were taken as starting point again. The dlf's of the remaining non-linked y j to linefits of incomplete strands were determined. Seeds of which the dlf was below the distance-threshold were linked. There was no threshold for the curvature in this step.
After finishing the linking steps, Y was transformed with an affine transformation. The transformation was determined using a linear least squares regression of ∑ ll ij 2 with singular value decomposition. The transformed compared distribution (Y k+1 ) forms the starting point for the next iteration. The transformation took place in the following order: From the second iteration on, the links were compared to the preceding iteration. If no differences were observed in the links the linking procedure was ended. A maximum of 6 iterations was allowed.

Parameters
In the linking procedure two parameters had considerable effect on the resulting links: First of all, the distance-threshold for the dlf and ll determined if links were allowed. Secondly the curvature-threshold strongly affected the outcome. Increasing the thresholds led to more links but also allowed more mismatches. Initial parameters were set by visual inspection of all 699 linked distributions. Inspection was performed using a rotatable 3D graph similar to the graphs displayed in figure 4. Cases in which a mismatch was present were identified. A subset of 24 cases with mismatches was created. The parameters were adjusted to minimize the number of mismatches in the subset. Subsequently the new settings were validated on the remaining cases by visual assessment. This process was repeated until the results were not further improved. The tuning of the parameters lead to values of 8 mm for the distance-threshold and 0.0125 mm −1 for the curvature-threshold. The visual inspection was performed in approximately 10 s per implant.

Analysis
The results of the linking procedure were compared to that of the initial linking only ( iteration 1, step I). All implants were visually inspected by one observer (HW). If the observer could improve the automatically generated links, the match was classified as mismatch; if the observer considered the result as possible mismatch, but he could not improve it, the match was classified as uncertain. For each implant, the mismatched and uncertain links, and involved strands were registered.

Computation
The linking procedure was implemented in Python code (Rossum et al 1998). All DICOM RT Plan objects were parsed to a SQLite database. Most scientific calculations were carried out with Numpy and Scipy (Oliphant 2007). For specific parts of the code, computation speed was greatly improved by using Numba (Continuum Analytics, Austin, TX, USA), a just-in-time compiler. Visualizations were realized with the Matplotlib plotting library (Hunter 2007). Analysis was performed using IPython (Perez and Granger 2007)

Results
An example of the linking procedure is depicted in figure 4. Step I shows the initial linking of seeds, resulting in multiple mismatches in this case. In step II strand constraints are introduced and links causing curvature-or distance-thresholds violation are excluded. The next step, III, shows the re-linking of seeds that were unlinked in step II, restricted by curvature-and distance-thresholds. In step IV, seeds close to linefits are linked, regardless of the curvature of the resulting strand. Than, the ll of all seeds was minimized by transforming Y with an affine transformation: T 1 [ ]. Next the linking procedure (step I-IV) was restarted on the transformed distribution. Notice that the transformation improved the initial results of step I. In this example step IV of iteration 2 shows that one seed was not linked: the distance (dlf ) between the residual seed and the strand missing a seed exceeded the distance-threshold.
Comparing the final result of the linking procedure with the results after the initial linking in step I, iteration 1, 233 out of 699 (33.3%) cases showed differences. This illustrates the added value of introducing the strand properties. The number of links after step I is by definition equal to the maximum number possible (100%). In step II the linked percentage decreases (∼97%), whereas in step III (∼99%) and IV (>99.8%) most seeds are linked again leaving a small fraction of links not feasible due to the thresholds.
Of all implants the observer classified 11 (1.6%) to contain a mismatch, 12 (1.7%) as uncertain. On the strand level 13 of 12860 (0.10%) strands were mismatched, 25 (0.19%) uncertain. Of 49897 seeds 17 (0.026%) were mismatched and 28 (0.050%) had an uncertain link. The results of the linking procedure, subdivided according to their final iteration, are shown in table 1. In this table the number of implants (second column) that reached their final linking are shown per iteration (row). The last three columns show the high percentage of linked seeds, the mean ll and the mean dlf of the registered implants, respectively. The mean ll after finishing the linking procedure was ± 3.4 0.8 mm (1 SD). The average dlf equaled ± 1.6 0.7 mm (1 SD). Of all seeds more than 99.8% was linked after the final iteration. In figure 5 the mean link-length per implant is visualized after each iteration. In iteration 1 and 2 all implants were incorporated, this number was much smaller in the last four iterations, shown in the top of figure 5. Notice, the reduction in mean link-length values from iteration 1 to 2, resulting from the transformation (equation (7)). On average, the linking procedure was performed in 42 ms per case.

Discussion
We present an automated method that allows tracking of seeds between intraoperatively acquired TRUS and CBCT. Our linking procedure improved in 33% of the cases the initial linking (iteration 1, step I) that is based on minimization of the cost matrix using the Hungarian Algorithm. Figure 4 exemplifies that the linking procedure made changes to the initial linking of step I. Especially the transformation applied after the first iteration improved the results as illustrated in figure 5. The linking procedure ended after the second (560 cases, 80%) or third iteration (124 cases, 18%) for most implants (table 1), which shows that the linking procedure was rapidly converging.
One observer visually inspected all results. Currently a metric that takes into account all considerations of a human observer, for example the classification of links as uncertain or mismatch, is not available. The linking procedure could be extended with a method that only presents cases with a risk of mismatches. Note that it is impossible to link seeds manually if only seed coordinates are presented. The visual inspection revealed only a small percentage of implants with uncertainties (1.7%) or mismatches (1.6%).
The linking method we present partly resembles the work presented by Chng et al (2011). Both methods use the properties of strands and make use of the Hungarian Algorithm to make an initial estimate of the links. There are however important differences between both methods. The work of Chng et al, S-reconstruction, was developed to compare seed positions of pretreatment plans to realized seed positions, localized on Day 0 CT. We primarily developed the linking procedure to compare intraoperatively acquired TRUS-and CBCT-based seed positions. The number of parameters that needs to be set was two in the current study instead of 12 in case of the S-reconstruction. Contrary to our method S-reconstruction relies on a fixed interseed distance. Our study showed 1.6% mismatches and 1.7% uncertainty while S-reconstruction showed mismatches in 26% of the cases. This might be caused by the fact that other types of data were analyzed. Finally, the linking procedure was performed in 42 ms per case, whereas the S-reconstruction took 5-10 s to complete. The number of implants (699) used to design and validate the procedure in the present work is approximately ten times larger than that for the design of S-reconstruction.
In future work we will use the linking procedure to analyze differences between intraoperatively obtained TRUS and CBCT seed positions. In that work we will be able to quantify the early migration of seeds based on two intraoperatively acquired image datasets. Currently the cause of migration is explained in several ways. McLaughlin et al (2006) suggested that anchoring of strands in the levator ani may explain migration. Steggerda et al (2015) hypothesized that strands are pushed out during voiding. The currently presented linking procedure may help to characterize the migration patterns and corresponding time trends and provide a better understanding.
Potentially the linking procedure can be extended to allow for incorporation in an intraoperative dynamic procedure (Polo et al 2010) allowing high accuracy registration and dosimetry prediction. As the presented procedure does not restrict seed spacing, the concept of strands might be extended to implantation channels, making it possible to study migration of loose seeds as well. Seed positions extracted from other data than TRUS or (CB)CT, like MRI, can be processed with the presented linking procedure. Other time spans could possibly be studied with adjustments of the threshold-parameters. For example migration from Day 0 to Day 30 could be investigated using the presented linking procedure.

Conclusion
A high accuracy, high speed seed linking procedure for intraoperatively acquired TRUS and CBCT datasets is presented and validated. This method has the potential to be generalized to other time spans and imaging modalities.