To identify potential needle trajectories passing through clustered seeds, a 3D line detection method based on the Hough transform was employed. The classical Hough transform converts a geometric entity (such as a point or line) in the Cartesian coordinate system into a parametric representation in a transformed “voting space”. Each seed position contributes a “vote” to all possible line parameters that could pass through it. Peaks in this voting space correspond to lines supported by multiple nearby seeds, representing potential needle directions.

While the traditional Hough transform operates in two dimensions for line detection in images, in this work, it was extended to three-dimensional seed coordinates using methods provided by Dalitz (28). The 3D Hough space was parameterized by a direction vector b→ describing the insertion orientation within the body and an anchor point A(x′,y′), which lies on a plane passing through the origin and orthogonal to b→. The coordinates (x′,y′) provide a two-dimensional parameterization of this plane and uniquely define the spatial position of the corresponding line. Each seed point P(x,y,z) in the Cartesian space was projected to a family of possible (x′,y′,b→) combinations in the Hough space using the forward transform (Equations 1, 2). The cumulative voting values were then used to identify dominant lines that best align multiple seed points, forming initial needle candidates.

The cumulative voting value represents the number of points that are approximately co-linear. Then, the inverse transformation (Equation 3) is used to obtain the Cartesian space representation of the anchor point. A two-point representation is adopted to facilitate constraint modeling. Specifically, the needle is parameterized by a shaft reference point P_S and a needle tip point P_T. The scalar parameter t denotes the minimal distance along the needle direction vector from the current querying point P_q to the anchor point. It is computed as in Equation 4. The needle tip point P_T is then obtained as in Equation 5. The needle shaft point P_S is defined by extending the needle by a fixed offset t_l, which corresponds to the physical length of the needle (Equation 6).

Conceptually, the Hough-based grouping replaces the manual visual alignment of seeds with an automatic “voting” mechanism, enabling consistent and reproducible identification of multi-seed trajectories before robotic adjustment.

The degree of discretization in the parameter space and the chosen partitioning strategy directly determine its ability to distinguish between lines. The computation time of the Hough transform is not significantly correlated with the number of points in the seed cloud but rather with the number of discretization units. In this study, icosahedral subdivision was used for angular discretization in 3D space. Under six levels of iteration, the angular resolution reached 1.3°, and the spatial resolution was set to 0.2 mm. These parameters provided a good balance between computation efficiency and spatial accuracy, resulting in a per-case computation time of less than 2 sec and a maximum deviation of 0.08 mm between the identified line and the actual seed positions. The deviation is sufficiently small that it does not affect the subsequent optimization steps.

The primary needle placements obtained from the Hough transform provide an initial geometric reference for potential insertion trajectories. However, these preliminary paths may still intersect anatomical obstacles or violate robotic workspace limits. To ensure clinical feasibility, an optimization-based refinement step was introduced to locally adjust each needle trajectory while maintaining alignment with its associated seed group.

The goal of this refinement is to minimize deviations between the planned seeds and their corresponding optimized needle axes, subject to both surgical and robotic constraints. The adjustment is performed within a limited spatial neighborhood to preserve the overall geometry of the Hough-derived trajectory. For each needle, the shaft point P_S is allowed to move within a cubic neighborhood (x0±r,y0±r,z0±r) around its initial position PS0(x0,y0,z0), while the tip P_T remains fixed. The optimization minimizes the sum of perpendicular distances from each seed to its corresponding needle axis, as expressed in Equations 7, 8.

This formulation ensures that each needle passes as close as possible to all assigned seed dwell positions while maintaining mechanical and anatomical feasibility. To achieve this, two categories of constraints are introduced—hard and soft constraints. Hard constraints represent non-negotiable safety or mechanical limits, while soft constraints encode desirable geometric or clinical conditions that may be relaxed slightly during optimization.

Before generating feasible needle trajectories, an initial set of planned seed positions was created using a geometry-based procedure. Random tumor shapes were generated using a multi-sphere aggregation approach. Each tumor was constructed by placing several spheres of equal radius at selected positions within the target volume. The centers of these spheres were used as the initial planned seed locations. The union of all spheres was then enveloped by a smooth outer boundary, producing continuous, irregular tumor-like geometries suitable for needle placement planning. The number and radii of the spheres were empirically selected to produce realistic tumor volumes while preserving the prescribed dose distribution. Each tumor shape contained between 20 and 30 seeds, reflecting typical clinical scenarios. The spheres’ centers of generated tumors were then used as initial seed input for the needle placement planning.

Twenty randomly shaped tumors were generated to simulate hepatic lesions. Needle path planning and subsequent evaluation were performed accordingly. Iodine-125 seeds (model 6711) with an activity of 0.5 mCi were used in the simulations.

The abdomen dataset, organ segmentation models, and needle placement modeling functions were obtained from a prior study by Yamada (29). These modeling functions were used consistently with the present implementation. Dose computation followed the TG-43 formalism, and slicerRT (30) was used for dose comparison and generation of DVH plots. Due to the finite resolution (approximately 3–4 mm in XY directions; 1.5 mm in Z direction) of the three-dimensional dose grid, extremely high dose regions in the immediate vicinity of radioactive sources were not fully resolved, resulting in a numerical truncation of the high-dose tail of the Planning Target Volume (PTV) DVH. The planning algorithm was implemented using Python, with 3D reconstruction and visualization handled using 3D Slicer (31). The experiments were performed on a system equipped with an Intel^® Core™ i7-7700HQ CPU @ 2.80 GHz processor.

A needle was considered feasible if it satisfied all geometric hard constraints and soft constraints described in Section 2.2. The seed adjustment distance was defined as the 3D Euclidean displacement between the nominal seed positions and their optimized locations. Gamma analysis was performed to compare the dose distributions before and after the needle trajectory optimization using a distance criterion of 2 mm and a dose criterion of 2% (global gamma normalization).

Phantom experiments were conducted to validate the practical feasibility and accuracy of the proposed needle placement planning method under robotic operation conditions. Our experimental platform is shown in Figure 2A, which includes a UR5 robotic arm (Universal Robots A/S, Odense, Denmark) equipped with a custom puncture end-effector and a silicone phantom simulating human anatomy, including artificial bones, organs, and tumor-like structures. Five needle insertion experiments were conducted. Each phantom was first scanned using CT (slice thickness, 0.3 mm; field of view, 144 * 124 mm) to obtain preoperative imaging data. Based on the three-dimensional reconstructions, the geometric shape of the tumor was determined. Following the inverse geometry-based initialization strategy described in Section 2.3, a sphere-based partitioning algorithm was employed to determine the initial seed positions, and the target seed dwell positions were predefined. Needle placement planning was then performed according to the methods described in Sections 2.1 and 2.2.

The UR5 robotic arm, equipped with a custom puncture end-effector capable of rotation and linear feed, executed the planned trajectories. The UR5 provides a nominal repeatability of ±0.1 mm, which is sufficient for the needle placement task considered in this study. Tool calibration and image-to-robot registration were performed following our previously published method (32–34). Specifically, the tool calibration error was within ±0.05-mm accuracy, and the image-to-robot registration achieved an average error below 0.5 mm, with high repeatability across repeated trials. After trajectory execution, seeds were manually inserted, and post-implant CT scans were acquired for evaluation.

Implanted seed positions were determined using thresholding followed by Gaussian mixture regression (GMR) to identify the seed centers. To evaluate the implantation accuracy, the spatial errors between the implanted seed positions and the planned seed locations were quantified. Axial errors were defined along the needle axis, with positive values indicating displacement toward the needle tip and negative values indicating displacement in the opposite direction. Radial errors were measured in the direction perpendicular to the needle axis.

We compared the proposed method with a conventional coplanar template approach on five cases consistent with Section 2.4, yielding a total of 10 cases. The coplanar template layout was generated manually by an experienced clinician following geometric heuristics (surface-normal insertion and parallel coplanar constraints). This manually designed baseline represents a commonly used clinical strategy and was used for comparison with the proposed automatic method.

To ensure a fair comparison, both planning approaches employed the same seed model and the same total number of seeds; only the spatial distribution of the seeds differed between the two methods. The resulting dose distributions were evaluated and compared using DVHs, dose field cross-sectional visualizations, and quantitative dosimetric metrics including the D_mean, V₁₅₀, and the conformity index (CI). The conformity index was computed using Paddick’s formulation, defined as shown in Equation 14.

where VPTV denotes the target volume, VPIV is the volume enclosed by the prescription isodose, and VPTV∩PIV represents their intersection.

As summarized in Table 1, across the 20 cases, the proposed method generated feasible needle sets with high reliability. Reliability was assessed based on seed reachability, the maximum adjustment distance applied to seeds, and the overall gamma pass rate, rather than the absolute number of needles generated per case. On average, 8.2 ± 2.1 feasible needles were produced per case, with a mean maximum seed adjustment distance of 0.6 ± 0.3 mm. The average gamma pass rate was 97.3% ± 2.7%. Not all seeds were reachable in every case: for example, Case C7 contained two missed seeds, consistent with Table 1. Figure 3 illustrates examples of optimization results under different constraints. A case was considered successful if all seeds were reachable, or if the displacement of all reachable seeds remained below 1.5 mm, resulting in an overall simulation success rate of 95%.

In addition to Gamma analysis, DVH plots were generated to compare dose distributions before and after needle trajectory optimization or a single representative case, as in Figure 4. The DVH curves showed minimal differences, indicating that the optimization process preserved the target dose distribution. Primary needle-path generation (Hough transform) required less than 1 sec per case, while trajectory optimization required an average of 43 sec per needle.

Figure 2B shows the planned needle placement for one representative case, while Figure 2C presents the actual needle insertion performed by the robot; the two results are consistent. Figure 2D shows the corresponding post-implant CT scan.

Across all insertions, the mean deviation of axial and radial errors was 0.11 mm and 0.713 mm, respectively. Axial deviations ranged from a maximum of +2.69 mm in the needle insertion direction to −2.44 mm in the opposite direction. Of the seeds, 63.3% exhibited an axial error below 1 mm, and 96.6% were within 2 mm. Radial deviations reached a maximum of 2.37 mm, with 76.6% of the seeds showing radial errors below 1 mm and 93.3% within 2 mm. The deviations of all seeds were analyzed along the axial and radial directions, and the results were summarized and visualized as violin plots in Figure 5.

Figure 6 is a representative instance, where the two planning approaches are shown in panels A and B, the isolevel plots of dose distributions are shown in panels C and D, and panel E presents the cross-sectional view of the difference in dose distributions.

The DVH result for a single representative case is shown in Figure 7, and similar trends were observed across all cases. Solid lines represent TPS planned (group TP) dose distributions from our method, while dashed lines correspond to the coplanar template manual planning approach (group MP). Blue curves indicate the PTV, and red curves indicate the organ at risk (OAR). The inset provides a zoomed-in view of the OAR DVH, with an expanded dose range along the x-axis while keeping the y-axis scale unchanged. All DVHs were normalized such that the prescribed dose was equal across plans, ensuring fair comparison. The DVH curves produced by our dose-driven needle planning demonstrate a sharper decline beyond the prescription dose.

The dosimetric comparison between the proposed and manual plans across 10 cases, as summarized in Table 2, showed that the proposed method improved dose homogeneity, with lower high-dose volumes (V₁₅₀ = 56.9% vs. 72.9%, p = 0.002; V₂₀₀ = 22.2% vs. 43.2%, p = 0.012) and reduced mean doses (D_mean = 247.6 ± 25.7 Gy vs. 320.3 ± 76.1 Gy, p = 0.019). Additionally, the proposed method achieved higher CI values (CI = 0.67 ± 0.1 vs. 0.45 ± 0.17, p = 0.011). All p-values were computed using a paired two-sided t-test.