METHOD DIRECTED TO MAGNETIC RESONANCE (MR) IMAGING SIMULATION

Abstract

The present invention describes method directed to magnetic resonance (MR) imaging simulation, said method comprising—partitioning 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 method directed to magnetic resonance (MR) imaging simulation, said method comprising partitioning 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.

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=t−t0: 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 3×3 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*(1−exp(−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(i−1), 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*(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. 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)
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(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).

Description

DESCRIPTION OF THE DRAWINGS

[0103] 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.

[0104] In FIG. 1 there is shown a general relationship for the signal amplitude against time and a though readout point.

[0105] For a traditional Bloch simulation, for each time point t, [0106] with dt=t(n)−t(n−1), then the following is performed:

Compute: E(dt), Eeq(dt), R(dt)

M(t+dt)=E(dt)R(dt)M(t)+Eeq(dt).

[0107] This known alternative is shown in FIG. 2.

[0108] 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=t−to, the following is performed:

[0109] Compute: E(h), Eeq(h), R(AREA)

M(t)=E(h)R(AREA)M(t0)+Eeq(h).

[0110] 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

[0111] [1] F. Bloch, “Nuclear induction”, Physical review, 70:460-474, 1946. [0112] [2] H. C. Torrey, “Bloch Equations with Diffusion Terms”, Physical review, 104:563-565, 1956. [0113] [3] A. Hazra, “Numerical Simulation of Bloch Equations for Dynamic Magnetic Resonance Imaging”, Ph.D. Dissertation, Georg-August-Universität Göttingen, 2016. [0114] [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. [0115] [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. [0116] [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.