-
Notifications
You must be signed in to change notification settings - Fork 1
/
optim_minBias.m
84 lines (63 loc) · 2.31 KB
/
optim_minBias.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
% Minimum Bias cost function
%
% Created by Arian Beqiri, King's College London, September 2016.
% Email: [email protected]
%
% This code is free under the terms of the MIT license.
function [fval] = optim_minBias(P_A,Ap,Amask,VOP,QG,targmask,idx,...
b1_act_scale,Nc,t_enc,lSARmax,wbSARmax,delta_1,delta_2,theta,...
A_amp,P_av,P_peak)
% Load in MLS phase distribution and counter
zt = load('z_tmp'); z = zt.z; counter = zt.counter;
A = Ap * P_A;
AA = Amask * P_A;
T = P_A;
% Compute pulse duration, RMS duration and TR
[rf_dur,t_rms] = SSFP_pulse_params(P_A,theta,delta_1,delta_2);
TR = max([rf_dur*2 rf_dur+t_enc]);
% SAR Scaling
S = b1_act_scale*(P_A^2)*t_rms/TR;
% Peak and average power drive constraints
drmax = (sqrt(P_peak/A_amp))/P_A;
dravg = sqrt(P_av*TR/((P_A.^2*t_rms)*A_amp));
% Most constraining power constraint
D = min([drmax dravg]);
sinit = ones(Nc,1);
mls_shim = zeros(Nc,1);
len_z = length(z);
while (abs(norm(mls_shim - sinit)) > 0.1)
sinit = mls_shim;
%%% CVX Optimisation %%%%%%%%%%%
cvx_begin quiet
variable y(Nc) complex;
minimise(norm(sum((AA*y).*(conj(z)))/len_z - T));
subject to
(y'*QG*y)*S <= wbSARmax; % Global SAR Constraint
SAR(y,VOP)*S <= (lSARmax + 0i); % Local SAR Constraint
max(abs(y)) <= D; % Power Constraint
cvx_end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
mls_shim = y;
mls_tmpsoln = A*mls_shim;
z = exp(1i*(angle(mls_tmpsoln(idx))));
err = norm((abs(mls_tmpsoln(idx))-T*targmask))/norm(T*targmask);
conv = abs(norm(mls_shim - sinit));
fprintf('Convergence = %1.6f\tError = %1.4f\n',conv,err);
% Solutions do not improve much after 1st MLS iteration so force
% only one iteration for these
if counter;break;end
end
disp(['Pulse Amplitude = ' num2str(P_A) ' TR = ' num2str(TR)])
% Compute Bias of solution
Bias = abs(P_A - mean(abs(mls_tmpsoln(idx))))/P_A*100;
Error = err;
% Set counter from zero to one in first function call
counter = 1;
% Plot bias convergence
Bias_conv = [zt.Bias_conv Bias];
plot(Bias_conv);xlabel('Iteration');ylabel('Percentage Bias');grid on
drawnow;hold on
% Cost function value
fval = TR + Bias^2;
save('z_tmp','z','counter','mls_shim','Bias_conv','S','Error')
end