A three-dimensional unwrapping method for complex fMRI phase data based on quality-guided

Abstract

A kind of three-dimensional unwrapping method of complex fMRI phase data based on quality guide belongs to biomedical signal processing field.The path guide algorithm based on phase quality map of the present application carries out spatial three-dimensional unwrapping to complex fMRI phase data, and solves the problem of losing spatial dimension brain function information of current fMRI phase unwrapping method.Because phase derivative variance (PDV) map has strong robustness in identifying phase image noise area, the present application uses PDV map as phase quality map, reliably guarantees the performance of three-dimensional phase unwrapping method, and effectively excavates more complete fMRI spatial information.Taking resting state fMRI data of a healthy subject as an example, the present application method, traditional complex division method and PRELUDE method are used to respectively unwrap phase data, and form complex fMRI data by combining with amplitude data.Using complex EBM algorithm, DMN component and AUD component are respectively extracted.The present application method has obvious improvement in correlation coefficient and effective voxel number, and the advantage is more obvious under high threshold, which can provide complete fMRI evidence for brain cognition and brain disease research.

Description

Technical Field

[0001] This invention belongs to the field of biomedical signal processing and relates to a quality-guided three-dimensional unwinding method for complex fMRI phase data. Background Technology

[0002] fMRI (functional Magnetic Resonance Imaging) has been widely used in research on brain function and brain diseases due to its advantages such as high resolution, non-invasiveness, and ease of acquisition. The raw data from fMRI is complex, including both amplitude and phase. However, because phase data is susceptible to noise interference, researchers typically only analyze the amplitude data, i.e., real-valued fMRI. Nevertheless, increasing research indicates that fMRI phase data contains unique information about brain function. Effective utilization of this data can help reveal more complete information about brain function, and it holds particular potential in the study of brain diseases.

[0003] fMRI phase data are obtained by calculating the arctangent of the real and imaginary parts of complex data, with a value range of (-π, π], where phase wrapping inevitably occurs. Specifically, there are wrapped phase ψ(x, y, z) and unwrapped phase. The following relationship must be satisfied:

[0004]

[0005] In the formula, (x,y,z) represents three-dimensional coordinates, W(·) represents the winding operation, and modulo... 2π (·) indicates a modulo operation with 2π as the modulus. (See appendix) Figure 1 The example illustrates a tangled phase where the phase map exhibits discontinuous jumps at ±π. Only continuous phase information can reflect continuous changes in brain function; therefore, it is necessary to untangle the fMRI phase data, hereinafter referred to as phase untangling. Existing fMRI phase untangling methods include the complex division method. This involves dividing the data from other time points using the complex data from the first time point, and then taking the phase from the result of the complex division; the resulting untangled phase is the untangled phase. The drawback of the complex division method is that it only utilizes the temporal dimension of the fMRI phase, neglecting the three-dimensional spatial information of the fMRI phase. To avoid the loss of spatial information about brain function, it is necessary to perform three-dimensional spatial phase untangling.

[0006] In fact, three-dimensional and two-dimensional spatial phase unwrapping algorithms already exist in other fields, mainly divided into global optimization algorithms and path-guided algorithms. However, these algorithms do not consider the data characteristics and physical meaning of fMRI, so they are not suitable for fMRI phase unwrapping. For example, the global optimization algorithm PRELUDE (Phase Region Expanding Labeler for Unwrapping Discrete Estimates) is the gold standard algorithm for three-dimensional phase unwrapping in the field of quantitative magnetic susceptibility imaging. This method is implemented through global partition boundary fusion. This processing method is prone to losing brain functional information near the phase transition boundary and cannot effectively identify and process residual points in the phase map. Path-guided algorithms are mainly used for two-dimensional phase unwrapping of synthetic aperture radar images, and can be further subdivided into quality-guided methods, branching methods, etc., with superior unwrapping performance. However, fMRI phase is a three-dimensional image, and spatial information has continuity. Directly applying two-dimensional phase unwrapping algorithms will inevitably destroy the original spatial continuity of brain functional information. Summary of the Invention

[0007] To address the aforementioned issues, this invention employs a two-dimensional quality-guided method with better unwrapping performance for three-dimensional extension, forming an effective three-dimensional phase unwrapping method for fMRI. Because the Phase Derivative Variance (PDV) map is highly robust in identifying noisy regions in phase images, this invention uses the PDV map as the phase quality map, reliably ensuring the performance of the three-dimensional phase unwrapping method and effectively extracting more complete fMRI spatial information.

[0008] The technical solution adopted in this invention is as follows:

[0009] Step 1: Input the entangled fMRI phase data in, For single-subject fMRI complex data, T is the number of time points, t=1,…,T,V=X×Y×Z, X,Y,Z are the length, width and height of the fMRI input data;

[0010] Step 2: Devoxel removal from the brain. Given the amplitude data of the fMRI complex dataset. Amplitude data at the first time point mean As a threshold, a binary mask is generated. when BM(v) = 1, otherwise BM(v) = 0. The original phase data P′ is masked using BM: P = BM·P′, yielding the phase data after de-brainwashing.

[0011] Step 3: Let t = 1;

[0012] Step 4: Extract the phase data at time point t. x=1,…,X, y=1,…,Y, z=1,…,Z;

[0013] Step 5: Calculate P using formulas (2) and (3). t 3D PDV diagram of (x,y,z)

[0014] Q t (x,y,z)=Q′ t (x,y,z) / H 3 #(2)

[0015]

[0016] In the formula, H is the size of the three-dimensional window H×H×H, and (x,y,z) are the coordinates of the central voxel of the three-dimensional window. Let x, y, and z represent the local partial derivatives in the x, y, and z directions, respectively, and calculate them according to the following formula (where x, y, and z are the local partial derivatives in the x, y, and z directions, respectively). For example: i,j,k are offset variables. These are the mean values ​​of the local partial derivatives in the x, y, and z directions within the three-dimensional window, calculated using formula (4):

[0017]

[0018] At this time, phase diagram P t Each voxel in (x,y,z) corresponds to a PDV diagram Q. t A value in (x,y,z). Q t The smaller the (x,y,z) value, the higher the phase quality of the voxel; in the mass-guided method, Q is used. t The (x,y,z) values ​​guide the direction of unwrapping. Since edge voxels do not require processing, they are assigned the value max(Qt), i.e., the minimum phase quality, to avoid selecting the unwrapping start point from the edge voxels.

[0019] Step 6: Based on Q t (x,y,z) Take the voxel coordinates (x0,y0,z0) with the highest phase quality in the PDV diagram as the unwrapping starting point and generate an unwrapping queue. At this time, the queue contains only one coordinate element. Extract the original phase value P corresponding to its six neighbor coordinates (x0+1,y0,z0), (x0-1,y0,z0), (x0,y0+1,z0), (x0,y0-1,z0), (x0,y0,z0+1), (x0,y0,z0-1). t (x0+1,y0,z0),P t (x0-1,y0,z0),Pt (x0,y0+1,z0),P t (x0,y0-1,z0),P t (x0,y0,z0+1),P t (x0, y0, z0-1), perform phase unwrapping operation; with P t Taking (x0+1, y0, z0) as an example, the untangling operation is as follows:

[0020] unwrappedP t (x0+1,y0,z0)=P t (x0+1,y0,z0)+2hπ#(5)

[0021] h = mod 2π [P t (x0+1,y0,z0)-P t (x0,y0,z0)]#(6)

[0022] The PDV value Q of the six neighboring voxels of coordinates (x0, y0, z0) t (x0+1,y0,z0),Q t (x0-1,y0,z0),Q t (x0,y0+1,z0),Q t (x0, y0-1, z0), Q t (x0,y0,z0+1),Q t (x0, y0, z0-1) are sorted from smallest to largest, and the corresponding phase quality is sorted from high to low. The six coordinates corresponding to them are added to the unwrapping queue. At this time, there are seven coordinate elements in the unwrapping queue.

[0023] Step 7: Dequeue the voxel coordinate element (x0, y0, z0) with the highest quality in the current unwrapping queue, take the point with the highest phase quality in the unwrapping queue as the new unwrapping starting point, unwrap the unwrapped voxels in its six neighboring areas according to formulas (5) and (6), and sort the coordinates of the unwrapped voxels in descending order of PDV value and add them to the unwrapping queue;

[0024] Taking (x0+1, y0, z0) as a new starting point as an example, since (x0, y0, z0) does not need to be unwrapped, only five voxels P need to be unwrapped. t (x0+2,y0,z0),P t (x0+1,y0+1,z0),P t (x0+1,y0-1,z0),P t (x0+1,y0,z0+1),P t(x0+1,y0,z0-1); after unwrapping, the five voxel coordinates are sorted according to the PDV value from small to large and added to the unwrapping queue, at this time there are 11 coordinate elements in the queue; during the unwrapping process, the order of the unwrapping queue is always arranged according to the phase quality from high to low, and the length of the unwrapping queue is first gradually lengthened, and then shortened with the decrease of the six-neighborhood unwrapped voxels;

[0025] Step 8: judge whether the queue is empty, if not empty, repeat step 7, traverse all voxel points, until all brain voxels are unwrapped, that is, the queue is empty; if the queue is empty, all voxels are unwrapped, and the next step is executed;

[0026] Step 9: if t≠T, let t=t+1, repeat steps 4-8 to complete the phase unwrapping of the next time point fMRI data; if t=T, all time point phase data are unwrapped, and the next step is entered;

[0027] Step 10: output the unwrapped phase data

[0028] The method for testing the effect of phase unwrapping is to combine the unwrapped phase data and the amplitude data into the unwrapped complex data After the preprocessing steps of head motion correction, spatial standardization and spatial smoothing, X' is obtained The entropy bound minimization algorithm (EBM) of ICA is used to decompose X' into a spatial map (SM) matrix and a time course (TC) matrix K is the model order of the EBM algorithm. The K rows of S contain K spatial activation component vectors k=1,…,K; the K columns of A contain K time course component vectors k=1,…,K. According to the reference template of the spatial activation component (From the literature "S. M. Smith, P. T. Fox et al., Correspondence of the brain's functional architecture during activation and rest, Proceedings of the National Academy of Sciences of the United States of America, vol. 106, no. 31, pp. 13040-13045, 2009"), two important spatial components for distinguishing patients and healthy people, default network (DMN) and auditory network (AUD), are selected. BRIEF DESCRIPTION OF DRAWINGS

[0029] Figure 1 is a phase winding schematic.

[0030] Figure 2 is a flowchart for implementing the present application.

[0031] Figure 3 , 4, 5 are respectively the threshold value of 0.5, 1, 1.5, the comparison chart of the DMN spatial activation components extracted by the present application and other methods. A is the DMN component reference, B is the traditional complex division method, C is the PRELUDE method, and D is the present application method.

[0032] Figure 6 , 7, 8 are respectively the threshold value of 0.5, 1, 1.5, the comparison chart of the AUD spatial activation components extracted by the present application and other methods. A is the AUD component reference, B is the traditional complex division method, C is the PRELUDE method, and D is the present application method. DETAILED DESCRIPTION

[0033] One specific embodiment of the present application will be described in detail below in combination with the technical scheme.

[0034] Embodiment 1

[0035] Existing single-subject resting-state complex fMRI data There are T = 145 whole brain scans, and the intracerebral voxel V = 64x64x30 = 122880. The following are the specific steps of the three-dimensional phase unwinding and extraction of the DMN component and the AUD component of the present project:

[0036] Step 1: input the fMRI phase data with winding ​For single-subject fMRI complex data, T = 145 is the number of time points, t = 1, ..., 145, V = 122880 = 64 × 64 × 30, 64, 64, 30 are the length, width and height of the fMRI input data;

[0037] Step 2: Devoxel removal from the brain. Given the amplitude data of the fMRI complex dataset. Amplitude data at the first time point mean As a threshold, a binary mask is generated. When M (1) If (v) ≥ 4435, BM(v) = 1; otherwise, BM(v) = 0. The original phase data P′ is masked using BM: P = BM·P′, resulting in the phase data after de-brainwashing.

[0038] Step 3: Let t = 1;

[0039] Step 4: Obtain the phase data at the current time point x=1,…,64, y=1,…,64, z=1,…,30;

[0040] Step 5: Using formulas (2) and (3), and taking H=3, calculate the 3D PDV diagram at the current time point. At this time, the phase diagram P t Each voxel in (x,y,z) corresponds to a phase mass value Q. t (x,y,z), Q t (x,y,z)∈[0.018,100], since edge voxels do not need to be processed, assign the value max(Q) to the edge voxels. t =100, which is the minimum phase quality, to avoid selecting the untangling starting point from the edge voxel;

[0041] Step 6: Based on Q t (x,y,z) Take the voxel coordinates (27,42,17) with the highest phase quality in the PDV diagram as the unwrapping starting point and generate an unwrapping queue. At this time, the queue contains only one coordinate element. Extract the original phase value P corresponding to its six neighbor coordinates (28,42,17), (26,42,17), (27,43,17), (27,41,17), (27,42,18), and (27,42,16) of this coordinate. t (28,42,17), P t (26,42,17), P t (27,43,17), P t (27,41,17), P t (27,42,18), Pt (27,42,16) use formula (5) (6) to do phase unwrapping operation, the PDV value Q of the six-neighborhood voxel of coordinate (27, 42, 17) t (28,42,17), Q t (26,42,17), Q t (27,43,17), Q t (27,41,17), Q t (27,42,18), Q t (27,42,16) is sorted from small to large, and the corresponding phase quality is from high to low. The corresponding six coordinates are added to the unwrapping queue, and at this time, there are seven coordinate elements in the unwrapping queue;

[0042] Seventh step: the coordinate element (27, 42, 17) with the highest quality in the current unwrapping queue is dequeued, and the point (28, 42, 17) with the highest phase quality in the unwrapping queue is taken as a new unwrapping starting point. The six-neighborhood undetangled voxel P t (29,42,17), P t (28,43,17), P t (28,41,17), P t (28,42,18), P t (28,42,16) is unwrapped according to formula (5) (6), and the voxel coordinates after this time are added to the unwrapping queue, and at this time, there are 11 coordinate elements in the queue, which are arranged in order of the phase quality of the PDV value from high to low to update the unwrapping queue;

[0043] Eighth step: if the unwrapping queue is not empty, repeat the seventh step. If the unwrapping queue is empty, all voxels are traversed at this time to obtain the phase space unwrapped data of the t time point Execute the next step;

[0044] Ninth step: if t≠145, let t=t+1, repeat the fourth step to the eighth step to complete the phase unwrapping of the fMRI data of the next time point; if t=145, all phase data of the time point is unwrapped, and the next step is entered;

[0045] Tenth step: output the unwrapped phase data

[0046] Eleventh step: the unwrapped phase data is normalized to (-π, +π), and the amplitude data M is combined to form the complex data after unwrapping:

[0047]

[0048] Twelfth step: perform head motion correction, spatial standardization, and spatial smoothing preprocessing steps to obtain preprocessed complex data

[0049] Step 13: Decompose the preprocessed data X″ using the EBM algorithm to obtain the spatial component matrix. Time process matrix A = {a1, a2, ..., a K The order of the EBM algorithm model is K = 80;

[0050] Step 14: Select SM components and perform phase correction and phase denoising. Using the Smith2009 reference templates for DMN and AUD components, and based on the principle of maximizing the correlation coefficient with the reference components, select the DMN and AUD components, and calculate the performance indicators using equations (7) and (8) (see appendix). Figures 3-8 ).

[0051] Comparative Example 1

[0052] The results obtained by the method of the present invention are compared with those obtained by the traditional multiple division method and the PRELUDE method.

[0053] The performance metrics being compared include:

[0054] 1. With reference template The absolute value of the Pearson correlation coefficient |ρ c |:

[0055] |ρ c |=|corr(s rwf ,s k )|#(7)

[0056] |ρ c The larger the value, the more similar it is to the reference template, and the higher the accuracy of the composition.

[0057] 2. Number of activated voxels within the template, V inside :

[0058] V inside =s k ∩s ref #(8)

[0059] V inside s k Falling on s ref The number of voxels within V. inside The larger the value, the larger the activated region and the more information it contains.

[0060] Appendix Figures 3-8 The results show a comparison between the method of this invention, the traditional multiple division method, and the PRELUDE method at different thresholds. When a threshold of 0.5 is selected, the DMN components of the method of this invention (see...) Figure 3 The correlation coefficient is similar to that of other methods, but the effective voxel count V insideThe proposed method outperforms the complex division method by 103.2% (3902 vs. 1920) and the PRELUDE method by 148.5% (3902 vs. 1570); the performance of the proposed method on the AUD component (see Figure 6 ) is more significant. Compared with the traditional complex division method, the correlation coefficient is improved by 75% (0.35 vs. 0.2), and the number of effective voxels V inside is improved by 345.5% (4491 vs. 1008). Compared with the PRELUDE method, the correlation coefficient is improved by 105.8% (0.35 vs. 0.17), and the number of effective voxels V inside is improved by 63.3% (4491 vs. 2750). As can be seen from Figures 4-8 , the advantages of the proposed method over the traditional complex division method and the PRELUDE method are more obvious at the other two higher thresholds. In particular, in the DMN component extracted by the proposed method, the continuous and complete disease marker brain region, the anterior cingulate cortex (ACC), is successfully detected. The proposed method can perform high-performance unwrapping on the phase data of a single-subject original complex fMRI, and extract spatial activation components with richer information and better quality from the new complex data constructed from the unwrapped phase and original amplitude data, thereby providing new evidence for brain cognition and brain disease research.

Claims

1. A method for three-dimensional unwrapping of complex fMRI phase data based on quality-guided, characterized in that, First step: input fMRI phase data with wrapping ; wherein, is single subject fMRI complex data, is the number of time points, , , , , is the length, width, and height of the fMRI input data; Step 2: Devoxelation of the brain; amplitude data of known fMRI complex data. Using the amplitude data at the first time point mean As a threshold, a binary mask is generated. ,when , ,otherwise ; Use BM to process the raw phase data Take cover: Obtain phase data outside the brain ; Step 3: Let t = 1; Fourth step: Extracting phase data at the first t time point , ; Step 5: Using the formula to calculate the three-dimensional PDV map : (1)； (2)； where is the three-dimensional window size, is the three-dimensional window center voxel coordinate; , , denote the local partial derivatives in , , directions, computed as is the shift variable, , , are the mean values of the local partial derivatives in , , directions, computed as ​ (3)； phase map PDV map one value; the smaller the value, the higher the phase quality of the voxel; in the quality-guided method, the value is used the value is the direction of the guided unwrapping; since the edge voxels do not need to be processed, the edge voxels are assigned i.e. the lowest phase quality, to avoid selecting unwrapping starting points from the edge voxels; Step 6: Based on the coordinates of the voxel with the highest phase quality in the PDV map Generate a phase unwrapping queue with the coordinate as the starting point; the queue contains only one coordinate element at this time; extract the six-neighbor coordinates of the coordinate the corresponding original phase value , and perform phase unwrapping operation; the unwrapping operation is as follows: (4)； (5)； The coordinates The PDV values of the hexagon voxels The seven coordinates are added to the unwrapping queue in order from small to large, and the corresponding phase quality from high to low. Step 7: the element of the highest quality voxel point coordinate in the current unwinding queue is dequeued, and the point with the highest phase quality of the unwinding queue is taken as a new unwinding starting point, and the six-neighborhood undetermined voxels are unwound according to the formula , and the undetermined voxel coordinates are sorted according to the PDV value from large to small, and added to the unwinding queue. ​ Step 8: judging whether the queue is empty, if not, repeating step 7, traversing all voxel points until all brain voxels are unwrapped, i.e.the queue is empty; if the queue is empty, all voxels are unwrapped, and the next step is performed; Step 9: If , let , repeat steps 4-8 to complete the phase unwrapping of the next time point fMRI data; if , then all time point phase data are unwrapped, and proceed to the next step. Step 10: Output unwrapped phase data .

2. A method for three-dimensional unwrapping of complex fMRI phase data based on quality-guided, as claimed in claim 1, wherein, Selection of parameters in the calculation of PDV map, Using formula calculate 3D PDV diagram : (1)； (2)； In the formula is a three-dimensional window of size, taken .

Images