-
Notifications
You must be signed in to change notification settings - Fork 7
/
Copy pathmagia_fit_2tcm_iterative.m
64 lines (52 loc) · 1.37 KB
/
magia_fit_2tcm_iterative.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
function [y,x_optim,e_min,vt] = magia_fit_2tcm_iterative(roi_tac,t_plasma,ca,cb,frames,varargin)
if(nargin==5)
lb = [0 0.001 0.01 0 0];
ub = [0.5 10 1 10 0.2];
elseif(nargin==6)
lb = varargin{1};
ub = [0.5 10 1 10 0.2];
elseif(nargin==7)
lb = varargin{1};
ub = varargin{2};
else
error('Too many input parameters given.');
end
if(isempty(help('optimoptions')))
opts = optimset('lsqcurvefit');
opts = optimset(opts,'Display', 'off');
opts = optimset(opts,'MaxIter',500);
else
opts = optimoptions('lsqcurvefit');
opts = optimoptions(opts,'Display', 'off');
opts = optimoptions(opts,'MaxIter',500);
end
%% Exclude plasma and blood measurements before injection
pos_idx = t_plasma >= 0;
t_plasma = t_plasma(pos_idx);
ca = ca(pos_idx);
cb = cb(pos_idx);
%% Create an artificial data point at time zero
if(t_plasma(1) > 0)
t_plasma = [0;t_plasma];
ca = [0;ca];
cb = [0;cb];
end
%% Fit the model to the data
t_pet = mean(frames,2);
if(size(roi_tac,2)~=1)
roi_tac = roi_tac';
end
fun = @(x,t) spline(t_plasma,magia_2tcm_3(x,t,ca,cb),t_pet);
n = 25;
e_min = 1e10;
for i = 1:n
x0 = unifrnd(lb,ub);
[x,e] = lsqcurvefit(fun,x0,t_plasma,roi_tac,lb,ub,opts);
if(e < e_min)
e_min = e;
x_optim = x;
end
end
y = spline(t_plasma,magia_2tcm_3(x_optim,t_plasma,ca,cb),t_pet);
vt = x_optim(2)*(1+x_optim(4));
end