-
Notifications
You must be signed in to change notification settings - Fork 2
/
Copy pathSinglePool_HMs.m
73 lines (53 loc) · 3.74 KB
/
SinglePool_HMs.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
%%% Generates single-pool heat maps for SF 4.
close all; clear all;
% Tissue and sequence parameters.
T1 = 1; T2 = 0.1; M0 = 1; Delta = 0; PC1 = 0 + Delta; PC2 = pi + Delta;
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]);
% Ground-truth signals for mcDESPOT.
SPGR_Data = SPGR_SP_SteadyState(FA_SPGR, TR_SPGR,'T1',T1,'M0',M0);
SSFP_Data_0 = SSFP_SP_SteadyState(FA_SSFP0, TR_SSFP, PC1,'T1',T1,'T2',T2,'M0',M0);
SSFP_Data_180 = SSFP_SP_SteadyState(FA_SSFP180, TR_SSFP, PC2,'T1',T1,'T2',T2,'M0',M0);
% Concatenate SSFP signals.
SSFP_Data = [SSFP_Data_0 ; SSFP_Data_180];
%% Direct sampling of cost-function in 2D. Modify for signal fit plots.
SNR = 30; Sigma = mean(SPGR_Data)/SNR;
Upper = [1.1 1.1 0.11]; Lower = [0.9 0.9 0.09];
% Pre-normalisation.
Data_O = [SPGR_Data ; SSFP_Data_0 ; SSFP_Data_180];
Steps = 49; nTrials = 100000;
T1_Vector = linspace(Lower(1),Upper(1),Steps); M0_Vector = linspace(Lower(2),Upper(2),Steps);
T2_Vector = linspace(Lower(3),Upper(3),Steps);
P_T1 = zeros(Steps,Steps,nTrials); P_T2 = zeros(Steps,Steps,nTrials);
P_T1T2 = zeros(Steps,Steps,nTrials);
for ii = 1:Steps
disp(['Outer Step Number: ', num2str(ii), '/', num2str(Steps), '.'])
for jj = 1:Steps
parfor nn = 1:nTrials
% Draw random parameter values from a uniform distribution.
T1_Rand = (Upper(1) - Lower(1)) .* rand(1,1) + Lower(1);
M0_Rand = (Upper(2) - Lower(2)) .* rand(1,1) + Lower(2);
T2_Rand = (Upper(3) - Lower(3)) .* rand(1,1) + Lower(3);
P_T1(ii,jj,nn) = (logpdf_SinglePool([T1_Vector(ii) M0_Vector(jj) T2_Rand], Lower, Upper, FA_SPGR, FA_SSFP0, FA_SSFP180, TR_SPGR, TR_SSFP, Data_O, Sigma));
P_T2(ii,jj,nn) = (logpdf_SinglePool([T1_Rand M0_Vector(ii) T2_Vector(jj)], Lower, Upper, FA_SPGR, FA_SSFP0, FA_SSFP180, TR_SPGR, TR_SSFP, Data_O, Sigma));
P_T1T2(ii,jj,nn) = (logpdf_SinglePool([T1_Vector(ii) M0_Rand T2_Vector(jj)], Lower, Upper, FA_SPGR, FA_SSFP0, FA_SSFP180, TR_SPGR, TR_SSFP, Data_O, Sigma));
end
end
end
%% Heat map analysis and plots.
coloraxis = ([0 1]);
figure(3); subplot(1,3,1); %MaxCF_T1 = max(exp(P_T1),[],3);
imagesc([min(M0_Vector) max(M0_Vector)],[min(T1_Vector) max(T1_Vector)],MaxCF_T1_NoiseInd,coloraxis); shading interp; xlabel('M_{0}','FontSize',20); ylabel('T_{1} (s)','FontSize',20); hold on;
plot(M0,T1,'.','Color',[0 0.5 0], 'MarkerSize', 22); get(gca, 'XTick'); set(gca, 'FontSize', 20); get(gca, 'YTick'); set(gca, 'FontSize', 20); xlim([0.95 1.05])
axis square
subplot(1,3,2); %MaxCF_T2 = max(exp(P_T2),[],3);
imagesc([min(T2_Vector) max(T2_Vector)],[min(M0_Vector) max(M0_Vector)],MaxCF_T2_NoiseInd,coloraxis); shading interp; xlabel('T_{2} (s)','FontSize',20); ylabel('M_{0}','FontSize',20); hold on;
plot(T2,M0,'.','Color',[0 0.5 0], 'MarkerSize', 22); get(gca, 'XTick'); set(gca, 'FontSize', 20); get(gca, 'YTick'); set(gca, 'FontSize', 20); xlim([0.095 0.105])
axis square
subplot(1,3,3); %MaxCF_T1T2 = max(exp(P_T1T2),[],3);
imagesc([min(T2_Vector) max(T2_Vector)],[min(T1_Vector) max(T1_Vector)],MaxCF_T1T2_NoiseInd,coloraxis); shading interp; xlabel('T_{2} (s)','FontSize',20); ylabel('T_{1} (s)','FontSize',20); hold on;
plot(T2,T1,'.','Color',[0 0.5 0], 'MarkerSize', 22); get(gca, 'XTick'); set(gca, 'FontSize', 20); get(gca, 'YTick'); set(gca, 'FontSize', 20); xlim([0.095 0.105])
axis square
colormap(flipud(magma));
hcb = colorbar('Position',[0.916360983719776,0.300794551645857,0.013774574153008,0.450624290578888],'FontSize',14); hcb.FontSize = 18;
colorTitleHandle = get(hcb,'Title'); titleString = '\eta';
set(colorTitleHandle,'String',titleString,'FontSize',34);