-
Notifications
You must be signed in to change notification settings - Fork 2
/
Copy pathmcDESPOT_MC.m
127 lines (105 loc) · 6.22 KB
/
mcDESPOT_MC.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
%%% Performs MC simulations for box plots in Figures 5-7, 9 and SF 2. %%%
close all; clear all;
tissuetype = {'HB','WML','INT','GML'};
acquisition = 'Bouhrara';
exchangeGT = 'on';
exchangefit = 'off';
inc_kFS = 'off';
noise = 'on';
switch exchangefit
case 'on'
Params = 6;
case 'off'
Params = 5;
end
Solution_Bouhrara = zeros(length(tissuetype),1,1000,Params);
Solution_Deoni = zeros(length(tissuetype),1,1000,Params);
Solution_Wood = zeros(length(tissuetype),1,1000,Params);
Solution_Zhang = zeros(length(tissuetype),1,1000,Params);
for qq = 1:length(tissuetype)
disp(['Tissue Type: ', tissuetype{qq}])
switch exchangeGT
case 'on'
switch tissuetype{qq}
case 'HB'
T1_F = 0.45; T1_S = 1.4; T2_F = 0.015; T2_S = 0.09; M0_F = 0.15; k_FS = 8; Delta = 0; PC1 = 0 + Delta; PC2 = pi + Delta;
case 'WML'
T1_S = 1; T1_F = 0.35; T2_S = 0.080; T2_F = 0.015; k_FS = 10; M0_F = 0.25; Delta = 0; PC1 = 0 + Delta; PC2 = pi + Delta;
case 'INT'
T1_S = 1.15; T1_F = 0.4; T2_S = 0.11; T2_F = 0.02; k_FS = 7.5; M0_F = 0.175; Delta = 0; PC1 = 0 + Delta; PC2 = pi + Delta;
case 'GML'
T1_S = 1.3; T1_F = 0.45; T2_S = 0.14; T2_F = 0.025; k_FS = 5; M0_F = 0.1; Delta = 0; PC1 = 0 + Delta; PC2 = pi + Delta;
end
case 'off'
switch tissuetype{qq}
case 'HB'
T1_F = 0.45; T1_S = 1.4; T2_F = 0.015; T2_S = 0.09; M0_F = 0.15; k_FS = 0; Delta = 0; PC1 = 0 + Delta; PC2 = pi + Delta;
case 'WML'
T1_S = 1; T1_F = 0.35; T2_S = 0.080; T2_F = 0.015; k_FS = 0; M0_F = 0.25; Delta = 0; PC1 = 0 + Delta; PC2 = pi + Delta;
case 'INT'
T1_S = 1.15; T1_F = 0.4; T2_S = 0.11; T2_F = 0.02; k_FS = 0; M0_F = 0.175; Delta = 0; PC1 = 0 + Delta; PC2 = pi + Delta;
case 'GML'
T1_S = 1.3; T1_F = 0.45; T2_S = 0.14; T2_F = 0.025; k_FS = 0; M0_F = 0.1; Delta = 0; PC1 = 0 + Delta; PC2 = pi + Delta;
end
end
switch acquisition
case 'Bouhrara'
TR_SPGR = 6.5e-3; TR_SSFP = 6.5e-3; FA_SPGR = deg2rad([2 4 6 8 10 12 14 16 18 20]); FA_SSFP180 = deg2rad([2 6 14 22 30 38 46 54 62 70]); FA_SSFP0 = deg2rad([2 6 14 22 30 38 46 54 62 70]);
case 'Teixeira'
TR_SPGR = 7e-3; TR_SSFP = 7e-3; FA_SPGR = deg2rad([6 8 10 12 14 16]); FA_SSFP180 = deg2rad([15 25 35 45 55 65]); FA_SSFP0 = deg2rad([25 55]);
case 'Deoni'
TR_SPGR = 5.6e-3; TR_SSFP = 4.4e-3; FA_SPGR = deg2rad([4 5 6 7 9 11 14 18]); FA_SSFP0 = deg2rad([12 16 19 23 27 34 50 70]); FA_SSFP180 = deg2rad([12 16 19 23 27 34 50 70]);
end
switch inc_kFS
case 'on'
k_FS = 0:4:20;
end
%% Perform stochastic region contraction.
Realisations = 1000; Trials = 40000; Iterations = 30; N = 50; Runs = 1;
delete(gcp('nocreate')); c = parcluster('local'); c.NumWorkers = 32; parpool(c, c.NumWorkers);
for ii = 1:length(k_FS)
disp(['kFS Trial Number: ', num2str(k_FS(ii))])
% Ground-truth signals for mcDESPOT.
SPGR_Data = SPGR_SteadyState(FA_SPGR, TR_SPGR,'T1_S',T1_S,'T1_F',T1_F,'M0_F',M0_F,'k_FS',k_FS(ii));
SSFP_Data_0 = SSFP_SteadyState(FA_SSFP0, TR_SSFP, PC1,'T1_S',T1_S,'T2_S',T2_S,'T1_F',T1_F,'T2_F',T2_F,'M0_F',M0_F,'k_FS',k_FS(ii));
SSFP_Data_180 = SSFP_SteadyState(FA_SSFP180, TR_SSFP, PC2,'T1_S',T1_S,'T2_S',T2_S,'T1_F',T1_F,'T2_F',T2_F,'M0_F',M0_F,'k_FS',k_FS(ii));
% Concatenate SSFP signals.
SSFP_Data = [SSFP_Data_0 ; SSFP_Data_180];
% Define SNR and define wrt mean SPGR signal.
SNR = 100; Sigma = mean(SPGR_Data)/SNR;
parfor tt = 1:Realisations
switch noise
case 'on'
% Add different noise to each element.
SPGR_Data_Noisy = zeros(length(SPGR_Data),1);
for mm = 1:length(SPGR_Data)
SPGR_Data_Noisy(mm) = SPGR_Data(mm) + (normrnd(0,Sigma));
end
SSFP_Data_Noisy = zeros(length(SSFP_Data),1);
for nn = 1:length(SSFP_Data)
SSFP_Data_Noisy(nn) = SSFP_Data(nn) + (normrnd(0,Sigma));
end
% Normalise signals and concatenate.
SPGR_Data_Norm = SPGR_Data_Noisy./mean(SPGR_Data_Noisy);
SSFP_Data_Norm = SSFP_Data_Noisy./mean(SSFP_Data_Noisy);
Data = [SPGR_Data_Norm ; SSFP_Data_Norm];
case 'off'
SPGR_Data_Norm = SPGR_Data./mean(SPGR_Data);
SSFP_Data_Norm = SSFP_Data./mean(SSFP_Data);
Data = [SPGR_Data_Norm ; SSFP_Data_Norm];
end
switch exchangefit
case 'off'
[Solution_Bouhrara(qq,ii,tt,:), ~, ~] = SRC_Sim_NoExB(Trials, Iterations, N, Runs, FA_SPGR, FA_SSFP0, FA_SSFP180, TR_SPGR, TR_SSFP, Data);
[Solution_Deoni(qq,ii,tt,:), ~, ~] = SRC_Sim_NoExD(Trials, Iterations, N, Runs, FA_SPGR, FA_SSFP0, FA_SSFP180, TR_SPGR, TR_SSFP, Data);
[Solution_Wood(qq,ii,tt,:), ~, ~] = SRC_Sim_NoExW(Trials, Iterations, N, Runs, FA_SPGR, FA_SSFP0, FA_SSFP180, TR_SPGR, TR_SSFP, Data);
[Solution_Zhang(qq,ii,tt,:), ~, ~] = SRC_Sim_NoExZ(Trials, Iterations, N, Runs, FA_SPGR, FA_SSFP0, FA_SSFP180, TR_SPGR, TR_SSFP, Data);
case 'on'
[Solution_Bouhrara(qq,ii,tt,:), ~, ~] = SRC_SimB(Trials, Iterations, N, Runs, FA_SPGR, FA_SSFP0, FA_SSFP180, TR_SPGR, TR_SSFP, Data);
[Solution_Deoni(qq,ii,tt,:), ~, ~] = SRC_SimD(Trials, Iterations, N, Runs, FA_SPGR, FA_SSFP0, FA_SSFP180, TR_SPGR, TR_SSFP, Data);
[Solution_Wood(qq,ii,tt,:), ~, ~] = SRC_SimW(Trials, Iterations, N, Runs, FA_SPGR, FA_SSFP0, FA_SSFP180, TR_SPGR, TR_SSFP, Data);
[Solution_Zhang(qq,ii,tt,:), ~, ~] = SRC_SimZ(Trials, Iterations, N, Runs, FA_SPGR, FA_SSFP0, FA_SSFP180, TR_SPGR, TR_SSFP, Data);
end
end
end
end