-
Notifications
You must be signed in to change notification settings - Fork 0
/
three_comp_fit.m
88 lines (61 loc) · 2.67 KB
/
three_comp_fit.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
function [A_fit, T2s_fit, resnorm] = three_comp_fit(signal, echo_times, B0)
%%-------------------------------------------------------------------------
%% Set default settings according to relaxation time
%%-------------------------------------------------------------------------
if B0 == 3
%-------------------------------------------------------------------------
% 3 T values
%-------------------------------------------------------------------------
T2s_MW_guess = double(10e-3);
T2s_EW_guess = double(40e-3);
T2s_AW_guess = double(60e-3);
ub_MW = double(25e-3); % Nam et al. NeuroImage 2015
ub_EW = double(150e-3); % Nam et al. NeuroImage 2015
ub_AW = double(150e-3); % Nam et al. NeuroImage 2015
lb_MW = double(3e-3); % Nam et al. NeuroImage 2015
lb_EW = double(25e-3); % Nam et al. NeuroImage 2015
lb_AW = double(25e-3); % Nam et al. NeuroImage 2015
elseif B0 == 7
%-------------------------------------------------------------------------
% 7 T values
%-------------------------------------------------------------------------
T2s_MW_guess = double(6e-3);
T2s_EW_guess = double(30e-3);
T2s_AW_guess = double(40e-3);
ub_MW = double(15e-3);
ub_EW = double(150e-3);
ub_AW = double(150e-3);
lb_MW = double(3e-3);
lb_EW = double(15e-3);
lb_AW = double(15e-3);
%-------------------------------------------------------------------------
end
A_MW_guess = double(0.1*abs(signal(1)));
A_EW_guess = double(0.5*abs(signal(1)));
A_AW_guess = double(0.5*abs(signal(1)));
guess = [A_MW_guess, T2s_MW_guess, A_EW_guess, T2s_EW_guess, A_AW_guess, T2s_AW_guess];
lb_A_MW = 0;
lb_A_EW = 0;
lb_A_AW = 0;
ub_A_MW = double(0.3*signal(1));
ub_A_EW = double(1.5*signal(1));
ub_A_AW = double(1.5*signal(1));
lb = [lb_A_MW, lb_MW, lb_A_EW, lb_EW, lb_A_AW, lb_AW];
ub = [ub_A_MW, ub_MW, ub_A_EW, ub_EW, ub_A_AW, ub_AW];
%%-------------------------------------------------------------------------
%% data fitting and analysis
%%-------------------------------------------------------------------------
[result, resnorm, res, flag] = lsqnonlin(@calc_diff,guess,double(lb),double(ub));
A_fit.MW = result(1);
T2s_fit.MW = result(2);
A_fit.EW = result(3);
T2s_fit.EW = result(4);
A_fit.AW = result(5);
T2s_fit.AW = result(6);
%%-------------------------------------------------------------------------
%% function definitions
%%-------------------------------------------------------------------------
function diff = calc_diff(guess)
diff = squeeze(double(signal)) - (guess(1)*exp(-echo_times/guess(2))+guess(3)*exp(-echo_times/guess(4))+guess(5)*exp(-echo_times/guess(6)));
end
end