A METHOD FOR MULTI-COMPONENT ANALYSIS ON MRI MEASUREMENT DATA

Abstract

It is an object of present invention to provide for a faster method of multi-component analysis. This object is achieved by a method for multi-component analysis on MRI measurement data, wherein a component is defined by one or more tissue component parameters among which preferably one is a T2 or T1 value. The method comprising steps of receiving the MRI measurement data, wherein the MRI measurement data comprises multiple signals corresponding to multiple voxels in an MRI image and wherein the MRI measurement data is acquired by means of a sequence encoding the one or more tissue component parameters; identifying components in the multiple voxels by minimizing a difference between the multiple signals and a linear combination of weighted simulated temporal signal evolutions, wherein different simulated temporal signal evolutions represent different components and are based on different values of the one or more tissue component parameters, and wherein the identification of the components is performed under the assumption that the possible components are the same for all of the multiple voxels and wherein a higher total number of components is penalized over a lower total number of components, and wherein the simulated temporal signal evolutions are weighted by a weight factor that is non-negative.

Claims

1. A method for multi-component analysis on MRI measurement data, wherein a component is defined by one or more tissue component parameters, the method comprising steps of receiving the MRI measurement data, wherein the MRI measurement data comprises multiple signals corresponding to multiple voxels in an MRI image and wherein the MRI measurement data is acquired by a sequence encoding the one or more tissue component parameters; identifying components in the multiple voxels by minimizing a difference between the multiple signals and a linear combination of weighted simulated temporal signal evolutions, wherein different simulated temporal signal evolutions represent different components and are based on different values of the one or more tissue component parameters, and wherein the identification of the components is performed under the assumption that the possible components are the same for all of the multiple voxels and wherein a higher total number of components is penalized over a lower total number of components, and wherein the simulated temporal signal evolutions are weighted by a weight factor that is non-negative.

2. The method according to claim 1, further comprising the step of creating a set comprising the simulated temporal signal evolutions.

3. The method according to claim 1, wherein the MRI measurement data is acquired by a multi-echo spin-echo acquisition.

4. The method according to claim 1, wherein the MRI measurement data is acquired by an MR fingerprinting sequence.

5. The method according to claim 1, wherein the components are at least one selected from the group of: myelin water, intra- and extra-cellular water, or free water.

6. The method according to claim 1, wherein the MRI measurement data is acquired with a sequence encoding for diffusion and wherein different components are identified at least partly based on diffusion values.

7. The method according to claim 1, wherein the components are at least one selected from the group of epithelium, lumen and stroma.

8. The method according to claim 1, further comprising the step of receiving a B1 map for the multiple voxels and taking into account a B1 value for the respective voxels when identifying the components.

9. The method according to claim 8, wherein the method further comprises: creating a set of simulated temporal signal evolutions for a range of B1 values and; determining a B1 value for a voxel and; based on the voxel's B1 value, selecting a part of the set of simulated temporal signal evolutions for use in the identification of the components.

10. The method according to claim 8, wherein the B1 map originates from a B1 measurement.

11. A method according to claim 1, wherein B1 values are estimated from the MRI measurement data.

12. A computer program comprising program code configured to perform the method according to claim 1.

13. A magnetic resonance imaging system comprising the computer program according to claim 12.

Description

BRIEF DESCRIPTION OF THE FIGURES

[0021] FIG. 1 diagrammatically shows a method according to embodiments of the invention and

[0022] FIG. 2 displays the flip angle variation in a possible sequence and

[0023] FIG. 3 displays the results of the multi-component analysis on the phantom data and

[0024] FIG. 4 displays the estimated component weights in in-vivo data and

[0025] FIG. 5 diagrammatically displays a magnetic resonance imaging system according to the invention.

DETAILED DESCRIPTION OF THE INVENTION

[0026] FIG. 1 diagrammatically shows a method according to embodiments of the invention. First, temporal signal evolutions are being simulated 100, wherein different simulated temporal signal evolutions represent different components and are based on different values of the one or more tissue component parameters. Depending on the tissue component parameters of interest (e.g. T1, T2, etc) and the expected values of the tissue component parameters, the temporal signal evolutions are simulated. The creation of the set of simulated temporal signal evolutions does not necessarily need to be repeated for each patient. Therefore, this creation could also be performed at a different location and/or by different means from where and by which the multi-component analysis is being performed. Preferably, temporal signal evolutions are being simulated for a variety of B1 values.

[0027] Then, in step 102, a B1 map is received. This B1 map provides B1 values related to the voxels of interest for the following multi-component analysis. The B1 map preferably is created from B1 measurement data acquired by an MM system. After that, MM measurement data is received 102. The MM measurement data comprises multiple signals corresponding to multiple voxels in an MRI image. The MRI measurement data is acquired by means of a sequence encoding the one or more tissue component parameters. This sequence could for example be a multi-echo spin-echo acquisition in case T2 would be a tissue component parameter of interest. Another option is inversion recovery multi echo spin echo, as described in Kim D, Wisnowski J L, Nguyen C T, Haldar J P. Relaxation-Relaxation Correlation Spectroscopic Imaging (RR-CSI): Leveraging the Blessings of Dimensionality to Map In Vivo Microstructure. arXiv:180605752 [eess] 2018 June Also the sequence could be an MR fingerprinting sequence. It the latter case there are multiple tissue component parameters that can be encoded within a single sequence. Then in step 106 components are identified in multiple voxels simultaneously. This is achieved by minimizing a difference between the multiple signals and a linear combination of weighted simulated temporal signal evolutions. Different simulated temporal signal evolutions represent different components and are based on different values of the one or more tissue component parameters. The identification of the components is performed under the assumption that the possible components are the same for all of the multiple voxels. A higher total number of components is penalized over a lower total number of components. The simulated temporal signal evolutions are weighted by a weight factor that is non-negative.

[0028] Below more information about the background of the method and preferred implementations will be provided.

[0029] In a multi-component signal model, the temporal signal evolution x.sub.j of a voxel j, can be written as

x.sub.j=Dc.sub.j+e.sub.j  (1)

where D is a set of simulated temporal signal evolutions, which in MRF is often called a dictionary, c.sub.j is a vector containing the weights for the different components and e.sub.j is a noise term. The measured temporal signal evolutions are generally complex, however, when the phase of the signal is determined, the signal can be rotated to the real axis resulting in a real vector x.sub.jϵ.sup.M*N of length M, where M is the number of time points of the (fingerprinting) sequence or the length of the signal after singular value decomposition (SVD) compression.

[0030] D∈.sup.M*N contains the signal evolutions for N different components. The measured temporal signal evolution is modelled as a non-negative linear combination of the simulated signal evolutions or dictionary signals. The weights of these different components are contained

[0031] in the vector c.sub.jϵ.sub.≥0.sup.M. Besides the non-negativity constraint, it can be assumed that the weight vector c.sub.j is sparse, thus the measured signal can be represented by a small number of components, representing a small number of tissue types. The weights for each component in Eq. (1) can for example be obtained by least squares minimisation. When the requirement is included that c is non-negative, the following non-negative least squares (NNLS) problem for each voxel j is obtained:

[00001]$\begin{array}{cc}\begin{array}{c}\min \\ c_j\in {ℝ}^N\end{array}{.Math.x_j-D{c}_j.Math.}_2^2& \left(2\right)\end{array}$

[0032] For a set of simulated temporal signal evolutions with a large number of components, this problem is highly under-determined and has infinitely many solutions. This formulation is similar to a compressed sensing problem. Therefore, if the solution vector is sparse, there are some theoretical guarantees that it can be recovered using a sparsity constraint. However, due to the high coherence of simulated temporal signal evolutions, a unique solution only exists for very sparse solutions. One sparsity promoting approach to solve this problem is the active set NNLS algorithm as proposed by Lawson C L, Hanson R J. Solving Least Squares Problems. SIAM; 1974, Chapter 23. The NNLS algorithm aims at finding a non-negative least squares solution.

[0033] Another approach to restrict the solution would be by means of a joint sparsity constraint. This could for example be implemented in the form of regularisation. This could for example be achieved by means of the following equation:

[00002]$\begin{array}{cc}\begin{array}{c}\min \\ c_j\in {ℝ}^N\end{array}{.Math.x_j-D{c}_j.Math.}_2^2+{\lambda }^2.Math.c_j{.Math.}_1^2& \left(2\right)\end{array}$

[0034] Where λ>0 is the regularisation parameter.

[0035] The joint sparsity constraint may be introduced by using the weights w.sub.i=∥C[i, :]∥.sub.2+ϵ(here ϵ=10.sup.−4) (i=1, . . . , N) to couple the different voxels.

[0036] There are multiple implementations foreseeable how the skilled person may exploit joint sparsity and non-negativity to arrive at a faster method for multi-component analysis. The invention is not limited to the specific examples mentioned in the description.

[0037] Other examples of how joint sparsity can be exploited are known from literature. Some examples of papers on the use of joint sparsity are provided below.

[0038] Duarte Alf, Sarvotham S, Baron D, Wakin M B, Baraniuk R G. Distributed Compressed Sensing of Jointly Sparse Signals. In: Proceedings of the 2005 Asilomar Conference on Signals, Systems, and Computers IEEE; 2005. p. 1537-1541. In this paper all the algorithms proposed are greedy and try to minimize the .sub.0 norm

[0039] Cotter S F, Rao B D, Kjersti Engan, Kreutz-Delgado K. Sparse Solutions to Linear Inverse Problems with Multiple Measurement Vectors. IEEE Transactions on Signal Processing 2005 July; 53(7):2477-2488. This paper refers to multiple measurements vector, meaning measurements with multiple sensors that are correlated and deals with enforcing the sparsity using the .sub.p norm where p≤1. This is a very common expression in sparsity enforcing algorithms, not only joint sparsity.

[0040] Tropp J A, Gilbert A C, Strauss M J. Algorithms for Simultaneous Sparse Approximation. Part I: Greedy Pursuit. Signal Processing 2006 March; 86(3):572-588.

[0041] Tropp J A. Algorithms for Simultaneous Sparse Approximation. Part II: Convex Relaxation. Signal Processing 2006 March; 86(3):589-602. Both the papers from Tropp et al. call the joint sparsity problem “simultaneous sparse approximation” and propose a greedy pursuit algorithm, which is a very different approach.

[0042] Blanchard J D, Cermak M, Hanle D, Jing Y. Greedy Algorithms for Joint Sparse Recovery. IEEE Transactions on Signal Processing 2014 April; 62(7):1694-1704 focuses on the use of greedy algorithms.

[0043] Further acceleration of the algorithm may be obtained by using SVD compression of the dictionary and temporal signal evolutions. Also pruning of the dictionary after the second iteration may accelerate the method. In the pruning all the unused dictionary components (∥C[i, :]∥=0) are removed, making the other iterations very fast. After 10 iterations the results become stable in general.

[0044] If B1 information is available for the voxels of interest, this information will be used for correcting the simulated temporal signal evolutions or the multiple measured signals 104. For example, based on a voxel's B1 value, only the part of the set of simulated temporal signal evolutions (sub-dictionary in MRF) that corresponds to this specific B1 value can be selected for use in the identification of the components. When taking into account B1 values, still it is assumed that there is only a limited number of components in the region of interest and hence only a limited number of tissue component parameter (combinations).

[0045] To test the proposed method, simulations were performed with a numerical phantom containing three different components. The relaxation times for the simulated components were chosen according to a three tissue brain model, where the measured MR signal is a combination of myelin water, intra- and extracellular water, and free water. The first component is in the range of myelin water (MW) with relaxation times (T1=67 ms and T2=13 ms), the second component in the range intra- and extracellular water (IEW) (T1=1000 ms and T2=100 ms) and the third component in the range of free water (FW) (T1=2000 ms and T2=500 ms).

[0046] Ten multi-component compositions were simulated, the first component had a weight of 10% in each composition, the other two components vary from 0% to 90%. Based on each of the ten multi-component mixtures the signal evolution was calculated, to each of these signals different independent Gaussian noise was added, resulting in a total of 10×10 simulated voxels with a signal to noise ratio of 50.

[0047] An MRF-FISP sequence was used according to Jiang Y, Ma D, Seiberlich N, Gulani V, Griswold M A. MR Fingerprinting Using Fast Imaging with Steady State Precession (FISP) with Spiral Readout: MRFingerprinting with FISP. Magnetic Resonance in Medicine 2015 December; 74(6):1621-1631. The MRF-FISP sequence of 200 time points was used for the simulations. FIG. 2 displays the flip angle variation in the sequence used and a constant repetition time of TR=15 ms. A logarithmically spaced dictionary containing 80 T1 values from 10 ms to 5 s and 80 T2 values from 10 ms to 5 s was computed with the extended phase graphs algorithm (EPG), with the restriction T2≤T1 resulting in 3240 dictionary atoms.

[0048] FIG. 3 displays the results of the multi-component analysis on the phantom data. On the left column the ground truth is being displayed. On the right side the results achieved with a method according to the invention are displayed. The root mean squared error (RMSE) is considered to be low.

[0049] Further, ten 3T measurements on healthy volunteers were performed to test a method according to embodiments of the invention. The 10 volunteers were scanned on a 3.0 T Philips Achieva scanner with Sequence 2 as given in FIG. 4, 400 with a Cartesian sampling pattern, a FOV of 240×240 mm.sup.2, in plane resolution of 1.25×1.25 mm.sup.2, and 10 mm slice thickness. The acquisition time for one slice was 337 s. The dictionary included different relative B1-inhomogeneity values, ranging from 0.75 to 1.26 with step size 0.03, leading to a dictionary size of 845580.

[0050] The average values calculated for white matter over the 10 volunteers was T1=898.1 ms and T2=53.2 ms for white matter and T1=1241.0 ms and T2=58.8 ms for grey matter. Only four components were matched by the method in the 3T measurements brain measurements. FIG. 4 displays the estimated component weights in in-vivo data 402 achieved with the flip angle variation displayed in 400. Component a and b are assumed to be related to white and grey matter respectively, the other two components are related to CSF or free water.

[0051] FIG. 5 diagrammatically displays a magnetic resonance imaging system 500 according to the invention. The method could be performed by a computer program product 502. This computer program product could in turn be comprised in a magnetic resonance imaging system. This magnetic resonance imaging system could be further configured to acquire the MM measurement data and/or the B1 map. The multi-component analysis could become part of the standard image reconstructions performed by the MM system. The MM system could further comprise a display configured for displaying the results of the multi-component analysis to a user.

[0052] The multi-component analysis could also be performed on a separate workstation or in the cloud or any other means known to the person skilled in the art.

[0053] Whilst the invention has been illustrated and described in detail in the drawings and foregoing description, such illustrations and description are to be considered illustrative or exemplary and not restrictive; the invention is not limited to the disclosed embodiments.