Autism risk assessment method and system based on brain image feature fusion

Abstract

The application discloses a self-injury risk assessment method and system based on brain image feature fusion, and relates to the technical field of medical image processing; by extracting the skull contour geometric centroid change information in the functional magnetic resonance sequence, the quantitative characterization of the head movement state is realized; combined with the high-frequency gradient discrete characteristics of the current sequence frame, the initial motion evaluation benchmark is dynamically corrected, and the adaptability of the motion artifact identification process to different image quality states is improved. Compared with the traditional fixed threshold correction or abnormal frame direct elimination mode, the application uses the adaptive exponential decay weight to continuously evaluate the sequence reliability, and uses the weighted interpolation reconstruction for the low reliability sequence frame, so that the artifact suppression and the time sequence continuity are considered. On this basis, by fusing the dynamic brain function connection features and the structural magnetic resonance morphological features, the accuracy and stability of the neural function abnormality evaluation result are improved.

Description

Technical Field

[0001] This invention relates to the field of medical image processing technology, specifically to a method and system for autism risk assessment based on the fusion of brain image features. Background Technology

[0002] Autism Spectrum Disorder (ASD) is a complex neurodevelopmental disorder characterized by impairments in social interaction, restricted interests and behaviors, and neurodevelopmental abnormalities. With the development of medical imaging technology, auxiliary assessment methods based on brain imaging features have gradually become a research focus. Functional magnetic resonance imaging (fMRI) can reflect the functional connectivity between different brain regions by detecting changes in blood oxygenation-dependent signals; structural magnetic resonance imaging (sMRI) can characterize anatomical features such as brain tissue morphology and cortical structure. By fusing fMRI and sMRI imaging features, a comprehensive analysis of the neurological functional status of the target individual can be conducted from two dimensions: brain functional activity and brain structural development, providing a data foundation for auxiliary risk assessment of autism.

[0003] However, during functional magnetic resonance imaging (fMRI) scans of children, the difficulty in maintaining a static position for extended periods and the resulting head movements can lead to positional changes and localized image signal anomalies between different time frames in a continuous fMRI sequence. Current techniques typically identify and correct these anomalies based on pre-defined motion criteria. However, the image signal quality during actual acquisition is affected by factors such as the scanning environment, the patient's condition, and changes in local tissue signals. This results in variations in image stability between different sequence frames. Using fixed evaluation criteria cannot accurately reflect the actual reliability of the current sequence frame, leading to inaccurate assessments of the degree of motion influence. Consequently, this affects the stability of subsequent brain functional connectivity feature extraction and multimodal image fusion analysis results. Summary of the Invention

[0004] To address the shortcomings of existing technologies, this invention provides a method and system for autism risk assessment based on brain imaging feature fusion.

[0005] To achieve the above objectives, the technical solution of the present invention is as follows:

[0006] In a first aspect, the present invention discloses an autism risk assessment method based on brain imaging feature fusion, comprising the following steps:

[0007] Acquire a continuous functional magnetic resonance imaging (fMRI) sequence of a target object under a functional task state, wherein the continuous functional magnetic resonance imaging sequence includes multiple sequence frames arranged according to timestamps;

[0008] Edge features are extracted from the current sequence frame and the reference sequence frame in the continuous functional magnetic resonance imaging (fMRI) sequence to obtain the corresponding geometric centroid of the skull contour, and the spatial displacement vector between the two is calculated based on the corresponding geometric centroid of the skull contour.

[0009] Calculate the voxel gradient distribution of the current sequence frame and extract the real-time high-frequency gradient dispersion of the current sequence frame;

[0010] The preset initial reference velocity is weighted and corrected using the real-time high-frequency gradient discreteness to generate a dynamic reference velocity.

[0011] Based on the spatial displacement vector and the dynamic reference velocity, the adaptive exponential decay weight of the current sequence frame is calculated, and it is determined whether the adaptive exponential decay weight is less than a preset artifact threshold.

[0012] If the judgment result is yes, then within the preset time window, extract adjacent sequence frames whose corresponding adaptive exponential decay weight is not less than the preset artifact threshold as interpolation reference frames, and use the interpolation reference frames to perform weighted interpolation reconstruction on the current sequence frame to generate artifact-free sequence frames.

[0013] Based on the dynamic brain functional connectivity matrix extracted from the continuous functional magnetic resonance imaging (fMRI) sequence updated by the artifact removal sequence frames, the dynamic brain functional connectivity matrix is ​​fused with the structural magnetic resonance morphological features of the target object to generate a neurological function abnormality evaluation index for the target object.

[0014] Secondly, this invention discloses an autism risk assessment system based on brain imaging feature fusion, comprising:

[0015] The image sequence acquisition module is used to acquire a continuous functional magnetic resonance sequence of a target object under a functional task state, wherein the continuous functional magnetic resonance sequence includes multiple sequence frames arranged according to timestamps.

[0016] The head motion feature extraction module is used to extract edge features from the current sequence frame and the reference sequence frame in the continuous functional magnetic resonance imaging sequence, obtain the corresponding geometric centroid of the skull contour, and calculate the spatial displacement vector between the two based on the corresponding geometric centroid of the skull contour.

[0017] The gradient state analysis module is used to calculate the voxel gradient distribution of the current sequence frame and extract the real-time high-frequency gradient dispersion of the current sequence frame.

[0018] The reference speed correction module is used to perform a weighted correction operation on the preset initial reference speed using the real-time high-frequency gradient dispersion to generate a dynamic reference speed.

[0019] The artifact weight evaluation module is used to calculate the adaptive exponential decay weight of the current sequence frame based on the spatial displacement vector and the dynamic reference velocity, and to determine whether the adaptive exponential decay weight is less than a preset artifact threshold.

[0020] The sequence compensation and reconstruction module is used to extract adjacent sequence frames with corresponding adaptive exponential decay weights not less than the preset artifact threshold within a preset time window when the artifact weight evaluation module determines that the current sequence frame is weighted and reconstructed using the interpolation reference frame to generate artifact-free sequence frames.

[0021] The feature fusion evaluation module is used to extract a dynamic brain functional connectivity matrix based on the continuous functional magnetic resonance imaging sequence updated by the artifact removal sequence frame, and to fuse the dynamic brain functional connectivity matrix with the structural magnetic resonance morphological features of the target object to generate a neurological function abnormality evaluation index for the target object.

[0022] Compared with the prior art, the beneficial effects of the present invention are as follows:

[0023] 1. This invention adaptively analyzes the changes in image state in functional magnetic resonance imaging (fMRI) sequences, and can dynamically adjust the motion interference evaluation method according to different acquisition states. This avoids the problem of inaccurate head movement artifact recognition caused by using fixed judgment criteria in the prior art, and improves the adaptability of fMRI image quality evaluation in hyperactive scenarios in children.

[0024] 2. This invention effectively corrects the sequence data affected by motion, thereby reducing the impact of involuntary head movements on the continuity of functional magnetic resonance imaging time series, avoiding the problem of missing brain functional activity information due to the deletion of abnormal data, and improving the completeness and reliability of subsequent brain functional connectivity feature analysis.

[0025] 3. This invention integrates the brain functional activity characteristics reflected by functional magnetic resonance imaging (fMRI) with the brain morphological characteristics reflected by structural magnetic resonance imaging (SMRI), thereby reducing the lack of information in evaluations based on a single imaging type, improving the stability of analysis of abnormal neurological states, and providing data support for autism risk assessment. Attached Figure Description

[0026] The disclosure of this invention is illustrated with reference to the accompanying drawings. It should be understood that the drawings are for illustrative purposes only and are not intended to limit the scope of protection of this invention. In the drawings, the same reference numerals refer to the same parts. Wherein:

[0027] Figure 1 This is a flowchart of the method of the present invention;

[0028] Figure 2 This is a data flow diagram of the present invention;

[0029] Figure 3 This is a system architecture diagram of the present invention. Detailed Implementation

[0030] It is readily understood that, based on the technical solution of this invention, those skilled in the art can propose various interchangeable structural methods and implementations without altering the essential spirit of the invention. Therefore, the following detailed embodiments and accompanying drawings are merely illustrative examples of the technical solution of this invention and should not be considered as the entirety of the invention or as limitations or restrictions on the technical solution of this invention.

[0031] Application Overview

[0032] This invention relates to the field of medical image processing and intelligent assisted assessment technology, and proposes an autism risk assessment method based on brain image feature fusion. This method addresses the problem that involuntary head movements during functional magnetic resonance imaging (fMRI) acquisition in children can easily generate motion artifacts, leading to distortion of brain functional image sequences, inaccurate extraction of functional connectivity features, and reduced reliability of risk assessment. It constructs an adaptive brain image processing and fusion analysis scheme for dynamic motion interference scenarios. Utilizing the image state information of the fMRI sequence itself, the motion evaluation benchmark is dynamically adjusted, enabling the head movement artifact recognition process to adapt to different acquisition states. Finally, by combining the artifact-free functional image features and structural image features, the method achieves the assessment of abnormal neurological functions.

[0033] Existing brain imaging-based methods for autism assessment typically rely on functional magnetic resonance imaging (fMRI) and structural magnetic resonance imaging (SMRI) data for brain functional connectivity analysis and morphological feature extraction. However, during actual data collection on children, the target subject is unlikely to remain still for extended periods, and continuous fMRI sequences are easily affected by head movements. Current techniques often employ fixed motion thresholds or uniform correction methods to identify abnormal image data, which struggles to adapt to changes in image signal quality at different time points, easily leading to errors in motion artifact recognition. Furthermore, simply removing abnormal frames disrupts the continuity of the fMRI time series, impacting dynamic brain functional connectivity analysis. Therefore, improving the adaptive capability of motion interference recognition while enhancing the accuracy of brain image feature analysis while maintaining image data continuity has become a pressing issue in this field.

[0034] To address the aforementioned technical problems, this invention provides a method for autism risk assessment based on brain image feature fusion. This method acquires continuous functional magnetic resonance imaging (fMRI) sequences of the target subject under functional task conditions. By performing edge feature analysis on different time-series frames, spatial displacement information reflecting head movement is obtained. Simultaneously, combined with the image change features within the current sequence frame, the initial motion evaluation parameters are dynamically adjusted to generate a motion evaluation benchmark suitable for the current image state. Based on the dynamically adjusted evaluation parameters, the degree of motion influence on each sequence frame is quantitatively analyzed, and for sequence frames with insufficient reliability, reconstruction processing is performed using neighboring reliable sequence data, thereby obtaining a more continuous, artifact-free fMRI sequence. After image quality optimization, a dynamic brain functional connectivity matrix is ​​further extracted and fused with the structural magnetic resonance morphological features of the target subject to generate an evaluation index for characterizing abnormal neurological functions.

[0035] Through the above technical solutions, this invention can reduce the adverse effects of children's hyperactivity on functional magnetic resonance imaging (fMRI) results, and improve the stability and adaptability of brain imaging data processing in complex motion environments. Simultaneously, by maintaining the integrity of the fMRI time series, it reduces information loss caused by abnormal data processing, improving the reliability of brain functional connectivity feature analysis. Furthermore, by combining multidimensional feature information from functional and structural images, the accuracy of evaluating abnormal neurological states can be improved, providing more stable and effective imaging data support for autism risk assessment.

[0036] like Figure 1 , Figure 2 As shown, this application discloses a method for autism risk assessment based on brain image feature fusion, including the following steps:

[0037] A1. Obtain the continuous functional magnetic resonance imaging (fMRI) sequence of the target object under the functional task state. The continuous functional magnetic resonance imaging (fMRI) sequence includes multiple sequence frames arranged according to timestamps.

[0038] The target population refers to subjects (such as children) who need autism risk screening; autistic individuals often exhibit specific abnormal activation of brain regions when processing social or emotional stimuli; functional task status characterizes the brain's microscopic hemodynamic response to specific external stimuli (such as facial expression recognition, auditory stimulation, etc.). Continuous functional magnetic resonance imaging (fMRI) sequences are core data reflecting the dynamic changes of blood oxygen level dependent (BOLD) signals in the brain over time, and can be directly acquired by 3.0T or higher superconducting magnetic resonance imaging (fMRI) scanners routinely configured in hospitals. Continuous fMRI sequences contain multiple sequence frames arranged with strict timestamps; this temporal order provides fundamental data support for subsequently capturing the dynamic changes of transient neural networks in the brain.

[0039] A2. Extract edge features from the current sequence frame and the reference sequence frame in the continuous functional magnetic resonance imaging sequence to obtain the corresponding geometric centroid of the skull contour, and calculate the spatial displacement vector between the two based on the corresponding geometric centroid of the skull contour.

[0040] Because subjects are prone to head movements during scanning, edge feature extraction techniques were employed to quantify this rigid motion. Specifically, the 3D Canny or Sobel operators were used to perform edge detection on the voxel matrix to remove grayscale interference from internal brain tissue. The extracted geometric centroid of the skull contour represents the absolute rigid position of the entire head in the 3D physical coordinate system. The spatial displacement vector can be obtained by performing spatial difference operations on the centroid coordinates of the current sequence frame and a reference sequence frame (such as the initial frame or the immediately preceding frame). This spatial displacement vector not only indicates the physical direction of head movement, but its magnitude also more intuitively reflects the macroscopic displacement amplitude of the head movement.

[0041] A3. Calculate the voxel gradient distribution of the current sequence frame and extract the real-time high-frequency gradient dispersion of the current sequence frame.

[0042] Voxel gradient distribution is calculated using spatial high-pass filtering and other processing techniques. The voxel gradient distribution characterizes the degree of drastic change in grayscale values ​​between adjacent voxels in an image, reflecting the image's edges and texture details. Real-time high-frequency gradient dispersion is obtained by calculating the statistical variance or the ratio of eigenvalues ​​of the structure tensor from the local gradients. Real-time high-frequency gradient dispersion is used to characterize the degree of microscopic blurring in the image caused by high-frequency micro-twitches of local muscles in the subject (such as jaw tremors due to tension).

[0043] A4. Use the real-time high-frequency gradient dispersion to perform a weighted correction operation on the preset initial reference velocity to generate a dynamic reference velocity.

[0044] The preset initial reference velocity serves as the anchor point for distinguishing normal physiological fluctuations from abnormal head movement artifacts. The preset initial reference velocity can be based on the natural head movement velocity (such as minute displacements caused by breathing and heartbeat) of healthy, age-matched subjects at rest. For example, it can be preset as follows: By introducing real-time high-frequency gradient dispersion and performing nonlinear weighted correction operations (such as logarithmic attenuation modulation), a dynamic reference velocity is generated. When severe high-frequency blurring occurs in the image, the reference velocity is automatically tightened, reducing the tolerance and thus triggering a more stringent artifact removal mechanism, achieving adaptive adjustment to individual scan state fluctuations.

[0045] A5. Based on the spatial displacement vector and the dynamic reference velocity, calculate the adaptive exponential decay weight of the current sequence frame, and determine whether the adaptive exponential decay weight is less than the preset artifact threshold.

[0046] By extracting the magnitude of the spatial displacement vector and processing it in conjunction with the dynamic reference velocity to achieve dimensionlessness, the adaptive exponential decay weight is calculated by substituting it into a negative exponential function with the base of the natural logarithm. This adaptive exponential decay weight directly characterizes the physical reliability of the current sequence frame signal. A preset artifact threshold is the critical line determining whether to retain or discard the current frame. This threshold can be pre-set based on the acceptance of numerous clinical experts and ROC curve sensitivity analysis (e.g., set to 0.65). Through this comparative judgment, it is possible to accurately identify which sequence frames have been severely contaminated by head movements, thereby achieving precise intervention and avoiding the drawback of indiscriminately discarding usable data in traditional methods.

[0047] A6. If the judgment result is yes, then within the preset time window, extract the adjacent sequence frames with corresponding adaptive exponential decay weights not less than the preset artifact threshold as interpolation reference frames, and use the interpolation reference frames to perform weighted interpolation reconstruction on the current sequence frames to generate artifact-free sequence frames.

[0048] Once the current frame is determined to be contaminated, a preset time window (e.g., a 6-second pre- and post-contaminated 3D spatiotemporal neighborhood) is defined, and high-weight (uncontaminated) interpolation reference frames are selected from this window. Subsequently, nonlinear spatial convolution and other interpolation techniques are used to reconstruct the current contaminated frame using weighted interpolation based on the structural information of these high-quality reference frames. The artifact-free sequence is the restored high-quality functional image. This process completely suppresses head movement artifacts while preserving the absolute temporal continuity of the functional magnetic resonance imaging (fMRI) sequence.

[0049] A7. Based on the continuous functional magnetic resonance imaging (fMRI) sequence updated by artifact removal sequence frames, extract the dynamic brain functional connectivity matrix, fuse the dynamic brain functional connectivity matrix with the structural magnetic resonance morphological features of the target object, and generate the neurological function abnormality evaluation index of the target object.

[0050] Pearson correlation coefficients of blood oxygenation signals in each brain region were calculated using the repaired sequences to generate a dynamic brain functional connectivity matrix. This matrix characterizes the intensity of information interaction between different brain regions during the task, i.e., functional abnormalities. Structural magnetic resonance imaging (MRI) morphological features, such as cortical thickness and sulcus fractal dimension, were extracted from additionally acquired T1 / T2-weighted anatomical images, characterizing the physical developmental state of the brain's hardware structure. Through feature fusion techniques such as principal component analysis and Mahalanobis distance calculation, the dynamic functional data and static structural data were concatenated and correlated to ultimately output a neurological functional abnormality evaluation index.

[0051] The Neurological Dysfunction Assessment Index (NDI) is a quantitative value that objectively reflects the deviation of the current brain state from a healthy neurodevelopmental baseline. As a high-dimensional objective reference data obtained through low-level image denoising and deep feature calculations, the NDI will be output to external display devices or medical information systems and provided to clinical medical personnel. Medical personnel can combine this assessment index with other non-imaging information, such as the target individual's clinical behavioral observations (e.g., ADOS scale scores) and family history of genetic diseases, to conduct a more comprehensive and scientific autism risk analysis and early screening.

[0052] Regarding the process of obtaining the spatial displacement vector mentioned in step A2, this application further discloses the following:

[0053] B1. Calculate the three-dimensional voxel gradient maps of the current sequence frame and the reference sequence frame respectively;

[0054] First, the input continuous functional magnetic resonance imaging (fMRI) sequence needs to be decomposed into data along the time dimension. The specific data processing methods are as follows:

[0055] Traverse the continuous functional magnetic resonance imaging (fMRI) sequence according to timestamp order, extract the frame data at the front of the processing queue as the current sequence frame, and extract the adjacent previous frame (or the first frame of the sequence) as the reference sequence frame according to the preset time window pointer.

[0056] After extracting the frame data, three-dimensional voxel gradient maps are calculated for both the extracted current sequence frame and the reference sequence frame. The three-dimensional voxel gradient map is a spatial matrix parameter that characterizes the degree of drastic change in grayscale values ​​between adjacent voxels within a three-dimensional medical image.

[0057] In terms of specific gradient processing methods, a three-dimensional edge detection operator (such as the three-dimensional Sobel operator or the three-dimensional Canny operator) is invoked to perform discrete convolution differentiation operations on the three-dimensional voxel matrices of the current sequence frame and the reference sequence frame in three orthogonal directions: horizontal (X-axis), vertical (Y-axis), and depth (Z-axis), thereby obtaining the corresponding three-dimensional voxel gradient map.

[0058] B2. Binarize the three-dimensional voxel gradient map according to the preset gradient threshold to generate the corresponding skull boundary mask map.

[0059] The preset gradient threshold is an objective benchmark for determining whether a voxel truly belongs to the rigid boundary of the skull. In practice, the preset gradient threshold is set based on the inherent prior gray-level contrast step standard between the skull / dura mater tissue and adjacent cerebrospinal fluid or scalp in standard functional magnetic resonance T2* weighted imaging (e.g., in a normalized gradient field, the preset gradient threshold is empirically set between 0.65 and 0.75). The skull boundary mask generated after binarization is a pure Boolean space matrix, where voxel points at drastically changing boundaries are marked as valid (e.g., set to 1), and the remaining background and internal tissue are marked as invalid (e.g., set to 0). By forcibly discretizing the complex continuous gradient field, a large amount of weak pseudo-edge noise is eliminated, greatly reducing the computational cost of subsequent calculations while locking a complete rigid three-dimensional skull contour boundary composed of valid voxel points.

[0060] B3. Perform a three-dimensional spatial coordinate integration and mean calculation on the voxel points at the boundary in the skull boundary mask image, and use the obtained three-dimensional spatial coordinate mean as the geometric centroid of the skull contour corresponding to the current sequence frame and the reference sequence frame, respectively.

[0061] The geometric centroid of the skull contour represents the spatial geometric center of gravity of the target object's head as a rigid structure in the absolute three-dimensional physical coordinate system of the magnetic resonance scanner. This reduces tens of thousands of complex contour voxels into a single, highly stable macroscopic anchor point. In terms of data processing, all marked valid boundary voxels in the corresponding skull boundary mask are traversed, and their absolute coordinate coefficients along the X, Y, and Z axes are spatially integrated. This integration is then divided by the total number of valid voxels (i.e., the mean of the integral) to obtain the coordinates of the final geometric centroid of the skull contour.

[0062] B4. Perform spatial difference operation based on the three-dimensional coordinates of the geometric centroid of the skull contour in the current sequence frame and the three-dimensional coordinates of the geometric centroid of the skull contour in the reference sequence frame to construct a spatial displacement vector from the reference sequence frame to the current sequence frame.

[0063] The spatial displacement vector characterizes the absolute macroscopic motion state of the target object's head between two sampling timestamps (the current sequence frame and the reference sequence frame). In terms of specific processing, a spatial vector difference operation is performed, with the geometric centroid coordinates of the skull contour in the reference sequence frame as the origin and the geometric centroid coordinates of the skull contour in the current sequence frame as the endpoint (i.e., subtracting the origin coordinates from the endpoint coordinates). The magnitude of the spatial displacement vector quantifies the physical absolute distance of the head movement, while the direction of the vector records the spatial trajectory of the motion.

[0064] Regarding the extraction process of real-time high-frequency gradient discreteness mentioned in step A3, this application further discloses the following:

[0065] C1. Perform three-dimensional multi-scale wavelet decomposition on the voxel matrix of the current sequence frame to extract the corresponding high-frequency detail subband matrix.

[0066] A three-dimensional discrete wavelet transform (3D-DWT) algorithm is used to deeply decompose the original voxel matrix of the current sequence frame. The specific data processing method is as follows: A preset wavelet basis function (e.g., the Daubechies 4 wavelet, which has good orthogonality and compact support) is called, and convolution downsampling operations are performed on the voxel matrix along the three orthogonal coordinate axes (X, Y, and Z axes) using preset low-pass and high-pass filters. Data passing through the three-dimensional low-pass filters is divided into low-frequency approximate components, which represent the overall macroscopic anatomical structure of the brain, with extremely gradual gray-level changes between adjacent voxels. Data passing through at least one-dimensional high-pass filters is divided into high-frequency detail components. Since high-pass filters only allow drastically changing signals to pass through, these components specifically represent microscopic local details with sharp gray-level abrupt changes, such as the edges of brain sulci and gyri, and the contours of microvessels.

[0067] After each single-scale decomposition, one low-frequency approximate component subband and seven high-frequency detail subbands with clear directions are output (including high-frequency subbands in horizontal, vertical, depth, and diagonal directions, such as LLH, LHL, HHL, etc.). These seven high-frequency detail subband matrices are spatially stitched and combined to finally output the high-frequency detail subband matrix. The high-frequency detail subband matrix is ​​a spatial matrix specifically representing the microscopic texture of brain tissue, the edges of complex sulci and gyri, and the boundaries of minute blood vessels.

[0068] C2. Calculate the partial derivatives of each voxel in the high-frequency detail subband matrix in the three orthogonal spatial dimensions, and construct the three-dimensional structure tensor matrix of the corresponding voxel based on the partial derivatives.

[0069] A three-dimensional spatial difference operator (e.g., a three-dimensional central difference template with a preset size of 3×3×3) is used to perform discrete differentiation operations on each voxel in the high-frequency detail subband matrix, thereby calculating the first-order spatial partial derivatives of the voxel in the three orthogonal spatial dimensions of horizontal (X-axis), vertical (Y-axis), and depth (Z-axis) (denoted as ). , , These partial derivative parameters constitute the three-dimensional gradient vector of the current voxel, used to characterize the drastic changes in high-frequency textures in different physical directions.

[0070] However, the gradient vector of a single pixel is highly susceptible to interference from magnetic resonance radio frequency pulse noise. To construct an extremely stable 3D structure tensor matrix that characterizes the local neighborhood structure, the following construction process is further performed:

[0071] First, perform a cross product operation between the current voxel's 3D gradient vector and its transpose to generate a 3×3 initial symmetric cross product matrix (whose internal elements include squared terms of partial derivatives). , , and cross-product terms , , );

[0072] Next, a preset three-dimensional Gaussian smoothing window (e.g., a three-dimensional Gaussian convolution kernel with a preset standard deviation of 1.5) is introduced, and the Gaussian smoothing window is used to perform local weighted convolution integrals on the corresponding elements of all initial symmetric outer product matrices in the current voxel and its preset spatial neighborhood.

[0073] The final 3×3 matrix output after neighborhood smoothing integration is the three-dimensional structure tensor matrix of the voxel.

[0074] The above process perfectly encapsulates the gray-level distribution pattern and edge structure consistency within the local neighborhood (e.g., determining whether the region is flat, has a unidirectional straight edge, or has multiple corners) into a mathematical tensor model, providing underlying data support for subsequent eigenvalue decomposition that is not only highly noise-resistant but also accurately reflects the degree of anisotropic distortion in the image.

[0075] C3. Perform eigenvalue decomposition on the three-dimensional structural tensor matrix to obtain the principal eigenvalue and minimum eigenvalue corresponding to each voxel, and generate local orientation distortion coefficients based on the ratio between the principal eigenvalue and the minimum eigenvalue of each voxel.

[0076] The three-dimensional structure tensor matrix generated in step C2 above is a The real symmetric positive semi-definite matrix is ​​used to reduce its dimensionality by solving its characteristic equation using algebraic analytical methods. The specific data processing methods are as follows:

[0077] Construct the characteristic polynomial equation of the three-dimensional structure tensor matrix (i.e., solve for the determinant). ,in Let be the three-dimensional structure tensor matrix of the current voxel. It is a third-order identity matrix. The three non-negative real eigenvalues ​​corresponding to the voxel matrix are calculated by directly solving the real roots of the cubic equation or by calling the preset Jacobi eigenvalue iteration algorithm for numerical approximation.

[0078] Subsequently, the three non-negative real-valued features extracted earlier are subjected to targeted extraction based on their numerical magnitude. The specific processing method is as follows:

[0079] Call a numerical sorting algorithm (such as quicksort) to sort the three feature values ​​in descending order according to their absolute numerical values, forming a one-dimensional feature value sequence (that is, satisfying...). After sorting, based on the sequence index, the largest value at the beginning of the sequence is defined as the principal eigenvalue of the voxel, and the smallest value at the end of the sequence is defined as the minimum eigenvalue of the voxel. Physically, the principal eigenvalue represents the structural energy of the voxel in the direction of the most drastic grayscale change (i.e., the main edge direction); the minimum eigenvalue represents the structural energy in the direction of the most gradual grayscale change.

[0080] Normal and clear edge voxels have extremely large principal eigenvalues ​​and extremely small minimum eigenvalues. However, when high-frequency micro-flickers occur, the edges are stretched and blurred along the direction of motion, resulting in a sharp decrease in principal eigenvalues ​​and an abnormal increase in minimum eigenvalues. Therefore, local orientation distortion coefficients are calculated based on the ratio between these two values ​​(e.g., normalization by dividing the difference by the sum, or by directly dividing the principal eigenvalue by the minimum eigenvalue, etc., using a pre-defined algorithm). The local orientation distortion coefficient is a quantitative indicator reflecting the deviation of micro-voxel points from their normal anatomical shape due to external motion interference. By comparing the dimensionality of the eigenvalues, the motion blur effect caused by head movement is directly converted into a specific mathematical penalty coefficient, greatly improving the algorithm's detection specificity for microscopic physical deformations.

[0081] C4. Obtain the local orientation distortion coefficients corresponding to all voxels in the high-frequency detail subband matrix, perform logarithmic summation on all local orientation distortion coefficients, and use the summation result as the real-time high-frequency gradient discreteness of the current sequence frame.

[0082] Finally, the local orientation distortion coefficients of the entire image generated in step C3 are traversed and obtained, and these coefficients are accumulated using logarithmic summation. Specifically, each local orientation distortion coefficient is first compressed using a logarithmic mapping with a preset base, and then the entire sequence is summed. The final output real-time high-frequency gradient dispersion is the overall image quality degradation index of the entire sequence frame at the microscopic high-frequency level. During the operation of the magnetic resonance scanner, a very small number of voxels may generate extreme noise distortion values ​​due to radio frequency pulse spikes. Logarithmic mapping has the nonlinear characteristic of high-value compression, which can effectively suppress the linear amplification interference of these isolated extreme noises on the global sum, ensuring that the final generated real-time high-frequency gradient dispersion can truly and stably reflect the true physiological tremor of the subject, thus providing a high-confidence parameter guarantee for the subsequent adaptive artifact suppression mechanism.

[0083] In one embodiment of the present invention, the logarithmic summation operation in step C4 above uses the following calculation formula to perform global nonlinear aggregation processing on all local direction distortion coefficients:

[0084] ;

[0085] in, The real-time high-frequency gradient dispersion of the current sequence frame is characterized by the calculated output. This parameter serves as the final macroscopic judgment index and is used to quantify the severity of high-frequency artifact contamination in the current sequence frame as a whole. Characterizes the total number of effective voxels extracted from the high-frequency detail subband matrix; Characterizing the first The local orientation distortion coefficient corresponding to the individual element, as a low-level physical input variable, records the tissue edge deformation information at the microscopic level. For the natural constant Logarithmic functions with base 0.

[0086] The 1 in the formula represents the preset anti-zero singularity constant. This parameter is necessary because in a perfectly clear and smooth brain region, where no motion blur occurs, the calculated local orientation distortion coefficients... The value may approach zero infinitely. In this case, directly performing logarithmic operations would cause underlying data overflow (i.e., trigger a loop). (The program crashed), adding this constant completely blocked the computational singularities, ensuring the absolute safety of the mathematical domain of the logarithmic function. This is a preset nonlinear scaling penalty coefficient. This parameter is an empirically weighted parameter pre-calibrated based on the fundamental signal-to-noise ratio (SNR) of the MRI machine model. For example, for a 3.0T high-field-strength device, due to its relatively low radio frequency white noise, the nonlinear scaling penalty coefficient can be... A relatively large penalty weight (e.g., 1.5) is preset to amplify extremely weak jitter distortion features; however, for 1.5T low field strength equipment, to avoid baseline noise being misjudged as head motion distortion, the following can be used: The default weight is a low weight (e.g., 0.8).

[0087] The purpose of employing the aforementioned logarithmic summation operation is that, during complex functional magnetic resonance imaging (fMRI) scans, a very small number of voxels may exhibit extreme local directional distortion coefficients due to external electromagnetic spikes or implant metal artifacts. The inherent high-value saturated convergence characteristic of logarithmic mapping—where the larger the independent variable, the more gradually the dependent variable grows—can forcibly suppress these isolated extreme noise parameters, preventing them from having a linear amplification effect during global summation. This results in a more accurate and efficient real-time high-frequency gradient dispersion output. It can accurately and stably reflect the overall high-frequency physiological tremor level of the subjects.

[0088] Regarding the generation process of the local directional distortion coefficients mentioned in step C3, this application further discloses the following:

[0089] C31. Extract the relative difference magnitude between the principal eigenvalue and the minimum eigenvalue, and normalize the relative difference magnitude based on the cumulative feature value of the principal eigenvalue and the minimum eigenvalue to generate the initial anisotropic features.

[0090] Using the principal and minimum eigenvalues ​​obtained from eigenvalue decomposition as basic data input, the data processing first employs a discrete subtraction operator to calculate the absolute difference between the principal and minimum eigenvalues, extracting their relative difference magnitude. Simultaneously, an addition operator is used to calculate the sum of the principal and minimum eigenvalues ​​(or all three eigenvalues ​​of the 3D structure tensor), yielding the accumulated feature quantity (i.e., the total gradient energy of the current voxel). Subsequently, a normalized division operation is performed, dividing the relative difference magnitude by the accumulated feature quantity to output the initial anisotropy feature. This initial anisotropy feature characterizes the purely geometrical selectivity and stretching of the local gradient direction of the current voxel before applying noise penalties.

[0091] C32. Extract the tensor traces of the main diagonal elements of the three-dimensional structure tensor matrix and use them as the local structure energy features of the current voxel.

[0092] Taking the three-dimensional structure tensor matrix constructed in the previous step C2 as input, the specific data processing method reads this... The three elements on the main diagonal of the tensor matrix (i.e., the neighborhood smoothing value of the sum of the squares of the partial derivatives of the three orthogonal dimensions) are summed to extract the tensor trace of the matrix as a local structural energy feature. The extracted local structural energy feature characterizes the overall spatial gradient intensity and texture richness contained in the current voxel and its spatial neighborhood. It will then be used as an objective energy scale to determine whether the current voxel is located at the edge of high-contrast real tissue (high-energy region) or in a smooth, unstructured background / cerebrospinal fluid region (low-energy region).

[0093] C33. A nonlinear weighting mapping function with monotonically increasing and high-value saturation convergence characteristics is used to perform weighting mapping on the local structural energy features to generate the corresponding structural confidence penalty term.

[0094] The local structural energy features extracted in step C32 are used as input variables. A pre-defined nonlinear mapping function (e.g., Sigmoi function, sigmoid saturation function, or exponential saturation mapping function) with monotonically increasing and high-value saturation convergence characteristics is invoked to perform function mapping calculations. The structure confidence penalty term generated after mapping processing represents the physical confidence that the current voxel contains real anatomical structural information. The structure confidence penalty term acts as a dynamic switch for subsequent feature modulation to filter out noise false positives in unstructured background regions.

[0095] In one specific embodiment of the present invention, the structural confidence penalty term is generated using the following nonlinear weight mapping function:

[0096] ;

[0097] in:

[0098] The output includes a structure confidence penalty term;

[0099] The local structural energy features of the current voxel extracted above (i.e., the trace of the three-dimensional structural tensor matrix).

[0100] The global baseline noise variance of the magnetic resonance device in the current scanning environment is obtained by extracting the background signal variance of the magnetic resonance scanner in the background area at the edge of the non-scanning field of view (FOV), which reflects the thermal noise level of the current hardware device in real time.

[0101] This is a preset sensor hardware reference sensitivity coefficient, typically ranging from 1.2 to 5.0.

[0102] The preset nonlinear modulation index (and );

[0103] For the natural constant Exponential operations with base 0.

[0104] The above formula introduces the concept of including The dynamic denominator makes the generation of the structure confidence penalty term no longer an isolated local operation, but directly controlled by global parameters of the current machine hardware operating state.

[0105] Limited (For example, set to 2 or 3), when local energy When the energy is at a low level (i.e., in pure noise background), the exponentiation of the ratio will cause it to rapidly approach 0, thus forming an extremely smooth noise dead zone that completely blocks false alarms. However, once the energy crosses the hardware noise baseline, the mapped value will burst with an extremely steep slope and quickly saturate to 1. The formula itself has the mathematical characteristics of being highly sensitive to weak edges and absolutely suppressing background noise.

[0106] Sensor hardware reference sensitivity coefficient Its specific value is not an arbitrary constant, but needs to be adaptively preset according to the main magnetic field strength and basic signal-to-noise ratio of the functional magnetic resonance scanner used on the target object. The core physical meaning of this parameter is to control the tolerance threshold to background thermal noise.

[0107] In specific application scenarios, the default example standards are as follows:

[0108] For high-field-strength superconducting magnetic resonance imaging (MRI) devices of 3.0T and above: Due to the extremely high signal-to-noise ratio of high-field-strength devices, the background random thermal noise generated by the RF coils... It is inherently very weak. To prevent the algorithm from being overly sensitive and misinterpreting extremely weak normal blood flow pulsations as tissue edges, it will... The default is a relatively small tightening value, for example, the default value. At this point, it is possible to keenly capture the true microscopic edges, while utilizing... Even a slight scaling can effectively suppress the already weak hardware thermal noise.

[0109] For 1.5T conventional field strength MRI equipment scenarios: Due to the relatively low signal-to-noise ratio of conventional field strength equipment, the overall image is often enveloped in a thick layer of salt-and-pepper noise or white noise, resulting in... The absolute value is too large. If If the value is too small, a large amount of high-energy noise will exceed the penalty threshold. Therefore, in such hardware scenarios, [the following will be considered:] The default is a larger tolerance penalty value, for example, the default value is... even By increasing The value artificially raises the threshold for judging the edge of real tissue, ensuring that only those brain gyri physical structures with extremely strong contrast (such as the gradient energy of the gray-white matter interface) can break through the suppression of the denominator, thus maintaining extremely high artifact removal specificity even under low signal-to-noise ratio hardware conditions.

[0110] C34. Perform feature modulation operation on the initial anisotropic features according to the structural confidence penalty term, and use the modulated output as the local orientation distortion coefficient.

[0111] The structural confidence penalty term generated in step C33 and the initial anisotropic features obtained in step C31 are used as two-dimensional input streams. In terms of specific data processing, a dot product feature modulation operation is performed, that is, the structural confidence penalty term is used as a weighting coefficient to multiply the initial anisotropic features. The modulated final output local orientation distortion coefficient rigorously characterizes the degree of true net orientation distortion of the current voxel caused by the subject's high-frequency head movement after excluding background noise interference. The local orientation distortion coefficient, as the final quantization index at the micro-voxel level, is directly used in the global logarithmic summation step in step C4 to synthesize the frame-level real-time high-frequency gradient dispersion.

[0112] Regarding the process of generating the dynamic reference velocity mentioned in step A4, this application further discloses the following:

[0113] D1. Obtain the global baseline gradient dispersion of the target object in the resting state;

[0114] Global baseline gradient dispersion is a scalar parameter characterizing the background level of the combination of the inherent microtexture of the brain and microphysiological noise (such as spontaneous vascular pulsation and weak respiratory elevation) of a target object in a calm state (resting state) without specific brain function task load. Global baseline gradient dispersion can be directly obtained by acquiring an additional resting-state fMRI sequence (e.g., a resting-state fMRI sequence with a duration of 1 to 2 minutes) before performing the formal task fMRI scan; the high-frequency gradient dispersion of each frame is extracted by applying the aforementioned wavelet decomposition and structural tensor calculation methods, and a time-domain arithmetic mean is performed to finally generate the global baseline gradient dispersion.

[0115] D2. Extract the dynamic deviation characteristics of the real-time high-frequency gradient dispersion relative to the global benchmark gradient dispersion, and generate the dispersion variation index.

[0116] The real-time high-frequency gradient dispersion obtained in step A3 and the global baseline gradient dispersion obtained in step D1 are used as dual-source data inputs. In terms of specific data processing methods, firstly, a difference operator is used to calculate the relative residual between the real-time dispersion and the global baseline dispersion, and the relative residual feature is extracted. Next, a preset physiological baseline fluctuation tolerance feature is introduced (this tolerance feature is pre-calibrated based on the 95% confidence interval of the resting-state dispersion of a large number of healthy children, i.e., 1.96 times the standard deviation; for example, the preset upper limit of the tolerance envelope is 0.12). The relative residual feature is compared with this upper limit distribution envelope through feature mapping, and the dispersion variation index is obtained through mapping as a dimensionless control quantity that purely characterizes the intensity of abnormal head movement distortion.

[0117] D3. The dispersion variation index is mapped using a preset dynamic adjustment function with asymmetric decay characteristics to generate an adaptive penalty weight. When the dispersion variation index indicates a positive deviation, the adaptive penalty weight exhibits a nonlinear accelerated decay distribution.

[0118] The dispersion variation index generated in the preceding step D2 is used as the independent variable and input into a preset dynamic adjustment function with asymmetric decay characteristics to perform nonlinear mapping calculations; specifically, an example formula for this dynamic adjustment function can be constructed as follows:

[0119] ;

[0120] in, The adaptive penalty weights characterize the output, and their values ​​are distributed in the interval [range]. The dimensionless scaling factor within is used to characterize the compression ratio of the reference velocity; The input is the dispersion variability index; and These are preset function control parameters.

[0121] The decay rate coefficient is preset based on the image signal-to-noise ratio (SNR) of the magnetic resonance imaging (MRI) device (e.g., preset for a 3.0T device). ); It is a non-linear acceleration exponent, and is forcibly preset to a constant greater than 1 (such as the preset value). ).

[0122] When the dispersion variation index At that time, the function outputs adaptive penalty weights. (i.e., no punishment is imposed); and once (Indicates positive deviation, i.e., abnormal head motion blur), due to The function, As the variation exponent increases, a nonlinear acceleration decay curve emerges, characterized by a gradual initial decrease followed by a steeper decrease. This asymmetric acceleration decay function provides the algorithm with a certain tolerance margin when the microscopic image is only slightly blurred. However, once the high-frequency blur intensifies, the judgment threshold is rapidly tightened with an exponentially decreasing weight, achieving a highly sensitive and severe penalty for serious artifacts.

[0123] D4. Apply adaptive penalty weights to the preset initial reference velocity and perform dynamic attenuation modulation, using the modulated result as the dynamic reference velocity.

[0124] The preset initial reference velocity is a rigid initiation threshold for judging head movements at the macroscopic kinematic level. Its preset standard is based on the baseline of the average normal physiological head movement velocity in healthy children of the same age during functional magnetic resonance imaging (fMRI) scans (e.g., for children, the initial reference velocity is empirically preset to...). In terms of data processing methods, a dot product modulation operation is performed, which involves applying the adaptive penalty weights generated in step D3 above. Compared with the preset initial reference speed Multiplication (i.e.) The product of the modulated outputs is the final dynamic reference speed. This serves as the dynamic denominator for calculating the adaptive exponential decay weight in subsequent steps.

[0125] When the microtexture quality is good (the adaptive penalty weight value approaches or equals 1), the original quality is maintained. Forgiveness speed threshold;

[0126] When the microscopic high-frequency texture is severely blurred (the adaptive penalty weight value approaches 0), the speed threshold is automatically compressed to a lower value (e.g., This makes the subsequent artifact removal mechanism extremely stringent.

[0127] Regarding the generation process of the dispersion variation index mentioned in step D2, this application further discloses the following:

[0128] D21. Extract the relative residual feature between the real-time high-frequency gradient dispersion and the global baseline gradient dispersion, and obtain the preset physiological baseline fluctuation tolerance feature; wherein, the physiological baseline fluctuation tolerance feature is used to characterize the normal physiological rhythm noise distribution boundary of the target object in the resting state.

[0129] The real-time high-frequency gradient dispersion obtained in step A3 and the global baseline gradient dispersion obtained in step D1 are used as inputs. In terms of specific data processing, a difference operator is used to subtract the global baseline gradient dispersion from the real-time high-frequency gradient dispersion, thereby calculating and extracting the relative residual feature between the two. The relative residual feature is a fundamental physical difference that characterizes the absolute shift in the micro-texture blurring of the current sequence frame relative to the inherent resting baseline of the target object.

[0130] Simultaneously, a pre-defined physiological baseline fluctuation tolerance feature is obtained. This feature is specifically designed to characterize the boundary of normal physiological rhythm noise distribution in the target subject at rest. The pre-defined standard for the physiological baseline fluctuation tolerance feature is determined based on the statistical distribution of dispersion of healthy subjects in the same age group under large-scale resting-state fMRI scans. For example, the standard deviation of the dispersion fluctuation of healthy samples is extracted, and the upper limit of the distribution envelope of this tolerance feature is empirically pre-defined as 1.96 times the standard deviation (i.e., the upper limit of the 95% confidence interval), serving as an objective standard to distinguish between normal physiological fluctuations and abnormal head tremors.

[0131] D22. Compare the relative residual feature quantity with the physiological baseline fluctuation tolerance feature through feature mapping to determine whether the relative residual feature quantity exceeds the upper limit distribution envelope of the physiological baseline fluctuation tolerance feature.

[0132] The relative residual characteristic calculated in the previous step D21 is compared with the upper limit distribution envelope of the physiological baseline fluctuation tolerance characteristic to determine the feature size. Based on the different determination results, the subsequent two sets of branch processing are executed.

[0133] D231. If the judgment result is yes, then extract the effective distortion component whose relative residual feature exceeds the upper limit distribution envelope, and perform logarithmic compression mapping on the effective distortion component. Output the feature after logarithmic compression mapping as the dispersion variation index.

[0134] When the relative residual characteristic exceeds the upper limit envelope, it indicates that the current subject has indeed experienced abnormal high-frequency tremors that exceed the normal physiological extreme.

[0135] In the data processing of this branch, a truncation subtraction operation is first performed, that is, the upper limit distribution envelope value is subtracted from the relative residual feature value to strip and extract the effective distortion component. This component represents the net distortion degree caused purely by abnormal motion. Next, in order to prevent sudden electromagnetic spike pulses from causing the value of this component to exceed the limit, the logarithmic function processor is called to perform logarithmic compression mapping processing on the effective distortion component. The dispersion variability index output after processing is the dimensionless control parameter that finally characterizes the intensity of abnormal head motion. By utilizing the high-value gradual increase characteristic of the logarithmic function, the linear amplification effect of extreme outlier noise is effectively suppressed, ensuring that the variability index can still output a stable and safe control variable to the subsequent weight modulation module when faced with extremely severe local image distortion.

[0136] D232. If the judgment result is negative, the preset distortion-free baseline constant will be output as the dispersion variation index.

[0137] When the relative residual feature does not exceed the upper limit envelope, it indicates that the current sequence frame has not been contaminated by abnormal motion.

[0138] In this step, a preset distortion-free baseline constant is directly invoked as the current dispersion variability index for forced output. The preset standard for this parameter is typically set to a neutral value (e.g., a constant of 0) based on the mathematical interface of the subsequent dynamic adjustment function. By outputting a zero value (or a specific neutral value), a complete nullification masking operation is performed on the current small residuals. When the subsequent dynamic adjustment function receives this distortion-free baseline constant, it maintains its original tolerance speed threshold without any penalty attenuation, achieving seamless filtering of normal physiological noise and logical self-consistency of the underlying data link.

[0139] The present invention further discloses the process of performing logarithmic compression mapping on the effective distortion components in step D231; specifically, the effective distortion components are subjected to logarithmic compression mapping by calling the following specific mathematical formula:

[0140] ;

[0141] in, The dispersion variability index of the output is calculated and directly output as a dimensionless control parameter to the subsequent dynamic adjustment function; The effective distortion component of the relative residual feature quantity extracted above exceeds the upper limit distribution envelope, which is the underlying physical input variable.

[0142] This parameter represents the standard deviation of the baseline dispersion fluctuation of the target object in a resting state. It directly uses this inherent parameter, which was used in the previous step to construct the tolerance feature, as an adaptive denominator to transform the effective distortion component into a normalized relative variability rate relative to the subject's own baseline noise level. The introduction of this parameter ensures that the compression algorithm can automatically adapt to the individual brain tissue contrast differences among different subjects.

[0143] The preset acceleration shape parameter is calibrated based on the nonlinear aggravation characteristics of abnormal head movements in disrupting functional connectivity signals in the initial stage. To improve sensitivity to early micro-tremors, it is forcibly preset to a constant slightly greater than 1 (e.g., preset...). ).

[0144] This is the preset structural compression penalty coefficient. The preset standard for this parameter depends on the maximum dynamic quantization range of the magnetic resonance imaging equipment in the current sequence. For example, if the image acquisition sequence is prone to producing high-amplitude transient motion artifact spikes, it is preset to a larger value (e.g., preset...). ), to enhance the compression curvature of the logarithmic function.

[0145] The above formula exhibits the characteristics of sensitive amplification at low values ​​and strict clamping at high values. Internally, through... and Normalization of division and The exponential acceleration effect allows for rapid nonlinear amplification of the effective distortion component just as it crosses the envelope (i.e., when an extremely weak but definite abnormal head movement occurs), ensuring extremely high artifact detection sensitivity. Conversely, when external interference (such as metal implant artifacts or electromagnetic spikes) causes... When extreme, off-the-charts values ​​occur, the outer parameterized logarithmic function... It will leverage its convergence properties to forcibly compress the out-of-the-charts values ​​into a safe and stable range.

[0146] The calculation process for the adaptive exponential decay weight mentioned in step A5 is further disclosed in this application as follows:

[0147] E1. Extract the spatial magnitude of the previously obtained spatial displacement vector and obtain the resampling time interval between adjacent sequence frames of the continuous functional magnetic resonance sequence; and perform dimensional transformation on the dynamic reference velocity based on the resampling time interval to generate a dynamic reference displacement threshold.

[0148] The system receives the previously constructed spatial displacement vector as input data and calls the Euclidean distance solver to perform L2 norm calculation (i.e., calculate the square root of the sum of squares of the coordinate components of the vector), thereby stripping away the directional attributes and extracting the spatial magnitude of the spatial displacement vector. The spatial magnitude of the spatial displacement vector represents the absolute physical distance of the target object's head during rigid motion between two consecutive frames.

[0149] Simultaneously, by reading the DICOM low-level file header information of the continuous functional magnetic resonance imaging (fMRI) sequence, the preset resampling time interval for that fMRI sequence during MRI scanning is parsed and obtained. The resampling time interval is an objectively existing physical scalar (e.g., typically 10 ... or The absolute time difference represents the time taken by the scanner to acquire two adjacent complete brain voxel matrices.

[0150] Then the multiplication operator is called to apply the dynamic reference velocity (whose dimensions are length / time, e.g.) ) and the repeated sampling time interval obtained above (its dimension is time, such as The multiplication operation is performed. Through multiplication, the time dimensions are successfully canceled out, and the output dynamic reference displacement threshold is perfectly converted into an absolute length in terms of physical dimensions (e.g., ...). The dynamic reference displacement threshold characterizes the maximum tolerable macroscopic rigid displacement physical limit within the time interval between two adjacent frames, under the constraint of the current subject's microscopic high-frequency tremor state.

[0151] E2. Calculate the ratio of the spatial modulus to the dynamic reference displacement threshold to obtain the dimensionless motion parameters;

[0152] The extracted spatial modulus is used as the numerator, and the converted dynamic reference displacement threshold is used as the denominator for division. Since both the numerator and denominator are unified to physical length dimensions, this division operation outputs a purely dimensionless motion parameter. This dimensionless motion parameter intuitively represents the rate at which the actual macroscopic displacement amplitude of the current sequence frame exceeds the system's dynamic tolerance threshold, providing a safe, uniform, and unaffected input variable for subsequent pure nonlinear mathematical exponential operations (i.e., calculating adaptive exponential decay weights).

[0153] E3. Using the natural constant as the base and the negative number of the dimensionless motion parameter as the exponential variable, the adaptive exponential decay weight of the current sequence frame is calculated.

[0154] The floating-point exponentiation function in the underlying mathematical function library is used to perform a nonlinear continuous mapping on the aforementioned dimensionless motion parameters. The specific data processing method involves using the inherent natural constants in mathematics (i.e., Euler's number). Its underlying preset precision is usually 100%. (and its floating-point extension) are used as a fixed basis, and the dimensionless motion parameters of the input are negativeed (i.e., multiplied by) It is forcibly converted to a negative number field, and then the negative number is used as an exponent variable to perform floating-point exponentiation. Its underlying core calculation logic can be represented by the following specific mathematical formula:

[0155] ;

[0156] In this formula, Characterizes the adaptive exponential decay weights of the calculated output; The base of the natural constant; The input is a dimensionless motion parameter.

[0157] The adaptive exponential decay weight generated after this step has a value strictly distributed in... A continuous floating-point control parameter within a closed interval. This adaptive exponential decay weight directly characterizes the physical reliability of the magnetic resonance signal in the current sequence frame or the signal confidence level unaffected by head motion artifacts. The adaptive exponential decay weight directly serves as the sole core decision factor in subsequent logic branches for determining whether to trigger the three-dimensional spatiotemporal interpolation reconstruction mechanism, and after triggering repair, it directly acts as the physical weight coefficient in the weighted interpolation operation.

[0158] The aforementioned negative exponential mapping perfectly matches the nonlinear degradation physical law of motion-induced time-dependent oxygenation (BOLD) signals in magnetic resonance imaging.

[0159] When the dimensionless motion parameters approach 0 (i.e., the actual spatial displacement is much smaller than the current tolerance threshold), the adaptive exponential decay weight is kept at a high level close to 1 very smoothly, ensuring that the high-quality uncontaminated image will not be damaged by any computational power.

[0160] When the dimensionless motion parameters gradually increase (i.e., when obvious head motion artifacts that exceed the baseline threshold occur), the negative exponential function will apply a non-linear accelerated decay penalty to the adaptive exponential decay weight in a "steep first, then slow" manner, causing the adaptive exponential decay weight of severely polluted frames to quickly drop to a penalty range close to 0.

[0161] This continuous and smooth nonlinear decay mechanism completely replaces the hard cut-off judgment of non-zero or 1 that deletes images greater than a threshold in traditional medical image processing. It ensures that the data weight transition between frames is absolutely smooth in the subsequent spatiotemporal interpolation reconstruction process, and enhances the temporal continuity and data fidelity of the underlying blood oxygenation signal in autism risk assessment.

[0162] Regarding the process of generating the artifact-free sequence frames mentioned in step A6, this application further discloses the following:

[0163] F1. Extract the adaptive exponential decay weights corresponding to the interpolation reference frame, and perform spatial normalization on the extracted adaptive exponential decay weights to generate spatial smoothing kernel weights for the current sequence frame.

[0164] The algorithm calls upon the adjacent, uncontaminated sequence frames selected in the preceding decision logic (i.e., interpolation reference frames, whose adaptive exponential decay weights are not less than a preset artifact threshold) and their respective adaptive exponential decay weights, and performs arithmetic normalization on this set of weight data. Specifically, it calculates the sum of the decay weights of all selected interpolation reference frames within the current time window, and then divides the decay weight of each interpolation reference frame by this sum. The spatial smoothing kernel weights output after division normalization are a set of relative contribution coefficients whose values ​​are strictly distributed between 0 and 1, and whose sum of weights of all reference frames within the current time window is strictly equal to 1.

[0165] Obtaining the spatial smoothing kernel weights transforms absolute physical reliability into a relative signal allocation ratio during subsequent image fusion and reconstruction. Forced normalization ensures that the total signal energy of each voxel remains physically conserved in the reconstructed image. This completely avoids overall image grayscale anomalies (such as abnormal brightness or local oversaturation) caused by direct weighted summation of multiple frames of data, fundamentally preventing the introduction of artificial, systematic baseline drift into the extremely weak blood oxygen level-dependent (BOLD) signal during the restoration process.

[0166] F2. Using spatial smoothing kernel weights, nonlinear convolution interpolation is performed on the voxel matrix of the current sequence frame and the voxel matrix of the interpolation reference frame at the corresponding voxel positions to generate artifact-free sequence frames.

[0167] The spatial smoothing kernel weights generated in the preceding F1 step are used as dynamic convolution kernel parameters. The voxel matrix of the current sequence frame and the voxel matrix of the multi-frame interpolation reference frame are used as the input for deep fusion of the underlying multidimensional image. Specifically, a three-dimensional voxel-level linear weighted interpolation algorithm is employed for data processing. The specific data processing methods are as follows:

[0168] Using the voxel matrix of the current contaminated sequence frame as a spatial coordinate template, for each three-dimensional spatial coordinate in the matrix... By calling pre-extracted multi-frame interpolation reference frames, a weighted summation calculation is performed according to the following standard:

[0169] ;

[0170] In this calculation formula, The characterization of the artifact-free sequence frame output by computation in coordinates Repair grayscale values ​​at the location; Characterizing the first The true grayscale values ​​of voxels at the same physical coordinates of each interpolation reference frame; Characterizing the generation of the preceding step F1, corresponding to the first step Spatial smoothing kernel weights for each interpolation reference frame; Characterizes the total number of interpolation reference frames extracted within a preset time window (this total number is usually preset according to the time window's tolerance, for example, a preset...). or frame).

[0171] At the underlying execution logic, due to severe head-movement contamination of the current sequence frame, the preservation of the original distorted pixel values ​​of the current frame is directly abandoned. Instead, targeted replacement is performed entirely using a high-quality reference frame that is temporally adjacent. Since the preceding steps have already forcibly guaranteed... (i.e., weighted normalization) This purely linear weighted summation ensures that the reconstructed voxel grayscale values ​​strictly fall within the grayscale envelope of surrounding high-quality frames, absolutely preventing secondary artifacts such as abnormal brightness or local oversaturation caused by mathematical function overshoot. Through this efficient voxel-level fusion, grayscale distortion caused by high-frequency micro-vibrations is completely covered, while perfectly restoring the smooth continuity of the functional magnetic resonance imaging (fMRI) sequence on the time axis with extremely low computational latency, providing pure underlying data support for subsequent extraction of dynamic brain network features.

[0172] Regarding the extraction process of the dynamic brain functional connectivity matrix mentioned in step A7, this application further discloses the following:

[0173] G1. Map the updated continuous functional magnetic resonance imaging (fMRI) sequence to a preset brain region topology map and extract the time-dependent blood oxygenation level signal sequence for each target brain region.

[0174] Pre-defined brain region topology atlases are pre-constructed template data used to deconstruct the three-dimensional continuous anatomical space of the brain into specific functional or anatomical regions of interest (ROIs). The preset standards for these brain region topology atlases are pre-defined based on the consistency of brain function or anatomical boundaries in a large sample of resting-state functional magnetic resonance imaging (fMRI) data from healthy individuals. For example, presets may use the Schaefer 200 functional area atlas containing 200 functional areas of the cerebral cortex, or the AAL (Automated Anatomical Labeling) anatomical atlas containing 116 anatomical brain regions.

[0175] In terms of specific processing methods, the updated continuous functional magnetic resonance imaging (fMRI) sequences are first subjected to spatial alignment and registration transformations to convert their voxel coordinate system to a standard space (such as MNI space). Then, a mask based on a pre-defined brain region topology is used to spatially segment and overlay the mapped sequences. Subsequently, for each segmented target brain region, all voxels contained within the mask of that region are traversed, and the grayscale values ​​of the voxels are averaged at each sampling timestamp to extract the time-dependent oxygenation level signal sequence for that target brain region. This time-dependent oxygenation level signal sequence rigorously characterizes the dynamic fluctuations in neuronal metabolism and oxygen consumption in the corresponding target brain region throughout the entire scanning time axis.

[0176] G2. By traversing the time-dependent blood oxygenation signal sequence through a sliding time window mechanism, the Pearson correlation coefficient of each target brain region in different time windows is calculated. The Pearson correlation coefficients in each time window are aggregated to generate a dynamic brain functional connectivity matrix.

[0177] Set a window with a fixed time length on the time axis, and traverse and extract the time-dependent blood oxygenation signal sequence along the time progression direction with a fixed step size.

[0178] Based on the dominant low-frequency fluctuation frequency of brain blood oxygenation level-dependent signals (usually concentrated in...) The window length and sliding step size can be preset. For example, the window length is empirically preset to 30 repeated sampling time intervals (about 60 seconds), and the sliding step size is preset to 1 sampling time interval (about 2 seconds). The preset window length of 60 seconds can ensure that there are enough signal periods in the window to accurately estimate the correlation, while avoiding the transient dynamic information being overwhelmed by the smoothing due to the window being too long.

[0179] Within each extracted time window, time-dependent blood oxygenation level signal subsequences for all target brain regions within that window are extracted. A binary statistical correlation operator is then used to calculate the Pearson correlation coefficient between each pair of target brain region signal sequences. The covariance of the two brain region signals is divided by the product of their respective standard deviations, outputting a value distributed across... Scalar values ​​within the interval. The Pearson correlation coefficient within each time window is used to quantify the temporal synchronization coupling strength of neural activity between any two brain regions within a specific time window.

[0180] Finally, collect the Pearson correlation coefficient matrices for all time windows calculated using the aforementioned iterative method (assuming there are a total of...). For each target brain region, a calculation is performed for each time window. The symmetric correlation matrix, the entire sequence is extracted (a time window).

[0181] In terms of specific data processing methods, time-domain tensor concatenation and aggregation operations are performed, that is... The size is The double matrices are stacked along the third dimension (time dimension) according to the order of the time windows, forming a matrix with dimension 1. The three-dimensional tensor is aggregated into a dynamic brain functional connectivity matrix that characterizes the dynamic properties of connectivity. The final output dynamic brain functional connectivity matrix records and characterizes the trajectory of the evolution of functional connectivity strength between target brain regions over time and the frequency of state switching, serving as a high-dimensional functional feature input for subsequent feature fusion and risk assessment.

[0182] Regarding the process of obtaining the structural magnetic resonance morphological features mentioned in step A7, this application further discloses the following:

[0183] G3. Acquire high-resolution structural magnetic resonance images of the target object, and perform three-dimensional surface mesh reconstruction on the high-resolution structural magnetic resonance images to generate inner and outer boundary meshes of the cortex.

[0184] High-resolution structural magnetic resonance imaging (MRI) is anatomical image data characterizing the three-dimensional, high-contrast anatomical structure of brain tissue. High-resolution structural MRI images can be acquired by directly scanning the head of the target subject using 3.0T or higher MRI equipment and a 3D T1-weighted magnetic resonance imaging (MRMI) sequence. The preset scanning parameters typically require an isotropic voxel resolution of [insert value here]. (i.e., the volume of a single voxel is) This ensures that the fine folds of the microscopic cortical structure do not produce a partial voxel effect due to excessively large voxels.

[0185] In terms of data processing, the original high-resolution structural images were first subjected to intensity non-uniformity correction (e.g., using the N4 bias field correction algorithm) and skull dissection to remove interference from non-brain tissues. Subsequently, a three-dimensional tissue segmentation algorithm (e.g., a Bayesian segmentation operator based on Markov random fields and prior brain atlases) was used to divide the brain anatomical structure into gray matter, white matter, and cerebrospinal fluid. Based on this, a three-dimensional deformable surface mesh reconstruction method (e.g., the Moving Cubes algorithm combined with topology correction and the surface deformation relaxation algorithm) was used to continuously deform and expand outward from the white matter / gray matter interface as the initial interface, thereby constructing the inner cortical boundary mesh (i.e., the white matter surface mesh, representing the gray-white matter interface) and the outer cortical boundary mesh (i.e., the pia mater surface mesh, representing the gray matter-cerebrospinal fluid interface) composed of a large number of three-dimensional triangular facets and vertices.

[0186] G4. Calculate the shortest Euclidean distance between corresponding vertices of the inner and outer boundary grids of the cortex, and generate the cortical thickness distribution matrix.

[0187] The cortical intracortical boundary mesh and the cortical outer cortical boundary mesh generated in step G3 above are used as the dual-mesh geometric input. In terms of specific data processing, every three-dimensional vertex on the cortical intracortical boundary mesh is traversed. A ray is emitted outward along the surface normal at the vertex, and the corresponding intersecting vertex or nearest projection point is searched on the outer boundary mesh of the cortex. Subsequently, the L2 norm operator is called to calculate the vector difference magnitude between two corresponding vertices, i.e., to solve for the shortest Euclidean distance in 3D space:

[0188] ;

[0189] After traversing and solving all vertices of the entire mesh (typically containing tens of thousands to hundreds of thousands of vertices), the calculated Euclidean distance values ​​of each vertex are arranged in spatial index order and output as a cortical thickness distribution matrix. The cortical thickness distribution matrix represents the spatial distribution of the local physical thickness of gray matter in the entire cortical surface.

[0190] G5. The fractal dimension of gyri folds is extracted based on the inner and outer boundary grid of the cortex, and the cortical thickness distribution matrix and the fractal dimension of gyri folds are combined to form structural magnetic resonance morphological features.

[0191] In addition to microscopic thickness indicators, geometric parameters characterizing the macroscopic complexity of the cerebral cortex's folding are further extracted. In terms of specific data processing methods, non-Euclidean geometric analysis is introduced, employing a three-dimensional box-counting algorithm to measure the spatial folding complexity of the cortical inner and outer boundary grids. Specifically, a series of preset parameters with different side lengths are used. A three-dimensional equidistant grid cube (e.g., a cube with a preset side length) Take values ​​in sequence }) Overlay the reconstructed cortical surface mesh and count the values ​​at each edge length The total number of 3D cubes containing curved mesh patches Subsequently, the data pairs The least squares linear regression was performed to fit the data, and the slope of the fitted line was extracted and output as the fractal dimension (FD) of the brain folds.

[0192] The fractal dimension of gyral folds is a dimensionless geometric scalar parameter that rigorously characterizes the density and nonlinear geometric complexity of the folding and embedding of the cerebral cortex sulci and gyri in three-dimensional space. The purpose of obtaining this parameter is that abnormal neural network development in the brains of children with autism is often accompanied by phenotypes of hyper-gyrification or under-gyrification of the cerebral cortex.

[0193] After the above extraction is completed, the feature stitching operator is invoked to vector-combine and stack the cortical thickness distribution matrix containing local detailed features and the gyral fold fractal dimension containing global geometric complexity along the feature dimension, ultimately outputting multi-scale structural magnetic resonance morphological features. This feature combination achieves complementary advantages and multi-scale deep fusion of microscopic cortical thickness and macroscopic fold morphology, constructing a high-dimensional morphological baseline characterizing the developmental state of the brain's anatomical hardware. This provides comprehensive and interference-resistant structural imaging feature support for subsequent multimodal joint stitching with the aforementioned dynamic brain functional connectivity matrix, and further assessment of autism risk.

[0194] Regarding the process of obtaining the neurological dysfunction assessment index mentioned in step A7, this application further discloses the following:

[0195] G6. The dynamic brain function connectivity matrix and structural magnetic resonance morphological features are concatenated into a multimodal joint feature vector, and the Mahalanobis distance between the multimodal joint feature vector and the preset health benchmark feature spatial distribution center is calculated.

[0196] The dynamic brain functional connectivity matrix containing time-dependent evolutionary information extracted in step G2 above (i.e. The dual-source input data consists of a tensor of the dimension or a feature matrix aggregated in the time domain, and the structural magnetic resonance morphological features of the brain's anatomical hardware development state extracted in step G5 above (i.e., a one-dimensional vector that splices together the cortical thickness distribution matrix and the fractal dimension of the brain folds).

[0197] In terms of specific data processing methods, firstly, a matrix vectorization operator is used to flatten the effective elements of the upper triangular region of the dynamic brain function connectivity matrix (excluding diagonal and redundant symmetric elements) into a one-dimensional functional feature vector. Then, a feature normalization operator (such as the Z-score operator) is used to perform cross-modal dimensional normalization on both the functional and structural morphological feature vectors. Finally, a feature concatenation operator is used to concatenate the normalized functional and structural feature vectors along the column dimension, generating a high-dimensional multimodal joint feature vector. This multimodal joint feature vector represents a panoramic composite neuroimaging atlas of the target brain in two core dimensions: dynamic functional network synergy and static anatomical morphological development.

[0198] Multimodal magnetic resonance imaging (MRI) data of a large number of healthy, developing children (Typical Development, TD) of the same age and sex are collected in advance (e.g., a sample size of no less than 500 cases). The multimodal joint feature vector of each healthy sample is extracted using the same method described above. The arithmetic mean vector of the feature vectors of this healthy sample set is calculated, and this mean vector is preset as the center of the spatial distribution of the health baseline feature vector (denoted as ). Simultaneously, the covariance matrix of the feature vectors of the healthy sample set is calculated and inverted, pre-defined as the inverse covariance matrix of the health baseline feature space (denoted as...). ).

[0199] The current multimodal joint feature vector of the target object (denoted as...) ) and the aforementioned preset center vector and the inverse covariance matrix The following matrix quadratic form algebraic operations are performed by calling the Mahalanobis distance calculation formula:

[0200] ;

[0201] Solve for the output Mahalanobis distance It is a dimensionless statistical distance scalar that represents the high-dimensional geometric distance of the target object's multimodal brain features from the typical developmental distribution center of the healthy group after considering the covariance relationship between the features. It serves as the underlying physical scalar for measuring the degree of abnormal deviation of the brain's nervous system.

[0202] G7. Use a preset nonlinear mapping function to convert Mahalanobis distance into a neurological function abnormality evaluation index with values ​​within a preset range.

[0203] Due to the calculated Mahalanobis distance It is a statistical value with a lower limit of 0 but an unbounded upper limit and is non-linear. In order to meet the needs of clinical medical personnel for intuitive image reading and risk classification, it needs to be bounded and non-linearly mapped to clinical semantics.

[0204] In terms of specific data processing methods, a pre-defined parameterized sigmoid nonlinear mapping function is invoked to perform numerical transformation of the Mahalanobis distance. The underlying calculation formula is shown below:

[0205] ;

[0206] In this mapping formula, The output of the computational analysis is a neurological function abnormality assessment index (with values ​​strictly limited to a preset range). Within, such as the preset interval is ); The preset upper limit constant of the risk index (such as the preset upper limit constant) ); The input value is the Mahalanobis distance calculated above.

[0207] This is the preset baseline distance offset threshold. This parameter is pre-calibrated based on the 95th percentile of the Mahalanobis distance in a healthy control group (i.e., the upper limit of the 95% confidence interval) (e.g., calibrated in a healthy sample database). When the Mahalanobis distance of the target object is equal to At that time, the neurological dysfunction evaluation index output by mapping was exactly 50 points, which served as a neutral anchor point to distinguish between low risk and high risk.

[0208] This is the preset mapping scaling sensitivity coefficient. This parameter is pre-calibrated based on the slope corresponding to the maximum Youden's Index of the ROC curve on the clinical validation set (e.g., preset). ), used to control the steepness of the mapping function near the critical point.

[0209] The neurological dysfunction evaluation index output by the mapping directly quantifies the probability level of the target object having autism or being on a high-risk developmental trajectory. This parameter, as an objective, non-invasive and highly interference-resistant quantitative indicator of low-level neuroimaging, is output to an external display terminal for clinicians to conduct comprehensive risk screening in conjunction with the behavioral ADOS scale.

[0210] The above processing smoothly maps the unbounded, high-dimensional, and abstract statistical Mahalanobis distance into a well-defined (0 to 100 points) percentage-based risk index that is semantically intuitive and conforms to clinical diagnostic logic. Meanwhile, the saturation at both ends and sensitivity in the middle of the sigmoid function avoids the distortion caused by extreme outlier data, and significantly amplifies the risk discrimination of borderline suspected cases, greatly enhancing the readability, persuasiveness, and practical value of this assessment method in early clinical screening scenarios.

[0211] like Figure 3 As shown, this application discloses an autism risk assessment system based on brain image feature fusion, including an image sequence acquisition module, a head movement feature extraction module, a gradient state analysis module, a baseline velocity correction module, an artifact weight evaluation module, a sequence compensation reconstruction module, and a feature fusion evaluation module; wherein:

[0212] The image sequence acquisition module is used to acquire the continuous functional magnetic resonance sequence of the target object under the functional task state. The continuous functional magnetic resonance sequence includes multiple sequence frames arranged according to timestamps.

[0213] The head motion feature extraction module is used to extract edge features from the current sequence frame and the reference sequence frame in the continuous functional magnetic resonance imaging sequence, obtain the corresponding geometric centroid of the skull contour, and calculate the spatial displacement vector between the two based on the corresponding geometric centroid of the skull contour.

[0214] The gradient state analysis module is used to calculate the voxel gradient distribution of the current sequence frame and extract the real-time high-frequency gradient dispersion of the current sequence frame.

[0215] The reference velocity correction module is used to perform weighted correction calculations on the preset initial reference velocity using real-time high-frequency gradient dispersion to generate a dynamic reference velocity.

[0216] The artifact weight evaluation module is used to calculate the adaptive exponential decay weight of the current sequence frame based on the spatial displacement vector and the dynamic reference base velocity, and to determine whether the adaptive exponential decay weight is less than the preset artifact threshold.

[0217] The sequence compensation and reconstruction module is used to extract adjacent sequence frames with corresponding adaptive exponential decay weights not less than a preset artifact threshold within a preset time window when the artifact weight evaluation module determines that the sequence frame is correct. The module then uses the interpolation reference frames to perform weighted interpolation reconstruction on the current sequence frame to generate artifact-free sequence frames.

[0218] The feature fusion evaluation module is used to extract the dynamic brain functional connectivity matrix based on the continuous functional magnetic resonance imaging sequence updated by the artifact removal sequence frame, and to fuse the dynamic brain functional connectivity matrix with the structural magnetic resonance morphological features of the target object to generate a neurological function abnormality evaluation index for the target object.

[0219] The technical scope of this invention is not limited to the content described above. Those skilled in the art can make various modifications and variations to the above embodiments without departing from the technical concept of this invention, and all such modifications and variations should fall within the protection scope of this invention.

Claims

1. An autism risk assessment method based on brain imaging feature fusion, characterized in that: Includes the following steps: Acquire a continuous functional magnetic resonance imaging (fMRI) sequence of a target object under a functional task state, wherein the continuous functional magnetic resonance imaging (fMRI) sequence includes multiple sequence frames arranged according to timestamps; Edge features are extracted from the current sequence frame and the reference sequence frame in the continuous functional magnetic resonance imaging (fMRI) sequence to obtain the corresponding geometric centroid of the skull contour, and the spatial displacement vector between the two is calculated based on the corresponding geometric centroid of the skull contour. Calculate the voxel gradient distribution of the current sequence frame and extract the real-time high-frequency gradient dispersion of the current sequence frame; The preset initial reference velocity is weighted and corrected using the real-time high-frequency gradient discreteness to generate a dynamic reference velocity. Based on the spatial displacement vector and the dynamic reference velocity, the adaptive exponential decay weight of the current sequence frame is calculated, and it is determined whether the adaptive exponential decay weight is less than a preset artifact threshold. If the judgment result is yes, then within the preset time window, extract adjacent sequence frames whose corresponding adaptive exponential decay weight is not less than the preset artifact threshold as interpolation reference frames, and use the interpolation reference frames to perform weighted interpolation reconstruction on the current sequence frame to generate artifact-free sequence frames. Based on the dynamic brain functional connectivity matrix extracted from the continuous functional magnetic resonance imaging (fMRI) sequence updated by the artifact removal sequence frames, the dynamic brain functional connectivity matrix is ​​fused with the structural magnetic resonance morphological features of the target object to generate a neurological function abnormality evaluation index for the target object.

2. The autism risk assessment method based on brain imaging feature fusion according to claim 1, characterized in that: The process of obtaining the spatial displacement vector is as follows: Calculate the three-dimensional voxel gradient maps of the current sequence frame and the reference sequence frame, respectively; The three-dimensional voxel gradient map is binarized according to a preset gradient threshold to generate a corresponding skull boundary mask map. The average value of the three-dimensional spatial coordinates of the voxel points at the boundary in the skull boundary mask image is calculated by integrating the three-dimensional spatial coordinates. The average value of the three-dimensional spatial coordinates is used as the geometric centroid of the skull contour corresponding to the current sequence frame and the reference sequence frame, respectively. Spatial difference operations are performed based on the three-dimensional coordinates of the centroid of the skull contour in the current sequence frame and the three-dimensional coordinates of the centroid of the skull contour in the reference sequence frame to construct the spatial displacement vector from the reference sequence frame to the current sequence frame.

3. The autism risk assessment method based on brain imaging feature fusion according to claim 1, characterized in that: The extraction process of the real-time high-frequency gradient discreteness is as follows: Perform three-dimensional multi-scale wavelet decomposition on the voxel matrix of the current sequence frame to extract the corresponding high-frequency detail subband matrix; Calculate the partial derivatives of each voxel in the high-frequency detail subband matrix in the three orthogonal spatial dimensions, and construct the three-dimensional structure tensor matrix of the corresponding voxel based on the partial derivatives; The three-dimensional structural tensor matrix is ​​decomposed into eigenvalues ​​to obtain the principal eigenvalues ​​and minimum eigenvalues ​​for each voxel. Local orientation distortion coefficients are generated based on the ratio between the principal eigenvalues ​​and the minimum eigenvalues ​​for each voxel. Obtain the local orientation distortion coefficients corresponding to all voxels in the high-frequency detail subband matrix, perform a logarithmic summation on all the local orientation distortion coefficients, and use the summation result as the real-time high-frequency gradient discreteness of the current sequence frame.

4. The autism risk assessment method based on brain imaging feature fusion according to claim 3, characterized in that: The process of generating the local directional distortion coefficient is as follows: Extract the relative difference magnitude between the principal feature value and the minimum feature value, and normalize the relative difference magnitude based on the cumulative feature value of the principal feature value and the minimum feature value to generate initial anisotropic features; Extract the tensor traces of the main diagonal elements of the three-dimensional structure tensor matrix and use them as the local structure energy features of the current voxel. A weighted mapping function with monotonically increasing and high-value saturation convergence characteristics is used to perform weighted mapping processing on the local structural energy features to generate the corresponding structural confidence penalty term. The initial anisotropic features are subjected to feature modulation operation based on the structural confidence penalty term, and the modulated output is used as the local orientation distortion coefficient.

5. The autism risk assessment method based on brain imaging feature fusion according to claim 1, characterized in that: The process of generating the dynamic reference velocity is as follows: Obtain the global baseline gradient dispersion of the target object in a resting state; Extract the dynamic deviation characteristics of the real-time high-frequency gradient dispersion relative to the global benchmark gradient dispersion, and generate the dispersion variation index; The dispersion variation index is mapped using a preset dynamic adjustment function to generate an adaptive penalty weight, wherein when the dispersion variation index indicates a positive deviation, the adaptive penalty weight exhibits an accelerated decay distribution. The preset initial reference velocity is dynamically attenuated and modulated using the adaptive penalty weights, and the modulated result is used as the dynamic reference velocity.

6. The autism risk assessment method based on brain imaging feature fusion according to claim 5, characterized in that: The process of generating the dispersion variation index is as follows: Extract the relative residual feature between the real-time high-frequency gradient dispersion and the global benchmark gradient dispersion, and obtain the preset physiological baseline fluctuation tolerance feature; wherein, the physiological baseline fluctuation tolerance feature is used to characterize the normal physiological rhythm noise distribution boundary of the target object in the resting state; The relative residual feature quantity is compared with the physiological baseline fluctuation tolerance feature by performing feature mapping to determine whether the relative residual feature quantity exceeds the upper limit distribution envelope of the physiological baseline fluctuation tolerance feature. If the judgment result is yes, then the effective distortion component of the relative residual feature quantity that exceeds the upper limit distribution envelope is extracted, and the effective distortion component is subjected to logarithmic compression mapping processing. The feature after logarithmic compression mapping processing is output as the dispersion variation index. If the judgment result is negative, the preset distortion-free baseline constant will be output as the dispersion variation index.

7. The autism risk assessment method based on brain imaging feature fusion according to claim 1, characterized in that: The calculation process for the adaptive exponential decay weight is as follows: Extract the spatial magnitude of the spatial displacement vector and obtain the resampling time interval between adjacent sequence frames of the continuous functional magnetic resonance sequence; perform dimensional transformation on the dynamic reference velocity based on the resampling time interval to generate a dynamic reference displacement threshold; Calculate the ratio of the spatial modulus to the dynamic reference displacement threshold to obtain the dimensionless motion parameters; Using the natural constant as the base and the negative of the dimensionless motion parameter as the exponential variable, the adaptive exponential decay weight of the current sequence frame is calculated.

8. The autism risk assessment method based on brain imaging feature fusion according to claim 1, characterized in that: The process of generating the artifact removal sequence frames is as follows: Extract the adaptive exponential decay weights corresponding to the interpolation reference frame, and perform spatial normalization on the extracted adaptive exponential decay weights to generate spatial smoothing kernel weights for the current sequence frame. Using the spatial smoothing kernel weights, convolution interpolation is performed on the voxel matrix of the current sequence frame and the voxel matrix of the interpolation reference frame at corresponding voxel positions to generate the artifact-free sequence frame.

9. The autism risk assessment method based on brain imaging feature fusion according to claim 1, characterized in that: The process of generating the neurological dysfunction evaluation index for the target object is as follows: The updated continuous functional magnetic resonance imaging (fMRI) sequence is mapped onto a preset brain region topology map, and time-dependent blood oxygenation level signal sequences of each target brain region are extracted. The time-dependent blood oxygenation signal sequence is traversed by a sliding time window mechanism, and the Pearson correlation coefficient of each target brain region is calculated in different time windows. The Pearson correlation coefficients in each time window are aggregated to generate the dynamic brain function connectivity matrix. Acquire a high-resolution structural magnetic resonance image of the target object, and perform three-dimensional surface mesh reconstruction on the high-resolution structural magnetic resonance image to generate a cortical inner and outer boundary mesh; Calculate the shortest Euclidean distance between corresponding vertices of the inner and outer boundary grids of the cortex to generate the cortical thickness distribution matrix; Based on the inner and outer boundary grids of the cortex, the fractal dimension of the gyri folds is extracted, and the cortical thickness distribution matrix and the fractal dimension of the gyri folds are combined to form the structural magnetic resonance morphological features. The dynamic brain function connectivity matrix and the structural magnetic resonance morphological features are concatenated into a multimodal joint feature vector, and the Mahalanobis distance between the multimodal joint feature vector and the preset health benchmark feature spatial distribution center is calculated. The Mahalanobis distance is converted into a neurological function abnormality evaluation index with values ​​within a preset range using a preset nonlinear mapping function.

10. An autism risk assessment system based on brain imaging feature fusion, characterized in that: include: The image sequence acquisition module is used to acquire a continuous functional magnetic resonance sequence of a target object under a functional task state, wherein the continuous functional magnetic resonance sequence includes multiple sequence frames arranged according to timestamps. The head motion feature extraction module is used to extract edge features from the current sequence frame and the reference sequence frame in the continuous functional magnetic resonance imaging sequence, obtain the corresponding geometric centroid of the skull contour, and calculate the spatial displacement vector between the two based on the corresponding geometric centroid of the skull contour. The gradient state analysis module is used to calculate the voxel gradient distribution of the current sequence frame and extract the real-time high-frequency gradient dispersion of the current sequence frame. The reference speed correction module is used to perform a weighted correction operation on the preset initial reference speed using the real-time high-frequency gradient dispersion to generate a dynamic reference speed. The artifact weight evaluation module is used to calculate the adaptive exponential decay weight of the current sequence frame based on the spatial displacement vector and the dynamic reference velocity, and to determine whether the adaptive exponential decay weight is less than a preset artifact threshold. The sequence compensation and reconstruction module is used to extract adjacent sequence frames with corresponding adaptive exponential decay weights not less than the preset artifact threshold within a preset time window when the artifact weight evaluation module determines that the current sequence frame is weighted and reconstructed using the interpolation reference frame to generate artifact-free sequence frames. The feature fusion evaluation module is used to extract a dynamic brain functional connectivity matrix based on the continuous functional magnetic resonance imaging sequence updated by the artifact removal sequence frame, and to fuse the dynamic brain functional connectivity matrix with the structural magnetic resonance morphological features of the target object to generate a neurological function abnormality evaluation index for the target object.

Images