Skip to content

Commit 544464f

Browse files
markmikkelsenmarkmikkelsen
committed
Refactor segmentation and QA metrics in Seg.m
Improved tissue segmentation by introducing a probability threshold to reduce partial volume effects and updated variable names for clarity. Enhanced MRI quality metrics calculation using more robust statistics (median, MAD) and added additional metrics (FBER, WM2MAX, SNR variants). Updated function signatures and internal logic for consistency in both Seg.m and GannetSegment.m.
1 parent 3a690e9 commit 544464f

File tree

2 files changed

+114
-81
lines changed

2 files changed

+114
-81
lines changed

‎CoRegStandAlone/Seg.m‎

Lines changed: 112 additions & 79 deletions
Original file line numberDiff line numberDiff line change
@@ -40,6 +40,8 @@
4040
vox = MRS_struct.p.vox(1);
4141
kk = 1;
4242
setup_spm = 1;
43+
prob_threshold = 0.9; % threshold to reduce partial volume effect and
44+
% improve accuracy of tissue probability maps
4345

4446
for ii = 1:length(MRS_struct.metabfile)
4547

@@ -50,13 +52,13 @@
5052

5153
% Check to see if segmentation has already been done (and all
5254
% probability tissue maps are present)
53-
tmp = {[T1dir filesep 'c1' T1name T1ext]
54-
[T1dir filesep 'c2' T1name T1ext]
55-
[T1dir filesep 'c3' T1name T1ext]
56-
[T1dir filesep 'c6' T1name T1ext]};
57-
filesExist = zeros(1,length(tmp));
58-
for jj = 1:length(tmp)
59-
filesExist(jj) = exist(tmp{jj}, 'file');
55+
tissue_maps = {[T1dir filesep 'c1' T1name T1ext]
56+
[T1dir filesep 'c2' T1name T1ext]
57+
[T1dir filesep 'c3' T1name T1ext]
58+
[T1dir filesep 'c6' T1name T1ext]};
59+
filesExist = zeros(1,length(tissue_maps));
60+
for jj = 1:length(tissue_maps)
61+
filesExist(jj) = exist(tissue_maps{jj}, 'file');
6062
end
6163
if ~all(filesExist)
6264
if setup_spm
@@ -88,89 +90,120 @@
8890
GM = [T1dir filesep 'c1' T1name T1ext];
8991
WM = [T1dir filesep 'c2' T1name T1ext];
9092
CSF = [T1dir filesep 'c3' T1name T1ext];
91-
air = [T1dir filesep 'c6' T1name T1ext];
93+
BG = [T1dir filesep 'c6' T1name T1ext];
9294

93-
GMvol = spm_vol(GM);
94-
WMvol = spm_vol(WM);
95-
CSFvol = spm_vol(CSF);
96-
airvol = spm_vol(air);
97-
98-
% Segmentation quality metrics (Chua et al. JMRI, 2009; Ganzetti et
99-
% al. Front. Neuroinform., 2016; Esteban et al. PLOS One, 2017)
95+
GM_vol = spm_vol(GM);
96+
WM_vol = spm_vol(WM);
97+
CSF_vol = spm_vol(CSF);
98+
BG_vol = spm_vol(BG);
99+
100+
% MRIQC image quality metrics (Esteban et al., 2017,doi:10.1371/journal.pone.0184661;
101+
% also see: Chua et al., 2009, doi:10.1002/jmri.21768; Ganzetti et al., 2016,
102+
% doi:10.3389/fninf.2016.00010)
100103
T1 = spm_vol(struc);
101104
T1_tmp = T1.private.dat(:,:,:);
102-
103-
WMvol_tmp = WMvol.private.dat(:,:,:);
104-
WMvol_tmp(WMvol_tmp < 0.9) = NaN;
105-
WMvol_thresh = WMvol_tmp .* T1_tmp;
106-
WMvol_thresh = WMvol_thresh(:);
107-
108-
GMvol_tmp = GMvol.private.dat(:,:,:);
109-
GMvol_tmp(GMvol_tmp < 0.9) = NaN;
110-
GMvol_thresh = GMvol_tmp .* T1_tmp;
111-
GMvol_thresh = GMvol_thresh(:);
112-
113-
airvol_tmp = airvol.private.dat(:,:,:);
114-
airvol_tmp(airvol_tmp < 0.9) = NaN;
115-
airvol_thresh = airvol_tmp .* T1_tmp;
116-
airvol_thresh = airvol_thresh(:);
117-
118-
MRS_struct.out.QA.CV_WM(ii) = std(WMvol_thresh, 'omitnan') / mean(WMvol_thresh, 'omitnan');
119-
MRS_struct.out.QA.CV_GM(ii) = std(GMvol_thresh, 'omitnan') / mean(GMvol_thresh, 'omitnan');
120-
MRS_struct.out.QA.CJV(ii) = (std(WMvol_thresh, 'omitnan') + std(GMvol_thresh, 'omitnan')) ...
121-
/ abs(mean(WMvol_thresh, 'omitnan') - mean(GMvol_thresh, 'omitnan'));
122-
MRS_struct.out.QA.CNR(ii) = abs(mean(WMvol_thresh, 'omitnan') - mean(GMvol_thresh, 'omitnan')) / ...
123-
sqrt(var(airvol_thresh, 'omitnan') + var(WMvol_thresh, 'omitnan') + var(GMvol_thresh, 'omitnan'));
124-
105+
106+
WM_vol_tmp = WM_vol.private.dat(:,:,:);
107+
WM_vol_tmp(WM_vol_tmp < prob_threshold) = 0;
108+
T1_WM = WM_vol_tmp .* T1_tmp;
109+
T1_WM = T1_WM(:);
110+
T1_WM = T1_WM(T1_WM > 0); % include only nonzero voxels
111+
112+
GM_vol_tmp = GM_vol.private.dat(:,:,:);
113+
GM_vol_tmp(GM_vol_tmp < prob_threshold) = 0;
114+
T1_GM = GM_vol_tmp .* T1_tmp;
115+
T1_GM = T1_GM(:);
116+
T1_GM = T1_GM(T1_GM > 0);
117+
118+
CSF_vol_tmp = CSF_vol.private.dat(:,:,:);
119+
CSF_vol_tmp(CSF_vol_tmp < prob_threshold) = 0;
120+
T1_CSF = CSF_vol_tmp .* T1_tmp;
121+
T1_CSF = T1_CSF(:);
122+
T1_CSF = T1_CSF(T1_CSF > 0);
123+
124+
BG_vol_tmp = BG_vol.private.dat(:,:,:);
125+
BG_vol_tmp(BG_vol_tmp < prob_threshold) = 0;
126+
T1_BG = BG_vol_tmp .* T1_tmp;
127+
T1_BG = T1_BG(:);
128+
T1_BG = T1_BG(T1_BG > 0);
129+
130+
head_vol_tmp = 1 - BG_vol.private.dat(:,:,:);
131+
head_vol_tmp(head_vol_tmp < prob_threshold) = 0;
132+
T1_head = head_vol_tmp .* T1_tmp;
133+
T1_head = T1_head(:);
134+
T1_head = T1_head(T1_head > 0);
135+
136+
MRS_struct.out.QA.CV.WM(ii) = mad(T1_WM,1) / median(T1_WM);
137+
MRS_struct.out.QA.CV.GM(ii) = mad(T1_GM,1) / median(T1_GM);
138+
MRS_struct.out.QA.CV.CSF(ii) = mad(T1_CSF,1) / median(T1_CSF);
139+
MRS_struct.out.QA.CJV(ii) = (mad(T1_WM,1) + mad(T1_GM,1)) / abs(median(T1_WM) - median(T1_GM));
140+
MRS_struct.out.QA.CNR(ii) = abs(median(T1_WM) - median(T1_GM)) / sqrt(std(T1_WM).^2 + std(T1_GM).^2 + std(T1_BG).^2);
141+
125142
T1_tmp = T1_tmp(:);
126143
n_vox = numel(T1_tmp);
127144
efc_max = n_vox * (1/sqrt(n_vox)) * log(1/sqrt(n_vox));
128145
b_max = sqrt(sum(T1_tmp.^2));
129-
MRS_struct.out.QA.EFC(ii) = (1/efc_max) .* sum((T1_tmp ./ b_max) .* log((T1_tmp + eps) ./ b_max));
130-
146+
MRS_struct.out.QA.EFC(ii) = (1/efc_max) .* sum((T1_tmp / b_max) .* log((T1_tmp + eps) / b_max));
147+
148+
MRS_struct.out.QA.FBER(ii) = median(abs(T1_head).^2) / median(abs(T1_BG).^2);
149+
MRS_struct.out.QA.WM2MAX(ii) = median(T1_WM) / prctile(T1_tmp, 99.95);
150+
151+
MRS_struct.out.QA.SNR.WM(ii) = median(T1_WM) / (std(T1_WM) * sqrt(numel(T1_WM) / (numel(T1_WM) - 1)));
152+
MRS_struct.out.QA.SNR.GM(ii) = median(T1_GM) / (std(T1_GM) * sqrt(numel(T1_GM) / (numel(T1_GM) - 1)));
153+
MRS_struct.out.QA.SNR.CSF(ii) = median(T1_CSF) / (std(T1_CSF) * sqrt(numel(T1_CSF) / (numel(T1_CSF) - 1)));
154+
MRS_struct.out.QA.SNR.total(ii) = mean([MRS_struct.out.QA.SNR.WM(ii) MRS_struct.out.QA.SNR.GM(ii) MRS_struct.out.QA.SNR.CSF(ii)]);
155+
156+
MRS_struct.out.QA.SNR_D.WM(ii) = median(T1_WM) / (sqrt(2 / (4 - pi)) * mad(T1_BG,1));
157+
MRS_struct.out.QA.SNR_D.GM(ii) = median(T1_GM) / (sqrt(2 / (4 - pi)) * mad(T1_BG,1));
158+
MRS_struct.out.QA.SNR_D.CSF(ii) = median(T1_CSF) / (sqrt(2 / (4 - pi)) * mad(T1_BG,1));
159+
MRS_struct.out.QA.SNR_D.total(ii) = mean([MRS_struct.out.QA.SNR_D.WM(ii) MRS_struct.out.QA.SNR_D.GM(ii) MRS_struct.out.QA.SNR_D.CSF(ii)]);
160+
131161
% Loop over voxels if PRIAM
132162
for kk = 1:length(vox)
133-
134-
voxmaskvol = spm_vol(cell2mat(MRS_struct.mask.(vox{kk}).fname(ii)));
135-
[a,b,c] = fileparts(voxmaskvol.fname);
136-
163+
164+
vox_mask_vol = spm_vol(cell2mat(MRS_struct.mask.(vox{kk}).fname(ii)));
165+
[a,b,c] = fileparts(vox_mask_vol.fname);
166+
137167
% GM
138-
O_GMvox.fname = fullfile(a, [b '_GM' c]);
139-
O_GMvox.descrip = 'MRS_voxel_mask_GM';
140-
O_GMvox.dim = voxmaskvol.dim;
141-
O_GMvox.dt = voxmaskvol.dt;
142-
O_GMvox.mat = voxmaskvol.mat;
143-
GM_voxmask_vol = GMvol.private.dat(:,:,:) .* voxmaskvol.private.dat(:,:,:);
144-
O_GMvox = spm_write_vol(O_GMvox, GM_voxmask_vol);
168+
GM_vox.fname = fullfile(a, [b '_GM' c]);
169+
GM_vox.descrip = 'MRS_voxel_mask_GM';
170+
GM_vox.dim = vox_mask_vol.dim;
171+
GM_vox.dt = vox_mask_vol.dt;
172+
GM_vox.mat = vox_mask_vol.mat;
173+
GM_vox_mask_vol = GM_vol.private.dat(:,:,:) .* vox_mask_vol.private.dat(:,:,:);
174+
GM_vox = spm_write_vol(GM_vox, GM_vox_mask_vol);
145175

146176
% WM
147-
O_WMvox.fname = fullfile(a, [b '_WM' c]);
148-
O_WMvox.descrip = 'MRS_voxel_mask_WM';
149-
O_WMvox.dim = voxmaskvol.dim;
150-
O_WMvox.dt = voxmaskvol.dt;
151-
O_WMvox.mat = voxmaskvol.mat;
152-
WM_voxmask_vol = WMvol.private.dat(:,:,:) .* voxmaskvol.private.dat(:,:,:);
153-
O_WMvox = spm_write_vol(O_WMvox, WM_voxmask_vol);
177+
WM_vox.fname = fullfile(a, [b '_WM' c]);
178+
WM_vox.descrip = 'MRS_voxel_mask_WM';
179+
WM_vox.dim = vox_mask_vol.dim;
180+
WM_vox.dt = vox_mask_vol.dt;
181+
WM_vox.mat = vox_mask_vol.mat;
182+
WM_vox_mask_vol = WM_vol.private.dat(:,:,:) .* vox_mask_vol.private.dat(:,:,:);
183+
WM_vox = spm_write_vol(WM_vox, WM_vox_mask_vol);
154184

155185
% CSF
156-
O_CSFvox.fname = fullfile(a, [b '_CSF' c]);
157-
O_CSFvox.descrip = 'MRS_voxel_mask_CSF';
158-
O_CSFvox.dim = voxmaskvol.dim;
159-
O_CSFvox.dt = voxmaskvol.dt;
160-
O_CSFvox.mat = voxmaskvol.mat;
161-
CSF_voxmask_vol = CSFvol.private.dat(:,:,:) .* voxmaskvol.private.dat(:,:,:);
162-
O_CSFvox = spm_write_vol(O_CSFvox, CSF_voxmask_vol);
186+
CSF_vox.fname = fullfile(a, [b '_CSF' c]);
187+
CSF_vox.descrip = 'MRS_voxel_mask_CSF';
188+
CSF_vox.dim = vox_mask_vol.dim;
189+
CSF_vox.dt = vox_mask_vol.dt;
190+
CSF_vox.mat = vox_mask_vol.mat;
191+
CSF_vox_mask_vol = CSF_vol.private.dat(:,:,:) .* vox_mask_vol.private.dat(:,:,:);
192+
CSF_vox = spm_write_vol(CSF_vox, CSF_vox_mask_vol);
163193

164194
% 3. Calculate a CSF-corrected i.u. value and output it to the structure
165195

166-
GMsum = sum(sum(sum(O_GMvox.private.dat(:,:,:))));
167-
WMsum = sum(sum(sum(O_WMvox.private.dat(:,:,:))));
168-
CSFsum = sum(sum(sum(O_CSFvox.private.dat(:,:,:))));
169-
170-
fGM = GMsum / (GMsum + WMsum + CSFsum);
171-
fWM = WMsum / (GMsum + WMsum + CSFsum);
172-
fCSF = CSFsum / (GMsum + WMsum + CSFsum);
173-
196+
GM_vox_n = GM_vox.private.dat(:,:,:);
197+
GM_sum = sum(GM_vox_n(GM_vox_n > prob_threshold));
198+
WM_vox_n = WM_vox.private.dat(:,:,:);
199+
WM_sum = sum(WM_vox_n(WM_vox_n > prob_threshold));
200+
CSF_vox_n = CSF_vox.private.dat(:,:,:);
201+
CSF_sum = sum(CSF_vox_n(CSF_vox_n > prob_threshold));
202+
203+
fGM = GM_sum / (GM_sum + WM_sum + CSF_sum);
204+
fWM = WM_sum / (GM_sum + WM_sum + CSF_sum);
205+
fCSF = CSF_sum / (GM_sum + WM_sum + CSF_sum);
206+
174207
MRS_struct.out.(vox{kk}).tissue.fGM(ii) = fGM;
175208
MRS_struct.out.(vox{kk}).tissue.fWM(ii) = fWM;
176209
MRS_struct.out.(vox{kk}).tissue.fCSF(ii) = fCSF;
@@ -247,7 +280,7 @@
247280
end
248281

249282
% Plot segmented voxel as a montage
250-
MRS_struct.mask.(vox{kk}).img_montage{ii} = PlotSegmentedVoxels(struc, voxoff, voxmaskvol, O_GMvox, O_WMvox, O_CSFvox);
283+
MRS_struct.mask.(vox{kk}).img_montage{ii} = PlotSegmentedVoxels(struc, voxoff, vox_mask_vol, GM_vox, WM_vox, CSF_vox);
251284

252285
if strcmp(MRS_struct.p.vendor,'Siemens_rda')
253286
[~,tmp,tmp2] = fileparts(MRS_struct.metabfile{ii*2-1});
@@ -348,13 +381,13 @@
348381
end
349382

350383

351-
function img_montage = PlotSegmentedVoxels(struc, voxoff, voxmaskvol, O_GMvox, O_WMvox, O_CSFvox)
384+
function img_montage = PlotSegmentedVoxels(struc, voxoff, vox_mask_vol, GM_vox, WM_vox, CSF_vox)
352385

353386
img_t = flipud(voxel2world_space(spm_vol(struc), voxoff));
354-
mask_t = flipud(voxel2world_space(voxmaskvol, voxoff));
355-
mask_t_GM = flipud(voxel2world_space(O_GMvox, voxoff));
356-
mask_t_WM = flipud(voxel2world_space(O_WMvox, voxoff));
357-
mask_t_CSF = flipud(voxel2world_space(O_CSFvox, voxoff));
387+
mask_t = flipud(voxel2world_space(vox_mask_vol, voxoff));
388+
mask_t_GM = flipud(voxel2world_space(GM_vox, voxoff));
389+
mask_t_WM = flipud(voxel2world_space(WM_vox, voxoff));
390+
mask_t_CSF = flipud(voxel2world_space(CSF_vox, voxoff));
358391

359392
w_t = zeros(size(img_t));
360393

‎GannetSegment.m‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -486,10 +486,10 @@
486486
end
487487

488488

489-
function img_montage = PlotSegmentedVoxels(struc, voxoff, voxmaskvol, GM_vox, WM_vox, CSF_vox, font_size_adj)
489+
function img_montage = PlotSegmentedVoxels(struc, voxoff, vox_mask_vol, GM_vox, WM_vox, CSF_vox, font_size_adj)
490490

491491
img_t = flipud(voxel2world_space(spm_vol(struc), voxoff));
492-
mask_t = flipud(voxel2world_space(voxmaskvol, voxoff));
492+
mask_t = flipud(voxel2world_space(vox_mask_vol, voxoff));
493493
mask_t_GM = flipud(voxel2world_space(GM_vox, voxoff));
494494
mask_t_WM = flipud(voxel2world_space(WM_vox, voxoff));
495495
mask_t_CSF = flipud(voxel2world_space(CSF_vox, voxoff));

0 commit comments

Comments
 (0)