-
Notifications
You must be signed in to change notification settings - Fork 0
/
calculateTwoComponentModel.m
135 lines (94 loc) · 5.69 KB
/
calculateTwoComponentModel.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
% ----------------------------------------------------------------------- %
% This function calculates and saves the values from the two-component %
% model (volume fractions) for a patient, using data from combined T2-DWI %
% with two TEs and two b-values. %
% ----------------------------------------------------------------------- %
function calculateTwoComponentModel(patientNr, patientDataFilePath, dicomFilePath, coRegistrationFilePath, resultsFilePath, noiseThreshold, useCoRegisteredImages)
% Find patient name string
patientNrString = getPatientNrString(patientNr);
% Write to output
fprintf('\nTwo-component model analysis of patient %s started %s\n', patientNrString, datetime(now,'ConvertFrom','datenum'));
% Load images
if useCoRegisteredImages
[im_shortTE_lowB,im_shortTE_highB,im_longTE_lowB,im_longTE_highB] = loadCoRegisteredImages(patientNrString,coRegistrationFilePath);
else
[im_shortTE_lowB,im_shortTE_highB,im_longTE_lowB,im_longTE_highB] = loadDicoms(patientNrString,dicomFilePath);
end
% Extract image size
[nRows, nCols, nSlices] = size(im_shortTE_lowB);
% Find x and y limits of the box ROI containing the prostate
[x_min,y_min,x_max,y_max] = findBoxROI(patientNr,patientDataFilePath);
% Specify b-values and TEs
b = [50e-3 700e-3]; % s/mm²
TE = [55 73]; % ms
% Specify fixed T2 values (globally optimized)
T2_slow = 45; % ms
T2_fast = 180; % ms
% Initialize maps
map_SF_fast = NaN(nRows, nCols, nSlices);
map_SI_0 = NaN(nRows, nCols, nSlices);
map_gof_sse = NaN(nRows, nCols, nSlices);
map_gof_rsquare = NaN(nRows, nCols, nSlices);
map_gof_dfe = NaN(nRows, nCols, nSlices);
map_gof_adjrsquare = NaN(nRows, nCols, nSlices);
map_gof_rmse = NaN(nRows, nCols, nSlices);
% Make variable that contains the number of excluded pixels
nExcludedPixels = 0;
% Start clock to measure time of model fit
tic
% Loops over slices, rows and columns to fit the two-component model
% for each voxel inside specified ROI ("box" containing the prostate)
for slice = 1:nSlices
for row = y_min:y_max % row = y
for col = x_min:x_max % col = x
% Check that the pixel values are above the threshold
% value, and ensure that the ADC and T2 values will be
% positive (increasing signal intensity with increasing
% b-value or TE is considered "unphysical")
if (im_longTE_highB(row,col,slice) > noiseThreshold) && ...
(im_shortTE_lowB(row,col,slice) > im_shortTE_highB(row,col,slice)) && ...
(im_longTE_lowB(row,col,slice) > im_longTE_highB(row,col,slice)) && ...
(im_shortTE_lowB(row,col,slice) > im_longTE_lowB(row,col,slice)) && ...
(im_shortTE_highB(row,col,slice) > im_longTE_highB(row,col,slice))
% Perform model fit using the 2x2 signal data
[fitresult, gof] = hybridfit_TCmodel(b, TE, [im_shortTE_lowB(row,col,slice) im_shortTE_highB(row,col,slice); im_longTE_lowB(row,col,slice) im_longTE_highB(row,col,slice)], [T2_slow T2_slow], [T2_fast T2_fast]);
% Fill maps with values from fit
map_SF_fast(row,col,slice) = fitresult.SF_fast;
map_SI_0(row,col,slice) = fitresult.SI_0;
map_gof_sse(row,col,slice) = gof.sse;
map_gof_rsquare(row,col,slice) = gof.rsquare;
map_gof_dfe(row,col,slice) = gof.dfe;
map_gof_adjrsquare(row,col,slice) = gof.adjrsquare;
map_gof_rmse(row,col,slice) = gof.rmse;
else
% Exclude pixel
nExcludedPixels = nExcludedPixels + 1;
end % if
end % for col
end % for row
end % for slice
% Write to output
fprintf('Time of model fit: %.f seconds.\n',toc);
% Make SF_slow map (SF_slow + SF_fast = 1)
map_SF_slow = 1 - map_SF_fast;
% Make value arrays from the maps, and remove NaN values
values_SF_slow = rmmissing(reshape(map_SF_slow,[],1));
values_SF_fast = rmmissing(reshape(map_SF_fast,[],1));
values_SI_0 = rmmissing(reshape(map_SI_0,[],1));
values_gof_sse = rmmissing(reshape(map_gof_sse,[],1));
values_gof_rsquare = rmmissing(reshape(map_gof_rsquare,[],1));
values_gof_dfe = rmmissing(reshape(map_gof_dfe,[],1));
values_gof_adjrsquare = rmmissing(reshape(map_gof_adjrsquare,[],1));
values_gof_rmse = rmmissing(reshape(map_gof_rmse,[],1));
% Save variables (and make directory if it does not already exist)
if ~isfolder([resultsFilePath patientNrString])
mkdir([resultsFilePath patientNrString])
end
if useCoRegisteredImages
save([resultsFilePath patientNrString '/resultsTwoComponentModel_coRegistered.mat'], 'map_SF_slow', 'map_SF_fast', 'map_SI_0', 'map_gof_sse', 'map_gof_rsquare', 'map_gof_dfe', 'map_gof_adjrsquare', 'map_gof_rmse');
save([resultsFilePath patientNrString '/resultsTwoComponentModel_coRegistered_allVariables.mat']);
else
save([resultsFilePath patientNrString '/resultsTwoComponentModel.mat'], 'map_SF_slow', 'map_SF_fast', 'map_SI_0', 'map_gof_sse', 'map_gof_rsquare', 'map_gof_dfe', 'map_gof_adjrsquare', 'map_gof_rmse');
save([resultsFilePath patientNrString '/resultsTwoComponentModel_allVariables.mat']);
end
end