A spinal cord diffusion tensor imaging fidelity denoising reconstruction method and device

Abstract

This application discloses a method and apparatus for high-fidelity denoising and reconstruction of spinal cord diffusion tensor imaging. The method includes: segmenting the pathological spinal cord based on T2-weighted structural images to generate preprocessed T2 structural images, spinal cord masks, and gray-white matter boundary masks; performing cross-modal high-fidelity denoising in the DWI domain of the original DWI data with prior constraints from the T2 structural images to generate a boundary-constrained anisotropic denoised DWI image sequence; performing Gibbs ringing removal, eddy current and EPI geometric distortion correction on the boundary-constrained anisotropic denoised DWI image sequence to obtain a corrected high-fidelity reconstructed DWI dataset; and obtaining a quantitative diffusion parameter image group and a whole-spinal cord symmetric positive definite diffusion tensor field based on the corrected DWI dataset. This application utilizes the anatomical prior information of high-resolution T2-weighted structural images to achieve pathologically adaptive spinal cord segmentation, independent of geometric assumptions about the normal spinal cord, and capable of adapting to various pathological morphological changes in patients with chronic spinal cord injury.

Description

Technical Field

[0001] This application relates to the field of image processing technology, specifically to a method and apparatus for fidelity-preserving denoising and reconstruction of spinal cord diffusion tensor imaging. Background Technology

[0002] Diffusion tensor imaging (DTI) is a magnetic resonance imaging technique based on the Brownian motion of water molecules. By applying diffusion-sensitive gradients in different directions, it can non-invasively assess the diffusion characteristics of water molecules within biological tissues. DTI can reflect the microstructural information of white matter fiber bundles and calculate parameters such as fractional anisotropy (FA), mean diffusion coefficient (MD), axial diffusion coefficient (AD), and radial diffusion coefficient (RD). It has significant value in the clinical diagnosis and scientific research of neurological diseases.

[0003] However, diffusion-weighted imaging (DWI) requires a strong diffusion-sensitive gradient, resulting in a typically low signal-to-noise ratio (SNR). Furthermore, geometric distortions, Gibbs ringing artifacts, and Rician noise distribution introduced by echo-planar imaging (EPI) sequences make improving DTI image quality a critical technical challenge. Compared to brain DTI, spinal cord DTI faces even greater technical challenges. The spinal cord's transverse diameter is typically only 8-15 mm, and the white matter fiber tracts it contains, such as the corticospinal tract (CST), have an even smaller area, further reducing the SNR. Simultaneously, the spinal cord is surrounded by cerebrospinal fluid (CSF), and the pulsation and flow of CSF significantly contaminate the diffuse signal; the presence of vertebral bony structures causes magnetic susceptibility artifacts, further degrading image quality.

[0004] Patients with chronic spinal cord injury often exhibit atypical morphological changes in their spinal cord due to long-term mechanical compression or pathological alterations. These changes include, but are not limited to, spinal cord swelling, localized atrophy, syringomyelia, and flattening of the cross-section. These pathological morphological changes frequently cause automatic segmentation methods based on normal spinal cord geometry assumptions (such as elliptical or circular cross-section assumptions) to fail in clinical applications. Inaccurate segmentation directly leads to subsequent image denoising being performed in incorrect areas, further affecting the accuracy of diffusion parameter estimation.

[0005] Currently, DTI image denoising methods mainly include the following categories:

[0006] (1) Non-Local Means (NLM) filtering: Weighted averaging is used to denoise similar blocks in the image. However, standard NLM uses isotropic search and smoothing strategies, which can easily blur the gray and white matter boundaries at the boundaries of slender structures such as the spinal cord, and does not take into account pathological morphological changes.

[0007] (2) Marchenko-Pastur principal component analysis (MP-PCA): This method separates the signal subspace and noise subspace through local PCA decomposition. Although effective, in spinal cord regions with low SNR, setting the PCA threshold is difficult, and it is easy to lose the true pathological signal while removing noise.

[0008] (3) Deep learning methods, such as Patch2Self self-supervised denoising, require a large amount of training data and the model has limited generalization ability to pathological data.

[0009] The aforementioned existing technologies have the following common defects: (a) they do not make full use of the anatomical prior information of high-resolution structural images (such as T2-weighted images); (b) the denoising process does not adaptively adjust to the pathological morphological changes of the spinal cord; and (c) they do not apply constraints at the tissue boundaries, resulting in blurred image boundaries after denoising, which affects subsequent fiber tracking and quantitative analysis.

[0010] Taking the widely used open-source spinal MRI processing tool Spinal Cord Toolbox (SCT) as an example, its deep learning segmentation method relies on training data of normal spinal cord morphology, and the segmentation accuracy decreases under pathological conditions (such as swelling and atrophy); its denoising module relies on external MP-PCA tools, adopts an isotropic smoothing strategy, and the segmentation results and denoising process are independent of each other, failing to achieve effective fusion of cross-modal prior information. Summary of the Invention

[0011] The purpose of this invention is to provide a method for fidelity-preserving denoising and reconstruction of spinal cord diffusion tensor imaging, which at least solves one of the aforementioned technical problems.

[0012] One aspect of the present invention provides a method for fidelity-preserving denoising reconstruction of spinal cord diffusion tensor imaging, the method comprising:

[0013] Step 1: Acquire raw DWI data and T2-weighted structural images of the spinal cord region of the subject;

[0014] Step 2: Perform pathological spinal cord segmentation based on T2-weighted structural images to generate preprocessed T2 structural images, spinal cord masks, and gray-white matter boundary masks;

[0015] Step 3: Using T2 structured image prior constraints, perform cross-modal fidelity denoising in the DWI domain of the original DWI data to generate a DWI image sequence after boundary constraint anisotropic denoising.

[0016] Step 4: Perform Gibbs ringing removal, eddy current and EPI geometric distortion correction on the DWI image sequence after boundary constraint anisotropic denoising to obtain a corrected and faithfully reconstructed DWI dataset;

[0017] Step 5: Obtain the quantitative diffusion parameter image set and the whole spinal cord symmetric positive definite diffusion tensor field based on the corrected DWI dataset.

[0018] Optionally, step 2 includes:

[0019] Step 21: Perform N4 off-field correction and nonlocal mean noise suppression on the acquired T2-weighted structure image to obtain the preprocessed T2-weighted structure image;

[0020] Step 22: The spinal cord region is extracted using Otsu adaptive threshold segmentation combined with the level set evolution algorithm to generate a spinal cord mask. During the execution of the level set evolution algorithm, no elliptical or circular geometric constraints are applied to the cross-section of the spinal cord.

[0021] Step 23: Within the area defined by the spinal cord mask, based on the gray-scale difference between the high signal of gray matter and the low signal of white matter in the T2-weighted structural image, the gray-white matter interface of the spinal cord is detected by gradient analysis combined with the region growing algorithm, and a gray-white matter boundary mask is generated.

[0022] Optionally, step 3 includes:

[0023] Step 31: Perform rigid registration and spatial resolution matching between the preprocessed T2-weighted structured image and the b=0 reference image in the original DWI data to obtain the cross-modal spatial mapping relationship;

[0024] Step 32: Combine the gray-white matter boundary mask and the cross-modal space mapping relationship to construct an anisotropic smoothing kernel, thereby obtaining a voxel-level orientation-dependent smoothing kernel, wherein the smoothing intensity along the fiber direction is greater than the smoothing intensity perpendicular to the fiber direction.

[0025] Step 33: Set a linear decay weight based on the distance from the voxel to the gray-white matter boundary to obtain a smoothing kernel with boundary constraints, wherein the amplitude of the smoothing kernel at the boundary is zero;

[0026] Step 34: Iteratively smooth the DWI images in each diffusion gradient direction within the spinal cord mask area until convergence, generating a boundary-constrained anisotropic denoised DWI image sequence.

[0027] Optionally, step 4 includes:

[0028] Step 41: Perform Gibbs ringing removal processing on the DWI image sequence after boundary constraint anisotropic denoising to obtain an intermediate DWI image sequence with ringing artifacts eliminated;

[0029] Step 42: Perform eddy current correction and EPI geometric distortion correction sequentially on the intermediate DWI image sequence after eliminating ringing artifacts to obtain the corrected and faithfully reconstructed DWI dataset.

[0030] Optionally, step 5 includes:

[0031] Step 51: The original diffusion tensor matrix is ​​obtained by fitting the corrected and faithfully reconstructed DWI dataset using the weighted linear least squares method;

[0032] Step 52: Symmetricize the original diffusion tensor matrix, and perform eigenvalue decomposition and non-negative truncation on the symmetric tensor to obtain the corrected eigenvalue and eigenvector matrix.

[0033] Step 53: Reconstruct the whole spinal cord symmetric positive definite diffusion tensor field using the correction eigenvalues ​​and eigenvector matrices;

[0034] Step 54: Calculate the anisotropy fraction map, mean diffusion coefficient map, axial diffusion coefficient map, and radial diffusion coefficient map based on the whole spinal cord symmetric positive definite diffusion tensor field. The anisotropy fraction map, mean diffusion coefficient map, axial diffusion coefficient map, and radial diffusion coefficient map constitute a quantitative diffusion parameter image group.

[0035] Optionally, before step 34, step 3 further includes:

[0036] Based on the spinal cord mask corresponding to the registered T2-weighted structural image, the continuous topological skeleton curve of the spinal cord center is extracted. The topological skeleton is segmented and truncated using the gray and white matter boundary mask as a constraint, dividing it into three types of topological units: the main segment of the corticospinal tract, the gray matter decussation segment, and the terminal segment of the spinal cord marginal fiber.

[0037] For each type of topological unit, a diffusion gradient source tracing mapping relationship is established. The original diffusion signal source of each voxel is traced back along the topological skeleton curve under different diffusion gradient directions, and a corresponding signal source tracing index set for the topological unit is generated.

[0038] Based on the signal source index set, a topology bundle shielding rule is established. The diffuse signal between different topology units is only allowed to propagate along the skeleton curve in the positive direction, blocking the interaction of transverse diffuse signals across units and generating a topology bundle isolation mask.

[0039] A topological-spatial dual-locked weight field is constructed by superimposing gray-white matter boundary space constraints on a topological beam isolation mask as a basis. Smoothing kernel weight superposition calculation is performed only within the same topological unit, and the voxel weights across topological units are independent of each other.

[0040] Step 34 includes:

[0041] The anisotropic smoothing kernel is corrected based on the topological-spatial dual-locked weight field, and then global iterative smoothing is performed to generate a DWI image sequence with boundary-constrained anisotropic denoising.

[0042] Optionally, the step of establishing a diffusion gradient source tracing mapping relationship for each type of topological unit, tracing the source of the original diffusion signal of each voxel under different diffusion gradient directions in reverse along the topological skeleton curve, and generating a corresponding topological unit-specific signal source tracing index set includes:

[0043] Establish a local curve coordinate system with the central axis of the spinal cord topological skeleton curve. For any voxel x within the spinal cord mask, solve for the shortest vertical projection point p from voxel x to the topological skeleton curve. Extract the skeleton curve arc length parameter s(p) corresponding to the projection point p, and construct the voxel-skeleton arc length mapping function. ;

[0044] Obtain the preset gradient domain diffusion signal source tracing operator;

[0045] Traverse all diffusion gradient directions, use the source tracing operator to calculate the source signal amplitude of the corresponding topological unit under each gradient on a voxel-by-voxel basis, and index and bind all voxels in the same topological unit whose source amplitude is higher than the topological signal basis threshold to generate a signal source tracing index set exclusive to that topological unit.

[0046] Mutual exclusion tags are applied to the signal tracing index sets of different topological units to ensure that there are no overlapping voxels in the index sets of different topological units, and the tracing index sets of all topological units are constructed.

[0047] Optionally, the preset gradient domain diffuse signal source tracing operator formula is as follows:

[0048] ;

[0049] in, The original diffuse signal of voxel x under the i-th diffusion gradient direction; gradient The skeletal arc length corresponding to the origin point of the original signal of the fiber; α is the total arc length of a single topological unit skeleton; α is the topological attenuation coefficient. This is the gradient domain diffuse signal source tracing operator.

[0050] Optionally, before performing eigenvalue decomposition and nonnegativity truncation on the symmetricized diffusion tensor matrix, step 5 further includes:

[0051] Based on the topological bundle isolation mask, all voxels within the spinal cord mask are assigned to corresponding fiber topological units, and the topological unit assignment labels for each voxel are obtained.

[0052] Using all symmetric diffusion tensors within the same topological unit as a sample set, solve for the basis matrix of the tensor specific to that topological unit, and perform a basis alignment transformation on the symmetric diffusion tensors of each voxel within the unit to obtain the basis aligned tensor.

[0053] By applying inter-element tensor orthogonality constraints to the basis alignment tensor and removing the interference projection components brought by the tensor basis of other topological elements, a pure tensor with cross-element diffusion interference removed is obtained.

[0054] Step 52 includes:

[0055] The pure tensor matrix after removing cross-unit diffusion interference is symmetricized, and the symmetricized tensor is subjected to eigenvalue decomposition and non-negative truncation to obtain the corrected eigenvalue and eigenvector matrix.

[0056] This application also provides a spinal cord diffusion tensor imaging fidelity denoising and reconstruction device, the spinal cord diffusion tensor imaging fidelity denoising and reconstruction device comprising:

[0057] Image acquisition module, which is used to acquire raw DWI data and T2-weighted structural images of the spinal cord region of the subject;

[0058] The mask acquisition module is used to perform pathological spinal cord segmentation based on T2-weighted structural images, and generate preprocessed T2 structural images, spinal cord masks, and gray-white matter boundary masks.

[0059] The DWI image sequence acquisition module is used to perform cross-modal fidelity denoising in the DWI domain of the original DWI data with T2 structural image prior constraints, and generate a DWI image sequence after boundary constraint anisotropic denoising.

[0060] The reconstruction module is used to perform Gibbs ringing removal, eddy current and EPI geometric distortion correction on the DWI image sequence after boundary constraint anisotropic denoising, so as to obtain a corrected and faithfully reconstructed DWI dataset.

[0061] The image group acquisition module is used to acquire a quantitative diffusion parameter image group and a whole spinal cord symmetric positive definite diffusion tensor field based on the corrected DWI dataset.

[0062] This application utilizes prior anatomical information from high-resolution T2-weighted structural images to achieve pathologically adaptive spinal cord segmentation. It does not rely on geometric assumptions about the normal spinal cord and can adapt to various pathological morphological changes (swelling, atrophy, cavitation, flattening, etc.) in patients with chronic spinal cord injury, significantly improving segmentation accuracy. Attached Figure Description

[0063] Figure 1 This is a schematic flowchart of a method for fidelity-preserving denoising and reconstruction of spinal cord diffusion tensor imaging according to an embodiment of this application.

[0064] Figure 2 This is a schematic flowchart illustrating the detailed process of pathological spinal cord segmentation according to an embodiment of this application.

[0065] Figure 3 This is a schematic diagram of an anisotropic smoothing kernel according to an embodiment of this application.

[0066] Figure 4 This is a graph of the boundary attenuation function according to an embodiment of this application.

[0067] Figure 5 This is a comparison of FA parameter diagrams before and after denoising according to an embodiment of this application (based on real data from patients with chronic spinal cord injury).

[0068] Figure 6 This is a comparison of fiber tracking results before and after denoising in an embodiment of this application (based on real data from patients with chronic spinal cord injury). Detailed Implementation

[0069] To make the objectives, technical solutions, and advantages of this application clearer, the technical solutions in the embodiments of this application will be described in more detail below with reference to the accompanying drawings. In the drawings, the same or similar reference numerals denote the same or similar elements or elements having the same or similar functions throughout. The described embodiments are some, but not all, embodiments of this application. The embodiments described below with reference to the accompanying drawings are exemplary and intended to explain this application, and should not be construed as limiting this application. All other embodiments obtained by those skilled in the art based on the embodiments of this application without creative effort are within the scope of protection of this application. The embodiments of this application will be described in detail below with reference to the accompanying drawings.

[0070] like Figure 1 The spinal cord diffusion tensor imaging fidelity denoising reconstruction method shown includes:

[0071] Step 1: Acquire raw DWI data and T2-weighted structural images of the spinal cord region of the subject;

[0072] Step 2: Perform pathological spinal cord segmentation based on T2-weighted structural images to generate preprocessed T2 structural images, spinal cord masks, and gray-white matter boundary masks;

[0073] Step 3: Using T2 structured image prior constraints, perform cross-modal fidelity denoising in the DWI domain of the original DWI data to generate a DWI image sequence after boundary constraint anisotropic denoising.

[0074] Step 4: Perform Gibbs ringing removal, eddy current and EPI geometric distortion correction on the DWI image sequence after boundary constraint anisotropic denoising to obtain a corrected and faithfully reconstructed DWI dataset;

[0075] Step 5: Obtain the quantitative diffusion parameter image set and the whole spinal cord symmetric positive definite diffusion tensor field based on the corrected DWI dataset.

[0076] In this embodiment, the subjects collected in step 1 are mainly patients with chronic spinal cord injury. These patients often have various pathological changes in their spinal cord, such as swelling, local atrophy, syringomyelia, and flattening of the cross-section.

[0077] The imaging hardware of this application uses a 3.0T medical magnetic resonance imaging system with an 8-channel cervical and thoracic spinal phased array coil, and is equipped with an MR post-processing workstation that can export NIfTI and gradient text; a 1.5T device can be used for acquisition, but its imaging signal-to-noise ratio is too low and it is not the preferred device.

[0078] T2-weighted structural image data acquisition rules:

[0079] Scanning sequence: Fast spin echo (FSE) sequence was used;

[0080] Core scanning parameters: repetition time TR=3000 ms, echo time TE=120 ms, imaging matrix 256×256, slice thickness 3 mm, slice spacing 0.3 mm; standardized positioning before acquisition: subject supine, cervical and thoracic sponge pads for support and fixation, headrest locked, laser uniformly locked the C1-T4 interval after three-plane pre-scanning, ensuring that T2 and DWI scan FOV and anatomical start and end layers are completely consistent.

[0081] Scan coverage: cervical spinal cord and upper thoracic spinal cord, corresponding anatomical segments C1-T4 level;

[0082] Image function: High-resolution anatomical images provide a clear overall outline of the spinal cord and the difference in gray-to-white matter grayscale contrast, serving as the core anatomical basis for subsequent adaptive pathological segmentation, gray-to-white matter boundary extraction, and cross-modal prior constraints.

[0083] Complete Composition and Acquisition Specifications of DWI Raw Data: DWI, or diffusion-weighted imaging raw data, is the foundation of DTI reconstruction. The complete dataset must include two types of images: b-values, gradient directions, and scan sequences, strictly adhering to the standards specified in the disclosure document.

[0084] Base reference image: at least one b=0s / The diffusion-free weighted reference image serves as the DWI baseline mode for subsequent rigid registration and spatial resolution matching of T2 structural images;

[0085] Multi-gradient diffusion-weighted images: Configure 15–64 b-value-weighted images with different diffusion gradient directions; the preferred clinical implementation uses 30 diffusion gradient directions.

[0086] The range of b-values ​​is uniformly set to 800-1000s / In the preferred embodiment, b = 1000s / ;

[0087] Scanning sequence: Single excitation echo planar imaging SS-EPI sequence;

[0088] DWI scan parameters: TR=8000ms, TE=85ms, imaging matrix 128×128, field of view FOV=220×220 The layer thickness is 3mm and the layer spacing is 0.3mm. One b=0 reference image is acquired simultaneously in both positive and negative phase encoding directions. The two phase encoding directions are opposite and the other scanning parameters are completely consistent, which are used to solve the magnetic susceptibility field map.

[0089] The implementation fully covers the numerical range boundaries of the claims. The lower limit scheme has 15 gradients and b=800, which is suitable for elderly patients who cannot tolerate long scans; the upper limit scheme has 64 gradients and b=1000, which is suitable for high-precision quantitative research. The TR / TE / slice thickness / FOV parameters remain completely unchanged in the three implementations.

[0090] The sequence naturally suffers from problems such as low signal-to-noise ratio (SNR), EPI geometric distortion, Gibbs ringing artifacts, cerebrospinal fluid pulsation signal contamination, and vertebral magnetic susceptibility artifacts. These are also the original data defects that need to be addressed in the subsequent denoising and distortion correction steps of this application.

[0091] After acquisition, the entire set of T2-weighted structural images and the complete set of raw DWI sequences are uniformly stored in the medical standard NIfTI image format, along with the scan parameters, gradient direction vectors, and b-value text information. The images are uniformly compressed into NIfTI.nii.gz format, with b-values ​​separated by commas on a single line, and gradient vectors in 3-column gx / gy / gz text. The folder is hierarchically structured as: subject number / T2, DWI subdirectory, and a four-level image verification log file is generated simultaneously for traceability.

[0092] Verification process:

[0093] (1) Verify that the T2 structural image and the DWI scan cover the spinal cord segments completely and there is no vertical offset; if the vertical offset of the vertebral body exceeds one complete slice, it is considered unqualified and a repositioning scan is required.

[0094] (2) Verify the integrity of the DWI data, confirm that it includes the complete b=0 reference image and all gradient direction weighted images, and that there are no missing tortuosities; the number of lines of gradient vector text must be equal to the number of DWI weighted images, and mismatch of line numbers indicates that the parameters are disordered;

[0095] (3) Check that the gradient vector file corresponds one-to-one with the number of diffusion-weighted images and that there are no disordered gradient parameters; if the spinal cord signal fluctuation of the b=0 image is greater than 20% and the signal loss of the spinal cord cross section exceeds 1 / 3, the image is deemed unqualified and needs to be re-sampled.

[0096] (4) Remove unqualified images with obvious motion artifacts, missing scan sections, or large areas of signal loss, and re-scan and acquire them.

[0097] Step 1 Output Results: After completing the acquisition and verification in this step, output two complete sets of original image datasets as input for Step 2:

[0098] Spinal cord anatomy dataset: raw images of T2-weighted FSE structural images of C1-T4 segments;

[0099] Diffusion imaging dataset: Complete SS-EPI sequence DWI raw data, including a single b=0 reference image and multiple gradient orientations at b=1000 s / Diffusion-weighted image, along with accompanying gradient vector and b-value parameter files.

[0100] In this embodiment, step 2 completes the entire process of anatomical segmentation using the T2-weighted structural image output from step 1. No geometric constraints of elliptical or circular cross-sections of the normal spinal cord are applied throughout the process. This can adapt to various pathological morphologies of the spinal cord in patients with chronic spinal cord injury, such as swelling, atrophy, cavitation, and flattening of the cross-section. It outputs three types of prior masks: standardized preprocessed T2 images, complete spinal cord region masks, and gray-white matter boundary masks, providing anatomical prior support for cross-modal boundary constraint denoising in step 3.

[0101] See Figure 2 In this embodiment, step 2 includes:

[0102] Step 21: Perform N4 off-field correction and nonlocal mean noise suppression on the acquired T2-weighted structure image to obtain the preprocessed T2-weighted structure image;

[0103] Specifically, the N4 bias correction operation: T2 weighted structural images are affected by the non-uniformity of the magnetic resonance radio frequency field, resulting in signal intensity shifts across the entire image domain. The N4 bias correction algorithm is used to equalize the signal across the entire domain of the original T2 image, eliminating the brightness deviation caused by the non-uniformity of the radio frequency field, unifying the signal benchmark of the spinal cord region, cerebrospinal fluid, and surrounding soft tissues, and avoiding the problems of missed segmentation and over-segmentation due to brightness shifts in subsequent segmentation.

[0104] The N4 correction uses B-spline order 3 and convergence threshold 1e-4; the NLM filter search window is 5×5×5, similarity threshold is 0.7, smoothing weight is 0.4, and the operation is only performed on the pre-defined region of the spinal cord mask in T2 images. Peripheral cerebrospinal fluid and vertebral bodies are not included in the filtering calculation.

[0105] The grayscale parameters mentioned above are all calculated based on the normalized grayscale values ​​of the preprocessed T2 image. The normalization method is to linearly map the signal intensity within the spinal cord mask to the [0,1] interval, with the minimum grayscale value within the spinal cord mask being 0 and the maximum grayscale value being 1. The similarity of the NLM filter uses the normalized Euclidean distance metric of image blocks. That is, for two image blocks of the same size, the sum of the squares of the grayscale differences of the corresponding pixels is calculated and divided by the total number of pixels in the block. Image blocks with a normalized distance less than the similarity threshold are included in the weighted average calculation. The smoothing weight is the normalized coefficient of the weight of each similar block, which is used to control the filtering intensity.

[0106] Non-local mean (NLM) noise suppression: T2 images after field correction contain random scanning noise. NLM filtering is performed to suppress the noise. This step only eliminates irrelevant random noise, fully preserving the texture and grayscale contrast between spinal cord gray matter, white matter, and cerebrospinal fluid, without blurring the subtle anatomical boundaries of the spinal cord.

[0107] Finally, after N4 bias correction and NLM noise suppression, a preprocessed T2-weighted structural image is obtained, which serves as the subsequent input base image.

[0108] Step 22: The spinal cord region is extracted using Otsu adaptive threshold segmentation combined with the level set evolution algorithm to generate a spinal cord mask. During the execution of the level set evolution algorithm, no elliptical or circular geometric constraints are applied to the cross-section of the spinal cord.

[0109] Specifically, this process consists of two steps: coarse segmentation and precise segmentation based on level set evolution. No assumptions are made regarding the geometry of the spinal cord cross-section throughout the process, making it suitable for various pathological spinal cord malformations.

[0110] Otsu adaptive thresholding coarse segmentation: Based on the global grayscale distribution of the preprocessed T2-weighted structural image, the optimal segmentation threshold is automatically calculated using the Otsu adaptive thresholding method to divide the high-signal region into the low-signal background. The coarse segmentation result includes the main spinal cord, surrounding cerebrospinal fluid (CSF), and adjacent soft tissue, resulting in a coarse segmentation region of interest containing the target spinal cord.

[0111] After coarse segmentation, a connected component filtering step is added to retain only connected regions with an axial cross-sectional pixel area between 30 and 300 pixels, automatically removing large connected regions such as cerebrospinal fluid and subcutaneous fat, thus reducing the amount of subsequent level set irrelevant calculations.

[0112] Level set evolution precise segmentation:

[0113] (1) Initialize the level set function within the region of interest obtained from the Otsu coarse segmentation. For the symbolic distance function, construct a segmentation optimization energy functional:

[0114] ;

[0115] in, The signal intensity of the T2 structure image is localized. , These represent the average signal strength within and outside the horizontal set, respectively. , , , This is the weight hyperparameter for the energy term.

[0116] Iteratively minimizing the energy functional driving the evolution of the level set, the optimization process relies solely on local grayscale consistency and edge alignment constraints to update the contour, without imposing standard geometric constraints such as ellipse or circle on the spinal cord cross-section throughout the process, and can naturally adapt to various pathological irregular morphologies such as spinal cord swelling, local atrophy, syringomyelia, and cross-sectional flattening.

[0117] In this embodiment, the preferred weight values ​​are: contour length weight μ∈[0.1,0.5], preferably 0.25; region area weight v∈[0.01,0.2], preferably 0.1; and grayscale fitting weight. = The value is fixed at 1; the upper limit of the level set iteration is set to 500 times. If the overlap rate of 99.9% is not reached, the iteration will be forcibly terminated to avoid infinite iteration loops of pathological cavity images.

[0118] All weight parameters of the aforementioned energy functional are adapted to the normalized grayscale image in the [0,1] interval. The original T2 image needs to be normalized before being input into the level set evolution process. The axial cross-sectional pixel area of ​​the aforementioned connected component screening is 30–300 pixels, which is adapted to the acquisition parameters of the 256×256 imaging matrix and 3mm layer thickness of the T2 structural image in this embodiment, corresponding to a physical area of ​​approximately 2.7–27mm². If structural images of other resolutions are used, the area threshold needs to be calculated proportionally according to the pixel physical size.

[0119] Level set evolution convergence determination: convergence is determined when the overlap rate of the contour region is higher than 99.9% in two consecutive iterations, the evolution is stopped, and the closed contour after convergence is extracted to generate a complete spinal cord mask that fits the patient's true pathological morphology.

[0120] Step 23: Within the area defined by the spinal cord mask, based on the gray-scale difference between the high signal of gray matter and the low signal of white matter in the T2-weighted structural image, the gray-white matter interface of the spinal cord is detected by gradient analysis combined with the region growing algorithm, and a gray-white matter boundary mask is generated.

[0121] Specifically, the anatomical grayscale features are based on the following: T2-weighted structural images have stable grayscale differentiation features, with spinal cord gray matter showing a relatively high signal and spinal cord white matter showing a relatively low signal. This grayscale difference is the core basis for gray and white matter boundary detection.

[0122] Boundary detection execution process:

[0123] (1) Limit the calculation range: only perform calculations within the output spinal cord mask area, shielding irrelevant backgrounds such as cerebrospinal fluid, vertebral body, and soft tissue outside the spinal cord, and reducing irrelevant gradient interference.

[0124] (2) Gradient analysis for initial edge detection: Gradient analysis (optional Canny edge detection, Gaussian smoothing kernel σ=1.0, gradient high threshold is 0.7 times the peak value of the gradient histogram, and low threshold is 0.4 times the high threshold) is used to calculate the local gray value amplitude of the T2 image. The gradient peak position is the potential interface between gray matter and white matter.

[0125] (3) Refined Boundary Correction for Region Growth: The region growth algorithm is executed using the gradient peak point as the seed point, growing bidirectionally along the gray-level gradient direction to distinguish between gray-matter connected regions and white-matter connected regions, and to remove isolated pseudo-edges and minor noise artifacts. The gray-level similarity threshold for region growth is set to 0.08. Growth is stopped immediately when the gray-level difference between the voxel to be grown and the seed exceeds the threshold or when it touches the outer boundary of the spinal cord mask, to prevent large-area overflow across gray and white matter. This similarity threshold is calculated based on the normalized T2 gray-level value in the aforementioned [0,1] interval, and the initial gray-level of the seed point is the normalized gray-level value corresponding to the gradient peak point.

[0126] Extract the interface contours of the gray and white matter regions, and binarize them to generate a gray-white matter boundary mask. The boundary mask precisely marks every gray-white matter boundary position.

[0127] After completing the entire process in step 2, output three types of standardized data:

[0128] 1. Preprocessed T2-weighted structural image: Equalized low-noise anatomical image after N4 field correction and NLM noise suppression;

[0129] 2. Patient-specific spinal cord mask: Generated through geometrically unconstrained level set evolution, adapting to various spinal cord pathological morphologies, precisely wrapping the entire spinal cord parenchyma region;

[0130] 3. Gray and white matter boundary mask: Based on T2 gray-level gradient and region growth, it completely marks the gray and white matter interface inside the spinal cord, providing a spatial constraint benchmark for subsequent boundary truncation smoothing kernel.

[0131] For extreme pathological samples with syringomyelia, severe atrophy, and flattened cross-sections, the level set evolution in this step does not force a closed contour size. If the minimum axial width of the segmented spinal cord mask is less than 5 pixels, a local grayscale secondary Otsu resegmentation is automatically initiated: First, the local resegmentation region is defined by extending 10 pixels outward from the failed region output by the level set. Otsu adaptive threshold calculation is re-executed within the local region, and the optimal segmentation threshold is solved only for the grayscale distribution of the local region. After the secondary segmentation result is combined with the original spinal cord mask, the level set evolution is re-executed, and the upper limit of iteration is adjusted to 200 times. Finally, a modified spinal cord mask adapted to the pathological morphology is generated to avoid missing small lesions.

[0132] Step 3 includes:

[0133] Step 31: Perform rigid registration and spatial resolution matching between the preprocessed T2-weighted structured image and the b=0 reference image in the original DWI data to obtain the cross-modal spatial mapping relationship;

[0134] In this embodiment, the b=0 image is a baseline DWI image without diffusion gradient, with a higher signal-to-noise ratio than the gradient-weighted DWI. It serves as the T2 structural registration reference, eliminating positional offset and resolution differences between the T2 anatomical image and the DWI diffusion image.

[0135] Affine rigid transformation is used to map the T2-weighted structural image to the b=0 image space. Only translation and rotation are performed, without introducing nonlinear deformation, to ensure the true correspondence of the spinal cord anatomical structure and avoid deformation from destroying the spatial position of the gray and white matter boundaries.

[0136] High-resolution T2 images are resampled to ensure that their voxel sizes are completely consistent with the original DWI data, achieving a one-to-one correspondence between T2 anatomical coordinates and DWI diffuse voxel coordinates, thus ensuring that subsequent gray and white matter boundary masks can be accurately mapped to all DWI gradient images.

[0137] Output mapping matrix: Stores the rigid registration transformation matrix and the resampled T2 anatomy map, serving as the cross-modal localization reference for subsequent smoothing kernel space constraints. In this embodiment, affine registration uses normalized mutual information (NMI) as the cost function, and the iterative convergence threshold is set to 10. -4Resampling uses trilinear interpolation; registration deviations at the fault boundaries of the upper edge of C1 and the lower edge of T4 are corrected separately to ensure anatomical alignment of the first and last segments; the same set of registration transformation matrices can be reused for DWI at any gradient level of 15 / 30 / 64.

[0138] Step 32: Combine the gray-white matter boundary mask and the cross-modal space mapping relationship to construct an anisotropic smoothing kernel, thereby obtaining a voxel-level orientation-dependent smoothing kernel, wherein the smoothing intensity along the fiber direction is greater than the smoothing intensity perpendicular to the fiber direction.

[0139] In this embodiment, a voxel-level orientation-dependent anisotropic smoothing kernel K(x) is constructed by combining a gray-white matter boundary mask. The smoothing intensity along the fiber direction is greater than that perpendicular to the fiber direction, and the anisotropic ratio r takes a value of 2-5.

[0140] For any voxel X within DWI, the local spinal cord long axis is extracted based on the registered T2 spinal cord contour. (i.e., the corticospinal tract fibers are approximately parallel), and the local voxel normal vector n(x) is calculated by combining the gray and white matter boundary mask.

[0141] Construction of the smoothing kernel covariance matrix: A two-dimensional anisotropic Gaussian kernel model is used, and the covariance expression is as follows:

[0142] ;

[0143] in, The standard deviation is smoothed parallel to the fiber direction. The standard deviation is the smoothed deviation perpendicular to the fiber direction; the ratio of the two is... The value of r is limited to the range of 2 to 5. In this embodiment, r=2.5 is preferred. =2.0 mm =0.8 mm. The three-dimensional Gaussian smoothing kernel discrete sampling size is fixed at 7×7×7 voxels, the kernel truncation threshold is 3 times the standard deviation, and the weight is directly set to 0 if it exceeds the range; the standard deviation is distributed along the long axis of the spinal cord and in two orthogonal directions perpendicular to the front and back / perpendicular to the left and right in the three-dimensional voxel coordinates, to fully realize three-dimensional anisotropic smoothing, not limited to two-dimensional planar calculation.

[0144] In alternative embodiments, this application further includes the following steps:

[0145] To address the inherent Rician noise distribution characteristics of spinal cord DWI images, noise bias correction is applied to the aforementioned anisotropic smoothing kernel. The expression for the corrected smoothing kernel is as follows:

[0146] ;

[0147] in The aforementioned voxel-level anisotropic Gaussian smoothing kernel, The standard deviation of background noise is taken from the standard deviation of signal intensity in the cerebrospinal fluid region within a 30-pixel range outside the spinal cord. (x) represents the b=0 reference image signal intensity corresponding to voxel x; when the voxel b=0 signal intensity is less than twice the noise standard deviation, the lower limit of the correction factor is clamped at 0.2 to avoid over-correction in low signal areas and the introduction of new artifacts.

[0148] Regional spinal cord long axis The calculation is performed as follows: Principal component analysis is performed on the cross-section of the spinal cord mask at each axial level, and the direction of the first principal component is extracted as the short axis direction of the spinal cord cross-section at that level. Combined with the longitudinal interlayer position of the spinal cord, the voxel-level three-dimensional fiber long axis vector is obtained by correcting the tangent direction of the three-dimensional skeleton. The three-dimensional skeleton is the single-pixel center curve output by the subsequent topological skeleton extraction step. The long axis direction corresponding to any voxel is taken as the skeleton tangent unit vector at the shortest projection point on the skeleton, ensuring that the long axis direction is completely consistent with the spinal cord fiber course.

[0149] The covariance matrix of the three-dimensional anisotropic Gaussian kernel is expanded to a 3×3 form, where the dimension parallel to the fiber's long axis is adopted. The remaining two perpendicular dimensions orthogonal to the major axis are both adopted. The two perpendicular directions are obtained through orthogonal basis decomposition, corresponding to the anterior-posterior and left-right directions of the spinal cord cross section, respectively, ensuring that the smooth nucleus stretches along the long axis of the fiber and converges in the vertical direction in three-dimensional space, which is consistent with the physical characteristics of white matter fiber diffusion.

[0150] Nuclear characteristics description: A larger standard deviation in the direction parallel to the fibers can effectively suppress random noise in low SNR areas of the spinal cord; a smaller standard deviation in the direction perpendicular to the fibers avoids excessive smoothing across the gray and white matter boundaries, which can cause tissue contrast loss. This is different from the circular isotropic smoothed nuclei used in NLM and MP-PCA.

[0151] like Figure 3 This is a schematic diagram of the anisotropic smooth nucleus in this embodiment, where (a) is a schematic diagram of the spinal cord cross-section, (b) is the isotropic circular smooth nucleus used in the existing NLM / MP-PCA, (c) is the anisotropic smooth nucleus (r=2.5) stretched along the fiber direction in this invention, and (d) is the final smooth nucleus after superimposed gray and white matter boundary truncation constraint.

[0152] Step 33: Set a linear decay weight based on the distance from the voxel to the gray-white matter boundary to obtain a smoothing kernel with boundary constraints, wherein the amplitude of the smoothing kernel at the boundary is zero;

[0153] Specifically, it iterates through all voxels within the spinal canopy, calculates the Euclidean distance d from each voxel X to the nearest gray-white matter boundary, and sets a preset distance threshold. =1.5mm.

[0154] Linear decay weighting function: using the linear decay formula:

[0155] ;

[0156] when =0 (the voxel falls on the gray-white matter boundary line). =0, the smoothing kernel amplitude is completely zero, no smoothing is done on the boundary signal, and the true diffusion difference at the gray and white matter boundary is fully preserved;

[0157] when The weight of the smoothing kernel increases linearly with distance, and the farther away from the boundary, the stronger the noise reduction and smoothing effect.

[0158] when ≥ Weight =1, the smooth kernel maintains its complete anisotropic morphology, has no attenuation constraint, and effectively suppresses deep spinal cord noise. For transitional mixed voxels that fall on both sides of the gray and white matter boundary, the minimum distance from both boundaries is used to calculate the linear attenuation weight, avoiding weight distortion caused by multiple boundary superpositions.

[0159] Weight fusion operation: converting linearly decaying weights The resulting anisotropic smoothing kernel K(x) is multiplied voxel-by-voxel to generate the final smoothing kernel with boundary protection constraints. Compared with hard thresholding and exponential decay schemes, linear decay has a smooth transition and does not produce weight abrupt artifacts at the threshold.

[0160] like Figure 4 Figure (a) shows the boundary decay function curve of this embodiment. Figure (a) compares the three weight change trends of linear decay, hard threshold and exponential. This invention adopts linear decay without abrupt changes. Figure (b) shows the two-dimensional weight visualization effect of gray and white matter boundary.

[0161] In alternative embodiments, this application further includes the following steps:

[0162] (1) Construction of cross-modal signal fidelity benchmark

[0163] Based on the registered T2-weighted structural image and the DWI b=0 reference image, normal white matter and normal gray matter voxels were screened within the spinal cord mask (the screening criteria were consistent with the gray-level threshold of the aforementioned topological unit division). A univariate linear regression was performed with the T2-normalized gray level of the normal voxels as the independent variable and the b0 signal intensity as the dependent variable. The linear mapping coefficients and intercepts were obtained, and a voxel-level signal prior benchmark was constructed.

[0164] ;

[0165] in, For voxel x Normalized grayscale value, This serves as the baseline value for the diffusion signal fidelity of the voxel. Voxels from pathological regions are not included in the linear regression calculation to ensure that the baseline signal conforms to the physical correspondence with normal tissue.

[0166] (2) Automatic identification and adaptive weight allocation of pathological regions: Based on T2 normalized grayscale, pathological regions of chronic spinal cord injury are automatically identified and divided into three categories according to the degree of grayscale deviation from the normal range. Voxel-level fidelity weights are allocated to each category. :

[0167] Normal tissue area: grayscale is within the normal white or gray matter range, fidelity weight. =0.2, prioritizing noise reduction strength;

[0168] Transitional zone of lesion: grayscale deviation from the normal range by 10%-20%, or located within 1 voxel of the lesion edge, with fidelity weighting. =0.5, balancing noise reduction and signal fidelity;

[0169] Core lesion area: grayscale greater than 0.8 (corresponding to edema, syringomyelia) or less than 0.25 (corresponding to atrophy, softening lesions) and ≥5 conjoints, with fidelity weighting. =0.8, prioritizing the preservation of pathological diffusion characteristics.

[0170] Step 34: Iteratively smooth the DWI images in each diffusion gradient direction within the spinal cord mask area until convergence, generating a boundary-constrained anisotropic denoised DWI image sequence.

[0171] Specifically, within the spinal cord mask, the DWI image for each diffusion gradient direction is independently and iteratively smoothed until the iteration converges, generating a DWI image sequence with boundary-constrained anisotropic denoising:

[0172] Processing range constraints: Smoothing operations are performed only within the binary spinal cord mask, and the extraspinal cerebrospinal fluid, vertebral bodies, and soft tissue regions are all masked. Irrelevant background noise does not participate in the convolution calculation, reducing invalid operations and avoiding peripheral artifacts from interfering with the spinal cord parenchyma signal.

[0173] Independent gradient processing: The DWI weighted images corresponding to the 15-64 diffusion gradient directions are all denoised separately through iterative processes. The smoothing results are not shared between different gradients, ensuring that the diffusion signal is independently preserved in each gradient dimension.

[0174] Each iteration uses anisotropic diffusion partial differential equations for updating, and the update formula is:

[0175] ;

[0176] in, Let be the diffuse signal intensity at voxel x in the t-th iteration, and Δt be the iteration step size, which is fixed at 0.15 in this embodiment. For divergence operators, For gradient operators, The constrained anisotropic smoothing kernel function is given for voxel x. During the iteration process, Neumann boundary conditions are used at the boundary of the spinal cord mask, that is, the gradient value of the boundary voxel is taken as the gradient value of the adjacent voxel inside the boundary to avoid distortion of the boundary signal.

[0177] In an alternative embodiment, each iteration uses an anisotropic diffusion equation with fidelity constraints to perform the update, and the update formula is:

[0178] ;

[0179] The first term is the anisotropic smoothing term after Rician noise correction, used to suppress random noise in the image; the second term is the cross-modal fidelity constraint term, which constrains the deviation between the denoised signal and the T2 prior benchmark, preventing the pathological signal from being over-smoothed; Δt is the iteration step size, which is fixed at 0.12 in this embodiment. For divergence operators, For gradient operators, The constrained smoothing kernel function after Rician noise correction for voxel x;

[0180] Iterative convergence criterion: The root mean square error (RMSE) of the global domain between two consecutive iterations is used as the criterion, and the convergence threshold is set to 10. -3 If the RMSE is less than the threshold for two consecutive iterations, the noise is considered to have converged, and the smoothing iteration is stopped.

[0181] The final output is a full gradient dimension boundary-constrained anisotropic denoised DWI image sequence: covering the b=0 reference image and all b-value weighted diffusion images, and each image is subjected to fiber-oriented smoothing + gray-white matter boundary protection denoising.

[0182] The accompanying intermediate parameter file includes the fiber orientation vectors of each voxel, the voxel-gray matter boundary distance matrix, the gradient-by-gradient iterative convergence log, and the cross-modal rigid registration transformation matrix, which can be used in subsequent steps.

[0183] Step 4 includes:

[0184] Step 41: Perform Gibbs ringing removal processing on the DWI image sequence after boundary constraint anisotropic denoising to obtain an intermediate DWI image sequence with ringing artifacts eliminated;

[0185] Correction was performed on Gibbs ringing artifacts in DWI images. These artifacts are concentrated at the junctions of high-contrast tissues such as spinal cord and cerebrospinal fluid, and gray matter and white matter. They cause distortion of diffuse signals in the junction areas and interfere with the acquisition of true voxel signals.

[0186] The correction is limited to the area inside the spinal cord mask, and only the image within the spinal cord parenchyma is corrected. Background areas such as cerebrospinal fluid, vertebral bodies, and peripheral soft tissues are not included in the calculation, so as to avoid secondary interference of background signals to the effective signals of the spinal cord.

[0187] Fourier domain ringing filtering was performed on DWI images for each diffusion gradient dimension to suppress oscillation artifacts caused by alternating light and dark areas at tissue edges. This eliminated Gibbs oscillation interference while preserving the true diffusion signals of the spinal cord gray and white matter boundaries and the fine corticospinal tracts. After processing, an intermediate DWI image sequence with ringing artifacts removed was generated and proceeded to the next correction step. Hanning window low-pass filtering was used in the Fourier domain, with a frequency cutoff threshold set to 0.75 of the image sampling Nyquist frequency. Fourier transforms and inverse transforms were performed independently along the layer, row, and column dimensions of the three-dimensional image, using 64×64×64 overlapping blocks to avoid new oscillation artifacts at block boundaries. The overlap size of adjacent blocks was 8 voxels. After block processing, the overlapping areas were fused using a linear weighting method, with higher weights closer to the block center and zero weights at the edges, eliminating splicing artifacts at block boundaries. Before correction, the signal within the spinal cord mask is uniformly normalized to the [0,1] interval. After correction, the original signal amplitude is restored to prevent over-correction of high-signal cerebrospinal fluid regions.

[0188] Step 42: Perform eddy current correction and EPI geometric distortion correction sequentially on the intermediate DWI image sequence after eliminating ringing artifacts to obtain the corrected and faithfully reconstructed DWI dataset.

[0189] For intermediate DWI image sequences after Gibbs ringing removal, eddy current correction and EPI geometric distortion correction were performed continuously. Both types of distortion are inherent imaging defects of the SS-EPI sequence: First, the diffusion gradient switching generates an eddy current field, causing translation and local distortion of DWI images in different gradient directions, which destroys the voxel spatial correspondence between multi-gradient images; Second, the vertebral bony structure brings local magnetic susceptibility differences, which induces geometric stretching and compression deformation in the EPI readout direction, resulting in distortion of the spinal cord fiber anatomy.

[0190] The correction operation employs a combination of field map correction and eddy current correction tools: First, a magnetic susceptibility field distribution is constructed based on the b=0 reference image to complete EPI geometric distortion correction and correct overall spinal cord deformation. The magnetic susceptibility field map is constructed using the phase difference method of positive and negative phase b0 images: Fourier transform and phase unwrapping are performed on two positive and negative phase b=0 images respectively, the phase difference between the two phase encoding directions is calculated, and the global magnetic susceptibility offset field is obtained by combining the readout direction bandwidth and pixel size conversion. The number of iterations for solving the field map is set to 30, the phase unwrapping uses the minimum cost flow algorithm, and the deformation fitting uses cubic B-spline interpolation. Only a single generation of the full field map based on the b=0 image is required, and all 15 / 30 / 64 gradient-weighted DWIs reuse the same magnetic susceptibility field map, eliminating the need for gradient-by-gradient reconstruction of the field map and significantly reducing the computational load.

[0191] Then, by combining the coding parameters of each diffusion gradient, the image shift and distortion caused by gradient eddy currents are eliminated frame by frame. The entire correction process strictly preserves the amplitude of the true diffuse signal after denoising in step 3, without compromising the fidelity of the gray-white matter boundary. The number of iterations for solving the field map is set to 30, and cubic B-spline interpolation is used for deformation fitting; only a single generation of the full field map based on b=0 is required, and all 15 / 30 / 64 gradient-weighted DWIs reuse the same magnetic susceptibility field map, eliminating the need for gradient-by-gradient reconstruction of the field map and significantly reducing the computational load. Eddy current correction uses a third-order polynomial to fit the gradient-induced spatial shift, and the shift coefficients are solved separately for each gradient.

[0192] After correction, the spatial alignment accuracy of all gradient images is uniformly verified to ensure that the voxel coordinates of the same anatomical location are matched in all DWI images. Finally, the complete corrected DWI dataset is output, which includes the b=0 reference image and the weighted DWI images of all 15-64 diffusion gradient directions, as the standard input data for diffusion tensor fitting in step 5.

[0193] Alignment quantization judgment criteria: The deviation of the three-dimensional coordinates of the same anatomical voxel in all gradient images is ≤0.5 voxels to be qualified; if the deviation of a single weighted DWI exceeds the threshold, the second fine correction of the eddy current is re-performed until the deviation requirement is met.

[0194] After all corrections are completed, a unified global signal amplitude calibration is performed: the average signal of gray matter within the spinal cord mask is used as a reference to correct the intensity drift between gradients, eliminate the global signal offset caused by eddy current and deformation corrections, and ensure that the amplitudes of diffuse signals between different gradients are comparable.

[0195] Step 5 includes:

[0196] The 3D tensor fitting employs overlapping 64×64×32 voxel blocks for computation, with a 2-voxel overlap region set at the block boundaries to eliminate tomographic artifacts. This method is suitable for full-size 3D DWI big data of the spinal cord and avoids computer memory overflow. The tensor calculation results in the overlapping regions are fused using a weighted average method, with the weights negatively correlated with the distance from the voxel to the block center, ensuring the continuity and smoothness of the tensor field at the block boundaries. For the 15 / 30 / 64 gradient datasets, only the dimensions of the equation system change; the fitting weights and convergence thresholds do not require modification.

[0197] Step 51: The original diffusion tensor matrix is ​​obtained by fitting the corrected and faithfully reconstructed DWI dataset using the weighted linear least squares method;

[0198] Specifically, tensor fitting operations are performed on each individual voxel within the spinal cord mask, and the fitting follows the diffusion signal fundamental model:

[0199] ;

[0200] in, Let be the intensity of the diffuse signal of the voxel under the i-th diffusion gradient direction; The voxel signal corresponding to the diffusion-weighted reference image with b=0; The diffusion weighting coefficient has a value range of 800 s / -1000 s / ; Let be the unit direction vector corresponding to the i-th diffusion gradient; D is the original diffusion tensor matrix to be solved.

[0201] A weighted least squares algorithm is used to solve for the multi-gradient signals simultaneously, reducing the interference of low signal-to-noise ratio voxels on the tensor solution, and obtaining the uncorrected original diffusion tensor matrix D for each voxel. Fitting calculations are performed only within the generated spinal cord mask area, and calculations are directly skipped for cerebrospinal fluid, vertebral body, and peripheral soft tissue regions, reducing unnecessary computation.

[0202] The weighting is assigned based on the reciprocal of the spinal cord SNR of each gradient image; the lower the SNR, the smaller the gradient weight. The signal-to-noise ratio of the image in the i-th diffusion gradient direction is calculated using the following formula:

[0203] ;

[0204] in, Let be the voxel average signal intensity of the spinal cord mask-covered region in the i-th gradient image. Let be the standard deviation of the signal intensity of the background region outside the spinal cord in the i-th gradient image. The background region is the cerebrospinal fluid region within a range of 30 pixels outward from the outer boundary of the spinal cord mask.

[0205] The fitting weights are inversely proportional to the signal-to-noise ratio. The fitting weights for the i-th gradient are: ;

[0206] Where N is the total number of diffusion gradient directions, and the sum of the weights of all gradients after normalization is 1.

[0207] The convergence threshold for iterative solutions of linear equation systems is set to 10. -4 The maximum iteration limit is 300 times. If the limit is exceeded, singular value decomposition will be automatically used as a fallback solution. Only voxels within the spinal cord mask participate in the calculation. Cerebrospinal fluid and bony regions are directly skipped from the fitting calculation.

[0208] Step 52: Symmetricize the original diffusion tensor matrix, and perform eigenvalue decomposition and non-negative truncation on the symmetric tensor to obtain the corrected eigenvalue and eigenvector matrix.

[0209] Specifically, the original tensor matrix is ​​affected by imaging noise and gradient error, and contains asymmetric components, which does not meet the symmetry requirements of biological diffusion physics. Therefore, a symmetry transformation is performed:

[0210] ;

[0211] To symmetricize the diffusion tensor matrix, the asymmetric noise components of the original matrix are eliminated, ensuring that the tensor satisfies the physical symmetry premise of water molecule dispersion, and providing input for subsequent eigenvalue decomposition.

[0212] Step 53: Reconstruct the whole spinal cord symmetric positive definite diffusion tensor field using the correction eigenvalues ​​and eigenvector matrices;

[0213] For symmetric matrices Eigenvalue decomposition yields three sets of ordered eigenvalues ​​that satisfy the following size relationship: ≥ ≥ Simultaneously, it is equipped with a feature vector matrix V;

[0214] Considering that noise can cause some voxel eigenvalues ​​to be negative (physically meaningless, as the dispersion coefficient of water molecules cannot be negative), a non-negativity truncation constraint is applied to all eigenvalues:

[0215] ;

[0216] in It is a very small positive number. =10 -6 After truncation, all eigenvalues ​​are greater than or equal to This ensures that the diffusion tensor possesses physical positive definiteness from the source. If the original characteristics... If the value is much less than 0 (less than -0.01), the voxel gradient signal is determined to have severe distortion. The abnormality is marked and replaced with the correction feature value of the neighboring voxel in the same unit to prevent discrete abnormal noise from appearing in the parameter map.

[0217] The selection rule for neighboring voxels is as follows: voxels within three topological units above and below the abnormal voxel along the fiber long axis, excluding voxels from other units; the alternative method is to take the median of the feature values ​​of all neighboring normal voxels, and the corresponding feature vector is normalized using the original feature vector of the abnormal voxel to ensure the continuity of the fiber direction.

[0218] Using truncated and corrected eigenvalues , Using the original eigenvector matrix V, reconstruct the diffusion tensor that satisfies the symmetric positive definite constraint:

[0219] ;

[0220] (·) represents a diagonal matrix, with the diagonal elements being the three corrected eigenvalues. This is a symmetrical positive definite diffusion tensor field that conforms to the physical laws of biological diffusion, covering all voxels within the spinal cord mask to form a continuous tensor field throughout the spinal cord.

[0221] Step 54: Calculate the anisotropy fraction map, mean diffusion coefficient map, axial diffusion coefficient map, and radial diffusion coefficient map based on the whole spinal cord symmetric positive definite diffusion tensor field. The anisotropy fraction map, mean diffusion coefficient map, axial diffusion coefficient map, and radial diffusion coefficient map constitute a quantitative diffusion parameter image group.

[0222] Based on the corrected eigenvalues, calculate the four core clinical quantitative parameters: FA, MD, AD, and RD.

[0223] Axial dispersion coefficient AD: represents the degree of dispersion along the principal fiber direction.

[0224] ;

[0225] Radial dispersion coefficient RD: represents the average dispersion degree perpendicular to the fiber direction.

[0226] ;

[0227] Mean diffusion coefficient (MD): Overall average diffusion level of voxels.

[0228] ;

[0229] Anisotropy fraction (FA): A core evaluation indicator reflecting the regularity of white matter fibers and pathological damage.

[0230] ;

[0231] in, For tensor trace;

[0232] After calculating the parameter values ​​for each voxel, four sets of quantitative diffusion parameter images—FA, MD, AD, and RD—are generated to form a complete set of quantitative diffusion parameter images.

[0233] Clinical standard range of parameters: FA∈[0,1], MD unit mm² / s ranges from 0 to 3×10. -3 The unit AD / RD is mm², ranging from 0 to 2.5 × 10⁻⁶. -3 The output parameter graph is uniformly linearly normalized to 0~255 grayscale for visualization, and the original floating-point tensor field retains the true physical values ​​without compression.

[0234] Grayscale mapping only performs linear scaling, without piecewise nonlinear transformation; when reading ROI delineation in radiology, the original floating-point physical values ​​are directly retrieved, and the visualization of grayscale will not cause deviation in the quantitative measurement of lesions.

[0235] like Figure 6 The images show a comparison of fiber tracking results before and after denoising in this embodiment, drawn based on real patient data of chronic spinal cord injury. (a) is a fiber tracking image after isotropic denoising using traditional MP-PCA, showing multiple breaks in the corticospinal tract, poor fiber continuity, and severe loss of fiber tracts in the local spinal cord injury area. (b) is a fiber tract reconstructed after topology + boundary dual constraint denoising in this application, which significantly improves fiber integrity, significantly reduces the number of tracking breakpoints, and clearly shows subtle changes in fiber course in the injury area.

[0236] In existing technologies, traditional anisotropic smoothing relies solely on the single-layer constraint of gray and white matter boundaries. The diffusion signals of the three types of topological units within the spinal cord—the main corticospinal tract, gray matter crossings, and marginal fiber terminals—interfere with each other. The smoothing process can smooth out local pathological diffusion gradients across fiber units, causing distortion of local FA and MD values ​​in spinal cord injuries.

[0237] The diffusion direction of water molecules within different fiber topological units varies significantly. The uniform smooth nucleus lacks unit differentiation, making it easy to treat the lesion diffusion signal of a certain unit as noise smoothing, thus losing the microscopic damage characteristics corresponding to tiny spinal cord cavitation and local atrophy.

[0238] Relying solely on the two-dimensional boundary constraints of gray and white matter without utilizing the continuous topological structure of the three-dimensional fibers of the spinal cord, cross-unit transverse signal crosstalk cannot be blocked, reducing the continuity of fiber tracking and the accuracy of parameter measurement.

[0239] Existing methods such as NLM, MP-PCA, and basic anisotropic denoising do not distinguish fiber topology partitions, and cannot achieve smooth noise reduction within the same unit or signal isolation across units, resulting in insufficient identification of subtle pathological damage.

[0240] For the reasons stated above, prior to step 34, step 3 further includes:

[0241] Based on the spinal cord mask corresponding to the registered T2-weighted structural image, the continuous topological skeleton curve of the spinal cord center is extracted. The topological skeleton is segmented and truncated using the gray-white matter boundary mask as a constraint, dividing it into three types of topological units: the main segment of the corticospinal tract, the gray matter decussation segment, and the terminal segment of the spinal cord marginal fibers. The three types of units are radially distributed along the cross-section of the spinal cord and extend continuously along the longitudinal direction of the skeleton without spatial overlap.

[0242] The quantification and determination rules are as follows:

[0243] (1) Gray scale benchmark: All judgments are based on the normalized T2 gray value in the range of [0,1], with the minimum gray value in the spinal cord mask being 0 and the maximum gray value being 1;

[0244] (2) Reference benchmark: Using the spinal cord topological skeleton curve as the central axis, calculate the shortest Euclidean distance from any voxel in the spinal cord mask to the central axis, and denot it as the voxel radial distance r(x); For each axial cross section, calculate the equivalent radius R of the spinal cord mask at that level, that is, the radius of a circle with the same area as the spinal cord at that level.

[0245] (3) Unit determination priority: The units are determined in the order of gray matter cross segment, corticospinal tract trunk segment, and peripheral fiber terminal segment. Voxels are assigned to only the units with the highest priority.

[0246] (4) Specific judgment conditions:

[0247] ① Gray matter cross segment: Voxels with T2 normalized gray level higher than 0.6 and radial distance r(x) in the interval [0, 1 / 3R] correspond to the butterfly-shaped structure region of spinal cord gray matter;

[0248] ②Corticospinal tract trunk segment: Voxels that are not identified as gray matter crossing segments, have a T2 normalized gray value of less than 0.4, and whose radial distance r(x) is in the range of [1 / 3R, 2 / 3R], correspond to the corticospinal tract course region of the white matter of the lateral columns of the spinal cord;

[0249] ③ Spinal cord marginal fiber terminal segment: Voxels that are not classified into the first two categories and whose distance from the outer boundary of the spinal cord mask is ≤2 voxels, corresponding to the fine fiber terminal region of the peripheral white matter of the spinal cord.

[0250] Digital labels are assigned to all voxels within the spinal cord mask to complete the partition marking, and finally, a three-dimensional continuous spinal cord topology skeleton curve file and a partition label map of the entire spinal cord topology unit are obtained.

[0251] Specifically, a three-dimensional morphological refinement operation is performed based on the spinal cord mask, the peripheral voxels are stripped, and a continuous three-dimensional central topological skeleton curve that fits the longitudinal fibers of the spinal cord is extracted; the three-dimensional morphological refinement iteration is limited to 200 times, and 3×3×3 structural elements are used for layer-by-layer erosion, retaining only the single-pixel central skeleton; the cervical spinal cord segment and thoracic spinal cord segment skeleton extraction share the same set of refinement parameters, without the need for segmental adjustment.

[0252] Using a gray-white boundary mask as a spatial truncation constraint, the skeleton curve is segmented for cutting.

[0253] Based on anatomical distribution, three types of non-overlapping topological regions are divided: the main segment of the corticospinal tract, the gray matter decussation segment, and the terminal segment of the spinal cord marginal fibers.

[0254] Quantitative classification criteria: 1. Cortical spinal cord trunk: T2 gray level below the normalized threshold of 0.4, a ring-shaped region with a radius of 1 / 3 to 2 / 3 from the center of the spinal cord; 2. Gray matter crossover: T2 gray level above the normalized threshold of 0.6, a region with a radius of 0 to 1 / 3 from the center; 3. Peripheral region: a region on the periphery of the outer boundary of the spinal cord with a radius of 0 to 2 voxels. The three types of regions are mutually exclusive and do not overlap.

[0255] Digital labels are assigned to all voxels within the spinal cord mask to complete the partition marking, and finally a three-dimensional continuous spinal cord topology skeleton curve file and a whole spinal cord topology unit partition label map are obtained.

[0256] For each type of topological unit, a diffusion gradient source tracing mapping relationship is established. The original diffusion signal source of each voxel is traced back along the topological skeleton curve under different diffusion gradient directions, and a corresponding signal source tracing index set for the topological unit is generated.

[0257] Establish a local curvilinear coordinate system using the topological skeleton. For any voxel X within the spinal cord mask, calculate the shortest vertical projection point p from the voxel to the skeleton and extract the projected arc length. Construct a mapping function ;

[0258] Loading gradient domain diffuse signal source tracing operator:

[0259] ;

[0260] Traverse all 15-64 diffusion gradients and calculate the source signal amplitude of the unit at each gradient on a voxel-by-voxel basis. ;

[0261] Grouped by topological unit, voxels with source amplitudes higher than the base threshold are bound to a dedicated index set; different unit indices are marked mutually exclusive, and there are no overlapping voxels. Finally, the voxel-skeleton arc length mapping matrix L(x) and the independent dedicated signal source index set (mutually exclusive marking) of each topological unit are obtained; and the three-dimensional mapping map of the full gradient source amplitude is obtained.

[0262] Based on the signal source index set, a topology bundle shielding rule is established. The diffuse signal between different topology units is only allowed to propagate along the skeleton curve in the positive direction, blocking the interaction of transverse diffuse signals across units and generating a topology bundle isolation mask.

[0263] Specifically, the beam splitting and shielding rules are set: diffuse signals are only allowed to propagate in the longitudinal direction along the topological skeleton, blocking lateral interaction across topological units;

[0264] Traverse all spinal cord voxels. If two voxels belong to different topological units, reset the smoothing weight between them to 0 and prohibit cross-unit weighting. Voxels within the same unit remain connected and can be smoothed.

[0265] A global binary isolation mask is generated, where a mask value of 1 represents a smooth region within the same cell and 0 represents an isolated region across cells, ultimately generating a topological bundle isolation mask. Transition voxels at the boundary between two types of topological cells are not directly set to 0 / 1; instead, a linear transition weight of 0-1 is used to gradually change the transition weight, avoiding artificial artifacts caused by abrupt changes in cell boundary smoothness. The mask is stored using uint8 integers, with 0 / 255 corresponding to isolated / connected regions, respectively.

[0266] A topological-spatial dual-locked weight field is constructed by superimposing gray-white matter boundary space constraints on a topological beam isolation mask as a basis. Smoothing kernel weight superposition calculation is performed only within the same topological unit, and the voxel weights across topological units are independent of each other.

[0267] Using topological beam separation mask as the underlying locking basis, the weight for permanent locking across cell regions is 0;

[0268] Within the same topological unit, the gray and white matter boundary linearly decays weight w(d): the closer the voxel is to the boundary, the lower the weight, and the weight returns to zero when the boundary d=0;

[0269] Two constraints are fused per voxel to generate a dual-locked weight field that simultaneously satisfies cross-cell isolation and gray-white matter boundary protection. After fusion, the weight values ​​are forcibly truncated within the interval [0,1], and when floating-point overflow occurs during calculation, they are uniformly clamped to the boundary value to prevent noise amplification caused by smoothing weights greater than 1.

[0270] The weight fusion adopts a product rule, that is, the final locking weight of any voxel = topological beam separation mask weight × gray-white matter boundary linear attenuation weight w(d); after fusion, the weight values ​​are forcibly truncated to the [0,1] interval, and values ​​outside the range are uniformly clamped to the corresponding boundary value to prevent the smoothing weight from being greater than 1 and causing signal amplification distortion.

[0271] In alternative embodiments, this application adopts the following solution:

[0272] The weight fusion adopts a product rule to construct a topology-space-signal triple-locked weight field. The final smoothed effective weight of the voxel after fusion is calculated by the following formula:

[0273] ;

[0274] in, The value of the topology bundle isolation mask at voxel x is the topology bundle isolation weight, which is used to block smoothing across topology units; The linear decay weight of the gray-white matter boundary corresponding to voxel x is the boundary constraint weight defined in step 33, and the value at the boundary is 0. The voxel x-adaptive fidelity weights range from 0.2 to 0.8; the calculated... Forced truncation to the [0,1] interval, values ​​outside the range are clamped to the interval boundary to avoid signal amplification and distortion caused by the smoothing weight exceeding 1.

[0275] Step 34 includes:

[0276] The anisotropic smoothing kernel is corrected based on the topological-spatial dual-locked weight field, and then global iterative smoothing is performed to generate a DWI image sequence with boundary-constrained anisotropic denoising.

[0277] Specifically, the smoothing kernel is modified by multiplying the basic anisotropic smoothing kernel with the double-locked weight field voxel by voxel.

[0278] In an alternative embodiment, the anisotropic smoothing kernel corrected for Rician noise is multiplied voxel-by-voxel by a triple-locked weighted field to complete the smoothing kernel correction.

[0279] Existing spinal cord DWI denoising methods generally assume that the noise is Gaussian distributed and directly use anisotropic smoothing or nonlocal mean filtering. However, the noise in DWI images is essentially Rician distributed, and the spinal cord itself has a small transverse diameter and a signal-to-noise ratio much lower than that of the brain. Rician noise in low-signal areas (spinal cord margins, lesion areas) will produce a significant positive signal bias. Direct smoothing will lead to systematic deviations in quantitative parameters such as FA and MD, and the measured values ​​will deviate from the true physiological level.

[0280] In the alternative embodiments of this application, based on the original dual geometric constraints, three types of designs are integrated: Rician noise adaptation, cross-modal signal fidelity benchmark, and pathological weights, to form a fidelity denoising framework with triple locking of topological isolation, spatial edge preservation, and signal fidelity preservation.

[0281] This alternative embodiment addresses the following problems in the existing field of spinal cord DTI denoising:

[0282] First, the noise model is not well matched. The existing scheme defaults to Gaussian noise, which does not adapt to the inherent Rician noise distribution of DWI. The positive bias of the signal in the low signal-to-noise ratio region of the spinal cord will lead to systematic deviations in the quantitative parameters.

[0283] Second, the cross-modal prior utilization level is low. The T2 high-resolution structural image is only used for geometric constraints such as spatial masking, without exploring the physical correlation at the signal level. The denoising process has no signal truth anchor point, which easily leads to overall signal drift.

[0284] Third, the smoothing strategy is rigid. The uniform smoothing intensity cannot distinguish between normal tissue and pathological areas. There is an inherent contradiction that the noise reduction in normal areas is insufficient and the pathological features in lesion areas are over-smoothed. In addition, the original geometric constraints only limit the smoothing range and cannot control the smoothing amplitude, resulting in a single constraint dimension.

[0285] This alternative embodiment directly embeds the Rician noise bias correction term into the anisotropic smoothing kernel, dynamically adjusts the smoothing kernel amplitude based on the voxel b0 reference signal intensity and the standard deviation of background noise, and sets a correction lower limit clamp for the extremely low signal region. This can correct the signal amplitude bias caused by Ricean noise from the bottom layer of the smoothing kernel, improve the physical authenticity of the diffuse signal in the low signal-to-noise ratio region, and reduce the systematic error in the subsequent calculation of diffuse quantitative parameters from the input level.

[0286] Simultaneously, based on the linear fitting of T2 normalized grayscale and b0 signal of normal tissue, a voxel-level signal prior benchmark is constructed as a fidelity constraint term for iterative denoising. The high-contrast grayscale features of T2 image are used to automatically divide regions of different pathological grades and assign gradient fidelity weights. The smoothing intensity is adjusted in reverse. This can both anchor the amplitude range after denoising through cross-modal signal priors, avoid the overall signal drift caused by iterative smoothing, and improve the signal consistency between images of different diffusion gradients, and adaptively reduce the effective smoothing intensity in pathological areas to retain the real diffusion features corresponding to lesions and avoid pathological information being mistakenly identified as noise filtering.

[0287] Finally, by using the product rule, the three different types of constraints—topological isolation, boundary attenuation, and fidelity adaptation—are integrated into a unified voxel-level weight field. This corrects the effective strength of the smoothing kernel in one go. The three types of constraints act synchronously throughout the entire smoothing iteration process, achieving a complete constraint closed loop from limiting the smoothing range to finely controlling the smoothing intensity. The constraint granularity is refined from the region level to the voxel level, while avoiding the accumulation of errors caused by step-by-step processing. This achieves precise protection of pathological signals while blocking cross-fiber unit signal crosstalk and preserving the contrast between gray and white matter boundaries, thus balancing noise reduction intensity and signal fidelity.

[0288] This alternative embodiment breaks through the common technical bias in the field. The industry generally assumes that T2-weighted images and DWI-weighted images have different imaging weights and signal mechanisms, and can only be used for spatial registration and tissue segmentation. This alternative embodiment not only expands the fidelity benchmark at the signal amplitude level based on the common physical characteristics of normal tissues, but also relies on high signal-to-noise ratio T2 structural images to identify lesions to guide DWI denoising, thus expanding the utilization dimension of cross-modal prior information.

[0289] On the other hand, existing technologies generally process noise correction, boundary protection, and denoising as independent steps in a series. This alternative embodiment integrates noise model adaptation, geometric topology constraints, and signal fidelity control into a single weight field and smoothing kernel for synchronous operation. It is a collaborative design of multi-dimensional constraints rather than a simple splicing of conventional technologies. It can simultaneously balance the three contradictory goals of noise reduction, boundary preservation, and pathological preservation, and solves the inherent problem in the field of spinal cord DTI denoising that blurs the boundary when reducing noise and retains noise when preserving the boundary.

[0290] The operation is performed only within the spinal cord mask area, and the DWI image is independently iterated and smoothed for each diffusion gradient direction. The 3D DWI image is processed in memory in blocks, with a single block size of 64×64×32 voxels. The block boundaries are smoothed by overlapping tomography to avoid block boundary artifacts. 15 / 30 / 64 gradient sequences can be iterated in blocks in parallel. The number of gradients only affects the total operation time, and the smoothing logic remains unchanged.

[0291] After each iteration, the root mean square error (RMSE) of the entire image is calculated. When the RMSE of two consecutive iterations is less than 10... -3 Once convergence is determined, smoothing is stopped.

[0292] like Figure 5 The images shown are comparisons of FA parameters before and after denoising in this embodiment. They are drawn based on images of real patients with chronic spinal cord injury. (a) is the original FA image without topological-spatial dual-constraint denoising, where the gray and white matter boundaries are blurred and small pathological gradients are smoothed. (b) is the FA image after processing by the method of this application, where the contrast between gray and white matter is significantly improved and microscopic damage features such as syringomyelia and local atrophy are completely preserved.

[0293] In this embodiment, a diffusion gradient source tracing mapping relationship is established for each type of topological unit. The original diffusion signal source of each voxel under different diffusion gradient directions is traced backward along the topological skeleton curve, generating a corresponding topological unit-specific signal source tracing index set, including:

[0294] Establish a local curve coordinate system with the central axis of the spinal cord topological skeleton curve. For any voxel x within the spinal cord mask, solve for the shortest vertical projection point p from voxel x to the topological skeleton curve. Extract the skeleton curve arc length parameter s(p) corresponding to the projection point p, and construct the voxel-skeleton arc length mapping function. ;

[0295] Specifically, a local curved coordinate system is established based on the central axis of the spinal cord topology, abandoning the traditional Cartesian straight-line distance calculation method and adapting to the curved and long fiber structure of the spinal cord.

[0296] Traverse every voxel X inside the spinal cord mask and find the shortest vertical projection point p from voxel X to the topological skeleton curve.

[0297] Calculate the arc length parameter s(p) corresponding to the projection point p along the topological skeleton curve. The starting point of the arc length is uniformly set to the origin of the skeleton at the C1 end of the cervical spinal cord.

[0298] Construct a one-to-one mapping function L(x)=s(p) to enable any voxel to quickly query its own position arc length along the fiber longitudinal direction.

[0299] Obtain the preset gradient domain diffusion signal source tracing operator;

[0300] In this embodiment, the preset gradient domain diffusion signal source tracing operator formula is as follows:

[0301] ;

[0302] in, The original diffuse signal of voxel x under the i-th diffusion gradient direction; gradient The skeletal arc length corresponding to the origin point of the original signal of the fiber; The total arc length of a single-segment topological unit skeleton; The topological attenuation coefficient is preferably in the range of 0.8-1.2. For tracing the source of diffuse signals in the gradient domain; To output the voxel X skeleton arc length. The preferred value is 1.0; the fixed threshold for the topological signal basis is 0.35 times the average amplitude of the global image source tracing. If the amplitude is lower than this value, it is directly judged as a cross-unit interference signal and is not included in the index set.

[0303] Gradient corresponds to the origin of the original signal in the fiber. Determined according to the following rules: for the i-th diffusion gradient direction Calculate the angle between the gradient and the tangent direction at each point on the topology skeleton curve, and take the skeleton point with the smallest angle as the signal source point in that gradient direction; if multiple points have the same angle, take the point whose arc length is closest to the midpoint of the arc length of the topology unit as the source point. Total arc length of a single topology unit skeleton. The total length of the curve from the start point to the end point of the corresponding skeleton segment of this unit is given by the topology attenuation coefficient α, which is preferably set to 1.0.

[0304] Based on the three categories of topological units, the total arc length of the skeleton for each category of units was calculated. ;

[0305] A fixed topological signal base threshold is used as the critical value for distinguishing between effective signals and weak interference signals.

[0306] Traverse all diffusion gradient directions, use the source tracing operator to calculate the source signal amplitude of the corresponding topological unit under each gradient on a voxel-by-voxel basis, and index and bind all voxels in the same topological unit whose source amplitude is higher than the topological signal basis threshold to generate a signal source tracing index set exclusive to that topological unit.

[0307] Specifically, iterate through each diffusion gradient direction in turn. ;

[0308] For the current gradient, traverse all voxels within the spinal cord mask and substitute them into the source tracer operator to calculate the source tracer magnitude of each voxel under this gradient. ;

[0309] Voxels are grouped by topological unit labels, and all voxels within the same unit whose traceability amplitude is greater than the base threshold are selected.

[0310] The selected high-amplitude reliable voxels are indexed and bound to generate a temporary index set for the topological unit and the gradient.

[0311] Mutual exclusion tags are applied to the signal tracing index sets of different topological units to ensure that there are no overlapping voxels in the index sets of different topological units, and the tracing index sets of all topological units are constructed.

[0312] Traverse the index set of all topological units and check if there are any overlapping voxels across units;

[0313] If a voxel appears in multiple cell index sets at the same time, only the index of its own topological cell is retained, and the voxel's tag in other cells is deleted to achieve mutual exclusion of indexes.

[0314] Add independent unit identifier labels to the index set of each type of topological unit to distinguish between the three types of indexes: corticospinal tract trunk, gray matter decussation, and marginal fiber terminal.

[0315] Integrate all mutually exclusive unit indexes to build a globally unified traceability index database.

[0316] Before performing eigenvalue decomposition and nonnegativity truncation on the symmetricized diffusion tensor matrix, step 5 further includes:

[0317] Based on the topological bundle isolation mask, all voxels within the spinal cord mask are assigned to corresponding fiber topological units, and the topological unit assignment labels for each voxel are obtained.

[0318] Specifically, iterate through all voxels within the spinal cord mask and read the cell labels of the topological bundle isolation mask;

[0319] Each voxel is assigned a unique topological unit affiliation label, which is divided into three categories: the main corticospinal tract trunk, the gray matter decussation segment, and the terminal segment of the spinal cord marginal fiber.

[0320] Remove cross-element boundary mixed voxels (belonging to two or more types of elements at the same time) and mark them separately as boundary transition voxels, and then perform smoothing constraint processing separately;

[0321] Generate a voxel-topological unit mapping table, recording the unit number to which each voxel belongs.

[0322] Using all symmetric diffusion tensors within the same topological unit as a sample set, solve for the basis matrix of the tensor specific to that topological unit, and perform a basis alignment transformation on the symmetric diffusion tensors of each voxel within the unit to obtain the basis aligned tensor.

[0323] Specifically, the original tensors are grouped according to three types of topological units, and all voxel tensors within the same unit constitute the unit's exclusive sample set;

[0324] For each set of tensor samples, perform principal component analysis (PCA) to extract the diffuse feature principal vectors of that unit and construct the tensor basis matrix specific to that topological unit. (k=1,2,3, corresponding to the three types of fiber units respectively);

[0325] The basis matrix solution process is as follows: First, each 3×3 symmetric diffusion tensor within a cell is vectorized into a 6-dimensional independent component vector (containing three diagonal elements and three independent off-diagonal elements), forming an n×6 sample matrix (n is the number of voxels in the cell); Principal component analysis is performed on the sample matrix to extract the first three principal component eigenvectors. Each eigenvector is reconstructed into a 3×3 symmetric matrix, and then orthogonalized to obtain three orthogonal unit basis vectors, which together constitute a 3×3 orthogonal basis matrix Bk; If the number of effective voxel samples in a certain topological cell is less than 10, the basis vectors of adjacent cells are automatically reused by distance-weighted fusion to avoid the basis matrix being singular and uncomputable.

[0326] PCA retains the first three principal eigenvectors to form a 3×3 basis matrix; if the number of effective voxel samples in a certain topological unit is less than 10, the basis of the adjacent unit is automatically reused for weighted substitution to avoid the basis matrix being singular and unable to be calculated.

[0327] The basis matrix dimension is consistent with the 3×3 dimension of the diffusion tensor, and only the most representative diffuse feature basis of the unit is retained, while the random noise component within the unit is filtered out.

[0328] Store the base of the corticospinal tract trunk separately Gray cross-substrate , edge fiber terminal base .

[0329] By applying inter-element tensor orthogonality constraints to the basis alignment tensor and removing the interference projection components brought by the tensor basis of other topological elements, a pure tensor with cross-element diffusion interference removed is obtained.

[0330] For any voxel, read its topological unit k and retrieve its corresponding dedicated basis. ;

[0331] Establish a tensor basis alignment transformation equation to project the original tensor onto the feature basis of this element:

[0332] ;

[0333] By completing the substrate alignment, the voxel tensor is fully mapped to the feature space of its own fiber unit, thus weakening the hybrid components brought by the foreign unit substrate.

[0334] Boundary transition mixed voxels are aligned using a weighted average of adjacent two-unit bases, with the weights determined by the voxel's distance from both sides of the boundary. This process ultimately yields the whole spinal cord base-aligned tensor field. (The tensor matrix of each voxel after unit basis alignment).

[0335] Constructing orthogonal constraints for topological elements: The basis matrices of different topological elements are mutually orthogonal. Iterative orthogonalization convergence threshold 10 -5 This ensures that the inner product of different unit substrates approaches zero infinitely, eliminating residual aliasing components across unit substrates.

[0336] Inter-unit basis orthogonalization was performed using the Gram-Schmidt orthogonalization method, processed sequentially in the order of gray matter cross segments, corticospinal tract trunk segments, and peripheral fiber terminal segments: first, the original basis of the first type of unit was preserved; then, the basis of the second type of unit was projected onto the orthogonal space of the first type of basis; finally, the basis of the third type of unit was projected onto the orthogonal space of the two types of basis; this process was iterated until the inner product of the F-norm of any two types of basis matrices was less than 10. -5 Convergence is determined at that time.

[0337] Project the basis alignment tensor onto the other two types of heterogeneous basis elements respectively, and extract the cross-element interference projection components. ;

[0338] Subtract all cross-element interference components from the alignment tensor to remove the diffuse aliasing signal introduced by foreign fiber elements:

[0339] ;

[0340] This means removing the interference between the diffusion of different topological units and retaining only the true diffusion characteristics of the current fiber unit.

[0341] In this formula, all operations are 3×3 matrix multiplications. This is the projection matrix corresponding to the heterogeneous basis, and this matrix is ​​related to the alignment tensor. Multiplying them yields the projection component of the tensor onto the basis space of the different elements, which is the cross-element diffuse interference component. Subtracting all the interference components of the different elements from the aligned tensor leaves only the effective component of the tensor in its own basis space, resulting in a pure tensor free of cross-element interference.

[0342] Step 52 includes:

[0343] The pure tensor matrix after removing cross-unit diffusion interference is symmetricized, and the symmetricized tensor is subjected to eigenvalue decomposition and non-negative truncation to obtain the corrected eigenvalue and eigenvector matrix.

[0344] This application also provides a spinal cord diffusion tensor imaging fidelity denoising and reconstruction device, the spinal cord diffusion tensor imaging fidelity denoising and reconstruction device comprising:

[0345] Image acquisition module, which is used to acquire raw DWI data and T2-weighted structural images of the spinal cord region of the subject;

[0346] The mask acquisition module is used to perform pathological spinal cord segmentation based on T2-weighted structural images, and generate preprocessed T2 structural images, spinal cord masks, and gray-white matter boundary masks.

[0347] The DWI image sequence acquisition module is used to perform cross-modal fidelity denoising in the DWI domain of the original DWI data with T2 structural image prior constraints, and generate a DWI image sequence after boundary constraint anisotropic denoising.

[0348] The reconstruction module is used to perform Gibbs ringing removal, eddy current and EPI geometric distortion correction on the DWI image sequence after boundary constraint anisotropic denoising, so as to obtain a corrected and faithfully reconstructed DWI dataset.

[0349] The image group acquisition module is used to acquire a quantitative diffusion parameter image group and a whole spinal cord symmetric positive definite diffusion tensor field based on the corrected DWI dataset.

[0350] This application has the following advantages:

[0351] Traditional SCT segmentation algorithms and conventional level set segmentation both presuppose that the spinal cord cross-section is an elliptical / circular geometric model. When patients with chronic spinal cord injury have swelling, cavitation, atrophy, or flattened cross-sections, the segmentation contour deviates from the real spinal cord boundary. Subsequent denoising and tensor calculations are performed in the wrong region, resulting in systematic parameter distortion. Conventional deep learning segmentation relies on a large number of normal spinal cord labeled samples, but it has poor generalization to pathological samples and cannot be applied to small clinical datasets.

[0352] This application employs Otsu coarse segmentation and unconstrained level set evolution. The energy functional relies solely on image grayscale consistency and edge term iteration, without imposing any cross-sectional geometric constraints. It adapts to the patient's individual pathological morphology based solely on the texture and contrast of the T2 image itself. It eliminates the need for pre-setting a standard spinal cord cross-sectional shape, possessing natural adaptability to various spinal cord pathological malformations. Furthermore, it does not depend on massive pathological annotation training datasets; precise spinal cord and gray-white matter mask extraction can be directly performed from small batches of patient images in primary care settings, avoiding end-to-end computational errors caused by segmentation region offsets from the source.

[0353] In mainstream spinal cord processing tools, segmentation and denoising are separate modules. T2 high-resolution anatomical information is only used for segmentation and is not passed to the DWI denoising stage. The gray-white matter boundary information cannot constrain the DWI smoothing process at all, and the anatomical prior of the segmentation output is wasted.

[0354] This application achieves rigid registration and resolution matching between T2 and b=0, and completely maps the gray and white matter boundary mask generated in step 2 onto all DWI gradient images. The anatomical boundary is used as the core basis for smoothing kernel attenuation constraints, enabling the segmentation result to directly drive denoising calculations. This approach establishes a link between T2 anatomical modality and DWI diffusion modality information, achieving integrated reuse of cross-modal information for segmentation and denoising. It fully utilizes the gray-scale differences between high-resolution T2 gray and white matter to compensate for the boundary blurring defects of low SNR DWI images, significantly improving the utilization rate of prior information.

[0355] Existing technologies such as NLM and MP-PCA all use circular isotropic smoothing cores, with consistent smoothing intensity in both perpendicular and parallel fiber directions. During smoothing, microscopic dispersion differences are smoothed across gray and white matter boundaries, and the FA gradient in the damaged area is eliminated by smoothing. However, simple strong boundary truncation will cause noise residue.

[0356] This application constructs a direction-dependent smoothing kernel, with the smoothing standard deviation of parallel fibers being greater than that of vertical fibers (ratio 2-5). A linear decay weight is set, with the smoothing weight set to zero at the boundary and gradually restored as the weight moves away from the boundary. This approach strengthens smoothing along the long axis of the corticospinal tract, effectively suppressing random noise in low signal-to-noise ratio regions of the spinal cord. The smoothing intensity is weakened along the vertical fiber direction, and combined with an adaptive boundary truncation mechanism, this fundamentally avoids blurring of the gray-white matter boundary, fully preserving the true diffusion differences between gray and white matter, thus balancing noise reduction and tissue contrast.

[0357] Existing noise reduction techniques rely solely on single-layer constraints at the two-dimensional gray-white matter boundaries, failing to distinguish between the three types of topological units: the corticospinal tract trunk, gray matter decussations, and marginal fibers. During smoothing, signals from different fiber bundles mix, and subtle pathological diffusion gradients such as local cavities and mild atrophy are masked by signals from surrounding normal fibers, making it easy to miss mild spinal cord injuries.

[0358] This application employs a complementary topological skeleton extraction, gradient source indexing, and topological bundle isolation mask to construct a topological-spatial dual-locked weight field. Smooth weighting is only allowed within the same topological unit; cross-unit voxel smoothing weights are directly set to zero, isolating longitudinal fiber signals from transverse gray matter signals. By superimposing three-dimensional fiber topological longitudinal constraints on two-dimensional gray-white matter boundary constraints, dual fidelity preservation is achieved, maintaining transverse boundaries and longitudinal fiber pathological gradients. It distinguishes between native diffuse signals from different fiber units and cross-unit interference clutter, avoiding the filtering out of local lesion signals as noise, thus improving the ability to identify small spinal cord injury lesions.

[0359] Existing smoothing techniques rely solely on Euclidean spatial distance to screen similar voxels, ignoring the curve distribution characteristics of long, curved fibers in the spinal cord. Voxels of the same fiber at long distances are judged as unrelated, while voxels of different fibers at close distances are incorrectly included in the smoothing process, introducing systematic calculation bias.

[0360] This application designs a curve coordinate system that maps to the skeletal arc length, and uses an exponential decay source tracing operator to trace the original source location of the signal along the fiber direction. It quantifies and distinguishes between effective signals within the same fiber and interference signals across units, generating a set of mutually exclusive unit indices. It overcomes the inherent limitations of planar straight-line distances and adapts to the anatomical morphology of spinal cord curved fibers. Signal source tracing is completed independently gradient by gradient, and the consistency is stronger when fitting tensors with multiple gradients, reducing the overall offset error of parameters such as MD and RD, and improving the stability of quantitative measurements.

[0361] Existing techniques directly symmetricize and decompose the original tensor after fitting it, resulting in the overlapping of diffuse substrates of different fiber topological units and mutual interference between gray and white matter tensor components. Consequently, the AD / RD parameters cannot truly reflect the degree of damage to a single fiber bundle.

[0362] This application adds a unit attribution labeling, PCA-specific base construction, base alignment, and orthogonal constraint interference removal process before tensor symmetry to eliminate projection interference components from other fiber units, obtaining a pure unit tensor before performing standard tensor solving. This approach can separate the inherent diffusion characteristics of different fiber bundles and eliminate parameter shifts caused by tensor aliasing between units; the pure tensor undergoes non-negative eigenvalue truncation, and the output FA, MD, AD, and RD more closely reflect the local microscopic pathological state of the spinal cord, resulting in higher reliability of quantitative assessment.

[0363] Compared to the technical solutions that use neural networks, this application employs imaging algorithms such as morphological operations, linear algebra, weighted least squares, and level set numerical optimization. It does not involve neural network training or inference and relies solely on the raw T2 and DWI data of a single patient.

[0364] No large-scale labeled training dataset is required; spinal cord images of any single patient can be reconstructed independently. It is compatible with different segments of the cervical and thoracic spinal cord and can also be extended to the analysis of brain CST Waller degeneration DTI, making it applicable to a wider range of scenarios and allowing for direct deployment and use in primary care clinics.

[0365] Most DTI processing tools perform full-image, full-domain operations, with a large amount of irrelevant background such as cerebrospinal fluid, vertebral bodies, and soft tissues participating in segmentation, denoising, and tensor calculations. This results in numerous redundant operations, long processing times, and continuous interference from peripheral magnetic susceptibility artifacts and cerebrospinal fluid pulsation noise in the calculation of spinal cord parenchyma.

[0366] This application outputs a spinal cord mask as a unified computational region constraint throughout the entire process. Segmentation, registration, denoising, artifact correction, and tensor fitting are all performed only within the mask, and background voxels are directly skipped from calculation. This eliminates a large amount of irrelevant background redundant calculations, shortening the overall processing time; it also completely shields the interference of background noise such as cerebrospinal fluid pulsation and vertebral magnetic susceptibility on the spinal cord parenchymal parameters, further improving the signal-to-noise ratio and quantitative stability of the reconstruction results.

[0367] Although the present invention has been described in detail above with general descriptions and specific embodiments, modifications or improvements can be made to it, which will be obvious to those skilled in the art. Therefore, all such modifications or improvements made without departing from the spirit of the present invention fall within the scope of protection claimed by the present invention.

Claims

1. A method for fidelity-preserving denoising and reconstruction of spinal cord diffusion tensor imaging, characterized in that, The spinal cord diffusion tensor imaging fidelity denoising reconstruction method includes: Step 1: Acquire raw DWI data and T2-weighted structural images of the spinal cord region of the subject; Step 2: Perform pathological spinal cord segmentation based on T2-weighted structural images to generate preprocessed T2 structural images, spinal cord masks, and gray-white matter boundary masks; Step 3: Using T2 structured image prior constraints, perform cross-modal fidelity denoising in the DWI domain of the original DWI data to generate a DWI image sequence after boundary constraint anisotropic denoising. Step 4: Perform Gibbs ringing removal, eddy current and EPI geometric distortion correction on the DWI image sequence after boundary constraint anisotropic denoising to obtain a corrected and faithfully reconstructed DWI dataset; Step 5: Obtain the quantitative diffusion parameter image set and the whole spinal cord symmetric positive definite diffusion tensor field based on the corrected DWI dataset.

2. The method for fidelity-preserving denoising and reconstruction of spinal cord diffusion tensor imaging as described in claim 1, characterized in that, Step 2 includes: Step 21: Perform N4 off-field correction and nonlocal mean noise suppression on the acquired T2-weighted structure image to obtain the preprocessed T2-weighted structure image; Step 22: The spinal cord region is extracted using Otsu adaptive threshold segmentation combined with the level set evolution algorithm to generate a spinal cord mask. During the execution of the level set evolution algorithm, no elliptical or circular geometric constraints are applied to the cross-section of the spinal cord. Step 23: Within the area defined by the spinal cord mask, based on the gray-scale difference between the high signal of gray matter and the low signal of white matter in the T2-weighted structural image, the gray-white matter interface of the spinal cord is detected by gradient analysis combined with the region growing algorithm, and a gray-white matter boundary mask is generated.

3. The method for fidelity-preserving denoising and reconstruction of spinal cord diffusion tensor imaging as described in claim 2, characterized in that, Step 3 includes: Step 31: Perform rigid registration and spatial resolution matching between the preprocessed T2-weighted structured image and the b=0 reference image in the original DWI data to obtain the cross-modal spatial mapping relationship; Step 32: Combine the gray-white matter boundary mask and the cross-modal space mapping relationship to construct an anisotropic smoothing kernel, thereby obtaining a voxel-level orientation-dependent smoothing kernel, wherein the smoothing intensity along the fiber direction is greater than the smoothing intensity perpendicular to the fiber direction. Step 33: Set a linear decay weight based on the distance from the voxel to the gray-white matter boundary to obtain a smoothing kernel with boundary constraints, wherein the amplitude of the smoothing kernel at the boundary is zero; Step 34: Iteratively smooth the DWI images in each diffusion gradient direction within the spinal cord mask area until convergence, generating a boundary-constrained anisotropic denoised DWI image sequence.

4. The method for fidelity-preserving denoising and reconstruction of spinal cord diffusion tensor imaging as described in claim 3, characterized in that, Step 4 includes: Step 41: Perform Gibbs ringing removal processing on the DWI image sequence after boundary constraint anisotropic denoising to obtain an intermediate DWI image sequence with ringing artifacts eliminated; Step 42: Perform eddy current correction and EPI geometric distortion correction sequentially on the intermediate DWI image sequence after eliminating ringing artifacts to obtain the corrected and faithfully reconstructed DWI dataset.

5. The method for fidelity-preserving denoising and reconstruction of spinal cord diffusion tensor imaging as described in claim 4, characterized in that, Step 5 includes: Step 51: The original diffusion tensor matrix is ​​obtained by fitting the corrected and faithfully reconstructed DWI dataset using the weighted linear least squares method; Step 52: Symmetricize the original diffusion tensor matrix, and perform eigenvalue decomposition and non-negative truncation on the symmetric tensor to obtain the corrected eigenvalue and eigenvector matrix. Step 53: Reconstruct the whole spinal cord symmetric positive definite diffusion tensor field using the correction eigenvalues ​​and eigenvector matrices; Step 54: Calculate the anisotropy fraction map, mean diffusion coefficient map, axial diffusion coefficient map, and radial diffusion coefficient map based on the whole spinal cord symmetric positive definite diffusion tensor field. The anisotropy fraction map, mean diffusion coefficient map, axial diffusion coefficient map, and radial diffusion coefficient map constitute a quantitative diffusion parameter image group.

6. The method for fidelity-preserving denoising and reconstruction of spinal cord diffusion tensor imaging as described in claim 5, characterized in that, Prior to step 34, step 3 further includes: Based on the spinal cord mask corresponding to the registered T2-weighted structural image, the continuous topological skeleton curve of the spinal cord center is extracted. The topological skeleton is segmented and truncated using the gray and white matter boundary mask as a constraint, dividing it into three types of topological units: the main segment of the corticospinal tract, the gray matter decussation segment, and the terminal segment of the spinal cord marginal fiber. For each type of topological unit, a diffusion gradient source tracing mapping relationship is established. The original diffusion signal source of each voxel is traced back along the topological skeleton curve under different diffusion gradient directions, and a corresponding signal source tracing index set for the topological unit is generated. Based on the signal source index set, a topology bundle shielding rule is established. The diffuse signal between different topology units is only allowed to propagate along the skeleton curve in the positive direction, blocking the interaction of transverse diffuse signals across units and generating a topology bundle isolation mask. A topological-spatial dual-locked weight field is constructed by superimposing gray-white matter boundary space constraints on a topological beam isolation mask as a basis. Smoothing kernel weight superposition calculation is performed only within the same topological unit, and the voxel weights across topological units are independent of each other. Step 34 includes: The anisotropic smoothing kernel is corrected based on the topological-spatial dual-locked weight field, and then global iterative smoothing is performed to generate a DWI image sequence with boundary-constrained anisotropic denoising.

7. The method for fidelity-preserving denoising and reconstruction of spinal cord diffusion tensor imaging as described in claim 6, characterized in that, The process of establishing a diffusion gradient source tracing mapping relationship for each type of topological unit, tracing the original diffusion signal source of each voxel under different diffusion gradient directions along the topological skeleton curve, and generating a corresponding topological unit-specific signal source tracing index set includes: Establish a local curve coordinate system with the central axis of the spinal cord topological skeleton curve. For any voxel x within the spinal cord mask, solve for the shortest vertical projection point p from voxel x to the topological skeleton curve. Extract the skeleton curve arc length parameter s(p) corresponding to the projection point p, and construct the voxel-skeleton arc length mapping function. ; Obtain the preset gradient domain diffusion signal source tracing operator; Traverse all diffusion gradient directions, use the source tracing operator to calculate the source signal amplitude of the corresponding topological unit under each gradient on a voxel-by-voxel basis, and index and bind all voxels in the same topological unit whose source amplitude is higher than the topological signal basis threshold to generate a signal source tracing index set exclusive to that topological unit. Mutual exclusion tags are applied to the signal tracing index sets of different topological units to ensure that there are no overlapping voxels in the index sets of different topological units, and the tracing index sets of all topological units are constructed.

8. The method for fidelity-preserving denoising and reconstruction of spinal cord diffusion tensor imaging as described in claim 7, characterized in that, The preset gradient domain diffuse signal source tracing operator formula is as follows: ； in, The original diffuse signal of voxel x under the i-th diffusion gradient direction; gradient The skeletal arc length corresponding to the origin point of the original signal of the fiber; α is the total arc length of a single topological unit skeleton; α is the topological attenuation coefficient. This is the gradient domain diffuse signal source tracing operator.

9. The method for fidelity-preserving denoising and reconstruction of spinal cord diffusion tensor imaging as described in claim 8, characterized in that, Before performing eigenvalue decomposition and nonnegativity truncation on the symmetricized diffusion tensor matrix, step 5 further includes: Based on the topological bundle isolation mask, all voxels within the spinal cord mask are assigned to corresponding fiber topological units, and the topological unit assignment labels for each voxel are obtained. Using all symmetric diffusion tensors within the same topological unit as a sample set, solve for the basis matrix of the tensor specific to that topological unit, and perform a basis alignment transformation on the symmetric diffusion tensors of each voxel within the unit to obtain the basis aligned tensor. By applying inter-element tensor orthogonality constraints to the basis alignment tensor and removing the interference projection components brought by the tensor basis of other topological elements, a pure tensor with cross-element diffusion interference removed is obtained. Step 52 includes: The pure tensor matrix after removing cross-unit diffusion interference is symmetricized, and the symmetricized tensor is subjected to eigenvalue decomposition and non-negative truncation to obtain the corrected eigenvalue and eigenvector matrix.

10. A fidelity-preserving, noise-reducing reconstruction device for spinal cord diffusion tensor imaging, characterized in that, The spinal cord diffusion tensor imaging fidelity denoising reconstruction device includes: Image acquisition module, which is used to acquire raw DWI data and T2-weighted structural images of the spinal cord region of the subject; The mask acquisition module is used to perform pathological spinal cord segmentation based on T2-weighted structural images, and generate preprocessed T2 structural images, spinal cord masks, and gray-white matter boundary masks. The DWI image sequence acquisition module is used to perform cross-modal fidelity denoising in the DWI domain of the original DWI data with T2 structural image prior constraints, and generate a DWI image sequence after boundary constraint anisotropic denoising. The reconstruction module is used to perform Gibbs ringing removal, eddy current and EPI geometric distortion correction on the DWI image sequence after boundary constraint anisotropic denoising, so as to obtain a corrected and faithfully reconstructed DWI dataset. The image group acquisition module is used to acquire a quantitative diffusion parameter image group and a whole spinal cord symmetric positive definite diffusion tensor field based on the corrected DWI dataset.

Images