A method and device for planning a minimally invasive radiofrequency ablation surgical path for L5 / S1 based on MR images

Abstract

The application provides a kind of L5 / S1 minimally invasive disc radiofrequency ablation surgery path planning method and device based on MR image, comprising the following steps: obtaining the MRI data of the object in DICOM format;Data is preprocessed;The path planning environment is modeled, and the data is input into the modeled model, to obtain the modeled path planning environment;According to the obtained path planning environment, the intervertebral disc protrusion to be ablated 3D sub-model is tetrahedral divided, and the effective ablation target point is selected;Candidate path is constructed;The optimal path is obtained by detecting the hard constraint condition and soft constraint condition of candidate path;The path planning environment and optimal path are output. The determined surgical path can reduce the risk of surgery and the incidence of complications while ensuring the effect of surgery. Additional cross-sectional angle, coronal plane angle, sagittal plane angle and other clinical information are provided, which improves work efficiency and decision quality.

Description

TECHNICAL FIELD

[0001] The application relates to the technical field of minimally invasive intervertebral disc radiofrequency ablation surgery, in particular to an L5 / S1 minimally invasive intervertebral disc radiofrequency ablation surgery path planning method and device based on MR images. BACKGROUND

[0002] In the field of minimally invasive medicine, the treatment technology for lumbar disc herniation is constantly evolving. Radiofrequency ablation, with its low invasiveness, fast recovery and high safety, has become an important treatment method for relieving pain caused by lumbar disc herniation. This technology precisely controls low-temperature or high-temperature energy to act on the herniated disc tissue, achieving local decompression and effectively relieving low back pain caused by nerve compression. The key to a successful operation is to accurately identify and locate the ablation target, and to safely and efficiently plan the operation path to guide the puncture needle to the target area.

[0003] However, the current radiofrequency ablation for lumbar disc herniation still faces many challenges in the implementation process. First, the planning of the operation path is highly dependent on the personal experience and subjective judgment of the doctor, which not only prolongs the operation time, but also increases the risk of operation and the incidence of complications due to lack of experience. Second, although traditional C-arm or CT-guided puncture operations can assist in positioning to some extent, they are limited by imaging accuracy and the complexity of anatomical structures, often making it difficult to achieve precise positioning of the compressed nerve root, increasing the risk of intraoperative nerve injury. In addition, individual differences in patients, variations in nerve root anatomy, and uncertainty in the location of nerve compression further increase the difficulty and uncertainty of the operation.

[0004] Therefore, Chinese patent CN118252614A proposes a transforaminal approach lumbar disc herniation radiofrequency ablation puncture path planning method based on MRI data of the object. Through 5 steps, the optimal path is finally obtained. This invention uses deep learning automatic segmentation, which can achieve the precision comparable to manual segmentation while greatly speeding up the segmentation speed. Through computer parallel computing and headless rendering, the running time of the algorithm is greatly accelerated. However, the path planning method of the radiofrequency ablation for lumbar disc herniation such as the above still has the following defects in practical use:

[0005] Firstly, in the existing path planning method for radiofrequency ablation for lumbar disc herniation, the doctor needs to operate on third-party commercial software (such as Mimics) when manually planning the path. Frequent switching between different software and standard differences affect the efficiency and accuracy of the doctor's work. And the segmentation method used in the comparative document is a general deep learning segmentation framework 3D U-Net, but its performance still has optimization space in the fine segmentation of the specific anatomical structure L5 / S1 lumbar vertebrae.

[0006] Second, the segmentation method used in the comparison documents can only display path planning results, but displays little clinical information, which seriously consumes the physician's energy.

[0007] Therefore, in view of the shortcomings of the existing technology, it is necessary to provide a method and device for L5 / S1 minimally invasive intervertebral disc radiofrequency ablation surgery path planning based on MR images to solve the shortcomings of the existing technology. Summary of the Invention

[0008] The purpose of the present invention is to avoid the shortcomings of the existing technology and provide a method and device for L5 / S1 minimally invasive intervertebral disc radiofrequency ablation surgery path planning based on MR images. This method and device have the advantages of high path planning accuracy, high visualization level, rich auxiliary decision-making information, and high execution efficiency.

[0009] In order to achieve the above-mentioned purpose, the present invention is realized by the following technical measures:

[0010] A device for minimally invasive radiofrequency ablation surgery path planning for L5 / S1 intervertebral disc based on MR images includes a data preprocessing module, a path planning environment modeling module, a target calibration module, an optimal path decision module, and a visualization module.

[0011] A method for minimally invasive radiofrequency ablation of the L5 / S1 intervertebral disc based on MR images is performed by the following steps:

[0012] Step 1: Obtain MRI data of the subject in DICOM format;

[0013] Step 2: Input the data into the data preprocessing module;

[0014] Step 3: Input the pre-processed data into the path planning environment modeling module;

[0015] Step 4: Input the modeled path planning environment into the target calibration module;

[0016] Step 5: Input the path planning environment and target into the optimal path decision module. If, after comprehensive search and analysis, the optimal path decision module fails to find any feasible path that meets the preset safety standards, it will switch to the target calibration module and recommend recalibrating the target position;

[0017] Step 6: Output the path planning environment and the optimal path to the visualization module, and display the cross-sectional angulation, coronal angulation, sagittal angulation, target coordinates, skin insertion point coordinates, puncture depth, Kambin triangle area and other auxiliary decision-making information of the optimal path.

[0018] Preferably, the above step 2 is performed by the following steps:

[0019] Step 2.1: Convert MRI from DICOM format to NIFTI format;

[0020] Step 2.2: resample the MR image to an isotropic resolution of 0.75*0.75*0.75;

[0021] Step 2.3: Normalize the grayscale values ​​of the MR images to the range of [0, 1].

[0022] Step 2.4: Input the processed NIFTI format MRI data into the path planning environment modeling module;

[0023] Preferably, the above step 3 is performed by the following steps:

[0024] Step 3.1. Input the preprocessed MR images into the WT-UNet model to segment all target structures (L5 vertebra, S1 vertebra, intervertebral disc, nerve root, dura mater, skin, ilium, and intervertebral disc herniation).

[0025] Step 3.2: Use the final model trained in step 3.1 to segment all objects in the MR image, and apply morphological operations (such as hole filling and removal of isolated false positive regions) to the model output to optimize the segmentation results.

[0026] Step 3.3: Use the Marching Cubes algorithm to reconstruct the binary images of each segmented structure into a 3D image. Smooth the reconstructed 3D model with 4 iterations and an iteration factor of 0.4 to complete the path planning environment modeling.

[0027] The content of the model training phase in step 3.1 is specifically as follows:

[0028] Step 3.1.1: Construct a 3D WT-UNet, using the 3D U-Net framework as the base and integrating the wavelet transform and inverse transform modules. The encoding and decoding parts of the model are divided into five stages of feature layers. Each stage of the encoder / decoder contains two convolutional structures, and the convolution stride is set to achieve pooling and upsampling operations.

[0029] Step 3.1.2: Replace all up / down sampling parts in the basic model in 3.1.1 with 3D inverse wavelet transform / wavelet transform to realize up / down sampling function and achieve lossless transmission of information.

[0030] Step 3.1.3: The model uses the cross entropy loss function and the soft threshold dice loss function as the total loss function for training the model, uses the Stochastic Gradient Descent optimizer for model training, and enables the deep supervision strategy during the model training phase;

[0031] Step 3.1.4: The initial learning rate for model training is set to 1e-2 and decreases with the number of training rounds, for a total of 1000 rounds. Five-fold cross-validation is used during the training process to select the optimal model parameters.

[0032] The specific contents of the wavelet transform / inverse wavelet transform in step 3.1.2 are as follows:

[0033] Step 3.1.2.1. Considering that the three-dimensional wavelet transform will increase the feature channels of the model by eight times, the image is subjected to one-dimensional wavelet transform in the XYZ axis direction to reduce the calculation time.

[0034] Step 3.1.2.2: Assume the model has C feature channels. For the multiplied feature channel number 8C, use a convolution module to reduce the channel number to C. Due to the characteristics of the wavelet transform, the original image is reduced to half the size of the original image after a wavelet transform, thus having the same effect as the traditional downsampling operation in terms of image size.

[0035] In step 3.1.2.3, for the inverse wavelet transform, the number of model feature channels, C, is increased to 8C through a convolution module before the inverse wavelet transform is performed. Similarly, in step 3.1.2.2, a one-dimensional inverse wavelet transform is performed on the image along the X, Y, and Z axes to reduce computational time. The downsampled image size after the inverse wavelet transform is restored to its pre-wavelet transform size, with the number of channels restored from 8C to C.

[0036] Preferably, the above step 4 is performed by the following steps:

[0037] Step 4.1: Tetrahedron-divide the 3D sub-model of the intervertebral disc herniation to be ablated.

[0038] Step 4.2: The doctor selects the ablation target by scaling, rotating, and moving the 3D sub-model to be ablated and using the cursor-assisted positioning function.

[0039] Step 4.3: After the doctor preliminarily selects the target, an estimated range of the ablation area with adjustable transparency is generated at the target (a sphere model with adjustable radius, an ellipsoid model with adjustable major and minor axes).

[0040] Step 4.4: The doctor evaluates the legitimacy of the target selection based on the visualized effective ablation area and confirms whether to define the currently selected point as the ablation target.

[0041] In step 4.1, the tetrahedron decomposition is specifically as follows:

[0042] Step 4.1.1: Obtain all vertices of the 3D model to form an initial point set.

[0043] Step 4.1.2. Construct a super tetrahedron large enough to contain all the initial points.

[0044] Step 4.1.3. Randomly insert the initial points into the super regular tetrahedron one by one. For each newly inserted point, traverse to find all tetrahedrons whose known circumscribed spheres contain the point (called influencing tetrahedrons). Delete the common faces of the influencing tetrahedrons and connect the new points with all vertices of the influencing tetrahedrons to form new tetrahedrons. Perform empty circumscribed sphere detection on the newly generated tetrahedrons to ensure that they meet the empty circumscribed sphere characteristics of the Delaunay tetrahedron partitioning. If a tetrahedron that does not meet the empty circumscribed sphere characteristics is detected, local optimization is required. This usually involves exchanging the diagonals of the tetrahedrons to form new tetrahedrons until all tetrahedrons meet the Delaunay conditions.

[0045] Step 4.1.4: Randomly divide the points on the edges, constraint lines, and constraint surfaces of the initial tetrahedron obtained in step 4.1.3. The distance between the sampled points and the known points cannot be less than 1, otherwise they will be discarded. Perform point interpolation according to step 4.1.3.

[0046] Step 4.1.5: Randomly sample points in the tetrahedron space obtained in step 4.1.4. The distance between the sampled points and the known points cannot be less than 1, otherwise they will be discarded. Then perform point interpolation according to step 4.1.3.

[0047] Step 4.1.6. Delete all tetrahedrons associated with the vertices of the super tetrahedron. The remaining tetrahedrons constitute the final Delaunay tetrahedron.

[0048] The construction of the super tetrahedron in step 4.1.2 is specifically as follows:

[0049] Step 4.1.2.1. Calculate the bounding box of the initial point set. The bounding box is a cuboid whose sides are parallel to the x, y, and z axes, and which can exactly contain all given spatial points. That is, the coordinates of each spatial point are located inside or on the surface of the cuboid.

[0050] Step 4.1.2.2: Construct a sphere with the center of the bounding box as the center and a diameter R equal to four times the length of the bounding box's diagonal.

[0051] Step 4.1.2.3, construct a regular tetrahedron. The side length of the regular tetrahedron is And the ball in step 4.1.2.2 is the inscribed ball of the regular tetrahedron.

[0052] In step 4.2, the cursor-assisted positioning is specifically as follows:

[0053] Step 4.2.1: Model transformation: transform the vertices of the 3D sub-model to be ablated from model space to world space to obtain the world coordinates of the vertex.

[0054] worldPosition=modelMatrix·localPosition;

[0055] Where worldPosition is the coordinate (vector) of the model vertex in world space. modelMatrix is ​​the model matrix (4×4 matrix) containing scale, rotation and translation information. localPosition is the coordinate (vector, usually expressed as homogeneous coordinates, i.e., adding an extra 1 as the fourth component) of the model vertex in model space (local space).

[0056] Step 4.2.2: View transformation, convert the 3D model from world space to observation space (camera space).

[0057] eyePosition=viewMatrix·worldPosition;

[0058] Where eyePosition is the coordinate (vector) of the model's vertex in the observation space (camera space). viewMatrix is ​​the view matrix (4×4 matrix) that contains the camera's position, orientation, and "up" direction information.

[0059] Step 4.2.3, projection transformation, converts the 3D model from the observation space to the clipping space, and finally maps it to the normalized device coordinate (NDC) space.

[0060] clipPosition=projectionMatrix·eyePosition;

[0061] Where projectionMatrix is ​​the projection matrix (4×4 matrix), which is a perspective projection.

[0062] Step 4.2.4: Clipping removes vertices outside the view frustum, which are not visible in the final image. Clipping is done by comparing with the six clipping planes (left, right, top, bottom, near, and far) of the clipping space and using the Sutherland-Hodgman vertex clipping algorithm.

[0063] Step 4.2.5, perspective division, converts the clip space coordinates to normalized device coordinates (NDC).

[0064]

[0065] Where clipPosition is the coordinate (vector) of the model vertex in clip space. ndcPosition is the coordinate (vector) of the model vertex in normalized device coordinate (NDC) space.

[0066] Step 4.2.6, Rasterization, converts geometric primitives (such as triangles) into pixels (or fragments) on the screen. Vertices are converted into fragments that cover the pixel area of ​​the screen, and corresponding attributes (such as color, depth, etc.) are generated for each fragment.

[0067] Step 4.2.7, depth test, is used to perform pixel discarding. Compare the fragment's depth value with the depth value at the corresponding position in the frame buffer. If the fragment's depth value is less than or equal to the value in the frame buffer, the decision is discarded.

[0068] Step 4.2.8, viewport transformation, maps normalized device coordinates (NDC) to screen coordinates.

[0069]

[0070] Where ScreenX is the horizontal screen coordinate of the point, and ScreenY is the vertical screen coordinate of the point. Where width and height are the width and height of the viewport, and viewportX and viewportY are the coordinates of the lower-left corner of the viewport on the screen (if the viewport fills the entire screen, these values ​​are usually 0 and 0).

[0071] Step 4.2.9. Get the screen coordinates of the mouse cursor.

[0072] Step 4.2.10. Calculate the distance between the mouse cursor screen coordinates and the 3D model vertex screen coordinates. If there is a vertex whose Euclidean distance to the cursor screen coordinates is less than 2, move the cursor to that point. If the number of vertices whose Euclidean distance is less than 2 is greater than 1, move the cursor to the nearest point.

[0073] Preferably, the above step 5 is performed by the following steps:

[0074] The ablation targets obtained in step 5.1 and step 4 are input into the optimal path decision module.

[0075] Step 5.2: Sample the surface of the peritoneal side of the 3D skin model to obtain candidate needle insertion points.

[0076] Step 5.3: The candidate path is abstracted as the connection between the ablation target and each candidate needle entry point.

[0077] Step 5.4: Check the hard constraints on the candidate paths and delete the paths that do not meet the hard constraints.

[0078] Step 5.5: Calculate the soft constraint values ​​of the candidate paths that meet the hard constraint conditions.

[0079] Step 5.6: Determine the optimal path based on the soft constraint values ​​of the candidate paths that meet the hard constraints.

[0080] The specific process of obtaining candidate needle insertion points in step 5.2 is as follows:

[0081] Step 5.2.1. Crop the 3D skin model, retaining only the local area on the peritoneal side.

[0082] Step 5.2.2: After cropping, all vertices of the triangular mesh of the peritoneal surface of the 3D skin model are added to the candidate needle insertion point set.

[0083] Step 5.2.3, calculate the area of ​​the triangle and the total area. For each triangle, assuming its vertices are A, B, C, then the area S i ,

[0084]

[0085] Total area S sum ,

[0086]

[0087] Where AB=AB, AC=AC.

[0088] Step 5.2.4: Generate a random number r, satisfying 0<=r<1.

[0089] Step 5.2.5. Select a triangular face. Traverse all faces in order and calculate their cumulative area ratio P i ,

[0090]

[0091] Find the one that satisfies P i-1 <=r <P i , then the i-th face is selected.

[0092] Step 5.2.6. Generate random numbers r1 and r2 in the range [0,1] and ensure that their sum is less than 1.

[0093] P j =P A +r1*Vector AB+r2*Vector AC.

[0094] Step 5.2.7. Repeat step 5.2.6 N times to generate N sampling points and add them to the candidate needle insertion point set.

[0095] N=5*S j ;

[0096] Step 5.2.8. Repeat steps 5.2.4 to 5.2.7 10,000 times.

[0097] In step 5.4, the hard constraint condition detection of the candidate path is specifically as follows:

[0098]

[0099] In step 5.5, the calculation of the soft constraint condition value of the candidate path is specifically as follows:

[0100]

[0101] The optimal path decision in step 5.6 is specifically as follows:

[0102] Step 5.6.1. Calculate the soft constraint values ​​of each candidate path and normalize them. Items a, c, and d are calculated using the following formula:

[0103]

[0104] Item b is as follows:

[0105]

[0106] Step 5.6.2: Calculate the risk value of each candidate path. The formula is as follows:

[0107]

[0108] ∑α i +β+r+δ=1;

[0109] Among them, R P is the total risk value of the puncture path, and its value range is [0,1]. The closer the value is to 0, the smaller the risk of the path is, and the closer it is to 1, the greater the risk of the path is. i is the key structure distance risk sub-objective function, a total of 5 key structures are included in the calculation, α i is the weight coefficient corresponding to the key structure distance. S P L is the Kambin triangle projection area risk sub-objective function, and β is the corresponding weight coefficient. P is the puncture depth risk sub-objective function, and γ is the corresponding weight coefficient. p is the angle between the path and the surface of the intervertebral disc to be ablated, and δ is the corresponding weight coefficient.

[0110] Preferably, the above step 5.6.2 is performed by the following steps:

[0111] Step 5.6.2.1. For the n evaluation indicators, denote them as A_1, A_2, ..., A_n.

[0112] Step 5.6.2.2: Construct a priority graph for each pair of indicators A_i and A_j (i≠j, 1≤i, j≤n).

[0113] If A_i is more important than A_j, mark the arrow from A_i to A_j in the priority diagram, and denote it as a_{ij}=1.

[0114] If A_j is more important than A_i, mark the arrow from A_j to A_i in the priority diagram, and denote it as a_{ij}=0.

[0115] If A_i and A_j are of equal importance, no arrow is drawn, but usually to avoid confusion, they can be regarded as a_{ij} = 0.5, and a_{ji} = 0.5 must be ensured.

[0116] Note: a_{ii} is always undefined or set to 0, since it is not compared to itself.

[0117] Step 5.6.2.3. Calculate the importance score. For each indicator A_i, calculate its importance score S_i:

[0118] S_i=∑(a_{ij}-a_{ji}), where j=1 to n and j≠i

[0119] This formula calculates the "net importance" of indicator Ai relative to all other indicators (taking into account comparisons in both positive and negative directions).

[0120] Step 5.6.2.4, normalize the weights. Calculate the sum of all indicator importance scores:

[0121] S_total=∑S_i

[0122] For each indicator A_i, calculate its normalized weight W_i:

[0123] W_i=S_i / S_total.

[0124] Advantages compared to existing technologies:

[0125] First, the improved WT-UNet deep neural segmentation network not only inherits the powerful feature extraction and contextual information integration capabilities of the 3D U-Net, but also creatively incorporates a wavelet transform module to alleviate information loss during upsampling and downsampling operations. Compared to the 3D U-Net, it achieves improvements in DICE coefficient, precision, and recall for both bony and non-bone tissue segmentation at the L5 / S1 level, laying the foundation for an accurate and reliable path planning environment.

[0126] The target calibration module's tetrahedral segmentation of the 3D model of the herniated intervertebral disc to be ablated and cursor-assisted positioning assist physicians in ensuring accurate identification and positioning of the target. At the same time, the preset ablation area estimation model with adjustable transparency assists physicians in judging the rationality of the target, reducing the subjectivity and uncertainty of human judgment.

[0127] The optimal path decision module not only considers the clinical constraints of conventional path planning but also incorporates the Kambin triangle area of ​​the L5 / S1 segment. Based on this comprehensive consideration of various factors, a priority diagram is constructed, converting the physician's subjective experience into indicator weights. The physician can then fine-tune these weights based on actual circumstances, thereby determining a surgical path that minimizes surgical risks and complications while ensuring surgical effectiveness.

[0128] Secondly, this patent proposes an integrated, high-performance, one-stop path planning solution by integrating a data preprocessing module, a path planning environment modeling module, a target calibration module, an optimal path decision module, and a visualization module. This solves the problem of physicians frequently switching between multiple software programs when planning paths. It also provides additional clinical information such as transverse, coronal, and sagittal angulation, improving overall work efficiency and decision-making quality. BRIEF DESCRIPTION OF THE DRAWINGS

[0129] Figure 1 This is a flow chart of a method for planning a surgical path for minimally invasive L5 / S1 intervertebral disc radiofrequency ablation based on MR images according to the present invention;

[0130] Figure 2 This is a schematic diagram of the system structure of an MR image-based L5 / S1 minimally invasive intervertebral disc radiofrequency ablation surgical path planning device of the present invention;

[0131] Figure 3 Schematic diagram of the 3D WT-UNet model structure of the present invention;

[0132] Figure 4 This is a schematic diagram of the path planning environment of the present invention;

[0133] Figure 5 Schematic diagram of the structure of the 3D sub-model of the intervertebral disc herniation to be ablated which is divided into tetrahedra according to the present invention;

[0134] Figure 6 This is an example diagram of candidate paths and optimal paths of the present invention; the black line in the figure is the candidate path, and the white line is the optimal path;

[0135] Figure 7 This is a schematic diagram of the indicator weight diagram of the present invention;

[0136] Figure 8 This is an example diagram of the optimal path of the present invention;

[0137] Figure 9 Comparison of DSC results segmented by different models. DETAILED DESCRIPTION

[0138] Example 1, refer to the attached Figures 1-8 A method for minimally invasive radiofrequency ablation of L5 / S1 intervertebral disc based on MR images includes the following steps:

[0139] Step 1: Obtain the MRI data of the subject in DICOM format.

[0140] Step 2: preprocess the data;

[0141] Preprocessing includes format conversion, resampling, and grayscale value unification, specifically including the following steps:

[0142] Step 2.1: Convert MRI from DICOM format to NIFTI format;

[0143] Step 2.2: resample the MR image to an isotropic resolution of 0.75*0.75*0.75;

[0144] Step 2.3: Normalize the grayscale values ​​of the MR images to the range of [0, 1].

[0145] Step 2.4: Obtain the processed MRI data in NIFTI format.

[0146] Step 3: Model the path planning environment and input the pre-processed data in step 2 into the model to obtain the modeled path planning environment;

[0147] The details of step 3 are as follows: Step 3.1, input the preprocessed MR image into the WT-UNet model to segment all target structures (L5 vertebra, S1 vertebra, intervertebral disc, nerve root, dura mater, skin, ilium, intervertebral disc herniation).

[0148] The content of the model training phase in step 3.1 is specifically as follows:

[0149] Step 3.1.1: Construct a 3D WT-UNet, using the 3D U-Net framework as the base and integrating the wavelet transform and inverse transform modules. The encoding and decoding parts of the model are divided into five stages of feature layers. Each stage of the encoder / decoder contains two convolutional structures, and the convolution stride is set to achieve pooling and upsampling operations.

[0150] Step 3.1.2: Replace all up / down sampling parts in the basic model in 3.1.1 with 3D inverse wavelet transform / wavelet transform to realize up / down sampling function and achieve lossless transmission of information.

[0151] Step 3.1.3: The model uses the cross entropy loss function and the soft threshold dice loss function as the total loss function for training the model, uses the Stochastic Gradient Descent optimizer for model training, and enables the deep supervision strategy during the model training phase;

[0152] Step 3.1.4: The initial learning rate for model training is set to 1e-2 and decreases with the number of training rounds, for a total of 1000 rounds. Five-fold cross-validation is used during the training process to select the optimal model parameters.

[0153] The specific contents of the wavelet transform / inverse wavelet transform in step 3.1.2 are as follows:

[0154] Step 3.1.2.1. Considering that the three-dimensional wavelet transform will increase the feature channels of the model by eight times, the image is subjected to one-dimensional wavelet transform in the XYZ axis direction to reduce the calculation time.

[0155] Step 3.1.2.2: Assume the model has C feature channels. For the multiplied feature channel number 8C, use a convolution module to reduce the channel number to C. Due to the characteristics of the wavelet transform, the original image is reduced to half its original size after a wavelet transform, thus having the same effect as traditional downsampling in terms of image size.

[0156] In step 3.1.2.3, for the inverse wavelet transform, the number of model feature channels, C, is increased to 8C through a convolution module before the inverse wavelet transform is performed. Similarly, in step 3.1.2.2, a one-dimensional inverse wavelet transform is performed on the image along the X, Y, and Z axes to reduce computational time. The downsampled image size after the inverse wavelet transform is restored to its pre-wavelet transform size, with the number of channels restored from 8C to C.

[0157] Step 3.2: Use the final model trained in 3.1 to segment all targets in the MR image, and perform morphological operations (such as hole filling and removal of isolated false positive areas) on the output of the model to optimize the segmentation results.

[0158] Step 3.3: Use the Marching Cubes algorithm to reconstruct the binary images of each segmented structure into a 3D image. Smooth the reconstructed 3D model with 4 iterations and an iteration factor of 0.4 to complete the path planning environment modeling.

[0159] Step 4: Based on the path planning environment obtained in step 3, the 3D sub-model of the intervertebral disc herniation to be ablated is tetrahedroned to select an effective ablation target; that is, the ablation target is obtained.

[0160] The “selection of effective ablation targets” includes the following steps:

[0161] Step 4.1: Tetrahedron-divide the 3D sub-model of the intervertebral disc herniation to be ablated.

[0162] In step 4.1, the tetrahedron decomposition is specifically as follows:

[0163] Step 4.1.1: Obtain all vertices of the 3D model to form an initial point set.

[0164] Step 4.1.2. Construct a super tetrahedron large enough to contain all the initial points.

[0165] The construction of the super tetrahedron in step 4.1.2 is specifically as follows:

[0166] Step 4.1.2.1. Calculate the bounding box of the initial point set. The bounding box is a cuboid whose sides are parallel to the x, y, and z axes, and which can exactly contain all given spatial points. That is, the coordinates of each spatial point are located inside or on the surface of the cuboid.

[0167] Step 4.1.2.2: Construct a sphere with the center of the bounding box as the center and a diameter R equal to four times the length of the bounding box's diagonal.

[0168] Step 4.1.2.3, construct a regular tetrahedron. The side length of the regular tetrahedron is And the ball in step 4.1.2.2 is the inscribed ball of the regular tetrahedron.

[0169] Step 4.1.3. Randomly insert the initial points into the super regular tetrahedron one by one. For each newly inserted point, traverse to find all tetrahedrons whose known circumscribed spheres contain the point (called influencing tetrahedrons). Delete the common faces of the influencing tetrahedrons and connect the new points with all vertices of the influencing tetrahedrons to form new tetrahedrons. Perform empty circumscribed sphere detection on the newly generated tetrahedrons to ensure that they meet the empty circumscribed sphere characteristics of the Delaunay tetrahedron partitioning. If a tetrahedron that does not meet the empty circumscribed sphere characteristics is detected, local optimization is required. This usually involves exchanging the diagonals of the tetrahedrons to form new tetrahedrons until all tetrahedrons meet the Delaunay conditions.

[0170] Step 4.1.4: Randomly divide the points on the edges, constraint lines, and constraint faces of the initial tetrahedron obtained in step 4.1.3. The Euclidean distance between the sampled points and the known points must be at least 1, otherwise they will be discarded. Perform point interpolation as in step 4.1.3.

[0171] Step 4.1.5: Randomly sample points within the tetrahedron space obtained in step 4.1.4. The Euclidean distance between the sampled points and the known points must not be less than 1, otherwise they will be discarded. Then, perform point interpolation according to step 4.1.3.

[0172] Step 4.1.6. Delete all tetrahedrons associated with the super tetrahedron vertices. The remaining tetrahedrons constitute the final Delaunay tetrahedron.

[0173] Step 4.2: The doctor selects the ablation target by scaling, rotating, and moving the 3D sub-model to be ablated and using the cursor-assisted positioning function.

[0174] In step 4.2, the cursor-assisted positioning is specifically as follows:

[0175] Step 4.2.1: Model transformation: transform the vertices of the 3D sub-model to be ablated from model space to world space to obtain the world coordinates of the vertex.

[0176] worldPosition=modelMatrix·localPosition;

[0177] Where worldPosition is the coordinate (vector) of the model vertex in world space. modelMatrix is ​​the model matrix (4×4 matrix) containing scale, rotation and translation information. localPosition is the coordinate (vector, usually expressed as homogeneous coordinates, i.e., adding an extra 1 as the fourth component) of the model vertex in model space (local space).

[0178] Step 4.2.2: View transformation, convert the 3D model from world space to observation space (camera space).

[0179] eyePosition=viewMatrix·worldPosition;

[0180] Where eyePosition is the coordinate (vector) of the model's vertex in the observation space (camera space). viewMatrix is ​​the view matrix (4×4 matrix) that contains the camera's position, orientation, and "up" direction information.

[0181] Step 4.2.3, projection transformation, converts the 3D model from the observation space to the clipping space, and finally maps it to the normalized device coordinate (NDC) space.

[0182] clipPosition=projectionMatrix·eyePosition;

[0183] Where projectionMatrix is ​​the projection matrix (4×4 matrix), which is a perspective projection.

[0184] Step 4.2.4: Clipping removes vertices outside the view frustum, which are not visible in the final image. Clipping is done by comparing with the six clipping planes (left, right, top, bottom, near, and far) of the clipping space and using the Sutherland-Hodgman vertex clipping algorithm.

[0185] Step 4.2.5, perspective division, converts the clip space coordinates to normalized device coordinates (NDC).

[0186]

[0187] Where clipPosition is the coordinate (vector) of the model vertex in clip space. ndcPosition is the coordinate (vector) of the model vertex in normalized device coordinate (NDC) space.

[0188] Step 4.2.6, Rasterization, converts geometric primitives (such as triangles) into pixels (or fragments) on the screen. Vertices are converted into fragments that cover the pixel area of ​​the screen, and corresponding attributes (such as color, depth, etc.) are generated for each fragment.

[0189] Step 4.2.7, depth test, is used to perform pixel discarding. Compare the fragment's depth value with the depth value at the corresponding position in the frame buffer. If the fragment's depth value is less than or equal to the value in the frame buffer, the decision is discarded.

[0190] Step 4.2.8, viewport transformation, maps normalized device coordinates (NDC) to screen coordinates.

[0191]

[0192] Where ScreenX is the horizontal screen coordinate of the point, and ScreenY is the vertical screen coordinate of the point. Where width and height are the width and height of the viewport, and viewportX and viewportY are the coordinates of the lower-left corner of the viewport on the screen (if the viewport fills the entire screen, these values ​​are usually 0 and 0).

[0193] Step 4.2.9. Get the screen coordinates of the mouse cursor.

[0194] Step 4.2.10. Calculate the distance between the mouse cursor screen coordinates and the 3D model vertex screen coordinates. If there is a vertex whose Euclidean distance to the cursor screen coordinates is less than 2, move the cursor to that point. If the number of vertices whose Euclidean distance is less than 2 is greater than 1, move the cursor to the nearest point.

[0195] Step 4.3: After the doctor initially selects the target, an estimated ablation area with adjustable transparency is generated at the target (a sphere model with adjustable radius, an ellipsoid model with adjustable major and minor semi-axes);

[0196] Step 4.4: The doctor evaluates the legitimacy of the target selection based on the visualized effective ablation area to confirm that the currently selected point is designated as the ablation target.

[0197] Step 5. Perform hard constraint checks and soft constraint checks on the ablation target in step 4 and the candidate entry points obtained by surface sampling on the peritoneal side of the 3D skin model to obtain the optimal path. That is, a candidate path is constructed based on the "ablation target in step 4" and the "candidate entry points obtained by surface sampling on the peritoneal side of the 3D skin model." The candidate path is abstracted as a line connecting the ablation target and each candidate entry point. Perform hard constraint checks and soft constraint checks on the candidate path to obtain the optimal path.

[0198] The step 5 is performed by the following steps:

[0199] Step 5.1, input the ablation target obtained in step 4;

[0200] Step 5.2: Sampling the surface of the peritoneal side of the 3D skin model to obtain candidate needle insertion points;

[0201] The specific process of obtaining candidate needle insertion points in step 5.2 is as follows:

[0202] Step 5.2.1. Crop the 3D skin model, retaining only the local area on the peritoneal side.

[0203] Step 5.2.2: After cropping, all vertices of the triangular mesh of the peritoneal surface of the 3D skin model are added to the candidate needle insertion point set.

[0204] Step 5.2.3, calculate the area of ​​the triangle and the total area. For each triangle, assuming its vertices are A, B, C, then the area S i ,

[0205]

[0206] Total area S sum ,

[0207]

[0208] Where AB=AB, AC=AC.

[0209] Step 5.2.4: Generate a random number r, satisfying 0<=r<1.

[0210] Step 5.2.5. Select a triangular face. Traverse all faces in order and calculate their cumulative area ratio P i ,

[0211]

[0212] Find the one that satisfies P i-1 <=r <P i , then the i-th face is selected.

[0213] Step 5.2.6. Generate random numbers r1 and r2 in the range [0,1] and ensure that their sum is less than 1.

[0214] P j =P A +r1*Vector AB+r2*Vector;

[0215] Step 5.2.7. Repeat step 5.2.6 N times to generate N sampling points and add them to the candidate needle insertion point set.

[0216] N=5*Si;

[0217] Step 5.2.8. Repeat steps 5.2.4 to 5.2.7 for a total of 50,000 times.

[0218] Step 5.3: The candidate path is abstracted as the connection between the target and each candidate injection point.

[0219] Step 5.4: Check the hard constraints on the candidate paths and delete those that do not meet the hard constraints.

[0220] Hard constraints include the following:

[0221] A. Avoid risky tissue structures in the L5-S1 segment (L5 vertebra, S1 vertebra, nerve roots, dura mater, ilium), avoid collision, and keep a safety margin of more than 0.5 mm;

[0222] B. The puncture distance is less than 100 mm of the ablation needle length;

[0223] C. The entry point of the puncture path must be located on the skin of the L5-S1 segment to limit the insertion range of the ablation needle;

[0224] D. To prevent the ablation needle from sliding on the intervertebral disc, the angle between the puncture path and the plane where the triangular facets of the 3D intervertebral disc model intersect should be greater than 20 degrees;

[0225] E. The puncture path must pass through the intervertebral foramen;

[0226] F. The puncture must be performed from the back.

[0227] Step 5.5: Calculate the soft constraint values ​​of the candidate paths that meet the hard constraint conditions.

[0228] In step 5.5, the calculation of the soft constraint condition value of the candidate path is specifically as follows:

[0229] A. Distances between the puncture path and key tissue structures (L5 vertebra, S1 vertebra, nerve roots, dura mater, ilium);

[0230] B. puncture depth;

[0231] C. Projected area of ​​Kambin triangle;

[0232] Step 5.6: The process of determining the optimal path based on the soft constraint values ​​of the candidate paths that meet the hard constraint conditions is as follows;

[0233] The optimal path decision in step 5.6 is specifically as follows:

[0234] Step 5.6.1. Calculate the soft constraint values ​​of each candidate path and normalize them. Items a, c, and d are calculated using the following formula:

[0235]

[0236] Item b is as follows:

[0237]

[0238] Step 5.6.2: Calculate the risk value of each candidate path. The formula is as follows:

[0239]

[0240] ∑α i +β+r+δ=1;

[0241] Among them, R P is the total risk value of the puncture path, and its value range is [0,1]. The closer the value is to 0, the smaller the risk of the path is, and the closer it is to 1, the greater the risk of the path is. i is the key structure distance risk sub-objective function, a total of 5 key structures are included in the calculation, α i is the weight coefficient corresponding to the key structure distance. S P L is the Kambin triangle projection area risk sub-objective function, and β is the corresponding weight coefficient. Pis the puncture depth risk sub-objective function, and γ is the corresponding weight coefficient. p is the angle between the path and the surface of the intervertebral disc to be ablated, and δ is the corresponding weight coefficient.

[0242] The above step 5.6.2 is performed as follows:

[0243] Step 5.6.2.1. For the eight evaluation indicators, denote them as A_1, A_2, ..., A_n;

[0244] Step 5.6.2.2: Construct a priority graph. For each pair of indicators A_i and A_j (i≠j, 1≤i, j≤n):

[0245] If A_i is more important than A_j, mark the arrow from A_i to A_j in the priority diagram as a_{ij}=1;

[0246] If A_j is more important than A_i, mark the arrow from A_j to A_i in the priority diagram, and denote it as a_{ij}=0;

[0247] If A_i and A_j are equally important, no arrow is drawn. However, to avoid confusion, it is usually considered that a_{ij} = 0.5, and a_{ji} = 0.5 must be ensured.

[0248] Note: a_{ii} is always undefined or set to 0 because it is not compared to itself;

[0249] Step 5.6.2.3. Calculate the importance score. For each indicator A_i, calculate its importance score S_i:

[0250] S_i=∑(a_{ij}-a_{ji}), where j=1 to n and j≠i,

[0251] This formula calculates the "net importance" of indicator A_i relative to all other indicators (taking into account the comparison in both positive and negative directions);

[0252] Step 5.6.2.4, normalize the weights. Calculate the sum of all indicator importance scores:

[0253] S_total=∑S_i

[0254] For each indicator A_i, calculate its normalized weight W_i:

[0255] W_i=S_i / S_total.

[0256] Step 6: Output the path planning environment and the optimal path to a display device, displaying auxiliary decision-making information such as the cross-sectional angulation, coronal angulation, sagittal angulation, target coordinates, skin insertion point coordinates, puncture depth, and Kambin triangle area of ​​the optimal path.

[0257] Example 2, an MR image-based L5 / S1 minimally invasive intervertebral disc radiofrequency ablation surgical path planning device, includes a data preprocessing module, a path planning environment modeling module, a target calibration module, an optimal path decision module, and a visualization module.

[0258] In the data preprocessing module, the MRI format is converted from DICOM to NIFTI format; and the MR image is resampled to an isotropic resolution, specifically 0.75*0.75*0.75.

[0259] The grayscale values ​​of MR images are normalized to the range of [0, 1].

[0260] Finally, the processed MRI data in NIFTI format are input into the path planning environment modeling module.

[0261] In the path planning environment modeling module, the preprocessed MR images are input into the WT-UNet model to segment all target structures (L5 vertebra, S1 vertebra, intervertebral disc, nerve root, dura mater, skin, ilium, and intervertebral disc herniation).

[0262] After the WT-UNet model segments the required targets in the MR image, it also applies morphological operations (such as hole filling and removal of isolated false positive areas) to the output to optimize the segmentation results.

[0263] The binary images of each segmented structure were reconstructed into a 3D image using the Marching Cubes algorithm. The reconstructed 3D model was smoothed with 4 iterations and an iteration factor of 0.4 to complete the path planning environment modeling.

[0264] The target point calibration module is used to model the obtained path planning environment, tetrahedron the 3D sub-model of the intervertebral disc herniation to be ablated, and select effective ablation targets, that is, to obtain the ablation targets.

[0265] The optimal path decision module performs hard constraint condition detection and soft constraint condition detection based on the ablation target and the surface sampling candidate needle entry points on the peritoneal side of the skin 3D model to obtain the optimal path.

[0266] The visualization module includes a display and a projector, which are used to display auxiliary decision information, including the cross-sectional angle, coronal angle, sagittal angle, target coordinates, skin insertion point coordinates, puncture depth, and Kambin triangle area of ​​the optimal path.

Claims

1. A method for minimally invasive radiofrequency ablation of L5 / S1 intervertebral disc based on MR images, characterized in that: The following steps are involved: Step 1: Obtain MRI data of the subject in DICOM format; Step 2: preprocess the data; Step 3: Model the path planning environment and input the pre-processed data in step 2 into the model to obtain the modeled path planning environment; Step 4: Based on the path planning environment obtained in step 3, the 3D sub-model of the intervertebral disc herniation to be ablated is tetrahedroned to select effective ablation targets; Step 5: Construct a candidate path based on the ablation target in Step 4 and the candidate insertion points obtained by sampling the peritoneal surface of the 3D skin model. The candidate path is abstracted as the line connecting the ablation target and each candidate insertion point. Hard and soft constraint conditions are tested on the candidate path to obtain the optimal path. The specific process of obtaining candidate needle insertion points is as follows: Step 5.2.

1. Crop the 3D skin model to retain only the local area on the peritoneal side. Step 5.2.2: All vertices of the triangular mesh of the cropped peritoneal surface of the 3D skin model are added to the candidate needle insertion point set; Step 5.2.3, calculate the area and total area of ​​the triangle; for each triangle, assuming its vertices are A, B, C, then the area S i , Total area S sum , Among them, AB=AB, AC=AC; Step 5.2.

4. Generate a random number r, satisfying 0<=r<1; Step 5.2.

5. Select a triangular face; traverse all faces in order and calculate their cumulative area ratio P i , Find the one that satisfies P i-1 <=r <P i , then the i-th face is selected; Step 5.2.

6. Generate random numbers r1 and r2 in the range [0,1] and ensure that their sum is less than 1. P j =P A +r1*Vector AB+r2*Vector AC; Step 5.2.7: Repeat step 5.2.6 N times to generate N sampling points and add them to the candidate needle insertion point set; N=5*S j ; Step 5.2.

8. Repeat steps 5.2.4 to 5.2.7 10,000 times. Step 6: Output the path planning environment and the optimal path, and display auxiliary decision information, which includes the cross-sectional angle, coronal angle, sagittal angle, target coordinates, skin insertion point coordinates, puncture depth, and Kambin triangle area of ​​the optimal path.

2. The method for minimally invasive radiofrequency ablation of L5 / S1 intervertebral disc based on MR images according to claim 1, characterized in that: Preprocessing includes format conversion, resampling, and grayscale value unification, specifically including the following steps: Step 2.1, convert MRI from DICOM format to NIFTI format; Step 2.2: resample the MR image to an isotropic resolution of 0.75*0.75*0.

75. Step 2.3: Normalize the grayscale values ​​of the MR images to the range of [0, 1]. Step 2.4: Obtain the processed MRI data in NIFTI format.

3. The method for minimally invasive radiofrequency ablation of L5 / S1 intervertebral disc based on MR images according to claim 1, characterized in that: The specific process of step 3 is as follows: Step 3.

1. Input the preprocessed MR image into the WT-UNet model to segment the target structures; the target structures include L5 vertebra, S1 vertebra, intervertebral disc, nerve root, dura mater, skin, ilium, and intervertebral disc herniation; Step 3.2: Use the final model trained in step 3.1 to segment all objects in the MR image, and apply morphological operations to the model output to optimize the segmentation results. Morphological operations include hole filling and removal of isolated false positive regions. Step 3.3: Reconstruct the binary image of each segmented structure into a 3D image using the Marching Cubes algorithm. Smooth the reconstructed 3D model with 4 iterations and an iteration factor of 0.4 to complete the path planning environment modeling. The content of the model training phase in step 3.1 is specifically as follows: Step 3.1.

1. Construct a 3D WT-UNet, using the 3D U-Net framework as the basic framework and integrating the wavelet transform module and the inverse transform module. The encoding and decoding parts of the model can be divided into five stages of feature layers. Each stage of the encoder / decoder contains two convolutional structures. The pooling and upsampling operations of the model are achieved by setting the convolution step size. Step 3.1.2: Replace all up / down sampling parts in the basic model in 3.1.1 with 3D inverse wavelet transform / wavelet transform to realize up / down sampling function and achieve lossless transmission of information. Step 3.1.3: The model uses the cross entropy loss function and the soft threshold dice loss function as the total loss function for training the model, uses the Stochastic Gradient Descent optimizer for model training, and enables the deep supervision strategy during the model training phase; Step 3.1.4: The initial learning rate of the model training is set to 1e-2 and decreases with the increase of the number of training rounds. The total number of training rounds is 1000. The training process uses five-fold cross-validation to select the optimal model parameters.

4. The method for minimally invasive radiofrequency ablation of L5 / S1 intervertebral disc based on MR images according to claim 1, characterized in that: Step 4, "Selecting effective ablation targets," includes the following steps: Step 4.1, tetrahedroning the 3D sub-model of the intervertebral disc herniation to be ablated; Step 4.2: Select the ablation target by scaling, rotating, and moving the 3D sub-model to be ablated and using the cursor-assisted positioning function; Step 4.3: After the target point is initially selected, a "preset ablation area estimation model" with adjustable transparency is generated at the target point; Step 4.4: The doctor evaluates the legitimacy of the target selection based on the visualized effective ablation area in the "preset ablation area estimation model" with adjustable transparency in 4.3, to confirm that the currently selected point is designated as the ablation target.

5. The method for planning a minimally invasive radiofrequency ablation surgical path for L5 / S1 intervertebral disc based on MR images according to claim 1, characterized in that: The step 5 is performed by the following steps: Step 5.1, input the ablation target obtained in step 4; Step 5.2: Sampling the surface of the peritoneal side of the 3D skin model to obtain candidate needle insertion points; Step 5.3: The candidate path is abstracted into the lines between the target and each candidate insertion point; Step 5.4: Check the hard constraints on the candidate paths and delete those that do not meet the hard constraints. Step 5.5: Calculate the soft constraint values ​​of the candidate paths that meet the hard constraint conditions. Step 5.6: Determine the optimal path based on the soft constraint values ​​of the candidate paths that meet the hard constraints.

6. The method for planning a minimally invasive radiofrequency ablation surgical path for L5 / S1 intervertebral disc based on MR images according to claim 3, characterized in that: The specific content of the "3D inverse wavelet transform / wavelet transform" constructed in step 3.1.2 is: Step 3.1.2.1, perform one-dimensional wavelet transform on the image in the X, Y, and Z axis directions to reduce the calculation time; Step 3.1.2.2: Assume that the number of feature channels of the model is C. For the multiplied feature channel number 8C, use a convolution module to reduce the channel number to C. Step 3.1.2.3: For the inverse wavelet transform, the number of model feature channels C needs to be increased to 8C through a convolution module, and then the inverse wavelet transform operation is performed; in the inverse wavelet transform, the image is subjected to a one-dimensional inverse wavelet transform in the XYZ axis direction to reduce the calculation time; the size of the downsampled image after the inverse wavelet transform will be restored to the size before the wavelet transform, and its number of channels will be restored from 8C to C.

7. The method for planning a minimally invasive radiofrequency ablation surgical path for L5 / S1 intervertebral disc based on MR images according to claim 4, characterized in that: In step 4.1, the tetrahedron decomposition is specifically as follows: Step 4.1.1, obtain all vertices of the 3D model to form an initial point set; Step 4.1.2: Construct a super tetrahedron large enough to contain all the initial points. Step 4.1.

3. Randomly insert the initial point into the super tetrahedron one by one. For each newly inserted point, traverse to find all tetrahedrons whose known circumscribed spheres contain the point. Delete the common faces of the influencing tetrahedrons and connect the new point to all vertices of the influencing tetrahedron to form a new tetrahedron. Perform an empty circumscribed sphere check on the newly generated tetrahedron to ensure that it meets the empty circumscribed sphere property of the Delaunay tetrahedron partitioning. If a tetrahedron that does not meet the empty circumscribed sphere property is detected, perform local optimization. Local optimization involves swapping the diagonals of tetrahedra to form new tetrahedra until all tetrahedra satisfy the Delaunay condition; Step 4.1.4: Randomly divide the points on the edges, constraint lines, and constraint surfaces of the initial tetrahedron obtained in step 4.1.

3. The Euclidean distance between the sampled points and the known points cannot be less than 1, otherwise they will be discarded. Perform point interpolation according to step 4.1.

3. Step 4.1.5: Randomly sample points in the tetrahedron space obtained in step 4.1.

4. The Euclidean distance between the sampled points and the known points cannot be less than 1, otherwise they will be discarded. Then perform point interpolation according to step 4.1.

3. Step 4.1.

6. Delete all tetrahedra associated with the supertetrahedron vertices. The remaining tetrahedra constitute the final Delaunay tetrahedron. The specific steps for constructing the super tetrahedron in step 4.1.2 are: Step 4.1.2.

1. Calculate the bounding box of the initial point set; the bounding box is a cuboid whose sides are parallel to the x, y, and z axes, and the cuboid contains all given spatial points, that is, the coordinates of each spatial point are located inside or on the surface of the cuboid; Step 4.1.2.

2. Construct a sphere with the center of the bounding box as the center and a diameter R equal to four times the length of the diagonal of the bounding box. Step 4.1.2.3, construct a regular tetrahedron; the side length of the regular tetrahedron is And the ball in step 4.1.2.2 is the inscribed ball of the regular tetrahedron.

8. The method for planning a minimally invasive radiofrequency ablation surgical path for L5 / S1 intervertebral disc based on MR images according to claim 4, characterized in that: The "cursor-assisted positioning" described in step 4.2 is specifically as follows: Step 4.2.1: Model transformation: transform the vertices of the 3D sub-model to be ablated from model space to world space to obtain the world coordinates of the vertex. worldPosition=modelMatrix·localPosition; Where worldPosition is the coordinate of the model vertex in world space; modelMatrix is ​​the model matrix containing scaling, rotation and translation information; localPosition is the coordinate of the model vertex in model space; Step 4.2.2: View transformation, converting the 3D model from world space to observation space; eyePosition=viewMatrix·worldPosition; Where eyePosition is the coordinate of the model vertex in the observation space; viewMatrix is ​​the view matrix containing the camera position, orientation and "up" direction information; Step 4.2.3: Projection transformation, converting the 3D model from the observation space to the clipping space, and finally mapping it to the normalized device coordinate space; clipPosition=projectionMatrix·eyePosition; Among them, projectionMatrix is ​​the projection matrix, which is perspective projection; Step 4.2.4: Clipping removes vertices outside the viewing frustum, which are not visible in the final image. Clipping is done by comparing with the six clipping planes of the clipping space and the Sutherland-Hodgman vertex clipping algorithm. Step 4.2.5, perspective division, convert the clip space coordinates to normalized device coordinates; Where clipPosition is the coordinate of the model vertex in the clip space; ndcPosition is the coordinate of the model vertex in the normalized device coordinate space; Step 4.2.6: Rasterization converts geometric primitives into pixels on the screen. Vertices are converted into fragments that cover the area of ​​screen pixels, and corresponding attributes are generated for each fragment. Step 4.2.7, depth test, used to perform pixel discarding; compare the depth value of the fragment with the depth value of the corresponding position in the frame buffer; if the depth value of the fragment is less than or equal to the value in the frame buffer, discard the judgment; Step 4.2.8: Viewport transformation, mapping normalized device coordinates to screen coordinates; Where ScreenX is the horizontal coordinate of the screen corresponding to the point, and ScreenY is the vertical coordinate of the screen corresponding to the point; where width and height are the width and height of the viewport, and viewportX and viewportY are the coordinates of the lower left corner of the viewport on the screen; Step 4.2.9, get the screen coordinates of the mouse cursor; Step 4.2.

10. Calculate the distance between the mouse cursor screen coordinates and the 3D model vertex screen coordinates. If there is a vertex whose Euclidean distance to the cursor screen coordinates is less than 2, move the cursor to that point. If the number of vertices whose Euclidean distance is less than 2 is greater than 1, move the cursor to the nearest point.

9. The method for planning a surgical path for minimally invasive L5 / S1 intervertebral disc radiofrequency ablation based on MR images according to claim 5, characterized in that: Hard constraints include the following: A. Avoid collision with risky tissue structures in the L5-S1 segment, and the distance between the puncture path and each structure must be greater than 0.5 mm to leave a safety margin; B. The puncture distance is less than 150 mm of the ablation needle length; C. The angle between the puncture path and the plane of the triangular facet of the 3D model of the intervertebral disc herniation must be greater than 20 degrees to prevent the ablation needle from slipping when entering the intervertebral disc; D. The puncture path must pass through the intervertebral foramen; The calculation of soft constraint condition value is as follows: A. The distance between the puncture path and key tissue structures; B. Puncture depth; C. Kambin triangle projection area; D. The angle between the puncture path and the plane where the triangular facet of the intervertebral disc 3D model intersects.

10. The method for minimally invasive radiofrequency ablation of L5 / S1 intervertebral disc based on MR images according to claim 5, characterized in that: The process of determining the optimal path in step 5.6 is as follows: Step 5.6.

1. Calculate the soft constraint values ​​of each candidate path and normalize them. Items a, c, and d are calculated using the following formula: Item b is as follows: Step 5.6.2: Calculate the risk value of each candidate path. The formula is as follows: Yes i +β+r+δ=1; Among them, R P is the total risk value of the puncture path, and its value range is [0,1]. The closer the value is to 0, the smaller the risk of the path is, and the closer it is to 1, the greater the risk of the path is. i is the key structure distance risk sub-objective function, a total of 5 key structures are included in the calculation, α i is the weight coefficient corresponding to the key structure distance; S P is the risk sub-objective function of Kambin triangle projection area, β is the corresponding weight coefficient; L P is the puncture depth risk sub-objective function, γ is the corresponding weight coefficient; A p is the angle between the path and the surface of the intervertebral disc to be ablated, and δ is the corresponding weight coefficient; The above step 5.6.2 is performed as follows: Step 5.6.2.

1. For n evaluation indicators, denote them as A_1, A_2, ..., A_n; Step 5.6.2.

2. Construct a precedence diagram. For each pair of indicators A_i and A_j (i≠j, 1≤i, j≤n): if A_i is more important than A_j, mark the arrow from A_i to A_j in the precedence diagram as a_{ij}=1. If A_j is more important than A_i, mark the arrow from A_j to A_i in the priority diagram, and denote it as a_{ij}=0; If A_i and A_j are equally important, no arrow is drawn. However, to avoid confusion, it is usually considered that a_{ij} = 0.5, and a_{ji} = 0.5 must be ensured. a_{ii} is always undefined or set to 0; Step 5.6.2.

3. Calculate the importance score. For each indicator A_i, calculate its importance score S_i: S_i=∑(a_{ij}-a_{ji}), where j=1 to n and j≠i; Step 5.6.2.4, normalize the weights; calculate the sum of all indicator importance scores: S_total=∑S_i, For each indicator A_i, calculate its normalized weight W_i: W_i=S_i / S_total.

Images