Medical image analysis method and system based on image processing

Abstract

The invention relates to the technical field of medical image detection, in particular to a medical image analysis method and system based on image processing, and the method comprises the following steps: collecting a diffusion tensor image, carrying out the filtering preprocessing, normalizing the main diffusion direction, extracting an energy super-threshold voxel, calculating the node distance offset, and carrying out the correlation analysis to generate an atlas. And interpolating and aligning the brain return boundary to obtain a coupling graph, and outputting an abnormal probability distribution graph in a classified manner. According to the method, fiber details are ensured by adopting anisotropic diffusion filtering noise reduction of a diffusion tensor image sequence, three-dimensional normalization and wavelet energy analysis are coordinated, node detection specificity is enhanced by dynamic threshold screening, space-time coupling is quantified by combining displacement Euclidean distance and time sequence relevance, and an individualized coupling model is established by interpolating and fusing brain return boundaries; the template rigid constraint is broken through, the support vector machine recognizes an abnormal mode, structural distortion and functional abnormality cross-scale detection is achieved, and the early lesion image marker recognition efficiency is improved.

Description

Technical Field

[0001] This invention relates to the field of medical image detection technology, and in particular to medical image analysis methods and systems based on image processing. Background Technology

[0002] The field of medical image detection technology involves a comprehensive technical system for the automatic identification, quantitative evaluation, and morphological analysis of medical images related to the nervous system using image acquisition, processing, and analysis methods. This technical field encompasses the acquisition of image data of the central and peripheral nervous systems through medical imaging equipment such as magnetic resonance imaging, computed tomography, and positron emission tomography, followed by preprocessing work such as image denoising, registration, and standardization. Further key steps include feature extraction, image segmentation, region labeling, and quantitative analysis of brain tissue, spinal cord, nerve roots, and other neural structures. The aim is to support a structured understanding and evaluation of the tissue structure, lesion areas, and functional changes of the nervous system, improving data processing efficiency and result objectivity. It is widely applied to the automated analysis of neurosurgical images related to brain diseases, spinal diseases, peripheral nerve injuries, and neurotumors.

[0003] Traditional medical image analysis methods refer to the quantitative and qualitative analysis of neurological medical image data through rule-based image processing programs or image recognition processes designed based on human experience. These methods rely on image grayscale distribution statistics, edge detection algorithms, region growing algorithms, template matching, or morphological processing to extract and identify targets such as brain tissue, spinal cord structures, tumor boundaries, hemorrhage areas, nerve roots, or blood vessels. Some methods incorporate anatomical knowledge of the brain and spine, setting spatial constraints and using symmetry and regional partitioning information to assist in determining target regions and expressing their features.

[0004] Traditional methods rely on manually designed grayscale statistics and morphological processing rules, which lack self-matching noise reduction capabilities for non-uniform anisotropic signals in diffusion tensor imaging, resulting in artifacts and breaks in fiber tract tracking results. Spatial constraints based on anatomical knowledge only consider static partition information and cannot resolve dynamic structural displacements caused by neurophysiological activities. Region growing algorithms ignore the temporal variation characteristics of functional magnetic resonance imaging when registering gyral boundaries, causing spatial misalignment in structure-function cross-modal mapping. Template matching methods are not sensitive enough to individualized brain region variations. Manual feature extraction processes are difficult to capture the dynamic reconstruction characteristics of neural networks under stimulus conditions, limiting the ability to quantitatively assess functional compensation mechanisms. Summary of the Invention

[0005] To address the technical problems of traditional methods relying on manually designed grayscale statistics and morphological processing rules, lacking self-matching noise reduction capabilities for non-uniform anisotropic signals in diffusion tensor imaging, leading to artifacts and breaks in fiber tract tracking results, and the inability to analyze dynamic structural displacements induced by neurophysiological activity due to static partitioning information in anatomical knowledge-based spatial constraints, region growing algorithms ignoring temporal variations in functional magnetic resonance imaging (fMRI) during gyral boundary registration, resulting in spatial misalignment in cross-modal structure-function mapping, insufficient sensitivity of template matching to individualized brain region variations, and the difficulty of capturing the dynamic reconstruction characteristics of neural networks under stimulus conditions in manual feature extraction processes, thus limiting the quantitative assessment of functional compensation mechanisms, this invention provides a medical image analysis method and system based on image processing. The technical solution is as follows:

[0006] On the one hand, a medical image analysis method based on image processing is provided, which includes:

[0007] S1: Acquire the original diffusion tensor image sequence through magnetic resonance equipment, call the anisotropic diffusion filtering algorithm to perform noise suppression operation on the original image, extract the tensor direction change value to construct the structure vector set, and generate the preprocessed image sequence;

[0008] S2: Based on the preprocessed image sequence, perform three-dimensional normalization processing on the main diffusion direction value, use wavelet transform to decompose the sequence and extract the frequency band energy value, mark the voxels of energy exceeding the dynamic threshold, and output the fiber node coordinate set;

[0009] S3: Based on the fiber node coordinate set, perform point-by-point Euclidean distance offset calculation on the node coordinates in the resting and stimulated states, and perform Pearson correlation analysis on the distance offset and functional magnetic resonance time series data to generate a node displacement-functional response correlation map.

[0010] S4: Based on the node displacement-functional response correlation map, perform three-dimensional spatial interpolation on the whole brain nodes, align the interpolation results with the coordinates of the brain gyri boundary points in the structural magnetic resonance image, and output an individualized structure-function coupling strength distribution map;

[0011] S5: Based on the individualized structure-function coupling strength distribution map, a support vector machine classifier is used to perform classification operations, compare with clinical diagnostic standard parameters, and output an abnormal probability distribution map.

[0012] As a further aspect of the present invention, the preprocessed image sequence includes a signal-to-noise ratio parameter, an anisotropic fractional matrix, and a tensor eigenvalue distribution; the fiber node coordinate set specifically includes a high-frequency energy coordinate set, a low-frequency energy coordinate set, and an energy gradient coordinate set; the node displacement-functional response correlation map includes a displacement correlation coefficient matrix, a functional activation threshold distribution, and a spatiotemporal coupling matrix; the individualized structure-functional coupling strength distribution map specifically refers to the sulcus-gyrus coupling coefficient, gray-white matter interface intensity value, and cortical thickness mapping value; and the abnormality probability distribution map includes a classification confidence distribution, a regional abnormality index matrix, and a pathological feature weight map.

[0013] As a further aspect of the present invention, the specific steps of S1 include:

[0014] S101: Acquire the original diffusion tensor image sequence through magnetic resonance equipment, extract the tensor direction feature values ​​and corresponding gray-level distribution matrices at multiple locations in the image sequence, establish the diffusion direction distribution data structure, and generate diffusion direction change values;

[0015] S102: Call the diffusion direction change value, input it into the anisotropic diffusion filtering algorithm, iteratively calculate the grayscale correction value to suppress noise, and perform boundary judgment based on the tensor direction continuity to restrict pixel smoothing operation across the boundary and generate diffusion direction filtering weight value.

[0016] S103: Adjust the grayscale matrix values ​​of multiple frames in the original image sequence according to the diffusion direction filter weight value, reconstruct the tensor structure and output it uniformly to generate a preprocessed image sequence.

[0017] As a further aspect of the present invention, the specific steps of S2 include:

[0018] S201: Based on the preprocessed image sequence, extract the main diffusion direction value, call the three-dimensional vector normalization method, calculate the vector magnitude of the main diffusion direction, process the components to a uniform scale, and generate normalized direction vector values.

[0019] The principal diffusion direction value refers to the eigenvector corresponding to the principal eigenvalue in the diffusion tensor;

[0020] S202: Based on the normalized direction vector value, the vector value sequence is decomposed by wavelet transform, frequency band energy value data is extracted, voxels with energy values ​​greater than the dynamic threshold are selected, and a voxel sequence exceeding the threshold is obtained.

[0021] The wavelet transform performs multi-scale decomposition processing on the three-dimensional image sequence to obtain the feature information of the image under different frequencies and scales.

[0022] The dynamic threshold is adjusted based on statistical features and is set to the mean plus twice the standard deviation, referencing the mean and standard deviation of voxel energy from 200 clinical samples.

[0023] S203: Based on the above-threshold voxel sequence, locate the spatial position of the corresponding voxel in the three-dimensional coordinates of the image, summarize the spatial coordinate data, and generate a fiber node coordinate set.

[0024] As a further aspect of the present invention, the specific steps of S3 include:

[0025] S301: Based on the fiber node coordinate set, calculate the Euclidean distance between node pairs in resting and stimulated states, compare the differences between node pairs and construct a spatial offset sequence to obtain the node distance offset sequence.

[0026] S302: Call the node distance offset sequence, match the node distance offset with the time period of the functional magnetic resonance time series signal, filter the valid data, calculate the correspondence between the two sets of data according to the Pearson formula, and obtain the node correlation coefficient set.

[0027] The Pearson formula measures the strength of the linear relationship between nodal displacement and its corresponding functional magnetic resonance time-series signal, revealing whether structural changes are accompanied by functional changes.

[0028] S303: Based on the set of node correlation coefficients, pair the original coordinates with the correlation coefficients, construct a color level mapping relationship, and combine the coordinates to draw the color value spectrum corresponding to the node, generating a node displacement-functional response correlation spectrum.

[0029] As a further aspect of the present invention, the specific steps of S4 include:

[0030] S401: Based on the node displacement-functional response correlation map, select node coordinates and response values, calculate the weight coefficients of the response values ​​in three-dimensional space, construct an interpolation model using a weighted average method to calculate the grid values, and generate a three-dimensional functional response interpolation value distribution.

[0031] S402: Call the three-dimensional function response interpolation value distribution, filter out interpolation points in the structural magnetic resonance image template, match boundary point coordinates based on Euclidean distance, map interpolation values ​​to the gyri boundary, and obtain the gyri boundary matching value dataset;

[0032] S403: Based on the brain gyrus boundary matching value dataset, traverse the interpolation change points, calculate and accumulate the boundary point value differences, normalize the cumulative change to form the corresponding intensity, and generate an individualized structure-function coupling intensity distribution map.

[0033] As a further aspect of the present invention, the weighting coefficient is calculated using the following formula:

[0034]

[0035] Among them, w ijkΔx represents the weight coefficients of 3D mesh points i, j, k. i Δy represents the displacement of the i-th node along the x-axis, in mm. j ||d represents the displacement of the j-th node along the y-axis, in mm. k || represents the magnitude of the displacement vector at the k-th node, in mm, δ m This represents the absolute deviation of the m-th response value from the reference value, expressed in Pa. The gradient magnitude of the m-th response value is represented by r, in Pa / mm. ijk This represents the coordinates of the current grid point i, j, k, in mm. m Represents the coordinates of the m-th node, in mm, R max The maximum interpolation influence radius is represented in mm and is obtained based on a backtracking experiment of real-time MRI interpolation error, ranging from 10 to 20 mm. N represents the total number of nodes, and u represents the pressure-displacement conversion coefficient, calibrated by brain tissue elastic modulus experiments, in mm / Pa, with a value range of 0.02-0.05 mm / Pa.

[0036] As a further aspect of the present invention, the specific steps of S5 include:

[0037] S501: Based on the individualized structure-function coupling strength distribution map, divide the region to extract the coupling strength, match the spatial coordinates to classify the data, encode the region label to convert the vector format, unify the scale to form a vector set, and obtain the regional coupling feature coefficient set.

[0038] S502: Call the set of regional coupling feature coefficients as input to the support vector machine, calculate the classification interval according to the training set labels, perform interval optimization on the distance between the sample and the classification surface, calculate the support vector boundary coefficient values, and obtain the regional classification boundary coefficient matrix;

[0039] The support vector machine uses a radial basis function kernel, and the parameters are optimized in the training set through cross-validation.

[0040] S503: Based on the region classification boundary coefficient matrix, extract the set of coupling feature coefficients of the individuals to be distinguished, calculate the classification probability, compare the kernel function mapping offset between categories, and generate an anomaly probability distribution map.

[0041] As a further aspect of the present invention, the support vector boundary coefficient values ​​are calculated using the following formula:

[0042]

[0043] Where, α r,j n represents the support vector boundary coefficient value of the j-th class in the r-th region. rγ represents the number of support vectors in the r-th region. r,i w represents the coupling feature coefficient of the i-th sample in the r-th region. r,i b represents the support vector weight corresponding to the i-th sample in the r-th region. r Let m represent the classification bias term obtained from training in the r-th region. r λ represents the total number of dimensions involved in the support vector weight calculation in the r-th region. r The weight regularization adjustment parameter represents the r-th region.

[0044] On the other hand, an image processing-based medical image analysis system is provided, which is used to execute the above-described image processing-based medical image analysis method. The system includes:

[0045] The image preprocessing module is used to acquire the original diffusion tensor image sequence through the magnetic resonance device, call the anisotropic diffusion filtering algorithm to perform noise suppression operation on the original image, output the preprocessed image sequence, and pass it to the fiber node extraction module.

[0046] The fiber node extraction module is used to call the preprocessed image sequence, perform three-dimensional vector normalization on the main diffusion direction value, decompose the sequence using wavelet transform and extract the frequency band energy value, mark voxels whose energy exceeds the dynamic threshold, generate a fiber node coordinate set, and pass it to the correlation analysis module.

[0047] The correlation analysis module is used to perform point-by-point Euclidean distance offset calculation on the resting and stimulated state node coordinates based on the fiber node coordinate set, perform Pearson correlation analysis on the distance offset, generate a node displacement-functional response correlation map, and pass it to the coupling modeling module.

[0048] The coupling modeling module is used to call the node displacement-functional response correlation map, perform three-dimensional spatial interpolation operation on the whole brain nodes, spatially register the interpolation results with the coordinates of the brain gyri boundary points of the structural magnetic resonance image, generate a structural-functional coupling strength distribution map, and transmit it to the anomaly detection module.

[0049] The anomaly detection module is used to input the structure-function coupling strength distribution map into the support vector machine classifier, perform classification operations on the brain region feature vectors, compare them with the gray matter density threshold and white matter integrity index in the clinical diagnostic standard parameters, and output an anomaly probability distribution map.

[0050] The beneficial effects brought about by the technical solution provided by the embodiment of the present invention include at least:

[0051] By employing anisotropic diffusion filtering preprocessing of diffusion tensor image sequences, noise interference is reduced while preserving detailed fiber orientation features. The synergistic effect of three-dimensional vector normalization and wavelet transform frequency band energy analysis, along with a dynamic threshold screening mechanism, enhances the specificity of fiber node detection. The calculation of Euclidean distance offset of node displacement in resting and stimulated states, combined with functional magnetic resonance temporal correlation analysis, quantifies the spatiotemporal coupling relationship between neural activity and structural deformation. Three-dimensional spatial interpolation is used to fuse gyral boundary coordinates to establish an individualized structure-function coupling model, breaking through the rigid constraints of traditional fixed templates for anatomical partitioning. A support vector machine classifier identifies abnormal patterns based on coupling strength distribution features, achieving cross-scale correlation detection between microstructural distortions and macroscopic functional abnormalities, and improving the efficiency of imaging biomarker identification for early neurodegenerative diseases. Attached Figure Description

[0052] Figure 1 This is a schematic diagram of the workflow of the present invention;

[0053] Figure 2 This is a system flowchart of the present invention. Detailed Implementation

[0054] The technical solution of the present invention is described below in conjunction with the accompanying drawings.

[0055] In the embodiments of the present invention, words such as "exemplarily" and "for example" are used to indicate examples, illustrations, or explanations. Any embodiment or design described as an "exemplary" in the present invention should not be interpreted as being preferred or advantageous over other embodiments or designs. Rather, the use of the word "exemplary" is intended to present concepts in a concrete manner. Furthermore, in the embodiments of the present invention, "and / or" can mean both or either of the two.

[0056] In the embodiments of the present invention, the terms "image" and "picture" may be used interchangeably. It should be noted that, when the distinction between them is not emphasized, their intended meanings are the same. The terms "of," "corresponding," and "corresponding" may be used interchangeably. It should be noted that, when the distinction between them is not emphasized, their intended meanings are the same.

[0057] In the embodiments of the present invention, sometimes a subscript such as W1 may be written as a non-subscript such as W1. When the difference is not emphasized, the meanings to be expressed are the same.

[0058] In order to make the technical problems, technical solutions and advantages to be solved by the present invention clearer, a detailed description will be given below with reference to the accompanying drawings and specific embodiments.

[0059] Please see Figure 1This invention provides a medical image analysis method based on image processing. The processing flow of this method may include the following steps:

[0060] S1: Acquire the original diffusion tensor image sequence through magnetic resonance equipment, call the anisotropic diffusion filtering algorithm to perform noise suppression operation on the original image, extract the tensor direction change value to construct the structure vector set, and generate the preprocessed image sequence;

[0061] S2: Based on the preprocessed image sequence, perform three-dimensional normalization on the main diffusion direction value, use wavelet transform to decompose the sequence and extract the frequency band energy value, mark the voxels of energy exceeding the dynamic threshold, and output the fiber node coordinate set;

[0062] S3: Based on the fiber node coordinate set, perform point-by-point Euclidean distance offset calculation on the node coordinates in resting and stimulated states, and perform Pearson correlation analysis on the distance offset and functional magnetic resonance time series data to generate a node displacement-functional response correlation map.

[0063] S4: Based on the nodal displacement-functional response correlation map, perform three-dimensional spatial interpolation on the whole brain nodes, align the interpolation results with the coordinates of the gyri boundary points of the structural magnetic resonance image, and output an individualized structure-function coupling strength distribution map;

[0064] S5: Based on the individualized structure-function coupling strength distribution map, a support vector machine classifier is used to perform classification operations, compare with clinical diagnostic standard parameters, and output an abnormal probability distribution map.

[0065] The preprocessed image sequence includes signal-to-noise ratio parameters, anisotropy fraction matrix, tensor eigenvalue distribution, fiber node coordinate set specifically includes high-frequency energy coordinate set, low-frequency energy coordinate set, and energy gradient coordinate set, node displacement-functional response correlation map including displacement correlation coefficient matrix, functional activation threshold distribution, and spatiotemporal coupling matrix, individualized structure-functional coupling strength distribution map specifically refers to sulcus-gyrus coupling coefficient, gray-white matter interface intensity value, and cortical thickness mapping value, and abnormality probability distribution map including classification confidence distribution, regional abnormality index matrix, and pathological feature weight map.

[0066] Specifically, the steps of S1 are as follows:

[0067] S101: Acquire the original diffusion tensor image sequence through magnetic resonance equipment, extract the tensor direction feature values ​​and corresponding gray-level distribution matrices at multiple locations in the image sequence, establish the diffusion direction distribution data structure, and generate diffusion direction change values;

[0068] The change in diffusion direction measures the degree of directional change by comparing it with the principal direction vector of the diffusion tensor, reflecting the spatial continuity and directional consistency characteristics of the tissue microstructure.

[0069] To acquire the original diffusion tensor image sequence using a magnetic resonance imaging (MRI) device, the scanning parameters must first be set, such as echo time (TE) of 80 ms, repetition time (TR) of 7000 ms, and diffusion weight b value set to 1000 s / mm. 2 During scanning, multi-directional diffusion encoding is performed on the head region, such as using 30 diffusion gradient directions. An image is generated for each direction, forming a complete diffusion tensor image sequence. Subsequently, for each voxel in a selected region (such as the corpus callosum of the brain) in the image sequence, the diffusion signal under the differential diffusion direction is extracted. A 3×3 diffusion tensor matrix D is constructed for each voxel by fitting the Stejskal-Tanner formula. Then, the feature vector of the main diffusion direction is obtained through tensor eigenvalue decomposition. The eigenvalues ​​(λ1, λ2, λ3) reflect the directional diffusion intensity. Correspondingly, the gray values ​​of each voxel are collected to form a gray-level distribution matrix G. i,j,k Combined with grayscale change rate:

[0070]

[0071] A three-dimensional structured distributed data structure is established to record the spatial correlation between the principal directions of tensors and the gray-level distribution in multiple directions. After obtaining the directional distribution, the directional change value V is calculated on a voxel-by-voxel basis. Δ Its form is the cosine difference of the angle θ between the current voxel and the principal directions of the adjacent voxels, i.e., V Δ = 1 - cos(θ), where If the principal direction of a certain voxel is Its adjacent voxels are Then V Δ =1-(0.8×0.6+0.4×0.6+0.4×0.5)=1-(0.48+0.24+0.20)=1-0.92=0.08, the closer this directional change value is to 0, the more consistent the direction. If V Δ A value >0.2 indicates a region of drastic directional change. A three-dimensional directional change map is generated for each voxel. In clinical scenarios, such as the detection of white matter lesions, abnormal lesions exhibit drastic directional changes. Therefore, by using the V of continuous voxels... Δ The trend of value changes can further help identify lesion areas. For example, after collecting directional change values ​​in a 10×10×10 voxel area, the following statistics are generated:

[0072] Table 1: Statistics of Diffusion Direction Changes

[0073]

[0074]

[0075] As shown in Table 1, the angle between voxel A3 and the principal direction of its neighboring voxels varies considerably, with the direction change value V.Δ Reaching 0.11, combined with subsequent tensor orientation maps, suspicious lesion points can be located. Spatial visualization of this orientation change map forms a diffusion direction structure map, which is output via three-dimensional coordinate mapping (x, y, z, V). Δ The quadruple is used as the input for subsequent anisotropic diffusion weight allocation and filtering determination.

[0076] S102: Call the diffusion direction change value, input it into the anisotropic diffusion filter algorithm, iteratively calculate the grayscale correction value to suppress noise, and perform boundary judgment based on the continuity of tensor direction to restrict pixel smoothing operation across the boundary and generate diffusion direction filter weight value.

[0077] The diffusion direction filtering weight value refers to the directional smoothing coefficient that is self-matched and adjusted according to the local structural features of the image during the anisotropic diffusion process;

[0078] After obtaining the change value of diffusion direction, the V of each voxel must first be calculated. Δ Values ​​are read point by point, and initial classification is performed based on preset directional consistency intervals. Specifically, the interval [0, 0.1] is determined as a segment with highly consistent direction, the interval (0.1, 0.2) is determined as a segment with slight directional transition, and values ​​exceeding 0.2 are determined as regions with abrupt directional changes. Subsequently, the spatial distribution relationship of regions with abrupt directional changes is retrieved in the voxel and its neighborhood to obtain a local directional change structure map of the voxel. This structure map is used as a preliminary spatial mask and input into the grayscale smoothing stage. Then, the grayscale values ​​of the corresponding voxels in the original image sequence are subjected to initial differential processing, and the difference value ΔG = |G| is calculated. i,j,k -G i:1,j,k |+|G i,j,k -G i,j:1,k |+|G i,j,k -G i,j,k:1 |, forming the initial gradient tensor matrix T g , for T g Gradient magnitudes are sorted, and the top 30% of the largest gradient values ​​are used to create noise enhancement masks. These masks are then subjected to a cross-logic AND operation with the aforementioned directional change mask to extract regions with drastic directional changes and significant grayscale abrupt changes. These regions serve as the boundary for diffusion filtering. Next, directional weighted smoothing is applied to the voxel grayscale values ​​outside these regions. The calculation method involves assigning smoothing weights to the grayscale values ​​in each directional neighborhood according to the direction cosine relationship. For example, if the current voxel's principal direction vector is... Its direction vector with a certain neighboring voxel is Then the smoothing weight coefficient w of the neighborhood voxel n for In the neighborhood, this value is normalized to a probability-based weighting factor. After that Using the weights, a weighted average is calculated for the gray levels of the neighborhood. Suppose there are 5 neighborhood voxels with gray levels G1 = 56, G2 = 59, G3 = 53, G4 = 60, and G5 = 54, and corresponding cosine weighting coefficients [0.90, 0.92, 0.87, 0.89, 0.91]. Then, after normalization... Direction weighted correction value Replace the original voxel grayscale with this correction value to complete one anisotropic diffusion iteration operation. To avoid cross-region smoothing, all values ​​V that belong to directional changes are excluded. Δ For regions >0.2, the smoothing operation is interrupted and the grayscale is not updated. The overall iteration process uses N=5 rounds. When the grayscale change value ΔG of a certain voxel is <1, it is considered to be converged, thereby obtaining the diffusion direction filtering weight values ​​of all voxels.

[0079] S103: Adjust the grayscale matrix values ​​of multiple frames in the original image sequence according to the diffusion direction filter weight values, reconstruct the tensor structure and output it uniformly to generate a preprocessed image sequence;

[0080] After completing the anisotropic diffusion filtering operation, the voxel gray values ​​in the original image sequence are uniformly corrected based on the calculated directional filtering weight values. First, the gray response vector corresponding to each voxel under the differentiated diffusion direction needs to be extracted. Combined with its main diffusion direction vector With the neighborhood direction vector, through the cosine-weighted coefficient vector Perform element-wise weighted correction, specifically correcting the grayscale vector. That is, the gray value in each diffusion direction is adjusted according to its directional smoothing coefficient. For example, if the voxel has a gray value of 56 in direction 1 with a weight of 0.9, and a gray value of 58 in direction 2 with a weight of 0.85, then the gray value correction values ​​for the first two items are 56 × 0.9 = 50.4 and 58 × 0.85 = 49.3, respectively. After the directional correction is completed, a new gray value vector is formed. Then, based on the corrected grayscale vector, the Stejskal-Tanner model is refitted to construct a new diffusion tensor matrix D. ′ The process involves calculating the updated principal diffusion direction and eigenvalue vector through tensor reconstruction to update the structure of each voxel. Finally, the voxel-corrected grayscale image frames are combined into a complete tensor image sequence, which is then output as a preprocessed image dataset. The weight adjustment step in this process uses the directional consistency value V... Δ Define the baseline segmentation if V Δ When <0.1, directional weight The value range is set to [0.85, 1.0]; 0.1 <V Δ When V ≤ 0.2, the value is set to [0.6, 0.85]. ΔIf W > 0.2, then W d If the value is fixed at 0 and does not participate in the weighting, it ensures that no data interference from regions with abrupt changes in direction is introduced during the reconstruction of the structural image, thus completing the tensor structure reconstruction and outputting a unified image sequence.

[0081] Specifically, the steps of S2 are as follows:

[0082] S201: Based on the preprocessed image sequence, extract the main diffusion direction value, call the three-dimensional vector normalization method, calculate the vector magnitude of the main diffusion direction, process the components to a uniform scale, and generate normalized direction vector values.

[0083] The principal diffusion direction value refers to the eigenvector corresponding to the principal eigenvalue in the diffusion tensor;

[0084] Based on the preprocessed image sequence, the main diffusion direction value is extracted. First, the preprocessed tensor image sequence is read. For each voxel position in the sequence, its corresponding diffusion tensor matrix D′ is called sequentially. This matrix consists of six independent components, namely Di, Dj, ... xx D yy D zz D xy D xz D yz Taking voxel number T001 as an example, its tensor matrix is ​​initially extracted as [D]. xx =0.0011, D yy =0.0009, D zz =0.0013, D xy =0.0002, D xz =0.0001, D yz =0.0003]mm 2 / s, then perform eigenvalue decomposition on the extracted tensor matrix, transforming the tensor matrix into its corresponding eigenvalue vector [λ1, λ2, λ3] and eigenvector matrix. Taking voxel T001 as an example, λ1 = 0.0015 mm was calculated. 2 / s, λ2=0.0010mm 2 / s, λ3=0.0008mm 2 / s, which corresponds to the eigenvector of the main diffusion direction. Subsequently, the feature vector of the main diffusion direction was analyzed. Calculate its vector magnitude L v The calculation formula is:

[0085]

[0086] The input value is:

[0087]

[0088] Next, a normalization operation is performed, calculating the normalized value of each component. The normalization formula is as follows:

[0089]

[0090] Substituting the values, we get:

[0091]

[0092] Thus, the normalized principal diffusion direction vector of voxel T001 is obtained as [0.7510, 0.5007, 0.4306]. This normalized direction vector is unified to the standard modulus scale, which is beneficial for subsequent sequence analysis. After completing the voxel normalization, the voxels in the entire image sequence are traversed sequentially. Assuming that the image sequence contains a total of 100,000 voxels numbered from T001 to T100,000, the above feature decomposition and normalization process is called one by one in all voxels to form the normalized direction vector set of the entire image sequence. For example, the following normalized direction vectors are obtained for voxels T002 to T005 respectively:

[0093] Table 2: Normalized Direction Vector Table

[0094]

[0095] As shown in Table 2, the normalized main diffusion direction vector values ​​of the differential voxels are obtained by calculation, and then the set of normalized direction vectors is used as the basic data input for subsequent directional consistency analysis and spatial energy decomposition processing.

[0096] S202: Based on the normalized direction vector value, the vector value sequence is decomposed by wavelet transform, frequency band energy value data is extracted, voxels with energy values ​​greater than the dynamic threshold are selected, and the super-threshold voxel sequence is obtained.

[0097] Wavelet transform uses the db4 wavelet basis function to perform a 3-level decomposition, and performs multi-scale decomposition processing on the three-dimensional image sequence to obtain the feature information of the image under different frequencies and scales.

[0098] The dynamic threshold is adjusted based on statistical characteristics. It is set to the mean plus twice the standard deviation, referencing the mean and standard deviation of voxel energy from 200 clinical samples.

[0099] According to the normalized direction vector values, first perform three-dimensional wavelet transform processing on the normalized direction vector sequence of each voxel. Specifically, the vector components of the voxel in the x, y, and z directions are used as independent signals for input, and the Daubechies4 wavelet basis function is used for 3-layer decomposition. Taking the normalized direction vector [0.7510, 0.5007, 0.4306] of voxel T001 as an example, the approximation coefficient cA1 = 0.72 and the detail coefficient cD1 = 0.031 are obtained after wavelet decomposition of its x-direction component 0.7510. The cA1 = 0.48 and cD1 = 0.0207 are obtained after decomposition of the y-direction component 0.5007, and the cA1 = 0.41 and cD1 = 0.0206 are obtained after decomposition of the z-direction component 0.4306. Subsequently, energy calculation is performed on the third-layer detail coefficients, and the energy value formula is where n is the number of detail coefficients. Taking the third-layer detail coefficients in the x direction of voxel T001 as an example, its coefficient sequence is [0.008, 0.005, -0.003], then the energy value calculation is E x = 0.008 2 + 0.005 2 + (-0.003) 2 = 0.000064 + 0.000025 + 0.000009 = 0.000098. Similarly, calculate the energy E in the y direction y = 0.000035, and the energy E in the z direction z = 0.000028. The total energy E total = 0.000098 + 0.000035 + 0.000028 = 0.000161. The dynamic threshold is set based on the statistics of 200 clinical samples. Assuming the sample energy mean μ = 0.00015 and the standard deviation σ = 0.00002, then the dynamic threshold T = μ + 2σ = 0.00015 + 2×0.00002 = 0.00019. When the E of voxel T001 total = 0.000161 < T, it is screened out, while the normalized direction vector of voxel T003 is [0.8123, 0.3822, 0.4389], and its total energy calculation is E total = 0.00025 > T, so it is retained. After traversing the voxels, a sequence of super-threshold voxels is formed.

[0100] S203: Based on the sequence of super-threshold voxels, locate the spatial positions of the corresponding voxels in the three-dimensional coordinates of the image, summarize the spatial coordinate data, and generate a fiber node coordinate set;

[0101] Based on the super-threshold voxel sequence, the three-dimensional coordinate index of each super-threshold voxel is first extracted from the image sequence metadata. Taking voxel T003 as an example, its original voxel index is (128, 96, 64). By reading the voxel spacing parameters (Δx = 1.0 mm, Δy = 1.0 mm, Δz = 2.0 mm) from the image header file, a physical space coordinate transformation is performed. The calculation method is x = i × Δx, y = j × Δy, z = k × Δz. Substituting the values, we get x = 128 × 1.0 = 128.0 mm, y = 96 × 1.0 = 96.0 mm, z = 64 × 2.0 = 128.0 mm, generating the coordinates. The coordinates (128.0, 96.0, 128.0) are used. The same calculation is performed on the voxel T004 index (130, 98, 66), resulting in (130.0, 98.0, 132.0). The voxel T005 index (135, 102, 70) is converted to (135.0, 102.0, 140.0). Then, spatial continuity is checked on the coordinate data, with the rule that the distance between adjacent voxels should not exceed 3mm. If the coordinates of voxel T006 are (136, 103, 71), and after conversion to (136.0, 103.0, 142.0), its Euclidean distance to the preceding voxel T005 is calculated as follows:

[0102]

[0103] To ensure continuity, if the distance between a voxel's coordinates and the previous node exceeds 3mm, interpolation compensation is triggered. The interpolated intermediate point coordinates are linear interpolations of adjacent coordinates. For example, voxel T007 coordinates (140, 105, 75) are transformed to (140.0, 105.0, 150.0), and its distance from the previous voxel T006 is:

[0104]

[0105] If the threshold is exceeded, two intermediate points are inserted between T006 and T007. The interpolation formula is x′=x0+t×(x1-x0), where t=0.33 and t=0.66. The calculated intermediate point 1 is:

[0106] (136.0 + 0.33 × 4.0, 103.0 + 0.33 × 2.0, 142.0 + 0.33 × 8.0) = (137.32, 103.66, 144.64);

[0107] Midpoint 2 is:

[0108] (136.0+0.66×4.0,103.0+0.66×2.0,142.0+0.66×8.0)=(138.64,104.32,147.28);

[0109] Finally, the original coordinates and interpolated coordinates are arranged in spatial order to form a continuous trajectory point sequence. The node coordinates are sorted by anatomical orientation and output as a fiber node coordinate set. The data format is [(128.0, 96.0, 128.0), (130.0, 98.0, 132.0), (135.0, 102.0, 140.0), (136.0, 103.0, 142.0), (137.32, 103.66, 144.64), (138.64, 104.32, 147.28), (140.0, 105.0, 150.0), ...], thus completing the reconstruction of the spatial topology of the fiber bundle.

[0110] Specifically, the steps of S3 are as follows:

[0111] S301: Based on the fiber node coordinate set, calculate the Euclidean distance between node pairs in resting and stimulated states, compare the differences between node pairs and construct a spatial offset sequence to obtain the distance offset sequence between nodes.

[0112] Based on the fiber node coordinate set, the node coordinates are first sorted according to anatomical orientation, following the order of left to right, front to back, and top to bottom. The three-dimensional coordinates are then sorted in ascending order along the x, y, and z axes, i.e., x-values ​​are compared first; if x-values ​​are the same, y-values ​​are compared; if y-values ​​are also the same, z-values ​​are compared. This results in a sorted list of nodes. For example, if node A is (128.0, 96.0, 128.0) and node B is (130.0, 98.0, 132.0), then A is listed before B. After sorting, for each pair of adjacent nodes in the sorted list, the Euclidean distance difference between the resting and stimulated states is calculated. The Euclidean distance formula is:

[0113]

[0114] Where (x) i ,y i , z i ), (x j ,y j , z j Let be the coordinates of node i and node j, respectively. Substitute these coordinates into the above formula to calculate the distance between node pairs for the resting state coordinate set and the stimulus state coordinate set, respectively. and This leads to the offset:

[0115]

[0116] Taking some nodes as an example, the distance between nodes A and B in their resting state is:

[0117]

[0118] Let A be the stimulus state. ′ With B ′ The coordinates are (128.5, 96.5, 128.5) and (130.2, 98.3, 132.1), then:

[0119]

[0120] The distance offset of the node relative to AB is:

[0121] Δd AB =|4.37-4.90|=0.53;

[0122] These distance offset values ​​are calculated sequentially for adjacent node pairs, and the results are arranged in order to form a spatial offset sequence. To ensure outlier removal, if the offset value of a node pair exceeds a set threshold of 10.0 mm, it is considered a positioning or registration error and is removed. This threshold is defined by the overall node offset mean and twice the standard deviation. For example, if the average offset in the sampling is 0.75 and the standard deviation is 2.5, then the maximum allowable offset threshold is set to 0.75 + 2 × 2.5 = 5.75, rounded up to 10.0 mm. The table below lists the distance offset results for some node pairs:

[0123] Table 3: Node Pair Distance Offset Table (Unit: mm)

[0124]

[0125] As shown in Table 3, the distance offset sequence is obtained by calculating the distance between node pairs in the resting and stimulated states and subtracting them.

[0126] S302: Call the node distance offset sequence, match the node distance offset with the time period of the functional magnetic resonance time series signal, filter the valid data, calculate the correspondence between the two sets of data according to the Pearson formula, and obtain the set of node correlation coefficients.

[0127] Pearson's formula measures the strength of the linear relationship between nodal displacement and its corresponding functional magnetic resonance time series signal, revealing whether structural changes are accompanied by functional changes.

[0128] The above-mentioned inter-node distance offset sequence is used to compare with the time-series signal data recorded by functional magnetic resonance imaging (fMRI). To ensure the effectiveness of the comparison, the node offset sequence and the time-series signal must first be aligned in the time dimension. Specifically, they are aligned by the timestamp of each sampling point. If there is a difference in the sampling frequency between the two, the node offset sequence time is interpolated to the fMRI time for matching. The estimated value of the node offset at each fMRI sampling point is calculated using linear interpolation. After matching, invalid data points are removed, i.e., time points with NaN or signal amplitude less than the set baseline value (such as 0.01% signal change). The remaining time period data are then correlated according to the Pearson correlation formula.

[0129]

[0130] Where x i Let y be the distance offset between nodes at time i. i Let x and y be the mean values ​​of the fMRI signal at the corresponding time point, respectively. The sampling sequence is set to a sample group of length 10. The node offset sequence is [0.53, 0.52, 0.21, 0.33, 0.40, 0.39, 0.31, 0.35, 0.45, 0.50], and the fMRI signal sequence is [0.004, 0.003, 0.001, 0.002, 0.003, 0.004, 0.003, 0.002, 0.004, 0.005]. Calculate their mean values ​​respectively.

[0131]

[0132] Calculate the numerator:

[0133] ∑(x i -x)(y i -y)=0.0008863;

[0134] Calculate the denominator:

[0135]

[0136] The correlation coefficients show that there is a clear linear correspondence between the structural offset of the current node pair and the time variation of the functional signal, and finally a set of correlation coefficients corresponding to each pair of nodes is obtained.

[0137] S303: Based on the set of node correlation coefficients, pair the original coordinates with the correlation coefficients, construct a color level mapping relationship, and combine the coordinates to draw the color value spectrum corresponding to the node, generating a node displacement-functional response correlation spectrum;

[0138] Based on the aforementioned set of node correlation coefficients, the coordinate values ​​corresponding to each pair of nodes are paired with their correlation coefficient values. A color mapping scheme is set, with a mapping range of -1 to 1, where -1 maps to dark blue, 0 maps to white, and 1 maps to dark red. The average position of each node or node pair is used as the drawing reference point. The coordinates of node pair AB are (128, 96, 128) - (130, 98, 132). Therefore, the drawing center point is:

[0139]

[0140] If the correlation coefficient of the node pair is 0.8, the corresponding color value is medium red. The center points of the node pairs are drawn in sequence and filled with the corresponding color values ​​in this way to form a three-dimensional point cloud color chart. The color points are superimposed in the original sorting structure under a unified coordinate system to construct the node displacement-functional response correlation map.

[0141] Specifically, the steps of S4 are as follows:

[0142] S401: Based on the nodal displacement-functional response correlation map, select nodal coordinates and response values, calculate the weight coefficients of the response values ​​in three-dimensional space, construct an interpolation model using a weighted average method to calculate the mesh values, and generate a three-dimensional functional response interpolation value distribution.

[0143] Based on the nodal displacement-functional response correlation map, the coordinate set of the three-dimensional node distribution is first obtained. Displacement information of each node under real-time loading conditions is collected, and the corresponding functional response value of each node is extracted, such as the local pressure response of brain tissue under specific stimuli. The node coordinates are recorded in Cartesian coordinates in millimeters, and the functional response value is expressed as pressure (Pa). Subsequently, for each node, the displacement changes along the x, y, and z axes are recorded, forming a displacement vector (Δx) in three-dimensional space. i Δy j d k Given the node coordinates and their response values, calculate the influence of each node on the surrounding grid points, i.e., its weight coefficient w. ijk In the calculation, assuming a functional map of a brain tissue slice, five sampling nodes are defined on it. These nodes are located at three-dimensional coordinates (10, 10, 10), (12, 11, 10), (15, 14, 12), (16, 15, 13), and (20, 18, 15), with corresponding response values ​​of 1500 Pa, 1700 Pa, 1800 Pa, 2000 Pa, and 2200 Pa, respectively. Within the space near the current computational grid point (13, 12, 11), the weight coefficients of these five points are calculated, and their weight values ​​are averaged with the response values ​​to obtain the interpolated value for that grid point. Here, Δx... i with Δy jLet |d| represent the displacement difference between corresponding nodes of the grid points in the x and y directions, respectively. k || represents the nodal displacement modulus, calculated from the length of the three-dimensional displacement vector. For example, if the displacement of a node in the three axes is (1.2, 0.8, 1.5) mm, then its modulus is... absolute deviation of the node's response value δ m The difference between the node value and the reference value is taken. If the reference value is set to 1600 Pa, then the deviation of the first node is |1500-1600|=100 Pa, and the gradient magnitude of the response value is... This can be calculated from the difference in response values ​​and the distance between adjacent nodes. For example, if the difference in response values ​​between node 1 and node 2 is 200 Pa, then...

[0144] The Euclidean distance between the points is but Then, the parameters in the weighting formula are substituted into the calculation:

[0145]

[0146] Among them, w ijk Δx represents the weight coefficients of 3D mesh points i, j, k. i Δy represents the displacement of the i-th node along the x-axis, in mm. j ||d represents the displacement of the j-th node along the y-axis, in mm. k || represents the magnitude of the displacement vector at the k-th node, in mm, δ m This represents the absolute deviation of the m-th response value from the reference value, expressed in Pa. The gradient magnitude of the m-th response value is represented by r, in Pa / mm. ijk This represents the coordinates of the current grid point i, j, k, in mm. m Represents the coordinates of the m-th node, in mm, R max The maximum interpolation influence radius is represented in mm and is obtained based on a backtracking experiment of real-time MRI interpolation error, ranging from 10 to 20 mm. N represents the total number of nodes, and u represents the pressure-displacement conversion coefficient, calibrated by brain tissue elastic modulus experiments, in mm / Pa, with a value range of 0.02-0.05 mm / Pa.

[0147] Compare grid point (13, 12, 11) with the first node (10, 10, 10): Δx i =13-10=3, Δy j =12-10=2,|Δx i ·Δy j |=6,||d k ||=2.08mm;

[0148] Molecular part:

[0149] δ1 = |1500-1600| = 100 Pa;

[0150]

[0151] The denominator term corresponds to the first node:

[0152] Assuming that after calculating for all 5 nodes, ∑≈0.0332+0.0285+0.0257+0.0231+0.0210=0.1315;

[0153] Distance Term

[0154] The distance ratio term (1 - 3.74 / 15) = 0.7507;

[0155] final

[0156] Similarly, calculate w for the 5 nodes. ijk After being normalized proportionally, the values ​​are used as weighting coefficients, which are then multiplied by the corresponding response values ​​to complete the interpolation.

[0157] Table 4: Three-dimensional nodal displacement and response parameters:

[0158]

[0159] Table 4 shows the displacement vector, response value, deviation, and gradient information recorded by the five nodes in three-dimensional space.

[0160] The results show that the interpolated values ​​of the functional response of the current three-dimensional grid points are calculated by the weighted average of the surrounding nodes, and the influence of the differentiated nodes is determined by the parameters. The above interpolation model can be used to continuously estimate the functional response at any location in space and form a three-dimensional functional response interpolation value distribution.

[0161] S402: Call the three-dimensional function response interpolation value distribution, filter out interpolation points in the structural magnetic resonance image template, match boundary point coordinates based on Euclidean distance, map interpolation values ​​to the gyri boundary, and obtain the gyri boundary matching value dataset;

[0162] After invoking the 3D functional response interpolation value distribution, all interpolation grid points within the structural magnetic resonance image template are first extracted. The spatial coordinates of each interpolation point are obtained, recorded in millimeters as ternary Cartesian coordinates. Simultaneously, the coordinate set of gyral boundary points obtained from image segmentation in the cerebral cortex region is located. Then, for each interpolation point, a distance calculation operation is performed between it and the gyral boundary point using Euclidean distance calculation. It is determined whether the minimum distance between the interpolation point and the gyral boundary point is less than a set mapping threshold. If the condition is met, it is determined that the interpolation point can be mapped to the boundary point. The mapping method is to assign the functional response value of the interpolation point to the nearest boundary point, completing one projection of the interpolation value to the gyral boundary. The above process is repeated to traverse all interpolation points and gyral boundary points, establishing a mapping table between interpolation values ​​and gyral boundary points. The calculation of Euclidean distance follows a 3D spatial formula. For example, if the coordinates of an interpolation point are (13.0, 12.0, 11.0) mm and the coordinates of a boundary point are (12.0, 12.5, 10.5) mm, then their Euclidean distance is... If the set mapping threshold is 2mm, the point is successfully mapped, and the interpolated value is assigned to the boundary point. If another interpolated point (18, 20, 14) is more than 2mm away from any boundary point, then the interpolated point does not participate in the mapping. The successfully mapped interpolated points form the set of response values ​​of the gyral boundary, which constitutes the gyral boundary matching value dataset. Each item in this dataset consists of the spatial coordinates of the gyral boundary point and its mapped response value.

[0163] S403: Based on the brain gyrus boundary matching value dataset, traverse the interpolation change points, calculate and accumulate the boundary point value difference, normalize the cumulative change to form the corresponding intensity, and generate an individualized structure-function coupling intensity distribution map.

[0164] Based on the constructed brain gyrus boundary matching value dataset, the interpolated response values ​​in the boundary point coordinate set are first traversed. For each point, response value changes in its local spatial region are detected. The detection method involves constructing a local spherical neighborhood (e.g., enclosing a region with a radius of 3 mm) centered on each boundary point. The remaining boundary points within this neighborhood are then selected, and their response values ​​are extracted. The difference between these points and the response value of the center point is calculated, and the absolute value is taken to obtain the set of response variations within the local region. This set is then summed to obtain the total local response variation value. This value is then normalized to the original response value of the point. Specifically, the normalization method is to divide the total response variation value by the original response value of the point, representing the intensity of the corresponding response change. This process is repeated to traverse the boundary points, forming a set of intensity values ​​corresponding to each boundary point. For example, if the response value of a brain gyrus boundary point is 2000 Pa, and the response values ​​of the other five points in the neighborhood are 1950 Pa, 1980 Pa, 2100 Pa, 2020 Pa, and 1900 Pa, then the corresponding response variation value is:

[0165] |2000-1950|+|2000-1980|+|2000-2100|+|2000-2020|+|2000-1900|=50+20+100+20+100=290Pa;

[0166] The total response change is 290 Pa, and the corresponding intensity value is 290 / 2000 = 0.145. That is, the corresponding intensity of the boundary point is 0.145. After generating the corresponding intensity value of the boundary point in this way, the value can be combined with the spatial location of the boundary point to construct a three-dimensional spatial intensity distribution map, thereby completing the generation of the individualized structure-function coupling intensity map.

[0167] Specifically, the steps of S5 are as follows:

[0168] S501: Based on the individualized structure-function coupling strength distribution map, the coupling strength is extracted by dividing the region, the spatial coordinates are matched and the data is classified, the region labels are encoded and converted into vector format, the vector set is formed by unifying the scale, and the set of regional coupling feature coefficients is obtained.

[0169] Based on the individualized structure-function coupling intensity distribution map, this step first requires acquiring structural images (such as T1-weighted MRI) and functional images (such as fMRI) of brain regions. By comparing the consistency between the functional connectivity intensity in space and the structural connectivity map, the coupling intensity value corresponding to each voxel in the brain region is calculated. For example, in a standard space of 10mm×10mm×10mm, the brain map is divided into 1000 voxel regions. Pearson correlation calculations of structural and functional signals are performed on each region to obtain the individual coupling intensity map. Then, a region division operation is performed on this distribution map to divide the entire brain map into predefined functional modules (such as 90 regions in the AAL template). The mean coupling intensity of voxels in each region is extracted as the representative coupling value of that region. When matching spatial coordinates for data classification, a unified mind map segmentation template needs to be used to align individual mind maps to the MNI standard space. Regions are located based on this spatial coordinate system, and region label codes are recorded (e.g., 01 for the prefrontal cortex, 02 for the parietal cortex, etc.). Then, one-hot encoding is performed on the region numbers; for example, region number 02 is converted to a vector format of [0, 1, 0, ..., 0]. When forming a vector set with a unified scale, the region vectors are Z-score standardized to have a distribution characteristic of 0 mean and 1 variance, facilitating equal-weight comparisons in the subsequent classification model. Finally, the region number vectors and coupling strength values ​​are combined into a vector of length n+1, where n is the dimension of the region label vector. For example, if 90 regions are used, each vector is 91-dimensional, forming the entire set of region coupling feature coefficients for the individual. For example, the coupling strength of a subject's prefrontal region is 0.65, the region code is 01, and the one-hot vector is [1, 0, 0, ..., 0]. After standardization, the coupling strength is (0.65-μ) / σ. If μ = 0.5 and σ = 0.1, then the coupling strength after transformation is 1.5, and the final vector of this region is [1, 0, 0, ..., 0, 1.5]. Finally, the set of region coupling feature coefficients is constructed.

[0170] S502: Call the set of regional coupling feature coefficients as input to the support vector machine, calculate the classification margin according to the training set labels, perform margin optimization on the distance between samples and classification surfaces, calculate the support vector boundary coefficient values, and obtain the regional classification boundary coefficient matrix;

[0171] Support Vector Machines use radial basis function kernels, and the parameters are optimized in the training set through cross-validation.

[0172] When using the set of region-coupled feature coefficients as input to the support vector machine, the region vectors of each individual need to be reorganized. Features of individuals within the same region are grouped into a sample matrix, with each row representing a single sample and each column representing the vector feature values ​​corresponding to the same region. When calculating the classification margin for the training set labels, each sample is compared to its corresponding category (e.g., patient vs. healthy). The classification boundary is maximized for each group of support vectors. During margin optimization, samples are projected into a high-dimensional kernel function mapping space, the distance difference between the sample and the classification surface is calculated, and the direction and position of the normal vector of the classification surface are adjusted. The weight values ​​w of the support vectors are then gradually adjusted. r,i , and the bias value b r This maximizes the margin between the support vector points and the separating hyperplane; subsequently, the support vector boundary coefficients are calculated using the following formula:

[0173]

[0174] Where, α r,j n represents the support vector boundary coefficient value of the j-th class in the r-th region. r γ represents the number of support vectors in the r-th region. r,i w represents the coupling feature coefficient of the i-th sample in the r-th region. r,i b represents the support vector weight corresponding to the i-th sample in the r-th region, and its value is generally set between -1 and 1. r Let m represent the classification bias term obtained from training in the r-th region. r λ represents the total number of dimensions involved in the support vector weight calculation in the r-th region. r The weight regularization adjustment parameter represents the r-th region, typically with a value of 0.01.

[0175] Taking region 5 as an example, there are n support vectors. r =4, the coupled feature value array is [1.5, 0.8, -0.2, 1.1], the corresponding weight array is [0.4, -0.3, 0.2, 0.5], and the classification bias is b. r =0.1, number of dimensions is m r =4,λ r =0.01, then calculate the numerator first:

[0176]

[0177] Calculate the denominator again:

[0178]

[0179] Substitute into the calculation:

[0180]

[0181] Therefore, the support vector boundary coefficient value for the j-th class in the r=5 region is α. 5,j ≈0.1598. Table 5 lists the specific parameter settings for the support vector samples:

[0182] Table 5: Support Vector Sample Feature Value Table

[0183]

[0184] As shown in Table 5, the sample support vector parameters are explicitly set, enabling the formula to calculate specific boundary coefficient values.

[0185] The advantage of this formula is that by performing absolute averaging of the residuals after normalization between the feature values ​​and weights of multiple support vectors within a region, the most representative region boundary parameters can be extracted while minimizing classification bias. These parameters can then be used for subsequent classification decisions, thereby enabling the classifier to have significant discriminative ability for specific regions.

[0186] S503: Based on the regional classification boundary coefficient matrix, extract the set of coupling feature coefficients of the individual to be distinguished, calculate the classification probability, compare the kernel function mapping offset between categories, and generate an anomaly probability distribution map.

[0187] When extracting the set of coupling feature coefficients for individuals to be discriminated based on the region classification boundary coefficient matrix, it is necessary to obtain the structure-function coupling vector set of the target individual in the predefined region and substitute it into the aforementioned support vector machine model. Based on the category boundary parameters, the matching score for each category corresponding to the individual is calculated. Specifically, this is done by performing a region-by-region dot product between the individual feature vector and the support vector boundary coefficient matrix, and then calculating the distance difference. When calculating the classification probability, the multi-class matching scores are normalized using the Softmax method. For example, if the score for class A is 2.5, the score for class B is 1.8, and the score for class C is 0.5, then the calculations are as follows:

[0188]

[0189] The final probability distribution is [0.612, 0.304, 0.083]. The class is determined to be A based on the maximum probability. When comparing the kernel function mapping offsets between classes, the target individual region features are mapped to the high-dimensional space of the kernel function, and then compared with the Euclidean distance of the class model center point. The magnitude and direction of the offset vector are recorded to calculate the degree of anomaly. If the offset direction is opposite to the normal class direction in most regions, and the magnitude exceeds a set threshold (e.g., the offset exceeds the mean plus 1.5 times the standard deviation), it is marked as a potential anomalous region. When generating the anomaly probability distribution map, a two-dimensional probability heatmap corresponding to the MNI space is generated by combining the region offset direction and magnitude. The color value reflects the probability level of each region being an anomalous region, with a value range of 0 to 1, and is output in heatmap form. For example, if the anomalous probabilities in the left parietal lobe and right occipital lobe are 0.74 and 0.66 respectively, the corresponding region color blocks on the corresponding spatial coordinate heatmap are highlighted red and orange areas respectively.

[0190] like Figure 2 As shown, the medical image analysis system based on image processing includes:

[0191] The image preprocessing module is used to acquire the original diffusion tensor image sequence through the magnetic resonance device, call the anisotropic diffusion filtering algorithm to perform noise suppression operation on the original image, output the preprocessed image sequence, and pass it to the fiber node extraction module.

[0192] The fiber node extraction module is used to call the preprocessed image sequence, perform three-dimensional vector normalization on the main diffusion direction value, decompose the sequence using wavelet transform and extract the frequency band energy value, mark the voxels whose energy exceeds the dynamic threshold, generate the fiber node coordinate set, and pass it to the correlation analysis module.

[0193] The correlation analysis module is used to perform point-by-point Euclidean distance offset calculation on the resting and stimulated state node coordinates based on the fiber node coordinate set, perform Pearson correlation analysis on the distance offset, generate a node displacement-functional response correlation map, and pass it to the coupling modeling module.

[0194] The coupling modeling module is used to call the node displacement-functional response correlation map, perform three-dimensional spatial interpolation on the whole brain nodes, spatially register the interpolation results with the coordinates of the brain gyri boundary points of the structural magnetic resonance image, generate a structural-functional coupling strength distribution map, and pass it to the anomaly detection module.

[0195] The anomaly detection module is used to input the structure-function coupling strength distribution map into the support vector machine classifier, perform classification operations on the brain region feature vectors, compare them with the gray matter density threshold and white matter integrity index in the clinical diagnostic criteria parameters, and output an anomaly probability distribution map.

[0196] The above description is merely a specific embodiment of the present invention, but the scope of protection of the present invention is not limited thereto. Any variations or substitutions that can be easily conceived by those skilled in the art within the technical scope disclosed in the present invention should be included within the scope of protection of the present invention. Therefore, the scope of protection of the present invention should be determined by the scope of the claims.

Claims

1. A medical image analysis method based on image processing, characterized in that, Includes the following steps: S1: Acquire the original diffusion tensor image sequence through magnetic resonance equipment, call the anisotropic diffusion filtering algorithm to perform noise suppression operation on the original image, extract the tensor direction change value to construct the structure vector set, and generate the preprocessed image sequence; S2: Based on the preprocessed image sequence, perform three-dimensional normalization processing on the main diffusion direction value, use wavelet transform to decompose the sequence and extract the frequency band energy value, mark the voxels of energy exceeding the dynamic threshold, and output the fiber node coordinate set; S3: Based on the fiber node coordinate set, perform point-by-point Euclidean distance offset calculation on the node coordinates in the resting and stimulated states, and perform Pearson correlation analysis on the distance offset and functional magnetic resonance time series data to generate a node displacement-functional response correlation map. S4: Based on the node displacement-functional response correlation map, perform three-dimensional spatial interpolation on the whole brain nodes, align the interpolation results with the coordinates of the brain gyri boundary points in the structural magnetic resonance image, and output an individualized structure-function coupling strength distribution map; S5: Based on the individualized structure-function coupling strength distribution map, a support vector machine classifier is used to perform classification operations, compare with clinical diagnostic standard parameters, and output an abnormal probability distribution map.

2. The medical image analysis method based on image processing according to claim 1, characterized in that, The preprocessed image sequence includes signal-to-noise ratio parameters, anisotropy fractional matrix, and tensor eigenvalue distribution. The fiber node coordinate set specifically includes a high-frequency energy coordinate set, a low-frequency energy coordinate set, and an energy gradient coordinate set. The node displacement-functional response correlation map includes a displacement correlation coefficient matrix, functional activation threshold distribution, and spatiotemporal coupling matrix. The individualized structure-functional coupling strength distribution map specifically refers to the sulcus-gyrus coupling coefficient, gray-white matter interface intensity value, and cortical thickness mapping value. The abnormality probability distribution map includes a classification confidence distribution, a regional abnormality index matrix, and a pathological feature weight map.

3. The medical image analysis method based on image processing according to claim 1, characterized in that, The specific steps of S1 include: S101: Acquire the original diffusion tensor image sequence through magnetic resonance equipment, extract the tensor direction feature values ​​and corresponding gray-level distribution matrices at multiple locations in the image sequence, establish the diffusion direction distribution data structure, and generate diffusion direction change values; The change value of diffusion direction measures the degree of directional change by comparing it with the principal direction vector of the diffusion tensor, reflecting the spatial continuity and directional consistency characteristics of the tissue microstructure. S102: Call the diffusion direction change value, input it into the anisotropic diffusion filtering algorithm, iteratively calculate the grayscale correction value to suppress noise, and perform boundary judgment based on the tensor direction continuity to restrict pixel smoothing operation across the boundary and generate diffusion direction filtering weight value. The diffusion direction filtering weight value refers to the directional smoothing coefficient that is self-matched and adjusted according to the local structural features of the image during the anisotropic diffusion process. S103: Adjust the grayscale matrix values ​​of multiple frames in the original image sequence according to the diffusion direction filter weight value, reconstruct the tensor structure and output it uniformly to generate a preprocessed image sequence.

4. The medical image analysis method based on image processing according to claim 3, characterized in that, The specific steps of S2 include: S201: Based on the preprocessed image sequence, extract the main diffusion direction value, call the three-dimensional vector normalization method, calculate the vector magnitude of the main diffusion direction, process the components to a uniform scale, and generate normalized direction vector values. The principal diffusion direction value refers to the eigenvector corresponding to the principal eigenvalue in the diffusion tensor; S202: Based on the normalized direction vector value, the vector value sequence is decomposed by wavelet transform, frequency band energy value data is extracted, voxels with energy values ​​greater than the dynamic threshold are selected, and a voxel sequence exceeding the threshold is obtained. The wavelet transform uses the db4 wavelet basis function to perform a 3-level decomposition, and performs multi-scale decomposition processing on the three-dimensional image sequence to obtain the feature information of the image under different frequencies and scales. S203: Based on the above-threshold voxel sequence, locate the spatial position of the corresponding voxel in the three-dimensional coordinates of the image, summarize the spatial coordinate data, and generate a fiber node coordinate set.

5. The medical image analysis method based on image processing according to claim 4, characterized in that, The specific steps of S3 include: S301: Based on the fiber node coordinate set, calculate the Euclidean distance between node pairs in resting and stimulated states, compare the differences between node pairs and construct a spatial offset sequence to obtain the node distance offset sequence. S302: Call the node distance offset sequence, match the node distance offset with the time period of the functional magnetic resonance time series signal, filter the valid data, calculate the correspondence between the two sets of data according to the Pearson formula, and obtain the node correlation coefficient set. The Pearson formula measures the strength of the linear relationship between nodal displacement and its corresponding functional magnetic resonance time-series signal, revealing whether structural changes are accompanied by functional changes. S303: Based on the set of node correlation coefficients, pair the original coordinates with the correlation coefficients, construct a color level mapping relationship, and combine the coordinates to draw the color value spectrum corresponding to the node, generating a node displacement-functional response correlation spectrum.

6. The medical image analysis method based on image processing according to claim 5, characterized in that, The specific steps of S4 include: S401: Based on the node displacement-functional response correlation map, select node coordinates and response values, calculate the weight coefficients of the response values ​​in three-dimensional space, construct an interpolation model using a weighted average method to calculate the grid values, and generate a three-dimensional functional response interpolation value distribution. S402: Call the three-dimensional function response interpolation value distribution, filter out interpolation points in the structural magnetic resonance image template, match boundary point coordinates based on Euclidean distance, map interpolation values ​​to the gyri boundary, and obtain the gyri boundary matching value dataset; S403: Based on the brain gyrus boundary matching value dataset, traverse the interpolation change points, calculate and accumulate the boundary point value differences, normalize the cumulative change to form the corresponding intensity, and generate an individualized structure-function coupling intensity distribution map.

7. The medical image analysis method based on image processing according to claim 6, characterized in that, The weighting coefficients are calculated using the following formula: Among them, w ijk Δx represents the weight coefficients of 3D mesh points i, j, k. i Δy represents the displacement of the i-th node along the x-axis, in mm. j ||d represents the displacement of the j-th node along the y-axis, in mm. k || represents the magnitude of the displacement vector at the k-th node, in mm, δ m This represents the absolute deviation of the m-th response value from the reference value, expressed in Pa. The gradient magnitude of the m-th response value is represented by r, in Pa / mm. ijk This represents the coordinates of the current grid point i, j, k, in mm. m Represents the coordinates of the m-th node, in mm, R max The maximum interpolation influence radius is represented in mm and is obtained based on a backtracking experiment of real-time MRI interpolation error, ranging from 10 to 20 mm. N represents the total number of nodes, and u represents the pressure-displacement conversion coefficient, calibrated by brain tissue elastic modulus experiments, in mm / Pa, with a value range of 0.02-0.05 mm / Pa.

8. The medical image analysis method based on image processing according to claim 6, characterized in that, The specific steps of S5 include: S501: Based on the individualized structure-function coupling strength distribution map, divide the region to extract the coupling strength, match the spatial coordinates to classify the data, encode the region label to convert the vector format, unify the scale to form a vector set, and obtain the regional coupling feature coefficient set. S502: Call the set of regional coupling feature coefficients as input to the support vector machine, calculate the classification interval according to the training set labels, perform interval optimization on the distance between the sample and the classification surface, calculate the support vector boundary coefficient values, and obtain the regional classification boundary coefficient matrix; The support vector machine uses a radial basis function kernel, and the parameters are optimized in the training set through cross-validation. S503: Based on the region classification boundary coefficient matrix, extract the set of coupling feature coefficients of the individuals to be distinguished, calculate the classification probability, compare the kernel function mapping offset between categories, and generate an anomaly probability distribution map.

9. The medical image analysis method based on image processing according to claim 8, characterized in that, The support vector boundary coefficients are calculated using the following formula: Where, α r,j n represents the support vector boundary coefficient value of the j-th class in the r-th region. r γ represents the number of support vectors in the r-th region. r,i w represents the coupling feature coefficient of the i-th sample in the r-th region. r,i b represents the support vector weight corresponding to the i-th sample in the r-th region. r Let m represent the classification bias term obtained from training in the r-th region. r λ represents the total number of dimensions involved in the support vector weight calculation in the r-th region. r The weight regularization adjustment parameter represents the r-th region.

10. A medical image analysis system based on image processing, characterized in that, The system is used to implement the image processing-based medical image analysis method according to any one of claims 1-9, and the system comprises: The image preprocessing module is used to acquire the original diffusion tensor image sequence through the magnetic resonance device, call the anisotropic diffusion filtering algorithm to perform noise suppression operation on the original image, output the preprocessed image sequence, and pass it to the fiber node extraction module. The fiber node extraction module is used to call the preprocessed image sequence, perform three-dimensional vector normalization on the main diffusion direction value, decompose the sequence using wavelet transform and extract the frequency band energy value, mark voxels whose energy exceeds the dynamic threshold, generate a fiber node coordinate set, and pass it to the correlation analysis module. The correlation analysis module is used to perform point-by-point Euclidean distance offset calculation on the resting and stimulated state node coordinates based on the fiber node coordinate set, perform Pearson correlation analysis on the distance offset, generate a node displacement-functional response correlation map, and pass it to the coupling modeling module. The coupling modeling module is used to call the node displacement-functional response correlation map, perform three-dimensional spatial interpolation operation on the whole brain nodes, spatially register the interpolation results with the coordinates of the brain gyri boundary points of the structural magnetic resonance image, generate a structural-functional coupling strength distribution map, and transmit it to the anomaly detection module. The anomaly detection module is used to input the structure-function coupling strength distribution map into the support vector machine classifier, perform classification operations on the brain region feature vectors, compare them with the gray matter density threshold and white matter integrity index in the clinical diagnostic standard parameters, and output an anomaly probability distribution map.

Images