Radiotherapy target region adaptive delineation method and system based on CBCT image

Abstract

The application provides a radiotherapy target area adaptive delineation method and system based on a CBCT image, and relates to the technical field of image processing.The method comprises the following steps: acquiring a reference positioning CT image and a current fraction CBCT image, performing multi-scale spatial frequency domain alignment and feature extraction processing, and obtaining an initial spatiotemporal registration feature manifold; performing high-dimensional topological reconstruction and feature affinity calculation on the initial spatiotemporal registration feature manifold, constructing a spatial interaction intensity mapping model, and converting the spatial interaction intensity mapping model into a local structure potential energy gradient distribution to obtain a local heterogeneity potential field gradient distribution matrix.Through multi-scale frequency domain registration, high-dimensional topological modeling, curvature analysis and dynamic deformation guidance, the application can effectively suppress CBCT image artifact noise and adapt to organ nonlinear deformation, automatically generate a radiotherapy target area contour without holes, reduce manual repair, and improve adaptive delineation accuracy and efficiency.

Description

Technical Field

[0001] This invention relates to the field of image processing technology, and in particular to a method and system for adaptive delineation of radiotherapy target areas based on CBCT images. Background Technology

[0002] In cone-beam computed tomography (CBCT) guided radiotherapy, existing methods typically employ non-rigid image registration techniques, such as B-spline-based free deformation (FFD) models or optical flow-based deformation registration methods, to adaptively map the target contour delineated on the reference positioning CT image onto the current fractionated CBCT image. These methods, by deriving the spatial deformation relationship between the two images and then transferring the reference target contour to the CBCT image via a deformation field, achieve, to some extent, fractional updating of the target area.

[0003] However, in practical clinical applications, the existing methods described above still have room for improvement in certain situations. CBCT images themselves have physical imaging characteristics such as significant scattering artifacts, relatively low soft tissue contrast, and high noise levels. Simultaneously, during fractionated treatment, the patient's anatomical structures may undergo local nonlinear deformations, such as changes in organ filling, tumor regression, and changes in the volume of adjacent hollow organs (e.g., rectum, bladder). When the target area is adjacent to regions with strong local heterogeneity, such as bone, gas, or soft tissue-bone interfaces, conventional uniform smoothing constraints are insufficient to simultaneously maintain the continuity of the deformation field and the fidelity of local geometry. Near interfaces with large deformation gradients, especially in topological transition regions from solid tissue to cavities or bone, the deformation vector field derived by existing methods sometimes exhibits non-physical local abrupt changes or directional disorder. This instability of the deformation field directly leads to the creation of minute holes, boundary breaks, or abnormal depressions in the mapped target area contour at corresponding locations, thereby disrupting the integrity of the target area contour.

[0004] For example, in online adaptive radiotherapy for prostate cancer, the rectal gas volume may change significantly between fractions. When the rectum is significantly dilated, it compresses the prostate target area posteriorly, causing localized deformation of a small area but with a noticeable gradient at the junction between the posterior border of the prostate and the anterior wall of the rectum. When using conventional B-spline deformation registration methods, this junction area, being a topological transition zone between soft tissue structures and air-containing cavities, is prone to local discontinuities in the deformation field. These discontinuities result in sporadic holes or micro-cracks in the mapped clinical target volume (CTV) contour at the edge adjacent to the rectum. Although these defects are small in area, they are quite noticeable on thin-slice CBCT images, usually requiring manual retouching by physicists or physicians on each image layer, which to some extent affects the efficiency of the online adaptive delineation process. Summary of the Invention

[0005] This invention provides a method and system for adaptive delineation of radiotherapy target areas based on CBCT images, which improves the automation level of adaptive delineation of radiotherapy fractionation target areas.

[0006] To solve the above-mentioned technical problems, the technical solution of the present invention is as follows:

[0007] Firstly, a method for adaptive delineation of radiotherapy target areas based on CBCT images, the method comprising:

[0008] Step 1: Acquire the reference positioning CT image and the current fractional CBCT image, perform multi-scale spatial frequency domain alignment and feature extraction processing to obtain the initial spatiotemporal registration feature manifold;

[0009] Step 2: Perform high-dimensional topological reconstruction and feature affinity calculation on the initial spatiotemporal registration feature manifold, construct a spatial interaction intensity mapping model and transform the spatial interaction intensity mapping model into a local structural potential energy gradient distribution to obtain the local heterogeneous potential field gradient distribution matrix;

[0010] Step 3: Perform curvature tensor analytical operations and spatiotemporal state projection on the local heterogeneous potential field gradient distribution matrix, extract local curvature eigenvalues, and solve for the extreme value feature nodes of the principal deformation axis, the inflection point feature nodes of the secondary deformation axis, and the topological transition region feature nodes; derive the initial node deformation vectors based on the projection coordinate difference of the three feature nodes, and extract the first initial deformation base point and the second initial deformation base point; construct the first inward convergent deformation guiding surface and the second outward expansion deformation guiding surface based on the first initial deformation base point and the second initial deformation base point.

[0011] Step 4: Calculate the dynamic boundary coupling coefficient based on the complementary normal overlap between the first inward convergent deformation guiding surface and the second outward expansion deformation guiding surface. Use the dynamic boundary coupling coefficient to perform temporal smoothing filtering on the initial node deformation vector to obtain the temporal deformation state evolution trajectory tensor.

[0012] Step 5: Perform multi-objective non-convex optimization and dynamic potential energy balance iteration on the temporal deformation state evolution trajectory tensor to obtain a multi-constraint dynamic equilibrium steady-state deformation field; perform spatial inverse mapping reconstruction and verification on the multi-constraint dynamic equilibrium steady-state deformation field to obtain a complete target region contour set without holes.

[0013] Secondly, a radiotherapy target volume adaptive delineation system based on CBCT images includes:

[0014] The extraction module is used to acquire reference positioning CT images and current fractional CBCT images, perform multi-scale spatial frequency domain alignment and feature extraction processing, and obtain the initial spatiotemporal registration feature manifold;

[0015] The mapping module is used to perform high-dimensional topological reconstruction and feature affinity calculation on the initial spatiotemporal registration feature manifold, construct a spatial interaction intensity mapping model and transform the spatial interaction intensity mapping model into a local structural potential energy gradient distribution, and obtain the local heterogeneous potential field gradient distribution matrix.

[0016] The module is used to perform curvature tensor analytical operations and spatiotemporal state projection on the gradient distribution matrix of the local heterogeneous potential field, extract local curvature eigenvalues, and solve for the extreme feature nodes of the principal deformation axis, the inflection point feature nodes of the secondary deformation axis, and the feature nodes of the topological transition region. Based on the difference of the projected coordinates of the three feature nodes, the initial node deformation vector is derived, and the first initial deformation base point and the second initial deformation base point are extracted. Based on the first initial deformation base point and the second initial deformation base point, the first inward convergent deformation guiding surface and the second outward expansion deformation guiding surface are constructed.

[0017] The calculation module is used to calculate the dynamic boundary coupling coefficient based on the complementary normal overlap between the first inward convergent deformation guiding surface and the second outward expansion deformation guiding surface, and to use the dynamic boundary coupling coefficient to perform temporal smoothing filtering on the initial node deformation vector to obtain the temporal deformation state evolution trajectory tensor.

[0018] The reconstruction module is used to perform multi-objective non-convex optimization and dynamic potential energy balance iteration on the temporal deformation state evolution trajectory tensor to obtain a multi-constraint dynamic equilibrium steady-state deformation field; and to perform spatial inverse mapping reconstruction and verification on the multi-constraint dynamic equilibrium steady-state deformation field to obtain a complete target region contour set without holes.

[0019] Thirdly, a computing device including a memory and a processor;

[0020] The memory stores one or more computer programs, the one or more computer programs including instructions; when the instructions are executed by the processor, the computing device performs the method as described in the first aspect.

[0021] Fourthly, a computer-readable storage medium for storing a computer program for performing the method as described in the first aspect.

[0022] The above-described solution of the present invention has at least the following beneficial effects:

[0023] This invention employs multi-scale spatial frequency domain separation and phase normalization processing to effectively suppress inherent scattering artifacts and high-frequency noise interference in CBCT images, compensating for the imaging defects of insufficient soft tissue contrast in CBCT. Combined with local mutual information entropy and gradient-weighted matching mechanisms, it enhances the feature correlation stability between CT and CBCT cross-modal images, improving the basic accuracy and robustness of registration under low-quality clinical images compared to a single grayscale registration method. Through high-dimensional topological reconstruction, curvature tensor analysis, and local potential field gradient modeling, it abandons the uniform smooth constraint mode and relies on spatial interaction strength and local heterogeneous potential fields to quantify the deformation resistance and structural differences in the interface regions of soft and hard tissues, bone, and air-filled cavities. By dividing multiple types of deformation feature nodes and constructing anisotropic complementary deformation guiding surfaces, it can adaptively adapt to nonlinear anatomical deformations such as tumor regression, organ filling changes, and tissue compression and expansion, effectively suppressing non-physical problems such as local abrupt changes in the deformation field and vector direction disorder. The system incorporates topological connectivity detection, boundary self-intersection correction, and adaptive hole closure mechanisms to specifically address the inherent problems of cracks, holes, and edge defects in the target area contour during topological transition regions. During deformation reconstruction, it takes into account both global topological constraints and local curvature morphology matching, ensuring that the updated target area contour closely matches the current anatomical structure, avoiding abnormal contour deformation, eliminating the need for clinicians to manually repair contour defects, and reducing the cost of manual intervention. Attached Figure Description

[0024] Figure 1 This is a schematic flowchart of a radiotherapy target area adaptive delineation method based on CBCT images provided in an embodiment of the present invention.

[0025] Figure 2 This is a schematic diagram of a radiotherapy target area adaptive delineation system based on CBCT images provided in an embodiment of the present invention. Detailed Implementation

[0026] Exemplary embodiments of the present disclosure will now be described in more detail with reference to the accompanying drawings. While exemplary embodiments of the present disclosure are shown in the drawings, it should be understood that the present disclosure may be implemented in various forms and should not be limited to the embodiments set forth herein. Rather, these embodiments are provided so that this disclosure will be thorough and complete, and will fully convey the scope of the disclosure to those skilled in the art.

[0027] like Figure 1 As shown, embodiments of the present invention propose a method for adaptive delineation of radiotherapy target areas based on CBCT images, the method comprising the following steps:

[0028] Step 1: Acquire the reference positioning CT image and the current fractional CBCT image, perform multi-scale spatial frequency domain alignment and feature extraction processing to obtain the initial spatiotemporal registration feature manifold;

[0029] Step 2: Perform high-dimensional topological reconstruction and feature affinity calculation on the initial spatiotemporal registration feature manifold, construct a spatial interaction intensity mapping model and transform the spatial interaction intensity mapping model into a local structural potential energy gradient distribution to obtain the local heterogeneous potential field gradient distribution matrix;

[0030] Step 3: Perform curvature tensor analytical operations and spatiotemporal state projection on the local heterogeneous potential field gradient distribution matrix, extract local curvature eigenvalues, and solve for the extreme value feature nodes of the principal deformation axis, the inflection point feature nodes of the secondary deformation axis, and the topological transition region feature nodes; derive the initial node deformation vectors based on the projection coordinate difference of the three feature nodes, and extract the first initial deformation base point and the second initial deformation base point; construct the first inward convergent deformation guiding surface and the second outward expansion deformation guiding surface based on the first initial deformation base point and the second initial deformation base point.

[0031] Step 4: Calculate the dynamic boundary coupling coefficient based on the complementary normal overlap between the first inward convergent deformation guiding surface and the second outward expansion deformation guiding surface. Use the dynamic boundary coupling coefficient to perform temporal smoothing filtering on the initial node deformation vector to obtain the temporal deformation state evolution trajectory tensor.

[0032] Step 5: Perform multi-objective non-convex optimization and dynamic potential energy balance iteration on the temporal deformation state evolution trajectory tensor to obtain a multi-constraint dynamic equilibrium steady-state deformation field; perform spatial inverse mapping reconstruction and verification on the multi-constraint dynamic equilibrium steady-state deformation field to obtain a complete target region contour set without holes.

[0033] In this embodiment of the invention, a multi-scale spatial frequency domain separation and phase normalization processing method is adopted, which can effectively suppress the inherent scattering artifacts and high-frequency noise interference of CBCT images and compensate for the imaging defects of insufficient soft tissue contrast in CBCT. Combined with local mutual information entropy and gradient weighted matching mechanism, the feature correlation stability of CT and CBCT cross-modal images is enhanced, which improves the basic accuracy and robustness of registration under low-quality clinical images compared with the single gray-scale registration method. Through high-dimensional topological reconstruction, curvature tensor analysis and local potential field gradient modeling, the uniform smooth constraint mode is abandoned. Relying on the spatial interaction strength and local heterogeneous potential field, the deformation resistance and structural differences of the interface regions of soft and hard tissues, bone, and air-containing cavities are quantified. By dividing multiple types of deformation feature nodes and constructing anisotropic complementary deformation guiding surfaces, it can adaptively adapt to nonlinear anatomical deformations such as tumor regression, organ filling changes, and tissue compression expansion, effectively suppressing non-physical problems such as local abrupt changes in deformation field and vector direction disorder. The system incorporates topological connectivity detection, boundary self-intersection correction, and adaptive hole closure mechanisms to specifically address the inherent problems of cracks, holes, and edge defects in the target area contour during topological transition regions. During deformation reconstruction, it takes into account both global topological constraints and local curvature morphology matching, ensuring that the updated target area contour closely matches the current anatomical structure, avoiding abnormal contour deformation, eliminating the need for clinicians to manually repair contour defects, and reducing the cost of manual intervention.

[0034] In a preferred embodiment of the present invention, step 1 includes:

[0035] Step 100: Acquire the reference positioning CT image and the current fractionated CBCT image; perform a three-dimensional discrete frequency domain orthogonal transform on the reference positioning CT image and the current fractionated CBCT image to separate the low-frequency energy baseband and the high-frequency detail band, specifically including:

[0036] The reference positioning CT image acquired by spiral CT during the radiotherapy positioning phase and the current fraction CBCT image acquired in real time by the current treatment fraction CBCT device are read. Both images are preprocessed: the voxel spacing is standardized to 1mm×1mm×1mm, blank and invalid areas at the edges are cropped, and the grayscale values ​​of both CT and CBCT images are normalized to [value missing]. HU ensures that both are three-dimensional volumetric data of the same patient and the same anatomical region. After preprocessing, the voxel matrix of the reference localization CT image is denoted as... The current voxel matrix of CBCT images is ,in These are the positive integer coordinates of the voxel in the transverse, coronal, and sagittal directions, respectively, corresponding to the voxel's index position in the image, used to determine the spatial location of each voxel. To achieve accurate frequency domain alignment and noise separation, a three-dimensional discrete Fourier transform (3DDFT) is performed on the two preprocessed images to convert the spatial domain grayscale signal into a frequency domain signal, separating the overall structure from local details and noise information. The specific transformation formula is as follows: ; ;

[0037] In the formula, It is the frequency domain signal of the reference positioning CT image. It is the frequency domain signal of the current fractional CBCT images. For spatial domain coordinates The frequency domain coordinates that correspond one-to-one have a range of values. These correspond to the frequency components in the transverse, coronal, and sagittal directions, respectively. The larger the value, the higher the frequency in the corresponding direction, and the finer the image details it carries. The total number of voxels in each dimension of the image after preprocessing. The total number of voxels in all three dimensions is consistent; It is the imaginary unit.

[0038] After the DFT is completed, frequency bands are divided according to the energy percentage threshold. First, the energy amplitude of each component in the frequency domain is calculated. ( (These are the real and imaginary parts of the frequency domain components, respectively); then calculate the total energy. The cumulative energy percentage is calculated by sorting the frequency domain coordinates from low to high. The low-frequency components with a cumulative energy percentage of ≤90% are classified as the low-frequency energy baseband (characterizing the overall anatomical structure, grayscale trend and large-scale morphology, which is the core basis for cross-modal registration). The high-frequency components with a cumulative energy percentage of ≤10% are classified as the high-frequency detail band (carrying tissue edges, textures, Gaussian noise and CBCT scattering artifacts), thus completing the frequency band separation.

[0039] Step 101: Perform cross-power spectrum phase normalization on the high-frequency detail bands to extract the cross-modal frequency domain phase difference distribution; perform inverse frequency domain spatial mapping transformation on the cross-modal frequency domain phase difference distribution to obtain the spatial cross-correlation response kernel; based on the spatial cross-correlation response kernel, perform sliding window local mutual information entropy calculation and local gradient magnitude matching weighted aggregation to obtain the cross-modal spatial correlation distribution map, specifically including:

[0040] Based on the high-frequency detail bands separated in step 100, and utilizing the consistency of the phase characteristics of CT and CBCT structures to avoid interference from grayscale amplitude differences, the cross-power spectrum of the high-frequency components of the two is first calculated to capture the phase correlation, as shown in the following formula: ;

[0041] The cross-power spectrum represents the high-frequency components of the reference positioning CT and the current sub-CBCT, used to capture the phase correlation between the high-frequency signals of the two. The high-frequency component frequency domain matrix of the reference CT image is the frequency domain data corresponding to the high-frequency detail frequency band of the CT obtained in step 100. The matrix dimension is consistent with the CT frequency domain matrix. The complex conjugate of the frequency domain matrix of the high-frequency components of the current fractional CBCT image is given by the following formula: It can preserve the phase difference between the high-frequency components of CT and CBCT, while eliminating the interference caused by the difference in grayscale amplitude between the two. These are frequency domain coordinates, consistent with the meaning in step 100; , These are the real and imaginary parts of the frequency domain component, respectively, with the same meaning as in step 100.

[0042] To highlight phase characteristics and reduce amplitude fluctuations, the cross-power spectrum is phase-normalized using the following formula: ;

[0043] In the formula The normalized cross-modal frequency domain phase difference distribution retains only the phase characteristics. The amplitude of the cross-power spectrum is calculated in the same way as the amplitude of the frequency domain component energy in step 100. After normalization, the amplitude of the cross-modal frequency domain phase difference distribution is unified to 1, which can reflect the spatial correspondence of high-frequency details in the two images. The smaller the phase difference, the higher the structural matching degree at the corresponding position. A three-dimensional inverse discrete Fourier transform (3DIDFT) is performed on the normalized cross-modal frequency domain phase difference distribution to map the frequency domain phase difference to the spatial domain, obtaining the spatial cross-correlation response kernel characterizing the local structural matching degree. The formula is as follows: ;

[0044] The magnitude of the cross-correlation response kernel value in the spatial domain directly reflects the coordinate relationship between the reference positioning CT and the current fractionated CBCT images. The degree of local structural matching at a given location; the larger the value, the higher the degree of matching. Based on this, a 3×3×3 three-dimensional sliding window (step size 1 voxel) is used to traverse the entire image, and two indicators are calculated simultaneously to quantify the local matching degree, one of which is the local mutual information entropy. This characterizes the dependency between CT and CBCT grayscale values ​​within the window; a larger value indicates a higher degree of matching. In the formula CT grayscale values ​​within the window With CBCT grayscale values The joint probability is obtained by combining gray levels within a statistical window. The number of occurrences is obtained by dividing by the total number of voxels in the window (27). The first is the normalized grayscale value, ranging from [-1000, 1000] HU, which is discretized into 256 levels during calculation. The second is the local gradient magnitude. By calculating the partial derivatives of voxels within the window in three directions, the three-dimensional spatial gradient magnitude is solved, thereby characterizing the clarity of tissue edges and the intensity of structural changes. The specific formula is as follows: ;

[0045] In the formula These are the partial derivatives of the voxels within the window in three directions. The gradient magnitude represents the grayscale value of the voxel within the window. The larger the gradient magnitude, the clearer the tissue edge and the more dramatic the structural changes.

[0046] The three indicators mentioned above are cascaded and fused to calculate the comprehensive matching confidence score for each sliding window position. Specifically, for voxel coordinates... The current sliding window centered on spatial cross-correlation response values As the basic matching response, this value reflects the degree of alignment of the underlying structure at that location after frequency domain phase normalization and inverse transform. Simultaneously, it incorporates the local mutual information entropy calculated for the current window. With local gradient magnitude Both are subjected to deviation standardization to eliminate the dimensional differences between the two and unify them to a comparable interval, resulting in the normalized local mutual information entropy. With normalized local gradient magnitude And all are in the interval [0,1]. The standardized values ​​are... and We perform weighted summation to construct a comprehensive confidence adjustment factor. : ;

[0047] Among them, the weighting coefficient =0.5, used to balance the contribution weights of gray-level statistical correlation and marginal structural significance to the overall confidence score, ensuring that both have equal influence on the final confidence score evaluation. This overall confidence score adjustment factor... As a weighting coefficient, the basic cross-correlation response value Make corrections to obtain the current window center position. The final overall match confidence score at the location: ;

[0048] It provides a basic matched response based on frequency domain phase, while Then, based on the correlation of grayscale distribution ( ) and salience of edge structure ( The confidence-weighted adjustment of the base response is performed using two dimensions. When the grayscale statistical correlation between two images within the window is high and their edge structures are clear... A score close to 1 indicates that the overall confidence score is fully preserved; conversely, a score lower than 1 indicates poor matching quality in any dimension. The corresponding decrease in confidence score effectively suppresses the overall confidence score, thereby reducing the risk of mismatch caused by CBCT artifacts or insufficient soft tissue contrast.

[0049] Following the above calculation process, the entire image is traversed using a 3×3×3 three-dimensional sliding window with a step size of 1 voxel, and the calculation is performed at each voxel position. Calculate the overall matching confidence score for each one. After traversal, the confidence scores of all voxel locations are calculated according to their three-dimensional coordinates. The axis from left to right The axis from front to back, The axes are arranged in an orderly manner from lightest to darkest to construct a cross-modal spatial correlation distribution map that is completely consistent with the spatial dimensions of the original image. The values ​​at each voxel position in this distribution map comprehensively reflect the joint verification results of phase matching degree, gray-level statistical correlation, and edge structure saliency.

[0050] Step 102: Based on the cross-modal spatial correlation distribution map, perform global peak extremum search and sub-pixel level offset interpolation to solve the initial rigid spatial transformation parameters and construct the initial rigid transformation alignment field; based on the initial rigid transformation alignment field, perform voxel-level trilinear interpolation resampling and alignment on the current CBCT image to obtain aligned image data, specifically including:

[0051] Traverse the generated steps 101 The cross-correlation response values ​​of each voxel were compared, and the response value of each voxel was recorded. By comparing them one by one, the voxel with the largest response value was selected, and its coordinates were located. This coordinate is the optimal rigid alignment center between the reference positioning CT and the current fractional CBCT, with the highest matching degree, and serves as the reference point for rigid registration.

[0052] To achieve sub-pixel level registration accuracy, Using this as the core, a 5×5×5 fitting window is set, with the window range covering... The voxel region is defined; using cubic polynomial interpolation, the cross-correlation response values ​​and corresponding coordinates of all voxels within the window are first collected. The coordinates and response values ​​are then substituted into the cubic polynomial formula for fitting. The fitting formula is: ; ( to (fit coefficients) Given a constant, a continuous function of the response value as a function of spatial coordinates is obtained. The first derivative of this continuous function is taken, and the result is set to 0. Solving this equation yields the coordinates of the extreme point, which is the sub-pixel-level optimal alignment coordinate. The sub-pixel-level optimal alignment coordinates and the reference point are then calculated. The difference in the three directions yields the sub-pixel precision translation amount. Simultaneously, the second derivative of the fitted continuous function is calculated to obtain the rotation amount. This reflects the spatial rotation angle of the CBCT relative to the CT. Combining the aforementioned translation and rotation, initial rigid transformation parameters are constructed to form an alignment field, specifically through a 3×3 rotation matrix. With three-dimensional translation vector Implementation, in which It is the product of the rotation matrices in three directions. Composed of translation amounts, substituted into the rigid transformation solution formula This allows us to obtain the initial rigid space transformation parameters, thereby achieving rigid alignment of CBCT to CT. The original three-dimensional coordinate vector of the voxels in the current fractional CBCT image, compared with that in step 100 The meanings are consistent; both refer to positive integer coordinates. Determined by the amount of rotation, where around Axis rotation matrix , around axis, (The same applies to the axis rotation matrix). Multiplying the three matrices together yields the overall rotation matrix. This transformation allows each voxel in the CBCT to be precisely mapped to its corresponding spatial position in the CT.

[0053] For the non-integer target voxel positions obtained after rigid transformation, voxel-level trilinear interpolation resampling is performed: First, eight adjacent actual voxels around the target position are determined, which are eight coordinates formed by rounding down and up in the x, y, and z directions of the target, and the gray values ​​of each voxel are read. Then, weights are assigned according to the normalized distance between the target position and each adjacent voxel. The weight is the product of the corresponding values ​​in the x, y, and z directions, and all weights are between 0 and 1 with a sum of 1. The closer the distance, the greater the weight. Finally, each weight is multiplied by the corresponding voxel gray value, and the sum is obtained to obtain the gray value of the target voxel. After resampling, aligned image data with the reference positioning CT is obtained, whose spatial coordinates are completely consistent with the CT, with grayscale fidelity and no distortion.

[0054] Step 103: Perform nonlocal mean denoising and edge-aware convolution operations on the aligned image data to extract a multi-scale texture feature vector set; map the multi-scale texture feature vector set to a high-dimensional tensor space, perform Riemannian manifold dimensionality reduction embedding processing, and project the high-dimensional features to a low-dimensional continuous metric space to obtain the initial spatiotemporal registration feature manifold, specifically including:

[0055] Nonlocal mean denoising (NLM) is performed on the aligned image data obtained in step 102 to suppress scattering artifacts and Gaussian noise while preserving tissue edges and texture details. An 11×11×11 search window and a 3×3×3 similarity block window are set with each voxel as the center. All voxels within the search window are traversed, and the Euclidean distance between the current voxel's similarity block and other voxel similarity blocks within the search window is calculated. Gaussian weights are then assigned based on the Euclidean distance. In the formula For the current voxel With candidate voxels Similar block weights, For Euclidean distance, =1.0 represents the Gaussian kernel bandwidth; the denoised grayscale value of the current voxel is calculated using a weighted average. ,in =1, ensuring the grayscale value amplitude is consistent with the original image. This is for aligned image data in candidate voxels The grayscale value at the specified location is obtained by rigid registration and trilinear interpolation resampling in step 102. After denoising the entire image voxels, the denoised image is obtained.

[0056] Edge-aware convolution operations were performed on the denoised images to construct three sets of 3D convolution kernels with the following parameters: 3×3×3 small-scale convolution kernels (focusing on extracting subtle texture features within soft tissues), 5×5×5 medium-scale convolution kernels (capturing the boundary features between the target area and normal tissue), and 7×7×7 large-scale convolution kernels (capturing the overall morphology and structural changes of organs). All three sets of convolution kernels were subject to x, y, and z-directional perceptual constraints (ensuring accurate capture of features and edges in three dimensions). Convolution was performed region-by-region along the entire image domain to generate three corresponding feature maps. Batch normalization was performed on each type of feature map. ,in This represents the average gray level of all voxels in the current feature map. For grayscale variance, (To avoid calculation errors when the variance is 0) These are learnable parameters; after batch normalization, the ReLU activation function is used to enhance the nonlinear representation of features, as shown in the formula. The operation logic is to retain all non-negative feature values, set negative feature values ​​to 0, filter effective features, suppress invalid features, and finally integrate the three sets of feature maps to generate a 256-dimensional multi-scale texture feature vector set. Each voxel corresponds to a 256-dimensional feature vector, covering multi-dimensional information such as grayscale, edge, and texture.

[0057] A 256-dimensional feature vector set is mapped to a high-dimensional tensor space, constructing a high-dimensional composite feature tensor with dimensions of image height × width × number of layers × 256 (height, width, and number of layers correspond to the image's x, y, and z 3D dimensions, and 256 is the feature dimension). Then, the 256-dimensional feature vector of each voxel is expanded to 1024 dimensions, fusing voxel 3D spatial coordinates, voxel grayscale values, and voxel gradient magnitudes during the expansion process. A Riemannian manifold dimensionality reduction algorithm is used with parameter kw=20, selecting the 20 most similar neighboring nodes for each voxel, with dimension dy=30, effectively reducing the high-dimensional features to 30 dimensions. First, the Euclidean distance between voxels is calculated, and the 20 neighboring nodes with the smallest distance are selected. A local linear expression is then fitted, first defining the formula for the local linear expression as follows: ,in The 1024-dimensional high-dimensional feature vector of the target voxel. The set of 20 neighboring nodes of the target voxel. For the first The 1024-dimensional feature vector of the neighboring nodes, For the first Neighboring nodes for the target voxel linear weights, The fitting error (must satisfy) Solving for linear weights The core of the least squares method is to minimize the sum of squares of the fitting errors, that is, to minimize the objective function. At the same time, a weight normalization constraint is added. By taking the partial derivative of the objective function and setting the partial derivative to 0, the linear weights are obtained, and the high-dimensional feature vector of the target voxel is accurately expressed as a linear combination of the feature vectors of 20 neighboring nodes. Then, nonlinear dimensionality reduction is performed to map the fitted high-dimensional feature vector to a high-dimensional Riemannian manifold space. The linear correlation constraint of the high-dimensional features is broken by manifold embedding technology. Relying on the bending characteristics of the Riemannian manifold, the high-dimensional features are nonlinearly projected onto a 30-dimensional low-dimensional continuous metric space, which both suppresses the redundancy of high-dimensional features and preserves the spatial topology and anatomical structure association between voxels. During this process, redundant features with variance <0.01 are removed, highly coupled features with correlation coefficient >0.9 are merged, and time labels are fused in stages to finally generate an initial spatiotemporal registration feature manifold with dimensions of image height × width × number of layers × 30.

[0058] This embodiment achieves high- and low-frequency band layering through three-dimensional discrete frequency domain orthogonal transformation. Combined with cross-power spectrum phase normalization and cross-modal phase difference analysis, it can specifically overcome the imaging defects of CBCT images, such as prominent scattering artifacts, high noise levels, and low soft tissue contrast. It also overcomes the limitations of single-spatial-domain registration being susceptible to noise interference, improving the feature matching stability of two heterogeneous modal images: reference positioning CT and fractional CBCT. Relying on the sliding window local mutual information entropy calculation and gradient weighted aggregation, and sub-pixel level interpolation rigid registration mechanism, it achieves high-precision solution for coarse registration of cross-modal images. Combined with trilinear interpolation resampling and non-local mean denoising, it eliminates image coordinate offset and gray-level misalignment while completely preserving key anatomical details such as organ boundaries and tissue junctions. By employing a feature construction method that combines multi-scale edge-aware convolutional feature extraction with Riemannian manifold dimensionality reduction embedding, the shallow grayscale features of images are upgraded to high-dimensional continuous spatiotemporal registration feature manifolds. This enhances the spatiotemporal correlation expression of anatomical structures, can adapt to subtle changes in anatomical morphology between different phases, and solves the problems of limited low-dimensional feature expression and insufficient nonlinear deformation representation capabilities.

[0059] In a preferred embodiment of the present invention, step 2 includes:

[0060] Step 200: Extract the voxel neighborhood spatial topological connectivity relationships in the initial spatiotemporal registration feature manifold to construct a K-nearest neighbor feature association graph; calculate the texture gradient cosine similarity and spatial Euclidean distance decay weights between adjacent nodes in the K-nearest neighbor feature association graph, and obtain the feature affinity weight matrix through Gaussian kernel function nonlinear mapping, specifically including:

[0061] Based on the initial spatiotemporal registration feature manifold, a global traversal is performed on all voxel nodes within the manifold. Using each target voxel node as the center, and with a preset nearest neighbor count K=20, the 20 voxel nodes spatially closest to the target voxel node are selected as nearest neighbors in a 3D spatial coordinate system using Euclidean distance. The spatial topological connections between the target voxel node and each nearest neighbor node are extracted, clarifying the connection paths, adjacency relationships, and spatial positional constraints between nodes. This constructs a K-nearest neighbor feature association graph with a complete spatial topological structure from the discrete voxel nodes in the initial spatiotemporal registration feature manifold. For any pair of adjacent voxel nodes in the K-nearest neighbor feature association graph, the texture gradient feature vectors corresponding to the two nodes are extracted, and the cosine similarity of the texture gradients of the two feature vectors is calculated. The closer this value is to 1, the higher the consistency of the texture features of the adjacent nodes. Simultaneously, the Euclidean distance between two adjacent voxel nodes in 3D space is calculated, and the Euclidean distance is converted into a spatial Euclidean distance decay weight using an exponential decay function, consistent with the Gaussian weight formula in step 103. The smaller the distance between voxel nodes, the larger the weight value; the larger the distance, the smaller the weight value. The texture gradient cosine similarity and the spatial Euclidean distance decay weight are weighted and fused. Using the fused features as input, a Gaussian kernel function is substituted to complete the nonlinear mapping transformation, resulting in smooth and stable weight values. The calculation results of all adjacent node pairs are arranged in an orderly manner, row by row, column by column, and layer by layer, according to the 3D voxel coordinates x from left to right, y from front to back, and z from light to dark. Finally, a feature affinity weight matrix that perfectly matches the spatial dimension of the initial spatiotemporal registration feature manifold is generated. Each element in the matrix corresponds to a set of feature association strengths between adjacent voxel nodes. The feature association strength is expressed using continuous floating-point values ​​in the interval [0,1]. The closer the value is to 1, the stronger the feature association between voxel nodes; the closer the value is to 0, the weaker the association.

[0062] Step 201: Using the feature affinity weight matrix as the driving parameter, perform a coupling strength field fitting operation based on local kernel density estimation, mapping the feature affinity weights to a multidimensional force transfer potential function between nodes, iteratively updating the parameters of the multidimensional force transfer potential function until the local energy distribution converges, and constructing a spatial interaction strength mapping model, specifically including:

[0063] Using the feature affinity weight matrix calculated in step 200 as the core driving parameter, the coupling strength field fitting operation based on local kernel density estimation is initiated. Taking each voxel node in the K-nearest neighbor feature association graph as the computation center, kernel density estimation is performed within the three-dimensional local neighborhood to construct the basic coupling strength field. The kernel density estimation formula is: ;

[0064] in, It is the local kernel density estimate at the current voxel location. These are the spatial coordinates of the current voxel node. The coordinates of the voxel nodes in the neighborhood. This represents the total number of valid nodes in the neighborhood. K is the kernel width parameter, with a value of 1.0. ) is a three-dimensional Gaussian kernel function. It is an index of voxels within the neighborhood.

[0065] The feature affinity weights obtained in step 200 As a weighting factor, it is fused with the kernel density estimation results to construct a multidimensional force transfer potential function, expressed as: ;

[0066] in, It is a multidimensional force transmission potential function, which quantitatively characterizes the spatial interaction between adjacent voxel nodes in three-dimensional vector form. The vector amplitude directly represents the interaction strength through the magnitude of the potential function. The larger the magnitude, the stronger the interaction between adjacent voxel nodes, and the smaller the magnitude, the weaker the interaction. The vector direction is determined by the spatial unit vector pointing from the current voxel node to the neighboring voxel nodes, and is used to uniquely characterize the spatial transmission direction of the interaction. The vector magnitude attenuation rate is characterized by the distance-related exponential attenuation term. The larger the attenuation term, the higher the action transmission efficiency, and the smaller the attenuation term, the lower the action transmission efficiency, thus achieving simultaneous quantification of the three indicators: amplitude, direction, and transmission efficiency.

[0067] Set the initial parameters and iteration step size. =0.01, energy convergence threshold And construct a local energy convergence function: ;

[0068] in It is the sum of the total energy in both the global and local domains. The number of neighboring nodes is 20. The potential function value for the force transfer obtained in the previous iteration is used. The iterative update process then begins, where the potential function parameters are updated along the energy decreasing direction in each iteration, and the current global and local energy is calculated. Calculate the energy difference between two adjacent iterations. ,in It is the first The global energy of the next iteration, when the difference is less than the convergence threshold. When the local energy distribution reaches a stable convergence state, the iteration stops. Through the above fitting and iteration process, a spatial interaction intensity mapping model that can completely describe the spatial interaction law of all element nodes within the initial spatiotemporal registration feature manifold is finally constructed.

[0069] Step 202: Perform point-by-point potential field integration on the voxel spatial coordinates of the initial spatiotemporal registration feature manifold using the spatial interaction intensity mapping model to obtain the scalar field of local structural potential energy distribution; perform discrete difference gradient operator convolution on the scalar field of local structural potential energy distribution to transform the scalar potential energy field into a spatial direction gradient vector field, obtaining the local structural potential energy gradient distribution; perform voxel gridded discrete rearrangement and matrix dimension alignment on the local structural potential energy gradient distribution to obtain the local heterogeneous potential field gradient distribution matrix, specifically including:

[0070] Based on the spatial interaction intensity mapping model constructed in step 201, the spatial coordinates of each voxel within the initial spatiotemporal registration feature manifold are traversed. Perform point-by-point potential field integration for each voxel position, with the integration interval limited to the area centered on the current voxel. 3D neighborhood The integral formula for the potential field is: ,in It is a spatial interaction intensity mapping model in spatial location The intensity value of the effect at that location; It is centered on the current voxel A three-dimensional integral neighborhood limits the range of integral calculations. It is a three-dimensional spatial element used to perform spatial domain integration operations; It is the local structural potential energy corresponding to the current voxel position. The local structural potential energy corresponding to the current voxel position is characterized by the value obtained by the potential field integration operation. The larger the value, the stronger the spatial interaction of the current voxel with the surrounding voxels, and the greater the total potential energy. The smaller the value, the weaker the spatial interaction of the current voxel with the surrounding voxels, and the smaller the total potential energy. This quantitatively characterizes the potential energy accumulation level of this position in the local spatial structure.

[0071] A discrete difference gradient operator with dimensions of 3×3×3 is used to perform spatial differentiation calculations on the scalar field of the local structural potential energy distribution, transforming the scalar potential energy field containing only numerical information into a spatial gradient vector field that simultaneously contains direction and magnitude information. The three-dimensional gradient calculation formula is as follows: ;

[0072] in, It is the transpose sign. The partial differentials in each coordinate axis direction are calculated using the central difference scheme. Taking the x-direction as an example, the calculation formula is: ;

[0073] In the formula The three-dimensional gradient vector at the current voxel position is used to characterize the spatial rate and direction of change of the local structural potential energy at that point. These are the partial differential results of the local structural potential energy in the three orthogonal directions of x, y, and z, respectively; The image voxel spacing is used to convert the voxel coordinate difference into the differential length of the real physical space; and These represent the local structural potential energy values ​​of the current voxel's preceding and following voxels in the x-direction. The local structural potential energy gradient distribution is obtained through spatial differentiation.

[0074] According to the standard voxel grid coordinate rules of three-dimensional radiotherapy images, the local structural potential energy gradient distribution is discretized and rearranged using voxel grids. The gradient vector data is then arranged in a regular three-dimensional coordinate order: x-axis from left to right, y-axis from front to back, and z-axis from light to dark. After the rearrangement is completed, a matrix dimension alignment operation is performed to ensure that the dimension and size of the gradient distribution data are completely consistent with the spatial grid structure and the initial spatiotemporal registration feature manifold. Finally, a standardized local heterogeneous potential field gradient distribution matrix is ​​generated.

[0075] In this embodiment, pointwise potential field integration can transform abstract spatial interactions into an intuitive scalar field of local structural potential energy distribution. After being converted into a gradient vector field by gradient operator convolution, the deformation resistance and structural change trend of the local region can be accurately quantified. The resulting local heterogeneous potential field gradient distribution matrix can provide standardized and highly reliable data support for subsequent deformation feature node extraction and deformation guidance surface construction, effectively suppressing the interference of CBCT image defects on subsequent calculations.

[0076] In a preferred embodiment of the present invention, step 3 includes:

[0077] Step 300a: Perform spatial second-order differential gradient analytical operation on the gradient distribution matrix of the local heterogeneous potential field to extract local curvature feature values ​​characterizing the degree of spatial curvature of the local potential field and the deformation resistance gradient; construct a spatial coordinate mapping grid based on the local curvature feature values, and perform an extremum search along the direction of potential energy gradient in the spatial coordinate mapping grid to locate the center position of the voxel region where the local potential field gradient magnitude is greater than the preset gradient aggregation threshold and the spatial connectivity satisfies the eight-neighbor continuity condition. This location is denoted as the extremum feature node of the principal deformation axis, specifically including:

[0078] Using the local heterogeneous potential gradient distribution matrix generated in step 202 as input data, perform spatial second-order differential gradient analytical operations on this matrix. First, perform analysis on the scalar field function corresponding to the local heterogeneous potential gradient distribution matrix. The second-order partial differentials are calculated in three-dimensional space along the three orthogonal directions x, y, and z using a central difference scheme. Based on the results of the second-order partial differentials, a three-dimensional curvature tensor is constructed. This tensor is a 3×3 symmetric matrix, specifically in the form of: ;

[0079] in, Let x, y, and z represent the second-order partial derivatives of the scalar field function along the principal axes, respectively. Represents the cross-second partial derivative of the scalar field function between the two principal axes; for the three-dimensional curvature tensor Perform eigenvalue decomposition to solve the equation eigenvalues ,in For feature vectors, The maximum absolute value among the three eigenvalues ​​is taken as the local curvature eigenvalue of the current voxel. The local curvature eigenvalue is used to quantitatively characterize the degree of curvature of the local potential field and the tissue deformation resistance gradient. The larger the absolute value of the local curvature eigenvalue, the higher the degree of curvature of the local potential field, and the greater the tissue deformation resistance in the corresponding region; the smaller the absolute value of the local curvature eigenvalue, the lower the degree of curvature of the local potential field, and the smaller the tissue deformation resistance in the corresponding region.

[0080] Based on local curvature eigenvalues, a spatial coordinate mapping grid is constructed that is completely consistent with the voxel grid size and spatial coordinate system of CBCT images. Within this grid, a global voxel extremum search is performed along the direction of potential gradient ascent, and the local potential field gradient magnitude is calculated for each voxel. = The gradient magnitude is compared with a preset gradient aggregation threshold. Compare, retain The voxels are selected. For each voxel, an eight-neighbor spatial connectivity test is performed. The current voxel, in 3D space, has at least six of its eight neighboring voxels (front, back, left, right, top, bottom, left-front-top, right-back-bottom) with gradient magnitudes greater than a threshold. A valid voxel is defined as one that has no breaks or isolated gaps between adjacent voxels, thus satisfying the eight-neighborhood continuity condition. All voxels that meet the condition are aggregated into a connected region, and the arithmetic mean of the x, y, and z coordinates of all voxels within this region is calculated. The spatial position corresponding to this average coordinate is defined as the extreme value feature node of the principal deformation axis.

[0081] Step 301a: Tracing contour lines along the boundary where the potential energy gradient sign of the spatial coordinate mapping grid is reversed, locating the grid boundary position where the spatial coordinates are at the interface of abrupt change in deformation resistance and where the positive and negative signs of the local curvature characteristic values ​​are continuously reversed, denoted as the secondary deformation axis inflection point feature node, specifically including:

[0082] Within the constructed spatial coordinate mapping mesh, contour tracing is performed along the boundary locations where the potential energy gradient sign reverses. A voxel-by-voxel traversal detection method is used to determine the sign change of the local curvature eigenvalue for each voxel, capturing the mesh boundary locations where the potential energy gradient direction reverses. Based on the local heterogeneous potential field gradient distribution matrix and the distribution characteristics of local curvature eigenvalues, preliminary region division is performed. Soft tissue regions correspond to local curvature eigenvalues ​​with absolute values ​​<0.3 and local potential field gradient amplitudes <0.3. Regions with a value <0.3; regions with a value ≤0.3 for solid tissue (such as muscle and glands) correspond to local curvature eigenvalues ​​with an absolute value <0.6 and local potential gradient amplitudes ≤0.3. Regions with a curvature characteristic value < 0.6; regions containing air cavities correspond to local curvature characteristic values ​​that are negative and have large absolute values, i.e., local curvature characteristic value < -0.6, and local potential field gradient amplitude. Regions with a curvature eigenvalue ≥ 0.6; skeletal regions correspond to local curvature eigenvalues ​​that are positive and have large absolute values, i.e., local curvature eigenvalues ​​≥ 0.6 and local potential field gradient amplitudes. Regions with a value ≥0.6.

[0083] By dividing the area as described above, the spatial distribution range of various tissues is clarified, thereby locating the tissue interfaces between soft tissue and solid tissue, soft tissue and air-containing cavities, and soft tissue and bone. These interfaces, due to significant differences in tissue density and hardness, inevitably exhibit abrupt changes in tissue deformation resistance. Grid boundary voxels whose local curvature feature values ​​at the interface location show a continuous reversal of positive and negative signs are selected. Centered on the currently detected voxel, within its 3×3×3 local neighborhood, there exist voxels with local curvature feature values ​​>0 and voxels with local curvature feature values ​​<0, and the voxels with positive and negative feature values ​​are distributed alternately and continuously (i.e., the curvature feature values ​​of any two adjacent voxels have opposite signs, and there are no three or more voxels with the same sign consecutively). These voxels are the boundary voxels at the interface. All selected boundary voxels of this type are compiled into a boundary voxel set. The arithmetic mean of the x, y, and z coordinates of all voxels in the set is calculated. The spatial position corresponding to the coordinates of the geometric center is defined as the inflection point feature node of the secondary deformation axis. This node is used to mark the key position where the deformation trend of the anatomical structure changes, and can reflect the abrupt change characteristics of the deformation direction at the interface.

[0084] Step 302a: Extract the centroid positions of connected regions whose spatial coordinates are distributed within the projection range of the line connecting the extreme feature nodes of the principal deformation axis and the inflection point feature nodes of the secondary deformation axis, and whose local curvature feature value fluctuation amplitude is lower than a preset smoothing threshold. These are denoted as topological transition region feature nodes, specifically including:

[0085] In the spatial coordinate mapping mesh, using the line connecting the extreme feature nodes of the principal deformation axis and the inflection point feature nodes of the secondary deformation axis as a reference, an orthogonal projection range of this line in the 3D imaging space is constructed, and all voxels within the projection range are extracted. For each voxel, the standard deviation of the local curvature feature value within a 3×3×3 local neighborhood is calculated. This standard deviation is used as the fluctuation amplitude of the local curvature feature value. The fluctuation amplitude is compared with a preset smoothing threshold of 0.15, and voxels with fluctuation amplitudes <0.15 are retained and aggregated into connected regions. Using the local curvature feature value as weight, the spatial centroid coordinates of this connected region are calculated using a weighted average method. The calculation formula is as follows: ;

[0086] in The coordinates of the characteristic nodes in the topological transition region; For the first Local curvature eigenvalues ​​of individual elements; For the first Spatial coordinates of individual elements; This represents the total number of voxels within the connected region. The spatial location corresponding to the centroid coordinates is defined as a topological transition region feature node, used to mark the center location of topological transition regions such as the boundary between soft and hard tissues, and the boundary between tissues and cavities.

[0087] Step 300b involves extracting the spatial projection coordinate sequences of the extreme feature nodes of the primary deformation axis, the inflection point feature nodes of the secondary deformation axis, and the feature nodes of the topological transition zone within consecutive scanning cycles. The temporal position difference of the spatial projection coordinate sequence is calculated, and deformation integrals are derived from the temporal position difference using the time interval parameters to obtain the initial node deformation vector. Specifically, this includes:

[0088] Spatial projection coordinate sequences of the extreme value feature nodes of the primary deformation axis, the inflection point feature nodes of the secondary deformation axis, and the topological transition zone feature nodes were extracted over T consecutive radiotherapy fraction scan cycles. For each type of feature node, the temporal position difference between adjacent fractions was calculated. This yields the coordinate changes in the x, y, and z directions, where... It is the first Spatial projection coordinates within each scan cycle, It is the first Spatial projection coordinates within each scan cycle. The temporal position difference is then compared with the time interval between adjacent scan cycles. Combining these, we perform deformation integral derivation on the coordinate changes, and the integral formula is: ;

[0089] in This is the initial node deformation vector; The total number of integrations is denoted as . The initial node deformation vectors corresponding to the three types of feature nodes are obtained through this integration. These vectors are used to quantitatively characterize the initial deformation direction and displacement magnitude of each feature node during the treatment cycle. The spatial direction of the vector points to the actual deformation offset direction of the anatomical tissue in the corresponding three-dimensional space. The larger the magnitude of the vector, the longer the spatial offset distance of the corresponding feature node and the more obvious the overall displacement change. The smaller the magnitude, the shorter the spatial offset distance of the node and the smoother the displacement change.

[0090] Step 301b involves performing a weighted centroid clustering operation based on the spatial distribution matrix of local curvature eigenvalues ​​to extract the coordinates of the two spatial anchor points with the smallest spatial deviation, denoted as the first initial deformation base point and the second initial deformation base point, respectively. Specifically, this includes:

[0091] Based on the spatial distribution matrix of local curvature eigenvalues, a weighted centroid clustering operation is performed, using the local curvature eigenvalues ​​as cluster weights to perform weighted clustering on the global voxel coordinates. This process iterates through all candidate cluster centers across the entire domain. Calculate the spatial deviation between each candidate center and its surrounding neighborhood voxels. The formula for calculating spatial deviation is: ;

[0092] in, The coordinates are those of the neighboring voxels. Total number of neighboring voxels. Spatial deviation. The smaller the value, the more stable the candidate center. The two candidate centers with the smallest spatial deviation are selected as the final spatial anchor point coordinates, and the first anchor point is defined as the first initial deformation base point. The second anchor point is defined as the second initial deformation base point. The two base points serve as reference positions for tissue contraction deformation and tissue expansion deformation, respectively.

[0093] Step 302b: By performing spatial vertex connection and patch meshing operations on the first initial deformation base point, the extreme feature node of the main deformation axis, and the inflection point feature node of the secondary deformation axis, a first inwardly converging deformation guiding surface is obtained. This first inwardly converging deformation guiding surface provides an inwardly contracting geometric boundary constraint for the rapid reduction of target tissue volume, with its corresponding normal pointing inwards towards the interior of the deformation region. By performing spatial vertex connection and patch meshing operations on the second initial deformation base point, the topological transition region feature node, and the extreme feature node of the main deformation axis, a second outwardly expanding deformation guiding surface is obtained. This second outwardly expanding deformation guiding surface provides an outwardly extending geometric boundary pushing constraint for the edema expansion of tissue surrounding the target region, with its corresponding normal pointing outwards towards the exterior of the deformation region. The first inwardly converging deformation guiding surface and the second outwardly expanding deformation guiding surface constitute an anisotropic deformation guiding pair with opposite spatial normals and complementary functions, specifically including:

[0094] Using the first initial deformation base point (With contraction-type reference anchor points and fixed coordinates as the core, first clarify...) The three-dimensional coordinates, and then Interconnected with the extreme characteristic nodes of the principal deformation axis respectively Secondary deformation axis inflection point feature node Establish a straight-line connection in space, forming two independent spatial lines, namely and and The two connecting lines and and The lines are connected to form a closed triangular topological framework, which clarifies the construction range and boundary of the first inward convergence guiding surface.

[0095] Based on this triangular topological framework, Delaunay triangulation operations are performed, firstly using... Using three basic vertices, initial triangular facets are constructed as the mesh base. Uniform sampling is performed within the triangular topological framework, with a sampling interval controlled at 0.5mm to ensure uniform distribution of sampling points, avoiding overly dense or sparse patterns. The minimum interior angle of each triangular facet is ensured to be no less than 30°, eliminating narrow, sharp, or deformed facets. Each facet is inspected, removing those with a side length ratio exceeding 1.5:1 and correcting those with interior angles greater than 120° or less than 30°, ensuring that the side length deviation of each triangular facet does not exceed 0.3mm. For each triangular facet, three vertices (A, B, and C) are selected. The normal vector of the plane containing the facet is calculated. First, two edge vectors of the facet are constructed, i.e., vectors from A to B and A to C. The original normal vector of the facet is obtained through vector cross product. Then, the original normal vector is divided by its own magnitude for normalization, resulting in a standardized unit normal vector, ensuring that the normal vector magnitude is uniformly 1.

[0096] Select all adjacent triangular faces that share a boundary edge with the current triangular facet. Calculate the unit normal vector of each adjacent facet using the same method described above. Calculate the angle between the current facet and the normal vectors of each adjacent facet using vector dot product. Combine this angle to calculate the curvature of the current facet. That is, take the average of all angles based on the angle between the current facet and the normal vectors of each adjacent facet, and then normalize by dividing by the number of adjacent facets to obtain the curvature value of the current facet. The larger the angle, the larger the calculated curvature value, indicating a more severe degree of surface curvature in that region. Calculate the area of ​​each triangular facet using Heron's formula. First, calculate the lengths of the three sides of the facet (i.e., the magnitudes of each side vector), then calculate the semi-perimeter (half the sum of the three side lengths). Substitute the side lengths and semi-perimeter into Heron's formula to obtain the area of ​​each facet. Sum the areas of all triangular facets to obtain the total area of ​​all facets in the region. Using the area of ​​each facet as a weight, the curvature of a single facet is multiplied by its corresponding area, and the sum is divided by the total area of ​​all faces in the region to obtain a quantified value of the local average curvature of the region. This ensures that the influence of larger facets on the local average curvature is more consistent with the actual surface shape.

[0097] A curvature abrupt change threshold of 0.1 is preset. The calculated local average curvature of each region is compared with this threshold one by one. If the local average curvature is greater than 0.1, it is identified as a curvature abrupt change region, indicating that the surface in this region is severely bent and the existing sampling point density is insufficient, requiring the addition of transitional sampling points. If the local average curvature is not greater than 0.1, it indicates that the surface shape is gentle and no additional sampling points are needed. When adding sampling points, the spacing between adjacent sampling points is controlled to not exceed 0.5 mm, maintaining consistency with the initial sampling spacing. A bicubic interpolation algorithm is used to smooth the junctions of adjacent facets. The curvature difference between adjacent facets is calculated. If the difference is greater than 0.05, transition sampling points are added at the junction to eliminate splicing gaps and ensure that the entire guiding surface is continuous and without abrupt protrusions, with a surface roughness not exceeding 0.02 mm. When the target tissue shrinks (such as tumor regression or passive shrinkage of normal tissue), the guiding surface restricts the irregular deformation of the tissue through its own geometric boundary. Based on the normal direction of each triangular facet, a centripetal constraint force is generated on the surrounding voxels. The amplitude of the force is controlled between 0.8 and 1.2 N, which is proportional to the amplitude of the facet normal vector, causing the tissue to shrink along a preset direction, avoiding boundary distortion and uneven shrinkage. The preset direction is the core direction of the target area pointing to the first initial deformation base point, thereby guiding the tissue to shrink along this directional trajectory. At the same time, the vertex coordinates are adjusted with each iteration to ensure that the coordinate deviation does not exceed 0.1 mm, ensuring that the shrinkage process is traceable and controllable, conforming to the physiological structure.

[0098] Using the second initial deformation base point (Using the expansion-type deformation reference anchor point) as the core, an outward expansion guiding surface is constructed, which is... Interconnected with the extreme nodes of the principal deformation axis Topology transition node Establish spatial straight-line connections to form and , and Two lines, and and The lines are connected to form a closed triangular topological framework, which clarifies the construction scope of the outward guiding surface, ensures that there is no spatial overlap with the inward guiding surface, and that the boundary is naturally connected.

[0099] Using the same segmentation standard as the inward guiding surface, , , Using basic vertices, initial triangular meshes are generated by uniformly sampling at 0.5mm intervals. Each facet is checked, and deformed meshes with interior angles less than 30° or greater than 120°, or side length ratios exceeding 1.5:1, are removed to ensure a standardized mesh topology without redundancy. A Gaussian filtering algorithm is used to smooth and optimize the triangular mesh. The angle between the normal vectors of adjacent faces is calculated by first obtaining the standardized unit normal vectors of two adjacent triangular faces through vector cross product and normalization, then calculating the dot product of the two unit normal vectors, and substituting this dot product into the inverse cosine function to obtain the angle between the normal vectors of two adjacent faces. This ensures the angle does not exceed 15°, eliminating abruptness at face joints and resulting in a smooth guide surface without obvious splicing marks.

[0100] Based on the coordinates of the three vertices of the facet, two edge vectors are constructed, and the original normal vector is obtained through vector cross product. After normalization, the unit normal vector is obtained, and the normal direction is adjusted to be away from the target core (completely opposite to the normal of the inward guiding surface). The vector dot product is verified (the dot product with the vector in the opposite direction to the target core is ≥0.9) to ensure that the normal is consistent and without deviation. When the tissue has expansion needs such as edema or organ filling, the guiding surface exerts an outward constraint and guidance effect on the surrounding voxels through its own geometric boundary, guiding the tissue to deform along the natural expansion direction. At the same time, the coordinates are checked every 5 sampling points to ensure that the coordinate deviation of the sampling points does not exceed 0.1mm, avoiding excessive expansion or local compression, and ensuring that the expansion process conforms to the physiological structural characteristics without abnormal twisting or deformation.

[0101] Both guide surfaces use a 0.5mm sampling interval and a standard inner angle of 30° to 120°. The inward guide surface uses bicubic interpolation smoothing, while the outward guide surface uses Gaussian filtering smoothing to ensure consistent surface flatness (surface roughness ≤ 0.02mm). The normal of the inward guide surface points towards the core of the target area, while the normal of the outward guide surface points away from the core of the target area, with an angle of 180° between their normals, ensuring no directional conflict. The inward guide surface covers the core of the target area and a 5mm radius around it, while the outward guide surface covers the surrounding expansion area of ​​the target area. The coordinate deviation at the boundary connection is ≤ 0.3mm, with no spatial overlap. All coordinate calculations, normal vector solutions, and curvature calculations use the same right-handed coordinate system (x-axis horizontal, y-axis vertical, z-axis depth). The sampling interval, surface calibration standard, and force amplitude (0.8 to 1.2N) remain consistent to ensure the traceability of the entire constraint process.

[0102] This embodiment extracts local curvature feature values ​​through spatial second-order differential gradient analysis, quantifying tissue deformation resistance and potential field curvature. This effectively distinguishes heterogeneous regions such as soft and hard tissues, cavities, and bones, improving the accuracy of positioning key deformation nodes. It locates extreme feature nodes along the main deformation axis, inflection point feature nodes along the secondary deformation axis, and feature nodes in the topological transition zone, comprehensively covering three key locations: the core deformation area, the trend reversal area, and the structural transition area. Calculating the initial node deformation vector based on a continuous, fractional coordinate sequence accurately reflects the temporal deformation pattern of the anatomical structure during the treatment cycle, enhancing the rationality and physical consistency of the deformation trend characterization. Using weighted centroid clustering to determine the initial deformation base point ensures stable base point positions and minimal spatial deviation, avoiding deviations in subsequent guide surface construction due to base point offset. Constructing a deformation guide surface that complements inward convergence and outward expansion simultaneously adapts to two common clinical deformation modes: tissue contraction and tissue expansion, suppressing abrupt changes and non-physical distortions in the deformation field. Anisotropic complementary deformation guidance pairs can provide clear geometric constraints for subsequent deformation field optimization, reduce problems such as contour holes, fractures, and abnormal depressions, and improve the automation level of adaptive delineation.

[0103] In a preferred embodiment of the present invention, step 4 includes:

[0104] Step 400: Extract the geometric center points of the first inward-converging deformation guiding surface and the second outward-expanding deformation guiding surface, along with their corresponding unit normal vectors. Calculate the cosine value of the spatial angle between the two unit normal vectors. Project the first inward-converging deformation guiding surface and the second outward-expanding deformation guiding surface onto the reference imaging plane and calculate the spatial projection overlap area ratio. Perform a weighted fusion operation on the spatial angle cosine value and the spatial projection overlap area ratio using a normalized mapping function to obtain the dynamic boundary coupling coefficient, specifically including:

[0105] Extract all triangular facets from the first inward-converging deformation guiding surface and the second outward-expanding deformation guiding surface, respectively. For each triangular facet, perform the extraction operation of its geometric center point and corresponding unit normal vector. For any triangular facet, select its three vertex coordinates and calculate the geometric center point of the facet using the arithmetic mean method. Calculate the geometric center points of all triangular facets on both guiding surfaces in turn, forming a set of geometric center points for the two guiding surfaces. For each triangular facet, construct two edge vectors based on its three vertex coordinates. Obtain the original normal vector through the vector cross product operation, and then normalize the original normal vector by its own magnitude to obtain the standardized unit normal vector corresponding to the facet. Extract the unit normal vectors of all triangular facets on both guiding surfaces, forming a set of unit normal vectors for the two guiding surfaces. Select any unit normal vector n1 from the set of unit normal vectors of the first inward-converging deformation guiding surface, and select the unit normal vector n2 of the corresponding position patch from the set of unit normal vectors of the second outward-expanding deformation guiding surface. Calculate the dot product of the two vectors. Since both vectors are standardized unit normal vectors with a magnitude of 1, the cosine of the spatial angle is... =n1×n2, calculate the unit normal vector of the corresponding position patch of the two guide surfaces one by one, and take the average value of all calculation results as the final cosine value of the spatial angle.

[0106] To calculate the spatial projection overlap ratio of the two guide surfaces projected onto the reference imaging plane, first determine the reference imaging plane (using a right-handed coordinate system consistent with the previous coordinate calculation, selecting the xy plane as the reference imaging plane). Project all triangular facets of the first inward-converging deformation guide surface and the second outward-expanding deformation guide surface onto this reference imaging plane. The projection process involves retaining the x and y coordinates of the geometric center point of each facet and setting the z coordinate to 0, thus obtaining the two-dimensional projection graphics of the two guide surfaces.

[0107] The pixel scanning method is used to calculate the overlap area and individual projected areas of two two-dimensional projected graphics. A suitable pixel resolution (0.5 mm) is set, and the two projected graphics are discretized into pixel matrices. All pixels are traversed, and it is determined whether a pixel belongs to both projected graphics simultaneously. The number of pixels belonging to the overlap region is counted, and the overlap area is calculated based on the pixel resolution. The number of pixels in each projected graphic is counted separately, and the projected areas of the first inward-converging deformation guiding surface and the second outward-expanding deformation guiding surface are calculated. The smaller of the two values ​​is taken as the reference area. The spatial projection overlap area ratio R is calculated. L =Overlapping area / Reference area. The dynamic boundary coupling coefficient is obtained by weighted fusion calculation of the ratio of the spatial angle cosine value and the spatial projection overlap area using a normalized mapping function. First, the spatial angle cosine value... The ratio of the spatial projection overlap area R is normalized respectively, i.e. =( - ) / ( - ), R L norm=(R L -R L min) / (R L max-R L min), where , These are the minimum and maximum values ​​of the cosine of the spatial angle, R. L min、R L `max` represents the minimum and maximum values ​​of the spatial projection overlap area ratio, respectively, and the normalized values ​​are all in the range [0,1]. Weighting coefficients ω1 and ω2 are set (ω1+ω2=1, where ω1=0.6, ω2=0.4, ensuring that the normal consistency weight is higher than the projection overlap weight). The weighted fusion calculation formula is: K... J =ω1× +ω2×R L norm, where K J This refers to the dynamic boundary coupling coefficient, whose value range is [0,1], K J The closer the value is to 1, the better the normal consistency of the two guiding surfaces, the higher the projection overlap, and the stronger the coupling.

[0108] Step 401: The dynamic boundary coupling coefficient is input as a temporal weight attenuation factor into a preset multi-order inertial low-pass filter to construct a historical state sliding window with the continuous scanning period as the time axis. A weighted integral operation of the historical deformation data across consecutive scan periods is performed on the initial node deformation vector. The damping intensity of the deformation transition between adjacent scan periods is adjusted in real time according to the dynamic boundary coupling coefficient to obtain the temporal deformation state evolution trajectory tensor, specifically including:

[0109] The dynamic boundary coupling coefficient K calculated in step 400 is used as the dynamic boundary coupling coefficient K. J As a time-series weight attenuation factor, it is input into a preset third-order inertial low-pass filter (with the core objectives of smoothing time-series deformation data and suppressing deformation noise, it is constructed using a cascaded recursive structure of three first-order inertial elements, and the transfer function of each first-order inertial element is G(s)=1 / (T)). v s+1), where G(s) is the transfer function and the time constant T vThe value is set to 0.2s, where s is the Laplace operator. The three cascaded components form a third-order inertial low-pass filter. The order is set to third order to ensure smooth filtering to eliminate temporal deformation fluctuations and to consider deformation response speed to match the sub-scanning cycle. Its core parameters are consistent with the previous smoothing parameters. A historical state sliding window is constructed with the continuous sub-scanning cycle as the time axis. The length of the sliding window is set to 5 sub-scanning cycles. The window stores the relevant data corresponding to each sub-scanning cycle. The initial node deformation vector is the previous reference deformation data (associated with the initial deformation base point in step 302b, which is the initial deformation input parameter). The dynamic boundary coupling coefficient is calculated in step 400. The deformation timestamp is the real-time time record of each sub-scanning. The window adopts a sliding update mechanism, that is, after each new sub-scanning is completed, the oldest sub-scanning data in the window is deleted and the latest sub-scanning data is added to ensure that the data in the window is always the latest deformation-related data of the 5 consecutive sub-scanning cycles.

[0110] Perform a weighted integral operation on the initial node deformation vector across historical deformation data, and let the i-th node in the sliding window... S each fraction (i S The initial nodal deformation vectors for nodes (=1,2,3,4,5) are: The corresponding dynamic boundary coupling coefficient is The time weighting coefficient is ( =(T- ) / T d ,in For the i-th S Each time stamp, T d (Given the time span of all intervals within the sliding window, ensuring that newer intervals have higher weights), the weighted integral calculation formula is as follows: = dt, where the integration interval is the time span of the sliding window, calculates the weighted integral deformation vector across multiple divisions. The damping strength of the deformation transition between adjacent divisions is adjusted in real time based on the dynamic boundary coupling coefficient. A base value of damping strength D0 is set (0.3 based on physiological tissue deformation characteristics), and the damping strength adjustment formula is Dp = D0 × K. J When K J The larger the value (the stronger the coupling between the two guiding surfaces), the greater the damping strength Dp, and the smoother the deformation transition between adjacent fractions; when K J The smaller the value (the weaker the coupling between the two guiding surfaces), the smaller the damping strength D, and the more obvious the deformation transition between adjacent fractions. This avoids excessive suppression when deformation uncertainty is high, preserving a reasonable deformation response. The weighted integral deformation vector, damping strength, and deformation state parameters of each fraction are integrated to construct a time-series deformation state evolution trajectory tensor, which has a dimension of 3×N. b ×Mb (3 corresponds to three-dimensional spatial coordinates, N) b N is the number of nodes. b (The number of times within the sliding window), each element in the tensor corresponds to the three-dimensional deformation state parameters of a certain time segment and a certain node, and completely records the deformation evolution process of the node in the continuous scanning cycle.

[0111] In this embodiment, by extracting the geometric center points and unit normal vectors of the two guiding surfaces, calculating the ratio of the spatial angle cosine to the spatial projection overlap area, and weighting and fusing them, the resulting dynamic boundary coupling coefficient can accurately quantify the coupling degree of the two guiding surfaces, effectively avoiding conflict in the constraint directions of the two guiding surfaces, and ensuring the coordination and accuracy of deformation constraints. Using the dynamic boundary coupling coefficient as a temporal weight attenuation factor, combined with a multi-order inertial low-pass filter and a historical state sliding window, a cross-segment weighted integral operation is performed on the initial node deformation vector. Simultaneously, the deformation transition damping intensity is dynamically adjusted. The constructed temporal deformation state evolution trajectory tensor can completely record the temporal evolution process of node deformation, effectively suppressing temporal deformation noise, avoiding excessive deformation transitions between adjacent segments, and ensuring the continuity and stability of the deformation process.

[0112] In a preferred embodiment of the present invention, step 5 includes:

[0113] Step 500: Obtain the temporal deformation state evolution trajectory tensor, construct a joint optimization objective function including data fidelity constraints and topological connectivity penalties, and perform gradient descent operation with the temporal deformation state evolution trajectory tensor as the initial state for iteration. During the iteration process, dynamically adjust the weight ratio of each constraint term in the joint optimization objective function according to the current deformation residual until the residual change between adjacent iteration steps meets the preset convergence threshold, thus obtaining the multi-constraint dynamic equilibrium steady-state deformation field, specifically including:

[0114] First, obtain the temporal deformation state evolution trajectory tensor constructed in step 401. Using this tensor as the initial state for iteration, construct a joint optimization objective function that includes data fidelity constraints and topological connectivity penalties. ,in For the deformation field parameter vector, Assign weights to data fidelity constraint terms. The weights are the topological connectivity penalty terms. The initial weights are set to Ensuring data fidelity is a top priority. Data fidelity constraints. Used to ensure consistency between the deformation field and the original deformation data. ,in This represents the total number of nodes. The initial deformation parameter in the temporal deformation state evolution trajectory tensor is used as a constraint term to ensure that the deformation field does not deviate from the original temporal deformation law by calculating the mean square error between the current deformation field parameter and the initial deformation parameter.

[0115] Topological connectivity penalty term Used to ensure that the topology of nodes does not break in a deformation field. ,in This represents the total number of edges connecting the nodes. , The deformation parameters of adjacent nodes are used as the constraint term. This term penalizes cases where the deformation differences between adjacent nodes are too large by calculating the mean square error of the deformation parameters of adjacent nodes, thus maintaining the connectivity of the topology. Using the temporal deformation state evolution trajectory tensor as the initial state for iteration, gradient descent is performed. The learning rate for gradient descent is set to 0.01 (balancing iteration speed and convergence stability). The iteration process is as follows: ,in For the first Deformation field parameters of the next iteration. For the first Deformation field parameters of +1 iteration For learning rate, For the objective function in The gradient at that point.

[0116] During the iteration process, the weight ratio of each constraint term in the joint optimization objective function is dynamically adjusted according to the current deformation residual. The calculation process of the current deformation residual is as follows: ,in For the first Deformation residual of the next iteration , The first sequence The objective function value of the next iteration. At the same time, increase the weight of the data fidelity constraint. Reduce the weight of the topological connectivity penalty term to 0.8. Up to 0.2, prioritize data fidelity; when At that time, , Restored to initial values ​​of 0.7 and 0.3, balancing data fidelity and topological connectivity. The preset convergence threshold is... During the iteration process, the residual change between adjacent iteration steps is continuously calculated. ,when When the iteration stops, the deformation field parameter vector X obtained at this time is the multi-constraint dynamic equilibrium steady-state deformation field. This deformation field not only conforms to the temporal deformation evolution law, but also maintains the connectivity of the topological structure.

[0117] Step 501: Perform 3D surface meshing discretization on the reference target area boundary features in the initial spatiotemporal registration feature manifold to obtain the initial target area mesh; perform spatial coordinate resampling operation on the multi-constraint dynamic equilibrium steady-state deformation field forward mapped to the initial target area mesh to obtain the reconstructed boundary mesh; perform topological connectivity detection and geometric self-intersection correction operation on the reconstructed boundary mesh to remove mesh breaks and abnormal normal flipping regions, and obtain the adaptive topological reconstruction target area contour set, specifically including:

[0118] The boundary features of the reference target region in the initial spatiotemporal registration feature manifold are obtained. These boundary features are three-dimensional continuous surfaces. Three-dimensional surface meshing is performed on them, using the same standard as the previous Delaunay triangulation. The sampling interval is set to 0.5 mm. The boundary features of the reference target region are uniformly sampled to generate the initial triangular mesh, i.e. the initial target region mesh. The interior angles of the facets of the initial target region mesh are controlled between 30° and 120°, and the side length ratio is ≤1.5:1 to ensure that the mesh shape is regular.

[0119] The multi-constraint dynamic equilibrium steady-state deformation field obtained in step 500 is positively mapped to the initial target area mesh, and spatial coordinate resampling is performed. Let the coordinates of a certain vertex in the initial target area mesh be P(x,y,z), and the deformation mapping function of the steady-state deformation field be T(X,Y,Z)=(X b ,Y b Z b ), (X b ,Y b Z b P(x,y,z) represents the deformed spatial coordinates. Substituting P(x,y,z) into the mapping function yields the resampled vertex coordinates. This mapping operation is then performed on all vertices of the initial target area mesh, and the mapped vertices are reconnected to form the reconstructed boundary mesh. Topological connectivity detection and geometric self-intersection correction are performed on the reconstructed boundary mesh: Topological connectivity detection uses a half-edge data structure, traversing all edges of the mesh, marking the number of associated faces for each edge, and determining whether there are isolated edges, isolated vertices, and mesh break regions. If so, transition vertices and edges are added at a sampling interval of 0.5mm to connect break regions and restore topological connectivity. Geometric self-intersection correction uses the ray method, traversing all triangular faces, and detecting self-intersection by judging the intersection of rays and faces. If self-intersection exists, the coordinates of the vertices of the intersecting faces are adjusted by a fine adjustment of 0.05mm (fine-tuning in the opposite direction of the face normal to ensure that the faces do not overlap after fine-tuning), while eliminating self-intersecting regions. Simultaneously, the normal vector dot product is used to verify and eliminate abnormal regions of normal vector flipping (if the angle between the normal vectors of a face and its adjacent face is >90°, it is judged as abnormal), and finally, the adaptive topological reconstruction target region contour set is obtained.

[0120] Step 502: Perform boundary loop hole detection on the adaptive topology reconstruction target region contour set. Traverse all boundary edges in the adaptive topology reconstruction target region contour set, identify and extract closed boundary loops that are connected end to end through a half-edge data structure, and mark the internal region enclosed by each closed boundary loop as a hole region with missing triangular facets. Specifically, this includes:

[0121] Boundary loop and hole detection is performed on the adaptive topology reconstruction target area contour set obtained in step 501. First, all boundary edges in the contour set are traversed. A half-edge data structure is used to uniquely mark each boundary edge, explicitly recording the coordinates of the starting vertex and the ending vertex of each edge, as well as the number of triangular faces to which the edge belongs (normal non-boundary edges belong to 2 triangular faces, and boundary edges belong to 1 triangular face). This is used to initially filter out all boundary edges and eliminate non-boundary edges to reduce detection redundancy. Based on the correlation characteristics of the half-edge data structure, closed loop extraction is performed on each of the filtered boundary edges. Taking any boundary edge as the starting edge, its starting vertex is recorded as A and its ending vertex as B. Then, the next boundary edge with B as the starting vertex is found (ensuring that the two edges are spatially continuous). This process is iterated until a boundary edge with A as the ending vertex is found, forming a closed boundary loop that is connected end to end without breaks. After extracting each closed loop, all boundary edges within the loop are marked to avoid repeated extraction.

[0122] For each extracted closed boundary loop, a hole determination is performed. The determination criteria are clearly defined in two aspects: first, the internal region enclosed by the closed loop must not be covered by any triangular facets, and the spatial extent of this region must be continuous with the surrounding target area contour mesh without significant disconnection; second, the normal direction of the closed loop must be consistent with the normal direction of the surrounding adjacent triangular facets, verified by the dot product of normal vectors (the dot product of the loop normal vector and the surrounding facet normal vectors ≥ 0.8, and the normal deviation not exceeding 37°), ensuring compliance with the overall geometric rules of the target area contour and avoiding misjudging normal contour boundaries as holes. For closed boundary loops that meet the above two determination criteria, the internal region enclosed by them is marked as a hole region lacking triangular facets. Simultaneously, key information for each hole region is recorded in detail, including the vertex sequence on the boundary loop (recorded in clockwise order), the total length of the loop, the coordinates of the hole center, and the hole area.

[0123] Step 503: For each hole region, obtain the vertex sequence on the boundary loop, calculate the edge length and dihedral angle between adjacent vertices on the boundary loop, select the current advancing edge according to the minimum interior angle criterion, estimate the coordinates of the new vertex using Laplace weighted average, and perform smoothing optimization on the new vertex. Connect the advancing edge and the new vertex sequentially to form a filled triangular patch. Iterate the optimization process until the hole region is completely covered by the new triangular patch, obtaining the new filled triangular patch, specifically including:

[0124] For each hole region marked in step 502, first obtain the vertex sequence on the boundary loop of the hole region, let the vertex sequence be U1, U2, ..., ( (The number of vertices in the boundary loop) Calculate the edge length and dihedral angle between adjacent vertices on the boundary loop (obtain the normal vectors corresponding to the two edge vectors through vector cross product, and then calculate the angle between the two normal vectors through vector dot product; this angle is the dihedral angle). Calculate the average value of the two adjacent dihedral angles corresponding to each edge on the boundary loop, and select the edge with the smallest average value as the current advancing edge to ensure that the filled triangular patches have a regular shape and avoid deformed patches.

[0125] The coordinates of the new vertex are estimated by Laplace weighted average, i.e., by selecting the two endpoints of the current advancing front edge. ( ), and the boundary ring vertices adjacent to these two endpoints. Set weighting coefficients ( , The weight is 0.3. (with a weight of 0.2), new vertex ( The coordinates of ) are calculated as follows: ; ; ;

[0126] The new vertex is smoothed using a Gaussian filtering algorithm. Centered on the new vertex, three adjacent vertices are selected. The average x, y, and z coordinates of these three adjacent vertices are calculated. Then, the new vertex's x, y, and z coordinates are weighted and fused with this average value. Specifically, the initial x, y, and z coordinates of the new vertex are multiplied by a weight of 0.7, and the corresponding average coordinates are multiplied by a weight of 0.3. The two sets of weighted values ​​are then summed to obtain the smoothed x, y, and z coordinates of the new vertex. This correction ensures the geometric continuity between the new vertex and its surrounding boundary vertices. A filled triangular facet is formed by connecting the current advancing edge to the new vertex and adding it to the mesh of the hole region. This process is repeated, continuously selecting new advancing edges, estimating new vertices, and constructing filled triangular facets, iteratively optimizing until the hole region is completely covered by the new triangular facets, with no blank areas. At this point, the new filled triangular facets corresponding to all hole regions are obtained.

[0127] Step 504 involves performing normal consistency calibration on the newly filled triangular facets and performing neighborhood matching with the curvature distribution of the surrounding original triangular facets. This ensures a continuous transition between the filled region and the original geometry of the adaptive topology reconstruction target region contour set, resulting in a complete target region contour set without holes. Specifically, this includes:

[0128] First, perform normal consistency calibration on the new filled triangular facets obtained in step 503. The calibration method is the same as the normal calibration step corresponding to the first initial deformation base point, that is, by combining vector cross product with normalization, the unit normal vector of each new filled triangular facet is obtained, and its normal direction is adjusted to ensure that it is consistent with the normal direction of the surrounding original triangular facets, with a normal deviation of no more than 5°, and to maintain consistency with the previous normal calibration standard. Perform neighborhood matching on the curvature distribution of the new filled triangular facets and the surrounding original triangular facets, calculate the local average curvature of each new filled triangular facet (the calculation method is the same as the local average curvature calculation method in step 302b), and at the same time calculate the local average curvature of the three surrounding original triangular facets. Calculate the curvature difference between the new filled facet and each surrounding original facet. If the curvature difference is >0.05, adjust the vertex coordinates of the new filled facet to correct its curvature until the curvature difference is ≤0.05. Through the above-mentioned normal consistency calibration and curvature neighborhood matching, the original geometric shape of the filled region and the adaptive topology reconstruction target region contour set is continuously transitioned without obvious splicing marks, and the surface roughness is ≤0.02mm; finally, all the original triangular patches and the optimized new filled triangular patches are integrated to obtain a complete target region contour set with no holes, topological connectivity and regular geometric shape.

[0129] In this embodiment, the multi-constraint dynamic equilibrium steady-state deformation field balances the fidelity of deformation data with the connectivity of the topological structure, effectively avoiding anomalies such as distortion and breakage in the deformation field, and ensuring that the deformation field conforms to the deformation law of physiological tissue. The adaptive topology reconstruction target area contour set can accurately match the steady-state deformation field, eliminate abnormal mesh regions, and ensure the topological integrity and geometric accuracy of the target area contour. A half-edge data structure is used to identify closed boundary loops and mark hole regions, effectively solving the problems of hole omission and misjudgment during target area contour reconstruction, ensuring the comprehensiveness and accuracy of hole detection, and providing accurate boundary references for hole filling. The advancing front edge is selected according to the minimum interior angle criterion, and new vertices are estimated and smoothed using a Laplacian weighted average. The resulting triangular facets are regular in shape and can accurately fit the hole boundaries, effectively avoiding problems such as deformity and breakage of the filled facets, and ensuring the effectiveness and regularity of hole filling. By using normal consistency calibration and curvature neighborhood matching, the filling area and the original target area contour geometry are continuously transitioned, eliminating splicing traces. The final result is a complete target area contour set without holes, which is topologically connected and geometrically regular, and can accurately restore the actual shape of the target area.

[0130] like Figure 2As shown, embodiments of the present invention also provide a radiotherapy target volume adaptive delineation system based on CBCT images, comprising:

[0131] The extraction module is used to acquire reference positioning CT images and current fractional CBCT images, perform multi-scale spatial frequency domain alignment and feature extraction processing, and obtain the initial spatiotemporal registration feature manifold;

[0132] The mapping module is used to perform high-dimensional topological reconstruction and feature affinity calculation on the initial spatiotemporal registration feature manifold, construct a spatial interaction intensity mapping model and transform the spatial interaction intensity mapping model into a local structural potential energy gradient distribution, and obtain the local heterogeneous potential field gradient distribution matrix.

[0133] The module is used to perform curvature tensor analytical operations and spatiotemporal state projection on the gradient distribution matrix of the local heterogeneous potential field, extract local curvature eigenvalues, and solve for the extreme feature nodes of the principal deformation axis, the inflection point feature nodes of the secondary deformation axis, and the feature nodes of the topological transition region. Based on the difference of the projected coordinates of the three feature nodes, the initial node deformation vector is derived, and the first initial deformation base point and the second initial deformation base point are extracted. Based on the first initial deformation base point and the second initial deformation base point, the first inward convergent deformation guiding surface and the second outward expansion deformation guiding surface are constructed.

[0134] The calculation module is used to calculate the dynamic boundary coupling coefficient based on the complementary normal overlap between the first inward convergent deformation guiding surface and the second outward expansion deformation guiding surface, and to use the dynamic boundary coupling coefficient to perform temporal smoothing filtering on the initial node deformation vector to obtain the temporal deformation state evolution trajectory tensor.

[0135] The reconstruction module is used to perform multi-objective non-convex optimization and dynamic potential energy balance iteration on the temporal deformation state evolution trajectory tensor to obtain a multi-constraint dynamic equilibrium steady-state deformation field; and to perform spatial inverse mapping reconstruction and verification on the multi-constraint dynamic equilibrium steady-state deformation field to obtain a complete target region contour set without holes.

[0136] It should be noted that this system is a system corresponding to the above method. All implementation methods in the above method embodiments are applicable to this embodiment and can achieve the same technical effect.

[0137] Embodiments of the present invention also provide a computing device, including: a processor and a memory storing a computer program, wherein the computer program, when executed by the processor, performs the method described above. All implementations in the above method embodiments are applicable to this embodiment and can achieve the same technical effects.

[0138] Embodiments of the present invention also provide a computer-readable storage medium storing instructions that, when executed on a computer, cause the computer to perform the method described above. All implementations in the above method embodiments are applicable to this embodiment and can achieve the same technical effects.

[0139] The above description represents the preferred embodiments of the present invention. It should be noted that those skilled in the art can make various improvements and modifications without departing from the principles of the present invention, and these improvements and modifications should also be considered within the scope of protection of the present invention.

Claims

1. A method for adaptive delineation of radiotherapy target areas based on CBCT images, characterized in that, The method includes: Step 1: Acquire the reference positioning CT image and the current fractional CBCT image, perform multi-scale spatial frequency domain alignment and feature extraction processing to obtain the initial spatiotemporal registration feature manifold; Step 2: Perform high-dimensional topological reconstruction and feature affinity calculation on the initial spatiotemporal registration feature manifold, construct a spatial interaction intensity mapping model and transform the spatial interaction intensity mapping model into a local structural potential energy gradient distribution to obtain the local heterogeneous potential field gradient distribution matrix; Step 3: Perform curvature tensor analytical operations and spatiotemporal state projection on the local heterogeneous potential field gradient distribution matrix, extract local curvature eigenvalues, and solve for the extreme value feature nodes of the principal deformation axis, the inflection point feature nodes of the secondary deformation axis, and the topological transition region feature nodes; derive the initial node deformation vectors based on the projection coordinate difference of the three feature nodes, and extract the first initial deformation base point and the second initial deformation base point; construct the first inward convergent deformation guiding surface and the second outward expansion deformation guiding surface based on the first initial deformation base point and the second initial deformation base point. Step 4: Calculate the dynamic boundary coupling coefficient based on the complementary normal overlap between the first inward convergent deformation guiding surface and the second outward expansion deformation guiding surface. Use the dynamic boundary coupling coefficient to perform temporal smoothing filtering on the initial node deformation vector to obtain the temporal deformation state evolution trajectory tensor. Step 5: Perform multi-objective non-convex optimization and dynamic potential energy balance iteration on the temporal deformation state evolution trajectory tensor to obtain a multi-constraint dynamic equilibrium steady-state deformation field; perform spatial inverse mapping reconstruction and verification on the multi-constraint dynamic equilibrium steady-state deformation field to obtain a complete target region contour set without holes.

2. The method for adaptive delineation of radiotherapy target areas based on CBCT images according to claim 1, characterized in that, Step 1 includes: Acquire reference positioning CT images and current fractionated CBCT images, perform three-dimensional discrete frequency domain orthogonal transformation on the reference positioning CT images and current fractionated CBCT images to separate the low-frequency energy baseband and the high-frequency detail band; Cross-power spectrum phase normalization is performed on the high-frequency detail bands to extract the cross-modal frequency domain phase difference distribution; inverse frequency domain spatial mapping transformation is performed on the cross-modal frequency domain phase difference distribution to obtain the spatial cross-correlation response kernel; based on the spatial cross-correlation response kernel, sliding window local mutual information entropy calculation and local gradient magnitude matching weighted aggregation are performed to obtain the cross-modal spatial correlation distribution map; Based on the cross-modal spatial correlation distribution map, global peak extremum search and sub-pixel level offset interpolation are performed to solve the initial rigid spatial transformation parameters and construct the initial rigid transformation alignment field. Based on the initial rigid transformation alignment field, voxel-level trilinear interpolation resampling and alignment are performed on the current CBCT images to obtain the aligned image data. Nonlocal mean denoising and edge-aware convolution operations are performed on the aligned image data to extract a multi-scale texture feature vector set. The multi-scale texture feature vector set is mapped to a high-dimensional tensor space, and Riemannian manifold dimensionality reduction embedding is performed to project the high-dimensional features onto a low-dimensional continuous metric space to obtain the initial spatiotemporal registration feature manifold.

3. The method for adaptive delineation of radiotherapy target areas based on CBCT images according to claim 2, characterized in that, Step 2 includes: Extract the voxel neighborhood spatial topological connectivity in the initial spatiotemporal registration feature manifold and construct the K-nearest neighbor feature association graph; calculate the texture gradient cosine similarity and spatial Euclidean distance decay weights between adjacent nodes in the K-nearest neighbor feature association graph, and obtain the feature affinity weight matrix through Gaussian kernel function nonlinear mapping; Using the feature affinity weight matrix as the driving parameter, a coupling strength field fitting operation based on local kernel density estimation is performed, mapping the feature affinity weights to the multidimensional force transfer potential function between nodes, iteratively updating the parameters of the multidimensional force transfer potential function until the local energy distribution converges, and constructing a spatial interaction strength mapping model. By performing point-by-point potential field integration on the voxel spatial coordinates of the initial spatiotemporally registered characteristic manifold using the spatial interaction intensity mapping model, a scalar field of local structural potential energy distribution is obtained. Then, the scalar field of local structural potential energy distribution is convolved with the discrete difference gradient operator to transform it into a spatial directional gradient vector field, thus obtaining the local structural potential energy gradient distribution. Finally, the local structural potential energy gradient distribution is subjected to voxel gridded discrete rearrangement and matrix dimension alignment to obtain the local heterogeneous potential field gradient distribution matrix.

4. The method for adaptive delineation of radiotherapy target areas based on CBCT images according to claim 3, characterized in that, Perform curvature tensor analytical operations and spatiotemporal state projection on the gradient distribution matrix of the local heterogeneous potential field, extract local curvature eigenvalues, and solve for the extreme feature nodes of the principal deformation axis, the inflection point feature nodes of the secondary deformation axis, and the feature nodes of the topological transition region, including: A spatial second-order differential gradient analytical operation is performed on the gradient distribution matrix of the local heterogeneous potential field to extract local curvature feature values ​​that characterize the degree of spatial curvature of the local potential field and the deformation resistance gradient. Based on the local curvature feature values, a spatial coordinate mapping grid is constructed. An extreme value search is performed in the spatial coordinate mapping grid along the direction of potential energy gradient ascent. The spatial coordinates are located at the center of the voxel region where the local potential field gradient amplitude is greater than the preset gradient aggregation threshold and the spatial connectivity satisfies the eight-neighborhood continuity condition. This location is recorded as the extreme feature node of the main deformation axis. By mapping the potential energy gradient of the grid along the spatial coordinates and performing contour tracing along the boundary where the sign of the grid potential energy gradient is reversed, the location of the grid boundary where the spatial coordinates are at the interface of abrupt change in the resistance to deformation and where the sign of the local curvature characteristic value is continuously reversed is located, and this location is recorded as the inflection point feature node of the secondary deformation axis. Extract the centroid location of the connected region whose spatial coordinates are distributed within the projection range of the line connecting the extreme feature node of the primary deformation axis and the inflection point feature node of the secondary deformation axis, and whose local curvature feature value fluctuation amplitude is lower than the preset smooth threshold. This location is denoted as the topological transition region feature node.

5. The method for adaptive delineation of radiotherapy target areas based on CBCT images according to claim 4, characterized in that, The initial node deformation vector is derived by the difference of the projected coordinates of the three feature nodes, and the first initial deformation base point and the second initial deformation base point are extracted. Based on the first initial deformation base point and the second initial deformation base point, a first inward convergent deformation guiding surface and a second outward expansion deformation guiding surface are constructed, including: Extract the spatial projection coordinate sequence of the extreme value feature nodes of the main deformation axis, the inflection point feature nodes of the secondary deformation axis, and the topological transition zone feature nodes in consecutive scanning cycles. Calculate the temporal position difference of the spatial projection coordinate sequence. Combine the time interval parameter to perform deformation integral deduction on the temporal position difference to obtain the initial node deformation vector. Based on the spatial distribution matrix of local curvature eigenvalues, a weighted centroid clustering operation is performed to extract the coordinates of the two spatial anchor points with the smallest spatial deviation, which are denoted as the first initial deformation base point and the second initial deformation base point, respectively. By performing spatial vertex connection and patch meshing operations on the first initial deformation base point, the extreme feature node of the main deformation axis, and the inflection point feature node of the secondary deformation axis, a first inwardly converging deformation guiding surface is obtained. The first inwardly converging deformation guiding surface provides an inward contraction geometric boundary traction constraint for the rapid reduction of the target tissue volume by pointing the corresponding normal to the interior of the deformation region. By performing spatial vertex connection and patch meshing operations on the second initial deformation base point, the topological transition region feature node, and the extreme feature node of the main deformation axis, a second outwardly expanding deformation guiding surface is obtained. The second outwardly expanding deformation guiding surface provides an outwardly extending geometric boundary pushing constraint for the edema expansion state of the tissue around the target region by pointing the corresponding normal to the exterior of the deformation region. The first inwardly converging deformation guiding surface and the second outwardly expanding deformation guiding surface constitute an anisotropic deformation guiding pair with opposite spatial normals and complementary functions.

6. The method for adaptive delineation of radiotherapy target areas based on CBCT images according to claim 5, characterized in that, Step 4 includes: The geometric center points of the first inward convergent deformation guide surface and the second outward expansion deformation guide surface and their corresponding unit normal vectors are extracted respectively. The cosine value of the spatial angle between the two unit normal vectors is calculated. The first inward convergent deformation guide surface and the second outward expansion deformation guide surface are projected onto the reference imaging plane to calculate the spatial projection overlap area ratio. The cosine value of the spatial angle and the spatial projection overlap area ratio are weighted and fused by the normalized mapping function to obtain the dynamic boundary coupling coefficient. The dynamic boundary coupling coefficient is used as a time-series weight attenuation factor input to a preset multi-order inertial low-pass filter. A historical state sliding window with the continuous scanning period as the time axis is constructed. The weighted integral operation of the historical deformation data across the initial node deformation vector is performed. The damping intensity of the deformation transition between adjacent periods is adjusted in real time according to the dynamic boundary coupling coefficient to obtain the time-series deformation state evolution trajectory tensor.

7. The method for adaptive delineation of radiotherapy target areas based on CBCT images according to claim 6, characterized in that, Step 5 includes: Obtain the temporal deformation state evolution trajectory tensor, construct a joint optimization objective function containing data fidelity constraint terms and topological connectivity penalty terms, perform gradient descent operation with the temporal deformation state evolution trajectory tensor as the initial state of iteration, and dynamically adjust the weight ratio of each constraint term in the joint optimization objective function according to the current deformation residual during the iteration process until the residual change between adjacent iteration steps meets the preset convergence threshold, and obtain the multi-constraint dynamic equilibrium steady-state deformation field; The reference target region boundary features in the initial spatiotemporal registration feature manifold are discretized using three-dimensional surface meshing to obtain the initial target region mesh. The multi-constraint dynamic equilibrium steady-state deformation field is positively mapped to the initial target region mesh, and spatial coordinate resampling is performed to obtain the reconstructed boundary mesh. Topological connectivity detection and geometric self-intersection correction are performed on the reconstructed boundary mesh to remove mesh breakage and normal flip anomaly regions, resulting in an adaptive topological reconstruction target region contour set. Boundary loop hole detection is performed on the adaptive topology reconstruction target area contour set. All boundary edges in the adaptive topology reconstruction target area contour set are traversed. Closed boundary loops that are connected end to end are identified and extracted through half-edge data structure. The internal area enclosed by each closed boundary loop is marked as the hole area of ​​missing triangular facet. For each hole region, obtain the vertex sequence on the boundary loop, calculate the edge length and dihedral angle between adjacent vertices on the boundary loop, select the current advancing edge according to the minimum interior angle criterion, estimate the coordinates of the new vertex by Laplace weighted average and perform smooth optimization on the new vertex, connect the advancing edge and the new vertex in sequence to form a filling triangle, iterate the optimization process until the hole region is completely covered by the new triangle, and obtain the new filling triangle. Perform normal consistency calibration on the newly filled triangular facets and perform neighborhood matching with the curvature distribution of the surrounding original triangular facets to ensure a continuous transition between the filled region and the original geometry of the adaptive topology reconstruction target region contour set, resulting in a complete target region contour set without holes.

8. A radiotherapy target volume adaptive delineation system based on CBCT images, wherein the system implements the method as described in any one of claims 1 to 7, characterized in that, include: The extraction module is used to acquire reference positioning CT images and current fractional CBCT images, perform multi-scale spatial frequency domain alignment and feature extraction processing, and obtain the initial spatiotemporal registration feature manifold; The mapping module is used to perform high-dimensional topological reconstruction and feature affinity calculation on the initial spatiotemporal registration feature manifold, construct a spatial interaction intensity mapping model and transform the spatial interaction intensity mapping model into a local structural potential energy gradient distribution, and obtain the local heterogeneous potential field gradient distribution matrix. The module is used to perform curvature tensor analytical operations and spatiotemporal state projection on the gradient distribution matrix of the local heterogeneous potential field, extract local curvature eigenvalues ​​and solve for the extreme value feature nodes of the principal deformation axis, the inflection point feature nodes of the secondary deformation axis and the feature nodes of the topological transition region; the initial node deformation vector is derived from the difference of the projected coordinates of the three feature nodes and the first initial deformation base point and the second initial deformation base point are extracted. Construct a first inward convergent deformation guide surface and a second outward expansion deformation guide surface based on the first initial deformation base point and the second initial deformation base point; The calculation module is used to calculate the dynamic boundary coupling coefficient based on the complementary normal overlap between the first inward convergent deformation guiding surface and the second outward expansion deformation guiding surface, and to use the dynamic boundary coupling coefficient to perform temporal smoothing filtering on the initial node deformation vector to obtain the temporal deformation state evolution trajectory tensor. The reconstruction module is used to perform multi-objective non-convex optimization and dynamic potential energy balance iteration on the temporal deformation state evolution trajectory tensor to obtain a multi-constraint dynamic equilibrium steady-state deformation field; and to perform spatial inverse mapping reconstruction and verification on the multi-constraint dynamic equilibrium steady-state deformation field to obtain a complete target region contour set without holes.

9. A computing device, characterized in that, Including memory and processor; The memory stores one or more computer programs, the one or more computer programs including instructions; when the instructions are executed by the processor, the computing device performs the method as described in any one of claims 1 to 7.

10. A computer-readable storage medium, characterized in that, The computer-readable storage medium is used to store a computer program for performing the method as described in any one of claims 1 to 7.

Images