Method for automatic evaluation of cervical spinal cord compression degree based on MRI image

Abstract

The present application relates to the technical field of image processing, in particular to a cervical vertebra spinal cord compression degree automatic evaluation method based on MRI image. The present application obtains a three-dimensional sparse graph structure composed of voxels as nodes, obtains a connected probability value of each node according to a one-dimensional initial restart distribution vector composed of a tissue density distribution of all nodes and a starting lumen coordinate and a terminal lumen coordinate; obtains a comprehensive search generation value of adjacent nodes according to a coordinate position distribution of the adjacent nodes and the terminal lumen coordinate, a tissue density, a connected probability value and a pixel calibration interval constant, executes a path search algorithm to obtain a three-dimensional spinal cord center coordinate sequence; and obtains a spinal cord compression degree according to a coordinate distribution in the three-dimensional spinal cord center coordinate sequence and a gray scale distribution of different voxels in a three-dimensional image matrix. The present application obtains accurate three-dimensional spinal cord center coordinates, and more objectively and accurately evaluates the spinal cord compression degree.

Description

Technical Field

[0001] This invention relates to the field of image processing technology, and specifically to an automatic assessment method for the degree of cervical spinal cord compression based on MRI images. Background Technology

[0002] In standard T2-weighted imaging, the cerebrospinal fluid, which acts as an isolation zone, presents a high-brightness water signal, while the spinal cord parenchyma and surrounding bones and intervertebral discs present a medium-to-low gray signal. When severe tissue degeneration occurs, the protruding tissue severely compresses the spinal canal, causing the local high-brightness cerebrospinal fluid to be completely drained. The medium-to-low gray soft tissue of the spinal cord comes into direct contact with the protruding tissue, which also presents a low gray signal, resulting in severe tissue gray-scale aliasing in the local image.

[0003] In existing technologies, graph search algorithms are introduced to find a continuous central line of the spinal canal in three-dimensional space as a reference, and then calculate the cross-sectional area of ​​the spinal canal to assess the degree of compression. However, conventional graph search algorithms rely solely on local grayscale features for path expansion, ignoring the fact that the blind ends of the nerve roots extending to both sides of the cervical spine are also filled with cerebrospinal fluid with high-brightness features. When the search wavefront encounters high resistance caused by the compression of diseased tissue in the main spinal canal, it is very easy to get stuck in the anatomical blind end, resulting in poor accuracy in assessing the degree of compression. Summary of the Invention

[0004] To address the problem that conventional image search algorithms, which rely solely on local grayscale features for path expansion, are prone to getting stuck in anatomical blind spots, leading to poor accuracy in assessing the degree of compression, this invention aims to provide an automatic assessment method for the degree of cervical spinal cord compression based on MRI images. The specific technical solution adopted is as follows: This invention proposes an automatic assessment method for the degree of cervical spinal cord compression based on MRI images, the method comprising: Obtain the three-dimensional image matrix of the cervical spinal cord and the pixel calibration spacing constant in different dimensions; Based on the grayscale characteristics of different voxels in the 3D image matrix, the tissue density of each voxel is obtained; based on the tissue density distribution of voxels in 2D slices at different positions on the Z-axis of the 3D image matrix, the coordinates of the starting lumen and the ending lumen are obtained. A three-dimensional sparse graph structure composed of voxels as nodes is obtained. Based on the tissue density distribution of all nodes, the state transition weights between different nodes are obtained. Based on the one-dimensional initial restart distribution vector composed of the starting lumen coordinates and the ending lumen coordinates, the preset restart probability constant, and the state transition weight distribution between different nodes, the connectivity probability value of each node is obtained. Based on the coordinate distribution of each node's neighboring nodes and the coordinates of the terminating lumen, tissue density, connectivity probability, and pixel calibration spacing constant, the comprehensive search cost of neighboring nodes is obtained. Based on the comprehensive search cost of different neighboring nodes between the starting and ending lumen coordinates, a path search algorithm is executed to obtain a three-dimensional discrete coordinate sequence and a three-dimensional spinal cord center coordinate sequence. Based on the coordinate distribution in the three-dimensional spinal cord center coordinate sequence, the pixel calibration spacing constant, and the gray-level distribution of adjacent voxels at corresponding positions in the three-dimensional image matrix for different center coordinates, the cross-sectional area of ​​the spinal cord canal at each center coordinate is obtained, and the degree of spinal cord compression is obtained.

[0005] Furthermore, the method for obtaining the tissue density includes: Obtain the maximum and minimum gray values ​​of all voxels in the 3D image matrix; The gray values ​​of each voxel are normalized based on the maximum and minimum gray values, and the normalization results are negatively correlated and mapped as the tissue density of each voxel.

[0006] Furthermore, the method for obtaining the starting lumen coordinates and the ending lumen coordinates includes: Based on the tissue density distribution of voxels in two-dimensional slices at different positions on the Z-axis of the three-dimensional image matrix, the main vertebral canal connectivity domain of the two-dimensional slices at the corresponding positions is obtained. For any given location, the mean value of the position coordinates of all voxels within the main spinal canal connected domain on the two-dimensional slice is obtained as the average coordinate; the relative distance between the position coordinates of different voxels within the main spinal canal connected domain and the average coordinate is obtained, and the position coordinate of the voxel corresponding to the minimum relative distance is used as the average corrected coordinate. The three-dimensional coordinates formed by the highest position of the Z-axis and the average corrected coordinates on the corresponding two-dimensional slice are used as the starting lumen coordinates; the three-dimensional coordinates formed by the lowest position of the Z-axis and the average corrected coordinates on the corresponding two-dimensional slice are used as the ending lumen coordinates.

[0007] Furthermore, the method for obtaining the main vertebral canal connectivity region includes: For any location, if the tissue density of a voxel in the two-dimensional slice is less than the preset low-resistivity cutoff constant, the corresponding voxel will be used as a candidate voxel. Obtain the connected domains formed by consecutive adjacent candidate voxels, and select the connected domain with the most candidate voxels among all connected domains as the main vertebral canal connected domain.

[0008] Furthermore, the method for obtaining the state transition weights includes: For any node, obtain the difference between the positive integer 1 and the tissue density of each node's corresponding neighboring nodes, calculate the sum of the difference result and the preset adjustment coefficient, and use it as the state transition weight of each node relative to its neighboring nodes. For nodes other than their neighbors, set the state transition weight of each node relative to the other nodes to 0.

[0009] Furthermore, the method for obtaining the connectivity probability value includes: Construct a one-dimensional initial restart distribution vector for the number of nodes, obtain the one-dimensional linear index positions of the starting and ending lumen coordinates, set the element values ​​of the corresponding positions of the one-dimensional initial restart distribution vector to preset probability values, set the element values ​​of other positions to 0, normalize the state transition weights between different nodes, and construct a global state transition sparse matrix between different nodes. Using the one-dimensional initial restart distribution vector as the initial probability distribution vector, we obtain the product matrix of the global state transition sparse matrix and the initial probability distribution vector. We obtain the difference between the positive integer 1 and the preset restart probability constant, and calculate the product of the difference result and the product matrix as the transition probability vector. We obtain the product of the preset restart probability constant and the one-dimensional initial restart distribution vector as the restart probability vector. We obtain the sum of the transition probability vector and the restart probability vector as the probability distribution vector for the next iteration. Obtain the difference vector between the probability distribution vector of the next iteration and the initial probability distribution vector. Calculate the sum of the absolute values ​​of all elements in the difference vector. If the sum of the absolute values ​​is greater than or equal to a preset convergence threshold, use the probability distribution vector of the next iteration as the new initial probability distribution vector. Obtain the probability distribution vector for the corresponding next iteration. Continue until the sum of the absolute values ​​is less than the preset convergence threshold. Map the probability distribution vector of the next iteration onto a three-dimensional space to obtain the connectivity probability value of each node.

[0010] Furthermore, the method for obtaining the comprehensive search cost value includes: For any two nodes, obtain the product of the difference in coordinates between the nodes in each dimension and the pixel calibration spacing constant, as the actual spacing in each dimension; obtain the sum of the squares of the actual spacings in all dimensions and take the square root, as the actual distance between the corresponding nodes. Based on the organization compactness, connectivity probability value, and actual distance between a node and its neighboring nodes, the single-step movement cost increment between a node and its neighboring nodes is obtained. Organization compactness and actual distance are positively correlated with the single-step movement cost increment, while connectivity probability value is negatively correlated with the single-step movement cost increment. Obtain the actual distance between the corresponding adjacent node and the coordinates of the termination cavity as the heuristic cost; obtain the cumulative movement cost from the starting cavity coordinates to the current node; obtain the sum of the cumulative movement cost, the heuristic cost of adjacent nodes, and the single-step movement cost increment as the comprehensive search cost of adjacent nodes.

[0011] Further, obtaining the three-dimensional discrete coordinate sequence and the three-dimensional spinal cord center coordinate sequence includes: The search origin is set at the starting lumen coordinates, and the ending lumen coordinates are set at the ending lumen coordinates. The search origin is added to the candidate node set. Each time, the node with the smallest comprehensive search cost is selected from the candidate node set for expansion. If the corresponding node is the ending point, the search ends. Otherwise, its adjacent unvisited nodes are added to the candidate node set, and the cumulative movement cost of the adjacent nodes is updated until the search ending point is reached. The search origin is backtracked from the search ending point to obtain a three-dimensional discrete coordinate sequence consisting of all the nodes passed through. A sliding filter is applied to the three-dimensional discrete coordinate sequence, and the filtered three-dimensional discrete coordinate sequence is used as the three-dimensional spinal cord center coordinate sequence.

[0012] Furthermore, the method for obtaining the cross-sectional area of ​​the spinal canal includes: An orthogonal cross-sectional mesh is constructed along the three-dimensional spinal cord center coordinate sequence. The orthogonal cross-sectional mesh is then subjected to grayscale interpolation using a three-dimensional image matrix interpolation to generate a two-dimensional orthogonal slice image. Obtain the grayscale variance of all pixels in each two-dimensional orthogonal slice image. If the grayscale variance is less than the preset variance threshold, set the cross-sectional area of ​​the spinal canal of the corresponding two-dimensional orthogonal slice image to 0. If the gray-level variance is greater than or equal to the preset variance threshold, the Otsu algorithm is used to adaptively segment the two-dimensional orthogonal slice image and perform binarization to obtain the connected component formed by adjacent pixels after binarization. A connected region with a central coordinate is selected as the closed lumen contour; the product of the pixel calibration spacing constants in two dimensions on the horizontal plane is obtained as the pixel calibration physical area; the product of the number of pixels in the closed lumen contour and the pixel calibration physical area is obtained as the cross-sectional area of ​​the spinal canal of the corresponding two-dimensional orthogonal slice image.

[0013] Furthermore, the method for obtaining the degree of spinal cord compression includes: Obtain the spinal cord lumen cross-sectional area sequence corresponding to the three-dimensional spinal cord center coordinate sequence, and obtain the minimum value of the element in the spinal cord lumen cross-sectional area sequence as the lumen cross-sectional area of ​​the most severe compression point; Within the neighborhood of the cross-sectional area of ​​the lumen at the most severe compression point, the maximum value in the sequence of all spinal cord lumen cross-sectional areas is obtained as the local healthy baseline area. The degree of spinal cord compression is obtained based on the cross-sectional area of ​​the lumen at the point of most severe compression and the local healthy baseline area. The cross-sectional area of ​​the lumen at the point of most severe compression is negatively correlated with the degree of spinal cord compression, while the local healthy baseline area is positively correlated with the degree of spinal cord compression.

[0014] The present invention also proposes an automatic assessment system for the degree of cervical spinal cord compression based on MRI images. The system includes a memory, a processor, and a computer program stored in the memory and executable on the processor. When the processor executes the computer program, it implements any of the steps of the automatic assessment method for the degree of cervical spinal cord compression based on MRI images.

[0015] The present invention has the following beneficial effects: This invention obtains a three-dimensional sparse graph structure composed of voxels as nodes. Based on the tissue density distribution of all nodes, it obtains the state transition weights between different nodes, reflecting the probability of transitions between nodes. Based on the one-dimensional initial restart distribution vector composed of the starting and ending lumen coordinates, the preset restart probability constant, and the state transition weight distribution between different nodes, it obtains the connectivity probability value of each node, reflecting the significance of the node being located within the main vertebral canal. Based on the coordinate distribution of each node's adjacent nodes and the ending lumen coordinates, tissue density, connectivity probability value, and pixel calibration interval constant... The algorithm obtains the comprehensive search cost value of adjacent nodes, executes a path search algorithm to obtain a three-dimensional discrete coordinate sequence, and obtains a three-dimensional spinal cord center coordinate sequence. This ensures that the search wavefront can maintain a searchable normal cost value even in areas with severe tissue gray-level overlap, guaranteeing the continuity of main pathway pathfinding. Based on the coordinate distribution in the three-dimensional spinal cord center coordinate sequence, the pixel calibration spacing constant, and the gray-level distribution of adjacent voxels at corresponding positions in the three-dimensional image matrix for different center coordinates, the cross-sectional area of ​​the spinal cord lumen at each center coordinate is obtained. This eliminates physiological curvature errors, quantifies the lumen size at each location, and obtains the degree of spinal cord compression. This invention obtains accurate three-dimensional spinal cord center coordinates, enabling a more objective and precise assessment of the degree of spinal cord compression. Attached Figure Description

[0016] Figure 1 A flowchart illustrating an automatic assessment method for cervical spinal cord compression based on MRI images, provided as an embodiment of the present invention; Figure 2 A flowchart illustrating a method for obtaining the starting lumen coordinates and the ending lumen coordinates according to an embodiment of the present invention; Figure 3 This is a flowchart illustrating a method for obtaining a comprehensive search cost value according to an embodiment of the present invention. Detailed Implementation

[0017] Unless otherwise defined, all technical and scientific terms used herein have the same meaning as commonly understood by one of ordinary skill in the art to which this invention pertains.

[0018] The following describes in detail, with reference to the accompanying drawings, a specific scheme for an automatic assessment method for the degree of cervical spinal cord compression based on MRI images provided by the present invention.

[0019] Please see Figure 1 The diagram illustrates a flowchart of an automatic assessment method for cervical spinal cord compression based on MRI images, according to an embodiment of the present invention. The specific method includes: Step S1: Obtain the three-dimensional image matrix of the cervical spinal cord and the pixel calibration spacing constant in different dimensions.

[0020] In an embodiment of the present invention, the system first receives the original three-dimensional image matrix output by the medical imaging device and simultaneously reads the medical digital imaging file header label associated with the matrix. After reading the file header label, the system parses it and extracts the pixel calibration spacing constant in the three-dimensional coordinate system. That is, when extracting the pixel calibration spacing constant, the system obtains the pixel calibration spacing of the image in the X-axis, Y-axis and Z-axis dimensions respectively, which helps to fix the multiplier reference for all subsequent distance and area conversion formulas. The unit is millimeters.

[0021] It should be noted that the three-dimensional image matrix is ​​composed of small cubes with volume. Each cube is a voxel, which contains positional information, namely the three-dimensional coordinates of X, Y, and Z and the corresponding gray value. The three-dimensional image matrix is ​​a discrete grid, and the coordinates of all voxels are positive integers in the matrix index. In clinical MRI, the head-to-tail orientation is usually defined as the Z-axis.

[0022] Step S2: Obtain the tissue density of each voxel based on the grayscale characteristics of different voxels in the three-dimensional image matrix; obtain the starting lumen coordinates and ending lumen coordinates based on the tissue density distribution of voxels in two-dimensional slices at different positions on the Z-axis of the three-dimensional image matrix.

[0023] Considering the differences in overall brightness range of MRI images from different devices, the absolute device grayscale values ​​are transformed into a consistent quantitative standard characterizing local tissue density, maintaining the physical continuity of anatomical tissue grayscale in three-dimensional space. Based on the grayscale characteristics of different voxels in the three-dimensional image matrix, the tissue density of each voxel is obtained.

[0024] Preferably, in one embodiment of the present invention, the method for obtaining tissue density includes: Obtain the maximum and minimum gray values ​​of all voxels in the 3D image matrix; The gray values ​​of each voxel are normalized based on the maximum and minimum gray values, and the normalization results are negatively correlated and mapped as the tissue density of each voxel.

[0025] It should be noted that, in one embodiment of the present invention, normalization is achieved by calculating the difference between the positive integer 1 and the normalization result. The larger the normalization result, the smaller the tissue density, and the value ranges from 0 to 1. In other embodiments of the present invention, normalization can be achieved by taking the reciprocal or an exponential function with the natural constant as the base. A negative correlation mapping is performed, in which a very small positive number with the same dimensions is added to the denominator when calculating the reciprocal. The value of the number is set according to the range of the denominator. After the negative correlation mapping, the larger the normalized result, the smaller the tissue density. In order to keep the value in the range of 0-1, the negative correlation mapping result is normalized. The specific means are well known to those skilled in the art and will not be described in detail here.

[0026] It should be noted that, in the embodiments of the present invention, the gray value of each voxel is normalized to the range of 0-1 based on the maximum and minimum gray values, and the maximum and minimum value normalization method is adopted. The method includes: Calculate the first gray level difference between the gray level value and the minimum gray level value of each voxel to obtain the gray level range between the maximum gray level value and the minimum gray level value. Considering that the gray level range may be 0 when the image is completely black or completely white, obtain the sum of the gray level range value and the preset adjustment coefficient of the smallest positive number as the gray level range adjustment value. The preset adjustment coefficient can be set according to the value range and dimension of the gray level range value, which will not be elaborated here, so that the denominator is always not 0 during the calculation. The ratio of the first grayscale difference to the grayscale range adjustment value is calculated as the normalization result; the specific method is a well-known technique to those skilled in the art and will not be elaborated here.

[0027] Based on this, the smaller the tissue density, the larger the gray value of each voxel, and the closer the voxel is to the cerebrospinal fluid region that presents a bright water signal.

[0028] The formula for tissue density is expressed as: ;in, Indicates the first Tissue density of individual elements; Indicates the first The grayscale value of a single individual; This represents the minimum grayscale value; Indicates the maximum grayscale value; This indicates the preset adjustment coefficient; Indicates the grayscale range; This indicates the grayscale range adjustment value; This indicates that the grayscale values ​​are based on the maximum and minimum grayscale values ​​for the first... The gray values ​​of individual pixels are normalized to their maximum and minimum values.

[0029] Tissue density reflects the compactness of the tissue corresponding to the voxel. Considering that the lumen is filled with bright cerebrospinal fluid when it is not compressed, the higher the gray value, the lower the resistance area and the lower the tissue density. On the other hand, for the surrounding bone or protruding lesion tissue, the lower the gray value, the higher the resistance area and the higher the tissue density. Therefore, analyzing the tissue density distribution of voxels in two-dimensional slices at different locations can help locate the actual lumen. Based on the tissue density distribution of voxels in two-dimensional slices at different locations on the Z-axis of the three-dimensional image matrix, the coordinates of the starting lumen and the ending lumen are obtained.

[0030] Preferably, in one embodiment of the present invention, the method for obtaining the starting lumen coordinates and the ending lumen coordinates is described in [reference needed]. Figure 2 It illustrates a flowchart of a method for obtaining the starting and ending lumen coordinates, including: Step S201: Based on the tissue density distribution of voxels in two-dimensional slices at different positions on the Z-axis of the three-dimensional image matrix, obtain the main spinal canal connected domain of the two-dimensional slice at the corresponding position.

[0031] Preferably, in one embodiment of the present invention, the method for obtaining the main spinal canal connectivity region includes: For any location, if the tissue density of a voxel in the two-dimensional slice is less than the preset low-resistivity cutoff constant, the corresponding voxel will be used as a candidate voxel. Obtain the connected domains formed by consecutive adjacent candidate voxels, and select the connected domain with the most candidate voxels among all connected domains as the main vertebral canal connected domain.

[0032] It should be noted that the lower the tissue density, the more likely it is to be a voxel in a low-resistance state, and the more it tends to be high-brightness cerebrospinal fluid, the more it represents being analyzed within the lumen. The fewer candidate voxels, the more likely it is to be an isolated artifact region. Therefore, the connected region with the most candidate voxels is taken as the main spinal canal connected region. The main spinal canal connected region helps to reflect the true anatomical cross-section of the spinal canal filled with cerebrospinal fluid.

[0033] It should be noted that, in one embodiment of the present invention, considering that the tissue density of cerebrospinal fluid is closer to 0, in order to screen out high-brightness cerebrospinal fluid regions for analysis, the preset low-resistance cutoff constant is set to 0.1 based on relevant historical experience; in other embodiments of the present invention, the preset low-resistance cutoff constant can be set according to specific circumstances, and will not be limited or elaborated here.

[0034] Step S202: For any position, obtain the mean value of the position coordinates of all voxels in the main spinal canal connected domain on the two-dimensional slice, and use it as the average coordinate; obtain the relative distance between the position coordinates of different voxels in the main spinal canal connected domain and the average coordinate, and use the position coordinates of the voxel corresponding to the minimum relative distance as the average corrected coordinate.

[0035] The average coordinates reflect the geometric center in the main spinal canal's connected domain. The further the geometric center deviates from the main spinal canal's connected domain, the more severe the curvature of the lumen is. The geometric center may fall on high-resistance tissue outside the lumen, making it impossible to accurately analyze the cerebrospinal fluid lumen region. Therefore, adjustments are needed to find a center point with a valid integer, i.e., the average corrected coordinates.

[0036] It should be noted that, in the embodiments of the present invention, the relative distance is obtained by existing distance calculation methods such as Euclidean distance or Manhattan distance. The smaller the relative distance, the closer the spatial positions are, and the more real the integer coordinate value closest to the average coordinate is obtained for subsequent analysis. The specific means are well known to those skilled in the art and will not be described in detail here.

[0037] Step S203: The three-dimensional coordinates formed by the highest position of the Z-axis and the average corrected coordinates on the corresponding two-dimensional slice are used as the starting lumen coordinates; the three-dimensional coordinates formed by the lowest position of the Z-axis and the average corrected coordinates on the corresponding two-dimensional slice are used as the ending lumen coordinates.

[0038] Considering that the two ends of the cervical spine are usually not subjected to severe compression, the true opening cross section of the lumen can be extracted most stably. Also, since the cervical spine scan is stacked layer by layer from the top of the head to the tailbone, the highest and lowest positions of the Z-axis correspond to the first and last layers of the scan, i.e. the two ends of the cervical spine. Therefore, the analysis is performed at the highest and lowest positions to obtain the coordinates of the starting lumen and the ending lumen.

[0039] Step S3: Obtain a 3D sparse graph structure composed of voxels as nodes. Based on the tissue density distribution of all nodes, obtain the state transition weights between different nodes. Based on the linear index positions of the starting and ending lumen coordinates, the preset restart probability constant, and the state transition weight distribution between different nodes and adjacent nodes, obtain the new probability distribution vector of the entire space after iteration, and construct the corresponding 3D spatial connectivity probability matrix.

[0040] The lower the tissue density of the brain and spinal cord in low-resistance areas, the greater the diffusion probability; conversely, the higher the tissue density in high-resistance compressed areas, the lower the diffusion probability. Therefore, based on the tissue density analysis of different nodes, the degree of transfer to adjacent nodes is reflected. A three-dimensional sparse graph structure composed of voxels as nodes is obtained, and the state transition weights between different nodes are obtained according to the tissue density distribution of all nodes.

[0041] Preferably, in one embodiment of the present invention, the method for obtaining the state transition weights includes: For any node, obtain the difference between the positive integer 1 and the tissue density of each node's corresponding neighboring nodes, calculate the sum of the difference result and the preset adjustment coefficient, and use it as the state transition weight of each node relative to its neighboring nodes. It should be noted that each node is connected to 26 neighboring nodes in the three-dimensional space, which are called adjacent nodes. Each adjacent node is analyzed in turn to obtain the state transition weight of each node relative to any adjacent node.

[0042] It should be noted that, in the embodiments of the present invention, in order to avoid the situation where the state transition weight is meaningless when the tissue compaction is 1, the preset adjustment coefficient is set to a very small positive number with consistent dimensions. Its value can be specifically set according to the value range of tissue compaction, such as 0.001, to retain weak spatial connectivity. The greater the tissue compaction of adjacent nodes, the stronger the positional resistance, the lower the probability of a node transferring to an adjacent node, and the smaller the state transition weight of a node relative to its adjacent nodes.

[0043] For nodes other than their neighbors, set the state transition weight of each node relative to the other nodes to 0.

[0044] Considering that the blind ends of nerve roots extending to both sides of the cervical spine in medical imaging are filled with cerebrospinal fluid and exhibit the same low resistance characteristics as the main spinal canal, and that the state transition weights between different nodes are similar, the restart distribution vector and the restart probability constant in the transition are combined to enable the main spinal canal region near the double-ended entrance and through to accumulate a more significant connectivity probability than the deep blind ends. Based on the one-dimensional initial restart distribution vector formed by the coordinates of the starting lumen and the ending lumen, the preset restart probability constant, and the state transition weight distribution between different nodes, the connectivity probability value of each node is obtained.

[0045] Preferably, in one embodiment of the present invention, the method for obtaining the connectivity probability value includes: Step 1: Construct a one-dimensional initial restart distribution vector for the number of nodes, obtain the one-dimensional linear index positions of the starting and ending lumen coordinates, set the element values ​​of the corresponding positions of the one-dimensional initial restart distribution vector to preset probability values, set the element values ​​of other positions to 0, and normalize the state transition weights between different nodes to form a global state transition sparse matrix between different nodes.

[0046] It should be noted that, based on existing technology, the conversion from 3D coordinates to 1D indexes is as follows: using the common row-by-row storage method, X-coordinates are the fastest, followed by Y-coordinates, and Z-coordinates are the slowest. The resulting 1D index formula is: Where, index represents the one-dimensional index value; X represents the X-coordinate value in three dimensions; This represents the Y-coordinate value in three dimensions; Represents the Z-coordinate value in three dimensions; Represents the width in the X dimension; Represents the height of the Y dimension; such as the X dimension in a three-dimensional matrix. 2, Y dimension The coordinates are 5, the Z dimension is 2, and the number of nodes is 2×5×2, which is 20 nodes; The mapping to a one-dimensional index is 9.

[0047] It should be noted that only by continuously injecting the probability from the two ends of the anatomically determined main vertebral canal, i.e., the start and end coordinates, can the system utilize this bidirectional water flow to penetrate the pressure zone and identify the true main channel. Therefore, in the embodiments of the present invention, the initial probability is equally divided by the one-dimensional linear index positions corresponding to the starting and ending lumen coordinates, and the preset probability value is set to 0.5.

[0048] It should be noted that, in one embodiment of the present invention, the normalization method is as follows: for any node, the sum of the state transition weights of the node relative to all its neighboring nodes is obtained as the overall weight level; the ratio of the state transition weight of the node relative to each neighboring node to the overall weight level is obtained and normalized.

[0049] It should be noted that the dimension of the one-dimensional initial restart distribution vector is the number of nodes × 1, and the dimension of the global state transition sparse matrix is ​​the number of nodes × the number of nodes.

[0050] Step 2: Using the one-dimensional initial restart distribution vector as the initial probability distribution vector, obtain the product matrix of the global state transition sparse matrix and the initial probability distribution vector. Obtain the difference between the positive integer 1 and the preset restart probability constant. Calculate the product of the difference result and the product matrix, and use it as the transition probability vector. Obtain the product of the preset restart probability constant and the one-dimensional initial restart distribution vector, and use it as the restart probability vector. Obtain the sum of the transition probability vector and the restart probability vector, and use it as the probability distribution vector for the next iteration.

[0051] The formula is expressed as: ;in, Indicates the first The probability distribution vector of the next iteration; This represents the preset restart probability constant; Represents the sparse matrix of global state transitions; Indicates the first The probability distribution vector of the next iteration; This represents a one-dimensional initial restart distribution vector.

[0052] It should be noted that for the initial iteration When =0, the one-dimensional initial restart distribution vector is used as the initial probability distribution vector for the next iteration; the dimension of the product matrix of the global state transition sparse matrix and the initial probability distribution vector is the number of nodes × 1.

[0053] It should be noted that, in one embodiment of the present invention, considering that if the preset restart probability constant is set too high, the probability will be concentrated near the starting point and unable to effectively penetrate the severely compressed area of ​​the main channel, in order to ensure that the probability of injection at both ends of the main spinal canal can be steadily concentrated on the compressed segment, the preset restart probability constant is set to a small positive number based on the average physical length of the cervical spine anatomy and the spatial resolution of the imaging equipment, such as 0.15 based on relevant historical experience; in other embodiments of the present invention, the size of the preset restart probability constant can be set according to specific circumstances, and will not be limited or elaborated here.

[0054] Step 3: Obtain the difference vector between the probability distribution vector of the next iteration and the initial probability distribution vector. Calculate the sum of the absolute values ​​of all elements in the difference vector. If the sum of the absolute values ​​is greater than or equal to the preset convergence threshold, use the probability distribution vector of the next iteration as the new initial probability distribution vector to obtain the probability distribution vector of the corresponding next iteration. Continue until the sum of the absolute values ​​is less than the preset convergence threshold. Map the probability distribution vector of the next iteration to the three-dimensional space to obtain the connectivity probability value of each node.

[0055] It should be noted that, in one embodiment of the present invention, if the probability distribution vectors of adjacent iterations change similarly, it indicates that the probability distribution in the space has reached a stable dynamic equilibrium, and the iteration loop is stopped. At this time, the elements in the difference vector are all close to 0, therefore the preset convergence threshold is set to a very small positive number, such as... The specific settings can be adjusted according to the specific circumstances, and will not be limited or elaborated here.

[0056] It should be noted that mapping the probability distribution vector of the next iteration onto three-dimensional space to obtain the connectivity probability value of each node is done by following the reverse spatial order rule when mapping from three dimensions to one dimension. , , Where, index represents the one-dimensional index value; X represents the X-coordinate value in three dimensions; This represents the Y-coordinate value in three dimensions; Represents the Z-coordinate value in three dimensions; Represents the width in the X dimension; The modulo operator represents the height in the Y dimension; mod represents the remainder. This indicates rounding down to the nearest integer.

[0057] Step S4: Based on the coordinate distribution of each node's neighboring nodes and the coordinates of the terminating lumen, tissue density, connectivity probability, and pixel calibration spacing constant, obtain the comprehensive search cost of neighboring nodes; based on the comprehensive search cost of different neighboring nodes between the starting lumen coordinates and the terminating lumen coordinates, execute the path search algorithm to obtain a three-dimensional discrete coordinate sequence, and obtain a three-dimensional spinal cord center coordinate sequence.

[0058] When the search algorithm wave encounters high pressure from tissue compression in the main spinal canal, it is very easy to deviate from the path to the blind end of the transverse nerve root with the same low resistance gray level. Therefore, combining tissue density and connectivity probability value helps to assign a smaller cost value to nodes with higher connectivity probability. The pixel calibration spacing constant represents the actual physical millimeter length, which helps to quantify the actual distance. Therefore, based on the coordinate position distribution of each node's adjacent nodes and the terminating lumen coordinates, tissue density, connectivity probability value, and pixel calibration spacing constant, the comprehensive search cost value of adjacent nodes is obtained.

[0059] Preferably, in one embodiment of the present invention, the method for obtaining the comprehensive search cost value is described in [reference needed]. Figure 3 It illustrates a flowchart of a method for obtaining a comprehensive search cost value, including: Step S301: For any two nodes, obtain the product of the difference in coordinates between the nodes in each dimension and the pixel calibration spacing constant, as the actual spacing in each dimension; obtain the sum of the squares of the actual spacings in all dimensions and take the square root, as the actual distance between the corresponding nodes.

[0060] It should be noted that the sum of the squares of the actual distances in all dimensions is obtained, and the square root is taken to represent the actual physical Euclidean straight-line distance in three-dimensional space.

[0061] Step S302: Based on the organization compactness, connectivity probability value, and actual distance between the node and its neighboring nodes, obtain the single-step movement cost increment between the node and its neighboring nodes. Organization compactness and actual distance are positively correlated with the single-step movement cost increment, while connectivity probability value is negatively correlated with the single-step movement cost increment.

[0062] It should be noted that the greater the tissue compactness, the greater the resistance and the greater the cost of movement; the greater the spatial connectivity probability, the easier the transfer and the lower the cost of movement; and the greater the actual distance, the greater the cost of movement. Therefore, tissue compactness and actual distance are both positively correlated with the increment of single-step movement cost, while connectivity probability is negatively correlated with the increment of single-step movement cost.

[0063] In one embodiment of the present invention, a first sum between tissue compactness and a preset first adjustment coefficient is obtained, a second sum between the connectivity probability value of adjacent nodes and a preset second adjustment coefficient is obtained, the ratio of the first sum and the second sum is calculated as a first cost coefficient, and the product of the first cost coefficient and the actual distance is obtained as the single-step movement cost increment.

[0064] The formula is expressed as: ;in, Indicates the first The node and the first The incremental cost of a single-step movement between adjacent nodes; Indicates the first Tissue compactness of adjacent nodes; Indicates the first The connectivity probability value of each adjacent node; This indicates the preset first adjustment coefficient; This indicates the preset second adjustment coefficient; Indicates the first The node and the first The actual distance between adjacent nodes.

[0065] Based on this, when the search wavefront touches a transverse nerve root branch, although the tissue density... Extremely small, exhibiting low-resistance cerebrospinal fluid, but due to its topological blindness, the connectivity probability value is... Even smaller, extremely small denominator This causes the cost after division to be amplified exponentially, thus transforming the low-resistance anatomical blind end into a high-cost barrier at the numerical calculation level, forcibly limiting the lateral spread of errors in the path, and the greater the increment in cost per step; conversely, in the spinal canal lesion segment where severe tissue compression occurs, although tissue density... The larger the value, the greater the local resistance; however, because the main road maintains a relatively high connectivity probability during the random walk, this resistance remains relatively high. The numerical advantage of the denominator effectively mitigated the surge in resistance.

[0066] It should be noted that, in the embodiments of the present invention, in order to prevent the numerator and denominator from being 0 and the formula from being meaningless, a preset first adjustment coefficient and a preset second adjustment coefficient are added. The preset first adjustment coefficient and the preset second adjustment coefficient are set to extremely small positive numbers, and their values ​​are specifically set according to the value range of the corresponding numerator or denominator, such as 0.01. This is not limited or elaborated here.

[0067] Step S303: Obtain the actual distance between the corresponding adjacent node and the coordinates of the termination cavity, as the heuristic cost; obtain the cumulative movement cost from the starting cavity coordinates to the current node; obtain the sum of the cumulative movement cost, the heuristic cost of adjacent nodes, and the single-step movement cost increment, as the comprehensive search cost of adjacent nodes.

[0068] It should be noted that, in the embodiments of the present invention, the cumulative movement cost of the starting cavity coordinates is set to zero; during the progressive expansion of the search, the cumulative movement cost of the currently traversed node is obtained through the priority queue popping mechanism of the graph search algorithm.

[0069] Based on this, the comprehensive search cost value can reflect the actual cumulative cost value from the starting lumen coordinates through the current node to the subsequent adjacent nodes, as well as the sum of the heuristic cost value predicted from the adjacent nodes to the terminating lumen coordinates. The larger the cost value, the larger the distance, the smaller the connectivity probability, the larger the tissue density, and the less likely it is to be located at the terminating lumen position.

[0070] The comprehensive search cost reflects the combined constraints of local tissue penetration resistance and global topological connectivity security during path expansion. A smaller comprehensive search cost indicates safer and smoother movement of a node to adjacent nodes, and a greater tendency to proceed within the main spinal canal. Based on the comprehensive search cost of different adjacent nodes between the starting and ending lumen coordinates, a path search algorithm is executed to obtain a three-dimensional discrete coordinate sequence, and a three-dimensional spinal cord center coordinate sequence is also obtained.

[0071] Preferably, in one embodiment of the present invention, obtaining a three-dimensional discrete coordinate sequence and obtaining a three-dimensional spinal cord center coordinate sequence includes: The search origin is set at the starting lumen coordinates, and the ending lumen coordinates are set at the ending lumen coordinates. The search origin is added to the candidate node set. Each time, the node with the smallest comprehensive search cost is selected from the candidate node set for expansion. If the corresponding node is the ending point, the search ends. Otherwise, its adjacent unvisited nodes are added to the candidate node set, and the cumulative movement cost of the adjacent nodes is updated until the search ending point is reached. The search origin is backtracked from the search ending point to obtain a three-dimensional discrete coordinate sequence consisting of all the nodes passed through. It should be noted that the search algorithm used is the A-Star graph search algorithm, and the specific methods are well known to those skilled in the art. Considering that when encountering severe image artifacts that cause a single-layer full tomography, the search algorithm may get stuck in infinite expansion and cause memory overflow, a maximum expansion node limit is preset. The value is set according to the specific situation, such as 20% of the total number of nodes. That is, if the number of nodes traversed exceeds the maximum expansion node limit, the search endpoint has not been reached, and the search process needs to be terminated for manual intervention and warning. The specific search algorithm is a well known technique to those skilled in the art, and will not be described in detail here.

[0072] A sliding filter is applied to the three-dimensional discrete coordinate sequence, and the filtered three-dimensional discrete coordinate sequence is used as the three-dimensional spinal cord center coordinate sequence.

[0073] It should be noted that, in the embodiments of the present invention, the method for obtaining the sliding filter is as follows: First, a sliding window is constructed to traverse the three-dimensional discrete coordinate sequence and obtain the mean of all three-dimensional discrete coordinates in the sliding window as the local smoothing coordinate; It should be noted that, in one embodiment of the present invention, the size of the sliding window is 5; in other embodiments of the present invention, the size of the sliding window can be set according to specific circumstances, and will not be limited or described in detail here.

[0074] The second step is to perform spatial interpolation based on tissue density at different location coordinates to obtain the fitted tissue density at local smooth coordinates. It should be noted that, in the embodiments of the present invention, a trilinear interpolation algorithm can be used for spatial interpolation. The specific means are well known to those skilled in the art and will not be described in detail here.

[0075] Third, if the fitted tissue density is greater than the preset out-of-bounds warning threshold, the length of the sliding window is decreased by a preset step size to obtain new local smooth coordinates within the sliding window, until the sliding window length reaches the preset minimum window length, or until the fitted tissue density of the local smooth coordinates is less than or equal to the preset out-of-bounds warning threshold; if the fitted tissue density of the local smooth coordinates is less than or equal to the preset out-of-bounds warning threshold, the three-dimensional discrete coordinates of the current node are replaced with local smooth coordinates; if the preset minimum window length is reached and the corresponding fitted tissue density is still greater than the preset out-of-bounds warning threshold, the three-dimensional discrete coordinates of the current node are retained.

[0076] Based on this, the three-dimensional spinal cord center coordinate sequence not only possesses the mathematical properties of spatial continuity, smoothness, and differentiability, but also strictly encloses the actual anatomical cavity, providing a reliable geometric framework for subsequent reconstruction of spatial orthogonal slices.

[0077] It should be noted that, in one embodiment of the present invention, considering that the greater the tissue density, the greater the resistance, and the more necessary it is to adjust the coordinates, the preset step size is set to 1, and the length of the sliding window is reduced by 1 each time for re-filtering; the preset out-of-bounds warning threshold is set to 0.8; in engineering practice, a sliding window length of 0 has no physical meaning and will lead to division by zero errors or memory out-of-bounds crashes, so the preset minimum window length is set to 1; in other embodiments of the present invention, the size of the preset out-of-bounds warning threshold can be set according to specific circumstances, and is not limited or elaborated here.

[0078] Step S5: Based on the coordinate distribution in the three-dimensional spinal cord center coordinate sequence, the pixel calibration spacing constant, and the gray-level distribution of adjacent voxels at corresponding positions in the three-dimensional image matrix for different center coordinates, obtain the cross-sectional area of ​​the spinal cord canal for each center coordinate and the degree of spinal cord compression.

[0079] The human cervical spine exhibits a distinct physiological lordosis. If the lumen contour is directly extracted on the global absolute horizontal plane, the tilted anatomical lumen will undergo elliptical stretching deformation on the two-dimensional plane, resulting in an inflated measured area. Therefore, analysis along the centerline trajectory and pixel calibration spacing constant is necessary to construct a cross-section that is strictly perpendicular to the actual local orientation of the lumen. Since the acquired image matrix only supports integer index reading, directly using the center coordinates to read grayscale will result in severely jagged tomographic slices, disrupting the smooth boundaries of soft tissues. Therefore, analyzing the grayscale distribution of adjacent voxels is crucial for more accurately quantifying the cross-sectional area of ​​the spinal canal.

[0080] Preferably, in one embodiment of the present invention, the method for obtaining the cross-sectional area of ​​the spinal canal includes: Step 1: Construct an orthogonal cross-sectional mesh along the three-dimensional spinal cord center coordinate sequence, and perform grayscale interpolation on the orthogonal cross-sectional mesh using three-dimensional image matrix interpolation to generate a two-dimensional orthogonal slice image; It should be noted that, in the embodiments of the present invention, the method for obtaining two-dimensional orthogonal slice images is as follows: First, for the starting center coordinates or ending center coordinates in the three-dimensional spinal cord center coordinate sequence, obtain the coordinate difference vector between the subsequent center coordinates and the starting center coordinates, and obtain the coordinate difference vector between the previous center coordinates and the ending center coordinates. The second step is to obtain the coordinate differences between the next central coordinate and the previous central coordinate in each dimension for other central coordinates in the three-dimensional spinal cord central coordinate sequence; multiply the coordinate differences in each dimension by the corresponding pixel calibration spacing constant to obtain the physical coordinate difference vector; normalize the physical coordinate difference vector of each central coordinate to use it as the unit tangent vector of each central coordinate. It should be noted that, in the embodiments of the present invention, in order to avoid the modulus value being 0, the modulus value and the sum of the minimum positive number are calculated to obtain the modulus value. The minimum positive number can be specifically set according to the range of the modulus value. The physical coordinate difference of each center coordinate is obtained by dividing it by the sum of the modulus value and the minimum positive number, that is, normalization is performed, and it is used as the unit tangent vector of each center coordinate.

[0081] The third step involves constructing a two-dimensional discrete grid with the center coordinates as the geometric center and the corresponding unit tangent vector as the normal vector. This yields a vertical cross-sectional grid formed by the two-dimensional discrete grid in three-dimensional space. The average gray value of the adjacent voxels of each grid point in each vertical cross-sectional grid is then used as the average gray level of each grid point, thus forming a two-dimensional orthogonal slice image.

[0082] Step 2: Obtain the grayscale variance of all pixels in each two-dimensional orthogonal slice image. If the grayscale variance is less than a preset variance threshold, set the cross-sectional area of ​​the spinal canal of the corresponding two-dimensional orthogonal slice image to 0. It should be noted that, in the embodiments of the present invention, if based on the 0-255 standard grayscale, a normal T2WI slice contains obvious bright cerebrospinal fluid and dark tissue, and its global variance is usually between 2000 and 3500. If a single-peak image of severe compression leading to cerebrospinal fluid drainage occurs, the variance usually drops sharply to below 500. Therefore, the preset variance threshold is set to 500. In other embodiments of the present invention, the preset variance threshold can be set according to specific circumstances, which is not limited or elaborated here.

[0083] Step 3: If the gray-level variance is greater than or equal to the preset variance threshold, the Otsu algorithm is used to adaptively segment the two-dimensional orthogonal slice image and perform binarization to obtain the connected components formed by adjacent pixels after binarization. The product of the pixel calibration spacing constants in two dimensions on the horizontal plane is obtained as the pixel calibration physical area; the product of the number of pixels within the closed lumen contour and the pixel calibration physical area is obtained as the cross-sectional area of ​​the spinal canal in the corresponding two-dimensional orthogonal slice image.

[0084] It should be noted that, in the embodiments of the present invention, the binarization processing method is as follows: if the gray value of a pixel is less than the Otsu threshold, the corresponding pixel value is set to a binary value of 0; conversely, the corresponding pixel value is set to a binary value of 1; the two-dimensional connected component analysis algorithm is called to obtain the connected component formed by adjacent pixels; the Otsu algorithm automatically determines the optimal threshold by finding the maximum inter-class variance, which is used as the Otsu threshold; the specific means are well known to those skilled in the art and will not be described in detail here.

[0085] Preferably, in one embodiment of the present invention, the method for obtaining the degree of spinal cord compression includes: Obtain the spinal cord lumen cross-sectional area sequence corresponding to the three-dimensional spinal cord center coordinate sequence, and obtain the minimum value of the element in the spinal cord lumen cross-sectional area sequence as the lumen cross-sectional area of ​​the most severe compression point; Within the neighborhood of the cross-sectional area of ​​the lumen at the most severe compression point, the maximum value in the sequence of all spinal cord lumen cross-sectional areas is obtained as the local healthy baseline area. It should be noted that, in one embodiment of the present invention, the method for obtaining the neighborhood range is as follows: based on the cross-sectional area of ​​the lumen at the most severe compression point, a range consisting of 10 points before and after the cross-sectional area of ​​the spinal cord lumen is selected; in other embodiments of the present invention, the size of the neighborhood range can be specifically set according to the specific situation, and will not be limited or described in detail here.

[0086] The degree of spinal cord compression is obtained based on the cross-sectional area of ​​the lumen at the point of most severe compression and the local healthy baseline area. The cross-sectional area of ​​the lumen at the point of most severe compression is negatively correlated with the degree of spinal cord compression, while the local healthy baseline area is positively correlated with the degree of spinal cord compression.

[0087] It should be noted that the cross-sectional area of ​​the lumen at the point of most severe compression represents the actual physical area of ​​the lumen remaining at the extreme point of the lesion; the local healthy baseline area represents the physical area of ​​the reference lumen under normal conditions near the lesion point; therefore, the larger the local healthy baseline area, the smaller the cross-sectional area of ​​the lumen at the point of most severe compression, indicating that the compression at the lesion point is more severe and the degree of spinal cord compression is greater than that under normal conditions.

[0088] In one embodiment of the present invention, considering that the local health reference area may be 0, the sum of the local health reference area and the minimum positive number is calculated. The minimum positive number can be specifically set according to the value range of the local health reference area, such as... The first ratio is obtained by dividing the cross-sectional area of ​​the lumen at the point of most severe compression by the sum of the local healthy baseline area and the smallest positive number. The difference between the positive integer 1 and the first ratio is calculated as the degree of spinal cord compression.

[0089] Based on this, the degree of spinal cord compression directly and objectively characterizes the relative severity of deformation in the target cervical spine region; the greater the degree of spinal cord compression, the greater the severity of deformation.

[0090] In summary, this invention obtains a three-dimensional sparse graph structure composed of voxels as nodes. Based on the tissue density distribution of all nodes and the one-dimensional initial restart distribution vector composed of the starting and ending lumen coordinates, the connectivity probability value of each node is obtained. Based on the coordinate distribution of each node's neighboring nodes and the ending lumen coordinates, tissue density, connectivity probability value, and pixel calibration spacing constant, the comprehensive search cost of neighboring nodes is obtained. A path search algorithm is then executed to obtain a three-dimensional spinal cord center coordinate sequence. Based on the coordinate distribution in the three-dimensional spinal cord center coordinate sequence and the grayscale distribution of different voxels in the three-dimensional image matrix, the degree of spinal cord compression is obtained. This invention obtains accurate three-dimensional spinal cord center coordinates, enabling a more objective and precise assessment of the degree of spinal cord compression.

[0091] It should be noted that, in other embodiments of the present invention, based on the same application concept as the automatic assessment method for cervical spinal cord compression degree based on MRI images provided in the embodiments of this application, an automatic assessment system for cervical spinal cord compression degree based on MRI images is also proposed. The system includes a memory, a processor, and a computer program stored in the memory and executable on the processor. When the processor executes the computer program, it implements the steps of an automatic assessment method for cervical spinal cord compression degree based on MRI images.

Claims

1. An automatic assessment method for the degree of cervical spinal cord compression based on MRI images, characterized in that, The method includes: Obtain the three-dimensional image matrix of the cervical spinal cord and the pixel calibration spacing constant in different dimensions; Based on the grayscale characteristics of different voxels in the 3D image matrix, the tissue density of each voxel is obtained; based on the tissue density distribution of voxels in 2D slices at different positions on the Z-axis of the 3D image matrix, the coordinates of the starting lumen and the ending lumen are obtained. A three-dimensional sparse graph structure composed of voxels as nodes is obtained. Based on the tissue density distribution of all nodes, the state transition weights between different nodes are obtained. Based on the one-dimensional initial restart distribution vector composed of the starting lumen coordinates and the ending lumen coordinates, the preset restart probability constant, and the state transition weight distribution between different nodes, the connectivity probability value of each node is obtained. Based on the coordinate distribution of each node's neighboring nodes and the coordinates of the terminating lumen, tissue density, connectivity probability, and pixel calibration spacing constant, the comprehensive search cost of neighboring nodes is obtained. Based on the comprehensive search cost of different neighboring nodes between the starting and ending lumen coordinates, a path search algorithm is executed to obtain a three-dimensional discrete coordinate sequence and a three-dimensional spinal cord center coordinate sequence. Based on the coordinate distribution in the three-dimensional spinal cord center coordinate sequence, the pixel calibration spacing constant, and the gray-level distribution of adjacent voxels at corresponding positions in the three-dimensional image matrix for different center coordinates, the cross-sectional area of ​​the spinal cord canal at each center coordinate is obtained, and the degree of spinal cord compression is obtained.

2. The automatic assessment method for cervical spinal cord compression based on MRI images according to claim 1, characterized in that, The method for obtaining tissue density includes: Obtain the maximum and minimum gray values ​​of all voxels in the 3D image matrix; The gray values ​​of each voxel are normalized based on the maximum and minimum gray values, and the normalization results are negatively correlated and mapped as the tissue density of each voxel.

3. The automatic assessment method for cervical spinal cord compression based on MRI images according to claim 1, characterized in that, The methods for obtaining the starting lumen coordinates and the ending lumen coordinates include: Based on the tissue density distribution of voxels in two-dimensional slices at different positions on the Z-axis of the three-dimensional image matrix, the main vertebral canal connectivity domain of the two-dimensional slices at the corresponding positions is obtained. For any given location, the mean value of the position coordinates of all voxels within the main spinal canal connected domain on the two-dimensional slice is obtained as the average coordinate; the relative distance between the position coordinates of different voxels within the main spinal canal connected domain and the average coordinate is obtained, and the position coordinate of the voxel corresponding to the minimum relative distance is used as the average corrected coordinate. The three-dimensional coordinates formed by the highest position of the Z-axis and the average corrected coordinates on the corresponding two-dimensional slice are used as the starting lumen coordinates; the three-dimensional coordinates formed by the lowest position of the Z-axis and the average corrected coordinates on the corresponding two-dimensional slice are used as the ending lumen coordinates.

4. The automatic assessment method for cervical spinal cord compression based on MRI images according to claim 3, characterized in that, The method for obtaining the main spinal canal connectivity region includes: For any location, if the tissue density of a voxel in the two-dimensional slice is less than the preset low-resistivity cutoff constant, the corresponding voxel will be used as a candidate voxel. Obtain the connected domains formed by consecutive adjacent candidate voxels, and select the connected domain with the most candidate voxels among all connected domains as the main vertebral canal connected domain.

5. The automatic assessment method for cervical spinal cord compression based on MRI images according to claim 1, characterized in that, The method for obtaining the state transition weights includes: For any node, obtain the difference between the positive integer 1 and the tissue density of each node's corresponding neighboring nodes, calculate the sum of the difference result and the preset adjustment coefficient, and use it as the state transition weight of each node relative to its neighboring nodes. For nodes other than their neighbors, set the state transition weight of each node relative to the other nodes to 0.

6. The automatic assessment method for cervical spinal cord compression based on MRI images according to claim 1, characterized in that, The method for obtaining the connectivity probability value includes: Construct a one-dimensional initial restart distribution vector for the number of nodes, obtain the one-dimensional linear index positions of the starting and ending lumen coordinates, set the element values ​​of the corresponding positions of the one-dimensional initial restart distribution vector to preset probability values, set the element values ​​of other positions to 0, normalize the state transition weights between different nodes, and construct a global state transition sparse matrix between different nodes. Using the one-dimensional initial restart distribution vector as the initial probability distribution vector, we obtain the product matrix of the global state transition sparse matrix and the initial probability distribution vector. We obtain the difference between the positive integer 1 and the preset restart probability constant, and calculate the product of the difference result and the product matrix as the transition probability vector. We obtain the product of the preset restart probability constant and the one-dimensional initial restart distribution vector as the restart probability vector. We obtain the sum of the transition probability vector and the restart probability vector as the probability distribution vector for the next iteration. Obtain the difference vector between the probability distribution vector of the next iteration and the initial probability distribution vector. Calculate the sum of the absolute values ​​of all elements in the difference vector. If the sum of the absolute values ​​is greater than or equal to a preset convergence threshold, use the probability distribution vector of the next iteration as the new initial probability distribution vector. Obtain the probability distribution vector for the corresponding next iteration. Continue until the sum of the absolute values ​​is less than the preset convergence threshold. Map the probability distribution vector of the next iteration onto a three-dimensional space to obtain the connectivity probability value of each node.

7. The automatic assessment method for cervical spinal cord compression based on MRI images according to claim 1, characterized in that, The method for obtaining the comprehensive search cost value includes: For any two nodes, obtain the product of the difference in coordinates between the nodes in each dimension and the pixel calibration spacing constant, as the actual spacing in each dimension; obtain the sum of the squares of the actual spacings in all dimensions and take the square root, as the actual distance between the corresponding nodes. Based on the organization compactness, connectivity probability value, and actual distance between a node and its neighboring nodes, the single-step movement cost increment between a node and its neighboring nodes is obtained. Organization compactness and actual distance are positively correlated with the single-step movement cost increment, while connectivity probability value is negatively correlated with the single-step movement cost increment. Obtain the actual distance between the corresponding adjacent node and the coordinates of the termination cavity as the heuristic cost; obtain the cumulative movement cost from the starting cavity coordinates to the current node; obtain the sum of the cumulative movement cost, the heuristic cost of adjacent nodes, and the single-step movement cost increment as the comprehensive search cost of adjacent nodes.

8. The automatic assessment method for cervical spinal cord compression based on MRI images according to claim 1, characterized in that, The process of obtaining the three-dimensional discrete coordinate sequence and the three-dimensional spinal cord center coordinate sequence includes: The search origin is set at the starting lumen coordinates, and the ending lumen coordinates are set at the ending lumen coordinates. The search origin is added to the candidate node set. Each time, the node with the smallest comprehensive search cost is selected from the candidate node set for expansion. If the corresponding node is the ending point, the search ends. Otherwise, its adjacent unvisited nodes are added to the candidate node set, and the cumulative movement cost of the adjacent nodes is updated until the search ending point is reached. The search origin is backtracked from the search ending point to obtain a three-dimensional discrete coordinate sequence consisting of all the nodes passed through. A sliding filter is applied to the three-dimensional discrete coordinate sequence, and the filtered three-dimensional discrete coordinate sequence is used as the three-dimensional spinal cord center coordinate sequence.

9. The automatic assessment method for cervical spinal cord compression based on MRI images according to claim 1, characterized in that, The method for obtaining the cross-sectional area of ​​the spinal canal includes: An orthogonal cross-sectional mesh is constructed along the three-dimensional spinal cord center coordinate sequence. The orthogonal cross-sectional mesh is then subjected to grayscale interpolation using a three-dimensional image matrix interpolation to generate a two-dimensional orthogonal slice image. Obtain the grayscale variance of all pixels in each two-dimensional orthogonal slice image. If the grayscale variance is less than the preset variance threshold, set the cross-sectional area of ​​the spinal canal of the corresponding two-dimensional orthogonal slice image to 0. If the gray-level variance is greater than or equal to the preset variance threshold, the Otsu algorithm is used to adaptively segment the two-dimensional orthogonal slice image and perform binarization to obtain the connected component formed by adjacent pixels after binarization. A connected region with a central coordinate is selected as the closed lumen contour; the product of the pixel calibration spacing constants in two dimensions on the horizontal plane is obtained as the pixel calibration physical area; the product of the number of pixels in the closed lumen contour and the pixel calibration physical area is obtained as the cross-sectional area of ​​the spinal canal of the corresponding two-dimensional orthogonal slice image.

10. The automatic assessment method for cervical spinal cord compression based on MRI images according to claim 1, characterized in that, The methods for obtaining the degree of spinal cord compression include: Obtain the spinal cord lumen cross-sectional area sequence corresponding to the three-dimensional spinal cord center coordinate sequence, and obtain the minimum value of the element in the spinal cord lumen cross-sectional area sequence as the lumen cross-sectional area of ​​the most severe compression point; Within the neighborhood of the cross-sectional area of ​​the lumen at the most severe compression point, the maximum value in the sequence of all spinal cord lumen cross-sectional areas is obtained as the local healthy baseline area. The degree of spinal cord compression is obtained based on the cross-sectional area of ​​the lumen at the point of most severe compression and the local healthy baseline area. The cross-sectional area of ​​the lumen at the point of most severe compression is negatively correlated with the degree of spinal cord compression, while the local healthy baseline area is positively correlated with the degree of spinal cord compression.

Images