MRI Simulation - GRE and SSFP
这里记录一下不同种类的 Sequence 得到的 M。
比较 GRE 与 refocussed-SSFP 的效果
主程序
首先是主程序：1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35T1 = 600; % ms.
T2 = 100; % ms.
TE = 0:2.5:10; % ms.
TR = 10; % ms.
flip = pi/3;
df = [-100:100]; %Hz
Sig = zeros(length(df),length(TE));
for n=1:length(TE)
for k=1:length(df)
[Msig,Mss] = sssignal(flip,T1,T2,TE(n),TR,df(k));
% [Mss,Msig] = gresignal(flip,T1,T2,TE(n),TR,df(k));
Sig(k,n)=Mss(1)+i*Mss(2);
end;
end;
% ===== Plot the Results ======
subplot(2,1,1);
plot(df,abs(Sig));
xlabel('Frequency (Hz)');
ylabel('Magnitude');
grid on;
legend('TE=0', 'TE=2.5', 'TE=5.0', 'TE=7.5', 'TE=10');
subplot(2,1,2);
plot(df,angle(Sig));
xlabel('Frequency (Hz)');
ylabel('Phase (radians)');
axis([min(df) max(df) -pi pi]);
grid on;
legend('TE=0', 'TE=2.5', 'TE=5.0', 'TE=7.5', 'TE=10');
比较结果、文档解释
上面这个是基于 Steady state signal 得到的情况，图片为：
这一情况被称为 balanced steady-state free-precession (SSFP), 这里sssignal.m据文档所述，已经计算好了 refocussed-SSFP。
The balanced SSFP sequence simply consists of alpha excitation pulses spaced TR apart. TR is usually very short, on the order of several milliseconds.
Refocussed-SSFP means that the imaging gradients are fully refocussed. But we haven’t put in any imaging gradients yet, so this is true!
将上面的sssignal.m函数直接替换成gresignal.m，得到下面的这张图：
在文档中是这样说的：
Notice that the signal level of gradient echo (GRE) signal is exactly the same as the mean refocussed-SSFP signal. That makes sense, because the GRE signal was just the average over many spins that had different amounts of phase twist.
说明一件事情：GRE signal 只是 SSFP 的升级版，区别在于 GRE 考虑很多个 spin 并取均值。每个 spin 不同就不同在他们的 dephase 相位是不同的。具体看下它们各自的代码。
具体程序比较
首先看sssignal.m：1
2
3
4
5
6
7function [Msig,Mss] = sssignal(flip,T1,T2,TE,TR,dfreq)
Rflip = yrot(flip);
[Atr,Btr] = freeprecess(TR-TE,T1,T2,dfreq);
[Ate,Bte] = freeprecess(TE,T1,T2,dfreq);
Mss = inv(eye(3)-Ate*Rflip*Atr) * (Ate*Rflip*Btr+Bte);
Msig = Mss(1)+1i*Mss(2);
这个就是很平淡无奇地计算 steady state. 但是下面的 GRE 考虑的东西很多：1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28function [Mss,Msig]=gresignal(flip,T1,T2,TE,TR,dfreq)
N = 100;
M = zeros(3,N);
phi = ([1:N]/N - 0.5) * 4 * pi % 先在 [-2pi, 2pi]之间获得均匀分布的 N 个phi。
for k=1:100
M(:,k) = gssignal(flip,T1,T2,TE,TR,dfreq,phi(k));
end
Mss = mean(M, 2);
Msig = Mss(1)+i*Mss(2);
% 把两个函数放在同一个框内了，在实际使用的时候不要放在同一个.m文件中
function [Mss,Msig]=gssignal(flip,T1,T2,TE,TR,dfreq,phi)
Rflip = yrot(flip);
[Atr,Btr] = freeprecess(TR-TE,T1,T2,dfreq);
[Ate,Bte] = freeprecess(TE,T1,T2,dfreq);
Atr = zrot(phi)*Atr; % TR末尾的dephase
% 下面的和sssignal一样的
Mss = inv(eye(3)-Ate*Rflip*Atr) * (Ate*Rflip*Btr+Bte);
Msig = Mss(1)+i*Mss(2);
% phi is the angle by which the magnetization is dephased at the end of the TR
% >>> Mss=gssignal(pi/3,600,100,2,10,0,pi/2)
% >>> Mss=[0.1248, 0.1129, 0.1965]'
% 如果加入了一个线性的Gz，那么跟随不同的z坐标，不同的M就会有不同的dephase
% 所以这里只是模拟了一下，假设给了一个dephase角度，该怎么算。
% 这里给的结论就是按 Rz(dephase)来算的。
在文档中的解释：
gresignal.m calculates the average magnetization over, say, 100 points with the phase varying uniformly from -2pi to 2pi.
也就是说这里区别在于多了一个 phi。这个就是模拟梯度磁场进来之后，对其中100个点（100个不同z坐标上的点）的相位影响，具体影响就在 Atr = zrot(phi)*Atr; 也就是按z轴上旋转了phi的角度。
这样一来，前后就都能串起来了。