L5 / S1 minimally invasive intervertebral disc radiofrequency ablation operation path planning method and device based on MR image

Abstract

The invention provides an L5 / S1 minimally invasive intervertebral disc radiofrequency ablation operation path planning method and device based on an MR image. The method comprises the following steps that MRI data of an object in a DICOM format are obtained; preprocessing the data; modeling a path planning environment, and inputting data into a modeled model to obtain a modeled path planning environment; according to the obtained path planning environment, performing tetrahedron subdivision on the to-be-ablated 3D sub-model of the protrusion of intervertebral disc, and selecting an effective ablation target spot; constructing a candidate path; performing hard constraint condition detection and soft constraint condition detection on the candidate path to obtain an optimal path; and outputting the path planning environment and the optimal path. And the determined operation path can reduce the operation risk and the occurrence rate of complications while ensuring the operation effect. Clinical information such as cross section angulation, coronal plane angulation and sagittal plane angulation is additionally provided, and working efficiency and decision quality are improved.

Description

Technical Field

[0001] The present invention relates to the technical field of minimally invasive radiofrequency ablation of intervertebral disc surgery, and in particular to a method and device for planning a path for L5 / S1 minimally invasive radiofrequency ablation of intervertebral disc surgery based on MR images. Background Art

[0002] In the field of minimally invasive medicine, the treatment technology for lumbar disc herniation continues to evolve. Among them, radiofrequency ablation has become an important treatment method for relieving pain caused by lumbar disc herniation due to its low trauma, fast recovery and high safety. This technology acts on the protruding disc tissue by precisely controlling low or high temperature energy to achieve local decompression, thereby effectively relieving low back pain caused by nerve compression. The key to the success of the operation lies in accurately identifying and locating the ablation target, as well as safely and efficiently planning the surgical path and guiding the puncture needle to the target area.

[0003] However, the current radiofrequency ablation of lumbar disc herniation still faces many challenges in its implementation. First, the planning of the surgical path is highly dependent on the doctor's personal experience and subjective judgment, which not only prolongs the operation time, but may also increase the surgical risk and the incidence of complications due to lack of experience. Secondly, although the traditional C-arm or CT-guided puncture operation can assist in positioning to a certain extent, it is often difficult to accurately locate the compressed nerve root due to the limitations of imaging accuracy and the complexity of the anatomical structure, increasing the risk of accidental injury to the nerve during surgery. In addition, individual differences among patients, variations in the anatomical position of the nerve roots, and uncertainty in the location of nerve compression have further increased the difficulty and uncertainty of the operation.

[0004] To this end, Chinese patent CN118252614A proposes a puncture path planning method for radiofrequency ablation of lumbar disc herniation via the intervertebral foramen, which uses the MRI data of the object and goes through 5 steps to finally obtain the optimal path. The invention uses deep learning automatic segmentation to achieve the same accuracy as manual segmentation and greatly speeds up the segmentation speed. The algorithm's running time is greatly accelerated through computer parallel computing and headless rendering. However, the above-mentioned path planning method for radiofrequency ablation of lumbar disc herniation still has the following defects in practical use:

[0005] First, in the existing path planning method for radiofrequency ablation of lumbar disc herniation, doctors need to operate on third-party commercial software (such as Mimics) when manually planning the path. Frequent switching between different software and differences in standards affect the doctors' work efficiency and operation accuracy. In addition, the segmentation method used in the comparison document is the general deep learning segmentation framework 3D U-Net, but its performance still has room for improvement in the fine segmentation of the specific anatomical structure L5 / S1 lumbar spine.

[0006] Second, the segmentation method used in the comparison document can only show the path planning results, but less 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 prior art and provide a method and device for L5 / S1 minimally invasive intervertebral disc radiofrequency ablation surgery path planning based on MR images. The 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 object, the present invention is implemented through the following technical measures:

[0010] A L5 / S1 minimally invasive intervertebral disc radiofrequency ablation surgery path planning device 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 L5 / S1 intervertebral disc based on MR images is performed by the following steps:

[0012] Step 1, obtaining MRI data of the object in DICOM format;

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

[0014] Step 3: Input the preprocessed 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 the optimal path decision module fails to find any feasible path that meets the preset safety standards after comprehensive search and analysis, it will go to the target calibration module and suggest recalibrating the target position.

[0017] Step 6: Output the path planning environment and the optimal path to the visualization module to display the cross-sectional angle, coronal angle, sagittal angle, target point 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, specifically 0.75*0.75*0.75;

[0021] Step 2.3, normalize the grayscale value of the MR image 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 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.

[0026] Step 3.3: Use the Marching Cubes algorithm to reconstruct the binary image of each structure into a 3D image. The reconstructed 3D model is smoothed, the number of iterations is set to 4, and the iteration factor is set to 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, build 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 convolution structures, and the pooling and upsampling operations of the model are realized by setting the convolution step size.

[0029] Step 3.1.2, replace all up / down sampling parts in the basic model of 3.1.1 with three-dimensional inverse wavelet transform / wavelet transform to realize the 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 to train the model, uses the Stochastic Gradient Descent optimizer for model training, and enables the deep supervision strategy during the training phase of the model.

[0031] 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, and the total number of training rounds is 1000. The training process uses five-fold cross validation to select the optimal parameters of the model;

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

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

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

[0035] 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. Similarly, in the inverse wavelet transform, the image is subjected to a one-dimensional inverse wavelet transform in the XYZ axis direction to reduce the time required for calculation. 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.

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

[0037] Step 4.1, tetrahedron-dividing 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 semi-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 current 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 balls contain the point (called influential 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 ball detection on the newly generated tetrahedrons to ensure that they satisfy the empty circumscribed ball characteristics of the Delaunay tetrahedron partitioning. If a tetrahedron that does not satisfy the empty circumscribed ball characteristics is detected, local optimization is required. This usually involves swapping the diagonals of the tetrahedrons to form new tetrahedrons until all tetrahedrons satisfy the Delaunay conditions.

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

[0046] Step 4.1.5: Randomly sample the points in the tetrahedral 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 insertion according to step 4.1.3.

[0047] Step 4.1.6. Delete all tetrahedrons associated with the super tetrahedron vertices. 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 the cuboid can just 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 the diameter R as 4 times the length of the diagonal of the bounding box.

[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, convert the vertices of the 3D sub-model to be ablated from the model space to the world space, and obtain the world coordinates worldPosition corresponding to 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 scaling, 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 observation space (camera space). viewMatrix is ​​the view matrix (4×4 matrix) that contains information about the camera's position, orientation, and "up" direction.

[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, remove the vertices outside the view frustum, which are not visible in the final image. The clipping operation is performed by comparing with the six clipping planes (left, right, top, bottom, near, far) of the clipping space and 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 the clip space. ndcPosition is the coordinate (vector) of the model vertex in the 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 covering 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 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, the judgment is discarded.

[0068] Step 4.2.8: Viewport transformation, mapping normalized device coordinates (NDC) to screen coordinates.

[0069]

[0070] Where ScreenX is the horizontal screen coordinate corresponding to the point, and ScreenY is the vertical screen coordinate 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 (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 screen coordinate whose Euclidean distance to the cursor screen coordinates is less than 2, move the cursor to that point; if the number of vertices with Euclidean distance 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] Step 5.1: The ablation target points obtained in step 4 are input into the optimal path decision module.

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

[0076] Step 5.3: The candidate paths are abstracted as lines between the ablation target and each candidate needle entry point.

[0077] Step 5.4: Check the hard constraint conditions of the candidate paths and delete the paths that do not meet the conditions.

[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 constraint conditions.

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

[0081] Step 5.2.1. Crop the skin 3D model and keep only the local area on the peritoneum side.

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

[0083] 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 ,

[0084]

[0085] Total area sum ,

[0086]

[0087] Among them, 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, 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 set of needle insertion points.

[0095] N=5*S j ;

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

[0097] In step 5.4, the candidate path hard constraint condition detection 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] In step 5.6, the optimal path decision is specifically:

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

[0103]

[0104] Item b is calculated according to the following formula:

[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. P is the risk sub-objective function of the Kambin triangle projection area, 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, denoted 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, a_{ij}=0.5 can be considered, 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 A_i relative to all other indicators (taking into account comparisons in both positive and negative directions).

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

[0121] S_total=∑S_i

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

[0123] W_i=S_i / S_total.

[0124] Advantages compared to the prior art:

[0125] First, the improved WT-UNet deep neural segmentation network not only inherits the powerful feature extraction and context information integration capabilities of 3D U-Net, but also creatively incorporates the wavelet transform module to alleviate the problem of information loss in up / down sampling operations. Compared with 3D U-Net, it has improved the DICE coefficient, precision, recall rate and other indicators in the segmentation of bony tissue and non-bone tissue in the L5 / S1 segment, 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 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 takes into account the Kambin triangle area of ​​the L5 / S1 segment. Based on comprehensive consideration of various factors, a priority diagram is constructed to convert the physician's subjective experience into the weight of the indicator, and the physician is allowed to fine-tune the weight according to the actual situation, so that the determined surgical path can minimize the surgical risk and the incidence of complications while ensuring the surgical effect.

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

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

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

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

[0132] Figure 4 A schematic diagram of the path planning environment of the present invention;

[0133] Figure 5 It is a schematic structural diagram of a 3D sub-model of a herniated intervertebral disc to be ablated which is divided into tetrahedrons according to the present invention;

[0134] Figure 6 : 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 is an example diagram of the optimal path of the present invention;

[0137] Fig. 9 This is a comparison diagram of DSC results segmented by different models. DETAILED DESCRIPTION

[0138] Example 1, refer to the attached Figure 1-8 , a L5 / S1 minimally invasive radiofrequency ablation surgical path planning method based on MR images, comprising the following steps:

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

[0140] Step 2: pre-process the data;

[0141] Preprocessing includes format conversion, resampling, and grayscale value unification, which specifically includes 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, specifically 0.75*0.75*0.75;

[0144] Step 2.3, normalize the grayscale value of the MR image to the range of [0,1];

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

[0146] Step 3: Model the path planning environment, and input the preprocessed 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, build 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 convolution structures, and the pooling and upsampling operations of the model are realized by setting the convolution step size.

[0150] Step 3.1.2, replace all up / down sampling parts in the basic model of 3.1.1 with three-dimensional inverse wavelet transform / wavelet transform to realize the 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 to train the model, uses the Stochastic Gradient Descent optimizer for model training, and enables the deep supervision strategy during the training phase of the model.

[0152] 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, and the total number of training rounds is 1000. The training process uses five-fold cross validation to select the optimal parameters of the model;

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

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

[0155] Step 3.1.2.2, assuming that the number of feature channels of the model is C, for the multiplied feature channel number 8C, a convolution module is used to reduce the number of channels to C. Due to the characteristics of wavelet transform, the image size of the original image is reduced to half of the original image after a wavelet transform, thus having the same effect as the traditional downsampling operation in terms of image size.

[0156] 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. Similarly, in the inverse wavelet transform, the image is subjected to a one-dimensional inverse wavelet transform in the XYZ axis direction to reduce the time required for calculation. 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.

[0157] Step 3.2: Use the final model trained in 3.1 to segment all targets in the MR image, and use 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 image of each structure into a 3D image. The reconstructed 3D model is smoothed, the number of iterations is set to 4, and the iteration factor is set to 0.4 to complete the path planning environment modeling.

[0159] Step 4: According to the path planning environment obtained in step 3, the 3D sub-model of the intervertebral disc herniation to be ablated is tetrahedron-divided, and an effective ablation target is selected; that is, the ablation target is obtained.

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

[0161] Step 4.1, tetrahedron-dividing 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 the cuboid can just 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 the diameter R as 4 times the length of the diagonal of the bounding box.

[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 balls contain the point (called influential 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 ball detection on the newly generated tetrahedrons to ensure that they satisfy the empty circumscribed ball characteristics of the Delaunay tetrahedron partitioning. If a tetrahedron that does not satisfy the empty circumscribed ball characteristics is detected, local optimization is required. This usually involves swapping the diagonals of the tetrahedrons to form new tetrahedrons until all tetrahedrons satisfy the Delaunay conditions.

[0170] Step 4.1.4: Randomly divide the edges, constraint lines, and constraint faces of the initial tetrahedral decomposition obtained in step 4.1.3. The Euclidean distance between the sampling points and the known points cannot be less than 1, otherwise they will be discarded. Insert the points according to step 4.1.3.

[0171] Step 4.1.5: Randomly sample the points in the tetrahedral 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 insertion 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, convert the vertices of the 3D sub-model to be ablated from the model space to the world space, and obtain the world coordinates worldPosition corresponding to 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 scaling, 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 observation space (camera space). viewMatrix is ​​the view matrix (4×4 matrix) that contains information about the camera's position, orientation, and "up" direction.

[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, remove the vertices outside the view frustum, which are not visible in the final image. The clipping operation is performed by comparing with the six clipping planes (left, right, top, bottom, near, far) of the clipping space and 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 the clip space. ndcPosition is the coordinate (vector) of the model vertex in the 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 covering 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 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, the judgment is discarded.

[0190] Step 4.2.8: Viewport transformation, mapping normalized device coordinates (NDC) to screen coordinates.

[0191]

[0192] Where ScreenX is the horizontal screen coordinate corresponding to the point, and ScreenY is the vertical screen coordinate 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 (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 screen coordinate whose Euclidean distance to the cursor screen coordinates is less than 2, move the cursor to that point; if the number of vertices with Euclidean distance less than 2 is greater than 1, move the cursor to the nearest point.

[0195] Step 4.3: After the doctor preliminarily selects the target point, an estimated range of the ablation area with adjustable transparency is generated at the target point (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 condition detection and soft constraint condition detection on the ablation target in step 4 and the candidate needle entry point obtained by surface sampling on the peritoneal side of the skin 3D model to obtain the optimal path; that is, construct a candidate path based on "ablation target in step 4" and "candidate needle entry point obtained by surface sampling on the peritoneal side of the skin 3D model", and the candidate path is abstracted as a connection between the ablation target and each candidate needle entry point; perform hard constraint condition detection and soft constraint condition detection 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 point obtained in step 4;

[0200] Step 5.2, sampling the surface of the peritoneum side of the skin 3D 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 skin 3D model and keep only the local area on the peritoneum side.

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

[0204] 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 ,

[0205]

[0206] Total area sum ,

[0207]

[0208] Among them, 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, 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 set of needle insertion points.

[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 paths are abstracted into lines between the target and each candidate injection point.

[0219] Step 5.4: Check the hard constraint conditions of the candidate paths, and delete the paths that do not meet the conditions;

[0220] Hard constraints include the following:

[0221] A. Avoid risky tissue structures (L5 vertebra, S1 vertebra, nerve root, dura mater, ilium) in the L5-S1 segment. No collision should occur, and the distance should be greater than 0.5 mm, leaving a safety margin.

[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 facet of the intervertebral disc 3D model intersects 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, calculating the soft constraint condition 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. The distance between the puncture path and key tissue structures (L5 vertebra, S1 vertebra, nerve root, dura mater, ilium);

[0230] B. puncture depth;

[0231] C. Kambin triangle projection area;

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

[0233] In step 5.6, the optimal path decision is specifically:

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

[0235]

[0236] Item b is calculated according to the following formula:

[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. P is the risk sub-objective function of the Kambin triangle projection area, 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 by the following steps:

[0243] Step 5.6.2.1, denoted as A_1, A_2, ..., A_n for the 8 evaluation indicators;

[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, then mark the arrow from A_i to A_j in the priority diagram, denoted 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, denoted as a_{ij}=0;

[0247] If A_i and A_j are equally important, no arrow is drawn, but usually to avoid confusion, a_{ij} = 0.5 can be considered, 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 weights. Calculate the sum of all indicator importance scores:

[0253] S_total=∑S_i

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

[0255] W_i=S_i / S_total.

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

[0257] Embodiment 2, a L5 / S1 minimally invasive intervertebral disc radiofrequency ablation surgery path planning device 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.

[0258] In the data preprocessing module, the MRI is converted from DICOM format 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 images, 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 image of each structure after segmentation is reconstructed into a 3D image using the Marching Cubes algorithm. The reconstructed 3D model is smoothed, the number of iterations is set to 4, and the iteration factor is set to 0.4 to complete the path planning environment modeling.

[0264] The target point calibration module is used to model the path planning environment according to 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, and the auxiliary decision information includes the cross-sectional angle, coronal angle, sagittal angle, target point 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, obtaining MRI data of the object in DICOM format; Step 2: pre-process the data; Step 3: Model the path planning environment, and input the preprocessed data in step 2 into the model to obtain the modeled path planning environment; Step 4: According to the path planning environment obtained in step 3, the 3D sub-model of the intervertebral disc herniation to be ablated is tetrahedron-divided to select effective ablation targets; Step 5: Construct a candidate path based on the "ablation target in step 4" and "candidate needle entry points obtained by surface sampling on the peritoneal side of the skin 3D model". The candidate path is abstracted as a line between the ablation target and each candidate needle entry point; Perform hard constraint condition detection and soft constraint condition detection on candidate paths to obtain the optimal path; 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 point 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, which specifically includes 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, specifically 0.75*0.75*0.75; Step 2.3, normalize the grayscale value of the MR image to the range of [0,1]; Step 2.4: Obtain 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, the preprocessed MR image is input into the WT-UNet model to segment the target structure; the target structure includes 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 targets in the MR image, and use morphological operations on the output of the model to optimize the segmentation results; the morphological operations include hole filling and removal of isolated false positive areas; Step 3.3, reconstruct the binary image of each structure after segmentation into a 3D image using the Marching Cubes algorithm; smooth the reconstructed 3D model, set the number of iterations to 4, and the iteration factor to 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 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, and each stage of the encoder / decoder contains two convolution structures, and the pooling and upsampling operations of the model are realized by setting the convolution step size; Step 3.1.2, replace all up / down sampling parts in the basic model of 3.1.1 with three-dimensional inverse wavelet transform / wavelet transform to realize the 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 to train the model, uses the Stochastic Gradient Descent optimizer for model training, and enables the deep supervision strategy during the training phase of the model. Step 3.1.4: The initial learning rate of the model training is set to 1e-2 and decreases with the increase in the number of training times. 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, tetrahedron-dividing the 3D sub-model of the intervertebral disc herniation to be ablated; Step 4.2, by scaling, rotating, and moving the 3D sub-model to be ablated, and using the cursor-assisted positioning function, select the ablation target; 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, so as to confirm that the current selected point is set as the ablation target.

5. The method for minimally invasive radiofrequency ablation of 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 point obtained in step 4; Step 5.2, sampling the surface of the peritoneum side of the skin 3D model to obtain candidate needle insertion points; Step 5.3, the candidate path is abstracted as the connection between the target and each candidate injection point; Step 5.4: Check the hard constraint conditions of the candidate paths, and delete the paths that do not meet the conditions; Step 5.5, calculating the soft constraint condition 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 constraint conditions.

6. The method for minimally invasive radiofrequency ablation of L5 / S1 intervertebral disc based on MR images according to claim 3, characterized in that: The specific contents of the "3D wavelet inverse transform / wavelet transform" constructed in step 3.1.2 are: Step 3.1.2.1, perform one-dimensional wavelet transform on the image in the XYZ axis direction to reduce the calculation time; Step 3.1.2.2, assuming that the number of feature channels of the model is C, for the multiplied feature channel number 8C, a convolution module is used to reduce the number of channels 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 the number of channels will be restored from 8C to C.

7. The method for planning a surgical path for L5 / S1 minimally invasive radiofrequency ablation of 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 regular tetrahedron point by point; for each newly inserted point, traverse to find all tetrahedrons whose known circumscribed balls contain the point; delete the common faces of the affected tetrahedrons, and connect the new point with all the vertices of the affected tetrahedron to form a new tetrahedron; perform empty circumscribed ball detection on the newly generated tetrahedron to satisfy the empty circumscribed ball property of Delaunay tetrahedral partitioning; if a tetrahedron that does not satisfy the empty circumscribed ball property is detected, local optimization is required; local optimization involves exchanging the diagonals of the tetrahedrons to form new tetrahedrons until all tetrahedrons satisfy the Delaunay condition; Step 4.1.4: Randomly divide the points in the edges, constraint lines, and constraint surfaces of the initial tetrahedral decomposition obtained in step 4.1.

3. The Euclidean distance between the sampling points and the known points cannot be less than 1, otherwise they will be discarded; perform point insertion according to step 4.1.3; Step 4.1.5: Randomly sample the points in the tetrahedral 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; and insert the points according to step 4.1.3; Step 4.1.6, delete all tetrahedrons associated with the super tetrahedron vertices, and the remaining tetrahedrons constitute the final Delaunay tetrahedron; The specific steps of 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 respectively, 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; take the center of the bounding box as the center of the sphere, and the diameter R is 4 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 surgical path for minimally invasive L5 / S1 intervertebral disc radiofrequency ablation 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, convert the vertices of the 3D sub-model to be ablated from the model space to the world space, and obtain the world coordinates worldPosition corresponding to 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 a 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 covering the screen pixel area, 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 screen coordinate corresponding to the point, and ScreenY is the vertical screen coordinate 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 screen coordinate whose Euclidean distance to the cursor screen coordinates is less than 2, move the cursor to that point; if the number of vertices with Euclidean distance less than 2 is greater than 1, move the cursor to the nearest point.

9. The method for planning a surgical path for L5 / S1 minimally invasive radiofrequency ablation of intervertebral disc based on MR images according to claim 5, characterized in that: The specific process of obtaining candidate needle insertion points in step 5.2 is as follows: Step 5.2.1, crop the skin 3D model and keep only the local area on the peritoneum side; Step 5.2.2, all vertices of the triangular mesh of the cropped peritoneal surface of the skin 3D model are added to the candidate point set of the needle insertion point; 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 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, 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 point set for needle insertion; N=5*S j ; Step 5.2.8, repeat steps 5.2.4 to 5.2.7 for a total of 10,000 times; Hard constraints include the following: A. It is forbidden to collide with the risk tissue structure of the L5-S1 segment, and the distance between the puncture path and each structure must be greater than 0.5mm to leave a safety margin; B. The puncture distance is less than the length of the ablation needle by 150 mm; C. The angle between the puncture path and the plane of the triangular facet where the 3D model of the intervertebral disc herniation intersects 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 plane where the puncture path and the triangular facet of the 3D model of the intervertebral disc intersect.

10. The method for planning a surgical path for L5 / S1 minimally invasive radiofrequency ablation of 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 according to the following formula: Item b is calculated according to the following formula: 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 its 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 the 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 by the following steps: Step 5.6.2.1, for n evaluation indicators, record them as A_1, A_2, ..., A_n; Step 5.6.2.2, construct a priority 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 priority diagram, denoted 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, denoted as a_{ij}=0; If A_i and A_j are equally important, no arrow is drawn, but usually to avoid confusion, a_{ij} = 0.5 can be considered, 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 weights; calculate the sum of all indicator importance scores: S_total=∑S_i, For each indicator A_i, calculate its normalized weight Wi: W_i=S_i / S_total.

Images