Methods for targeting magnetic resonance (MR) imaging simulations

Abstract

A plurality of transmitters arranged in a predetermined space and transmitting transmission signals for wirelessly transmitting power 1 The present invention describes a method for magnetic resonance (MR) imaging simulation, which includes dividing a pulse sequence into portions with RF excitation and portions without RF excitation, and for each portion without RF active, called a gradient portion, performing a compression of the gradient portion, then expressing the signal driving the gradient as an accumulation of area up to the readout time, and then performing the simulation.

Description

[Technical Field]

[0001] The present invention relates to a method directed to magnetic resonance (MR) imaging simulation. [Background technology]

[0002] The present invention is generally directed to advanced MRI simulations that do not require the use of an actual MRI scanner. More particularly, the present invention is directed to a method for magnetic resonance (MR) imaging simulation that can advantageously analyze a particular hypothetical volume of an object using a preferred computer model. Summary of the Invention [Means for solving the problem]

[0003] The latter object is achieved by a method directed to magnetic resonance (MR) imaging simulation, said method comprising: - dividing the pulse sequence into a part with RF excitation and a part without RF excitation; - for each part that does not have RF activity, called gradient part, carry out a compression of the gradient part, and then express the signal that drives the gradient as an accumulation of the area up to the readout time; -Then run the simulation Includes. The present invention is directed to a method that focuses on the portions that do not have RF excitation, which further means that the method includes isolating these non-RF portions, which becomes possible when RF transmission is turned off. Furthermore, as can be seen from the above, the method according to the invention involves compressing the gradient portion and then expressing the signal driving the gradient as an accumulation of area up to the readout time, which means that the time compression of the sequence is defined not as a time signal but as an accumulation of area up to the readout time (see (A,h) pair below), which is a novel feature and concept compared to known methods currently in use. [Brief explanation of the drawings]

[0004] [Figure 1] The general relationship of signal amplitude versus time and readout points is shown. [Figure 2] This known alternative is shown below. [Figure 3] The implications of one embodiment of the present invention are illustrated using the proposed Bloch simulation. [Figure 4] As shown in FIG. 3, one possible continuation of an embodiment according to the invention is shown. DETAILED DESCRIPTION OF THE INVENTION

[0005] In the following, a representation of an embodiment of the present invention compared with the state of the art will be described with reference to the accompanying drawings.

[0006] Figure 1 shows the general relationship of signal amplitude and readout points with respect to time. In a traditional Bloch simulation, for each time t, let dt = t(n) - t(n-1), then:

number

[0007] Figure 2 illustrates this known alternative.

[0008] Figure 3 shows the implications of one embodiment according to the present invention using the proposed Bloch simulation. According to the present invention, for each readout at time t, if h = t - t0, the following is executed.

Number

[0009] Figure 4 shows one possible continuation of an embodiment according to the present invention as shown in Figure 3. Figure 4 shows a diagram of Riemann integration. Here, it should be noted that the area is calculated as the sum of the individual areas.

[0010] <Description of the Present Invention and Related Technical Background> <Area - Based Sequence Compression for High - Speed MRI Simulation> <Basis of MRI Simulation> MRI technology is based on the visual interpretation of the time evolution of the magnetization of specific atomic protons (spins) in soft tissue when excited by a smaller fluctuating magnetic field in the presence of a strong, constant magnetic field (main magnetic field).

[0011] Therefore, MRI is based on the implicit use of spin dynamics, and the development of magnetization in the presence of an external fluctuating magnetic field is modeled at a macroscopic scale by the Bloch equations [1]. The simulation of MRI phenomena depends on solving the Bloch equations (or upgraded derivative versions incorporating other physics such as Bloch - Torrey [2] for diffusion).

[0012] The basic form of the Bloch equation for MRI purposes is as follows.

Number

[0013] <Fundamentals of MRI> MRI is based on a series of events.

[0014] equilibrium : Usually, a strong constant main magnetic field B0 is applied in the z direction to align the spins in the initial state M0z (equilibrium) and cause precession motion according to the magnetic rotation ratio. Regarding z above and below, the z direction is considered the direction of the main magnetic field.

[0015] RF excitation : A short radio frequency (RF) magnetic field with transverse (x and y) magnetic components (Bx and By) at the frequency of the precession motion is applied to perturb the magnetization of the spins and make it a specific combination of Mx, My, and Mz.

[0016] Encoding and Retrieval : After the RF is turned off, a z-component gradient magnetic field (Bz that varies with position) is applied to generate so-called "spatial encoding." As a result, each voxel has a precession motion related to its position, causing a difference in the magnetization of Mx and My (transverse magnetization), which is later used to generate an image.

[0017] Due to the combined effect of these transverse magnetizations, an RF magnetic field with x and y components is radiated, which is received by a coil using Faraday's law and converted into an electrical signal S(t) = S(Mx(t), My(t)). Here, S is a function of the receiving coil sensitivity.

[0018] relief : Since no (which may be complete or partial) external magnetic field is applied, the transverse (xy) magnetization decays naturally by T2, while the Mz magnetization recovers towards M0z with a constant T1.

[0019] MRI acquires the time signal S(t) using precise timing and a sequence of repetitions of these basic events. These signals are decomposed into images using knowledge of encoding, and the type of signal applied depends more or less on T1 and T2, resulting in a natural contrast based on tissue characteristics.

[0020] <Details of MRI Simulation> Assuming that MRI operates under a constant strong static magnetic field B0 in the z direction and Bz(t) = B0 + ΔBz(t) (where B0 is the static magnetic field value and ΔBz(t) is the variation due to external excitation or inhomogeneity and non - linearity), a typical representation for the simulation of the Bloch equations is based on the rotating coordinate system of reference [ref4], and the precession effect of B0 is removed from Bz to simplify the formulation .

[0021] The Bx and By magnetic fields are due to RF excitation, while the Bz magnetic field is due to other effects such as the effect of the gradient magnetic field and inhomogeneities due to physical effects (magnetic susceptibility, chemical shift) or hardware imperfections. Therefore, when the effect of the B0 magnetic field is removed, the z - component of the magnetic field is as follows.

Equation

[0022] In this regard, it should also be noted that in more advanced cases the gradient maps may also be time-varying as follows: Gx(r,t), Gy(r,t), Gz(r,t).

[0023] The Bloch equations are ordinary differential equations (ODEs) [3] and can be expressed as dM / dt = F(M, B, T1, T2, M0z), where F is a function and the time and position dependence is implicit.

[0024] The ODE can be solved numerically by time discretization, where time is divided into small segments where the B field is assumed to be constant, given initial conditions, and the solution for the magnetization at each time step can be calculated iteratively by time integration from the previous solution.

[0025] There are other well-known methods that can be applied [3], for example, to solve the differential equation using an explicit forward difference scheme where dt = t(n) - t(n-1) and M(t) at t(n-1) is known:

number

[0026] Alternatively, an analytical solution can be obtained for a constant magnetic field time step by applying a matrix exponential. Given the structure of the Bloch equations, this operation can be decomposed as a series of rotations (dependent on the B field) and exponential decays [4, 5], which can be expressed in compact form as

number

number

number

[0027] In Cartesian components, the effect of the exponential decay operations E(dt) and Eeq(dt) during the interval dt is:

number

[0028] The computation of R(dt) involves a series of rotations, the angle of which depends on the B-field during the interval dt. A typical representation [4] decomposes R(dt) into a series of rotations:

number

number

number

number

[0029] Other representations and approximations of the same operations are possible [6], and these may be computationally advantageous in certain situations.

[0030] In either case, numerical Bloch simulations rely on time discretization into small intervals with a constant B field, and numerical or analytical time integration based on the magnetization at the beginning of the interval. Note that this is an iterative procedure, where each new solution requires the calculation of the previous solution.

[0031] When the magnetic field varies, the accuracy of the approximation depends on the time step, with shorter intervals resulting in smaller errors. However, each step must be applied to potentially millions of voxels, making the simulation expensive. In situations where the magnetic field remains constant for long periods of time, a solution can be obtained in a single (larger) time step without loss of accuracy [4].

[0032] Note that in most MRI simulations, the only time points of interest where magnetization is needed are during the acquisition steps (called READOUTS) or at the beginning and end of the sequence. These are only a small fraction of the time points that actually need to be calculated for an accurate simulation.

[0033] For practical purposes, such as simulating a three-dimensional field of view (FOV) containing biological tissue, spatial discretization is also applied, dividing the three-dimensional FOV into small voxels. A voxel is a small volume containing multiple nuclei with the same tissue properties and is assumed to experience the same external magnetic field due to their proximity. Therefore, B, M, T1, T2, and M0z are assumed to be constant for each voxel to which a coordinate r and volume are assigned. The magnetization of this voxel can be simulated using the procedure described above. In relation to the above, it should be noted that the present invention also covers the case where, in practice, there is more than one tissue type within a single voxel.

[0034] A typical Bloch simulation for MRI is based on the numerical approach presented above, where time is discretized finely, discretizing the time domain into small steps.

[0035] <Embodiments of the present invention and related details> The present invention relates to a method involving area-based sequence compression. According to the present invention, we propose applying an alternative representation of the sequence. As a preprocessing step, the sequence is divided into parts, separating those with RF excitation, and thus time-varying Bx, By, and Bz magnetic fields, from those without RF excitation (the encoding, acquisition, and relaxation parts of the sequence), where only the time-varying Bz magnetic field is present, i.e., Bx(t) = By(t) = 0. This can be done by identifying and grouping time steps with non-zero signals driving the transmit coils. If the sequence is not yet represented using time discretization, time discretization can be applied in advance.

[0036] For each section without RF activity, preferably called the GRADIENT section, the signal driving the gradient can be expressed as the accumulation of area A up to the READOUT time, where A=F(amplitude, h, t) (F is a function), and h=t(readout)-t0, where t0 is the initial time for this section, and only these value pairs (A, h) at the time of interest (READOUTS) are saved. The area A can be calculated (analytically or numerically using Riemann integrals) by integrating the signal over time.

[0037] Based on the above, according to one specific embodiment of the present invention, hereinafter referred to as the first embodiment aspect, compression of the gradient part is performed on the fly by pre-computing the (A,h) pair, where A is the signal driving the gradient as accumulated area up to the readout time t(readout) and h is h=t(readout)-t0, where t is time, or by generating the phase and time.

[0038] In line with the above, according to yet another specific embodiment of the present invention, the signal A driving the gradient as cumulative area up to the readout point is expressed as the integral of the signal amplitude over time, preferably for i within the discrete interval [1, readout]:

number

[0039] Expressing the sequences in each GRADIENT section as a set of AREAS(A) and PERIODS(h), a simulation can be performed. Note that the product gamma*A*Gradient gives the phase produced by the corresponding signal, and therefore the rotation due to the Bz field corresponds to:

number

[0040] Furthermore, in this case, given the initial magnetization (M0) at the start of the GRADIENT portion (t0), it is possible to calculate the magnetization solution at each of the READOUT points using a single step consisting of rotation and exponential decay applied to the initial magnetization M0 at the start of the GRADIENT sequence portion.

[0041] Based on the above, according to one specific embodiment of the present invention, the method is carried out by: for each reading at time t (h=t-t0), calculating the attenuation and rotation of: E(h), R(PHI), Eeq(h)

[0042] Applying these to M(t0), we obtain the solution at read time t as follows:

number

[0043] In this case, we only have the z component of the magnetic field, so the rotation is about the z axis:

number

[0044] Therefore, MREADOUT(t)=E(h)R(PHI)M(t0)+Eeq(h) can be expressed in Cartesian components as follows:

number

[0045] Note that without compression, the solution must instead be calculated in an iterative manner.

number

[0046] Furthermore, note that the iterative procedure is equivalent to dividing the signal duration dt into piecewise constant steps and applying a rotation of the area associated with the small steps. R(dt) = (Rho) where:

number

[0047] In a sense, the traditional (iterative) approach is conceptually similar to applying a Riemann integral with individual rotations and decays at each step.

[0048] Furthermore, in the following, a series of implications and advantages of the concepts according to the embodiments of the first aspect of the present invention disclosed above are provided.

[0049] 1) For each point of interest (READOUTS and segment endpoints), a precomputation is required to reformulate the sequences gx(t), gy(t) and gz(t) as a set of area-duration pairs (A, h). This is very efficient and cheap, since it only needs to be done once per sequence, and only needs to be done once per sequence.

[0050] 2) The accuracy of area calculation can be controlled arbitrarily. It is accurate in analytical integration, and can improve the discretization of the time signal in numerical integration without affecting the performance of the actual simulation.

[0051] 3) Instead of having to apply the numerical simulation to many time steps using dt discretization, rotation and damping only need to be applied once for each point of interest (READOUTS). This is particularly relevant for realistic waveforms (trapezoidal) and more advanced encodings, such as spirals, which employ sinusoidal waveforms and require finer discretization for high accuracy. This significantly improves performance, especially since realistic simulations may require applying these operations to millions of voxels.

[0052] 4) The solution for each READOUT depends only on the initial solution M(t0) at the start of the GRADIENT portion and its (A, h). Therefore, the solution for each READOUT can be calculated in parallel using multiple computer machines or computer cores, which further increases the potential for performance improvement. Based on this, according to one specific embodiment, once the (A, h) pairs are calculated in advance, the simulations for each readout instant are performed in parallel.

[0053] 5) The RF section must still be simulated in the conventional way. Also, the sections must be simulated sequentially. The initial M(t0) of each GRADIENT section is only available if the previous section has been simulated.

[0054] 6) The effect of magnetic field inhomogeneity or other effects can be accounted for as a constant (or time-varying, if applicable) signal and can be included as gamma*h*ΔBz(r) as above for each READOUT.

[0055] According to an aspect of the second embodiment of the present invention, an alternative representation of the magnetization, and therefore of the Bloch equations, is applied. Instead of representing the system in Cartesian coordinates, it is proposed to use polar coordinates, which represent the transverse magnetization as a time-dependent phasor. Thus, according to one embodiment of the present invention, the method also includes representing the transverse magnetization as a time-dependent phasor, and therefore representing the system in polar coordinates.

[0056] Using the magnitude (Mm) and phase (φ), the following can be applied:

number

[0057] There is a clear relationship between the Cartesian and polar components:

number

[0058] There are several advantages to using this representation: As with aspects of the first embodiment described above, when simulating a time point with no active RF (Bx(t)=By(t)=0), the calculations for simulating step dt are reduced to:

number

[0059] Based on the above, according to yet another embodiment of the present invention, the method includes avoiding calculating an exponential function by accumulating time when a time point is not the time point of interest, preferably using the following requirement:

number

[0060] Therefore, according to yet another embodiment, the method comprises tracking the spatial derivative of the phase with respect to the spatial direction as an additional variable, preferably according to:

number

[0061] This is useful in several relevant cases, namely to capture phase variations over voxels, to remove spurious echoes (which requires phase integrals over voxels), or to calculate scattering effects, which requires their derivatives [6].

[0062] The following provides conceptual speculations and advantages according to aspects of this second embodiment of the invention disclosed above.

[0063] 1) This approach can be useful in cases where sequence compression according to the first embodiment described above is not applicable. For example, in models involving motion or flow [5], the voxel position changes over time. Instead of applying a more complex approach, the generated phase can be easily accumulated even if the voxel position r changes over time. This motion can be tracked very accurately with small time steps, generating many time points that are not points of interest. Whenever magnetization needs to be calculated (point of interest), the accumulated phase and time steps can be used to generate a signal.

[0064] 2) Performance Improvement: When the approach according to the first aspect is not applicable, the computation is significantly reduced since there is no need to compute exponential functions or rotations (which require evaluating sin and cos), which are much more expensive operations than simple multiplication and addition.

[0065] 3) By tracking the phase over time, the effects of phase changes on a voxel can be easily accounted for. In Cartesian coordinates, the actual phase is lost and can only be obtained in the [0,2π] or [-π,π] interval. This creates difficult-to-handle phase wrapping, which can limit accurate calculation of scattering and other phase effects.

[0066] It should be understood that, according to the present invention, both of the concepts and approaches described above can be applied independently or can be combined.

[0067] When combined, the phase accumulation is obtained directly from the pre-calculated area A for each signal.

number

[0068] In this case, all that is required is to calculate the READOUT time. For each READOUT,

number

[0069] <Reference materials> [1] F. Bloch, “Nuclear induction”, Physical review, 70:460-474, 1946. [2] H.C. Torrey, “Bloch Equations with Diffusion Terms”, Physical review, 104:563-565, 1956. [3] A. Hazra, “Numerical Simulation of Bloch Equations for Dynamic Magnetic Resonance Imaging”, Ph.D. Dissertation, Georg-August-Universitat Gottingen, 2016. [4] C.G. Xanthis, I.E. Venetis, A.V. Chalkias, A.H. Aletras, “MRISIMUL: a GPU-based parallel approach to MRI simulations”, IEEE Trans. Med. Imaging 33: 607-616, 2014. [5] C.G. Xanthis, I.E. Venetis, A.H. Aletras, “High performance MRI simulations of motion on multi-GPU systems”, J. Cardiovasc. Magn. Reson. 16:48-62, 2014. [6] T. H. Jochimsen, A. Schafer, R. Bammer, and M. E. Moseley, “Efficient simulation of magnetic resonance imaging with Bloch-Torrey equations using intra-voxel magnetization gradients,” J. Magn. Reson., 180:29-38, 2006.

Claims

1. 1. A method directed to magnetic resonance (MR) imaging simulation, comprising: Dividing a predetermined pulse sequence into an RF excitation portion in which an RF pulse is applied and a non-RF excitation portion including a gradient period in which a gradient magnetic field is applied without applying an RF pulse; determining, for each gradient period, the duration from the start (t0) of the gradient period to the end (readout) of the gradient period (h=t(readout)-t0 (t is time)) and a parameter (A) indicative of the cumulative effect of the amplitude of the gradient magnetic field applied during the gradient period over time; using the determined duration (h) and a parameter (A) indicative of the cumulative effect to calculate jointly the rotation and relaxation of the magnetization vector during the gradient period, thereby simplifying the calculation of the change in magnetization over time during the gradient period; calculating the change in magnetization over time for the RF excitation portion and using the simplified calculation for the gradient period of the non-RF excitation portion, and performing the magnetic resonance (MR) imaging simulation by combining these; A method comprising:

2. 2. The method of claim 1, wherein the determined duration (h) and a pair of parameters (A) indicative of a cumulative effect are pre-calculated before the magnetic resonance (MR) imaging simulation is performed, or are determined on the fly by sequentially generating phase changes and elapsed times corresponding to the gradient periods while the magnetic resonance (MR) imaging simulation is being performed.

3. 3. The method of claim 1 or 2, The calculation of the change in magnetization over time during the gradient period is performed by executing, for each readout at time t (h=t-t0), Compute: E(h), R(PHI), Eeq(h)MREADOUT (t) = E(h) R(PHI) M(t0) + Eeq(h); where M is the magnetization vector, t is time, E and E are the effects of the relaxation constants T and T during time t, and R(PHI) is the rotation of the phase by the total cumulative phase (PHI), M = [Mx, My, Mz]', where [ ]′ denotes transposition, and E(h) and R(PHI) are 3x3 matrices, E(h) = diag(exp(-h / T2), exp(-h / T2), exp(-h / T1)). where R(PHI) has the rotation and Eeq(h) is a vector Eeq(h) = [0,0,M0z*(1-exp(-h / T1))]', where M0z is the equilibrium magnetization, however the method is applied.

4. 4. A method according to any one of claims 1 to 3, wherein the signal A driving the ramp period as a cumulative area up to the readout time is expressed as the integral of the signal amplitude over time, dividing the ramp period into discrete intervals (dt) and numerically approximated by A = Sum_i(dt * amplitude(t_i)), where t_i = t(i), dt = t(i) - t(i-1).

5. 4. The method of claim 3, wherein the total accumulated phase (PHI) is a phase PHI(r) that depends on the position r, as follows: The product gamma*A*Gradient gives the phase caused by the corresponding signal, and the rotation due to the Bz field, which is the z-component field due to gradients and other effects but not the RF excitation, corresponds to PHI(r) = gamma*(h * ΔBz(r) + Ax * Gx(r) + Ay * Gy(r) + Az * Gz(r)), where r is the coordinate of the "spin" or proton, ΔBz(r) represents the spatial distribution of the magnetic field nonlinearity, Gx(r), Gy(r) and Gz(r) are the spatial distributions of the z-component of the magnetic field (Bz) due to the gradient coils, and Ax, Ay and Az are the AREAS calculated by the time integral of the gradient signals gx(t), gy(t) and gz(t), respectively.

6. 6. The method of claim 5, wherein the calculation of the magnetization at each readout time point, which is the end point of each gradient period, is performed in parallel for multiple readout time points, when the pair of the duration (h) and parameter (A) indicating the cumulative effect is pre-calculated before the magnetic resonance (MR) imaging simulation is performed.

7. 7. A method according to any one of claims 1 to 6, which also comprises expressing the transverse magnetization as a time-dependent phasor and therefore expressing the system in polar coordinates.

8. 8. The method of claim 7, comprising using the relationships: Mxy=Mm exp(iφ), where Mm is magnitude, φ is phase, and i is the complex unit; Mx=real(Mxy)=Mm*cos(φ), My=imag(Mxy)=Mm*sin(φ), Mz=Mz.

9. 9. The method of claim 8, wherein, for a simulation interval dt, the phase (φ), transverse magnetization magnitude (Mm), and longitudinal magnetization (Mz) updates are calculated based on φ(t+dt) = φ(t) - gamma*dt*(ΔBz(r)+gx(t=ti)*Gx(r)+gy(t=ti)*Gy(r)+gz(t=ti)*Gz(r))Mm(t+dt) = Mm(t) * exp(-dt / T2(r))Mz(t+dt) = M0z + (Mz(t) - M0z)*exp(-dt / T1(r)), where ti is the start time t (t=ti) for identifying the gradient field value for each time step dt.

10. 10. A method according to claim 8 or 9, comprising, if a time point is not of interest, calculating as follows: If READOUT: φ-= gamma*dt*(ΔBz(r) + gx(t=ti) * Gx(r) + gy(t=ti) * Gy(r) + gz(t=ti) * Gz(r)) Mm = Mm(t0) * exp(-h / T2(r)) Mz = M0z + (Mz(t0) - M0z )* exp(-h / T1(r)), else: φ-= gamma*dt*(ΔBz(r) + gx(t=ti) * Gx(r) + gy(t=ti) * Gy(r) + gz(t=ti) * Gz(r)) h += dt requirement to avoid calculating the exponential function by accumulating time.

11. 11. A method according to any one of claims 7 to 10, comprising tracking spatial derivatives of the phase with respect to spatial directions as additional variables by dφ / dx-= gamma*dt*(dΔBz(r) / dx + gx(t=ti) * dGx(r) / dx)dφ / dy-= gamma*dt*(dΔBz(r) / dy + gy(t=ti) * dGy(r) / dy)dφ / dz-= gamma*dt*(dΔBz(r) / dz + gz(t=ti) * dGz(r) / dz).

12. 12. A method according to any one of claims 7 to 11, wherein the time dependence of a voxel position ri=r(t=ti) of a voxel is calculated for tracking its movement. If READOUT:φ-= gamma*dt*(ΔBz(ri) + gx(t=ti) * Gx(ri) + gy(t=ti) * Gy(ri) + gz(t=ti) * Gz(ri)) Mm = Mm(t0) *exp(-h / T2) Mz = M0z + (Mz(t0) - M0z )*exp(-h / T1) else:φ-= gamma*dt*(ΔBz(ri) + gx(t=ti) * Gx(ri) + gy(t=ti) * Gy(ri) + gz(t=ti) * Gz(ri)) h += dt and dφ / dx-= gamma*dt*(dΔBz(ri) / dx + A method involving gx(t=ti) * dGx(ri) / dx)dφ / dy-= gamma*dt*(dΔBz(ri) / dy + gy(t=ti) * dGy(ri) / dy)dφ / dz-= gamma*dt*(dΔBz(ri) / dz + gz(t=ti) * dGz(ri) / dz).

Images