-
Notifications
You must be signed in to change notification settings - Fork 2
/
a01_DEMO_Linear.m
128 lines (100 loc) · 3.79 KB
/
a01_DEMO_Linear.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
% Quick demo using the MATLAB implementation of the CLRRM.
% Some data is simulated and then the LRRM and CLRRM are fitted to it
% When running this script, please select the "change folder" option if prompted
%%
clearvars
addpath('./mfiles')
% Directory to save data which will be used by Julia for CNRRM
outDir = './data';
if ~exist(outDir,'dir')
mkdir(outDir);
end
%% Simulate some data
rng(12345); % Set arbitrary seed for reproducibility
% Simulated reference region properties
KTransRR = 0.1;
kepRR = 1.0;
veRR = KTransRR ./ kepRR;
% Tissue of interest properties
KTrans = 0.25;
kep = 0.625;
ve = KTrans ./ kep;
% Temporal Resolution (in seconds)
tRes = 5;
% Concentration-Noise Ratio (lower value -> more noise)
CNR = 10;
nRep = 1000; % Number of replications (copies) in which noise is added
% Define the time-span
tEnd = 10; % Duration of acquisition (in minutes - most other variables are in seconds)
tInject = 60; % Time of injection (in seconds)
initialTRes = 0.1; % Simulate at 0.1 s initially
t = initialTRes:initialTRes:tEnd*60;
t = t'/60; % Convert time from seconds into minutes
% Generate the AIF - use the injection time defined earlier
Cp = ParkerAif(t,t(tInject./initialTRes));
% Simulate the reference regions concentration-time curve
Crr = TrapzKety(Cp,[KTransRR kepRR],t)';
% Simulate the tissue of interest's concentration-time curve
Ct = TrapzKety(Cp,[KTrans kep],t)';
% Add noise to Ct
sigmaC = max(Ct(:))/CNR; % Determine sigma for the noise
Ct = repmat(Ct,[1 nRep]); % Make the copies of Ct
Ct = Ct + sigmaC * randn(size(Ct)); % Add gaussian noise to each Ct copy
% Downsample
t = downsample(t,tRes/initialTRes);
Ct = downsample(Ct,tRes/initialTRes);
Crr = downsample(Crr,tRes/initialTRes);
%% Do the fitting
% The implemented LRRM and CLRRM require that Ct be nT-by-nV matrix, where
% nT is number of frames, and nV is number of voxels
% The other two inputs are Crr (concentration in RefRegion) and t (time, in minutes)
tic
[estParamsLRRM, resnormLRRM] = LRRM(Ct,Crr,t);
runtimeLRRM=toc;
tic
[estParamsCLRRM, resnormCLRRM, estKepRR, stdKepRR] = CLRRM(Ct,Crr,t);
runtimeCLRRM=toc;
% The format of estParams is:
% [ KTrans/KTransRR , ve/veRR , kep ]
% So, if we knew KTransRR and veRR, then:
% estKTrans = KTransRR * estParams(:,1);
% estVe = veRR * estParams(:,2);
% Otherwise, the relative parameters themselves could be useful on their own
%% Compute percent errors
% Percent error in relative KTrans
errEstKTrans_LRRM = PercentError(estParamsLRRM(:,1), KTrans/KTransRR);
% Percent error in relative Ve
errEstVe_LRRM = PercentError(estParamsLRRM(:,2), ve/veRR);
% Percent error in kep
errEstKep_LRRM = PercentError(estParamsLRRM(:,3), kep);
errEstKTrans_CLRRM = PercentError(estParamsCLRRM(:,1), KTrans/KTransRR);
errEstVe_CLRRM = PercentError(estParamsCLRRM(:,2), ve/veRR);
errEstKep_CLRRM = PercentError(estParamsCLRRM(:,3), kep);
%% Plot the percent errors to compare the two models
modelOrder = {'LRRM','CLRRM'};
figure
boxplot([errEstKTrans_LRRM errEstKTrans_CLRRM])
rLine = refline([0 0]);
set(rLine,'LineStyle',':')
set(gca,'XTick',[1:2], 'XTickLabel',modelOrder)
ylim([-50 50])
ylabel('Percent error in relative K^{Trans}')
title('Percent error in relative K^{Trans}')
figure
boxplot([errEstVe_LRRM errEstVe_CLRRM])
rLine = refline([0 0]);
set(rLine,'LineStyle',':')
set(gca,'XTick',[1:2], 'XTickLabel',modelOrder)
ylim([-20 20])
ylabel('Percent error in relative v_e')
title('Percent error in relative v_e')
figure
boxplot([errEstKep_LRRM errEstKep_CLRRM])
rLine = refline([0 0]);
set(rLine,'LineStyle',':')
set(gca,'XTick',[1:2], 'XTickLabel',modelOrder)
ylim([-100 100])
ylabel('Percent error in k_{ep}')
title('Percent error in k_{ep}')
%% Save original data for Julia
save(fullfile(outDir,'demoData4Julia.mat'),'Ct','Crr','t','KTrans','KTransRR','ve','veRR','kep','kepRR');