-
Notifications
You must be signed in to change notification settings - Fork 2
/
e01_x_simAnalyzerStatic.m
95 lines (81 loc) · 3.29 KB
/
e01_x_simAnalyzerStatic.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
% Simulation Analyzer using static kepRef estimates.
% Note: This script is UNNECESSARY. Comments in code will be minimal.
% The code here was not used in the paper, but (I think) it is interesting nonetheless.
% This script uses the CLRRM with a range of kepRef values in order to see
% the effect they have on the estimates.
% In other words, this script tries to answer the question:
% "How bad will the CLRRM estimates be if the kepRef was over/underestimated by +/-50%?"
% Run time: ~470 seconds
%% Setup
clearvars
fclose('all')
addpath('./mfiles');
refName = 'refY';
% Choices are: 'refY', 'refWS', 'refP'
% The range of kepRef estimates
% If 'refY' is used, then this range goes from -50% to +50% from the true kepRef (which is 1.0 /min)
refList = {0.5,0.6,0.7,0.8,0.85,0.9,0.925,0.95,0.975,1.0,1.025,1.05,1.075,1.1,1.15,1.2,1.3,1.4,1.5};
%% Simulation Properties
simProp = SimProperties(refName);
CNR = simProp.CNR; % Range of Contrast-Noise Ratio to simulate
nVox = simProp.nVox; % Number of replications for each CNR
TRes = simProp.TRes; % Temporal resolutions (in seconds)
tDuration = simProp.tDuration; % Duration of DCE Acquisition (in minutes)
% Pharmacokinetics of tumour tissue
Kt = simProp.Kt;
ve = simProp.ve;
kep = simProp.kep;
% Pharmacokinetics of reference tissue
KtRR = simProp.KtRR;
veRR = simProp.veRR;
kepRR = simProp.kepRR;
relKt = Kt/KtRR;
relVe = ve/veRR;
%% Folder setup
dirLocation = DefaultFolders();
dataDir = fullfile(dirLocation.sim,simProp.name,'rawData');
outDir = fullfile(dirLocation.sim,simProp.name);
outFile = fullfile(dirLocation.results,['e01-simResultsStatic-' simProp.name '.csv']);
if ~exist(dirLocation.results)
mkdir(dirLocation.results)
end
%% Initialize CSV with header
hdr=['KepRef,AIF,CNR,TemporalRes,TrueRelKt,TrueRelVe,TrueKep,errKt,errVe,errKep,stdErrKt,stdErrVe,stdErrKep, meanResid, stdResid'];
outID = fopen(outFile, 'w+');
fprintf(outID, '%s\n', hdr);
%% Process the simulated data
tic
curAif = 'Default';
for indCNR = 1:length(CNR)
inFile = ['Sim-CNR-' num2str(CNR(indCNR)) '.mat'];
inFile = fullfile(dataDir,inFile);
load(inFile); % Load noisy simulation data
for indTRes = 1:length(TRes)
curTRes = TRes(indTRes);
% Downsample the high temporal resolution (1s) data to desired resolution
t = downsample(T, curTRes);
Ct = downsample(CtNoisy, curTRes);
Crr = downsample(CrrClean, curTRes);
for i=1:length(refList)
curRef = refList{i};
[pkParamsCL, residCL, meanRefKep, stdRefKep, p] = CLRRM(Ct, Crr, t, curRef);
ktError = 100*(pkParamsCL(:,1)-relKt)/relKt;
veError = 100*(pkParamsCL(:,2)-relVe)/relVe;
kepError = 100*(pkParamsCL(:,3)-kep)/kep;
meanKtErr = mean(ktError);
meanVeErr = mean(veError);
meanKepErr = mean(kepError);
stdKtErr = std(ktError);
stdVeErr = std(veError);
stdKepErr = std(kepError);
meanResid = mean(residCL(:));
stdResid = std(residCL(:));
outLine = {curRef,curAif,curCNR,curTRes,relKt,relVe,kep,meanKtErr,meanVeErr,...
meanKepErr,stdKtErr,stdVeErr,stdKepErr,meanResid,stdResid};
fprintf(outID,'%f,%s,%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,%f\n', outLine{:});
end
end
end
toc
fclose(outID);
fclose('all')