Skip to content

Commit 6b5a42c

Browse files
Remi-GauRemi-Gau
committed
trasnfer source code from CPP_SPM
1 parent ee398fb commit 6b5a42c

File tree

5 files changed

+343
-0
lines changed

5 files changed

+343
-0
lines changed

‎src/roi/checkRoiOrientation.m‎

Lines changed: 17 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,17 @@
1+
% (C) Copyright 2021 CPP BIDS SPM-pipeline developers
2+
3+
function [sts, images] = checkRoiOrientation(referenceImage, imagesToCheck)
4+
5+
% referenceImage - better if fullfile path
6+
% imagesToCheck - better if fullfile path
7+
8+
% TODO
9+
% - make it possible to pass several images in imagesToCheck at one
10+
11+
images = char({referenceImage; imagesToCheck});
12+
13+
% check if files are in the same space
14+
% if not we reslice the ROI to have the same resolution as the data image
15+
sts = spm_check_orientations(spm_vol(images));
16+
17+
end

‎src/roi/createRoi.m‎

Lines changed: 207 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,207 @@
1+
% (C) Copyright 2021 CPP BIDS SPM-pipeline developers
2+
3+
function mask = createRoi(type, varargin)
4+
%
5+
% Returns a mask to be used as a ROI by ``spm_summarize``.
6+
% Can also save the ROI as binary image.
7+
%
8+
% USAGE::
9+
%
10+
% mask = createROI(type, varargin);
11+
%
12+
% mask = createROI('mask', roiImage, volumeDefiningImage, saveImg = false);
13+
% mask = createROI('sphere', sphere, volumeDefiningImage, saveImg = false);
14+
% mask = createROI('intersection', roiImage, sphere, volumeDefiningImage, saveImg = false);
15+
% mask = createROI('expand', roiImage, sphere, volumeDefiningImage, saveImg = false);
16+
%
17+
% See the ``demos/roi`` to see examples on how to use it.
18+
%
19+
% :param type: ``'mask'``, ``'sphere'``, ``'intersection'``, ``'expand'``
20+
% :type type: string
21+
% :param roiImage: fullpath of the roi image
22+
% :type roiImage: string
23+
% :param sphere:
24+
% ``sphere.location``: X Y Z coordinates in millimeters
25+
% ``sphere.radius``: radius in millimeters
26+
% :type sphere: structure
27+
% :param volumeDefiningImage: fullpath of the image that will define the space
28+
% (resolution, ...) if the ROI is to be saved.
29+
% :type volumeDefiningImage: string
30+
% :param saveImg: Will save the resulting image as binary mask if set to
31+
% ``true``
32+
% :type saveImg: boolean
33+
%
34+
% :returns:
35+
%
36+
% mask - structure for the volume of interest adapted from ``spm_ROI``
37+
%
38+
% mask.def - VOI definition [sphere, mask]
39+
% mask.rej - cell array of disabled VOI definition options
40+
% mask.xyz - centre of VOI {mm} (for sphere)
41+
% mask.spec - VOI definition parameters (radius for sphere)
42+
% mask.str - description of the VOI
43+
% mask.roi.size - number of voxel in ROI
44+
% mask.roi.XYZ - voxel coordinates
45+
% mask.roi.XYZmm - voxel world coordinates
46+
% mask.global.hdr - header of the "search space" where the roi is
47+
% defined
48+
% mask.global.img
49+
% mask.global.XYZ
50+
% mask.global.XYZmm
51+
%
52+
%
53+
54+
if islogical(varargin{end})
55+
saveImg = varargin{end};
56+
else
57+
saveImg = false;
58+
end
59+
60+
if saveImg
61+
switch type
62+
case {'sphere', 'mask'}
63+
volumeDefiningImage = varargin{2};
64+
case {'intersection', 'expand'}
65+
volumeDefiningImage = varargin{3};
66+
end
67+
end
68+
69+
switch type
70+
71+
case 'sphere'
72+
73+
specification = varargin{1};
74+
75+
mask.def = type;
76+
mask.spec = specification.radius;
77+
mask.xyz = specification.location;
78+
79+
if size(mask.xyz, 1) ~= 3
80+
mask.xyz = mask.xyz';
81+
end
82+
83+
mask = spm_ROI(mask);
84+
mask.roi.XYZmm = [];
85+
86+
case 'mask'
87+
88+
roiImage = varargin{1};
89+
90+
mask = struct('XYZmm', []);
91+
mask = defineGlobalSearchSpace(mask, roiImage);
92+
93+
% in real world coordinates
94+
mask.global.XYZmm = returnXYZm(mask.global.hdr.mat, mask.global.XYZ);
95+
96+
assert(size(mask.global.XYZmm, 2) == sum(mask.global.img(:)));
97+
98+
locationsToSample = mask.global.XYZmm;
99+
100+
mask.def = type;
101+
mask.spec = roiImage;
102+
[~, mask.roi.XYZmm, j] = spm_ROI(mask, locationsToSample);
103+
104+
mask.roi.XYZ = mask.global.XYZ(:, j);
105+
106+
mask = setRoiSizeAndType(mask, type);
107+
108+
case 'intersection'
109+
110+
roiImage = varargin{1};
111+
mask = createRoi('mask', roiImage);
112+
113+
specification = varargin{2};
114+
mask2 = createRoi('sphere', specification);
115+
116+
locationsToSample = mask.global.XYZmm;
117+
118+
[~, mask.roi.XYZmm] = spm_ROI(mask2, locationsToSample);
119+
120+
mask = setRoiSizeAndType(mask, type);
121+
122+
case 'expand'
123+
124+
roiImage = varargin{1};
125+
specification = varargin{2};
126+
127+
% take as radius step the smallest voxel dimension of the roi image
128+
hdr = spm_vol(roiImage);
129+
dim = diag(hdr.mat);
130+
radiusStep = min(abs(dim(1:3)));
131+
132+
while true
133+
mask = createRoi('intersection', roiImage, specification);
134+
mask.roi.radius = specification.radius;
135+
if mask.roi.size > specification.maxNbVoxels
136+
break
137+
end
138+
specification.radius = specification.radius + radiusStep;
139+
end
140+
141+
mask = setRoiSizeAndType(mask, type);
142+
143+
end
144+
145+
if saveImg
146+
saveRoi(mask, volumeDefiningImage);
147+
end
148+
149+
end
150+
151+
function mask = defineGlobalSearchSpace(mask, image)
152+
153+
mask.global.hdr = spm_vol(image);
154+
mask.global.img = logical(spm_read_vols(mask.global.hdr));
155+
156+
[X, Y, Z] = ind2sub(size(mask.global.img), find(mask.global.img));
157+
158+
% XYZ format
159+
mask.global.XYZ = [X'; Y'; Z'];
160+
mask.global.size = size(mask.global.XYZ, 2);
161+
162+
end
163+
164+
function XYZmm = returnXYZm(transformationMatrix, XYZ)
165+
% "voxel to world transformation"
166+
XYZmm = transformationMatrix(1:3, :) * [XYZ; ones(1, size(XYZ, 2))];
167+
end
168+
169+
function saveRoi(mask, volumeDefiningImage)
170+
171+
checkDependencies('marsbar');
172+
173+
switch mask.def
174+
175+
case 'sphere'
176+
177+
[~, mask.roi.XYZmm] = spm_ROI(mask, volumeDefiningImage);
178+
mask = setRoiSizeAndType(mask, mask.def);
179+
180+
radius = mask.spec;
181+
center = mask.xyz;
182+
183+
descrip = sprintf('%0.1fmm radius sphere at [%0.1f %0.1f %0.1f]', radius, center);
184+
label = sprintf('sphere_%0.0f-%0.0f_%0.0f_%0.0f', radius, center);
185+
186+
otherwise
187+
188+
label = [mask.def '-roiMcRoiFace'];
189+
descrip = [mask.def '-roiMcRoiFace'];
190+
191+
end
192+
193+
% use the marsbar toolbox
194+
roiObject = maroi_pointlist(struct('XYZ', mask.roi.XYZmm, ...
195+
'mat', spm_get_space(volumeDefiningImage), ...
196+
'label', label, ...
197+
'descrip', descrip));
198+
199+
roiName = sprintf('%s.nii', label);
200+
save_as_image(roiObject, fullfile(pwd, roiName));
201+
202+
end
203+
204+
function mask = setRoiSizeAndType(mask, type)
205+
mask.def = type;
206+
mask.roi.size = size(mask.roi.XYZmm, 2);
207+
end

‎src/roi/getRetinoProbaAtlas.m‎

Lines changed: 68 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,68 @@
1+
% (C) Copyright 2021 CPP BIDS SPM-pipeline developers
2+
3+
function [maxProbaFiles, roiLabels] = getRetinoProbaAtlas()
4+
%
5+
% Loads the volumetric data from the
6+
% Probabilistic Maps of Visual Topography in Human Cortex
7+
%
8+
% maxProbaVol: 4D volume of ROIs in MNII space with
9+
%
10+
% - left hemisphere as maxProbaVol(:,:,:,1)
11+
% - right hemisphere as maxProbaVol(:,:,:,2)
12+
%
13+
% DOI 10.1093/cercor/bhu277
14+
% PMCID: PMC4585523
15+
% PMID: 25452571
16+
% Probabilistic Maps of Visual Topography in Human Cortex
17+
%
18+
% If the data is not present in the ``cpp_spm/atlas/ProbAtlas_v4`` of the repo, it will be
19+
% downloaded and unzipped
20+
21+
README_URL = 'http://scholar.princeton.edu/sites/default/files/napl/files/readme.txt';
22+
ATLAS_URL = 'http://scholar.princeton.edu/sites/default/files/napl/files/probatlas_v4.zip';
23+
24+
ATLAS_FOLDER = fullfile( ...
25+
fileparts(mfilename('fullpath')), ...
26+
'..', '..', 'atlas');
27+
28+
if ~exist(fullfile(ATLAS_FOLDER, 'ProbAtlas_v4'), 'dir')
29+
30+
mkdir(ATLAS_FOLDER);
31+
32+
try
33+
urlwrite(ATLAS_URL, 'probatlas_v4.zip');
34+
unzip('probatlas_v4.zip', ATLAS_FOLDER);
35+
catch
36+
system(sprintf('curl -L %s -o probatlas_v4.zip', ATLAS_URL));
37+
unzip('probatlas_v4.zip', ATLAS_FOLDER);
38+
end
39+
40+
end
41+
42+
% urlwrite(README_URL, fullfile(ATLAS_FOLDER, 'ProbAtlas_v4', 'README.txt'));
43+
% urlread(README_URL)
44+
45+
maxProbaFiles = spm_select('FPListRec', ...
46+
fullfile(ATLAS_FOLDER, 'ProbAtlas_v4'), ...
47+
'^max.*.nii$');
48+
49+
if size(maxProbaFiles, 1) < 2
50+
51+
gunzip(fullfile(pwd, 'atlas', 'ProbAtlas_v4', 'subj_vol_all', 'max*.nii.gz'));
52+
53+
maxProbaFiles = spm_select('FPListRec', ...
54+
fullfile(ATLAS_FOLDER, 'ProbAtlas_v4'), ...
55+
'^max.*.nii$');
56+
57+
end
58+
59+
if size(maxProbaFiles, 1) < 2
60+
error('no atlas present');
61+
end
62+
63+
% the label file was created manually to be easier to load into SPM
64+
roiLabels = spm_load(spm_select('FPListRec', ...
65+
fullfile(ATLAS_FOLDER), ...
66+
'^ROI_labels_ProbAtlas_v4.csv$'));
67+
68+
end

‎src/roi/resliceRoiImages.m‎

Lines changed: 30 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,30 @@
1+
% (C) Copyright 2021 CPP BIDS SPM-pipeline developers
2+
3+
function reslicedImages = resliceRoiImages(referenceImage, imagesToCheck)
4+
5+
% TODO
6+
% - make prefix more flexible
7+
% - make it possible to pass several images in imagesToCheck at one
8+
% - allow option to binarize output?
9+
10+
% check if files are in the same space
11+
% if not we reslice the ROI to have the same resolution as the data image
12+
[sts, images] = checkRoiOrientation(referenceImage, imagesToCheck);
13+
if sts == 1
14+
reslicedImages = imagesToCheck;
15+
16+
else
17+
18+
flags = struct( ...
19+
'mean', false, ...
20+
'which', 1, ...
21+
'prefix', 'r');
22+
23+
spm_reslice(images, flags);
24+
25+
basename = spm_file(imagesToCheck, 'basename');
26+
reslicedImages = spm_file(imagesToCheck, 'basename', [flags.prefix basename]);
27+
28+
end
29+
30+
end

‎src/roi/thresholdToMask.m‎

Lines changed: 21 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,21 @@
1+
% (C) Copyright 2021 CPP BIDS SPM-pipeline developers
2+
3+
function roiName = thresholdToMask(sourceImage, threshold)
4+
5+
% TODO
6+
% allow threshold to be inferior than, greater than or both
7+
8+
hdr = spm_vol(sourceImage);
9+
img = spm_read_vols(hdr);
10+
11+
roi_hdr = hdr;
12+
basename = spm_file(roi_hdr.fname, 'basename');
13+
roi_hdr.fname = spm_file(roi_hdr.fname, 'basename', [basename '-mask']);
14+
15+
img = img > threshold;
16+
17+
spm_write_vol(roi_hdr, img);
18+
19+
roiName = roi_hdr.fname;
20+
21+
end

0 commit comments

Comments
 (0)