-
Notifications
You must be signed in to change notification settings - Fork 3
/
run_GageRR.m
executable file
·108 lines (86 loc) · 3.98 KB
/
run_GageRR.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
clear all; close all; clc;
%% Set up time components and determine AIF with Cp10 function
time=0/60:6.4/60:5.4; time=time';
Cp=Cp10(time);
%% Simulate 30 values for ktrans and ve
ktrans=.1*rand(30,1)+0.25;
ve=.1*rand(30,1)+0.4;
%% Generate enhancement curves
for j=1:30
Ct_tumor(:,j)=myToftsCt2(ktrans(j),ktrans(j)/ve(j),time,Cp)'; %#ok<*SAGROW>
end
%% Generate enhancement curve for reference region and add white gaussian noise
S_RR=myToftsCt2(.1,1,time,Cp)';
S_RR=awgn(S_RR,40,'measured');
%% Define guesses for non-linear analysis
lb=zeros(3,1); lb=lb';
ub=ones(3,1)*10; ub=ub';
%% Set initial SNR
%% Set up info for semi-quant analysis
NumBaselineImages=10;
%% Loop through analysis changing SNR by 1 each time
mySNR=[10,20];
% mySNR=5:1:40; UNCOMMENT THIS IF YOU WANT TO RUN THE ENTIRE RANGE
Total_Reps=10;
%Totla_Resp=1000; ; UNCOMMENT THIS IF YOU WANT TO RUN THE ENTIRE RANGE
for SNR=1:length(mySNR)
%% Define number of reptitions to test 30 curves for 3 days by gage analysis
for Reps=1:Total_Reps
%%Create 3 simulated enhancement curves by adding white gaussian noise to
%%original enhancement curve
for q=1:3
for r=1:30
S_toi=awgn(Ct_tumor(:,r),mySNR(SNR),'measured');
[maxValue,index] = max(S_toi);
MER(r,q) = maxValue;
TTP(r,q) = time(index);
iauc64(r,q) = trapz(time(1:11),S_toi(1:11));
slope(r,q) = MER(r,q)./TTP(r,q);
% Non-Negative Linear
B=fitdcemri(S_toi,S_RR,time,'nonneg');
RKtrans1_linear_nonNeg(r,q)=B(1);
x0=abs(randn(3,1));
%x0=B(1:3); Uncomment if you want to use initial guesses
%for NRRM from LRRM to make NRRM*
B = fitdcemri(S_toi,S_RR,time,x0,lb,ub,'RRM');
RKtrans_nonlinear_nonNeg(r,q)=B(1);
B=fitdcemri(S_toi,S_RR,time,'lsq');
RKtrans2_linear_lsq(r,q)=B(1);
end
end
m=1:30;
Parts=[m,m,m]; Parts=Parts(:);
operators=ones(1,90)';
%%Perform gage analysis for semi-quant analysis
observations_MER=MER(:);
observations_TTP=TTP(:);
observations_iauc64=iauc64(:);
observations_slope=slope(:);
T=gagerr( observations_MER,{Parts,operators},'printtable','off','printgraph','off');
MER_Repeatibility(Reps,SNR)=T(2,2);
T=gagerr( observations_TTP,{Parts,operators},'printtable','off','printgraph','off');
TTP_Repeatibility(Reps,SNR)=T(2,2);
T=gagerr( observations_iauc64,{Parts,operators},'printtable','off','printgraph','off');
iauc64_Repeatibility(Reps,SNR)=T(2,2);
T=gagerr( observations_slope,{Parts,operators},'printtable','off','printgraph','off');
slope_Repeatibility(Reps,SNR)=T(2,2);
%%Perform gage analysis for linear methods
T=gagerr( RKtrans1_linear_nonNeg(:),{Parts,operators},'printtable','off','printgraph','off');
Nonneg_Repeatibility(Reps,SNR)=T(2,2);
T=gagerr( RKtrans2_linear_lsq(:),{Parts,operators},'printtable','off','printgraph','off');
Lsq_Repeatibility(Reps,SNR)=T(2,2);
%%Perform gage analysis for nonlinear method
T=gagerr( RKtrans_nonlinear_nonNeg(:),{Parts,operators},'printtable','off','printgraph','off');
Nonlinear_Repeatibility(Reps,SNR)=T(2,2);
end
end
%% Create table and displat results
X(:,1)=median(MER_Repeatibility);
X(:,2)=median(iauc64_Repeatibility);
X(:,3)=median(TTP_Repeatibility);
X(:,4)=median(slope_Repeatibility);
X(:,5)=median(Nonneg_Repeatibility);
X(:,6)=median(Lsq_Repeatibility);
X(:,7)=median(Nonlinear_Repeatibility);
X=[mySNR',X];
array2table(X,'VariableNames',{'SNR','MER','iauc64','TTP','slope','NonNeg','Lsq','NonLinear'})