Skip to content

Commit 0e64d9c

Browse files
Remi-GauRemi-Gau
committed
add small function to quickly vizualize data content of several ROIs
1 parent 931340a commit 0e64d9c

File tree

6 files changed

+278
-6
lines changed

6 files changed

+278
-6
lines changed

‎docs/source/roi.rst‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -4,7 +4,9 @@ Region of interest
44
.. automodule:: src.roi
55

66
.. autofunction:: createRoi
7+
.. autofunction:: getPeakCoordinates
78
.. autofunction:: keepHemisphere
9+
.. autofunction:: plotDataInRoi
810
.. autofunction:: thresholdToMask
911
.. autofunction:: renameNeuroSynth
1012
.. autofunction:: resliceRoiImages

‎src/roi/getPeakCoordinates.m‎

Lines changed: 17 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -3,7 +3,18 @@
33
% This function gets the coordinates of a peak within a specified region of
44
% interest.
55
%
6-
% [worldCoord, voxelCoord, maxVal] = getPeakCoordinates(dataImage, roiImage, criticalT)
6+
% USAGE::
7+
%
8+
% [worldCoord, voxelCoord, maxVal] = getPeakCoordinates(dataImage, roiImage, threshold)
9+
%
10+
% :param dataImage:
11+
% :type dataImage: path
12+
%
13+
% :param roiImage:
14+
% :type roiImage: path
15+
%
16+
% :param threshold: threshold above which peak must be found
17+
% :type threshold: numerical
718
%
819
%
920
% (C) Copyright 2021 CPP ROI developers
@@ -14,13 +25,13 @@
1425

1526
args.addRequired('dataImage', isFile);
1627
args.addRequired('roiImage', isFile);
17-
args.addOptional('criticalT', 0, @isnumeric);
28+
args.addOptional('threshold', 0, @isnumeric);
1829

1930
args.parse(varargin{:});
2031

2132
dataImage = args.Results.dataImage;
2233
roiImage = args.Results.roiImage;
23-
criticalT = args.Results.criticalT;
34+
threshold = args.Results.criticalT;
2435

2536
voxelCoord = nan(1, 3);
2637
worldCoord = nan(1, 3);
@@ -31,11 +42,11 @@
3142
return
3243
end
3344

34-
if maxVal < criticalT
45+
if maxVal < threshold
3546
warning('getPeakCoordinates:noMaxBeyondThreshold', ...
3647
'No max value found beyond threshold %f in image:\n %s\n', ...
3748
dataImage, ...
38-
criticalT);
49+
threshold);
3950
return
4051
end
4152

@@ -51,7 +62,7 @@
5162
warning('getPeakCoordinates:severalMaxBeyondThreshold', ...
5263
'Several equal max value found beyond threshold %f in image:\n %s\n', ...
5364
dataImage, ...
54-
criticalT);
65+
threshold);
5566
return
5667
end
5768

‎src/roi/plotDataInRoi.m‎

Lines changed: 156 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,156 @@
1+
function figHandle = plotDataInRoi(varargin)
2+
%
3+
% Creates a figure showing a histogram of the content of a set of ROIs used on
4+
% a set of data files.
5+
%
6+
% USAGE::
7+
%
8+
% figHandle = plotDataInRoi(dataImages, roiImages)
9+
%
10+
% :param dataImages:
11+
% :type dataImages: path or cellstr of paths
12+
%
13+
% :param roiImages:
14+
% :type roiImages: path or cellstr of paths
15+
%
16+
%
17+
% EXAMPLE::
18+
%
19+
% mask1 = fullfile(pwd, 'label-V1_mask.nii')
20+
% mask2 = fullfile(pwd, 'label-V2_mask.nii')
21+
%
22+
% data1 = fullfile(pwd, 'label-0001_beta.nii')
23+
% data2 = fullfile(pwd, 'label-0002_beta.nii')
24+
%
25+
% mask = cellstr(cat(1, mask1, mask2));
26+
% data = cellstr(cat(1, data1, data2));
27+
%
28+
% plotDataInRoi(data, mask);
29+
%
30+
%
31+
% (C) Copyright 2022 CPP ROI developers
32+
33+
defaultNbBins = 100;
34+
35+
isFile = @(x) iscellstr(x) || exist(x, 'file') == 2;
36+
37+
args = inputParser;
38+
39+
args.addRequired('dataImages', isFile);
40+
args.addRequired('roiImages', isFile);
41+
42+
args.parse(varargin{:});
43+
44+
dataImages = args.Results.dataImages;
45+
roiImages = args.Results.roiImages;
46+
47+
if ischar(dataImages)
48+
dataImages = {dataImages};
49+
end
50+
51+
if ischar(roiImages)
52+
roiImages = {roiImages};
53+
end
54+
55+
nbRois = numel(roiImages);
56+
nbData = numel(dataImages);
57+
58+
if nbRois == 1 || nbData == 1
59+
[rows, cols] = optimizeSubplotNumber(nbRois * nbData);
60+
else
61+
rows = nbRois;
62+
cols = nbData;
63+
end
64+
65+
figHandle = figure('position', [50 50 300 * rows 300 * cols]);
66+
67+
%% collect info to adapt the graphs later on
68+
69+
% y scale
70+
maxVox = [];
71+
72+
% x axis
73+
MIN = [];
74+
MAX = [];
75+
76+
% use the same number of bins for all graphs
77+
% based on the minimum number of unique values across all the datasets
78+
nbBins = [];
79+
80+
% to plot all the modes
81+
modes = [];
82+
83+
subplotList = [];
84+
85+
idxSubplot = 1;
86+
87+
for iRoi = 1:nbRois
88+
89+
for iData = 1:nbData
90+
91+
data{idxSubplot} = spm_summarise(spm_vol(dataImages{iData}), roiImages{iRoi});
92+
93+
[~, bins] = hist(data{idxSubplot}, defaultNbBins);
94+
MAX(end + 1) = max(bins);
95+
MIN(end + 1) = min(bins);
96+
97+
% modes and nbBins work better on rounded values
98+
modes(end + 1) = mode(round(data{idxSubplot}));
99+
nbBins(end + 1) = numel(unique(round(data{idxSubplot})));
100+
101+
subplotList(iRoi, iData) = idxSubplot;
102+
103+
idxSubplot = idxSubplot + 1;
104+
105+
end
106+
107+
end
108+
109+
nbBins = min(nbBins);
110+
111+
for i = 1:numel(data)
112+
tmp = hist(data{i}, nbBins);
113+
maxVox(end + 1) = max(tmp);
114+
end
115+
116+
%% plot histogram and mode
117+
118+
idxSubplot = 1;
119+
120+
for iRoi = 1:nbRois
121+
122+
for iData = 1:nbData
123+
124+
subplot(rows, cols, idxSubplot);
125+
126+
hold on;
127+
128+
hist(data{idxSubplot}, nbBins);
129+
130+
plot([modes(idxSubplot) modes(idxSubplot)], ...
131+
[0 max(maxVox)], ...
132+
'--r', ...
133+
'linewidth', 1);
134+
135+
bf = bids.File(roiImages{iRoi});
136+
title(['roi: ' bf.entities.label]);
137+
138+
axis([min(MIN) max(MAX) 0 max(maxVox)]);
139+
140+
idxSubplot = idxSubplot + 1;
141+
142+
end
143+
144+
end
145+
146+
for i = 1:cols
147+
subplot(rows, cols, subplotList(end, i));
148+
xlabel('intensities');
149+
end
150+
151+
for i = 1:rows
152+
subplot(rows, cols, subplotList(i, 1));
153+
ylabel('nb voxels');
154+
end
155+
156+
end

‎src/utils/optimizeSubplotNumber.m‎

Lines changed: 9 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,9 @@
1+
function [m, n] = optimizeSubplotNumber(mn)
2+
%
3+
% (C) Copyright 2022 CPP ROI developers
4+
5+
% Optimizes the number of subplot to have on a figure
6+
n = round(mn^0.4);
7+
m = ceil(mn / n);
8+
9+
end

‎tests/test_getPeakCoordinates.m‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -22,6 +22,8 @@ function test_getPeakCoordinates_basic()
2222
assertEqual(voxelCoord, [28 8 24]);
2323
assertElementsAlmostEqual(maxVal, 1.6212, 'absolute', 1e-3);
2424

25+
delete('*hemi-L_space-MNI_atlas-wang_label-V1v_mask.*');
26+
2527
end
2628

2729
function value = thisDir()

‎tests/test_plotRoi.m‎

Lines changed: 92 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,92 @@
1+
function test_suite = test_plotRoi %#ok<*STOUT>
2+
% (C) Copyright 2022 CPP ROI developers
3+
try % assignment of 'localfunctions' is necessary in Matlab >= 2016
4+
test_functions = localfunctions(); %#ok<*NASGU>
5+
catch % no problem; early Matlab versions can use initTestSuite fine
6+
end
7+
initTestSuite;
8+
end
9+
10+
function test_test_plotDataInRoi_basic()
11+
12+
close all;
13+
14+
mask1 = createDummyMask(1);
15+
data1 = createDummyData(1);
16+
17+
plotDataInRoi(data1, mask1);
18+
19+
delete *.nii;
20+
21+
end
22+
23+
function test_test_plotDataInRoi_many_rois_and_data()
24+
25+
mask1 = createDummyMask(1);
26+
mask2 = createDummyMask(2);
27+
mask3 = createDummyMask(3);
28+
mask4 = createDummyMask(4);
29+
data1 = createDummyData(1);
30+
data2 = createDummyData(2);
31+
data3 = createDummyData(3);
32+
data4 = createDummyData(4);
33+
34+
mask = cellstr(cat(1, mask1, mask2, mask3, mask4));
35+
data = cellstr(cat(1, data1, data2, data3, data4));
36+
37+
plotDataInRoi(data, mask);
38+
39+
delete *.nii;
40+
41+
end
42+
43+
function M = mat()
44+
M = ...
45+
[-3 0 0 84; ...
46+
0 3 0 -116; ...
47+
0 0 3 -56; ...
48+
0 0 0 1];
49+
50+
end
51+
52+
function mask = createDummyMask(idx)
53+
54+
MIN = randi([1, 40 / 2], 1);
55+
MAX = randi([MIN + 5, 40], 1);
56+
57+
maskHdr = struct( ...
58+
'fname', ['label-' num2str(idx) '_mask.nii'], ...
59+
'dim', [40, 40, 40], ...
60+
'dt', [spm_type('uint8') spm_platform('bigend')], ...
61+
'mat', mat(), ...
62+
'pinfo', [1 0 0]', ...
63+
'descrip', 'mask');
64+
65+
maskVol = false(40, 40, 40);
66+
maskVol(MIN:MAX, MIN:MAX, MIN:MAX) = true;
67+
68+
spm_write_vol(maskHdr, maskVol);
69+
70+
mask = fullfile(pwd, maskHdr.fname);
71+
72+
end
73+
74+
function data = createDummyData(idx)
75+
76+
dataHdr = struct( ...
77+
'fname', ['label-' num2str(idx) '_beta.nii'], ...
78+
'dim', [40, 40, 40], ...
79+
'dt', [spm_type('float32') spm_platform('bigend')], ...
80+
'mat', mat(), ...
81+
'pinfo', [1 0 0]', ...
82+
'descrip', 'data');
83+
84+
dataVol = nan(40, 40, 40);
85+
dataVol(2:39, 2:39, 2:39) = zeros(38, 38, 38);
86+
dataVol(3:38, 3:38, 3:38) = randn() * 10 + randn(36, 36, 36) * 10;
87+
88+
spm_write_vol(dataHdr, dataVol);
89+
90+
data = fullfile(pwd, dataHdr.fname);
91+
92+
end

0 commit comments

Comments
 (0)