Vascular interventional operation multi-view image spatial expansion and navigation method

Abstract

The invention is applicable to the technical field of vascular interventional operations, and provides a vascular interventional operation multi-view image spatial expansion and navigation method, which comprises the following steps of: extracting a blood vessel center line; performing spatial expansion on the DSA image; carrying out elastic energy minimization constraint optimization on a spatial expansion result; and performing man-machine interaction optimization based on AR / VR. According to the blood vessel center line extraction method based on the MRA image, the probability optimization ray tracing method based on ray beam geometry and the elastic energy minimization constraint optimization algorithm, real-time three-dimensional guide wire reconstruction and accurate blood vessel deformation prediction are achieved through the multi-view image fusion technology, and the real-time three-dimensional guide wire reconstruction and the accurate blood vessel deformation prediction are achieved. The problem that the DSA image is difficult to quickly and accurately map to the three-dimensional space in the vascular interventional operation is solved. By means of the method, the safety and effect of the interventional operation are improved, high-timeliness and high-precision navigation is achieved, and a revolutionary breakthrough is brought to the field of the vascular interventional operation.

Description

Technical Field

[0001] The present invention belongs to the technical field of vascular interventional surgery, and particularly relates to a method for expanding and navigating the multi-view imaging space of vascular interventional surgery. Background Art

[0002] In vascular interventional surgery centered around guidewires and catheters, a vascular access is established through percutaneous puncture and the diagnosis and treatment are completed under image guidance, significantly reducing the surgical trauma and providing a safer and more effective treatment option for patients with cardiovascular and cerebrovascular peripheral vascular diseases. With the miniaturization and intelligence of interventional devices and the continuous expansion of the indications for endovascular treatment, the role of the intraoperative precise navigation system has become increasingly crucial. In modern interventional surgery, doctors mainly rely on a two-dimensional digital subtraction angiography (DSA) system for navigation. It relies on the real-time fluoroscopic images provided by a C-arm X-ray machine and combines selective angiography techniques to display the target blood vessels. However, this navigation method highly depends on the doctor's clinical experience and spatial imagination ability, and has limitations in terms of real-time performance and accuracy. For example, it requires doctors to develop a three-dimensional perception of two-dimensional projection images through long-term practice, and at the same time, it is difficult to accurately locate the spatial position of complex vascular lesions and manipulate surgical instruments. To overcome these limitations, modern interventional navigation technologies are developing towards multi-modal image fusion, real-time three-dimensional reconstruction, and intelligent assisted decision-making. The fusion navigation system integrates various imaging data such as CT, MRI, and DSA, and uses artificial intelligence algorithms to provide more accurate anatomical structure information. This system not only supports intraoperative real-time path planning and instrument positioning, but also significantly improves the safety and accuracy of the surgery.

[0003] In specific implementation, the existing DSA three-dimensional reconstruction mainly adopts rotational angiography (3D-RA) technology. This process includes: after injecting the contrast agent, the C-arm rotates at a fixed angular velocity to collect a sequence of continuous projection images within a range of about 180 - 200 degrees (usually 50 - 120 frames); subsequently, these two-dimensional projections are reconstructed into three-dimensional volume data through a cone-beam CT reconstruction algorithm (such as the Feldkamp algorithm). In the data acquisition stage, a mask scan is performed to obtain the background image, then the contrast agent is injected for contrast scanning, and a pure vascular projection sequence is obtained after digital subtraction; in the image reconstruction stage, the system matrix is calculated and projection mapping is performed to obtain a three-dimensional image. However, this process requires the patient to hold their breath and requires precise control of the contrast agent injection timing to ensure sufficient vascular visualization throughout the acquisition process. This not only increases the burden on the patient but may also cause motion artifacts due to physiological activities (such as respiratory movement and cardiac pulsation), affecting the reconstruction accuracy. Secondly, due to the need to process a large amount of projection data, the computational complexity is high, and there are a large number of matrix operations and iterative optimization processes, resulting in the reconstruction time often taking dozens of seconds or even longer. This time delay in acquisition and reconstruction makes the system unable to respond in real time to the dynamic changes in vascular morphology, severely restricting its application in interventional surgery navigation. In addition, the distribution and dilution process of the contrast agent in the blood vessels will cause contrast differences between different projection frames, affecting the accurate positioning of the vascular boundary. At the level of the reconstruction algorithm, traditional cone-beam CT reconstruction is limited by discrete sampling and interpolation errors, with limited spatial resolution and prone to generating artifacts in the vascular overlap area. Especially for small blood vessels and complex bifurcation structures, the reconstruction quality often cannot meet the requirements of precise navigation. At the same time, system geometric calibration errors and mechanical motion errors will also reduce the reconstruction accuracy, and these accumulated errors directly affect the accuracy of navigation.

[0004] In view of the above problems, the present invention proposes a multi-view image space expansion and navigation method for vascular interventional surgery. Summary of the Invention

[0005] The object of the present invention is to provide a multi-view image space expansion and navigation method for vascular interventional surgery, aiming to solve the problems raised in the above background technology.

[0006] The object of the present invention is achieved through the following technical solutions:

[0007] A multi-view image space expansion and navigation method for vascular interventional surgery includes the following steps:

[0008] Step 1: Extract the vascular centerline, including:

[0009] Step 11: Preprocess the MRA image using an improved bilateral filtering method;

[0010] Step 12: Extraction of vascular anisotropic features;

[0011] Step 13: Construction of Riemann metric;

[0012] Step 14: Calculation of the shortest path;

[0013] Step 15: Optimization of the centerline;

[0014] Step 2: Spatial expansion of DSA images, including:

[0015] Step 21: System geometric calibration;

[0016] Step 22: Image feature extraction;

[0017] Step 23: Construction of the ray bundle probability model;

[0018] Step 24: Vascular surface reconstruction;

[0019] Step 25: Spatial expansion of the guide wire based on the ray bundle probability model;

[0020] Step 3: Optimization of the spatially expanded result with minimum elastic energy constraint, including:

[0021] Step 31: Initialization and discretization;

[0022] Step 32: Construction of the energy function;

[0023] Step 33: Solution of the energy function;

[0024] Step 4: Human-computer interaction optimization based on AR / VR.

[0025] Furthermore, the specific steps of Step 1 include:

[0026] Step 11: Preprocessing of MRA images using an improved bilateral filtering method; including the following steps:

[0027] Step 111: Dynamic range adjustment;

[0028] The low gray value part of the image is expanded by logarithmic transformation, and the transformation method is as follows:

[0029] g(x,y) = C·log[1 + f(x,y)] (1)

[0030] where g(x,y) is the gray value of the image after adjustment; C is the scale ratio coefficient; the term 1 + f(x,y) is used to avoid taking the logarithm of zero;

[0031] Step 112: Bilateral filtering;

[0032] The image after dynamic range adjustment is processed by bilateral filtering to obtain the image h(x, y). The bilateral filter BF is as follows:

[0033]

[0034] where p is the current pixel; s is the spatial range of the filtering window; I is the input image;

[0035] Gs is the spatial distance weight:

[0036]

[0037] where q is the reference pixel in the filtering window;

[0038] Gr is the pixel value weight:

[0039]

[0040] where σ r is the gray value similarity parameter, which controls the attenuation rate of the pixel value weight;

[0041] Wq is the sum of the weights of each pixel value in the filtering window, which is used for weight normalization:

[0042]

[0043] where σ s is the spatial distance parameter, which controls the attenuation rate of the spatial weight;

[0044] In the flat area, the spatial distance weight Gs in the filter dominates the filtering effect; in the edge area, the edge information will be protected;

[0045] Step 12: Extracting the anisotropic features of blood vessels; including the following steps:

[0046] Step 121: Calculating the structure tensor;

[0047]

[0048] where J_ρ(x) is the structure tensor; T is the transpose operator; ρ is the integral scale parameter; G_ρ is the Gaussian kernel function; represents the gradient vector of the image, which contains the partial derivatives in the x, y, and z directions and is used to capture the intensity change information of the blood vessel edge;

[0049] Step 122: Feature decomposition;

[0050] J_ρ(x) = λ 1 e 1 e 1 T + λ 2 e2 e 2 T + λ 3 e 3 e 3 T (7)

[0051] where λ 1 , λ 2 , λ 3 are the eigenvalues of the structure tensor, arranged from largest to smallest; e 1 , e 2 , e 3 are the eigenvectors corresponding to λ 1 , λ 2 , λ 3 respectively; e 1 serves as the principal direction vector v(x), representing the direction of the most significant local structure change and parallel to the blood vessel's orientation; e 2 and e 3 define the plane perpendicular to the blood vessel's orientation;

[0052] Step 123: Calculate the anisotropy measure;

[0053] κ(x) = (λ 1 - λ 2 ) / (λ 1 + λ 2 + λ 3 ) (8)

[0054] where κ(x) represents the anisotropy measure, used to characterize the directional intensity of the local structure, with a value range of [0, 1];

[0055] Step 13: Riemann metric construction; includes the following steps:

[0056] Step 131: Define the basic metric tensor;

[0057] M(x) = exp(-γI(x)) · [εI + (1 - ε)v(x)v(x) T (9)

[0058] where M(x) is the basic metric tensor; γ is the intensity sensitivity parameter, with a value range of (0, 1); v(x) is the blood vessel principal direction vector; I(x) represents the intensity value of the image at position x;

[0059] Step 132: Introduce the adaptive weight and curvature tensor;

[0060] The adaptive weight w(x) is as follows:

[0061] w(x) = exp(-β|κ(x) - κ o |2 ) (10)

[0062] where β is an adaptive weight parameter with a value range of (0, 1); κ(x) is the anisotropy measure; κ 0 is the reference anisotropy value;

[0063] The curvature tensor K(x) is as follows:

[0064] K(x) = ∑ i κ i (x)t i (x)t i (x) T (11)

[0065] where i is all the directions of the vector at point x; κ i (x) represents the curvature value at point x; t i (x) represents the tangent vector at point x;

[0066] Step 133: Obtain the comprehensive metric tensor;

[0067] M_final(x) = w(x)M(x) + (1 - w(x))I + αK(x) (12)

[0068] where α is the curvature weight coefficient;

[0069] Step 14: Calculate the shortest path;

[0070] Based on the constructed Riemann metric, use the improved Fast Marching method to calculate the shortest path; introduce the anisotropic Eikonal equation:

[0071]

[0072] where M_final(x) is the metric tensor; F(x) is the scalar velocity field; T(x) is the shortest arrival time function from the starting point to any point x; is the gradient of T(x), representing the shortest path direction;

[0073] Discretize Equation (13) based on difference approximation to obtain the discretized anisotropic Eikonal equation:

[0074] max{(D - x T) T M_final -1 (D - x T ),(D + x T)T M_final -1 (D + x T)} = F(x) 2 (14)

[0075] where D - x T represents the backward difference; D + x T represents the forward difference; F(x) is the scalar velocity field; meanwhile, the max operation is introduced to handle the directionality of feature propagation;

[0076] Solving equation (14) gives the shortest arrival time function T(x); then starting from the end point x, trace backward along the gradient direction of T(x) to solve:

[0077]

[0078] Obtain the parametric curve representation γ(t), t ∈ [0, 1] of the shortest path; regard the parametric curve as the approximate centerline of the blood vessel;

[0079] Step 15: Centerline optimization;

[0080] Based on the variational principle, regularization theory, and blood vessel physical constraints, use the energy minimization method to perform morphological and energy optimization on the approximate centerline of the blood vessel obtained in Step 14; the optimization function in the energy minimization method includes a data term, a smoothing constraint term, and a topological constraint term;

[0081] The data term E_data is defined as:

[0082] E_data = ∫||M_final(x)|| 2 ds (16)

[0083] where the comprehensive metric tensor M_final(x) includes the directionality of the blood vessel, image intensity, and local structure features;

[0084] The smoothing constraint term E_smooth is defined as:

[0085] E_smooth = ∫||c”(s)|| 2 ds (17)

[0086] where c”(s) is the second derivative of the curve;

[0087] The topological constraint term E_topo is defined as:

[0088] E_topo = ∑ i μ i ψ(w i ) (18)

[0089] where μ i is the adaptive weight, ψ(·) is the penalty function, and w i is the local topological feature considering the anisotropic metric;

[0090] Integrate the data term, the smoothing constraint term, and the topological constraint term to obtain the objective optimization function:

[0091] E = E_data + αE_smooth + βE_topo (19)

[0092] where α and β are the weight coefficients in the energy field, which control the intensities of the smoothing constraint and the topological constraint respectively;

[0093] In the optimization strategy, a phased optimization method with progressive constraints is selected; specifically: in the first stage, the data term is optimized, and the random gradient descent method is used to explore limitedly on the premise of ensuring the matching degree between the centerline and the actual blood vessel characteristics; in the second stage, the smoothing term is optimized, and the quasi-Newton method is used to enhance the smoothness; in the third stage, the topological constraint term is optimized to check the consistency between the centerline and the blood vessel topology.

[0094] Furthermore, the specific steps of step 2 include:

[0095] Step 21: System geometric calibration; including the following steps:

[0096] Step 211: Intrinsic parameter calibration;

[0097] The pinhole camera model is used to describe the projection characteristics of the C-arm, and the focal length f, the principal point coordinates (cx, cy), the pixel size (dx, dy), and the distortion parameters (k1, k2, p1, p2) are determined; a special calibration board containing metal marker points with known spacings is used in the calibration process; by collecting the projection images of the calibration board at different angles, the correspondence between the three-dimensional coordinates of the marker points and their two-dimensional projections is established, and the correspondence is used to construct a non-linear optimization problem:

[0098] min_{K,D}∑ i ∑ j ||π i (K,D,P j ) - p ij || 2 (20)

[0099] where K is the camera intrinsic parameter matrix, D is the distortion parameter vector, π i (·) represents the projection function at the i-th viewing angle, P j is the three-dimensional coordinate of the j-th marker point, and p ij is the corresponding observed projection coordinate;

[0100] Solve the optimization problem through the Levenberg - Marquardt algorithm to obtain the camera internal parameter matrix and the distortion parameter vector;

[0101] Step 212: Extrinsic parameter calibration;

[0102] Establish a world coordinate system and select a certain feature point on the calibration board as the origin; for each projection angle, the pose of the C - arm is described by a 4×4 rigid - body transformation matrix T:

[0103] T = [Rt; 0 1] (21)

[0104] where R is a 3×3 rotation matrix and t is a 3×1 translation vector;

[0105] Use the marked points on the calibration board to construct an optimization problem based on the reprojection error:

[0106] min_T ∑ j ||π(K, D, T, P j ) - p j || 2 (22)

[0107] where K is the camera internal parameter matrix, D is the distortion parameter vector, and p j is the three - dimensional coordinate of the j - th marked point;

[0108] Solve the optimization problem for each projection angle separately and consider the constraints of the C - arm mechanical motion to improve the accuracy of pose estimation;

[0109] Step 213: Uncertainty modeling;

[0110] For the rigid - body transformation matrix T, its uncertainty is represented by a Gaussian distribution on the Lie algebra se(3):

[0111] ξ = log(T -1 T_true) ∼ N(0, Σ_pose) (23)

[0112] where ξ is the perturbation vector on the Lie algebra se(3); T_true is the actual rigid - body transformation matrix; N is the Gaussian distribution; Σ_pose is the covariance matrix, which is estimated by collecting data from multiple repeated calibrations, analyzing the statistical characteristics of the reprojection error, and considering the kinematic constraints of the mechanical system;

[0113] Step 22: Image feature extraction; including the following steps:

[0114] Step 221: Image pre - processing;

[0115] First, a multi-scale enhancement strategy is adopted for image quality enhancement, including adaptive histogram equalization to improve local contrast; applying anisotropic diffusion filtering to suppress noise and preserve the edge information of the guide wire:

[0116]

[0117] where \(I\) represents the image grayscale value function; \(t\) represents the diffusion time parameter; represents the image gradient; \(div(\cdot)\) represents the gradient operator; \(g(x)\) represents the edge stopping function, and \(g(x)\) adopts an exponential form: \(g(x)=\exp(-x 2 / \kappa 2 ); the parameter \(\kappa\) is adaptively determined according to the statistical characteristics of the image gradient; the filtering process is iterated until the convergence criterion is met:

[0118] \(\left\|\frac{I(t + 1)-I(t)}{I(t)}\right\|<\varepsilon\quad(25)

[0119] where \(\varepsilon\) is a positive anisotropic parameter, used to determine whether the relative change between two adjacent iteration results is small enough. If the relative change is less than \(\varepsilon\), it is considered that the iteration converges and the filtering stops. Therefore, \(\varepsilon\) represents the maximum relative change threshold allowed for the image during the iteration process, reflecting the degree of preservation of image details during the filtering process;

[0120] Finally, an interested region mask is constructed based on the projection position of the MRA vascular centerline to limit the region for subsequent processing;

[0121] Step 222: Guide wire enhancement and segmentation;

[0122] In the guide wire enhancement stage, a multi-scale line structure detection method is adopted; at each scale \(\sigma\), the Hessian matrix of the image is calculated; by analyzing the eigenvalues \(\lambda 1 、\lambda 2 (|\lambda 1 |\leq|\lambda 2 |), a linear metric is constructed:

[0123] V(x,y,\sigma)=\exp\left(-\frac{\lambda 1 2}{2\alpha 2}\lambda 2 2 \right)\left(1-\exp\left(-\lambda 1 2 +\frac{\lambda 2 2}{2\beta 2}\right)\right)\quad(26)

[0124] where V(x, y, σ) is the linear structure response value at scale σ; the parameter α controls the sensitivity of line structures, and β controls the background suppression intensity; this linear metric is calculated at multiple scales and the maximum response is taken:

[0125] V(x, y) = max_{σ∈[σmin,σmax]} V(x, y, σ) (27)

[0126] where V(x, y) is the maximum value of the multi-scale linear structure response;

[0127] For guide wire segmentation, the OTSU-based method is used to adaptively determine the threshold, and then morphological operations are used to repair broken points and remove noise;

[0128] Step 223: Feature point extraction;

[0129] First, calculate the distance transform of the binary image and extract the ridges of the distance transform, and apply B-spline interpolation to achieve sub-pixel accuracy positioning; then calculate the local curvature and detect the local extreme points of the curvature, and combine the vascular anatomical structure to screen feature points, and uniformly sample to supplement feature points in the straight section; finally, perform feature point evaluation;

[0130] Step 23: Ray bundle probability model construction;

[0131] For a point p on the image plane, its corresponding ray is represented as:

[0132] R(t) = O + tD + ε_sys(t) + ε_rand(t) (28)

[0133] where R(t) is the three-dimensional spatial position of the ray at parameter t; O is the ray origin (camera optical center); D is the ray direction vector; ε_rand(t) is the random error; ε_sys(t) is the systematic error, which is directly related to pose uncertainty and intrinsic parameter uncertainty;

[0134] The error caused by pose uncertainty is expressed as:

[0135] ε_sys(t) = J_pose(t)ξ (29)

[0136] where J_pose(t) is the Jacobian matrix that describes how pose perturbations affect the ray position;

[0137] The error caused by intrinsic parameter uncertainty is expressed as:

[0138] ε_intrinsic(t) = J_intrinsic(t)η (30)

[0139] where ε_intrinsic(t) is the error caused by the uncertainty of the camera's internal parameters; J_intrinsic(t) is the Jacobian matrix of the influence of the internal parameter perturbation on the ray position; η is the internal parameter perturbation vector;

[0140] The random error ε_rand(t) is used to describe other sources of uncertainty and is modeled as distance-dependent Gaussian noise:

[0141] ε_rand(t) ~ N(0, σ_base 2 + kt 2 ) (31)

[0142] where the parameters σ_base and k are determined by analyzing the statistical distribution of the reconstruction error;

[0143] An integrated weight w is assigned to each ray bundle to improve the model reliability:

[0144] w = w_geo·w_img (32)

[0145] where w_geo is the geometric weight; w_img is the image feature weight;

[0146] Finally, each feature point corresponds to a complete probabilistic ray bundle description:

[0147] {O, D, ∑_pose, ∑_intrinsic, σ_base, k, w} (33)

[0148] where ∑_pose is the covariance matrix of the pose uncertainty, ∑_intrinsic is the covariance matrix of the internal parameter uncertainty, and w is the integrated weight;

[0149] Step 24: Vascular surface reconstruction;

[0150] First, reconstruct the vascular surface S from the DSA image, and then construct the distance field function A probability model is adopted to soften the traditional hard constraints:

[0151]

[0152] where P(x ∈ S) is the probability that the point x belongs to the vascular surface S; σ_s is the spatial uncertainty parameter;

[0153] Introduce constraints based on the normal vector and principal curvature to construct an anisotropic probability distribution:

[0154]

[0155] where \(P(x|S)\) is the conditional probability that the position, direction, and curvature features of point \(x\) satisfy the vascular morphological features under the condition of the given vascular surface \(S\); \(\sigma_r\) is the radial uncertainty parameter; \(n(x)\) is the normal vector; \(d\) is the direction vector; \(\sigma_n\) is the normal vector uncertainty parameter; \(\kappa\) 1 (x) is the principal curvature at \(x\); \(\kappa\) 2 (x) is the secondary curvature at \(x\); \(\sigma_{\kappa}\) is the curvature uncertainty parameter;

[0156] Step 25: Spatial expansion of the guidewire based on the ray bundle probability model;

[0157] Under the unified probability optimization framework, the spatial expansion problem of the guidewire is transformed into a maximum a posteriori probability estimation:

[0158] \(P(X|I,S)\propto P(I|X)P(X|S)P(X)\ (36)\)

[0159] where \(P(X|I,S)\) is the maximum a posteriori probability that the position, direction, and curvature features of point \(x\) satisfy the vascular morphological features under the conditions of the given image \(I\) and vascular surface \(S\); \(P(I|X)\) is the observation likelihood based on the ray bundle probability model; \(P(X|S)\) is the softened vascular surface constraint; \(P(X)\) is the guidewire shape prior;

[0160] After taking the logarithm of Equation (36), the energy function is obtained:

[0161]

[0162] where \(E(X)\) is the energy function, and \(X\) represents the three-dimensional coordinates of the guidewire to be solved; each term in the energy function corresponds to a specific constraint or objective in the reconstruction process:

[0163] \(\sum\) i (X i -\(\mu\) i ) T \(\sum\) i -1 (X i -\(\mu\) i ) / 2 is the projection consistency term, where \(\mu\) i represents the position predicted based on the ray bundle model; \(\sum\) i is the corresponding uncertainty matrix; \(X\) i is the three-dimensional coordinate of the \(i\)-th discrete point on the guidewire;

[0164] is the vascular surface constraint term, where is the distance constraint term, represents the distance from the point \(X\) on the guidewire i to the vascular surface; \((n(X i ) T di ) 2 / 2σ_n 2 is the direction constraint term, where n(X i ) is the normal vector of the blood vessel surface at point X i , and d i is the direction vector of the guide wire at this point. (n(X i ) T d i ) 2 ensures that the direction of the guide wire adapts to the blood vessel trend; (κ 1 (X i ) 2 +κ 2 (X i ) 2 ) / 2σ_k 2 is the curvature term, where κ 1 (X i ) is the principal curvature of the guide wire at the i-th point, and κ 2 (X i ) is the secondary curvature of the guide wire at the i-th point;

[0165] is the shape prior term;

[0166] λ 3 ||X(t)-X(t - 1)|| 2 / Δt 2 is the temporal consistency term;

[0167] By minimizing the energy function, among all possible spatial expansion results, a solution that best conforms to all known information and physical constraints is found, and this solution is the three-dimensional coordinates and morphology of the guide wire.

[0168] Furthermore, the specific steps of step 3 include:

[0169] Step 31: Initialization and discretization;

[0170] First, the guide wire is discretized into N nodes: X = {x 1 , x 2 ,..., x n}, where X represents the set of discretized nodes, and x 1 to x n represent N discrete nodes, and each node contains information on position (position coordinates in three-dimensional space) and direction (the tangent vector at this point). A piecewise cubic Hermite interpolation function is used to construct a continuous representation:

[0171] X(s) = ∑ i H i (s)x i + ∑ i H′i (s)d i (38)

[0172] where \(X(s)\) is the parametric representation of the guide wire, and \(s\) is the arc length parameter; \(\sum\) i is the corresponding uncertainty matrix; \(H\) i (s) and \(H'\) i (s) are the Hermite basis function and its derivative respectively; \(d\) i represents the tangent vector at the node; \(x\) i represents the position coordinates of the discrete nodes;

[0173] Step 32: Construct the energy function;

[0174] \(E_{total}=w\) 1 \(E_{elastic}+w\) 2 \(E_{bending}+w\) 3 \(E_{contact}+w\) 4 \(E_{data} (39)\)

[0175] where \(E_{total}\) is the weighted sum of each energy function, representing the total energy of the guide wire;

[0176] \(E_{elastic}\) is the elastic potential energy, modeled based on Hooke's law, describing the tensile and compressive deformations of the guide wire in the axial direction to ensure the physical rationality of the guide wire length:

[0177]

[0178] where \(EA\) is the axial stiffness; \(L\) is the total length of the guide wire; \(s\) is the arc length parameter;

[0179] \(E_{bending}\) is the bending energy, proportional to the square of the curvature of the guide wire, describing the lateral bending deformation of the guide wire to characterize the lateral deformation and restricting the excessive bending of the guide wire to ensure smooth shape:

[0180]

[0181] where \(EI\) is the bending stiffness;

[0182] \(E_{contact}\) is the contact energy, adopting a soft constraint and allowing small penetration, describing the interaction between the guide wire and the blood vessel wall to ensure the movement of the guide wire in the blood vessel:

[0183] \(E_{contact}=\int\) 0 Lk(s)·max(0,d_min - d(X(s),C(s))) 2 ds (42)

[0184] Where C(s) is the blood vessel centerline; k(s) is the contact stiffness coefficient that controls the contact force between the guide wire and the blood vessel wall; d_min is the minimum allowable distance between the guide wire and the blood vessel wall; d(X(s), C(s)) is the distance from the guide wire position to the catheter centerline;

[0185] E_data is a data item constructed based on the least squares criterion, which describes the consistency between the reconstruction result and the observed data and ensures that the reconstruction result conforms to the actual observation:

[0186] E_data = ∑ i w i ||X(s i ) - X_obs(s i )|| 2 (43)

[0187] Where w i is the weight coefficient of the data item; X(s i ) is the optimized guide wire position; X_obs(s i ) is the observed guide wire position;

[0188] The weight coefficients w 1 to w 4 are used to balance the relative importance of each item;

[0189] Step 33: Solve the energy function;

[0190] Solve according to the Euler - Lagrange equation derived from the variational principle of nature:

[0191]

[0192] Where ε is the strain; κ is the curvature; t and n are the tangential and normal unit vectors respectively; f_contact is the integrand of the contact energy, that is, the contact force; f_data is the summation factor in the data item summation formula, that is, the binding force generated by the data item;

[0193] Introduce an adaptive control mechanism to make the high - strain region have greater stiffness by dynamically adjusting the stiffness parameter:

[0194] EA(s) = EA 0 ·(1 + αε(s) 2 ) (45)

[0195] Where EA(s) is the position - dependent axial stiffness; EA 0 is the basic axial stiffness; αε(ε) is the strain - related stiffness adjustment factor.

[0196] Compared with the prior art, the beneficial effects of the present invention are:

[0197] The present invention proposes an innovative multi - perspective imaging space expansion and navigation method for vascular interventional surgery. This method combines probability optimization and physical constraint technologies, and improves the quality and efficiency of vascular interventional surgery through real - time, precise, and interactive navigation path planning. Specifically, it is embodied as follows:

[0198] In the preoperative preparation stage, the vascular centerline extraction technology based on MSA images is used to provide accurate reconstruction of the vascular anatomical structure, helping doctors better understand the spatial position and morphological characteristics of the lesion site.

[0199] In the intraoperative treatment stage, the probability - optimized ray - tracing technology based on beam geometry is adopted to achieve precise positioning of the guide wire and catheter, enabling doctors to master the instrument position in real - time, and thus more precisely control the surgical process. At the same time, the optimization of elastic energy minimization constraint ensures the physical rationality of the navigation path and avoids unnecessary vascular injuries. In addition, the AR / VR display system provides doctors with an immersive three - dimensional visual experience, intuitively showing the spatial relationship between the surgical instruments and blood vessels, further improving the accuracy and safety of the surgery.

[0200] In terms of postoperative effects, by improving the navigation accuracy, the risk of surgical complications such as vascular perforation and dissection is significantly reduced. At the same time, the operation time is significantly shortened, reducing the radiation exposure of patients and medical staff, and realizing a faster, safer, and more efficient surgical process. Brief Description of the Drawings

[0201] Figure 1 It is the flowchart of the method of the present invention.

[0202] Figure 2 It is the flowchart of step 1 in the present invention.

[0203] Figure 3 It is the flowchart of step 2 in the present invention. Detailed Embodiments

[0204] For a clearer understanding of the technical features, objectives, and beneficial effects of the present invention, the technical solutions of the present invention are described in detail below, but it should not be construed as a limitation on the scope of implementation of the present invention.

[0205] The following describes the specific implementation of the present invention in detail with reference to specific embodiments.

[0206] An embodiment of the present invention provides a multi - perspective imaging space expansion and navigation method for vascular interventional surgery, and its flowchart is as Figure 1 shown, and the method includes the following steps:

[0207] Step 1: Extract the vascular centerline;

[0208] In an interventional surgical navigation system, accurately extracting the vascular centerline is a crucial step to ensure navigation accuracy. Traditional DSA image processing methods often struggle to handle complex situations such as vascular overlap and uneven contrast, especially when dealing with small blood vessels and lesion areas, where both extraction accuracy and stability face severe challenges. Therefore, this method adopts a vascular centerline extraction algorithm based on anisotropic diffusion. By introducing directional information and local structural features, this algorithm effectively maintains the continuity and integrity of blood vessels, providing a reliable anatomical structure reference for subsequent 3D reconstruction and navigation optimization. In the entire navigation system, the output of this algorithm directly affects the accuracy of subsequent ray beam geometry reconstruction and the constraint conditions of elastic energy optimization, and is an important guarantee for navigation accuracy. The specific process is as follows:

[0209] Step 11: Preprocess the MRA image using an improved bilateral filtering method; including the following steps:

[0210] Step 111: Dynamic range adjustment;

[0211] The purpose of dynamic range adjustment is to improve image contrast while minimizing the loss of feature details. The low gray value part of the image is expanded through logarithmic transformation to display more details. The transformation method is as follows:

[0212] g(x,y) = C·log[1 + f(x,y)] (1)

[0213] Where g(x,y) is the gray value of the adjusted image; C is the scale factor; the term 1 + f(x,y) is to avoid taking the logarithm of zero.

[0214] Step 112: Bilateral filtering;

[0215] Bilateral filtering takes into account both distance factors and the influence of pixel value differences. The image after dynamic range adjustment is processed using bilateral filtering to obtain the image h(x,y). The bilateral filter BF is:

[0216]

[0217] Where p is the current pixel point; s is the spatial range of the filtering window; I is the input image;

[0218] Gs is the spatial distance weight:

[0219]

[0220] Where q is the reference pixel point in the filtering window;

[0221] Gr is the pixel value weight:

[0222]

[0223] where σ r is the grayscale similarity parameter, which controls the attenuation rate of the pixel value weight;

[0224] Wq is the sum of the weights of each pixel value within the filtering window and is used for weight normalization:

[0225]

[0226] where σ s is the spatial distance parameter, which controls the attenuation rate of the spatial weight;

[0227] In flat regions, the Gr values of each pixel point in the filter are similar, and the spatial distance weight Gs dominates the filtering effect. In edge regions, the Gr values on the same side of the edge are similar and much larger than the Gr values on the other side of the edge. At this time, the weights of the pixel points on the other side have little influence on the filtering result, and the edge information will be protected.

[0228] Step 12: Extracting vascular anisotropic features;

[0229] In this stage, the direction information of blood vessels is captured by constructing a local structure tensor. It includes the following steps:

[0230] Step 121: Calculating the structure tensor;

[0231]

[0232] where J_ρ(x) is the structure tensor; T is the transpose operator; ρ is the integral scale parameter, which determines the range of local feature statistics. In vascular analysis, ρ is usually chosen to be slightly larger than the diameter of the blood vessel to ensure that the direction information of the blood vessel can be accurately captured. G_ρ is the Gaussian kernel function. represents the gradient vector of the image, which contains the partial derivatives in the x, y, and z directions and is used to capture the intensity change information of the blood vessel edge. The gradient value is large at the blood vessel edge and small inside the blood vessel.

[0233] Step 122: Feature decomposition;

[0234] J_ρ(x) = λ 1 e 1 e 1 T + λ 2 e 2 e 2 T + λ 3 e 3 e 3 T (7)

[0235] where λ 1 、λ 2, λ 3 is the eigenvalue of the structure tensor, arranged from largest to smallest; e 1 , e 2 , e 3 are respectively the eigenvectors corresponding to λ 1 , λ 2 , λ 3 e1 can be used as the principal direction vector v(x), representing the direction with the most significant local structure change. This direction is usually parallel to the blood vessel's orientation to accurately describe the local direction of the blood vessel. e 2 and e 3 correspond to smaller eigenvalues, which define the plane perpendicular to the blood vessel's orientation.

[0236] Step 123: Calculate the anisotropy measure;

[0237] κ(x) = (λ 1 - λ 2 ) / (λ 1 + λ 2 + λ 3 ) (8)

[0238] where κ(x) represents the anisotropy measure, used to characterize the directional intensity of the local structure. The value range is [0, 1]. When approaching 1, it indicates strong directionality (such as a straight segment of a blood vessel), and when approaching 0, it indicates an isotropic structure (such as a blood vessel bifurcation).

[0239] Step 13: Construct the Riemann metric;

[0240] The Riemann metric originates from differential geometry and is used to describe the distance measure on a surface or manifold. In the scenario of blood vessel centerline extraction, the Riemann metric combines the directionality of the blood vessel, image intensity, and local structure features to define the "distance" in the image space, thereby more accurately extracting the blood vessel centerline. It includes the following steps:

[0241] Step 131: Define the basic metric tensor;

[0242] M(x) = exp(-γI(x)) · [εI + (1 - ε)v(x)v(x) T (9)

[0243] where \(M(x)\) is the basic metric tensor; \(\gamma\) is the intensity sensitivity parameter that controls the influence of the image intensity on the metric tensor. In angiography images, a larger \(\gamma\) value helps to highlight high-contrast regions. \(\varepsilon\) is the anisotropy parameter with a value range of \((0, 1)\) that controls the anisotropy degree of the basic metric tensor in the vessel direction and the perpendicular direction. \(v(x)\) is the main vessel direction vector. \(I(x)\) represents the intensity value of the image at position \(x\): in the angiography region, the intensity value is high; in the background region, the intensity value is low. The basic metric tensor \(M(x)\) takes into account the vessel directionality (\(v(x)v(x)\) T ) and the image intensity (\(I(x)\)).

[0244] Step 132: Introduce the adaptive weight and the curvature tensor;

[0245] The adaptive weight \(w(x)\) takes into account the local structural features of the vessels and can automatically adjust the metric strategy at straight vessels or bifurcated vessels;

[0246] \(w(x)=\exp(-\beta|\kappa(x)-\kappa 0 | 2 )\ (10)

[0247] where \(\beta\) is the adaptive weight parameter with a value range of \((0, 1)\) that controls the anisotropy degree of the metric tensor in the vessel direction and the perpendicular direction. \(\kappa(x)\) is the anisotropy metric. \(\kappa 0 is the reference anisotropy value, usually determined according to the anisotropy metric of typical vessel segments.

[0248] The curvature tensor \(K(x)\) takes into account the bending characteristics of the vessels.

[0249] \(K(x)=\sum i \kappa i (x)t i (x)t i (x) T \ (11)

[0250] where \(i\) is all directions of the vector at point \(x\); \(\kappa i (x)\) represents the curvature value at point \(x\) and describes the bending degree of the vessel at this point. \(t i (x)\) represents the tangent vector at point \(x\), and multiple \(t i (x)\) describe the bending trajectory of the vessel.

[0251] Step 133: Obtain the comprehensive metric tensor;

[0252] \(M_{final}(x)=w(x)M(x)+(1 - w(x))I+\alpha K(x)\ (12)

[0253] where α is the curvature weight coefficient that controls the influence degree of the curvature tensor K(x), and this value needs to be appropriately increased when dealing with curved blood vessels. The comprehensive metric tensor can represent the centerline of blood vessels with higher integrity, higher accuracy, and higher robustness at different scales.

[0254] Step 14: Calculate the shortest path;

[0255] Based on the constructed Riemann metric, use the improved Fast Marching method to calculate the shortest path. The original Fast Marching method solves the isotropic Eikonal equation to obtain the shortest arrival time function T(x) from the starting point to any point x, where F(x) is usually a scalar velocity field. To improve the accuracy of path selection, this method uses the metric tensor M_final(x) to replace the scalar velocity field F(x), and considers forward and backward differences as well as adding blood vessel direction information, and finally obtains the introduced anisotropic Eikonal equation:

[0256]

[0257] where is the gradient of T(x), representing the shortest path direction;

[0258] Since the computer cannot directly solve continuous partial differential equations, under the premise of maintaining numerical stability and meeting the causality requirements, equation (13) is discretized based on difference approximation for efficient calculation on the grid structure. The obtained discretized anisotropic Eikonal equation:

[0259] max{(D - x T) T M_final -1 (D - x T),(D + x T) T M_final -1 (D + x T)}=F(x) 2 (14)

[0260] where D - x T represents the backward difference, D + x T represents the forward difference, and F(x) is the scalar velocity field. At the same time, the max operation is introduced to handle the directionality of feature propagation.

[0261] Solving equation (14) can obtain the shortest arrival time function T(x). Then, starting from the end point x, trace backward along the gradient direction of T(x):

[0262]

[0263] The parametric curve representation of the shortest path γ(t), t ∈ [0, 1] is obtained. This parametric curve is regarded as the approximate centerline of the blood vessel. This step ensures that the path always advances along the optimal direction inside the blood vessel.

[0264] Step 15: Centerline optimization;

[0265] Based on the variational principle, regularization theory, and physical constraints of blood vessels, this step uses the energy minimization method to perform morphological and energy optimization on the approximate centerline of the blood vessel obtained in Step 14. The optimization function in the energy minimization method includes a data term, a smoothing constraint term, and a topological constraint term.

[0266] The data term is defined as:

[0267] E_data = ∫||M_final(x)|| 2 ds (16)

[0268] where the comprehensive metric tensor M_final(x) includes the directionality of the blood vessel, image intensity, and local structural features. The data term characterizes the matching degree between the centerline and the actual blood vessel features, ensuring that the centerline is located at the optimal position inside the blood vessel.

[0269] The smoothing constraint term is defined as:

[0270] E_smooth = ∫||c”(s)|| 2 ds (17)

[0271] where c”(s) is the second derivative of the curve. The integral form ensures global smoothness, and the quadratic form ensures global non-negativity. The smoothing constraint term is analogous to elastic energy, preventing the centerline from oscillating violently, ensuring its geometric continuity, and suppressing the influence of noise and local perturbations.

[0272] The topological constraint term is defined as:

[0273] E_topo = Σ i μ i ψ(w i ) (18)

[0274] where μ i is the adaptive weight, ψ(·) is the penalty function, and w i is the local topological feature considering the anisotropic metric. The topological constraint term maintains the topological structure of the blood vessel network, controls the connection relationship of branch points, and ensures anatomical correctness.

[0275] Integrate the data term, the smoothing constraint term, and the topological constraint term to obtain the target optimization function:

[0276] E = E_data + αE_smooth + βE_topo (19)

[0277] where α and β are weight coefficients in the energy field, controlling the strengths of the smoothing constraint and the topological constraint respectively. At the vascular bifurcation, β needs to be appropriately increased to maintain topological correctness.

[0278] In terms of the optimization strategy, a phased optimization method with progressive constraints is selected. Specifically: in the first stage, the data term is mainly optimized, and the stochastic gradient descent method is used to explore the optimization function to a limited extent on the premise of ensuring the matching degree between the centerline and the actual vascular characteristics; in the second stage, the smoothing term is mainly optimized, and the quasi - Newton method is used to enhance the smoothness; in the third stage, the topological constraint term is mainly optimized to check the consistency between the centerline and the vascular topology.

[0279] This improved algorithm for extracting vascular centerlines based on anisotropy significantly improves the extraction accuracy and robustness by introducing directional information and a multi - level optimization strategy. It can accurately capture the trend of complex vascular structures and ensure the smoothness of the centerline while maintaining topological correctness. In addition, the algorithm also has strong anti - noise ability and adaptability to vascular curvature. In practical applications, the anisotropy parameters can be adjusted according to specific vascular characteristics, and an appropriate optimization strategy can be selected in combination with the computational resource requirements.

[0280] Step 2: Perform spatial expansion on the DSA image;

[0281] In the interventional surgery navigation system, to connect the two - dimensional DSA projection image with the three - dimensional space reconstruction, a probability - optimized ray - tracing method based on beam geometry is adopted. The core of this method lies in its ability to systematically handle practical problems such as projection geometry uncertainty, image noise, and motion artifacts, thereby providing more reliable three - dimensional reconstruction results. In the entire navigation system, this algorithm receives the upstream vascular centerline information as a spatial constraint and provides an initial spatial position estimate for the downstream elastic energy optimization, which is the core link to achieve accurate navigation. The specific process is as follows:

[0282] Step 21: System geometric calibration;

[0283] In the wire three - dimensional reconstruction system based on probability - optimized ray - tracing of beam geometry, establishing an accurate geometric imaging model is the basis for achieving high - precision reconstruction. It includes the following steps:

[0284] Step 211: Intrinsic parameter calibration;

[0285] Internal parameter calibration aims to obtain the internal geometric parameters of the C-arm imaging system. The pinhole camera model is used to describe the projection characteristics of the C-arm, and the following parameters need to be determined: focal length f, principal point coordinates (cx, cy), pixel size (dx, dy), and distortion parameters (k1, k2, p1, p2). A special calibration board containing metal fiducial points with known spacings is used in the calibration process. By acquiring the projection images of the calibration board at different angles, the correspondence between the three-dimensional coordinates of the fiducial points and their two-dimensional projections is established. Using these correspondences, a non-linear optimization problem is constructed:

[0286] min_{K,D}Σ i ∑ j ||π i (K,D,P j )-p ij || 2 (20)

[0287] where K is the camera internal parameter matrix, D is the distortion parameter vector, and π i (·) represents the projection function at the i-th viewing angle, P j is the three-dimensional coordinate of the j-th fiducial point, and p ij is the corresponding observed projection coordinate.

[0288] This optimization problem can be solved by the Levenberg-Marquardt algorithm to obtain the camera internal parameter matrix and the distortion parameter vector.

[0289] Step 212: Extrinsic parameter calibration;

[0290] Extrinsic parameter calibration aims to determine the spatial pose of the C-arm at different projection angles. We establish a world coordinate system and select a certain feature point of the calibration board as the origin. For each projection angle, the pose of the C-arm can be described by a 4×4 rigid body transformation matrix T:

[0291] T = [Rt; 0 1] (21)

[0292] where R is a 3×3 rotation matrix and t is a 3×1 translation vector.

[0293] To obtain an accurate pose estimate, we also use the fiducial points on the calibration board to construct an optimization problem based on the reprojection error:

[0294] min_{T}∑ j ||π(K,D,T,P j )-p j || 2 (22)

[0295] where K is the camera internal parameter matrix, D is the distortion parameter vector, and p jis the three-dimensional coordinate of the j-th marker point;

[0296] This optimization problem needs to be solved separately for each projection angle, and the constraints of the C-arm mechanical movement, such as the isocenter rotation constraint, are considered to improve the accuracy of pose estimation.

[0297] Step 213: Uncertainty modeling;

[0298] Uncertainty modeling is the innovation of this method in establishing a complete geometric imaging model. Due to the inevitable errors in the system calibration process, these error sources are diverse, including marker point detection errors, mechanical movement errors, image noise, etc. To systematically describe these uncertainties, this method adopts a probability framework based on Lie group-Lie algebra. For the rigid body transformation matrix T, its uncertainty can be represented by a Gaussian distribution on the Lie algebra se(3):

[0299] ξ = log(T -1 T_true) ∼ N(0, ∑_pose) (23)

[0300] where ξ is the perturbation vector on the Lie algebra se(3); T_true is the actual rigid body transformation matrix; N is the Gaussian distribution; ∑_pose is the covariance matrix, which can be estimated by collecting data from multiple repeated calibrations, analyzing the statistical characteristics of the reprojection error, and considering the kinematic constraints of the mechanical system.

[0301] Step 22: Image feature extraction;

[0302] The goal of image feature extraction is to accurately extract the position information of the guidewire from the DSA sequence and quantify the related uncertainties. Since the MRA vascular centerline has been obtained, we can use this prior information to improve the accuracy and efficiency of feature extraction. The entire feature extraction process includes the following steps:

[0303] Step 221: Image preprocessing;

[0304] First, a multi-scale enhancement strategy is adopted to enhance the image quality, including adaptive histogram equalization (CLAHE) to improve the local contrast. To suppress noise while maintaining the guidewire edge information, anisotropic diffusion filtering is applied:

[0305]

[0306] where I represents the image gray value function; t represents the diffusion time parameter; represents the image gradient; div(·) represents the gradient operator; g(x) represents the edge stopping function, and g(x) adopts an exponential form: g(x) = exp(-x 2 / κ 2) The parameter κ is adaptively determined according to the statistical characteristics of the image gradient. The filtering process is iterated until the convergence criterion is met:

[0307] ||I(t + 1) - I(t)|| / ||I(t)|| < ε(25)

[0308] where ε is a positive anisotropic parameter used to determine whether the relative change between two adjacent iteration results is small enough. If the relative change is less than ε, the iteration is considered to have converged and the filtering stops. Therefore, ε represents the maximum relative change threshold allowed for the image during the iteration process, reflecting the degree of preservation of image details during the filtering process;

[0309] Finally, an interested region (ROI) mask is constructed based on the projection position of the MRA vessel centerline to limit the area for subsequent processing.

[0310] Step 222: Guidewire enhancement and segmentation;

[0311] In the guidewire enhancement stage, a multi-scale line structure detection method is adopted. At each scale σ, the Hessian matrix of the image is calculated. By analyzing the eigenvalues λ 1 、λ 2 (|λ 1 | ≤ |λ 2 |), a linear metric is constructed:

[0312] V(x, y, σ) = exp(-λ 1 2 / 2α 2 λ 2 2 )(1 - exp(-λ 1 2 + λ 2 2 / 2β 2 )) (26)

[0313] where V(x, y, σ) is the linear structure response value at scale σ; the parameter α controls the sensitivity of the line structure, and β controls the background suppression intensity. This linear metric is calculated at multiple scales and the maximum response is taken:

[0314] V(x, y) = max_{σ∈[σmin,σmax]} V(x, y, σ) (27)

[0315] where V(x, y) is the maximum value of the multi-scale linear structure response;

[0316] For guidewire segmentation, the threshold is adaptively determined using the OTSU-based method, and then morphological operations are used to repair breakpoints and remove noise to obtain a more accurate guidewire image.

[0317] Step 223: Feature point extraction;

[0318] First, calculate the distance transform of the binary image and extract the ridge of the distance transform, and apply B-spline interpolation to achieve sub-pixel accuracy positioning; then calculate the local curvature and detect the local extreme points of the curvature, and combine the vascular anatomical structure to screen the feature points, and use uniform sampling to supplement the feature points in the straight section; finally, perform feature point evaluation to ensure that they meet the requirements of subsequent processing.

[0319] Step 23: Ray bundle probability model construction;

[0320] This step unifies the system geometric uncertainty and the image feature uncertainty into the ray bundle model. For a point p on the image plane, its corresponding ray can be expressed as:

[0321] R(t) = O + tD + ε_sys(t) + ε_rand(t) (28)

[0322] where R(t) is the three-dimensional spatial position of the ray at parameter t; O is the ray origin (camera optical center); D is the ray direction vector; ε_rand(t) is the random error; ε_sys(t) is the system error, which is directly related to the pose uncertainty and the intrinsic parameter uncertainty.

[0323] The error caused by the pose uncertainty can be expressed as:

[0324] ε_sys(t) = J_pose(t)ξ (29)

[0325] where J_pose(t) is the Jacobian matrix that describes how the pose perturbation affects the ray position.

[0326] The error caused by the intrinsic parameter uncertainty can be expressed as:

[0327] ε_intrinsic(t) = J_intrinsic(t)η (30)

[0328] where ε_intrinsic(t) is the error caused by the camera intrinsic parameter uncertainty; J_intrinsic(t) is the Jacobian matrix of the influence of the intrinsic parameter perturbation on the ray position; η is the intrinsic parameter perturbation vector;

[0329] In addition, the random error ε_rand(t) is used to describe other sources of uncertainty and is modeled as distance-related Gaussian noise:

[0330] ε_rand(t) ~ N(0,σ_base 2 +kt 2 ) (31)

[0331] Among them, the parameters σ_base and k can be determined by analyzing the statistical distribution of the reconstruction error.

[0332] To improve the model reliability, a comprehensive weight w is considered to be assigned to each ray bundle:

[0333] w = w_geo · w_img (32)

[0334] Among them, the geometric weight w_geo reflects the projection angle and the spatial distribution uniformity of the rays; the image feature weight w_img is determined based on the local contrast of the image gradient magnitude, the segmentation confidence of the linear metric value, and the feature point stability of the local curvature change.

[0335] Finally, each feature point corresponds to a complete probabilistic ray bundle description:

[0336] {O, D, ∑_pose, ∑_intrinsic, σ_base, k, w} (33)

[0337] Among them, ∑_pose is the covariance matrix of the pose uncertainty, ∑_intrinsic is the covariance matrix of the intrinsic parameter uncertainty, and w is the comprehensive weight;

[0338] Step 24: Vessel surface reconstruction;

[0339] First, reconstruct the vessel surface S from the DSA image, and then construct the distance field function Considering the uncertainty in the reconstruction process, this method uses a probabilistic model to soften the traditional hard constraints:

[0340]

[0341] Among them, P(x ∈ S) is the probability that the point x belongs to the vessel surface S; σ_s is the spatial uncertainty parameter;

[0342] To better describe the local characteristics of the vessel, this method further introduces constraints based on the normal vector and the principal curvature to construct an anisotropic probability distribution:

[0343]

[0344] Among them, P(x|S) is the conditional probability that the position, direction, and curvature characteristics of the point x satisfy the vessel morphological characteristics under the condition of the given vessel surface S; σ_r is the radial uncertainty parameter; n(x) is the normal vector; d is the direction vector; σ_n is the normal vector uncertainty parameter; κ 1 (x) is the principal curvature at x; κ 2 (x) is the secondary curvature at x; σ_κ is the curvature uncertainty parameter;

[0345] By integrating these constraints, a more refined and accurate vascular surface reconstruction result is obtained.

[0346] Step 25: Spatial expansion of the guidewire based on the ray bundle probability model;

[0347] Under the unified probability optimization framework, the spatial expansion problem of the guidewire is transformed into a maximum a posteriori probability estimation:

[0348] P(X|I,S) ∝ P(I|X)P(X|S)P(X) (36)

[0349] where P(X|I,S) is the maximum a posteriori probability that the position, direction, and curvature characteristics of point x satisfy the vascular morphological characteristics under the condition of the given image I and vascular surface S; P(I|X) is the observation likelihood based on the ray bundle probability model; P(X|S) is the softened vascular surface constraint; P(X) is the prior of the guidewire shape.

[0350] To solve the problem of information loss in projecting from the known two-dimensional guidewire projection information to the unknown three-dimensional guidewire structure, this method proposes an energy minimization framework, which provides a systematic method to handle this ill-posedness. The energy function E(X) can be understood as a measure of the "goodness" of the reconstruction result, where X represents the three-dimensional coordinates of the guidewire to be solved. After logarithmizing Equation (36), the energy function is obtained:

[0351]

[0352] where E(X) is the energy function and X represents the three-dimensional coordinates of the guidewire to be solved; each term in the energy function corresponds to a specific constraint or objective in the reconstruction process:

[0353] 1) Projection consistency term: ∑ i (X i - μ i ) T ∑ i -1 (X i - μ i ) / 2; this term ensures that the reconstruction result is consistent with the two-dimensional DSA image observation. μ i represents the position predicted based on the ray bundle model; Σ i is the corresponding uncertainty matrix; X i is the three-dimensional coordinate of the i-th discrete point on the guidewire; the smaller this term is, the more consistent the reconstruction result is with the actual observation.

[0354] 2) Vascular surface constraint term: This term ensures that the reconstructed guidewire is located inside the blood vessel and follows the local geometric characteristics of the blood vessel. Among them is the distance constraint term, indicating the point X on the guide wire i to the distance from the blood vessel surface; (n(X i )) T d i ) 2 / 2σ_n 2 is the direction constraint term, n(X i ) is the normal vector of the blood vessel surface at point X i where d i is the direction vector of the guide wire at this point, (n(X i )) T d i ) 2 ensures that the direction of the guide wire is adapted to the blood vessel trend; (κ 1 (X i )) 2 +κ 2 (X i )) 2 ) / 2σ_κ 2 is the curvature term, where κ 1 (X i )) is the principal curvature of the guide wire at the i-th point, κ 2 (X i )) is the secondary curvature of the guide wire at the i-th point.

[0355] 3) Shape prior term: This term ensures the smoothness of the reconstruction result and avoids physically unrealistic mutations.

[0356] 4) Temporal consistency term: λ 3 ||X(t)-X(t - 1)|| 2 / Δt 2 ; This term ensures the continuity of the reconstruction results between consecutive frames and avoids jumps.

[0357] By minimizing this energy function, we are actually looking for an optimal three-dimensional configuration that simultaneously satisfies: conforming to the two-dimensional observation data, conforming to the blood vessel anatomical constraints, satisfying the physical properties of the guide wire, and maintaining temporal continuity. From a probabilistic perspective, energy minimization is equivalent to maximizing the posterior probability P(X|I,S). Each term in the energy function corresponds to a factor in the probability model: the projection consistency term corresponds to the observation likelihood P(I|X), the blood vessel surface constraint term corresponds to the conditional probability P(X|S), and the shape prior and temporal consistency terms correspond to the prior probability P(X). Therefore, solving the energy minimization problem is to find a solution that best conforms to all known information and physical constraints among all possible spatial extension results. This solution is the required three-dimensional coordinates and morphology of the guide wire. During the solution process, the weight parameter λ in the energy function can be adjustedi To control the relative importance of various constraints, so as to obtain the most suitable reconstruction results in different application scenarios.

[0358] This probability-optimized ray tracing algorithm based on beam geometry can not only effectively handle the uncertainty of projection geometry, but also significantly improve the stability and accuracy of reconstruction by introducing anatomical structure constraints and temporal consistency constraints. The three-dimensional spatial coordinates output by the algorithm provide a reliable initial estimate for subsequent elastic energy optimization, which is the key technical support for realizing precise navigation.

[0359] Step 3: Optimize the spatial expansion result with elastic energy minimization constraints;

[0360] In the interventional surgical navigation system, it is crucial to ensure the accuracy and physical rationality of navigation. To achieve this goal, we combine the three-dimensional coordinates of the guidewire catheter after spatial expansion with the blood vessel centerline for elastic energy minimization constraint optimization. This step aims to overcome the limitations of traditional pure geometric reconstruction methods, that is, ignoring the physical properties of the guidewire catheter, which may lead to reconstruction results violating the laws of material mechanics. By introducing elastic energy minimization constraints, we can model the guidewire catheter as an elastic body with specific physical properties, and at the same time use the blood vessel centerline as a spatial constraint to achieve a physically reasonable and anatomically accurate three-dimensional reconstruction. In the entire navigation system, this algorithm receives the upstream blood vessel centerline information and the initial spatial position estimate of ray tracing, and outputs the final navigation result, which is the guarantee for realizing precise navigation.

[0361] It should be particularly noted that the energy minimization method is used in both the extraction of the blood vessel centerline and the spatial expansion based on DSA images. However, the elastic energy minimization constraint optimization focuses more on reflecting the essential characteristics of the guidewire as a continuous elastic body and is more dedicated to ensuring that the reconstruction results conform to the physical properties of the guidewire and the blood vessel anatomical structure. These two optimization methods are complementary rather than repetitive.

[0362] The specific process is as follows:

[0363] Step 31: Initialization and discretization;

[0364] First, discretize the guidewire into N nodes: X = {x 1 , x 2 ,..., x n}, where X represents the set of discrete nodes, and x 1 to x n represent N discrete nodes, and each node contains information on position (position coordinates in three-dimensional space) and direction (tangent vector at this point). Use a piecewise cubic Hermite interpolation function to construct a continuous representation:

[0365] X(s) = ∑ i Hi (s)x i +∑ i H′ i (s)d i (38)

[0366] where X(s) is the parametric representation of the guidewire, s is the arc length parameter; ∑ i is the corresponding uncertainty matrix; H i (s) and H′ i (s) are the Hermite basis function and its derivative respectively; d i represents the tangent vector at the node; x i represents the position coordinates of the discrete nodes.

[0367] Step 32: Construction of the energy function;

[0368] E_total = w 1 E_elastic + w 2 E_bending + w 3 E_contact + w 4 E_data (39)

[0369] where E_total is the weighted sum of each energy function, representing the total energy of the guidewire;

[0370] E_elastic is the elastic potential energy, modeled based on Hooke's law, describing the tensile and compressive deformations of the guidewire in the axial direction to ensure the physical rationality of the guidewire length:

[0371]

[0372] where EA is the axial stiffness; L is the total length of the guidewire; s is the arc length parameter

[0373] E_bending is the bending energy, proportional to the square of the curvature of the guidewire, describing the transverse bending deformation of the guidewire to characterize the transverse deformation and restricting the excessive bending of the guidewire to ensure smooth shape:

[0374]

[0375] where EI is the bending stiffness;

[0376] E_contact is the contact energy, using soft constraints and allowing for small penetrations, describing the interaction between the guidewire and the vessel wall to ensure the movement of the guidewire within the blood vessel:

[0377] E_contact = ∫ 0 L k(s)·max(0, d_min - d(X(s), C(s)))2 ds (42)

[0378] where C(s) is the blood vessel centerline; k(s) is the contact stiffness coefficient that controls the contact force between the guide wire and the blood vessel wall; d_min is the minimum allowable distance between the guide wire and the blood vessel wall; d(X(s), C(s)) is the distance from the guide wire position to the catheter centerline;

[0379] E_data is a data item constructed based on the least squares criterion, which describes the consistency between the reconstruction result and the observed data and ensures that the reconstruction result conforms to the actual observation:

[0380] E_data = ∑ i w i ||X(s i ) - X_obs(s i )|| 2 (43)

[0381] where w i is the weight coefficient of the data item; X(s i ) is the optimized guide wire position; X_obs(s i ) is the observed guide wire position;

[0382] The weight coefficients w 1 to w 4 are used to balance the relative importance of each item.

[0383] Step 33: Solve the energy function;

[0384] Solve according to the Euler - Lagrange equation derived from the variational principle of nature:

[0385]

[0386] where ε is the strain; κ is the curvature; t and n are the tangential and normal unit vectors respectively; f_contact is the integrand of the contact energy, that is, the contact force; f_data is the summation factor in the data item summation formula, that is, the constraint force generated by the data item. This equation balances the internal force (the first two terms) and the external force (the last two terms).

[0387] In actual calculations, an adaptive control mechanism is also introduced to improve the stability and efficiency of the algorithm. By dynamically adjusting the stiffness parameter, the high - strain region has a greater stiffness, which can better handle large - deformation situations:

[0388] EA(s) = EA 0 ·(1 + αε(s) 2 ) (45)

[0389] where EA(s) is the position - dependent axial stiffness; EA0 is the basic axial stiffness; αε(s) is the strain-dependent stiffness adjustment factor.

[0390] This constraint optimization method based on elastic energy minimization not only ensures that the reconstruction result conforms to the material properties of the guide wire, but also guarantees the safety of navigation through contact constraints. In particular, the introduction of the adaptive control mechanism enables the algorithm to effectively handle the deformation problem of the guide wire under complex vascular morphologies.

[0391] Step 4: Optimization of human-computer interaction based on AR / VR;

[0392] The core motivation for introducing AR / VR technology for human-computer interaction optimization is to solve the spatial cognition obstacles and operation limitations brought by traditional two-dimensional display interfaces. Traditional navigation systems usually project the three-dimensional vascular structure and the positions of guide wires and catheters onto a two-dimensional plane for display. This representation method is difficult to intuitively convey complex spatial relationships, resulting in doctors having to constantly perform spatial transformation and reconstruction mentally, which undoubtedly increases the cognitive load and may affect the surgical efficiency and safety. The introduction of AR / VR technology provides doctors with a fully immersive three-dimensional visual experience, presenting complex anatomical structures and navigation information in the most natural way and significantly reducing the difficulty of spatial cognition. In the entire intraoperative navigation system, the AR / VR human-computer interaction optimization plays a key integration role. It is not just a simple display terminal, but an important link that intelligently integrates and spatially maps the outputs of Steps 1-3.

[0393] Through high-precision spatial positioning and real-time registration technologies, the system can accurately superimpose the three-dimensional reconstruction results of the vascular centerline and guide wires and catheters onto the patient's actual anatomical position. This intuitive spatial correspondence enables doctors to more accurately understand the navigation path and more precisely control the movement of instruments. At the same time, the system can also display the optimal path suggestions obtained based on elastic energy optimization in real time, providing strong support for doctors' operation decisions.

[0394] The above are only the preferred embodiments of the present invention. It should be noted that for those skilled in the art, without departing from the concept of the present invention, several modifications and improvements can still be made, and these should also be regarded as the protection scope of the present invention, and these will not affect the implementation effect of the present invention and the practicability of the patent.

Claims

1. A multi-view image space expansion and navigation method for vascular intervention surgery, characterized in that: The following steps are involved: Step 1: Extract the blood vessel centerline, including: Step 11: Preprocess the MRA images using an improved bilateral filtering method; Step 12: vascular anisotropy feature extraction; Step 13: Riemann metric construction; Step 14: Calculate the shortest path; Step 15: Centerline optimization; Step 2: Spatial expansion of the DSA image, including: Step 21: System geometry calibration; Step 22: Image feature extraction; Step 23: Building a ray beam probability model; Step 24: Blood vessel surface reconstruction; Step 25: Guidewire spatial expansion based on the beam probability model; Step 3: Perform elastic energy minimization constraint optimization on the spatial expansion results, including: Step 31: Initialization and discretization; Step 32: Energy function construction; Step 33: solving the energy function; Step 4: Optimize human-computer interaction based on AR / VR.

2. The multi-view image space expansion and navigation method for vascular intervention surgery according to claim 1, characterized in that: The specific steps of step 1 include: Step 11: Preprocess the MRA image using an improved bilateral filtering method; including the following steps: Step 111: dynamic range adjustment; The low gray value part of the image is expanded by logarithmic transformation. The transformation method is as follows: g(x,y)=C·log[1+f(x,y)] (1) Where g(x,y) is the adjusted grayscale value; C is the scale factor; the term 1+f(x,y) is used to avoid taking the logarithm of zero; Step 112: bilateral filtering; The image after dynamic range adjustment is processed by bilateral filtering to obtain the image h(x,y). The bilateral filter BF is: Where p is the current pixel; s is the spatial range of the filter window; I is the input image; Gs is the spatial distance weight: Where q is the reference pixel in the filter window; Gr is the pixel value weight: where σ r is the gray value similarity parameter, which controls the decay rate of pixel value weight; Wq is the weighted sum of each pixel value in the filter window, which is used for weight normalization: where σ s is the spatial distance parameter, which controls the decay rate of the spatial weight; In flat areas, the spatial distance weight Gs in the filter dominates the filtering effect; in edge areas, edge information will be protected; Step 12: Extracting vascular anisotropy features; including the following steps: Step 121: Calculate the structure tensor; Where J_ρ(x) is the structure tensor; T is the transpose operator; ρ is the integral scale parameter; G_ρ is the Gaussian kernel function; Represents the gradient vector of the image, including partial derivatives in the x, y, and z directions, which are used to capture the intensity change information of the blood vessel edge; Step 122: feature decomposition; J_ρ(x)=λ1e1e1 T +λ2e2e2 T +λ3e3e3 T (7) Among them, λ1, λ2, and λ3 are the eigenvalues ​​of the structure tensor, arranged from large to small; e1, e2, and e3 are the eigenvectors corresponding to λ1, λ2, and λ3 respectively; e1 is the main direction vector v(x), indicating the direction in which the local structure changes most significantly, parallel to the direction of the blood vessels; e2 and e3 define planes perpendicular to the direction of the blood vessels; Step 123: Calculate anisotropy metric; κ(x)=(λ1-λ2) / (λ1+λ2+λ3) (8) Where κ(x) represents the anisotropy metric, which is used to characterize the directional strength of the local structure and has a value range of [0,1]; Step 13: Riemann metric construction; including the following steps: Step 131: define the basic metric tensor; M(x)=exp(-γI(x))·[εI+(1-ε)v(x)v(x) T ] (9) Where M(x) is the basic metric tensor; γ is the intensity sensitivity parameter, ranging from (0,1); v(x) is the main direction vector of the blood vessel; I(x) represents the intensity value of the image at position x; Step 132: Introduce adaptive weights and curvature tensors; The adaptive weight w(x) is as follows: w(x)=exp(-β|κ(x)-κ0| 2 ) (10) Where β is the adaptive weight parameter, with a value range of (0,1); κ(x) is the anisotropy measure; κ0 is the reference anisotropy value; The curvature tensor K(x) is as follows: K(x)=∑ i κ i (x)t i (x)t i (x) T (11) where i is the direction of all vectors at point x; κ i (x) represents the curvature value at point x; t i (x) represents the tangent vector at point x; Step 133: Obtain a comprehensive metric tensor; M_final(x)=w(x)M(x)+(1-w(x))I+αK(x) (12) Where α is the curvature weight coefficient; Step 14: Calculate the shortest path; Based on the constructed Riemann metric, the improved Fast Marching method is used to calculate the shortest path; the anisotropic Eikonal equation is introduced: Where M_final(x) is the metric tensor; F(x) is the scalar velocity field; T(x) is the shortest arrival time function from the starting point to any point x; is the gradient of T(x), indicating the direction of the shortest path; The discretization of equation (13) based on the difference approximation gives the discretized anisotropic Eikonal equation: in represents reverse differential; represents forward difference; F(x) is the scalar velocity field; at the same time, the max operation is introduced to process the directionality of feature propagation; Solving equation (14) yields the shortest arrival time function T(x); then starting from the end point x, trace back along the gradient direction of T(x) to obtain the solution: Obtaining a parameter curve representation γ(t), t∈[0,1] of the shortest path; regarding the parameter curve as an approximate center line of the blood vessel; Step 15: Centerline optimization; Based on the variational principle, regularization theory and vascular physical constraints, the energy minimization method is used to optimize the morphology and energy of the vascular approximate center line obtained in step 14; the optimization function in the energy minimization method includes data terms, smoothing constraint terms and topological constraint terms; The data item E_data is defined as: E_data=∫||M_final(x)|| 2 ds (16) The comprehensive metric tensor M_final(x) contains the directionality, image intensity and local structural features of the blood vessels; The smooth constraint E_smooth is defined as: E_smooth=∫||c”(s)|| 2 ds (17) Where c”(s) is the second derivative of the curve; The topological constraint term E_topo is defined as: E_topo=Σ i m i ψ(w i ) (18) where μ i is the adaptive weight, ψ(·) is the penalty function, and w i It is the local topological feature that considers the anisotropic metric; Integrating data terms, smoothness constraints and topology constraints, we get the target optimization function: E=E_data+αE_smooth+βE_topo (19) Among them, α and β are weight coefficients in the energy field, which control the strength of smoothness constraint and topology constraint respectively; In terms of optimization strategy, a constrained progressive phased optimization method is selected; specifically: in the first stage, the data items are optimized, and the stochastic gradient descent method is used to make the optimization function explore to a limited extent while ensuring the degree of matching between the centerline and the actual vascular characteristics; in the second stage, the smoothness items are optimized, and the quasi-Newton method is used to enhance the smoothness; in the third stage, the topological constraint items are optimized to check the consistency between the centerline and the vascular topology.

3. The multi-view image space expansion and navigation method for vascular intervention surgery according to claim 1, characterized in that: The specific steps of step 2 include: Step 21: System geometry calibration; including the following steps: Step 211: internal reference calibration; The pinhole camera model is used to describe the projection characteristics of the C-arm, and the focal length f, principal point coordinates (cx, cy), pixel size (dx, dy) and distortion parameters (k1, k2, p1, p2) are determined. The calibration process uses a special calibration plate, which contains metal markers with known spacing. By collecting projection images of the calibration plate at different angles, the correspondence between the three-dimensional coordinates of the markers and their two-dimensional projections is established, and the correspondence is used to construct a nonlinear optimization problem: min_{K,D}Σ i S j ||p i (K,D,P j )-p ij || 2 (20) Where K is the camera intrinsic parameter matrix, D is the distortion parameter vector, π i (·) represents the projection function at the i-th viewing angle, P j is the three-dimensional coordinate of the jth marker point, p ij are the corresponding observation projection coordinates; Solve the optimization problem through the Levenberg-Marquardt algorithm to obtain the camera intrinsic parameter matrix and distortion parameter vector; Step 212: external parameter calibration; Establish a world coordinate system and select a feature point of the calibration plate as the origin; for each projection angle, the position of the C-arm is described by a 4×4 rigid body transformation matrix T: T=[R t;0 1] (21) Where R is a 3×3 rotation matrix and t is a 3×1 translation vector; Using the markers on the calibration board, we construct an optimization problem based on the reprojection error: min_{T}∑ j ||π(K,D,T,P j )-p j || 2 (22) Where K is the camera intrinsic parameter matrix, D is the distortion parameter vector, and p j is the three-dimensional coordinate of the jth marker point; The optimization problem is solved for each projection angle separately and the constraints of the C-arm mechanical motion are considered to improve the accuracy of pose estimation; Step 213: uncertainty modeling; For the rigid body transformation matrix T, its uncertainty is represented by the Gaussian distribution on the Lie algebra se(3): ξ=log(T -1 T_true)~N(0,Σ_pose) (23) where ξ is the perturbation vector on the Lie algebra se(3); T_true is the true rigid body transformation matrix; N is the Gaussian distribution; ∑_pose is the covariance matrix, which is estimated by collecting data from repeated calibrations, analyzing the statistical characteristics of the reprojection error, and considering the kinematic constraints of the mechanical system; Step 22: Image feature extraction; including the following steps: Step 221: image preprocessing; First, a multi-scale enhancement strategy is used to enhance the image quality, including adaptive histogram equalization to improve local contrast; anisotropic diffusion filtering is applied to suppress noise and maintain the edge information of the guidewire: Where I represents the image gray value function; t represents the diffusion time parameter; represents the image gradient; div(·) represents the gradient operator; g(x) represents the edge stop function, and g(x) is in exponential form: g(x) = exp(-x 2 / κ 2 ); the parameter κ is adaptively determined according to the statistical characteristics of the image gradient; the filtering process is iterative until the convergence criterion is met: ||I(t+1)-I(t)|| / ||I(t)||<ε(25) where ε is the specified anisotropy parameter; Finally, a region of interest mask is constructed based on the projection position of the MRA vascular centerline to limit the area for subsequent processing; Step 222: Guidewire enhancement and segmentation; The multi-scale line structure detection method is used in the guidewire enhancement stage. At each scale σ, the Hessian matrix of the image is calculated. By analyzing the eigenvalues ​​λ1 and λ2 of the Hessian matrix (|λ1|≤|λ2|), a linear metric is constructed: V(x,y,σ)=exp(-λ1 2 / 2a 2 λ2 2 )(1-exp(-λ1 2 +λ2 2 / 2b 2 )) (26) Where V(x,y,σ) is the linear structure response at scale σ; parameter α controls the sensitivity of the linear structure, and β controls the strength of background suppression; this linear metric is calculated at multiple scales and the maximum response is taken: V(x,y)=max_{σ∈[σmin,σmax]}V(x,y,σ) (27) Where V(x,y) is the maximum value of the multi-scale linear structure response; Guidewire segmentation uses an OTSU-based method to adaptively determine the threshold, and then uses graphic morphological operations to repair breakpoints and remove noise; Step 223: feature point extraction; First, the distance transform of the binary image is calculated and the ridge of the distance transform is extracted. The B-spline interpolation is applied to achieve sub-pixel precision positioning. Then, the local curvature is calculated and the local extreme point of the curvature is detected. The feature points are screened in combination with the vascular anatomical structure. Uniform sampling is used to supplement the feature points in the straight section. Finally, the feature points are evaluated. Step 23: Building a ray beam probability model; For a point p on the image plane, the corresponding ray is expressed as: R(t)=O+tD+ε_sys(t)+ε_rand(t) (28) Where R(t) is the 3D spatial position of the ray at parameter t; O is the ray starting point; D is the ray direction vector; ε_rand(t) is the random error; ε_sys(t) is the systematic error, which is directly related to the pose uncertainty and the intrinsic parameter uncertainty; The error caused by pose uncertainty is expressed as: ε_sys(t)=J_pose(t)ξ (29) Where J_pose(t) is the Jacobian matrix that describes how the pose perturbation affects the ray position; The error caused by internal parameter uncertainty is expressed as: ε_intrinsic(t)=J_intrinsic(t)η (30) Where ε_intrinsic(t) is the error caused by the uncertainty of the camera intrinsic parameters; J_intrinsic(t) is the Jacobian matrix of the influence of the intrinsic parameter perturbation on the ray position; η is the intrinsic parameter perturbation vector; The random error ε_rand(t) is used to describe other sources of uncertainty and is modeled as distance-dependent Gaussian noise: ε_rand(t)~N(0,σ_base 2 +kt 2 ) (31) The parameters σ_base and k are determined by analyzing the statistical distribution of the reconstruction error; Assign a comprehensive weight w to each ray bundle to improve model reliability: w=w_geo·w_img (32) Where w_geo is the geometric weight; w_img is the image feature weight; Finally, each feature point corresponds to a complete probabilistic ray beam description: {O,D,∑_pose,∑_intrinsic,σ_base,h,w} (33) Where ∑_pose is the covariance matrix of pose uncertainty, ∑_intrinsic is the covariance matrix of intrinsic uncertainty, and w is the comprehensive weight; Step 24: Blood vessel surface reconstruction; First, the vascular surface S is reconstructed from the DSA image, and then the distance field function is constructed Use probabilistic models to soften traditional hard constraints: Where P(x∈S) is the probability that point x belongs to the blood vessel surface S; σ_s is the spatial uncertainty parameter; Introduce constraints based on normal vectors and principal curvatures to construct anisotropic probability distribution: Where P(x|S) is the conditional probability that the position, direction and curvature characteristics of point x satisfy the vascular morphological characteristics given the vascular surface S; σ_r is the radial uncertainty parameter; n(x) is the normal vector; d is the direction vector; σ_n is the normal vector uncertainty parameter; κ1(x) is the principal curvature at x; κ2(x) is the secondary curvature at x; σ_κ is the curvature uncertainty parameter; Step 25: Guidewire spatial expansion based on the beam probability model; In a unified probabilistic optimization framework, the guidewire spatial extension problem is transformed into a maximum a posteriori probability estimation: P(X|I,S)∝P(I|X)P(X|S)P(X)(36) Where P(X|I,S) is the maximum a posteriori probability that the position, direction and curvature characteristics of point x satisfy the vascular morphological characteristics given the image I and the vascular surface S; P(I|X) is the observation likelihood based on the ray beam probability model; P(X|S) is the softened vascular surface constraint; P(X) is the guidewire shape prior; After logarithmization of equation (36), we get the energy function: Where E(X) is the energy function, X represents the three-dimensional coordinates of the guidewire to be solved; each term in the energy function corresponds to a specific constraint or goal in the reconstruction process: ∑ i (X i -μ i ) T Σ i -1 (X i -μ i ) / 2 is the projection consistency term, where μ i represents the position predicted based on the ray beam model; ∑ i is the corresponding uncertainty matrix; X i is the three-dimensional coordinate of the i-th discrete point on the guidewire; is the vessel surface constraint term, where is the distance constraint, Indicates point X on the guidewire i Distance to the blood vessel surface; (n(X i ) T d i ) 2 / 2σ_n 2 is the direction constraint term, n(X i ) is the blood vessel surface at point X i The normal vector at d i is the direction vector of the guide wire at this point, (n(X i ) T d i ) 2 Ensure that the guidewire direction is consistent with the direction of the blood vessel; (κ1(X i ) 2 +κ2(X i ) 2 ) / 2σ_κ 2 is the curvature term, where κ1(X i ) is the principal curvature of the guidewire at the i-th point, κ2(X i ) is the minor curvature of the guidewire at the i-th point; is the shape prior; is the temporal consistency term; By minimizing the energy function, a solution that best meets all known information and physical constraints is found among all possible spatial expansion results. The solution is the three-dimensional coordinates and shape of the guidewire.

4. The multi-view image space expansion and navigation method for vascular intervention surgery according to claim 1, characterized in that: The specific steps of step 3 include: Step 31: Initialization and discretization; First, the guide wire is discretized into N nodes: X = {x1, x2, ..., x n }, where X represents the discretized node set, x1 to x n Represents N discrete nodes, each of which contains position and direction information; a piecewise cubic Hermite interpolation function is used to construct a continuous representation: X(s)=Σ i H i (s)x i +S i H′ i (s)d i (38) Where X(s) is the parameterized representation of the guidewire, s is the arc length parameter; ∑ i is the corresponding uncertainty matrix; H i (s) and H′ i (s) are the Hermite basis functions and their derivatives respectively; d i represents the tangent vector at the node; x i Represents the position coordinates of discrete nodes; Step 32: Energy function construction; E_total=w1E_elastic+w2E_bending+w3E_contact+w4E_data (39) Where E_total is the weighted sum of each energy function, representing the complete energy; E_elastic is the elastic potential energy, which is modeled based on Hooke's law to describe the tensile and compressive deformation of the guidewire in the axial direction, ensuring the physical rationality of the guidewire length: Where EA is the axial stiffness; L is the total length of the guidewire; s is the arc length parameter; E_bending is the bending energy, which is proportional to the square of the curvature of the guidewire. It describes the lateral bending deformation of the guidewire and characterizes the lateral deformation, limiting the excessive bending of the guidewire and ensuring a smooth shape: Where EI is the bending stiffness; E_contact is the contact energy, which uses soft constraints to allow micro-penetration and describes the interaction between the guidewire and the vessel wall to ensure the movement of the guidewire in the vessel: E_contact=∫ o L k(s)·max(0,d_min-d(X(s),C(s))) 2 ds (42) Where C(s) is the centerline of the vessel; k(s) is the contact stiffness coefficient, which controls the contact force between the guidewire and the vessel wall; d_min is the minimum allowable distance between the guidewire and the vessel wall; d(X(s), C(s)) is the distance from the guidewire position to the centerline of the catheter; E_data is a data item, constructed based on the least squares criterion, which describes the consistency between the reconstruction result and the observed data, ensuring that the reconstruction result is consistent with the actual observation: E_data=S i w i ||X(s i )-X_obs(s i )|| 2 (43) where w i is the weight coefficient of the data item; X(s i ) is the optimized guidewire position; X_obs(s i ) is the observed guidewire position; The weight coefficients w1 to w4 are used to balance the relative importance of each item; Step 33: solving the energy function; Solve it according to the Euler-Lagrange equation derived from the variational principle of nature: Where ε is strain, which indicates the degree of local deformation; κ is curvature, which indicates the degree of local bending; t and n are tangent and normal unit vectors, respectively; f_contact is the integrated term of contact energy, i.e., contact force; f_data is the summation factor in the summation of data terms, i.e., the constraint force generated by the data terms; An adaptive control mechanism is introduced to dynamically adjust the stiffness parameters to make the high strain area have greater stiffness: EA(s)=EA0·(1+αε(s) 2 ) (45) where EA(s) is the position-dependent axial stiffness; EA0 is the basic axial stiffness; and αε(s) is the strain-dependent stiffness adjustment factor.

Images