WO2024231819A1 - Off-resonant encoded analytical parameter quantification using multi-dimensional linearised equations - Google Patents

Abstract

A magnetic resonance imaging method is proposed for quantification of one or more system parameters The method comprises: performing (11, 12) a bSSFP data acquisition of an object of the system to obtain a plurality of image volumes, a respective image volume corresponding to an effective radio frequency phase increment; obtaining (13) voxel-wise a set of multi-dimensional bSSFP signal values for a set of voxels of the plurality of image volumes with a corresponding radio frequency phase increment; constructing (14) voxel-wise a set of multi-dimensional bSSFP profiles for the set of voxels from the set of multi- dimensional bSSFP signal values; applying (15) voxel-wise a mathematical transformation on the multi-dimensional bSSFP profiles to obtain a set of multi-dimensional characteristic coefficients or a set of approximated multi-dimensional characteristic coefficients for the set of voxels; using (16) voxel-wise the set of multi-dimensional characteristic coefficients or approximated multi-dimensional characteristic coefficients to determine a set of parameter bases for the set of voxels; applying (17) one or more analytical solution functions onto the set of parameter bases to quantify voxel-wise the one or more system parameters.

Description

OFF-RESONANT ENCODED ANALYTICAL PARAMETER QUANTIFICATION USING MULTI-DIMENSIONAL LINEARISED EQUATIONS

FIELD OF THE INVENTION

The present invention relates to a method to decode balanced steady-state free precession (bSSFP) profiles analytically for multi-parametric quantitative mapping. The bSSFP solution derived from the Bloch equations is according to the present invention reshaped using mathematical rigor definitions and identities, formulating the bSSFP equation in a complex Fourier series. Then a Fourier transform or another mathematical transformation is used to exploit the Fourier series expression for analytical parameter quantification. This complete and compact representation of bSSFP enables rapid and robust quantification of multiple magnetic resonance imaging (MRI) parameters, including the 1 , T₂ (longitudinal and transverse relaxation times, respectively), proton density, and B_o inhomogeneity.

BACKGROUND OF THE INVENTION

Quantitative MRI has initiated a significant improvement of tissue characterisation and diagnostic potential for clinical applications. Quantitative relaxometry as well as quantitative susceptibility imaging (employing quantitative B_o maps) are successfully utilised for early detection of pathological changes and as a diagnostic tool for multiple sclerosis, tumours, Alzheimer’s and Parkinson’s disease. Simultaneous quantification of multiple tissue parameters leads to spatially coherent resolved quantitative maps as well as to a significantly accelerated and more efficient data acquisition compared to a sequential application of singleparameter mapping techniques.

For simultaneous multi-parameter quantification, magnetic resonance fingerprinting (MRF) methods can be used, and they have been intensively investigated recently. While MRF relies on dictionary-based reconstruction methods, MIRACLE (“Motion-Insensitive Rapid Configuration Relaxometry” by Damien Nguyen et al., Magnetic Resonance in Medicine 00:00-00 (2016)) or PLANET (“PLANET: An Ellipse Fitting Approach for Simultaneous T and T₂ Mapping Using Phase-Cycled Balanced Steady-State Free Precession” by Yulia Shcherbakova et al., Magnetic Resonance in Medicine 79:711-722 (2018)) are based on iterative analytical model reconstruction for the simultaneous quantification of 1 and T₂ times using the phase-cycled (PC) balanced steady-state free precession (bSSFP) sequences.

PC-bSSFP, which refers to multiple acquisitions with different linear phase increments of the radio frequency (RF) pulse, provides images with the highest signal-to-noise ratio per unit time. While relaxometry usually exploits the exponential decay of the magnitude data to extract the relaxation parameters, for PC-bSSFP in particular, the signal phase is important for the simultaneous assessment of T and T₂ parameters and is an indispensable source of information.

Recent methods for simultaneous T and T₂ quantification using PC-bSSFP dismiss valuable parts of the phase information, leading to an artificially induced instability in the model fit, or they hamper the extension from simultaneous 1 and T₂ quantification to additional parameter quantifications, such as (banding free) proton density and robust magnetic field inhomogeneity (B_o) quantification.

BRIEF DESCRIPTION OF THE INVENTION

The present invention aims to overcome at least some of the above-identified problems. More specifically, the present invention proposes a novel approach to quantify multiple parameters using phase-cycled bSSFP profiles.

According to a first aspect of the present invention, there is provided a magnetic resonance imaging method for quantification of one or more system parameters for a singlecompartment and/or a multi-compartment system as recited in claim 1.

The proposed method, which according to one example embodiment can be referred to as an Off-Resonant encoded Analytical parameter quantification using Complex Linearised Equations (ORACLE), decodes bSSFP profiles analytically, by using the entirety of encoded complex-valued information. ORACLE allows for the quantification of 1 , T₂, T₂ (which in the present description is called “effective transverse relaxation time”), proton density (M_o), diffusion, magnetic field inhomogeneity (B_o) of single-compartment or multi-compartment systems and/or inhomogeneities of the excitation field

as well as a clear visual separation of tissues by chemical shift. Additionally, ORACLE allows for the quantification of proton density fractions of multi-compartment systems. By exploiting the analytical solution for PC-bSSFP signal equation it replaces lengthy iterative algorithms or limited dictionary-sized reconstruction with ultra-rapid and robust parameter quantification. The proposed method makes no assumptions, and no information is dismissed.

According to a second aspect of the present invention, there is provided a computer program product comprising instructions for implementing the steps of the method when loaded and run on a computing apparatus or an electronic device.

According to a third aspect of the present invention, there is provided a magnetic resonance imaging apparatus for quantification of one or more system parameters for a single-compartment and/or a multi-compartment system as recited in claim 15. The apparatus is thus configured to implement the method according to the first aspect of the present invention.

Other aspects of the invention are recited in the dependent claims attached hereto.

BRIEF DESCRIPTION OF THE DRAWINGS

The invention will now be described in more detail with reference to the attached drawings, in which:

• Figure 1 a block diagram schematically illustrating an example MRI system where the teachings of the present invention can be applied;

• Figure 2 illustrates a geometric interpretation of a PC-bSSFP signal equation and duality between magnetisation space and Fourier space for 18 phase cycles and 4 modes, with the definition ₂ = A >TE + < >_rest;

• Figure 3 shows a flowchart illustrating the proposed MRI method for quantification of one or more system parameters of a single-compartment system and/or a multi-compartment system according to an example embodiment of the present invention;

• Figure 4 shows a flowchart of 1 , T₂, M₀,and Af quantification using the proposed ORACLE method;

• Figures 5a to 5c show first test results for the proposed ORACLE method;

• Figures 6a to 6c show second test results for the proposed ORACLE method; and

• Figure 7 shows quantitative maps obtained in 4 healthy subjects at 3 T (V1-V4) for the proposed ORACLE method.

DETAILED DESCRIPTION OF THE INVENTION

It should be noted that the figures are provided merely as an aid to understanding the principles underlying the invention and should not be taken as limiting the scope of protection sought. As used herein, unless otherwise specified the use of the ordinal adjectives ’’first”, ’’second”, “third”, etc. to describe a common object, merely indicate that different instances of like or different objects are being referred to, and are not intended to imply that the objects so described must be in a given sequence, either temporally, spatially, in ranking, or in any other manner. As utilised herein, “and/or” means any one or more of the items in the list joined by “and/or”. As an example, “x and/or y” means any element of the three-element set {(x), (y), (x, y)}. In other words, “x and/or y” means “one or both of x and y.” As another example, “x, y, and/or z” means any element of the seven-element set {(x), (y), (z), (x, y), (x, z), (y, z), (x, y, z)}. In other words, “x, y and/or z” means “one or more of x, y, and z.” Furthermore, the term “comprise” is used herein as an open-ended term. This means that the object encompasses all the elements listed, but may also include additional, unnamed elements. Thus, the word “comprise” is interpreted by the broader meaning “include”, “contain” or “comprehend”.

Figure 1 shows a block diagram of an MRI system or apparatus 1 schematically illustrating some functional elements that are useful for understanding the teachings of the present invention. The system 1 comprises a host module or unit 2, which is configured to interact with the user 3. For this purpose, the host module comprises a user interface through which the user is able to access one or more protocols and edit them. The user can also start a scan through the host module 2. After the scan finishes, the user 3 receives image data saved in an image database in the host module. Once the scan is started, the protocols are fed into a measurement system 4. Here the one or more protocols are translated into the scanner language. The measurement system is further configured to control the operation of an MRI scanner 5, which is configured to receive its input data from the measurement system 4 to perform the scan, and which is further configured to communicate with the measurement system to inform it if for example specific absorption rate (SAR) limits are exceeded, and/or send any error messages, etc. According to the present invention, the MRI scanner 5 is configured to perform a bSSFP data acquisition to obtain a plurality of image volumes, each volume comprising a set of voxels. A respective image volume corresponds to an effective radio frequency phase increment. The scanner 5 is configured to send any collected data to a reconstruction system 6, either directly or through the measurement system 4. The reconstruction system 6 is configured to implement reconstruction methods. For the calculation process in the reconstruction system, a calculation unit 7 may be used. The calculation unit is in this example configured to also implement the proposed ORACLE method according to the present invention. Instead of being a stand-alone unit, the calculation unit 7 could instead be part of the reconstruction system. Once the data have been reconstructed in a certain way, the data are transferred to the host module 2. Here the user can access the acquired data either as raw data or as reconstructed image data, reconstructed 1 maps, reconstructed T₂ maps, reconstructed B_o maps, etc. and visualise them. It is to be noted that one or more elements shown in Figure 1 could be merged to form a multifunctional module or unit. Furthermore, data communication between various elements can take place wirelessly or through wired connections. The proposed ORACLE method is next explained in more detail. The steady-state signal equation of the PC-bSSFP sequence yields:

The signal equation depends on the proton density M_o, the repetition time TR, the echo —TR time TE, the relaxation Ej = e ^Tj for (j = 1 and 2), the imaginary identity i = V-l, the radio frequency (RF) excitation angle a, the linear phase cycle increment of the RF pulse p, and the off-resonance with respect to the reference frequency

= - (AB₀ - S_csB₀). Every (Larmor) frequency of the respective proton is said to be off-resonant when they deviate from the reference frequency of the MRI scanner (on-resonant frequency). Those protons with deviating frequencies are in off-resonance. This leads to a linear shift of the bSSFP magnitude profile with respect to the corresponding phase cycle increments as well as a global rotation of the bSSFP profile in complex plane. B_o describes the main magnetic field strength and Z1B_O its local deviations, y describes the gyromagnetic ratio (of ¹H hydrogen protons), and <5_CS the chemical shift for the respective type of magnetic group and/or substance. 1 and T₂ are longitudinal and transverse relaxation times, respectively. _rest is a parameter, bundling different effects on the absolute phase, e.g., eddy currents and arbitrary absolute phase reference of coil elements.

We define Fourier modes of PC-bSSFP profiles in the following. Typically, each periodic and differentiable function can be represented as an infinite discrete sum of higher harmonics. Because the signal equation of PC-bSSFP (i.e., Equation 1) is periodic for integer multiples of the phase cycle increments p p + 2nn in the complex space, its complex Fourier representation yields:

where c_n e <C is the n^th mode i.e., the n^th higher harmonic. M₊(<p_k) is in the magnetisation space and represents one point of the complex-valued bSSFP profile. The complex bSSFP profile for a respective voxel is an array or a collection or sequence of signal values of the same or corresponding voxels across all the image volumes (or across a set of image volumes). In the present description, a voxel refers to a data point having a certain volume (e.g., a rectangular cuboids). An image is the arrangement of all voxels (e.g., tiny rectangular cuboids) due to their location. A two-dimensional (2D) image can be visualised as the arrangement of voxels (rectangular cuboids) in a matrix grid (similar to a rectangular chess board) with finite slice thickness. This shapes a 2D image volume. A three-dimensional (3D) image can be visualised as a stack of those 2D matrix grids in the third dimension, resulting in image volumes containing at least two slices. Thus, in the present description, voxels cover data points of both 2D and 3D objects, i.e., 2D images and 3D volumes. To assess the modes of the periodic signal by an ordinary inverse Fourier transform, it is required to sample the linear phase cycle increments p_k linearly and equidistantly with e.g., p_k = [0,2TT[. This results 271 in the parametrisation of p_k = — (k - 1), with N_PC the number of sampled phase cycle Mp_C increments and k = [1, 7V_PC] e N. The modes (which in this example are complex values) can then be extracted via the discrete Fourier transform, and in this case in particular via the inverse discrete Fourier transform:

The asymptotic limit between discrete Fourier transforms and Fourier series for PC- bSSFP is explained next. Inserting the Fourier series of Equation 2 into the Fourier transform of Equation 3 yields the following result:

For big enough N_PC (approximately N_PC > 7) the term

= N_PC8_{m n} approximates the Kronecker-delta 8_{m n}. The modes of the Fourier series are approximated with higher N_PC increasingly well c_m = c_m using a discrete ordinary Fourier transform along the phase cycle dimension.

We will next explain simultaneous and analytical assessment of AB₀, M_o, T and T₂. In the following, the mathematical steps to obtain an analytical framework for ZIB₀, M_o, T and T₂ quantification are briefly discussed. Equation 1 can be analytically reshaped by exploiting geometric series and the Euler formula in a combinatorial fashion. Then the PC-bSSFP signal can be formulated as a Fourier series such that its symbolic representation is shaped coherently to the structure of Equation 2. A complete and irreducible analytical equation for all the complex Fourier series modes of a PC-bSSFP signal can be found. Since the lowest harmonics i.e., the central modes of the complex PC-bSSFP spectrum {c_₁₍ c₀, c }, exhibit the highest signal as well the best asymptotical approximation, they are used. This is similar to the duality between the k-space and the image space, see Figure 2, where most of the information and signal is encoded in the centre of k-space. Hence, selecting the central modes is similar to a wavelet transformation, filtering noise contributions, which are encoded at higher frequencies. Equations relating to ORACLE and in particular relating to a single-compartment scenario are explained next. Equation 1 can be reshaped because it is a superimposed signal of many different echo trains, which are described with the help of a partition method and phase graph formalism. In order to decompose the individual echo contributions to the emergent bSSFP signal, the first step aims to perform a Fourier series expansion. The so- obtained Fourier coefficients (or modes) are the individual components of the contributing echoes. This decomposition allows ambiguities to be avoided in contrast to the superimposed PC-bSSFP signal, and it leads to better classification of each signal contribution, which leads to higher flexibilities and more potential post-processing methods and far easier quantification itself. This increased flexibility and method potential increase dramatically the chances of finding very easy and robust equations by this expansion, because simply more options of combinations exist. It is a combinatorial approach, aiming to find simple and unambiguous linear combinations for the solution of signal profile by increasing the mathematical flexibility without loss of accuracy and without any approximation. A series expansion is in general a change in the mathematical representation of the same signal evolution, for which a solution might be found easier. Examples are the series expansion of exponential and trigonometric functions, coordinate transformation from cartesian to radial which is beneficial for example if a system exhibits radial symmetry (e.g., polar coordinates are better for describing circles than cartesian coordinates). Additionally, a Fourier series expansion is beneficial because 1) the PC-bSSFP signal profiles are periodic and hence exhibit periodic symmetry which is exploited, 2) only elementary operations are used which are computationally very robust and not affected by numerical errors 3) the function is mapped to trigonometric functions which are very well understood and are numerically robust 4) Fourier transform is linear, leading to no error propagation/amplification and are very well understood. The derivation of Fourier coefficients (or modes) is done by exploiting mathematical identities and definitions like geometrical series, generating functions and relations between real and complex discrete Fourier series.

The irreducible configuration space basis in ORACLE representation for all PC-bSSFP systems yields:

with A the free induction decay (FID) and D the echo amplitude of the PC-bSSFP phase-graph configurations. Equation 5 may thus be considered as a simplified version of Equation 1 and is in the magnetisation space. All echo paths (FIDs and echoes) are falling into one single superimposed signal path, due to the boundary condition of bSSFP sequences, i.e. , balanced gradients, z (z being its complex conjugate) is an (complex) attenuation factor for the n^th refocused echo path. Imposing complexity of z e C is an efficient mathematical description and compression of information, which simultaneously describes both phase evolution (transversal rotation) and magnitude attenuation (|z| < 1) for the respective PC-bSSFP steady state echo path. This is corroborated by the £)-weighted configuration echo sum of Equation 5. The sum of the echo configurations starts at n = 1, i.e., echoes are at least refocused once by an RF refocusing pulse. The sum of the -weighted FID configurations starts at n = 0 since an FID can also be generated by applying one single excitation pulse without applying any refocusing pulses. The explicit bSSFP FID and echo factors of the configuration space yield:

2 is significantly influenced by the reference phase and hence it is just a

E1E2 -l+(Fi-Fz)cos(a)’ correction parameter for rotations of the complex plane due to references, eddy currents, etc. P_o and q are only necessary for the derivation of some solution functions but not for the parameter quantification itself. By applying the discrete Fourier transform given in Equation 3, the above Fourier series expansion given in Equation 5 transforms linearly into its dual Fourier space with the central Fourier modes (defining a system of equations):

C-L = Az, c₀ = A, c_ = Dz. (8)

By using basic arithmetic, every ORACLE parameter can be determined via (complex) division:

In the above example, A, z and D collectively form a parameter basis for solution functions as explained below. In this example, the equations to obtain these parameters are the same for each voxel, but the parameter values (complex values) typically vary from voxel to voxel. After the above ORACLE parameters are quantified, the next step aims to map the parameters to the desired sequence-independent system parameters, i.e. , biomarkers, {A,D,z} -> {T₁₍ T₂,B₀, M_O} via analytical solution functions. An elegant approach is proposed by taking the ratio f of the FID and echo factors according to Equations 6 and 7 to derive the analytical T₂ solution function:

Solving for E₂ and T₂ requires only elementary operations and yields:

resulting in a computationally ultra-rapid and easy T₂ quantification without loss of accuracy, compared to iterative or least-square-fit algorithms.

The derivation of the T solution function exploits the relation between r and q in dependence of E₂ and E

Derivations for the M_o and B_o solution functions are performed analogously by exploiting the definition of the ORACLE parameters:

cotg), (13) 1 = 4(z) = mod_27I(y S_o + 3_csB₀)TR). (14)

As explained above, the analytical solution functions use one or more reference values, which may be at least one of the following: one or more repetition times, one or more radio frequency excitation angles, a gyromagnetic ratio, a magnetic field strength, one or more reference frequencies, one or more chemical shifts, one or more nuclear magnetic resonance spectra, one or more gradient strengths, one or more gradient timings, and one or more gradient lengths. Furthermore, the analytical solution functions may be formulated by using at least one of the following elementary mathematical operations: equalling, addition, subtraction, multiplication, division, square root operation, exponential function, logarithmic function, trigonometric function, atan2 function, complex conjugation and absolute value operation of one or more complex-valued numbers. The operations are optionally performed more than once for a respective analytical solution function, and/or the operations use one or more parameters from a respective parameter basis and one or more reference values.

It is to be noted that the solution functions of the ORACLE formalism are by far not restricted to Fourier transforms. The formalism opens also the possibility to perform a leastsquares fitting in the dual magnetisation space. Here Equation 5 can be reshaped by exploiting geometrical series:

This proves that ORACLE describes the duality between magnetisation space and Fourier space inherently. This leads to higher flexibility. A slower and less convex leastsquares fitting approach leads to additional application possibilities like parameter quantification below the asymptotic limit of the Fourier transform, banding artifact removal with minimum three PC-bSSFP data points and more. Independently of whether A,D and z are determined with Fourier modes or with a least-squares fitting, the solution functions for the magnetisation space are 100% equal to the Fourier space solution functions. If it is desired to perform a least-squares fitting in the magnetisation space, then it is convenient to do a Fourier transform estimation for the start value problem. The above demonstrates the universality of the ORACLE formalism with abundant possibilities for the method development.

ORACLE extension to multi-compartment systems is next explained in more detail. Multi-compartment systems in MRI are systems in which more than one compartment/substance is present. Typically, in clinical MRI, the most abundant multicompartment system is a water-fat compartment. Since water is characterised by different T_lt T₂ and chemical shift than fatty tissues, this leads to a different signal behaviour of fat than water according to Equation 1. In general, if two substances with different 1 , T₂ and/or chemical shift are present in a system, it is called a multi-compartment system.

In multi-compartment systems, the signal of both substance x and substance y are superimposed. This is because electromagnetic signals (recorded by the MRI scanner due to induction) obey superposition principle. The superimposed signal for a two-compartment system can be written as:

while signal profiles for a multi-compartment system are obtained from the following equation:

where f is the (proton density) fraction of substance y. For instance, f = 0.4 means that 40%

-TR —TR of substance y and 60% of substance x is present in the system. E_{l x} = e^{Tl x} and E_2>x = e^{T2 x} are longitudinal and transverse relaxation parameters and 5_{CS X} is the chemical shift of substance x. Analogous considerations apply to substance y. Explicitly for Equation 5, the superimposed signal for two compartments looks like:

It is to be noted that the present invention in the present embodiment uses superposition principle to model signals of multi-compartment systems. In mathematical detail every electromagnetic signal obeys superposition principle. This is a fundamental property in physics (described by Maxwell equations). Thus, equation M₊

e^±m(pk + De link cⁿ e^=Fm*’^fc) employs superposition principle. Physically it is impossible to model (electro-magnetic) signals of multi-compartment systems without the use of superposition principle. Single-compartment systems are obeying superposition principle as well, but since all other signal sources are zero it does not change the single-compartment signal itself. In other words, a+0+0+0+0=a.

For convenience, the fraction is substituted into A_x,D_x,A_y and D_y, which then yields:

Multi-compartment systems for the bSSFP sequence lead typically to so-called asymmetries, which have until now been an unsolved problem. However, the ORACLE representation, in particular the complex parametrisation, enables the decomposition of multicompartments, where asymmetries are no longer an obstacle.

Equation 17 is dependent on six unknown parameters i.e., {A_x,D_x,z_x,A_y,D_y,z_y . These six parameters can be determined by calculating also six modes via Fourier transform of the sampled PC-bSSFP profile:

D ^y z^y~² ■ (18)

This non-linear system of equations can be solved for (A_x, D_x, z_x, A_y, D_y, z_y) analytically by the application of the so-called Buchenberger algorithm. Once (A_x, D_x,z_x,A_y, D_y,z_y) are obtained, they can be plugged into the single-compartment solution functions of Equations 11- to 14 for substance x and y independently. This represents a straightforward map or projection from multi-compartment systems to single-compartment systems. Then the additional (protondensity) fraction can be obtained as follows:

We explain next a multiplet special case for the multi-compartment system. If the substance, e.g., fat, has not only one chemical shift but several chemical shifts, i.e., it is a multiplet system, then the system can be further extended to incorporate the spectroscopic information of the fat spectrum:

The spectral information is completely summarised as effective (complex) weighting factors in the following form:

where

In the above equation, <5_cs,_y (n) is the chemical shift of each spectral peak of substance y and h_n is the corresponding (normalised) spectroscopic amplitude. Substance x is typically water (for clinical MRI) and hence has only one spectral peak i.e. , no spectral weighting factors are necessary for substance x. The analytical solution of the ORACLE parameter of Equation 19 is analogously performed via the Buchenberger algorithm (e.g., with the software package Mathematica®).

The special case of same relaxation is next explained for the multi-compartment system. If approximations are performed, like substance x and y are only different for the chemical shift, but both relaxation parameters are mutually substantially the same, the system of equations can be further reduced to a mathematical system of four modes:

where z denotes the complex conjugation of z and Acp = y(<5_cs,y ^> ^cs,x)B₀TR. This system is for instance suitable for water-acetone systems, where 1 and T₂ times are sufficiently equal and the chemical shift of acetone with respect to water is known (also y, B₀ and TR are always known). Then the fraction can be obtained (additionally to the above-mentioned system parameters) as follows:

Conclusively since ORACLE forms a complete, fundamental, irreducible and simple formalism, it enables the until now unsolved asymmetry problem to be solved to quantify proton density fractions. This is of particular interest for fat fraction quantification, to enable the detection of pathologies, such as fatty liver diseases, atherosclerosis, and the corresponding medical diagnosis, which helps in the decision management for the treatment of the disease.

The flow chart of Figure 3 summarises the above-described magnetic resonance imaging post-processing method for the quantification of one or more system parameters of a single-compartment system and/or a multi-compartment system. In the present example, by the word system is understood a tissue of a patient under examination.

In step 11 , MRI data acquisition is carried out by sampling the system response in k- space and by applying a radio frequency phase increment. In step 12, the acquired k-space data are transformed into image space to obtain a set or sequence of image volumes, which in this example are three-dimensional image volumes and consisting of voxels. In step 13, complex signal values are obtained voxel-wise for the voxels of the image volumes, i.e., voxel by voxel or individually or separately for different voxels, with a corresponding radio frequency phase increment. In step 14, complex-valued bSSFP profiles are obtained voxel-wise from the complex signal values obtained in step 13. More specifically, a respective bSSFP profile is obtained as an array or a collection of complex signal values obtained for the same (i.e., corresponding or matching in space) voxel position across the set of image volumes. Thus, a complex-valued bSSFP profile is obtained per voxel and per image volume. In step 15, a discrete Fourier transform (or inverse discrete Fourier transform) or another mathematical transformation more broadly, such as a Laplace transform, is applied voxel-wise onto the bSSFP profiles to obtain Fourier modes or characteristic coefficients more broadly. In step 16, the Fourier modes, which are complex values, are used vowel-wise to determine parameter bases for the voxels of the image volumes. The respective parameter basis comprises or consists of a set of complex values. Multi-variant polynomials based on monomials of integer power (c_m = c=i ^Ac ^zc^m> while complex-conjugated with a negative mode number being obtained as c_m* = ' c=₁ D_c* z_c ^m) are used to determine the parameter bases. The number of monomials is at most four for the single-compartment system and at most eight for the multicompartment system for the respective parameter basis. In step 17, analytical solution functions are applied voxel-wise to the parameter bases to quantify voxel-wise system parameters. In step 18, parameter maps are generated from the quantified system parameters. More specifically, a respective parameter map is generated for a respective image volume, where the respective parameter map is a collection, or an array of quantified system parameters obtained for the respective image volume.

In the above process, steps 12 to 18 may be implemented by the calculation unit 7, while step 11 may be carried out by the scanner 5. The aim of steps 11 to 14 is to obtain bSSFP profiles for the voxels. A bSSFP profile may be understood to be a system response of a bSSFP sequence to a given radio frequency phase increment. The system response refers to the underlying tissue properties (or system parameters) from an MRI point of view, which may be microstructure, 1 , T₂, T₂ , B_o, diffusion, fraction, etc. Depending on the implementation, the input of ORACLE may be formed by the set of image volumes, a set of complex signal values or a set of complex-valued bSSFP profiles. The analytical ORACLE solution functions are important components to ORACLE and contribute to the novelty of the present invention. Furthermore, the analytical ORACLE solution functions depend on one or more data acquisition parameters (i.e. , reference values for the solution functions) and also on the ORACLE parameter basis.

Figure 4 shows a flowchart or pipeline of 1 , T₂, M_o, and Af quantification using ORACLE and where output maps are also visualised. A) A plurality of bSSFP images are acquired by using RF phase cycle increments p = [0,2TT[ for the image acquisition. B) The bSSFP signal m₊ for each voxel are assigned to the respective RF phase cycle increment to obtain a voxel-wise bSSFP profile. C) A discrete Fourier transformation (DFT), or another suitable transformation, is applied on the bSSFP profile m₊(<p) in dependence of the sampled ^-values to obtain (bSSFP) modes c_n. For low numbers of sampled RF phase cycle increments N the ^-periodic (DFT) modes c_n may exhibit aliasing effects (N = 5). Aliasing of a signal occurs if the sampling rate employed to sample a certain signal is not high enough to determine all its frequency or Fourier coefficients in an accurate or unambiguous way (Nyquist theorem). Since we know that we are acquiring bSSFP signals, we can exploit this a priori knowledge to correct for any aliasing effects, which would not be possible if we would not have a priori knowledge of the origin of the signal. D) If aliasing occurs it can be completely corrected by mapping DFT modes to bSSFP modes using an exact correction formula. DFT and bSSFP modes are sufficiently equal if the profile was not undersampled (N = 20). E) For simplicity, auxiliary definitions can be calculated. F) Auxiliary definitions and calculated bSSFP modes are inserted into the analytical solution function. The analytical solutions are encoding the complete bSSFP-profile information without the necessity of any iterative methods. G) Using steps A-F for each voxel delivers aliasing free and coregistered 1 , T₂, M_o, complex

T sum, — and 6 = 2n fTR maps in ultra-rapid reconstruction time. The term “coregistered” in T2 this context means that two different images are perfectly or substantially perfectly matched by the location of the pixels/voxels. In fact, this is a major advantage of bSSFP images.

Figure 5a shows quantitative maps that were obtained in 4 healthy subjects at 3 T (V1- V4). 1 maps were obtained with MP2RAGE and phase-cycled bSSFP data. The proton density (PD) and T /T₂ maps were obtained by utilising the ORACLE method. The squares in the first row indicate the selected ROIs for local WM/GM/CSF value determination. A = T !T₂ values in the human brain decreases more than 550% from ambient tissue to cerebrospinal fluid (CSF), visible as dark regions in the 4^th row. This is an indication that A maps can be employed for CSF suppression or as a CSF filter, which is an important feature in neuroimaging for e.g. lesion detection. Figure 5b shows histograms for ORACLE and MP2RAGE 1 values. Figure 6c shows the determined peak values visualised as a Bland- Altman (BA) plot across all volunteers. The differences are calculated by taking the difference between the reference method and bSSFP. The straight line indicates the mean difference (bias) and the dashed lines the mean difference plus/minus the 1.96-fold standard deviation. The coefficient of variation, defined as the ratio of standard deviation and mean value of the BA plot, of CV=2.9% indicated a low variability and a high precision among the 10 volunteers.

Figure 6a shows quantitative maps obtained in 4 healthy subjects at 3 T (V1-V4). T₂ maps were obtained by using multi-echo spin-echo (ME-SE) and bSSFP data. Additionally, the complex sum (CS) images were obtained based on phase-cycled bSSFP data. Figure 6b shows histograms for ORACLE and ME-SE T₂ values. Figure 6c shows the determined peak values visualised as a Bland-Altman plot across all volunteers. The differences are calculated by taking the difference between the reference method and bSSFP. The straight line indicates the mean difference (bias) and the dashed lines the mean difference plus/minus the 1.96-fold standard deviation. The coefficient of variation, defined as the ratio of standard deviation and mean value of the BA plot, of CV=3.9% indicated a low variability and a high precision among the 10 volunteers.

Figure 7 shows quantitative maps obtained in 4 healthy subjects at 3 T (V1-V4). Off- resonance maps Af were obtained with phase cycled bSSFP and dual-echo gradient-echo data. The differences of the off-resonance maps are visualised in the 3^rd row. The mean absolute difference between both maps is 2.3 Hz with a standard deviation of 3.9 Hz among all 10 volunteers.

In view of the above, one aspect of the present invention relates to a magnetic resonance imaging post-processing method for quantification of one or more system parameters of a single-compartment system and/or a multi-compartment system. The method comprises the steps of: receiving data relating to a balanced steady-state free precession (bSSFP) data acquisition of an object of the system to generate a plurality of image volumes consisting of voxels or data points, a respective image volume corresponding to an effective radio frequency phase increment; obtaining voxel-wise a set of multi-dimensional bSSFP signal values for a set of voxels of the plurality of image volumes with a corresponding radio frequency phase increment; constructing voxel-wise a set of multi-dimensional bSSFP profiles for the set of voxels from the set of multi-dimensional bSSFP signal values; applying voxel-wise a mathematical transformation on the multi-dimensional bSSFP profiles to obtain a set of multi-dimensional characteristic coefficients or a set of approximated multi-dimensional characteristic coefficients for the set of voxels; using voxel-wise the set of multi-dimensional characteristic coefficients or approximated multi-dimensional characteristic coefficients to determine a set of parameter bases for the set of voxels; and applying one or more analytical solution functions onto the set of parameter bases to quantify voxel-wise the one or more system parameters.

At least some of the method steps can be considered as computer-implemented steps. The invention thus also relates to a non-transitory computer program product comprising instructions for implementing the steps or at least some of the steps of the method when loaded and run on computing means of a computing device, such as the calculation unit 7.

While the invention has been illustrated and described in detail in the drawings and foregoing description, such illustration and description are to be considered illustrative or exemplary and not restrictive, the invention being not limited to the disclosed embodiment. Other embodiments and variants are understood, and can be achieved by those skilled in the art when carrying out the claimed invention, based on a study of the drawings, the disclosure and the appended claims. New embodiments may be obtained by combining any of the teachings above. For example, all the features of the invention described in connection with the method can be used to characterise the apparatus or system that can implement the method.

In the claims, the word “comprising” does not exclude other elements or steps, and the indefinite article “a” or “an” does not exclude a plurality. The mere fact that different features are recited in mutually different dependent claims does not indicate that a combination of these features cannot be advantageously used. Any reference signs in the claims should not be construed as limiting the scope of the invention.

Claims

1. A magnetic resonance imaging method for quantification of one or more system parameters of a single-compartment system and/or a multi-compartment system, the method comprising: performing (11 , 12) a balanced steady-state free precession (bSSFP) data acquisition of an object of the system to obtain a plurality of image volumes, each volume comprising a set of voxels, a respective image volume corresponding to an effective radio frequency phase increment; obtaining (13) voxel-wise sets of multi-dimensional bSSFP signal values for the sets of voxels of the plurality of image volumes with a corresponding radio frequency phase increment; constructing (14) voxel-wise sets of multi-dimensional bSSFP profiles for the sets of voxels from the sets of multi-dimensional bSSFP signal values; applying (15) voxel-wise a mathematical transformation on the multi-dimensional bSSFP profiles to obtain sets of multi-dimensional characteristic coefficients or sets of approximated multi-dimensional characteristic coefficients for the sets of voxels; using (16) voxel-wise the sets of multi-dimensional characteristic coefficients or approximated multi-dimensional characteristic coefficients to determine sets of parameter bases for the sets of voxels; and applying (17) one or more analytical solution functions onto the sets of parameter bases to quantify voxel-wise the one or more system parameters.

2. The method according to claim 1 , wherein the mathematical transformation is a Fourier transform or a Laplace transform, the multi-dimensional bSSFP signal values are complex bSSFP signal values, the multi-dimensional bSSFP profiles are complex-valued bSSFP profiles and the multi-dimensional characteristic coefficients are complex-valued Fourier modes.

3. The method according to claim 1 or 2, wherein multi-variant polynomials based on monomials of integer power are used to determine the sets of parameter bases.

4. The method according to claim 3, wherein the number of monomials is at most four for the single-compartment system and at most eight for the multi-compartment system for a respective parameter basis.

5. The method according to any one of the preceding claims, wherein the parameter bases comprise complex numbers.

6. The method according to any one of the preceding claims, wherein the one or more analytical solution functions are formulated by using at least one of the following elementary mathematical operations: equalling, addition, subtraction, multiplication, division, square-root operation, exponential function, logarithmic function, trigonometric function, atan2 function, complex conjugation and absolute value operation of one or more complex-valued numbers.

7. The method according to claim 6, wherein the operations are performed more than once for a respective analytical solution function, and/or wherein the operations use one or more parameters from a respective parameter basis and one or more reference values.

8. The method according to any one of the preceding claims, wherein the one or more analytical solution functions use one or more reference values, and wherein the one or more reference values are at least one of the following: one or more repetition times, one or more radio frequency excitation angles, a gyromagnetic ratio, a magnetic field strength, one or more reference frequencies, one or more chemical shifts, one or more nuclear magnetic resonance spectra, one or more magnetic field gradient strengths, one or more magnetic field gradient timings, and one or more magnetic field gradient lengths.

9. The method according to any one of the preceding claims, wherein a respective multidimensional bSSFP profile for a respective voxel is defined as follows:

wherein A_c, D_c, z_c and its complex conjugate z_c are two-dimensional numbers, A_c, D_c, z_c collectively form a respective parameter basis, c denotes a compartment index such that N_c = 1 indicates a single-compartment system and N_c > 1 indicates a multi-compartment system, M₊ is a two-dimensional number for a respective radio frequency phase cycle increment p_k, and wherein the sign in the exponents depends on the handedness of a reference coordinate system for two-dimensional numbers.

10. The method according to any one of the preceding claims, wherein the one or more multi-dimensional characteristic coefficients are Fourier modes given as:

wherein c_m denotes a Fourier mode, m denotes a mode number and N_PC denotes the number of images volumes, wherein modes with zero and/or positive mode number are obtained as follows:

while complex-conjugated modes with ()* being the complex conjugation operator with a negative mode number are obtained as follows:

11 . The method according to any one of the preceding claims, wherein a system of equations linking the multi-dimensional characteristic coefficients to the parameter bases are solved by using a Buchenberger algorithm.

12. The method according to any one of the preceding claims, wherein the method further comprises generating (18) a set of system parameter maps from the quantified system parameters, wherein a respective system parameter map is an array of the quantified system parameters across voxels of a respective image volume.

13. The method according to any one of the preceding claims, wherein the one or more system parameters are one or more of the following: longitudinal relaxation time 1 , transverse relaxation time T₂, effective transverse relaxation time T₂*> proton density M_o, magnetic field inhomogeneity B_o, proton density fraction, diffusion and excitation field inhomogeneity

14. A computer program product comprising instructions for implementing the following steps when loaded and run on a computing apparatus: receiving data relating to a balanced steady-state free precession (bSSFP) data acquisition of an object of the system to obtain a plurality of image volumes, each volume comprising a set of voxels, a respective image volume corresponding to an effective radio frequency phase increment; obtaining voxel-wise sets of multi-dimensional bSSFP signal values for the sets of voxels of the plurality of image volumes with a corresponding radio frequency phase increment; constructing voxel-wise sets of multi-dimensional bSSFP profiles for the sets of voxels from the sets of multi-dimensional bSSFP signal values; applying voxel-wise a mathematical transformation on the multi-dimensional bSSFP profiles to obtain sets of multi-dimensional characteristic coefficients or sets of approximated multi-dimensional characteristic coefficients for the sets of voxels; using voxel-wise the sets of multi-dimensional characteristic coefficients or approximated multi-dimensional characteristic coefficients to determine sets of parameter bases for the sets of voxels; and applying one or more analytical solution functions onto the sets of parameter bases to quantify voxel-wise the one or more system parameters.

15. An imaging apparatus (1 ) for carrying out a magnetic resonance imaging method for quantification of one or more system parameters of a single-compartment system and/or a multi-compartment system, the imaging apparatus (1 ) comprising means (7) for: receiving data relating to a balanced steady-state free precession (bSSFP) data acquisition of an object of the system to obtain a plurality of image volumes, each volume comprising a set of voxels, a respective image volume corresponding to an effective radio frequency phase increment; obtaining voxel-wise sets of multi-dimensional bSSFP signal values for the sets of voxels of the plurality of image volumes with a corresponding radio frequency phase increment; constructing voxel-wise sets of multi-dimensional bSSFP profiles for the sets of voxels from the sets of multi-dimensional bSSFP signal values; applying voxel-wise a mathematical transformation on the multi-dimensional bSSFP profiles to obtain sets of multi-dimensional characteristic coefficients or sets of approximated multi-dimensional characteristic coefficients for the sets of voxels; using voxel-wise the sets of multi-dimensional characteristic coefficients or approximated multi-dimensional characteristic coefficients to determine sets of parameter bases for the sets of voxels; and applying one or more analytical solution functions onto the sets of parameter bases to quantify voxel-wise the one or more system parameters.