Functional magnetic resonance imaging method based on factor analysis and density peak automatic clustering

Abstract

The application discloses a functional magnetic resonance imaging method based on factor analysis and density peak automatic clustering, and belongs to the technical field of image data processing, and comprises the following steps: acquiring M brain fMRI images of a subject during a task experiment and preprocessing the M brain fMRI images into an initial data matrix; reducing the initial data matrix into a dimension-reduced matrix based on a factor analysis algorithm and a Bootstrap bootstrap method, and performing automatic clustering on time series in the dimension-reduced matrix to obtain clustering clusters based on an improved DPC algorithm; and mapping the clustering clusters back to a brain structure diagram to obtain a functional activity imaging diagram. The application firstly proposes to combine the factor analysis and the Bootstrap bootstrap method to filter and reduce high-dimensional fMRI data, select the most important feature data, greatly reduce the calculation amount of a subsequent clustering algorithm, and solve the shortcoming that fMRI image data cannot be clustered due to a large amount of data; and then the improved DPC algorithm is used for automatic clustering, the difficulty of manually selecting a clustering center in the existing DPC clustering is overcome, and the influence of human factors on the number of clusters is avoided.

Description

Technical Field

[0001] This invention relates to the field of image data processing technology, and in particular to a functional magnetic resonance imaging method based on factor analysis and automatic density peak clustering. Background Technology

[0002] Functional magnetic resonance imaging (fMRI) is a novel, non-invasive technique for detecting and imaging brain functional activity, providing clear anatomical and functional images. Using fMRI, the activity of multiple brain regions can be observed simultaneously, exploring the relationships between different functional areas, which is highly beneficial for studying higher brain functions. fMRI has made significant progress in basic research on normal human brain functional areas (vision, hearing, smell, motor, sensory, and language, etc.). fMRI plays a crucial role in studying the brain activity corresponding to the human motor, auditory, language, perceptual, emotional, and memory systems. Furthermore, fMRI has extensive clinical applications in neurosurgery, neurology, psychiatry, and mental disorders.

[0003] During fMRI data acquisition, hundreds of 3D brain function activation images can be obtained, with hundreds of thousands of voxels in each image. However, active voxels are only a small fraction of these images, while a large number of inactive voxels remain. If all voxels are used as input data for cluster analysis, it will consume a lot of time and computing resources and can easily lead to memory overflow.

[0004] Furthermore, when performing cluster analysis on fMRI data, clustering methods include K-means clustering, fuzzy clustering, and hierarchical clustering. However, these methods all have certain drawbacks. For example, K-means clustering requires pre-specifying the number of clusters K, and the choice of initial cluster centers has a significant impact on the clustering results; fuzzy clustering is particularly sensitive to initialization parameters, easily gets trapped in local minima, and requires prior knowledge of the number of clusters, which undermines the algorithm's adaptability and unsupervised nature; hierarchical clustering cannot undo splitting operations, cannot swap objects between clusters, and the choice of split points directly affects the clustering results, making it unsuitable for large datasets and having high time complexity.

[0005] In 2014, Alex Rodriguez et al. proposed a novel density-based clustering algorithm in Science, namely Density Peak Clustering (DPC). This algorithm proposes the following characteristics of cluster centers: cluster centers have high local density and are relatively far from other points with even higher density. Based on these characteristics, the DPC algorithm calculates the local density and relative distance of points, draws a decision map (an image composed of local density on the X-axis and relative distance on the Y-axis), manually selects cluster centers based on the decision map, and then classifies non-cluster centers to achieve the final clustering goal. However, manually selecting cluster centers based on the decision map is not ideal for complex brain function datasets, where the distribution of points on the decision map is also complex. Directly using local density and distance values ​​is not ideal, and it is difficult to accurately select the correct cluster centers manually.

[0006] Definitions: p-value: In statistics, researchers typically use hypothesis testing to assess differences or effects between two or more groups. In hypothesis testing, researchers calculate a statistic based on sample data and then calculate the probability of obtaining that statistic or a more extreme outcome under the null hypothesis (i.e., no difference or effect). This probability is the p-value, also known as the p-parameter value. The p-value is commonly used in significance testing. The criterion for significance testing is usually based on comparing the significance level (α, such as 0.05 or 0.01) with the calculated p-value. If p ≤ α, the null hypothesis is rejected, and the result is considered statistically significant. Summary of the Invention

[0007] The purpose of this invention is to provide a functional magnetic resonance imaging method based on factor analysis and density peak automatic clustering, which solves the problems mentioned above, such as the large amount of fMRI image data leading to large computational load of clustering algorithms, and the difficulties in parameter selection, time consumption, dependence on individual experience, and difficulty in manually selecting cluster centers of traditional clustering algorithms.

[0008] To achieve the above objectives, the technical solution adopted by this invention is as follows: a functional magnetic resonance imaging method based on factor analysis and automatic density peak clustering, comprising the following steps: S1. Acquire and preprocess M brain fMRI images of the subject during the task experiment, and label them as I1~I2. M The m-th image is I. m , will I m Divide into U voxels and read I m The gray values ​​corresponding to each voxel are used to obtain an initial U×M data matrix V, where the u-th row vector v in V is... u Given a time series of voxel u, the element v in the u-th row and m-th column of V. umFor voxel u in I m In the grayscale values, 1≤u≤U, 1≤m≤M; S2, based on factor analysis and Bootstrap, reduces the dimensionality of the initial matrix V to obtain an N×M dimensionality-reduced matrix S, where N is the total number of voxels after dimensionality reduction, and the element in the i-th row of S is labeled as the time series x. i ,but , 1≤i≤N; S3, based on the improved DPC algorithm, performs automatic functional area clustering on S, including S31~S35; S31, calculate each x i Local density ρ i and distance δ i ; S32, ρ i and δ i Normalize to the interval [0,1] to obtain the corresponding normalized value ρ. i_norm and δ i_norm And calculate ρ i and δ i absolute deviation ρ i_abs and δ i_abs ; S33, calculate x according to the following formula i The corresponding decision value γ i , forming a sequence ; , S34, to Press γ i Rearrange the values ​​in descending order to obtain the sorted vector. , where γ g for For the g-th element, substitute the sorted vector into the Matlab fitting function to fit a cubic polynomial, and label the value of the cubic polynomial in g as f(g). Compare each g in turn. If γ g >f(g), will γ g γ before rearrangement i time series x i As cluster centers, a total of n were identified. c Cluster centers; S35, based on the distance clustering algorithm, the time series of the dimensionality-reduced matrix S is clustered according to the cluster centers to obtain n c One cluster; S4 maps the clusters back to the brain structure map to obtain a functional activity imaging map.

[0009] Preferably, in S1, the preprocessing involves: performing head motion correction on the fMRI image, normalizing the image to an EPI template, spatial smoothing, and a 1 / 128Hz high-pass filter, then reading it into Matlab, resulting in the initial matrix V: .

[0010] Preferably, S2 includes S21 to S23: S21, the initial matrix V is reduced in dimensionality using factor analysis algorithm to obtain a U×k dimensionality-reduced matrix V1. U×k k is the number of factors selected based on contribution rate in the factor analysis algorithm; S22, for V1 U×k For each column vector, nonparametric statistics are performed using the Bootstrap algorithm to obtain the p-parameter value corresponding to each element in the vector. V1 U×k Set p parameter values ​​less than p < 0.05 to 1, and set all other p parameter values ​​to 0, to obtain a 0-1 matrix V2. U×k The u-th row corresponds to voxel u; S23, process V2 sequentially U×k In each row, if V2 U×k If at least one element in row u is 1, then voxel u is selected. A total of N voxels are selected. The row vectors corresponding to these N voxels are extracted from V to form an N×M dimension reduction matrix S.

[0011] As a preferred option, the factor selection based on contribution rate specifically selects factors with a contribution rate greater than 4%.

[0012] As a preferred option, in S31, x is calculated. i Local density ρ i Including Sa1~Sa3: Sa1, generates the index set I of S. s I s ={1,2,⋯,N}, where j is I s The number of elements is not equal to i. Sa2, calculate x i The Euclidean distances with other time series in S are T in total, where x i With the j-th time series x in S j The Euclidean distance is d ij Arrange the T Euclidean distances in ascending order as d1~d2. T Calculate the cutoff distance d c d c =d sda(T×2%) sda(∙) means rounding to the nearest integer; Sa3, the local density ρ is calculated according to the following formula. i ; , In the formula, e is the natural constant; In S32, calculate x i distance δ i Specifically: If x i If the local density is maximized, then d ij The maximum value is taken as δ i ; If x i If no local density is found, then search for all local densities greater than x. i Time series and x i The Euclidean distance is taken as the minimum value of the Euclidean distance. i .

[0013] As a preferred option, in S32, ρ i and δ i absolute deviation ρ i_abs and δ i_abs Calculate according to the following formulas respectively: , , In the formula, mean(∙) is used to calculate the average value.

[0014] Preferably, in S35, the distance in the distance clustering algorithm includes Euclidean distance.

[0015] Compared with the prior art, the advantages of the present invention are as follows: (1) A method for dimensionality reduction of fMRI data based on factor analysis algorithm and bootstrap method is proposed. First, the initial data matrix V of U×M is reduced to a dimensionality reduction matrix V1 of U×k by factor analysis algorithm. U×k k is the number of factors selected based on contribution rate in the factor analysis algorithm, and then the Bootstrap method is used to analyze V1. U×k For each column, generate a 0-1 matrix V2 U×k Based on V2 U×k The time series of voxels are selected from the initial data matrix V, and finally an N×M dimension-reduced matrix S is obtained. This method can select the most important activation voxels, greatly reducing the computational load of clustering algorithms and solving the problem that large amounts of fMRI image data cannot be clustered.

[0016] (2) Based on the above dimensionality reduction results, a density peak clustering algorithm was introduced for fMRI images. This algorithm overcomes the problems of difficult parameter selection, large time overhead and individual experience dependence of traditional clustering algorithms, and does not require prior information.

[0017] (3) An automatic clustering method is proposed based on the density peak clustering algorithm (DPC). Due to the complex distribution of cluster center decision points, DPC directly uses the local density ρ i and distance δ i This is not ideal, therefore the present invention uses ρ i and δ i The normalized values ​​ρ of the two i_norm and δ i_norm The absolute deviation ρ between the two i_abs and δ i_abs Calculate x i The corresponding decision value γ i The sequence is constructed, fitted into a cubic polynomial using a Matlab fitting function, and then the cluster centers are automatically detected according to the detection method of this invention. This method overcomes the difficulty of manually selecting cluster centers in DPC clustering and avoids the influence of human factors on the number of clusters.

[0018] In summary, this invention introduces factor analysis and automatic density peak clustering into functional magnetic resonance imaging, which can not only identify brain activity areas and reduce computational load, but also automatically complete clustering, avoiding the influence of human factors on clustering. Attached Figure Description

[0019] Figure 1 This is a flowchart of the present invention; Figure 2 This is a functional activity imaging diagram obtained in Example 3. Detailed Implementation

[0020] The present invention will be further described below with reference to the embodiments and accompanying drawings.

[0021] Example 1: See Figure 1 A functional magnetic resonance imaging method based on factor analysis and automatic density peak clustering includes the following steps: S1. Acquire and preprocess M brain fMRI images of the subject during the task experiment, and label them as I1~I2. M The m-th image is I. m , will I m Divide into U voxels and read I m The gray values ​​corresponding to each voxel are used to obtain an initial U×M data matrix V, where the u-th row vector v in V is... u Given a time series of voxel u, the element v in the u-th row and m-th column of V. um For voxel u in I m In the grayscale values, 1≤u≤U, 1≤m≤M; S2, based on factor analysis and Bootstrap, reduces the dimensionality of the initial matrix V to obtain an N×M dimensionality-reduced matrix S, where N is the total number of voxels after dimensionality reduction, and the element in the i-th row of S is labeled as the time series x. i ,but , 1≤i≤N; S3, based on the improved DPC algorithm, performs automatic functional area clustering on S, including S31~S35; S31, calculate each x i Local density ρ i and distance δ i ; S32, ρ i and δ i Normalize to the interval [0,1] to obtain the corresponding normalized value ρ. i_norm and δ i_norm And calculate ρ i and δ i absolute deviation ρ i_abs and δ i_abs ; S33, calculate x according to the following formula i The corresponding decision value γ i , forming a sequence ; , S34, to Press γ i Rearrange the values ​​in descending order to obtain the sorted vector. , where γ g for For the g-th element, substitute the sorted vector into the Matlab fitting function to fit a cubic polynomial, and label the value of the cubic polynomial in g as f(g). Compare each g in turn. If γ g >f(g), will γ g γ before rearrangement i time series x i As cluster centers, a total of n were identified. c Cluster centers; S35, based on the distance clustering algorithm, the time series of the dimensionality-reduced matrix S is clustered according to the cluster centers to obtain n c One cluster; S4 maps the clusters back to the brain structure map to obtain a functional activity imaging map.

[0022] Example 2, see Figure 1 This embodiment is more specific than embodiment 1: In S1, the preprocessing involves: performing head motion correction on the fMRI image, normalizing the image to an EPI template, spatial smoothing, and a 1 / 128Hz high-pass filter, then reading it into Matlab to obtain the initial matrix V: .

[0023] S2 includes S21~S23: S21, the initial matrix V is reduced in dimensionality using factor analysis algorithm to obtain a U×k dimensionality-reduced matrix V1. U×k k is the number of factors selected based on contribution rate in the factor analysis algorithm; and specifically, selecting factors based on contribution rate means selecting factors with a contribution rate greater than 4%. S22, for V1 U×k For each column vector, nonparametric statistics are performed using the Bootstrap algorithm to obtain the p-parameter value corresponding to each element in the vector. V1 U×k Set p parameter values ​​less than p < 0.05 to 1, and set all other p parameter values ​​to 0, to obtain a 0-1 matrix V2. U×k The u-th row corresponds to voxel u; S23, process V2 sequentially U×k In each row, if V2 U×k If at least one element in row u is 1, then voxel u is selected. A total of N voxels are selected. The row vectors corresponding to these N voxels are extracted from V to form an N×M dimension reduction matrix S.

[0024] In S31, calculate x i Local density ρ i Including Sa1~Sa3: Sa1, generates the index set I of S. s I s ={1,2,⋯,N}, where j is I s The number of elements is not equal to i. Sa2, calculate x i The Euclidean distances with other time series in S are T in total, where x i With the j-th time series x in S j The Euclidean distance is d ij Arrange the T Euclidean distances in ascending order as d1~d2. T Calculate the cutoff distance d c d c =d sda(T×2%) sda(∙) means rounding to the nearest integer; Sa3, the local density ρ is calculated according to the following formula. i ; , In the formula, e is the natural constant; In S32, calculate x i distance δ i Specifically: If x i If the local density is maximized, then d ij The maximum value is taken as δ i ; If x i If no local density is found, then search for all local densities greater than x. i Time series and x i The Euclidean distance is taken as the minimum value of the Euclidean distance. i .

[0025] In S32, ρ i and δ i absolute deviation ρ i_abs and δ i_abs Calculate according to the following formulas respectively: , , In the formula, mean(∙) is used to calculate the average value.

[0026] In S35, the distance in the distance clustering algorithm includes Euclidean distance.

[0027] Example 3: See Figure 1 and Figure 2 Based on Examples 1 and 2, in order to illustrate the process and effects of the present invention, this example uses block experiment fMRI images of finger movement tasks for analysis. The specific steps are as follows:

[0028] The fMRI experimental mode for the task block experiment: While the subject was being scanned, a computer generated visual stimuli according to the E-prime program (www.pstnet.com / e-prime) and projected them onto a screen above the subject's head. When the subject saw a small cross on the screen, it indicated a rest period. When the letter "R" appeared on the screen, it indicated grasping the right finger. When the letter "L" appeared, it indicated grasping the left finger. The rest period and the task period alternated for 20 seconds each. A total of 160 whole-brain images, or 160 fMRI images, were generated, representing the time allotted for each hand movement.

[0029] S1, acquire the above 160 fMRI images and preprocess them. The preprocessing is as follows: (1) Perform head motion correction on the fMRI images; (2) Normalize the fMRI images to the SPM2 EPI standard template; (3) Spatial smooth the data using a 12mm half-width at half-maximum (WHM); (4) Filter out low-frequency physiological noise components caused by breathing, etc., using a high-pass filter with a frequency of 1 / 128Hz; (5) Use Matlab to read the fMRI images processed by (1), (2), and (3) above into the initial matrix V. ; The initial data matrix V has U rows and M columns, and the vector v in the u-th row of V is... u Given a time series of voxel u, the element v in the u-th row and m-th column of V. um For voxel u in I m The grayscale values ​​in the data are 1≤u≤U, 1≤m≤M. In this embodiment, the total number of voxels is U=41×48×23=45264. S2, based on factor analysis and Bootstrap, reduces the dimensionality of the initial matrix V. Specifically, it first reduces the dimensionality using the factor analysis algorithm from S21, selecting five factors with a contribution rate greater than 4%, resulting in a U×k dimensionality-reduced matrix V1. U×k Where U is 45264 and k is 5; then generate a 0-1 matrix according to S22, and select voxels from V according to S23 to obtain an N×M dimension reduction matrix S.

[0030] S3, perform automatic functional area clustering on S, the same as step S3 in Example 1; S4. Map the clusters back to the brain structure map. Specifically, using the BrainNet Viewer toolkit in Matlab, the functionally related brain activity areas of the left and right hands identified by the clusters are projected back onto the brain structure glass image template to obtain a functional activity imaging map, such as... Figure 2 As shown.

[0031] Figure 2 It includes five perspectives of the brain, with "L" on the left representing the left hemisphere and "R" on the right representing the right hemisphere. From... Figure 2 It can be seen that the activity areas of left and right hand movements identified by the present invention mainly cover the left and right motor cortexes, and the red and orange colors indicate that the identified left and right motor cortexes belong to different clusters, which also indicates that the brain activity signals of the left and right hands can be correctly distinguished by the invention method.

[0032] The above description is merely a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions, and improvements made within the spirit and principles of the present invention should be included within the protection scope of the present invention.

Claims

1. A functional magnetic resonance imaging method based on factor analysis and automatic density peak clustering, characterized in that, Includes the following steps: S1. Acquire and preprocess M brain fMRI images of the subject during the task experiment, and label them as I1~I2. M The m-th image is I. m , will I m Divide into U voxels and read I m The gray values ​​corresponding to each voxel are used to obtain an initial U×M data matrix V, where the u-th row vector v in V is... u Given a time series of voxel u, the element v in the u-th row and m-th column of V. um For voxel u in I m In the grayscale values, 1≤u≤U, 1≤m≤M; S2, based on factor analysis and Bootstrap, reduces the dimensionality of the initial matrix V to obtain an N×M dimensionality-reduced matrix S, where N is the total number of voxels after dimensionality reduction, and the element in the i-th row of S is labeled as the time series x. i ,but , 1≤i≤N; S3, based on the improved DPC algorithm, performs automatic functional area clustering on S, including S31~S35; S31, calculate each x i Local density ρ i and distance δ i ; S32, ρ i and δ i Normalize to the interval [0,1] to obtain the corresponding normalized value ρ. i_norm and δ i_norm And calculate ρ i and δ i absolute deviation ρ i_abs and δ i_abs ; S33, calculate x according to the following formula i The corresponding decision value γ i , forming a sequence ; ， S34, to Press γ i Rearrange the values ​​in descending order to obtain the sorted vector. , where γ g for For the g-th element, substitute the sorted vector into the Matlab fitting function to fit a cubic polynomial, and label the value of the cubic polynomial in g as f(g). Compare each g in turn. If γ g >f(g), will γ g γ before rearrangement i time series x i As cluster centers, a total of n were identified c Cluster centers; S35, based on the distance clustering algorithm, the time series of the dimensionality-reduced matrix S is clustered according to the cluster centers to obtain n c One cluster; S4 maps the clusters back to the brain structure map to obtain a functional activity imaging map.

2. The functional magnetic resonance imaging method based on factor analysis and automatic density peak clustering according to claim 1, characterized in that, In S1, the preprocessing involves: performing head motion correction on the fMRI image, normalizing the image to an EPI template, spatial smoothing, and a 1 / 128Hz high-pass filter, then reading it into Matlab to obtain the initial matrix V: 。 3. The functional magnetic resonance imaging method based on factor analysis and automatic density peak clustering according to claim 1, characterized in that, S2 includes S21~S23: S21, the initial matrix V is reduced in dimensionality using factor analysis algorithm to obtain a U×k dimensionality-reduced matrix V1. U×k k is the number of factors selected based on contribution rate in the factor analysis algorithm; S22, for V1 U×k For each column vector, nonparametric statistics are performed using the Bootstrap algorithm to obtain the p-parameter value corresponding to each element in the vector. V1 U×k Set p parameter values ​​less than p < 0.05 to 1, and set all other p parameter values ​​to 0, to obtain a 0-1 matrix V2. U×k The u-th row corresponds to voxel u; S23, process V2 sequentially U×k In each row, if V2 U×k If at least one element in row u is 1, then voxel u is selected. A total of N voxels are selected. The row vectors corresponding to these N voxels are extracted from V to form an N×M dimension reduction matrix S.

4. The functional magnetic resonance imaging method based on factor analysis and automatic density peak clustering according to claim 3, characterized in that, The selection of factors based on contribution rate specifically involves selecting factors with a contribution rate greater than 4%.

5. The functional magnetic resonance imaging method based on factor analysis and automatic density peak clustering according to claim 1, characterized in that, In S31, calculate x i Local density ρ i Including Sa1~Sa3: Sa1, generates the index set I of S. s I s ={1,2,⋯,N}, where j is I s The number of elements is not equal to i. Sa2, calculate x i The Euclidean distances with other time series in S are T in total, where x i With the j-th time series x in S j The Euclidean distance is d ij Arrange the T Euclidean distances in ascending order as d1~d2. T Calculate the cutoff distance d c d c =d sda(T×2%) sda(∙) means rounding to the nearest integer; Sa3, the local density ρ is calculated according to the following formula. i ; ， In the formula, e is the natural constant; In S32, calculate x i distance δ i Specifically: If x i If the local density is maximized, then d ij The maximum value is taken as δ i ; If x i If no local density is found, then search for all local densities greater than x. i Time series and x i The Euclidean distance is taken as the minimum value of the Euclidean distance. i .

6. The functional magnetic resonance imaging method based on factor analysis and automatic density peak clustering according to claim 1, characterized in that, In S32, ρ i and δ i absolute deviation ρ i_abs and δ i_abs Calculate according to the following formulas respectively: ， ， In the formula, mean(∙) is used to calculate the average value.

7. The functional magnetic resonance imaging method based on factor analysis and automatic density peak clustering according to claim 1, characterized in that, In S35, the distance in the distance clustering algorithm includes Euclidean distance.

Images