-
Notifications
You must be signed in to change notification settings - Fork 3
/
fitdcemri.m
executable file
·473 lines (387 loc) · 15.7 KB
/
fitdcemri.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
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
function pars = fitdcemri(varargin)
% DCE_MRI_FIT performs fitting to estimate DCE MRI parameters for the following models:
%
% 1) Linear Reference Region Model
% pars = FITDCEMRI(toi,rr,time,'lsq')
%
% 2) Non-Linear Reference Region Model
% pars = FITDCEMRI(toi,rr,time,x0,lb,ub,'RRM')
%
% 3) Linear Tofts Model
% pars = FITDCEMRI(Ctoi,Cp,time)
%
% 4) Non-linear Tofts Model
% pars = FITDCEMRI(Ctoi,Cp,time,x0,lb,ub,'Tofts')
%
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Linear Reference Region Model (LRRM):
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% pars = fitdcemri(toi,rr,time,'fit')
% /::INPUTS::\
% toi: a column vector of either delta R1 values or Contrast Agent(CA)
% Concentration over time for the Tissue of Interest
% rr: a column vector of either delta R1 values or CA Concentration
% over time for the Reference Region
% time: a column vector of time in minutes corresponding to toi and
% rr
% 'fit': a string designating the desired fitting algorithm to be
% used. Can be:
% 'robust': uses robustfit
% 'lsq': uses linear least squared fitting (mldivide)
% 'nonneg': uses linear least squared fitting with a
% non-negative constraint
% 'lasso': uses lasso
% /::OUTPUT::\
% pars: a vector with the following components
% pars(1)= Relative Ktrans (RKtrans)
% pars(2)= kep for the tissue of interest, kepTOI
% pars(3)= kep for the reference region, kepRR
% pars(4)= rsquare for the fitting
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Non-Linear Reference Region Model (NRRM):
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% pars = fitdcemri(toi,rr,time,x0,lb,ub,'RRM')
% /::INPUTS::\
% toi: a column vector of either delta R1 values or Contrast Agent(CA)
% Concentration over time for the Tissue of Interest
% rr: a column vector of either delta R1 values or CA Concentration
% over time for the Reference Region
% time: a column vector of time in minutes corresponding to toi and
% rr
% x0: a vector of best guesses for parameters from the non-linear
% least squared fitting (Ex:[RKtrans_guess,kepTOI_guess,kepRR_guess])
% lb: vector of lowerbounds for the parameters (same form as above)
% ub: vector of upperbounds for the parameters (same form as above)
% 'RRM': This string indicates that you want the non-linear Reference
% Region Model (versus the non-linear Tofts Model)
% /::OUTPUT::\
% pars: a vector with the following components
% pars(1)= Relative Ktrans (RKtrans)
% pars(2)= kep for the tissue of interest, kepTOI
% pars(3)= kep for the reference region, kepRR
% pars(4)= rsquare for the fitting
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Linear Toft's Model (LTM):
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% pars = fitdcemri(Ctoi,Cp,time)
% /::INPUTS::\
% Ctoi: a column vector of Contrast Agent(CA) Concentration over time
% for the Tissue of Interest.
% NOTE: Must be concentration, not delta R1 (deltaR1 = r1*Conc
% where r1 is the relaxivity of the CA)
% Cp: an Arterial Input Function (AIF) for the tissue in question in
% the form of a column vector equal in size to toi
% time: a column vector of time in minutes corresponding to toi and
% Cp
% /::OUTPUT::\
% pars: a vector with the following components
% pars(1)= Ktrans
% pars(2)= kep for the tissue of interest, kepTOI
% pars(3)= rsquare for the fitting
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Non-Linear Toft's Model (NLTM):
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% pars = fitdcemri(Ctoi,Cp,time,x0,lb,ub,'Tofts')
% /::INPUTS::\
% Ctoi: a column vector of Contrast Agent(CA) Concentration over time
% for the Tissue of Interest.
% NOTE: Must be concentration, not delta R1 (deltaR1 = r1*Conc
% where r1 is the relaxivity of the CA)
% Cp: an Arterial Input Function (AIF) for the tissue in question in
% the form of a column vector equal in size to toi
% time: a column vector of time in minutes corresponding to toi and
% Cp
% x0: a vector of best guesses for parameters from the non-linear
% least squared fitting (Ex:[Ktrans_guess,kepTOI_guess])
% lb: vector of lowerbounds for the parameters (same form as above)
% ub: vector of upperbounds for the parameters (same form as above)
% 'Tofts': This string indicates that you want the non-linear Tofts
% Model (versus the non-linear RRM)
% /::OUTPUT::\
% pars: a vector with the following components
% pars(1)= Ktrans
% pars(2)= kep for the tissue of interest, kepTOI
% pars(3)= rsquare for the fitting
% Authors:
% Joey DeGrandchamp Julio Cardenas-Rodriguez
% University of Arizona University of Arizona
% jdegrandchamp@email.arizona.edu cardenaj@email.arizona.edu
%
% www.cardenaslab.org/resources
% v2.0 09/24/2015
% Check that the number of arguments is reasonable
narginchk(3,7);
if nargin == 4 %% LRRM
toi = varargin{1};
rr = varargin{2};
time = varargin{3};
fit = varargin{4};
% Make sure all vectors are column, not row
if size(toi,1) == 1
toi = toi';
elseif size(rr,1) == 1
rr = rr';
elseif size(time,1) == 1
time = time';
end
% Make sure data vectors are same size
if size(toi,1) ~= size(time,1) || size(toi,1) ~= size(rr,1) || size(rr,1) ~= size(time,1)
error(message('dce_mri_fit:toi,rr,and time vectors must all be same size'));
end
% Set up LRRM (see references)
y=toi;
x1=rr;
x2=cumtrapz(time,rr);
x3=cumtrapz(time,toi);
X=[x1,x2,-x3];
% Switch to users desired fitting algorithm
switch fit
case {'robust','Robust','ROBUST'}
% Estimate parameters using robust method
%c
B = robustfit(X,y);
B = B(2:end);
case {'lsq','LSQ','least-squared'}
% Estimate parameters using LLSQ method (mldivide)
B=X\y;
case {'nonneg','Nonneg','NONNEG','non-neg'}
% Estimate parameters using LLSQ method with non-negative constraint
B=lsqnonneg(X,y);
case {'lasso','Lasso','LASSO'}
% Estimate parameters using lasso
[b,FitInfo]=lasso(X,y);
[~,ii]=min(FitInfo.MSE);
B=b(:,ii);
otherwise
error('dce_mri_fit:Invalid fit option for the LRRM');
end
% Store each parameter of output pars (see reference) ##
pars(1,1)=B(1); %RKtrans;
pars(2,1)=B(3); %kepTOI;
pars(3,1)=(B(2))*B(1)^-1; %kepRR;
% Calculate predicted curve for the TOI
prediction=X*B; %#ok<*NBRAK>
% Calculate the the RMSE for the predicted curve and allocate to pars
pars(4,1)=rsquare(y,prediction);
elseif nargin == 5 || nargin == 6
error('dce_mri_fit:Invalid number of arguments');
elseif nargin == 7 %% Non-linear fits (RRM/Tofts)
model = varargin{7};
switch model
case {'RRM','NLRRM','NRRM','rrm','nlrrm','nrrm'}
toi = varargin{1};
rr = varargin{2};
time = varargin{3};
x0 = varargin{4};
lb = varargin{5};
ub = varargin{6};
% Make sure all vectors are column, not row
if size(toi,1) == 1
toi = toi';
elseif size(rr,1) == 1
rr = rr';
elseif size(time,1) == 1
time = time';
end
% Make sure data vectors are same size
if size(toi,1) ~= size(time,1) || size(toi,1) ~= size(rr,1) || size(rr,1) ~= size(time,1)
error('dce_mri_fit:toi,rr,and time vectors must all be same size');
end
% Set options for non-linear LSQ fitting
S=optimset; S.Algorithm='trust-region-reflective'; S.Display='off';
S.TolFun=1e-3; S.TolX=1e-3;
S.MaxIter=1000;
% Perform the fitting (the function "rrm" is in the nested functions below)
[B,~,~] = lsqcurvefit(@rrm,x0,[time,rr],toi,lb,ub,S);
% Store each parameter (see reference)
pars(1,1)=B(1); %RKtrans
pars(2,1)=B(2); %kepTOI
pars(3,1)=B(3); %kepRR
pred=rrm(B,[time,rr]);
pars(4)=rsquare(toi,pred);
case {'Tofts','tofts','NLTM','nltm'}
toi = varargin{1};
Cp = varargin{2};
time = varargin{3};
x0 = varargin{4};
lb = varargin{5};
ub = varargin{6};
% Make sure all vectors are column, not row
if size(toi,1) == 1
toi = toi';
elseif size(Cp,1) == 1
Cp = Cp';
elseif size(time,1) == 1
time = time';
end
% Make sure data vectors are same size
if size(toi,1) ~= size(time,1) || size(toi,1) ~= size(Cp,1) || size(Cp,1) ~= size(time,1)
error(message('dce_mri_fit:toi,Cp,and time vectors must all be same size'));
end
% Set options for non-linear LSQ fitting
S=optimset; S.Algorithm='trust-region-reflective'; S.Display='off';
S.TolFun=1e-3; S.TolX=1e-3;
S.MaxIter=1000;
% Perform the fitting (the function "rrm" is in the nested functions below)
[B,~,~] = lsqcurvefit(@tofts,x0,[time,Cp],toi,lb,ub,S);
% Store each parameter (see reference)
pars(1,1)=B(1); %Ktrans
pars(2,1)=B(2); %kepTOI
pred=tofts(B,[time,Cp]);
pars(3,1)=rsquare(toi,pred);
otherwise
error('dce_mri_fit:Must choose ''Tofts'' or ''RRM'' for last option');
end
elseif nargin ==3 %% LTM
toi = varargin{1};
Cp = varargin{2};
time = varargin{3};
% Make sure all vectors are column, not row
if size(toi,1) == 1
toi = toi';
elseif size(Cp,1) == 1
Cp = Cp';
elseif size(time,1) == 1
time = time';
end
% Make sure data vectors are same size
if size(toi,1) ~= size(time,1) || size(toi,1) ~= size(Cp,1) || size(Cp,1) ~= size(time,1)
error(message('dce_mri_fit:toi,Cp,and time vectors must all be same size'));
end
y=toi;
x1=cumtrapz(time,Cp);
x2=cumtrapz(time,toi);
X=[x1,-x2,ones(size(y))];
% Estimate parameters using a LLSQ with non-negative constraint
[B,~,~]=lsqnonneg(X,y);
% Calculate each parameter of output pars (see reference) ##
pars(1,1)=B(1); %Ktrans;
pars(2,1)=B(2); %kep;
% Calculate predicted curve for the TOI
prediction=X*B; %#ok<*NBRAK>
% Calculate the the r-square for the predicted curve and allocated to pars
pars(3,1)=rsquare(y,prediction);
end
end
%% ####### Internal Functions ##########
function [r2, rmse] = rsquare(y,f,varargin)
% Compute coefficient of determination of data fit model and RMSE
%
% [r2 rmse] = rsquare(y,f)
% [r2 rmse] = rsquare(y,f,c)
%
% RSQUARE computes the coefficient of determination (R-square) value from
% actual data Y and model data F. The code uses a general version of
% R-square, based on comparing the variability of the estimation errors
% with the variability of the original values. RSQUARE also outputs the
% root mean squared error (RMSE) for the user's convenience.
%
% Note: RSQUARE ignores comparisons involving NaN values.
%
% INPUTS
% Y : Actual data
% F : Model fit
%
% OPTION
% C : Constant term in model
% R-square may be a questionable measure of fit when no
% constant term is included in the model.
% [DEFAULT] TRUE : Use traditional R-square computation
% FALSE : Uses alternate R-square computation for model
% without constant term [R2 = 1 - NORM(Y-F)/NORM(Y)]
%
% OUTPUT
% R2 : Coefficient of determination
% RMSE : Root mean squared error
%
% EXAMPLE
% x = 0:0.1:10;
% y = 2.*x + 1 + randn(size(x));
% p = polyfit(x,y,1);
% f = polyval(p,x);
% [r2 rmse] = rsquare(y,f);
% figure; plot(x,y,'b-');
% hold on; plot(x,f,'r-');
% title(strcat(['R2 = ' num2str(r2) '; RMSE = ' num2str(rmse)]))
%
% Jered R Wells
% 11/17/11
% jered [dot] wells [at] duke [dot] edu
%
% v1.2 (02/14/2012)
%
% Thanks to John D'Errico for useful comments and insight which has helped
% to improve this code. His code POLYFITN was consulted in the inclusion of
% the C-option (REF. File ID: #34765).
if isempty(varargin); c = true;
elseif length(varargin)>1; error 'Too many input arguments';
elseif ~islogical(varargin{1}); error 'C must be logical (TRUE||FALSE)'
else c = varargin{1};
end
% Compare inputs
if ~all(size(y)==size(f)); error 'Y and F must be the same size'; end
% Check for NaN
tmp = ~or(isnan(y),isnan(f));
y = y(tmp);
f = f(tmp);
if c; r2 = max(0,1 - sum((y(:)-f(:)).^2)/sum((y(:)-mean(y(:))).^2));
else r2 = 1 - sum((y(:)-f(:)).^2)/sum((y(:)).^2);
if r2<0
% http://web.maths.unsw.edu.au/~adelle/Garvan/Assays/GoodnessOfFit.html
warning('Consider adding a constant term to your model') %#ok<WNTAG>
r2 = 0;
end
end
rmse = sqrt(mean((y(:) - f(:)).^2));
end
function C_toi = rrm(beta,X)
time=X(:,1);
C_rr=X(:,2);
RKtrans= beta(1); %ktrans_toi / ktrans_rr ; RKtrans
kepRR= beta(3); %ktrans_rr / ve_rr; KepRR
kepTOI= beta(2); %ktrans_toi / ve_toi; KepTOI
R4= RKtrans*(kepRR-kepTOI);
C_toi=RKtrans.*C_rr + R4.*convolution(C_rr,exp(-kepTOI.*time));
end
function C_toi = tofts(beta,X)
time=X(:,1);
Cp = X(:,2);
Ktrans = beta(1);
kepTOI = beta(2);
C_toi = Ktrans*convolution(Cp,exp(-kepTOI.*time));
end
function c = convolution(a, b, shape)
%CONVOLUTION MODIFIED BY JULIO CARDENAS, MAY OF 2011.
% SAME THAN CONV BUT RETURN A VECTOR FROM a(1) to a(end), not the central
% section as described for the usual convolution function.
%
if ~isvector(a) || ~isvector(b)
error(message('MATLAB:conv:AorBNotVector'));
end
if nargin < 3
shape = 'full';
end
if ~ischar(shape)
error(message('MATLAB:conv:unknownShapeParameter'));
end
% compute as if both inputs are column vectors
[rows,~]=size(a);
c = conv2(a(:),b(:),shape);
c=c(1:rows);
% restore orientation
if shape(1) == 'f'
if length(a) > length(b)
if size(a,1) == 1 %row vector
c = c.';
end
else
if size(b,1) == 1 %row vector
c = c.';
end
end
else
if size(a,1) == 1 %row vector
c = c.';
end
end
end