A method, device, medium, and program product for image feature point separation based on intensity terrain model and 3D image reconstruction based on applicator.

Abstract

This specification provides an image feature point separation method, device, medium, and program product based on an intensity terrain model and an applicator-based 3D image reconstruction, relating to the field of intelligent healthcare. The separation method includes: acquiring a target connected region of the image to be processed, the target connected region including multiple bright feature points adhered due to imaging artifacts or spatial adjacency of objects; mapping the two-dimensional coordinates and grayscale intensity values ​​of each pixel within the target connected region to a 3D intensity terrain model; clustering the pixels in the 3D intensity terrain model based on a set dynamic grayscale height threshold and grayscale connectivity constraints; for any two pixels, if there exists at least one connected path and the grayscale intensity values ​​of all pixels on that path are greater than or equal to the dynamic grayscale height threshold, then the two pixels are assigned to the same sub-cluster; extracting the pixel with the largest grayscale intensity value in each sub-cluster as an independent feature point representing that sub-cluster.

Description

Technical Field

[0001] This invention relates to the field of intelligent healthcare, and more specifically, to a method, device, medium, and program product for image feature point separation based on an intensity terrain model and three-dimensional image reconstruction based on an applicator. Background Technology

[0002] Interstitial brachytherapy is an important treatment for tumors in gynecology, head and neck, and other areas. In this treatment, multiple metal implantation needles (appenders) are precisely implanted into the tumor target area. Then, guided by medical imaging (such as CT), the three-dimensional spatial channel of each implantation needle is reconstructed in the treatment planning system (i.e., appender reconstruction) to plan the position and timing of the radiation source. This allows for increasing the radiation dose to the tumor target area while minimizing the radiation dose to surrounding normal tissues, thereby improving the local tumor control rate and reducing treatment-related toxicities.

[0003] For traditional intracavitary brachytherapy, reconstruction is relatively easy due to the fixed geometry of the applicator, its relatively stable position within the body's natural cavities, fewer imaging artifacts, and the fact that the catheter structure is usually clearly identifiable in images. However, in interstitial brachytherapy, the metal implantation needles, after penetrating the tumor tissue, easily produce obvious metallic streaks in CT images, affecting the accuracy of visual judgment of radiation source localization. Furthermore, due to the thin diameter and thin walls of the implantation needles, they are prone to local bending after removal of the core under the influence of muscle tissue, gauze packing expansion, or external forces. For larger tumors, due to anatomical constraints, multiple implantation needles often need to be arranged in a crisscross pattern to achieve ideal dose distribution, further increasing the difficulty of reconstruction. Under these complex conditions, the process is susceptible to subjective factors such as operator experience and visual fatigue, making it difficult to guarantee the consistency of reconstruction results. Therefore, there is an urgent need to develop an automatic reconstruction algorithm that can adapt to the above complex situations to improve the accuracy of implantation needle reconstruction and clinical efficiency.

[0004] Currently, researchers have introduced physical sensing-based positioning technologies, such as electromagnetic tracking (EMT), for the digital reconstruction of applicators. This technology integrates miniature electromagnetic sensors into the applicator catheter or guide core, acquiring the sensor's three-dimensional position and orientation information in real time within a known magnetic field environment constructed by an electromagnetic field generator. This information is then mapped to the patient's image coordinate system, thereby achieving automatic reconstruction of the applicator's true spatial path. The drawbacks include the need for additional hardware and disposable or reusable sensors, resulting in higher system costs. It is also sensitive to metallic environments and electromagnetic interference, and increases the complexity of clinical procedures.

[0005] In addition, some scholars have attempted to use deep learning methods to achieve automatic digitization and reconstruction of applicators for high-dose-rate brachytherapy (HDR-BT). Such methods typically model the identification and centerline extraction of the applicator in CT (or MRI) images as an automatic segmentation or regression problem. The typical process includes preprocessing the image data (such as normalization, cropping, and data augmentation), using convolutional neural networks (commonly U-Net and its variants) to automatically segment the applicator region to obtain a three-dimensional voxel mask, and then generating a continuous source channel geometric model through skeletonization, centerline extraction, or direct regression of key path points, combined with curve fitting. This type of method has a high degree of automation, strong adaptability to complex geometric structures and multi-catheter scenarios, and can achieve millimeter-level or even sub-millimeter-level reconstruction accuracy under good image quality conditions. However, its limitations lie in its strong reliance on large-scale, high-quality labeled data, and the model's generalization ability may be affected by changes in imaging equipment, scanning parameters, and applicator type. At the same time, errors in the output of deep learning models are not easily identified intuitively, lack clear physical constraints and interpretability, and still require strict quality control and manual review in clinical applications. Summary of the Invention

[0006] This invention aims to at least solve one of the technical problems existing in the prior art. To this end, this invention provides a method, apparatus, medium, and program product for image feature point separation based on an intensity terrain model and 3D image reconstruction based on an applicator. Compared to traditional automatic reconstruction methods based on thresholds or region growing, this application does not rely on the assumption of donor boundary continuity, effectively addressing issues such as boundary blurring or missing image signals caused by adhesion of multiple metal implants and severe artifacts. It significantly reduces the risk of missegmentation and channel mismatch, especially in clinical scenarios with small implant spacing and frequent crossings. The introduction of RANSAC significantly enhances robustness to outliers, automatically removing abnormal non-source path centroids caused by factors such as metal artifacts, thereby improving the stability and geometric consistency of channel fitting.

[0007] Compared to template library registration methods, this method does not rely on specific applicator geometry models, making it suitable for non-standardized or personalized catheter insertion and offering better versatility. Furthermore, compared to deep learning methods, this method does not require large-scale labeled data, and its algorithm flow has clear physical and geometric interpretability, facilitating clinical quality control and manual review. It is more feasible and safer in real-world clinical environments with limited sample size or poor image quality.

[0008] The first aspect of this application discloses an image feature point separation method based on an intensity terrain model, comprising: S101, Obtain the target connected region of the image to be processed, the target connected region includes multiple bright feature points that are stuck together due to imaging artifacts or object spatial adjacency; S102, map the two-dimensional coordinates and grayscale intensity values ​​of each pixel in the target connected region into a three-dimensional intensity terrain model, wherein the grayscale intensity value is the elevation; S103, based on the set dynamic grayscale height threshold, the pixels in the three-dimensional intensity terrain model are clustered according to the grayscale connectivity constraint. For any two pixels, if there is at least one connected path and the grayscale intensity value of all pixels on the path is greater than or equal to the dynamic grayscale height threshold, then the two pixels are assigned to the same sub-cluster. S104. Extract the pixel with the largest gray intensity value in each sub-cluster as the independent feature point representing that sub-cluster, and obtain the candidate set of independent feature points.

[0009] In some embodiments, the method further includes: If the number of independent feature point candidate sets is less than a preset minimum value and the geometric size of the target connected region is greater than a first size threshold, it is determined to be over-merging and a first supplementary strategy is executed to increase the number of feature points; if the number of independent feature point candidate sets is equal to a preset intermediate value and the combined size of the target connected region is greater than a second size threshold, it is determined that there is a missing peak and a second supplementary strategy is executed to increase the number of feature points. Output the set of independent feature point coordinates after the supplementary strategy is executed.

[0010] The second aspect of this application discloses a method for reconstructing a three-dimensional image of an applicator, comprising: S201, acquire CT medical images containing metal implanted needles, perform standardized preprocessing on the medical images to eliminate equipment differences and enhance tissue contrast; S202, perform pixel-level analysis on the preprocessed single-frame image, identify initial candidate high-brightness points, and use the initial candidate high-brightness points as seed points to construct corresponding connected regions using the region growing method to obtain multiple connected regions; S203, based on the geometric dimensions of the connected region, the connected region is divided into a first type of region for characterizing the cross-section of a single implantation needle, and a second type of region for characterizing the cross-sections of multiple implantation needles that are stuck together due to artifacts or adjacency. S204, for each of the second type of regions, the independent needle point coordinates of multiple implantation needle cross sections are separated and extracted using the method disclosed in the first aspect of this application; S205, integrate the centroid coordinates of the first type of region and the coordinates of all independent needle points separated from the second type of region to obtain a global candidate needle point set; S206, an improved random sampling consensus algorithm is used to perform spatial path clustering and fitting on the global candidate needle point set to obtain the spatial trajectory of the applicator.

[0011] A third aspect of this application discloses a computer device, comprising: a memory and a processor; the memory is used to store a computer program; the processor executes the computer program to implement the steps of the above-described method.

[0012] The fourth aspect of this application discloses a computer-readable storage medium having a computer program stored thereon, which, when executed by a processor, implements the steps of the above-described method.

[0013] The fifth aspect of this application discloses a computer program product, including a computer program that, when executed by a processor, implements the steps of the above-described method.

[0014] The complete technical path of this application: First, artifact suppression and structural enhancement are performed on the original CT images to reduce the interference of metal artifacts on the identification of high-density structures, while preserving the main geometric features of the metal implant needle. Then, candidate high-density voxel sets are extracted through thresholding and connected component analysis, and unsupervised clustering is used to group spatially overlapping or adherent voxels. Even when precise segmentation of individual catheter boundaries is not possible, the possible centroid locations of each independent applicator channel on the cross-section are identified. This strategy transforms the complex problem of relying on precise boundary segmentation into a centroid estimation problem based on statistical features, thus significantly improving stability in complex artifact environments. Finally, the axially ordered set of centroids is input into the RANSAC algorithm to fit a straight-line or curve model to the spatial orientation of the applicator channel, achieving robust reconstruction of the applicator path.

[0015] This application has the following beneficial effects: 1. This application addresses the issue that in real-world scenarios where needle points are very close or adjacent, complex metallic artifacts can occur due to the interplay of various factors such as beam hardening, photon starvation, and scattering effects caused by limitations in CT hardware, resulting in image adhesion. Traditional region growing algorithms struggle to effectively separate these artifacts. To address this problem, this application proposes a "multi-stage clustering and point separation strategy based on an intensity terrain model." This strategy treats pixel grayscale as terrain elevation and utilizes a dynamic "horizontal plane" (dynamic grayscale height threshold) for watershed-style clustering. This separates multiple needle points that are physically independent but connected at grayscale valleys due to artifacts into different sub-clusters. Furthermore, through vertex correction, it automatically removes spurious extrema caused by noise and fills in any missed needle points due to threshold settings, outputting 3D coordinates of independent needle points that more closely approximate their true coordinates. By constructing and segmenting grayscale intensity terrain, multiple implanted needle sections that are physically independent but connected at grayscale valleys due to metallic artifacts are separated into different sub-clusters, thus overcoming the difficulty of separation in such adhesion scenarios using traditional thresholding or region growing methods.

[0016] 2. In the needle point detection stage, the proposed solution employs a multi-stage strategy to achieve multi-step localization from coarse to fine: Preliminary screening and region growing: Candidate regions are quickly identified through neighborhood difference analysis, and then pixels are aggregated through region growing to efficiently distinguish between single-needle points and multi-needle adhesion regions.

[0017] Targeted adhesion handling: For adhesion areas, the above-mentioned intensity terrain model clustering is adopted, supplemented by a "multi-stage vertex correction" mechanism (such as supplementing the missing third high point), which can adaptively deal with the problem of missing needle points or excessive merging caused by changes in artifact intensity or local gray-scale anomalies, ensuring that the feature point set extracted in a single CT slice has higher integrity and accuracy.

[0018] 3. This application applies an improved random sampling consensus algorithm to 3D point cloud fitting, enhancing the geometric consistency and stability of 3D spatial path fitting. First, through multiple iterations and interior point screening, it can automatically and effectively remove abnormal discrete points caused by metal artifacts or residual points from the 2D detection stage, ensuring the purity of the fitted base data. Second, the algorithm specifically includes an optimization step for repeated points under the same Z coordinate, which can further solve the problem of extremely close adhesion points that were not completely separated in the 2D level at the 3D level, achieving a clear division of point cloud affixes. Finally, through local line segment fitting (based on multiple 2D plane projections) and a secondary allocation strategy for remaining points, it can better approximate the slight bending morphology that the implanted needle may have in the tissue, and the fitting results are more consistent with the actual clinical situation. Attached Figure Description

[0019] To more clearly illustrate the technical solutions in the embodiments of the present invention, the accompanying drawings used in the description of the embodiments will be briefly introduced below. Obviously, the accompanying drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0020] Figure 1 This is a schematic diagram of the method flow provided in the first aspect of the present invention; Figure 2 This is a schematic diagram of the method flow provided in the second aspect of the present invention; Figure 3 This is a flowchart of CT image stripe artifact correction provided in an embodiment of the present invention; Figure 4 This is a schematic diagram of needle point detection based on neighborhood difference provided in an embodiment of the present invention, wherein, Figure 4 A is a schematic diagram of the 16-neighborhood of pixel P. Figure 4 B is a schematic diagram of the bounding box of the target region in the pelvic CT image. Figure 4 C is a schematic diagram of the detection box for a multi-target region in a pelvic CT image; Figure 5 This is a schematic diagram of the target detection grayscale image and candidate box provided in an embodiment of the present invention; Figure 6 This is a schematic diagram of the intensity terrain model provided in an embodiment of the present invention; Figure 7 This is a schematic diagram illustrating the limitations of the initial point cloud provided in this embodiment of the invention compared to traditional methods; wherein, Figure 7 A represents the original three-dimensional point cloud distribution of the implanted needles. Figure 7 B represents the traditional overall linear fitting effect; Figure 8 This is a schematic diagram illustrating the outlier preprocessing effect and the fitting and clustering results of the improved RANSAC algorithm provided in this embodiment of the invention; wherein, Figure 8 A represents the point cloud distribution after outlier removal. Figure 8 B represents the clustering and fitting results of the three-dimensional paths of the eight implantation needles based on the improved RANSAC algorithm; Figure 9 This is an iterative convergence evaluation curve provided in an embodiment of the present invention; wherein, Figure 9 A represents the curve showing the change in the average distance of the interior point set during the iteration process. Figure 9 B is the curve showing the change in the length of the interior point set during the iteration process; Figure 10 This is an in-point distance distribution diagram in the fitting accuracy analysis provided in this embodiment of the invention; Figure 11These are comparison images before and after window width and window level processing, provided in an embodiment of the present invention. Figure 11 A represents the state before processing. Figure 11 B represents the processed version; Figure 12 This is a result image of manual annotation provided in an embodiment of the present invention. The right side is an enlarged view of the left center annotation area. The result of manual annotation is represented by red dots, and the annotations 1-8 accurately mark the center or prominent position. Figure 13 This is a result image of automatic annotation using the method of this embodiment provided by the present invention, wherein the right side is an enlarged view of the left central annotation area; Detailed Implementation To enable those skilled in the art to better understand the present invention, the technical solutions in the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings.

[0021] In some processes described in the specification, claims, and accompanying drawings of this invention, multiple operations appearing in a specific order are included. However, it should be clearly understood that these operations may not be executed in the order they appear in this application, or may be executed in parallel. The sequence numbers of the operations, such as 101, 102, etc., are merely used to distinguish different operations and do not represent any execution order. Furthermore, these processes may include more or fewer operations, and these operations may be executed sequentially or in parallel. It should be noted that the descriptions such as "first," "second," etc., in this application are used to distinguish different messages, devices, modules, etc., and do not represent a sequential order, nor do they limit "first" and "second" to different types.

[0022] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0023] To more clearly describe the specific process and results of the automatic reconstruction method in this application, which combines artifact suppression, centroid clustering, and random sample consensus (RANSAC) fitting, compared to traditional automatic reconstruction methods based on thresholds or region growing, this embodiment uses hospital clinical data from high-dose-rate brachytherapy with Ir-192 as an example to illustrate this application. However, this description is only one use case and is not intended to limit the scope of protection of this data processing method. The specific exemplary description is as follows: Figure 2This is a schematic flowchart of a three-dimensional image reconstruction method for an applicator provided by an embodiment of the present invention. Specifically, the method includes the following steps: 1. Database This application utilizes hospital clinical data from high-dose-rate brachytherapy using Ir-192. Interstitial needles from Elekta Medical, made of stainless steel with an outer diameter of 1.5 mm, were used for implantation. Implantation was performed precisely under CT image guidance after anesthesia. The scrub nurse recorded the length and number of implanted needles. Pelvic CT images were acquired using Siemens Caredose 4D technology and imported into the Oncentra Brachy treatment planning system. An experienced gynecologic oncologist delineated the target area and critical organs, including the rectum, bladder, and sigmoid colon. An experienced medical physicist used the TPS applicator reconstruction module to mark the tip and end of each needle on the 3D CT images and measure their length for needle reconstruction. The treatment plan aimed to ensure that D90 ≥ 100% of the prescribed dose. The specific prescribed dose depends on the treatment target. The DVH parameters for brachytherapy mainly include D2cc, D1cc, and D0.1cc, representing the radiation dose received by the highest 2cc, 1cc, and 0.1cc volumes, respectively. The dose distribution of the entire radiotherapy regimen is calculated using the EQD2 formula.

[0024] In some embodiments, the terms "subject," "test subject," or "sample to be tested" as used herein refer to any animal (e.g., a mammal), including but not limited to humans, non-human primates, rodents, etc., which will become the recipient of a particular treatment. Generally, the terms "subject" and "patient" are used interchangeably in this application when referring to human subjects. Preferably, the subject is a human.

[0025] In some embodiments, the applicator refers to an implantation needle, which is a hollow container or catheter inserted into a natural cavity or between tumor tissues before radiotherapy to provide a channel and positioning support for the subsequent delivery of the radiation source. The radiation source is delivered into the applicator to deliver a high dose of radiation to the tumor, thus precisely delivering the radiation source to the tumor site.

[0026] 2. Dataset Processing To address the complex challenges of imaging metal implants in CT images, such as stripe artifacts, noise interference, and multi-needle adhesion, this application proposes a phased, layer-by-layer data processing and needle point detection framework. This method progressively narrows the search space and enhances detection stability through image preprocessing, candidate needle point detection, region growing, clustering, and dynamic correction. Its flowchart is shown below. Figure 3 As shown.

[0027] S201: Acquire CT medical images containing metal implanted needles, perform standardized preprocessing on the medical images to eliminate equipment differences and enhance tissue contrast; In some specific embodiments, during the image preprocessing stage of CT imaging, the original pixel values... The original pixel values ​​lack a unified physical meaning, and their values ​​are significantly affected by scanning equipment, reconstruction algorithms, and imaging parameters. To eliminate grayscale differences caused by different scanning conditions and ensure the comparability and physical interpretability of subsequent grayscale-based region growing and clustering operations, this application linearly converts the original pixel values ​​into a standardized Hounsfield Unit (HU) representation based on the RescaleSlope and RescaleIntercept parameters in the DICOM header file. First, CT image data is read, and its two-dimensional pixel matrix Iraw(x,y) and related metadata information (including pixel spacing, slice thickness, RescaleSlope and RescaleIntercept parameters, etc.) are extracted. The original pixel values ​​of CT images have significant characteristics: First, there is no unified standard for original pixel values ​​under different equipment and scanning protocols, and the numerical range varies greatly (usually 16-bit integers), making direct cross-equipment quantitative analysis impossible; second, the original pixel values ​​only reflect the detector signal intensity and lack interpretable physical meaning in clinical diagnosis, failing to directly correspond to the density characteristics of human tissue. To solve the above problem, the original pixel values ​​need to be converted into standardized Hounsfield units (HU)—this is a universal physical quantity for characterizing tissue density in CT images. The conversion formula is shown in Formula 1 below:

[0028] Here, S represents RescaleSlope and B represents RescaleIntercept, both of which are determined by the calibration parameters during device scanning and can be extracted from DICOM metadata.

[0029] After HU conversion, although pixel values ​​have physical meaning, the HU value range for different human tissues (such as bones, soft tissues, and needles) can still be thousands. Direct visualization will compress most details, making it impossible to clearly distinguish the grayscale differences between the needle point and surrounding tissues. To enhance the contrast of the target area (needle point and surrounding tissues), the image needs to be transformed by window width (W) and window level (L). This transformation maps the range of HU values ​​of interest to the [0,1] interval, highlighting the grayscale levels of key tissues. The mapping formula is shown in Formula 2 below:

[0030] Furthermore, to address the unavoidable stripe artifacts in CT imaging of metal implanted needles, morphological operations (iterative operations of dilation and erosion) are employed for correction to reduce artifact interference while maintaining the high brightness features and edge smoothness of the needle point, thereby providing a more stable input for subsequent detection.

[0031] S202: Perform pixel-level analysis on the preprocessed single-frame image to identify initial candidate high-brightness points. Using the initial candidate high-brightness points as the gray-scale similarity criterion, construct corresponding connected regions using the region growing method to obtain multiple connected regions. In some specific embodiments, the method for identifying initial candidate high-brightness points includes: needle point detection based on its high-brightness performance in CT images. Since the starting end of the implanted needle is a solid, integrated structure with a sharp, conical shape, its physical size is significantly smaller than other parts of the needle body. During the candidate point generation stage, a multi-scale neighborhood local extremum detection method is first employed. Specifically, different thresholds are used for detection based on the area of ​​the needle tip and the needle body; for example, a small threshold is used for the smaller needle tip, and a large threshold is used for the larger needle body. In this embodiment, as an example, a 16-neighborhood can be used, such as... Figure 4 Figure A shows a real-world illustration of the 16-neighborhood of pixel P. While not limited to 16 neighbors, the 16-neighborhood offers a more comprehensive and accurate detection range compared to the 8-neighborhood. Red dots represent the current pixel coordinates P, blue dots indicate that the grayscale threshold difference with the center pixel meets the screening criteria, and yellow dots indicate that the grayscale difference exceeds the set threshold. Theoretically, in the absence of metal artifacts, all neighboring pixels should meet the screening criteria and appear blue; however, due to metal stripe artifacts and local highlight diffusion, some neighboring pixels may exhibit abnormal grayscale deviations. To avoid the erroneous rejection of real pinpoints during the screening process, this application introduces a tolerance parameter T to improve the robustness and completeness of pinpoint recognition under complex artifact conditions. Additionally, it should be noted that... Figure 4 In A, the blue dots are not all located on the circular line surrounding pixel P; in fact, as shown... Figure 4 The blue dot A can lie on the line. By limiting the neighborhood range and grayscale intensity threshold, false extrema points in non-target areas are filtered out, thus initially identifying candidate regions for the needle point. Let the intensity of the current pixel be... Its neighborhood set is A pixel is considered a candidate highlight pixel when it meets the conditions of Formula 3 below:

[0032] in, This represents the number of neighboring points that meet the criteria, with a threshold of T=6 to ensure robustness even when scanning at different angles causes changes in the cross-sectional size of the needle point. The initial candidate point set is thus obtained, as shown in Formula 4 below:

[0033] However, since calcification points and hemostatic coils can also form local extremes, further screening and correction are required in subsequent treatments.

[0034] In some specific embodiments, such as Figure 4 As shown in Figure A, when filtering high-brightness points, the red dot is the current point. The value of this point is compared with the values ​​of the surrounding 16 points in two concentric circles, and the number of times the current point is greater than the surrounding values ​​is recorded. When this number exceeds a set threshold, the current point is determined to be a high-brightness point. More specifically, in this embodiment, two concentric circles can be selected as representatives of the surrounding points, and the selected threshold for the number of comparisons can be 6. However, the number of circles and the threshold are not fixed; they are mainly determined based on experimental results. When using two concentric circles and a threshold of 6, the selected points basically meet the characteristics of high-brightness points. Simultaneously, because a maximum of 16 comparisons and threshold judgments are performed, small noise points and points in flat areas are filtered out. The threshold of 6 used here is not high compared to the maximum comparison count of 16, ensuring that small bright spots can also be successfully selected. However, this also brings disadvantages, namely, the number of identified bright spots will be more than the actual number, and some points surrounding actual high-brightness points will also be identified as bright spots. This problem is addressed by excluding and merging bright spots in subsequent processing. Specifically, as shown... Figure 4 B and Figure 4 Figure C shows a schematic diagram of the bounding box and detection box of the target region identified by the above-mentioned method for identifying candidate high-brightness points in a real pelvic CT image.

[0035] In some specific embodiments, the growth method of connected regions specifically includes: using candidate points As seed points, based on the gray-level similarity criterion The corresponding connected regions are constructed using the region growing method. (As shown in Formula 5 below):

[0036] in, Represents the image domain.

[0037] S203: Based on the geometric dimensions of the connected region, the connected region is divided into a first type of region for characterizing the cross-section of a single implantation needle, and a second type of region for characterizing the cross-sections of multiple implantation needles that are stuck together due to artifacts or adjacency. In some specific embodiments, candidate points can be classified based on the bounding box (BBox) of the growth region: Small areas refer to the first type of areas mentioned above ( The centroid can be directly calculated as the needle point position using the following formula 6:

[0038] in, Indicates the area The total number of pixels within.

[0039] Large areas refer to the second type of areas mentioned above ( Since it may contain multiple adhered pin points, clustering separation (such as K-means or DBSCAN) is required to group the pixels within a region into multiple sub-regions, each sub-region corresponding to an independent pin point. For example... Figure 5 As shown in the figure, the highlighted areas marked 1, 3, and 4 are relatively large, indicating that they may contain multiple adjacent pinpoint responses. To accurately identify and separate each individual pinpoint location, further clustering calculations are needed for the pixels within these large areas.

[0040] S204: For each of the second type of regions, the independent needle point coordinates of multiple implantation needle cross sections are separated and extracted using the method disclosed in the first aspect of this application; In some embodiments, when multiple needle points are adjacent in a large area, image adhesion often occurs due to the limited resolution of CT scans and the strong metallic stripe artifact effect, making it difficult for traditional region growing algorithms to effectively separate them. To address this problem, this application proposes a "multi-stage clustering and point separation strategy for intensity terrain models." Its core idea is to analogize the distribution of high-brightness points in an image to mountain peaks in terrain, and to distinguish adjacent peaks by setting a "level" parameter (minimumvalue level, minvalue). Unlike classic DBSCAN or mean shift clustering methods, this application's method emphasizes local gray-level transition characteristics rather than a global density threshold, thus better adapting to complex gray-level distributions caused by metallic artifacts.

[0041] like Figure 6 As shown, the gray intensity distribution of the needle point area is analogous to three-dimensional terrain. The four peaks marked in the figure (1, 2, 3, 4) correspond to four independent local gray-level maxima, i.e., four potential needle point locations. Peaks 1 and 3 represent the main needle points; peaks 2 and 4 are another set of actual needle points. Due to their proximity to other needle points during helical scanning, the photon flux of peaks 2 and 4 is abnormally reduced due to ray attenuation and beam hardening effects, resulting in lower HU values ​​and forming secondary peaks. Therefore, it is necessary to filter and separate the main and secondary peaks by setting the "horizontal plane" parameter (minvalue). Specifically, as shown... Figure 1As shown, the image feature point separation method based on an intensity terrain model disclosed in the first aspect of this application includes: S101: Obtain the target connected region of the image to be processed, the target connected region including multiple bright feature points that are stuck together due to imaging artifacts or spatial adjacency of objects; S102: Map the two-dimensional coordinates and grayscale intensity values ​​of each pixel within the target connected region to a three-dimensional intensity terrain model, wherein the grayscale intensity value is the elevation. S103: Based on the set dynamic grayscale height threshold, cluster the pixels in the three-dimensional intensity terrain model based on grayscale connectivity constraints. For any two pixels, traverse all possible connected paths between the two pixels. If there is at least one connected path and the grayscale intensity value of all pixels on the path is greater than or equal to the dynamic grayscale height threshold, then the two pixels are assigned to the same sub-cluster. In some embodiments, the dynamic grayscale height threshold is adaptively adjusted based on the local grayscale intensity distribution within the target connected region to distinguish between the primary grayscale peak formed by the real physical object and the secondary grayscale peak formed by imaging artifacts.

[0042] S104: Extract the pixel with the largest gray intensity value in each sub-cluster as the independent feature point representing that sub-cluster, and obtain the candidate set of independent feature points.

[0043] In some embodiments, the method further includes: If the number of independent feature point candidate sets is less than a preset minimum value and the geometric size of the target connected region is greater than a first size threshold, it is determined to be over-merging and a first supplementary strategy is executed to increase feature points; if the number of independent feature point candidate sets is equal to a preset intermediate value and the combined size of the target connected region is greater than a second size threshold, it is determined that there are missing peaks and a second supplementary strategy is executed to increase feature points; both the first and second supplementary strategies select local gray-level extreme points that do not overlap with existing independent feature points and do not satisfy the neighborhood connectivity condition as supplementary vertices through neighborhood exclusive search in the original set of high-brightness points.

[0044] Output the set of independent feature point coordinates after the supplementary strategy is executed.

[0045] S205: Integrate the centroid coordinates of the first type of region and the coordinates of all independent needle points separated from the second type of region to obtain a global candidate needle point set; In some more specific embodiments, the specific steps of S204 and S205 described above include: (1) Let the set of pixel coordinates of all candidate highlight points within the BBox be . As shown in Formula 7 below:

[0046] in, Indicates the first The coordinates of each pixel, and the corresponding CT grayscale intensity (HU value), are denoted as... .

[0047] Introducing a clustering operator based on gray-scale height constraints As shown in Formula 8 below:

[0048] Satisfy the following formula 9:

[0049] in, The grayscale height threshold. This represents a sub-cluster, corresponding to a potential pivot point.

[0050] Gray-scale connectivity constraints: For any two points If a connected path exists This ensures that the grayscale intensity of all points on the path is not lower than the threshold. That is, the following formula 10:

[0051] but and They were assigned to the same cluster. .

[0052] Here, This can be considered as the "horizontal plane height" in a grayscale terrain model: when At higher levels, retaining only the peak area is beneficial for separating needle points that are obviously adhered but have independent peak values; when At lower levels, adjacent bright areas are connected by gray-level valleys. Compared with traditional clustering methods based on Euclidean distance, this strategy is more in line with the imaging physical characteristics of metal needle points in CT images, which are characterized by steep local gray-level changes and peak dominance, and can avoid excessive fragmentation of a single needle point.

[0053] (2) Peak Extraction After completing the clustering, for each subset The pixel with the highest grayscale intensity is extracted as the peak point of the cluster, as shown in Formula 11 below:

[0054] in, Represents pixels CT grayscale intensity (HU value) at the location.

[0055] This yields the candidate set of needle points, as shown in Formula 12 below:

[0056] Each of them This corresponds to an independent potential needle point location.

[0057] (3) Vertex supplementation If the number of representative high-brightness points obtained after clustering is less than 2, and the corresponding bounding box (BBox) size is greater than or equal to 3×3 pixels, it indicates that the issue may be due to a grayscale height threshold. If the setting is too small, multiple pin points will be excessively merged. In this case, a local highest point q that does not overlap with the existing candidate points and does not satisfy neighborhood connectivity should be selected from the original candidate high point set V in the BBox through neighborhood exclusive search as a supplementary vertex to restore the excessively merged pin points.

[0058] (4) Search for the third highest highlight If the candidate point set obtained after clustering If the size is 2 (meaning two representative high-brightness points have been obtained), and the corresponding bounding box (BBox) size is greater than or equal to 20×20 pixels, it indicates that a third independent peak may exist in this region. This is because in a larger BBox, multiple pinpoints are more likely to stick together due to metallic artifacts, resulting in only two main peaks being detected during clustering. In this case, it is necessary to select from the original candidate high-brightness point set... In this process, a neighbor-exclusive search is used to select a point that is compatible with the existing candidate points. Local maxima that do not overlap and do not satisfy neighborhood connectivity This is used to add new high highlights to restore the missed needle points.

[0059] (5) Merging of high-brightness points If the final region contains only two candidate high-brightness points and Then, the square of its Euclidean distance is further calculated as shown in Formula 13 below:

[0060] Two points will be merged into one high point if one of the following conditions is met: If both points are located at the starting point of the needle, and ; If the two points are not located at the starting point of the needle, and .

[0061] This segmented constraint effectively balances the imaging differences caused by changes in the physical dimensions of the needle tip. The merged high-brightness area... Take the centroids of the two points, as shown in Formula 14 below:

[0062] The final set of needle points after candidate point screening, clustering, and correction is denoted by the following formula 15:

[0063] Each element represents a corrected needle point location. To achieve subsequent 3D reconstruction, the needle point coordinates in the image coordinate system need to be transformed to the patient's physical coordinate system.

[0064] Let the image coordinates be... The index of the current slice is According to the DICOM standard, physical coordinates The rigorous calculation formula is as follows, where P x The calculation is shown in Formula 16, P y The calculation is shown in Formula 17, P z The calculation is shown in Formula 18:

[0065]

[0066]

[0067] in: (S r , S c PixeSpacing (row / column spacing); (Slice thickness); That is, the three-dimensional coordinates of the top left corner of the slice in physical space; The orientation angle derived from ImageOrientationPatient describes the orientation of the image plane in three-dimensional space.

[0068] In the code implementation, to simplify calculations, a linear mapping as shown in the following formula is used (the negative sign is used to adapt to a specific clinical view coordinate system):

[0069]

[0070] In a single CT slice, when a large bounding box contains multiple spatially highly cohesive pinpoints (e.g., four pinpoints exist, but the algorithm only identifies two), a two-dimensional detection strategy relying solely on local clustering and gray-level extremum discrimination may theoretically fail to fully recover all real pinpoints. This phenomenon mainly stems from factors such as metal artifacts, partial volume effects, and enhanced gray-level continuity between pinpoints. These factors weaken the separability between pinpoints, causing multiple real pinpoints to appear as a single or a limited number of significant extrema in the local feature space.

[0071] In this embodiment, in the above preprocessing method, the image data values ​​are cropped by WinWidth=3600 and WinCenter=2200, and the values ​​in the window (set by window width and window level) are retained. Figure 11 These are comparison images before and after window width and position adjustments. Figure 11 A represents the state before processing. Figure 11 B represents the processed image. By comparing the images, it can be seen that the high-brightness areas are more concentrated and easier to segment after processing.

[0072] From an algorithm design perspective, while it is formally feasible to forcibly recover more pin points within a single slice, for example by lowering the threshold or introducing a higher-order sub-clustering strategy, such methods often significantly increase the risk of false detection in practical applications, especially in areas with severe artifacts or high noise levels. This may lead to the introduction of false pin points, which could adversely affect the subsequent pin trajectory reconstruction process.

[0073] Based on a comprehensive consideration of detection stability and overall reconstruction reliability, the method in this application intentionally constrains the number of needle points recovered in a single slice, but does not require the complete recovery of all needle points in extreme adhesion cases. Furthermore, robust optimization is used to further reconstruct the data, so as to achieve dual verification of local feature extrema and global geometric consistency, and collaboratively optimize the real data of the implanted needle path.

[0074] S206: An improved random sampling consensus algorithm is used to perform spatial path clustering and fitting on the global candidate needle point set to obtain the spatial trajectory of the applicator.

[0075] In some embodiments, the step of using an improved random sampling consensus algorithm to perform spatial path clustering and fitting on the global candidate needle point set includes: a: In a single iteration, two discrete points are randomly selected from the global candidate needle point set, and straight line models of the three projection planes XY, YZ, and XZ constructed based on the image coordinate system are fitted respectively, and the geometric deviations of all points in the global candidate needle point set to the straight line model are calculated. b: Filter out the set of inner points according to the preset first deviation threshold, and optimize the duplicate points with the same Z coordinate in the set of inner points, retain the point with the smallest set deviation, and eliminate the adhesion interference at the intersection of multiple needles. c: Repeat steps a and b above for a preset number of times to generate the set of interior points obtained in each iteration, and obtain the total set of interior points; d: From the total set of interior points, extract the set with the largest number of points as the target insertion needle cluster, and remove the points contained in the target insertion needle cluster from the global candidate needle point set. Repeat this process until a preset number of trajectory clusters are reached. e: The remaining candidate points in the global candidate needle point set that are not included in any cluster are then reassigned based on their geometric deviations from each target needle cluster using a second geometric deviation threshold, thus completing the feature point set of the curved needle path. Furthermore, a coordinate system transformation occurs during the fitting process. Specifically, the aforementioned XYZ coordinate system is based on the CT image coordinate system, while the three-dimensional intensity terrain coordinate system constructed by the centroid in a single CT image is based on a two-dimensional pixel system. The conversion between the two can be performed using conventional methods from existing technologies.

[0076] In some embodiments, the geometric deviation includes: calculating the vertical distances from the discrete point to the fitted lines of the XY plane, the YZ plane, and the XZ plane respectively, and then performing a linear weighted combination of the three vertical distances to obtain the comprehensive deviation value of the discrete point relative to the linear model. That is, the geometric deviation serves as the basis for subsequent interior point selection.

[0077] In some more specific embodiments, the algorithm implementation principle of S206 based on improved RANSAC for robust fitting of the three-dimensional path of the implantation needle includes: To address the discrete 3D coordinate point set obtained after single-slice needle point detection, and considering the common problems in clinical scenarios such as multi-needle cross-distribution, local bending deformation of single needles, and difficulty in distinguishing adhesion points and insignificant needle tract centroid features due to spatial adhesion of multiple needles, this application proposes an improved Random Sample Consensus (RANSAC) algorithm to achieve spatial path clustering and robust fitting of implanted needles. This algorithm overcomes the dependence of traditional fitting methods on centroid discrimination, appropriately fitting the bending morphology of implanted needles through local segmentation, and effectively separating adhesion point clouds by combining distance filtering and duplicate point optimization, while simultaneously solving the problem of outlier interference caused by metal artifacts. The algorithm first sets the number of iterations to 10,000. In each iteration, two points are randomly selected to fit straight lines on the XY, YZ, and XZ two-dimensional planes in three-dimensional space. The bending trajectory of the insertion needle is approximated using a local segmental fitting method. The combined distance from all discrete points to the fitted straight line is calculated, and inlier points are filtered using a distance threshold of 0.8. Simultaneously, the length and number distribution characteristics of the inlier point set in each iteration are recorded. Then, the Z-axis is repeated on the filtered inlier point set. The adhesion points of the coordinates are optimized. Even when the adhesion points are not effectively distinguished and the centroid features of the needle path are unclear, the points with the smallest distance are retained and the other points with larger distances are removed based on the parameters of the number and length of the inner points. This eliminates the fitting interference caused by the cross adhesion of multiple needles and achieves effective separation of the adhesion points. Then, all iteration results are sorted from the smallest to the largest number of inner points. The number and length of the implanted needles are used as termination conditions (as mentioned above, the number and length of implanted needles in each CT image are known). The cluster with the largest number of inner points is extracted in turn and the clustered points are removed. This ensures that even in complex scenarios where three implanted needles are distributed crosswise and one of them is accompanied by local curvature, accurate clustering can still be completed based on the length and number parameters. Finally, the remaining un-clustered points are assigned to the nearest implanted needle cluster through a secondary allocation strategy (distance threshold 2.5) to further complete the feature point set of the curved needle path. After finally achieving the 3D path fitting of the implanted needle, the fitting effect was verified by filtering out the minimum distance threshold of 0.03, calculating the mean distance, counting the number of points in each cluster (length distribution), and drawing a distance histogram (50 bins). This verified the robustness of the algorithm in scenarios where the adhesion point distinction is poor and the centroid features are not significant.

[0078] In some embodiments, the specific implementation details of the above algorithm include: In the preprocessing stage, the discrete three-dimensional coordinate point set of multiple implanted needles (including a small number of outliers introduced by metal artifacts) output by the needle point detection is read first. Two discrete points are randomly selected, and the linear models of the three two-dimensional planes XY, YZ, and XZ are fitted by the least squares method respectively. The perpendicular distance from all points to each two-dimensional fitted line is calculated and combined into the equivalent distance in three-dimensional space. A distance threshold of 0.8 is set to filter the initial inner point set. At the same time, duplicate points under the same Z coordinate are optimized, and points that are far away from the model among the duplicate Z coordinate points are removed (only the closest points are kept). After 10,000 iterations of sampling-fitting-screening (ensuring the algorithm converges to the optimal solution), all iteration results are sorted according to the number of inlier sets. The model with the largest number of inlier sets is selected step by step, and intersection points are removed. Finally, the non-intersection point set corresponding to the inserted needles is clustered. For the remaining un-clustered points, secondary allocation is performed based on their three-dimensional equivalent distance to each needle fitting path (threshold 2.5). Finally, the needle point set is visualized in three dimensions. By calculating the average distance of the inlier set, the length distribution of the inlier set, and the distance histogram during the iteration process, the fitting effect is quantitatively analyzed and verified.

[0079] In some more specific embodiments, the core steps of the above algorithm are implemented through the following steps: (1) Least squares three-dimensional linear fitting The least squares method was used to fit a three-dimensional line to two randomly selected points. Linear regressions were then performed on the XY, YZ, and XZ planes to solve for the linear model parameters. The fitting formula is as follows:

[0080] in, These are the parameters for fitting the XY plane. These are the parameters for fitting the YZ plane. These are the fitting parameters for the XZ plane. The Euclidean distance from a point to a 3D line is obtained by weighting the distances from the three planes, and the calculation formula is:

[0081] in, , , These are the perpendicular distances from the point to the three fitted planes.

[0082] (2) Improve the core RANSAC process Two discrete points are randomly selected to construct an initial straight-line model. Local fitting primitives are used to adapt to the local bending morphology of the implanted needles. The distance from all points to the model is calculated, and the set of inner points with a distance less than 0.8 mm is selected. To address the problem of ineffective differentiation of adhesion points and difficulty in identifying the centroid of the needle path in multi-needle intersection scenarios, the adhesion and overlapping points under the same Z coordinate are optimized. By relying on the length and number parameters of the inner points, only the closest points are retained to achieve separation of adhesion points, breaking through the dependence of traditional methods on the centroid differentiation. After 10,000 iterations of sampling-fitting-selection, the optimal model with the largest number of inner points is retained. Point sets are extracted sequentially from the iteration results to complete the non-intersection clustering of 8 implanted needles. Even in complex scenarios where three implanted needles are intersecting and one of them is bent, the overall bending needle path can still be accurately approximated through feature selection using length and number parameters. For the remaining unclustered points, a secondary allocation is performed based on their distance to the fitted path (≤ 2.5 mm), further improving the integrity of the curved needle path fitting and the robustness of the adhesion point handling.

[0083] (3) Quantitative and Visual Analysis Quantitative analysis: Calculate the distance from each set of points to the fitted path, remove redundant points with a distance less than or equal to 0.03 mm, and calculate the average distance as an evaluation index of path fitting accuracy; at the same time, statistically analyze the changing patterns of the number of internal point sets and the average distance during the iteration process to analyze the convergence of the algorithm. b. 3D visualization: A 3D coordinate system is constructed using the Matplotlib library to perform color clustering visualization of the feature point set of 8 implantation needles. At the same time, discrete points of the same needle point set are connected to generate a 3D path fitting curve for the implantation needles. c. Result saving: The feature point set of each implanted needle is saved as a CSV file to provide digital coordinate data for subsequent radiotherapy planning.

[0084] In some embodiments, the aforementioned preset thresholds are obtained through training with training set samples or set based on experience. They can be specific thresholds or ranges, and their specific forms are not specifically limited in this embodiment.

[0085] In some embodiments, the method further includes: converting the spatial trajectory of each implantation needle fitted by S206 into data with clear geometric relationships and clinical significance, and encoding it into an RTPLAN file conforming to the DICOM standard.

[0086] In some more specific embodiments, the process is as follows: First, the spatial trajectory of each implantation needle is represented as a three-dimensional physical coordinate sequence in the target TPS coordinate system, and the trajectory points are oriented to be consistent according to spatial continuity constraints. Then, based on the spatial positional relationship of the implantation needles or a preset numbering rule, a corresponding channel topology relationship is established for each implantation needle, and a source generator centerline path model is constructed.

[0087] Furthermore, based on the applicator centerline path model, path coordinate information that can characterize the geometric shape of the applicator is extracted, including path point sequence, path direction and effective length parameters, and a spatial coordinate system consistent with the original image reconstruction results is maintained, thereby forming a unified geometric data chain from image detection results to applicator path model, avoiding cumulative errors caused by multiple coordinate transformations.

[0088] Finally, the system uses the DICOM encoding module to organize and encapsulate the applicator path coordinates, channel number, effective length, and related geometric parameters according to the DICOM RTPLAN data structure, automatically generating an RTPLAN file that conforms to the DICOM standard requirements. This file can be directly read and used by the target TPS, eliminating the need for manual reconstruction in the traditional workflow and reducing reconstruction errors caused by subjective judgment and visual fatigue among clinical staff.

[0089] In some embodiments, the experimental results and analysis of specific experiments verified by the above methods are as follows: 1. Limitations of Initial Point Clouds and Traditional Methods The original distribution of the discrete point cloud of the implanted needles is as follows: Figure 7 As shown in Figure A, all three-dimensional coordinate points are uniformly presented in gray. Eight implanted needles are visible and intersecting in space. One of them is accompanied by local bending and deformation. At the same time, there are problems such as point cloud overlap caused by multiple needles sticking together and outlier scattered points caused by metal artifacts. The centroid features of the needle path are not significant, and the spatial paths of each implanted needle cannot be directly distinguished by visual means.

[0090] To verify the limitations of traditional fitting methods, a single straight line was fitted to all points using the global least squares method. The results are as follows: Figure 7 As shown in B, the red line represents the traditional overall fitting result. It can be seen that this method can only generate a single fitted line, which cannot adapt to the local curvature of the implanted needle. Furthermore, the adhesion and cross distribution of multiple needles cause the fitted line to deviate significantly from the actual needle path. At the same time, it cannot distinguish the spatial paths of different implanted needles, highlighting the adaptation defects of the traditional method in complex clinical scenarios.

[0091] 2. Outlier Preprocessing Effect Outlier removal was performed on the original point cloud based on a distance threshold (0.8). The preprocessing result is as follows: Figure 8As shown in Figure A, the blue dots represent the retained valid implanted needle feature points, while the light gray dots represent outliers that were removed (mostly caused by CT metal artifacts or needle detection errors). After preprocessing, the point cloud is free from noise interference, and the core features (multiple needle crossings, local bending, and point cloud adhesion) are fully preserved, providing a clean data foundation for subsequent improvements to the RANSAC algorithm's clustering and bending fitting.

[0092] 3. Improve the fitting and clustering results of the RANSAC algorithm. After removing outliers, an improved RANSAC algorithm was used to perform three-dimensional path clustering and fitting of the implanted needles. The results are as follows: Figure 8 As shown in Figure B, the algorithm successfully clustered the mixed discrete 3D coordinate points into 8 groups, each corresponding to the spatial feature points of 8 implantation needles. Furthermore, based on the feature points of the same implantation needle, a continuous spatial path can be effectively fitted, with no obvious clustering errors or outlier interference.

[0093] From a spatial distribution perspective, the paths of different implanted needles are crisscrossed in three-dimensional space, highly consistent with clinical implantation strategies, verifying the algorithm's ability to cluster complex distribution sets of needle points. Feature points of the same implanted needle are closely distributed along the fitted path, with some needle tracks exhibiting slight curvature. This indicates that the algorithm can still effectively identify the spatial continuity of needle points under metal artifact interference and adapt to local bending deformations of the implanted needles, providing a reliable geometric basis for subsequent three-dimensional reconstruction of the applicator.

[0094] Furthermore, a small number of isolated outliers not initially clustered are visible in the image. These outliers are mostly caused by metal artifacts in CT images or needle point detection. The algorithm effectively merges these outliers into the nearest needle cluster through a secondary allocation strategy (distance threshold 2.5), further improving the completeness and robustness of the path fitting. The algorithm not only achieves accurate clustering under multi-needle interleaving distribution, but also completes accurate needle path fitting based on in-point length and quantity parameters in complex scenarios where three needles are interleaved, one of which is locally bent, adhesion points are not effectively distinguished, and needle centroid features are insignificant. This effectively eliminates path fitting bias caused by multi-needle adhesion and difficulty in distinguishing centroids, ensuring consistency between the fitting results and the actual spatial morphology of the needles. Overall, this visualization clearly demonstrates the advantages of the improved RANSAC algorithm in multi-needle clustering, 3D line fitting, and outlier handling, providing strong support for accurate 3D reconstruction of clinical brachytherapy applicators.

[0095] The algorithm not only achieves accurate clustering under the condition of multiple needles interlaced, but also completes the accurate fitting of the needle path by relying on the parameters of inner point length and number in complex scenarios where three implantation needles are interlaced, one of which is accompanied by local bending, the adhesion points are not effectively distinguished, and the needle path centroid features are not significant. It effectively adapts to the local bending and deformation characteristics of implantation needles in clinical practice, eliminates the path fitting deviation caused by multiple needle adhesion and difficulty in distinguishing the centroid, and ensures the consistency between the fitting results and the actual spatial morphology of the implantation needles.

[0096] 4. Evaluation of iterative convergence like Figure 9 As shown in Figure A, in the early stages of iteration (the first 2000 samplings), the average distance from the inlier set to the fitted line fluctuates wildly and is relatively high, indicating that the random sampling model has not yet converged to a stable cluster of points. As the number of iterations increases, the average distance generally converges to a lower range, the proportion of high-value points decreases significantly, and it stabilizes at a low level of 0.3~0.6 mm, indicating that after the algorithm converges, the inlier set fits the fitted line well, and the path fitting accuracy is stable and reliable.

[0097] like Figure 9 As shown in B, the number of interior points increases stepwise with the iteration index, eventually converging to approximately 40-45 points. This indicates that during the iteration process, the algorithm can gradually find a better fitting model and effectively identify more valid points belonging to the same insertion needle, demonstrating good convergence and robustness.

[0098] Based on the quantitative analysis results of the two figures, the improved RANSAC algorithm can effectively converge to the optimal fitting model during the iteration process while maintaining high fitting accuracy, providing a reliable guarantee for the three-dimensional path reconstruction of the implanted needle.

[0099] 5. Fitting accuracy analysis Distance statistics were performed on the 8 sets of needle points fitted by the improved RANSAC algorithm. The distance distribution from each point to the fitted path of the corresponding implanted needle is shown below. Figure 10 As shown in the figure, the distances are mainly concentrated in the range of 0–0.8 mm, peaking in the range of 0.7–0.8 mm, indicating that the vast majority of points closely fit the fitted path. Statistically, the average fitting error is approximately 0.4–0.5 mm, far below the accuracy threshold for applicator reconstruction in clinical brachytherapy (typically 1.0 mm), verifying that the algorithm can still achieve high-precision 3D path fitting under the interference of metal artifacts, meeting the accuracy requirements for clinical applications.

[0100] In addition, with Figure 4 B and Figure 4 Based on the area shown in C, the results of manual annotation are as follows: Figure 12 As shown, the results of automatic annotation are as follows: Figure 13As shown, it is similar to Figure 12 Manual annotation ( Figure 12 In the comparison (where the number on the top is 6 and the number on the bottom is 4), the speed is faster while maintaining comparable accuracy; it should be noted that... Figure 13 Different colored boxes are used to represent them. Figure 13 There is partial adhesion between the edges of multiple frames of different colors (for example, in...). Figure 13 The color boxes from left to right are green, yellow, red, bright purple, red, blue, purple, and yellowish-green. However, due to partial adhesion, the red box (the third one from the left) is not fully displayed in the image. But in actual clinical software, a vector graphic is displayed, and selecting different boxes after zooming in works without issue. (Comparison) Figure 12 and Figure 13 It can be seen that after the screenshot Figure 12 The visual effect is clearer and better, which is reflected in Figure 13 The rectangles in the image partially overlap or are glued together. The reason for this is: Figure 12 For images displayed in commercial medical TPS systems, due to the combined effects of hardware graphics processing capabilities and the image filtering and resampling algorithms built into the planning system, artifacts related to metal applicators are weakened to a certain extent at the display level, thus achieving a better visual observation effect (not that artifacts are essentially eliminated, but rather that they are weakened at the display level). Figure 13 The display results are processed on a general PC platform. Due to the limitations of ordinary PC hardware performance and image compression processing, the visual display effect of this application is inferior to that of commercial TPS systems, especially in the overlapping areas of multiple rectangles. However, this difference is mainly reflected in the image display level and does not significantly affect the accuracy of the applicator center positioning based on the centroid calculation of the rectangular area. The results are quite consistent with manual annotation and are faster.

[0101] 6. Algorithm Implementation Code Description This algorithm can be developed using different languages, and various programming languages ​​can achieve the goals set forth in this core concept. In this embodiment, it is primarily developed using Python. In practical applications, the Pandas library is mainly used for reading and preprocessing coordinate data, the NumPy library for performing matrix operations and other numerical calculations, and the Matplotlib library for visualizing the 3D space. Based on the practical needs of fitting the 3D path of the implanted needle, the algorithm code is built with four core functional modules: a least-squares fitting function, a repeat point detection function, an improved RANSAC iterative function, and a 3D visualization function. The specific implementation idea is as follows: In the data reading phase, Pandas was used to read a text file containing the 3D coordinates of discrete pinpoints, extracting the X, Y, and Z coordinate data columns to prepare the foundational data for subsequent fitting and clustering calculations. In the function definition phase, the `linefit2dpoint` function was written to fit the projected lines of the XY, YZ, and XZ two-dimensional planes using the least squares method, providing a foundation for subsequent 3D spatial line fitting. Simultaneously, the `find_duplicates_with_indices` function was written to specifically detect duplicate coordinate points in the Z-axis direction, calculating the distance from duplicate points to the fitted line and filtering them, retaining the points with the smallest distance and removing redundant points with larger distances to reduce the interference of duplicate coordinates on the fitting. In the iterative clustering phase, 10,000 iterations were set. In each iteration, two sample points were randomly selected, a 2D projected line was fitted, and then the inner point set was filtered based on a distance threshold of 0.8 mm, thus completing 8 iterations. The algorithm performs the following steps: Initial clustering of the root needle point set; Remaining point allocation stage: For discrete points not included in the initial cluster, calculate the distance from each root needle to the fitted line, and assign points with a distance ≤ 2.5 mm to the corresponding needle point set to complete the feature point data for each root needle; Visualization and saving stage: Use Matplotlib's 3D plotting function to present the clustering results and path fitting effect, and export the feature point set of each root needle to a CSV file for easy subsequent data verification and reuse; Quantitative analysis stage: Statistically analyze the changes in the number of internal points and average distance during the iteration process, and plot a histogram of the distance distribution from points to the fitted path to evaluate the convergence and fitting accuracy of the algorithm.

[0102] The core logic and parameter settings of the algorithm code correspond one-to-one with the algorithm principles described above, ensuring that the experimental process and results can be repeatedly verified. At the same time, the code has strong practicality and extensibility, and can flexibly adjust key parameters such as distance threshold and number of iterations according to different needs of actual clinical application scenarios, adapting to different clinical data processing scenarios.

[0103] A third aspect of this application discloses a computer device, which may include: one or more processors and one or more memories; wherein the memories store computer-readable code that, when run by the one or more processors, can perform the methods described above.

[0104] The processor in this embodiment can be an integrated circuit chip with signal processing capabilities. The processor can be a general-purpose processor, a digital signal processor (DSP), an application-specific integrated circuit (ASIC), an off-the-shelf programmable gate array (FPGA), or other programmable logic devices, discrete gate or transistor logic devices, or discrete hardware components. It can implement or execute the methods, operations, and logic block diagrams disclosed in this embodiment. The general-purpose processor can be a microprocessor or any conventional processor, and can be based on an x86 or ARM architecture.

[0105] In general, the various exemplary embodiments of this disclosure can be implemented in hardware or dedicated circuitry, software, firmware, logic, or any combination thereof. Some aspects can be implemented in hardware, while others can be implemented in firmware or software that can be executed by a controller, microprocessor, or other computing device. When aspects of embodiments of this disclosure are illustrated or described as block diagrams, flowcharts, or using some other graphical representation, it will be understood that the blocks, apparatuses, systems, techniques, or methods described herein can be implemented as non-limiting examples in hardware, software, firmware, dedicated circuitry or logic, general-purpose hardware or controllers or other computing devices, or some combination thereof.

[0106] For example, the methods or apparatus according to embodiments of this disclosure can also be implemented using the architecture of a computing device. A computing device may include a bus, one or more CPUs, read-only memory (ROM), random access memory (RAM), a communication port connected to a network, input / output components, a hard disk, etc. Storage devices in the computing device, such as ROM or hard disk, may store various data or files used for processing and / or communication of the methods provided in this disclosure, as well as program instructions executed by the CPU. The computing device may also include a user interface. Of course, the architecture described above is merely exemplary, and one or more components of the computing device described above may be omitted as needed when implementing different devices.

[0107] This invention also provides a computer-readable storage medium storing computer-readable instructions. When the computer-readable instructions are executed by a processor, the methods disclosed in this embodiment can be performed. The computer-readable storage medium in this embodiment can be volatile memory or non-volatile memory, or may include both volatile and non-volatile memory. Non-volatile memory can be read-only memory (ROM), programmable read-only memory (PROM), erasable programmable read-only memory (EPROM), electrically erasable programmable read-only memory (EEPROM), or flash memory. Volatile memory can be random access memory (RAM), which is used as an external cache. By way of example, but not limitation, many forms of RAM are available, such as static random access memory (SRAM), dynamic random access memory (DRAM), synchronous dynamic random access memory (SDRAM), double data rate synchronous dynamic random access memory (DDRSDRAM), enhanced synchronous dynamic random access memory (ESDRAM), synchronous interconnected dynamic random access memory (SLDRAM), and direct memory bus random access memory (DRRAM). It should be noted that the memory of the methods described in this application is intended to include, but is not limited to, these and any other suitable types of memory.

[0108] This disclosure also provides a computer program product or system, including a computer program that, when executed by a processor, implements the steps of the above-described method. The computer program product or computer program is used to execute an image feature point separation system based on an intensity terrain model. The system includes: A connected region acquisition module is used or configured to acquire a target connected region of an image to be processed, the target connected region including multiple bright feature points that are stuck together due to imaging artifacts or spatial adjacency of objects; A terrain model generation module is used or configured to map the two-dimensional coordinates and grayscale intensity values ​​of each pixel within the target connected region to a three-dimensional intensity terrain model, wherein the grayscale intensity value is the elevation. The pixel clustering module is used or configured to cluster pixels in the three-dimensional intensity terrain model based on grayscale connectivity constraints according to a set dynamic grayscale height threshold. For any two pixels, if there is at least one connected path and the grayscale intensity value of all pixels on the path is greater than or equal to the dynamic grayscale height threshold, then the two pixels are assigned to the same sub-cluster. The feature point extraction module is used or configured to extract the pixel with the largest gray intensity value in each sub-cluster as the independent feature point representing that sub-cluster, thus obtaining a candidate set of independent feature points.

[0109] In some embodiments, the system further includes: The feature point correction module is used or configured to determine that if the number of independent feature point candidate sets is less than a preset minimum value and the geometric size of the target connected region is greater than a first size threshold, it is considered over-merging and a first supplementary strategy is executed to increase feature points; if the number of independent feature point candidate sets is equal to a preset intermediate value and the combined size of the target connected region is greater than a second size threshold, it is considered that there are missing peaks and a second supplementary strategy is executed to increase feature points. The feature point coordinate set output module is used or configured to output the independent feature point coordinate set after the supplementary strategy is executed.

[0110] In some embodiments, a computer program product or computer program is used to perform an applicator three-dimensional image reconstruction system, the system comprising: A medical image acquisition module is used or configured to acquire CT medical images containing metal implanted needles, perform standardized preprocessing on the medical images, eliminate device differences, and enhance tissue contrast. The connected region construction module is used or configured to perform pixel-level analysis on the preprocessed single-frame image, identify initial candidate high-brightness points, and construct corresponding connected regions using the initial candidate high-brightness points as seed points and a region growing method to obtain multiple connected regions. A connected region processing module is used or configured to divide the connected region into a first type of region for characterizing the cross-section of a single implantation needle, and a second type of region for characterizing the cross-sections of multiple implantation needles that are stuck together due to artifacts or adjacency, based on the geometric dimensions of the connected region. The feature point separation module is used or configured to, for each of the second type of regions, use the method disclosed in the first aspect of this application to separate and extract the independent needle point coordinates of multiple implantation needle cross sections; The candidate needle point extraction module is used or configured to integrate the centroid coordinates of the first type of region and the coordinates of all independent needle points separated from the second type of region to obtain a global candidate needle point set; The applicator spatial trajectory generation module is used or configured to perform spatial path clustering and fitting on the global candidate needle point set using an improved random sampling consensus algorithm to obtain the spatial trajectory of the applicator.

[0111] The exemplary embodiments of this disclosure described in detail above are merely illustrative and not restrictive. Those skilled in the art will understand that various modifications and combinations can be made to these embodiments or their features without departing from the principles and spirit of this disclosure, and such modifications should fall within the scope of this disclosure.

Claims

1. A method for separating image feature points based on an intensity terrain model, characterized in that, include: Obtain the target connected region of the image to be processed, the target connected region including multiple bright feature points that are stuck together due to imaging artifacts or spatial adjacency of objects; The two-dimensional coordinates and grayscale intensity values ​​of each pixel within the target connected region are mapped to a three-dimensional intensity terrain model, where the grayscale intensity value is the elevation. Based on the set dynamic grayscale height threshold, the pixels in the three-dimensional intensity terrain model are clustered according to the grayscale connectivity constraint. For any two pixels, if there is at least one connected path and the grayscale intensity value of all pixels on the path is greater than or equal to the dynamic grayscale height threshold, then the two pixels are assigned to the same sub-cluster. Extract the pixel with the largest gray intensity value in each sub-cluster as the independent feature point representing that sub-cluster, and obtain the candidate set of independent feature points.

2. The image feature point separation method based on intensity terrain model according to claim 1, characterized in that, The method further includes: If the number of independent feature point candidate sets is less than a preset minimum value and the geometric size of the target connected region is greater than a first size threshold, it is determined to be over-merging and a first supplementary strategy is executed to increase the number of feature points; if the number of independent feature point candidate sets is equal to a preset intermediate value and the combined size of the target connected region is greater than a second size threshold, it is determined that there is a missing peak and a second supplementary strategy is executed to increase the number of feature points. Output the set of independent feature point coordinates after the supplementary strategy is executed.

3. The image feature point separation method based on intensity terrain model according to claim 1, characterized in that, The dynamic grayscale height threshold is adaptively adjusted based on the local grayscale intensity distribution within the target connected region to distinguish between the primary grayscale peak formed by the real physical object and the secondary grayscale peak formed by imaging artifacts.

4. The image feature point separation method based on intensity terrain model according to claim 2, characterized in that, Both the first and second supplementary strategies select local gray-level extreme points that do not overlap with existing independent feature points and do not satisfy the neighborhood connectivity condition as supplementary vertices by performing a neighborhood exclusive search in the original set of high-brightness points.

5. A method for reconstructing a three-dimensional image of an applicator, characterized in that, include: Acquire CT medical images containing metal implanted needles, and perform standardized preprocessing on the medical images to eliminate device differences and enhance tissue contrast; Pixel-level analysis is performed on the preprocessed single-frame image to identify initial candidate high-brightness points. Then, using the initial candidate high-brightness points as seed points, a region growing method is used to construct corresponding connected regions, resulting in multiple connected regions. Based on the geometric dimensions of the connected regions, the connected regions are divided into a first type of region for characterizing the cross-section of a single implantation needle, and a second type of region for characterizing the cross-sections of multiple implantation needles that are stuck together due to artifacts or adjacency. For each of the second type of regions, the independent needle point coordinates of multiple implantation needle cross sections are separated and extracted using the method described in any one of claims 1-4; By integrating the centroid coordinates of the first type of region and the coordinates of all independent needle points separated from the second type of region, a global candidate needle point set is obtained; An improved random sampling consensus algorithm is used to perform spatial path clustering and fitting on the global candidate needle point set to obtain the spatial trajectory of the applicator.

6. The method for reconstructing a three-dimensional image of an applicator according to claim 5, characterized in that, The step of using an improved random sampling consensus algorithm to perform spatial path clustering and fitting on the global candidate needle point set includes: a: In a single iteration, two discrete points are randomly selected from the global candidate needle point set, and straight line models of the three projection planes XY, YZ, and XZ constructed based on the image coordinate system are fitted respectively, and the geometric deviations of all points in the global candidate needle point set to the straight line model are calculated. b: Filter out the set of inner points according to the preset first deviation threshold, and optimize the duplicate points with the same Z coordinate in the set of inner points, retain the point with the smallest set deviation, and eliminate the adhesion interference at the intersection of multiple needles. c: Repeat steps a and b above for a preset number of times to generate the set of interior points obtained in each iteration, and obtain the total set of interior points; d: From the total set of interior points, extract the set with the largest number of points as the target insertion needle cluster, and remove the points contained in the target insertion needle cluster from the global candidate needle point set. Repeat this process until a preset number of trajectory clusters are reached. e: The remaining candidate points in the global candidate needle point set that have not been included in any cluster are reassigned based on their geometric deviation from each target difference needle cluster using a second geometric deviation threshold to complete the feature point set of the curved needle path.

7. The method for reconstructing a three-dimensional image of an applicator according to claim 6, characterized in that, The geometric deviation includes: calculating the perpendicular distances from the discrete point to the fitted lines of the XY plane, the YZ plane, and the XZ plane, respectively, and then performing a linear weighted combination of the above three perpendicular distances to obtain the comprehensive deviation value of the discrete point relative to the linear model. That is, the geometric deviation serves as the basis for subsequent interior point selection.

8. A computer device, characterized in that, The device includes: a memory and a processor; the memory is used to store a computer program; the processor executes the computer program to implement the steps of the method according to any one of claims 1-7.

9. A computer-readable storage medium, characterized in that, It stores a computer program that, when executed by a processor, implements the steps of the method as described in any one of claims 1-7.

10. A computer program product, comprising a computer program, characterized in that, When executed by a processor, the computer program implements the steps of the method described in any one of claims 1-7.

Images