Skip to content

Commit bfbf490

Browse files
Remi-GauRemi-Gau
committed
allow to threshold image with with cluster size
1 parent a813876 commit bfbf490

File tree

7 files changed

+223
-41
lines changed

7 files changed

+223
-41
lines changed

‎src/atlas/labelClusters.m‎

Lines changed: 39 additions & 31 deletions
Original file line numberDiff line numberDiff line change
@@ -1,25 +1,59 @@
1-
function outputImage = labelClusters(sourceImage, peakThreshold, extendThreshold)
1+
function outputImage = labelClusters(varargin)
2+
%
3+
% Returns a binary mask for an image after applying voxel wise threshold and
4+
% an optional cluster size threshold
5+
%
6+
% USAGE::
7+
%
8+
% outputImage = thresholdToMask(inputImage, peakThreshold, clusterSize)
9+
%
10+
% :param inputImage:
11+
% :type inputImage: path
12+
% :param peakThreshold:
13+
% :type peakThreshold: float
14+
% :param clusterSize:
15+
% :type clusterSize: integer >= 0 (Default: 0)
16+
%
17+
% :returns: - :outputImage: (string)
18+
%
19+
% See also getClusters, sortAndLabelClusters
220
%
321
% Adapted from:
422
% https://en.wikibooks.org/wiki/SPM/How-to#How_to_remove_clusters_under_a_certain_size_in_a_binary_mask?
523
%
624
% (C) Copyright 2021 CPP ROI developers
725

8-
hdr = spm_vol(sourceImage);
26+
default_clusterSize = 0;
27+
28+
isFile = @(x) exists(x, 'file') == 2;
29+
isPositive = @(x) isinteger(x) && X >= 0;
30+
31+
p = inputParser;
32+
33+
addRequired(p, 'inputImage', isFile);
34+
addRequired(p, 'peakThreshold', @isnumeric);
35+
addOptional(p, 'clusterSize', default_clusterSize, isPositive);
36+
37+
parse(p, varargin{:});
938

10-
[l2, num] = getClusters(hdr, peakThreshold);
39+
inputImage = p.Results.inputImage;
40+
peakThreshold = p.Results.peakThreshold;
41+
clusterSize = p.Results.clusterSize;
42+
43+
[l2, num] = getClusters(inputImage, peakThreshold);
1144

1245
if ~num
1346
warning('No clusters found.');
1447
return
1548
end
1649

17-
vol = sortAndLabelClusters(l2, num, extendThreshold);
50+
vol = sortAndLabelClusters(l2, num, clusterSize);
1851

1952
% Write new image with cluster laebelled with their voxel size
20-
p = bids.internal.parse_filename(sourceImage);
53+
p = bids.internal.parse_filename(inputImage);
2154
p.suffix = 'dseg'; % discrete segmentation
2255

56+
hdr = spm_vol(inputImage);
2357
bidsFile = bids.File(p);
2458
hdr.fname = spm_file(hdr.fname, 'filename', bidsFile.filename);
2559

@@ -31,29 +65,3 @@
3165
bids.util.jsonencode(json.filename, json.content);
3266

3367
end
34-
35-
function [l2, num] = getClusters(hdr, peakThreshold)
36-
data = spm_read_vols(hdr);
37-
[l2, num] = spm_bwlabel(double(data > peakThreshold), 26);
38-
end
39-
40-
function vol = sortAndLabelClusters(l2, num, extendThreshold)
41-
42-
% Extent threshold, and sort clusters by their extent
43-
% Label corresponds to their rank
44-
[n, ni] = sort(histc(l2(:), 0:num), 1, 'descend');
45-
vol = zeros(size(l2));
46-
n = n(2:end);
47-
ni = ni(2:end) - 1;
48-
ni = ni(n >= extendThreshold);
49-
n = n(n >= extendThreshold);
50-
for i = 1:length(n)
51-
vol(l2 == ni(i)) = i;
52-
end
53-
54-
fprintf('Selected %d clusters (out of %d) in image.\n', length(n), num);
55-
for i = 1:length(n)
56-
fprintf('Cluster label %i ; size: %i voxels\n', i, n(i));
57-
end
58-
59-
end

‎src/roi/thresholdToMask.m‎

Lines changed: 40 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -1,26 +1,56 @@
1-
function outputImage = thresholdToMask(inputImage, threshold)
1+
function outputImage = thresholdToMask(varargin)
22
%
3-
% TODO
4-
% allow threshold to be inferior than, greater than or both
3+
% Returns a binary mask for an image after applying voxel wise threshold and
4+
% an optional cluster size threshold
5+
%
6+
% USAGE::
57
%
8+
% outputImage = thresholdToMask(inputImage, peakThreshold, clusterSize)
69
%
10+
% :param inputImage:
11+
% :type inputImage: path
12+
% :param peakThreshold:
13+
% :type peakThreshold: float
14+
% :param clusterSize:
15+
% :type clusterSize: integer >= 0 (Default)
16+
%
17+
% :returns: - :outputImage: (string)
718
%
819
% (C) Copyright 2021 CPP ROI developers
920

10-
hdr = spm_vol(inputImage);
11-
img = spm_read_vols(hdr);
21+
% TODO
22+
% allow threshold to be inferior than, greater than or both
23+
24+
default_clusterSize = 0;
25+
26+
isFile = @(x) exist(x, 'file') == 2;
27+
isPositive = @(x) isnumeric(x) && x >= 0;
28+
29+
p = inputParser;
1230

31+
addRequired(p, 'inputImage', isFile);
32+
addRequired(p, 'peakThreshold', @isnumeric);
33+
addOptional(p, 'clusterSize', default_clusterSize, isPositive);
34+
35+
parse(p, varargin{:});
36+
37+
inputImage = p.Results.inputImage;
38+
peakThreshold = p.Results.peakThreshold;
39+
clusterSize = p.Results.clusterSize;
40+
41+
% create output name
1342
p = bids.internal.parse_filename(inputImage);
1443
p.suffix = 'mask';
15-
1644
bidsFile = bids.File(p);
17-
45+
hdr = spm_vol(inputImage);
1846
hdr.fname = spm_file(hdr.fname, 'filename', bidsFile.filename);
19-
img = img > threshold;
20-
spm_write_vol(hdr, img);
21-
2247
outputImage = hdr.fname;
2348

49+
[l2, num] = getClusters(inputImage, peakThreshold);
50+
vol = sortAndThresholdClusters(l2, num, clusterSize);
51+
52+
spm_write_vol(hdr, vol);
53+
2454
json = bids.derivatives_json(outputImage);
2555
bids.util.jsonencode(json.filename, json.content);
2656

‎src/utils/getClusters.m‎

Lines changed: 21 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,21 @@
1+
function [l2, num] = getClusters(inputImage, peakThreshold)
2+
%
3+
% USAGE::
4+
%
5+
% [l2, num] = getClusters(inputImage, peakThreshold)
6+
%
7+
% :param sourceImage: string
8+
% :type sourceImage: string
9+
% :param peakThreshold: string
10+
% :type peakThreshold: string
11+
%
12+
% See also labelClusters, sortAndLabelClusters
13+
%
14+
% (C) Copyright 2021 CPP ROI developers
15+
16+
hdr = spm_vol(inputImage);
17+
vol = spm_read_vols(hdr);
18+
19+
[l2, num] = spm_bwlabel(double(vol > peakThreshold), 26);
20+
21+
end

‎src/utils/rootDir.m‎

Lines changed: 12 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,12 @@
1+
function rootDir = rootDir()
2+
%
3+
% USAGE::
4+
%
5+
% rootDir = rootDir()
6+
%
7+
%
8+
% (C) Copyright 2021 CPP ROI developers
9+
10+
rootDir = fullfile(fileparts(mfilename('fullpath')), '..', '..');
11+
12+
end

‎src/utils/sortAndLabelClusters.m‎

Lines changed: 34 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,34 @@
1+
function vol = sortAndLabelClusters(l2, num, clusterSize)
2+
%
3+
% Creates a binary image after:
4+
% - Sorting clusters by their extent
5+
% - Removes cluster with smaller than clusterSize
6+
% - Labelling cluster with their rank in the sorting
7+
%
8+
% USAGE::
9+
%
10+
% vol = sortAndLabelClusters(l2, num, clusterSize)
11+
%
12+
%
13+
% See also: getClusters, labelClusters
14+
%
15+
% (C) Copyright 2021 CPP ROI developers
16+
17+
% refactor with sortAndThresholdClusters
18+
19+
[n, ni] = sort(histc(l2(:), 0:num), 1, 'descend');
20+
vol = zeros(size(l2));
21+
n = n(2:end);
22+
ni = ni(2:end) - 1;
23+
ni = ni(n >= clusterSize);
24+
n = n(n >= clusterSize);
25+
for i = 1:length(n)
26+
vol(l2 == ni(i)) = i;
27+
end
28+
29+
fprintf('Selected %d clusters (out of %d) in image.\n', length(n), num);
30+
for i = 1:length(n)
31+
fprintf('Cluster label %i ; size: %i voxels\n', i, n(i));
32+
end
33+
34+
end

‎src/utils/sortAndThresholdClusters.m‎

Lines changed: 15 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,15 @@
1+
function vol = sortAndThresholdClusters(l2, num, clusterSize)
2+
%
3+
% (C) Copyright 2020 CPP ROI developers
4+
5+
[n, ni] = sort(histc(l2(:), 0:num), 1, 'descend');
6+
vol = zeros(size(l2));
7+
n = n(2:end);
8+
ni = ni(2:end) - 1;
9+
ni = ni(n >= clusterSize);
10+
n = n(n >= clusterSize);
11+
for i = 1:length(n)
12+
vol(l2 == ni(i)) = 1;
13+
end
14+
15+
end

‎tests/test_thresholdToMask.m‎

Lines changed: 62 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,62 @@
1+
function test_suite = test_thresholdToMask() %#ok<*STOUT>
2+
%
3+
% (C) Copyright 2021 CPP ROI developers
4+
5+
try % assignment of 'localfunctions' is necessary in Matlab >= 2016
6+
test_functions = localfunctions(); %#ok<*NASGU>
7+
catch % no problem; early Matlab versions can use initTestSuite fine
8+
end
9+
initTestSuite;
10+
end
11+
12+
% TODO refactor set up and teardown
13+
14+
function test_thresholdToMask_default()
15+
16+
image = 'visual motion_association-test_z_FDR_0.01.nii.gz';
17+
18+
inputImage = setUp(image);
19+
20+
peakThreshold = 5;
21+
outputImage = thresholdToMask(inputImage, peakThreshold);
22+
23+
% check that we have certain number voxels in the mask
24+
vol = spm_read_vols(spm_vol(outputImage));
25+
assertEqual(sum(vol(:) > 0), 2008);
26+
27+
teardown();
28+
29+
end
30+
31+
function test_thresholdToMask_cluster()
32+
33+
image = 'visual motion_association-test_z_FDR_0.01.nii.gz';
34+
35+
inputImage = setUp(image);
36+
37+
peakThreshold = 5;
38+
clusterSize = 10;
39+
outputImage = thresholdToMask(inputImage, peakThreshold, clusterSize);
40+
41+
% check that we have certain number voxels in the mask
42+
vol = spm_read_vols(spm_vol(outputImage));
43+
assertEqual(sum(vol(:) > 0), 1926);
44+
45+
teardown();
46+
47+
end
48+
49+
function inputImage = setUp(image)
50+
51+
gzImage = spm_file(fullfile(rootDir(), 'demos', 'roi', 'inputs', image), 'cpath');
52+
[~, basename] = spm_fileparts(gzImage);
53+
gunzip(gzImage, pwd);
54+
55+
inputImage = renameNeuroSynth(fullfile(pwd, basename));
56+
57+
end
58+
59+
function teardown()
60+
delete('*.nii');
61+
delete('*.json');
62+
end

0 commit comments

Comments
 (0)