Method directed to magnetic resonance (MR) imaging simulation

Abstract

The present invention describes method directed to magnetic resonance (MR) imaging simulation, said method comprisingpartitioning a pulse sequence in parts that have RF excitation and parts that do not have RF excitation; for each of the parts that do not have the RF active, called gradient parts, performing compression of the gradient parts and then representing signals driving the gradients as an accumulation of area until a readout time point; and then performing the simulation.

Claims

1. A computer-implemented method comprising magnetic resonance (MR) imaging simulation, said method comprising: partitioning, by at least one computer, a pulse sequence in parts that have RF excitation and parts that do not have RF excitation; for each of the parts that do not have the RF active, called gradient parts, using the at least one computer for performing compression of the gradient parts and then representing signals driving the gradients as an accumulation of area until a readout time point; resolving, by the at least one computer, the representing signals into an image using one or more simulations; and presenting, by the at least one computer, the image on a user interface.

2. The method according to claim 1, wherein the compression of the gradient parts is performed either by a pre-computation of (A, h) pairs where A are signals driving the gradients as an accumulation area until a readout time point, t(readout), and where h is h=t(readout)t0, with t being time, or on the fly by accruing phase and time.

3. The method according to claim 1, wherein the method is performed by performing the following: For each readout at time t, with h=tt0:
Compute:E(h),R(PHI),Eeq(h)
M.sub.READOUT(t)=E(h)R(PHI)M(t0)+Eeq(h) where M is magnetization vector, t is time, E and Eeq are the effects of relaxation constants T1 and T2 during the time t, and R(PHI) is the rotation due to accumulation of the phase and wherein the following is applied: M=[Mx, My, Mz] where [ ] indicates transpose, and E(h) and R(PHI) are 33 matrices according to the following:
E(h)=diag(exp(h/T2), exp(h/T2), exp(h/T1)) and where R(PHI) has the rotations and Eeq(h) is a vector with Eeq(h)=[0, 0, M0z*(1exp(h/T1))] with M0z being the equilibrium magnetization.

4. The method according to claim 1, wherein signals A driving the gradients as an accumulation area until the readout time point are presented as the integral of the signal amplitude over time, preferably numerically approximated, such as by A=Sum_i(dt*amplitude (t_i)), with t_i=t(i), dt=t(i)t(i1), for i in the discrete interval [1, readout].

5. The method according to claim 4, wherein the product gamma*A*Gradient provides the phase accrued due to corresponding signals, and where the rotation due to Bz fields, Bz fields being magnetic fields in the z component which are due to gradient fields and other effects but not 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 field non-linearities, Gx(r) Gy(r) and Gz(r) are the spatial distribution of the z component of the magnetic field (Bz) due to gradient coils, and Ax, Ay and Az are the AREAS computed by time integration for each of the gradient signals gx(t), gy(t) and gz(t).

6. The method according to claim 5, wherein the simulation for each readout time point is performed in parallel once the (A, h) pairs have been pre-computed.

7. The method according to claim 1, said method also comprising representing transverse magnetization as a time dependent phasor and thus representing a system in polar coordinates.

8. The method according to claim 7, wherein the method involves using the following relationship: if Mxy=Mm exp(i ) where Mm is magnitude and is phase, and i is the complex unit, then
Mx=real(Mxy)=Mm cos()
My=imag(Mxy)=Mm sin()
Mz=Mz.

9. The method according to claim 7, wherein computation of the simulation of the interval dt is 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)).

10. The method according to claim 8, wherein the method involves avoiding computing the exponentials by accumulating the time if a time point is not a time point of interest, preferably by use of the following requirements: 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.

11. The method according to claim 7, wherein the method involves tracking of a spatial derivate of the phase with respect to the spatial direction, preferably as extra variables according to the following:
d/dx=gamma*dt*(dBz(r)/dx+gx(t=ti)*dGx(r)/dx)
d/dy=gamma*dt*(dBz(r)/dy+gy(t=ti)*dGy(r)/dy)
d/dz=gamma*dt*(dBz(r)/dz+gz(t=ti)*dGz(r)/dz).

12. The method according to claim 7, wherein the method involves a time dependence of the position of the voxel, ri=r(t=ti), in order to track voxel 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(r))
Mz=M0z+(Mz(t0)M0z)*exp(h/T1(r)) 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*(dBz(ri)/dx+gx(t=ti)*dGx(ri)/dx)
d/dy=gamma*dt*(dBz(ri)/dy+gy(t=ti)*dGy(ri)/dy)
d/dz=gamma*dt*(dBz(ri)/dz+gz(t=ti)*dGz(ri)/dz).

Description

DESCRIPTION OF THE DRAWINGS

(1) Below there is described a representation of one embodiment according to the present invention, as compared to state of the art today, and with reference to the attached figures.

(2) In FIG. 1 there is shown a general relationship for the signal amplitude against time and a though readout point.

(3) For a traditional Bloch simulation, for each time point t, with dt=t(n)t(n1), then the following is performed:
Compute: E(dt),Eeq(dt),R(dt)
M(t+dt)=E(dt)R(dt)M(t)+Eeq(dt).

(4) This known alternative is shown in FIG. 2.

(5) In FIG. 3 there is shown the implication one embodiment according to the present invention, and with the proposed Bloch simulation. According to the present invention, for each Readout at time t, with h=tto, the following is performed:
Compute: E(h),Eeq(h),R(AREA)
M(t)=E(h)R(AREA)M(t0)+Eeq(h).

(6) In FIG. 4 there is shown one possible continuation of the embodiment according to the present invention, as shown in FIG. 3. In FIG. 4 there is shown the illustration of a Riemann integral. Here it should be noted that the area is computed as the summation of individual areas.

REFERENCES

(7) [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-Universitt Gttingen, 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.