CN118351406A - A globally optimized supervised multimodal image fusion method - Google Patents

Abstract

The present invention discloses a globally optimized supervised multimodal image fusion method, which preprocesses brain image data of each modality and performs feature extraction; preprocesses the obtained features of each modality to obtain the feature matrix of each modality after preprocessing; calculates the independence of the independent components of each modality; and solves the optimal solution of the mixing matrix when maximizing the independence of the independent components of each modality, the correlation between the mixing matrices of each modality, and the correlation between the mixing matrices of each modality and reference information. The global optimal solution is obtained by maximizing the independence of the independent components and maximizing the sum of the squares of the correlations between the independent components of each modality and between the independent components and the reference information. By introducing working memory scores as reference information, the synergistic covariant components that are related to clinical indicators and relatively independent in space are specifically identified. Thus, the interactive relationship between cognitive ability and multimodal neuroimaging in mental illness is deeply explored.

Description

A globally optimized supervised multimodal image fusion method

Technical Field

The present invention relates to medical image analysis, and in particular to a globally optimized supervised multimodal image fusion method.

Background technique

In recent years, the rapid development of neuroimaging technology has enabled researchers to use a variety of imaging methods to collect neuroimages of the same subject in different modalities, thereby providing information about the brain's anatomical structure or function in a multi-dimensional manner. For example, common structural imaging techniques include structural magnetic resonance imaging (sMRI), which can be used to evaluate changes in the local concentration or volume of gray matter and white matter on each voxel, thereby reflecting the corresponding anatomical changes. In the field of functional imaging, functional magnetic resonance imaging (fMRI) is based on blood oxygen level-dependent signals to reflect the activity of brain neurons in task or resting states. All of the above imaging methods have their own advantages and limitations. Using a single imaging mode alone cannot simultaneously obtain complete information on brain function and structure. Although the multimodal imaging data of the scanned object coexist, they are often analyzed separately in actual operations, or the analysis results of the two magnetic resonance imaging data are simply compared or correlated, and the complementary information between the modalities cannot be fully utilized. Multimodal data fusion analysis methods can integrate complementary information between modal data, provide a deeper understanding of the pathological mechanisms associated with neuropsychiatric diseases, and thus improve diagnostic and prediction accuracy.

Supervised fusion methods are more goal-oriented because they use reference information to guide fusion analysis, which can accurately locate specific components of interest from large and complex data sets, effectively improving the targeting of feature extraction. pICA-R (pICA with references) is a method that uses genes as reference information to guide multimodal fusion. However, the reference information in this method is limited to the distribution of spatial components, such as specific brain regions or specific sites, and cannot cover the indicators of individual subjects, such as various behavioral scores of subject samples. In addition, the fusion effect is heavily dependent on the accuracy of genetic information. MCCAR+jICA (multi-set CCA with reference+joint ICA), while maintaining the original basic performance of joint source separation, can further detect covariant multimodal features that are significantly correlated with reference information. However, this method optimizes the objective function in stages, which may cause the correlation extracted in the first stage to be lost in the second stage analysis.

Summary of the invention

Purpose of the invention: In view of the above shortcomings, the present invention provides a globally optimized supervised multimodal image fusion method with improved accuracy.

Technical solution: To solve the above problems, the present invention adopts a globally optimized supervised multimodal image fusion method, which includes the following steps:

(1) Preprocessing brain image data of each modality and performing feature extraction to obtain features of each modality;

(2) Preprocessing the obtained modal features to obtain the feature matrix of each modality after preprocessing;

(3) The independence of the independent components of each mode calculated based on the characteristic matrix of each mode, the mixing matrix of each mode and the weight matrix of each mode;

(4) Finding the optimal solution of the mixing matrix that maximizes the independence of the independent components of each mode, the correlation between the mixing matrices of each mode, and the correlation between the mixing matrices of each mode and the reference information;

(5) Based on the mixing matrix obtained in step (4), calculate the covariant component and mixing matrix of the multimodal information related to the reference information.

Furthermore, the independence of the independent components of each modality in step (3) is described by the information entropy function H(·), and the calculation formula for maximizing the independence of the independent components of each modality is:

U _k =W _k ·X _k +W _k0 ;

Among them, E(.) is the expected value; is the probability density function; Y _k and U _k are the intermediate variables of the calculation; W _k is the unmixing matrix of the mixing matrix _Ak ; _Ak is the mixing matrix of mode k; W _k0 is the bias weight matrix of mode k; X _k is the feature matrix of mode k after preprocessing.

Furthermore, the solution formula for maximizing the independence of the independent components of each mode, the correlation between the mixing matrices of each mode, and the correlation between the mixing matrices of each mode and the reference information in step (4) is:

Among them, m _k (k＝1,2) is the number of independent components of mode k; A _1i is the i-th column of the mixing matrix of mode 1; A _2j is the j-th column of the mixing matrix of mode 2; ref _h represents the h-th column of the reference information matrix; Corr(A _1i , A _2j ) represents the correlation between A _1i and A _2j ; Corr(A _1i , ref _k ) represents the correlation between A _1i and ref _h ; Corr(A _2j , ref _h ) represents the correlation between A ₂ j and ref _h ; α ₁ , α ₂ , and α ₃ represent regularization parameters.

Furthermore, the gradient descent algorithm is used to iterate the solution formula until convergence, and the optimal solution of the mixing matrix is obtained:

The iteration rules of the unmixing matrix of each mode satisfy:

Where _λk is the learning rate for independent component analysis of mode k; I is the identity matrix; and T is the transpose of the matrix.

Furthermore, the step of maximizing the correlation between the mixing matrices of each mode and the correlation between the mixing matrices of each mode and the reference information in step (4) specifically includes:

(41) Calculate the correlation between the mixing matrix of each mode and the reference information, and select the triplet A _1i , A _2j and ref _k with the largest correlation as the constrained items;

(42) Taking the partial derivatives of the sum of squared correlations with respect to A _1i and A _2j respectively and setting them equal to 0, we obtain the partial differential equation containing A _1i and A _2j ;

(43) Setting an initial point, iteratively updating the A _1i and A _2j , and ensuring that the value of the correlation square sum increases until the correlation square sum converges, and taking the A _1i and A _2j at this time as the optimal solution of this iteration process;

(44) Update A ₁ and A ₂ according to the optimal solution in step (43), and repeat steps (41), (42) and (43) to obtain the optimal solutions A _1i and A _2j .

Furthermore, the iterative update rule for iteratively updating A _1i and A _2j in step (43) satisfies:

Among them, Std(·) is the standard deviation; Cov(·) is the covariance; var(·) is the variance; They respectively represent the average values of A _1i , A _2j , and ref _k ; η ₁₁ , η ₁₂ , η ₂₁ , and η ₂₂ are all descent steps obtained by the gradient descent algorithm; λ _c1 and λ _c2 are learning rates.

Furthermore, each modality image in step (1) includes functional magnetic resonance image fMRI and structural magnetic resonance image sMRI. The modality features extracted from the functional magnetic resonance image fMRI include fractional low-frequency fluctuation amplitude; the modality features extracted from the structural magnetic resonance image sMRI include gray matter volume. In step (2), the preprocessing of each modality feature obtained includes: matrixing, centering, whitening and dimension reduction using principal component analysis method.

Beneficial effects: Compared with the prior art, the significant advantage of the present invention is that it maximizes the independence of independent components and simultaneously maximizes the sum of squares of correlations between independent components of each modality and between independent components and reference information, thereby obtaining a global optimal solution. By introducing working memory scores as reference information, synergistic covariant components that are related to clinical indicators and relatively independent in space are specifically identified. This allows for in-depth exploration of the interactive relationship between cognitive ability and multimodal neuroimaging in mental illness, providing strong support for understanding the pathological mechanisms of mental illness and formulating personalized treatment plans.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG1 is a schematic diagram of a process of a multimodal brain image fusion method of the present invention;

FIG2(a) is a schematic diagram showing the comparison between the absolute values of the correlation between the target components of each modality and the reference information estimated by MCCA, jICA, MCCA+jICA, separate ICA, PICA, MCCAR, MCCAR+jICA and the fusion method of the present invention and the actual absolute values of the correlation between the target components of each modality and the reference information;

FIG2( b ) is a schematic diagram showing the comparison of the accuracy of the mixed matrix and independent components of MCCA, jICA, MCCA+jICA, separate ICA, PICA, MCCAR, MCCAR+jICA and the fusion method of the present invention relative to the average of the absolute values of the true correlation between the target components of each modality and the reference information;

FIG3( a ) is a schematic diagram showing the accuracy of MCCA, jICA, MCCA+jICA, separate ICA, PICA, MCCAR, MCCAR+jICA and the supervised PICA method of the present invention relative to the independent components averaged over 14 noise levels;

FIG3( b ) is a schematic diagram showing the accuracy of the MCCA, jICA, MCCA+jICA, separate ICA, PICA, MCCAR, MCCAR+jICA, and supervised PICA methods relative to the mixing matrix averaged over 14 noise levels;

FIG4( a ) is a spatial view of the multimodal joint covariant components associated with working memory scores detected by the fusion method of the present invention;

Figure 4(b) shows the intergroup differences of each modality between the schizophrenia group and the healthy control group;

Figure 4(c) is a fitted graph of the correlation between each modality and working memory score;

Figure 4(d) is the correlation fitting diagram between each mode.

Detailed ways

As shown in FIG1 , a globally optimized supervised multimodal image fusion method in this embodiment includes the following steps:

S1: Preprocess the brain imaging data of each modality and extract its features.

Each modality can be a resting state functional brain imaging fMRI, structural magnetic resonance imaging sMRI, diffusion tensor imaging dMRI and other brain imaging technologies. In this embodiment, two modalities, fMRI and sMRI, are used, but in actual applications, multiple modalities can be selected and implemented as needed.

Although this embodiment only takes the multimodal data fusion of schizophrenia as an example for illustration, the protection scope of the present invention is not limited thereto. The data of each modality can cover the multimodal data of all mental disorders, including but not limited to mental illnesses such as schizophrenia and depression, and multimodal data of non-mental diseases such as bipolar disorder.

In this embodiment, the fractional amplitude fALFF of low-frequency fluctuations is extracted from the resting-state functional brain image fMRI as its feature. The preprocessing and feature extraction of fMRI are specifically as follows: in the MATLAB 2019 environment, standard preprocessing based on statistical parameter mapping is performed on the data as follows: (1) head motion correction, the image at each time point is adjusted to ensure the accurate correspondence of brain structure during the analysis process; (2) inter-layer time correction, the acquisition time difference between different layers (slices) in the fMRI image is adjusted; (3) standardization to the Montreal MNI standard space and resampling; preferably, resampling to 3×3×3mm; (4) regressing out 6 head motion parameters and white matter (WM), gray matter (GM), cerebrospinal fluid (CSF); (5) using 8mm full width at half maximum (full width half maximum) (6) The fractional amplitude of low-frequency fluctuations (fALFF) was calculated by dividing the sum of the amplitude values in the low-frequency power range of 0.01 to 0.08 Hz by the sum of the amplitudes of the entire detectable power spectrum.

In this embodiment, the gray matter volume GMV is extracted from the structural magnetic resonance image sMRI as its feature. The preprocessing and feature extraction of sMRI are specifically as follows: (1) using the unified segmentation method in SPM12, the sMRI data is standardized to the Montreal MNI standard space and resampled; preferably, the data is resampled to 3×3×3 mm; (2) the data is segmented into gray matter, white matter and cerebrospinal fluid, and the result is output as gray matter volume; (3) the gray matter is smoothed using an 8 mm half-maximum full-width Gaussian filter kernel; and (4) the abnormal values of the subjects are detected to ensure the correct segmentation.

After the above preprocessing of the fMRI and sMRI data of each subject, the features of each modality are converted into a 3D matrix. For fMRI, the features that can be used include the fractional low-frequency fluctuation amplitude feature. For sMRI, the modality features that can be used include the gray matter volume feature.

S2: Preprocessing of each modal feature, including matrixing, centering, whitening and dimension reduction using principal component analysis.

The matrix processing process can be as follows: the 3D feature matrix of each modality of each subject obtained after the processing of step S1 is stretched into a row vector (1×L), where L is the number of voxels, and the 3D matrices of all subjects are transformed to obtain an N×L matrix X _k (k＝1,2), where N represents the number of subjects and k represents the number of modalities. The above transformation is performed on both fMRI and sMRI modalities to obtain the feature matrices of the two modalities.

Next, the feature matrix is centered and whitened to eliminate the translation effect of the data so that it has the same variance and is uncorrelated with each other, ensuring that the feature values of the two modes are in the same range.

Finally, the principal component analysis method is used to reduce the dimension of the above feature matrix to remove redundant information and irrelevant features, thereby reducing computational complexity and reducing noise interference. This helps to retain the main direction of data change and improve the generalization performance of the model.

S3: Based on the reduced-dimensional features of each modal obtained in step S2, the reference information matrix is added, and the independence of each modal component, the correlation between the modalities, and the correlation between the components and the reference information are maximized. The gradient descent algorithm is used for cyclic iteration until convergence.

Specifically, the maximization is achieved according to the following formula:

In formula (1), H(·) is the information entropy function that ensures the independence of independent components of each mode; m _k (k＝1,2) represents the number of independent components of mode k; A _1i represents the i-th column of the mixing matrix of mode 1; A _2j represents the j-th column of the mixing matrix of mode 2; ref _k represents the k-th column of the reference information matrix; Corr(A _1i ,A _2j ) represents the correlation between A _1i and A _2j ; Corr(A _1i ,ref _k ) represents the correlation between A _1i and ref _k ; Corr(A _2j ,ref _k ) represents the correlation between A _2j and ref _k ; α ₁ , α ₂ , and α ₃ represent the control and The regularization parameter of the weights between them and ensures the convergence of the algorithm optimization process.

Specifically, the independent component analysis of each modal feature is performed in parallel according to the following formula:

in, and They represent the probability density functions of _Y1 and _Y2 respectively; E(.) is the expected value; H(.) is the entropy function; _Wk (k＝1,2) is the inverse of the mixing matrix _Ak , that is, the unmixing matrix of _Ak ; _Wk0 (k＝1,2) is the bias weight matrix of mode k; _Xk (k＝1,2) represents the feature matrix of mode k after preprocessing.

For solving ICA, the widely used Infomax algorithm is adopted here.

Solving the maximization problem of formula (1) may specifically include: calculating the correlation of each column of A ₁ , A ₂ and ref, and selecting the triples A _1i , A _2j and ref _k with the largest correlation as constrained items. Taking partial derivatives of the sum of squares of correlation with respect to A _1i and A _2j , respectively, and setting them equal to 0, a partial differential equation containing A _1i and A _2j is obtained. Since the sum of squares of correlation in the objective function is a quadratic equation of A _1i and A _2j , the partial differentials of the sum of squares of correlation with respect to A _1i and A _2j are linear functions of A _1i and A _2j , respectively, so an approximate solution can be obtained. Setting an initial point, iteratively updating A _1i and A _2j through formula (1), and ensuring that the value of the sum of squares of correlation increases, when the convergence criterion of the objective function is met, the iteration stops, and A _1i and A _2j at this time are taken as the optimal solution of this iteration process. According to the optimal solution obtained in each iteration process, A ₁ and A ₂ are updated, and the above steps are repeated until a stable solution is obtained, which is the optimal solution of formula (1).

The reference information is used to guide multimodal fusion, including but not limited to working memory scores, cognitive scores, symptom scores, and genes at a certain site. The reference information can be selected according to the actual research purpose.

S4: Step S3 can obtain multimodal covariant components and mixing matrices that are significantly correlated with the reference information, thereby achieving globally optimized supervised multimodal brain image fusion.

Figure 2(a) compares the absolute value of the correlation between each modal target component and the reference information estimated by each method with the absolute value of the true correlation between each modal target component and the reference information.

In this embodiment, MCCA refers to multi-group canonical correlation analysis; jICA refers to joint independent component analysis; separate ICA refers to independent component analysis of a single modality; PICA refers to parallel independent component analysis; MCCAR refers to supervised multivariate canonical correlation analysis; MCCAR+jICA refers to supervised multivariate canonical correlation analysis + joint independent component analysis; supervised PICA refers to supervised parallel independent component analysis. In addition, truecorrelation refers to true correlation.

As shown in Figure 2(a), the method proposed in the embodiment of the present invention can still extract independent components related to the reference information when the true correlation is small. As the true correlation increases, the method proposed in the embodiment of the present invention shows strong robustness in extracting independent components.

FIG2( b ) exemplarily shows the accuracy of each method relative to the mixing matrix _Ak and the independent components _Sk which are averaged over the absolute values of the correlations between the target components of each modality and the reference information.

Precision refers to the correlation between the mixing matrix and independent components obtained by each method and the true mixing matrix and components. The higher the correlation between the two, the higher the precision.

It can be seen from Figure 2(b) that the method proposed in the embodiment of the present invention is optimal in terms of the accuracy of the mixing matrix _Ak and the independent component _Sk compared to other supervised multimodal fusion methods, which illustrates the importance and necessity of using a global optimization solution method.

FIG3(a) exemplarily shows the accuracy of the averaged independent component _Sk in each method relative to 14 noise levels (peak signal-to-noise ratio PSNR = [-1110]). FIG3(b) exemplarily shows the accuracy of the averaged mixing matrix _Ak in each method relative to 14 noise levels (peak signal-to-noise ratio PSNR = [-11 10]). It can be seen from FIG3(a) and FIG3(b) that the method proposed in the embodiment of the present invention has a relatively stable ability to extract independent components under the influence of different noises, and is optimal in terms of the accuracy of the mixing matrix _Ak and the independent component _Sk .

As shown in Table 1, the comparison results of the ability of each method to extract independent components significantly correlated with reference when the peak signal-to-noise ratio PSNR = 7 are shown. Among them, real corr is the real correlation; A~A ^* is the accuracy of the mixed matrix extracted by the method proposed in the example of the present invention; S~S ^* is the accuracy of the independent components extracted by the method proposed in the example of the present invention.

Table 1 Comparison of the ability of each method to extract independent components significantly related to the reference information

The following is a detailed description of the test performed by the method proposed in the examples of the present invention.

Generation of simulated data in this example: Two data sets with the same dimensions as fMRI and sMRI data were simulated using the simTB toolbox. Each data set contained 6 brain networks with different distributions and two random mixing matrices.

The source of real data in this example: Function Biomedical Informatics Research Network (FBIRN) Phase III study data set, which contains 294 subjects (147 schizophrenia patients and 147 normal controls with matching demographic information), and the working memory scores of all subjects have been evaluated using the CMINDS cognitive assessment system. These data were used to test the method proposed in the embodiment of the present invention. The working memory score was used as reference information.

Figures 4(a), 4(b), 4(c) and 4(d) exemplarily show the results of the embodiment of the present invention on the data set FBIRN, where the IC is an independent component. Among them, Figure 4(a) is a spatial view of the multimodal covariant component related to the working memory score detected by the method proposed in the embodiment of the present invention. It can be seen that the location of the brain areas and the degree of activation contained therein are basically consistent with existing studies; Figure 4(b) is a group difference diagram of each modality; Figure 4(c) is a correlation fitting diagram of each modality and the working memory score Workingmemory; Figure 4(d) is a fitting diagram of the correlation between modalities. In summary, it can be seen that the method proposed in the present invention can separate independent components with inter-group differences, inter-modality correlations and significant correlations with working memory scores.

Claims (9)

1. A globally optimized supervised multimodal image fusion method, characterized in that it comprises the following steps: (1) Preprocessing brain image data of each modality and performing feature extraction to obtain features of each modality; (2) Preprocessing the obtained modal features to obtain the feature matrix of each modality after preprocessing; (3) The independence of the independent components of each mode calculated based on the characteristic matrix of each mode, the mixing matrix of each mode and the weight matrix of each mode; (4) Finding the optimal solution of the mixing matrix that maximizes the independence of the independent components of each mode, the correlation between the mixing matrices of each mode, and the correlation between the mixing matrices of each mode and the reference information; (5) Based on the mixing matrix obtained in step (4), calculate the covariant component and mixing matrix of the multimodal information related to the reference information. 2. The supervised multimodal image fusion method according to claim 1 is characterized in that the independence of the independent components of each modality in step (3) is described by an information entropy function H(·), and the calculation formula for maximizing the independence of the independent components of each modality is: U _k =W _k ·X _k +W _k0 ; Among them, E(.) is the expected value; is the probability density function; Y _k and U _k are the intermediate variables of the calculation; W _k is the unmixing matrix of the mixing matrix _Ak ; _Ak is the mixing matrix of mode k; W _k0 is the bias weight matrix of mode k; X _k is the feature matrix of mode k after preprocessing. 3. The supervised multimodal image fusion method according to claim 2 is characterized in that the solution formula for maximizing the independence of the independent components of each modality, the correlation between the mixing matrices of each modality, and the correlation between the mixing matrices of each modality and the reference information in step (4) is: Among them, m _k (k＝1,2) is the number of independent components of mode k; A _1i is the i-th column of the mixing matrix of mode 1; A _2j is the j-th column of the mixing matrix of mode 2; ref _h represents the h-th column of the reference information matrix; Corr(A _1i ,A _2j ) represents the correlation between A _1i and A _2j ; Corr(A _1i ,ref _k ) represents the correlation between A _1i and ref _h ; Corr(A _2j ,ref _h ) represents the correlation between A _2j and ref _h ; α ₁ , α ₂ , and α ₃ represent regularization parameters. 4. The supervised multimodal image fusion method according to claim 3 is characterized in that the gradient descent algorithm is used to iterate the solution formula until convergence to obtain the optimal solution of the mixing matrix: The iteration rules of the unmixing matrix of each mode satisfy: Where _λk is the learning rate for independent component analysis of mode k; I is the identity matrix; and T is the transpose of the matrix. 5. The supervised multimodal image fusion method according to claim 4, characterized in that the step of maximizing the correlation between the mixing matrices of each modality and the correlation between the mixing matrices of each modality and the reference information in step (4) specifically comprises: (41) Calculate the correlation between the mixing matrix of each mode and the reference information, and select the triplet A _1i , A _2j and ref _k with the largest correlation as the constrained items; (42) Taking the partial derivatives of the sum of squared correlations with respect to A _1i and A _2j respectively and setting them equal to 0, we obtain the partial differential equation containing A _1i and A _2j ; (43) Setting an initial point, iteratively updating the A _1i and A _2j , and ensuring that the value of the correlation square sum increases until the correlation square sum converges, and taking the A _1i and A _2j at this time as the optimal solution of this iteration process; (44) Update A ₁ and A ₂ according to the optimal solution in step (43), and repeat steps (41), (42) and (43) to obtain the optimal solutions A _1i and A _2j . 6. The supervised multimodal image fusion method according to claim 5, characterized in that the iterative update rule for iteratively updating A _1i and A _2j in step (43) satisfies: Among them, Std(·) is the standard deviation; Cov(·) is the covariance; var(·) is the variance; They respectively represent the average values of A _1i , A _2j , and ref _k ; η ₁₁ , η ₁₂ , η ₂₁ , and η ₂₂ are all descent steps obtained by the gradient descent algorithm; λ _c1 and λ _c2 are learning rates. 7. The supervised multimodal image fusion method according to claim 1 is characterized in that each modality image in step (1) comprises a functional magnetic resonance image (fMRI) and a structural magnetic resonance image (sMRI). 8. The supervised multimodal image fusion method according to claim 7 is characterized in that the modal features extracted from the functional magnetic resonance image (fMRI) include fractional low-frequency fluctuation amplitude; and the modal features extracted from the structural magnetic resonance image (sMRI) include gray matter volume. 9. The supervised multimodal image fusion method according to claim 1 is characterized in that the preprocessing of each modal feature obtained in step (2) includes: matrixing, centering, whitening and dimensionality reduction using a principal component analysis method.