-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathDictionary_function_CRLB.m
141 lines (110 loc) · 4.09 KB
/
Dictionary_function_CRLB.m
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
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
%%% Simulates the steady-state signal for a bSSFP sequence with multiple
%%% single-band (1B) and multi-band (MB) pulses.
function Mss_Total = Dictionary_function_CRLB(Theta,flips,TR,Dur,dphi,Delta,TBW,nMB,n1B,b1sqrd,G_1B,G_2B,G_3B,wloc_1B,wloc_2B,wloc_3B)
T1x = [Theta(1),Theta(2),Theta(3)*1e-3];
T2x = [Theta(4),Theta(5)]*1e-3;
f = Theta(9);
M0b = Theta(6); M0f = (1-M0b); M0 = Theta(7)*[M0f M0b*(1-f) M0b*f];
K = Theta(8);
%% Parameters for signal calculation.
B1ms_1B = [0 b1sqrd{1}(1,2) 0]; B1ms_2B = b1sqrd{1}; B1ms_3B = b1sqrd{2};
df_1B = [0 0 0]; df_2B = [0 0 Delta]; df_3B = [-Delta 0 Delta];
Gamma = 267.52219; nd = length(B1ms_1B);
W_3B = pi*Gamma^2*B1ms_3B.*G_3B;
W_2B = pi*Gamma^2*B1ms_2B.*G_2B;
W_1B = pi*Gamma^2*B1ms_1B.*G_1B;
%% Evolution during RF - flip and saturate. Two separate terms needed for 1B and MB.
WPulse_1B = zeros(3,3); WPulse_2B = zeros(3,3); WPulse_3B = zeros(3,3);
for ii=1:nd
D = 2*pi*df_1B(ii)/wloc_1B;
WPulse_1B = WPulse_1B + W_1B(ii)*([[-1 0 0];[0 -1 D];[0 D -D^2]]);
end
for ii=1:nd
D = 2*pi*df_2B(ii)/wloc_2B;
WPulse_2B = WPulse_2B + W_2B(ii)*([[-1 0 0];[0 -1 D];[0 D -D^2]]);
end
for ii=1:nd
D = 2*pi*df_3B(ii)/wloc_3B;
WPulse_3B = WPulse_3B + W_3B(ii)*([[-1 0 0];[0 -1 D];[0 D -D^2]]);
end
% Assume RF is rotation around x-axis. Phase doesn't matter here.
% Pick the sum(alpha) assuming central is the only non-zero frequency band.
T = RF_Rot(sum(flips),0);
AMat_1B = blkdiag(T,expm(WPulse_1B*Dur));
AMat_2B = blkdiag(T,expm(WPulse_2B*Dur));
AMat_3B = blkdiag(T,expm(WPulse_3B*Dur));
%% Evolution during free precession.
% Modify T2 as per Bieri et al.
R1a = 1/T1x(1); R2a = 1/T2x(1);
T_RFE = 1.20*(Dur/TBW);
Eta = 0.68-0.125*(1+(T_RFE/TR))*(R1a/R2a);
R2a_mod = (1-Eta*(T_RFE/TR))*R2a;
% R2a_mod = R2a;
M0f = M0(1);
M0s = M0(2)+M0(3); % Total semisolid fraction.
f = M0(3)/M0s; % Fraction of semisolid pool that is dipolar coupled.
% Reset f value from NaN.
if M0s == 0; f = 0; end
T1 = [T1x(1:2) T1x(2) T1x(3)]; % Replicate semisolid Zeeman T1.
R1 = 1./T1;
% Relaxation/exchange matrices.
E = diag([-R2a_mod -R2a_mod -R1]);
KK = blkdiag(0,0,[[-K*M0s K*M0f K*M0f];[K*(1-f)*M0s -K*M0f 0];[K*f*M0s 0 -K*M0f]],0);
M0tot = M0f+M0s;
KK = KK/M0tot; %%<--- Shaihan: these should not have units of magnetization.
O = diag([-1i*dphi/TR 1i*dphi/TR 0 0 0 0]);
Lambda = E+KK+O; % Matrix exponential of this will form Smat.
C = [0 0 R1(1)*M0f R1(2)*(1-f)*M0s R1(3)*f*M0s 0].';
Dphase = diag([-1 -1 1 1 1 1]); % Phase alternation matrix. %Dphase = diag([0 0 1 1 1 1]); % Spoiling matrix.
%% Perform eigenvector decomposition. Use augmented matrix approach.
LC = cat(1,[Lambda C],zeros(1,7));
AlphaPrime_1B = blkdiag(AMat_1B,1);
AlphaPrime_2B = blkdiag(AMat_2B,1);
AlphaPrime_3B = blkdiag(AMat_3B,1); % Post-expm, so 1 is added rather than 0.
Dphase = blkdiag(Dphase,1);
RDS_1B = AlphaPrime_1B * Dphase * expm(LC*TR);
RDS_2B = AlphaPrime_2B * Dphase * expm(LC*TR);
RDS_3B = AlphaPrime_3B * Dphase * expm(LC*TR);
%%% Signal calculation. %%%
% First reference point.
X = RDS_2B * (RDS_1B)^n1B * (RDS_3B)^nMB * (RDS_1B)^n1B * (RDS_2B)^(nMB-1);
[V,~,~] = eig(X);
% Create array to hold all state vectors.
Mstates = zeros([size(X,1) 2*(n1B+nMB)]);
Mstates(:,1) = V(1:end,end)/V(end,end);
Counter = 2;
% 2-band period.
for ii = 2:nMB
Mstates(:,Counter) = RDS_2B * Mstates(:,Counter-1);
Counter=Counter+1;
end
% 1-band period.
for ii = 1:n1B
Mstates(:,Counter) = RDS_1B * Mstates(:,Counter-1);
Counter=Counter+1;
end
% 3-band period.
for ii = 1:nMB
Mstates(:,Counter) = RDS_3B * Mstates(:,Counter-1);
Counter=Counter+1;
end
% 1-band period.
for ii = 1:n1B
Mstates(:,Counter) = RDS_1B * Mstates(:,Counter-1);
Counter=Counter+1;
end
Mss_Total = Mstates(1,:);
%%% Deals with rotation. %%%
function Tap = RF_Rot(a,p)
Tap = zeros([3 3]);
Tap(1) = cos(a/2).^2;
Tap(2) = exp(-2*1i*p)*(sin(a/2)).^2;
Tap(3) = -0.5*1i*exp(-1i*p)*sin(a);
Tap(4) = conj(Tap(2));
Tap(5) = Tap(1);
Tap(6) = 0.5*1i*exp(1i*p)*sin(a);
Tap(7) = -1i*exp(1i*p)*sin(a);
Tap(8) = 1i*exp(-1i*p)*sin(a);
Tap(9) = cos(a);
end
end