-
Notifications
You must be signed in to change notification settings - Fork 2
/
Copy pathdsc_fitGamma2listConc.m
80 lines (68 loc) · 3.86 KB
/
dsc_fitGamma2listConc.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
function [signalsConcFit,time] = dsc_fitGamma2listConc(FnameMatFileConcTime,outFname, limitCropInjTime)
%UNTITLED3 Summary of this function goes here
% time: time in seconds (need to have a regular sampling)
% -------------------------------------------------------------------------
% default parameter values
% -------------------------------------------------------------------------
if nargin<3
limitCropInjTime = false;
end
% -------------------------------------------------------------------------
% load data
% -------------------------------------------------------------------------
load(FnameMatFileConcTime, 'concAllSignals', 'acqTimeRegrid', 'injTime', 'TE', 'r2GdInBlood', 'firstPassStartTime', 'firstPassEndTime');
% duration to keep after injection
upperBoundCropDelay = 1.3*(firstPassEndTime - injTime)/1000;
lowerBoundCropDelay = 4.0*(firstPassEndTime - firstPassStartTime)/1000; % with respect to first pass end
% time in seconds starting at 0
time = (acqTimeRegrid - repmat(min(acqTimeRegrid,[],2), 1, size(acqTimeRegrid,2)))/1000;
injTime = (injTime - min(acqTimeRegrid(:)))/1000;
firstPassEndTime = (firstPassEndTime - min(acqTimeRegrid(:)))/1000;
% -------------------------------------------------------------------------
% crop signals after bolus because of big differences in steady-stade
% signal
% -------------------------------------------------------------------------
[ ~, upperCropBound] = min( abs( time(1,:) - (injTime+upperBoundCropDelay) ) );
[ ~, lowerCropBound] = min( abs( time(1,:) - (firstPassEndTime-lowerBoundCropDelay) ) );
if limitCropInjTime
[ ~, injRep] = min( abs( time(1,:) - injTime ) );
lowerCropBound = max([lowerCropBound, injRep]);
end
concAllSignalsCrop = concAllSignals(:, lowerCropBound:upperCropBound);
timeCrop = time(:, lowerCropBound:upperCropBound);
% -------------------------------------------------------------------------
% fit gamma curve to each signal in array
% -------------------------------------------------------------------------
signalsConcFit = zeros(size(concAllSignals));
progressbar = waitbar(0,['Fitting ' num2str(size(concAllSignals,1)) ' signals...']);
for i_sig=1:size(concAllSignals,1)
% fit cropped signal
signalCropFit = dsc_fitGamma(timeCrop(i_sig, :), concAllSignalsCrop(i_sig, :));
% extend signal before and after cropping
signalsConcFit(i_sig, upperCropBound:end) = signalCropFit(end)*ones(1,size(concAllSignals,2)-(upperCropBound-1));
signalsConcFit(i_sig, 1:lowerCropBound) = signalCropFit(1)*ones(1,lowerCropBound);
signalsConcFit(i_sig, lowerCropBound:upperCropBound) = signalCropFit;
waitbar(i_sig/size(concAllSignalsCrop,1),progressbar, [num2str(round(100*i_sig/size(concAllSignalsCrop,1),2)) ' % of signals processed...'])
end
delete(progressbar);
% -------------------------------------------------------------------------
% calculate relative metrics
% -------------------------------------------------------------------------
PR = 100*(1 - min(exp(-r2GdInBlood*TE*concAllSignals/1000), [], 2));
BAT = zeros(size(concAllSignals,1),1);
TT = zeros(size(concAllSignals,1),1);
TTP = zeros(size(concAllSignals,1),1);
rTTP = zeros(size(concAllSignals,1),1);
rBV = zeros(size(concAllSignals,1),1);
rBF = zeros(size(concAllSignals,1),1);
rMTT = zeros(size(concAllSignals,1),1);
for i_sig=1:size(concAllSignals,1)
[BAT(i_sig), TT(i_sig), TTP(i_sig), rTTP(i_sig)] = dsc_calculate_BAT_TT_TTP_rTTP(signalsConcFit(i_sig, :),time(i_sig, :),injTime);
[rBV(i_sig), rBF(i_sig), rMTT(i_sig)] = dsc_calculate_rBV_rBF_rMTT(signalsConcFit(i_sig, :), time(i_sig, :));
end
% -------------------------------------------------------------------------
% save results
% -------------------------------------------------------------------------
% matfile
save([outFname '.mat'], 'signalsConcFit', 'injTime', 'TE', 'r2GdInBlood', 'time', 'acqTimeRegrid', 'PR', 'BAT', 'TT', 'TTP', 'rTTP', 'rBV', 'rBF', 'rMTT');
end