Hello everyone,
I am trying to extract peak coordinates from given clusters in an 3D fMRI group image, e.g. a thresholded ICA z-score map.
Using MATLAB/SPM I am able to extract a single peak coordinate for each cluster by numbering the clusters with spm_bwlabel(), listing the voxel values per cluster and searching for the maximum value. However, I'd like to extract more than only one peak that might exist within each cluster (e.g. the SPM results viewer displays up to 3 peaks by default).
Did I miss an obvious solution or an SPM function, respectively? I could imagine to, within a main cluster, search for "subclusters" and again compute their peak values - but there might be an easier solution?
Any help is appreciated! :)
Best regards, Leon
The function you might be looking for is spm_max().
A minimal working example:
% load image
vol = spm_vol('some file.nii');
[Y,XYZ] = spm_read_vols(vol);
% voxel space
dim = size(Y);
[R,C,P] = ndgrid(1:dim(1),1:dim(2),1:dim(3));
RCP = [R(:)';C(:)';P(:)'];
% thresholding
ind = Y > 4;
Y = Y(ind);
RCP = RCP(:,ind);
XYZ = XYZ(:,ind);
% get maxima
[N,Z,M,A,~] = spm_max(Y(:),RCP);
All information regarding clusters and local maxima are now stored in the variabels N, Z, M and A. (See "help spm_max".) You can now easily print results tables with arbitrary numbers of peaks per cluster.
Minimal working example:
% create results table
resTbl = {'cluster-level','voxel-level','';'k','val','x,y,z (mm)'};
curRow = 3;
for iClust = 1:max(A)
% get associated voxels
idxClust = A == iClust;
% sort by value
[B,I] = sort(Z(idxClust),'descend');
% cluster extent
n = N(find(idxClust,1));
resTbl{curRow,1} = n;
% iterate over peaks
for iPeak = 1:min(length(I),Inf)
% value
resTbl{curRow,2} = B(iPeak);
% coordinate
tmp = XYZ(1:3,idxClust);
xyz = tmp(:,I(iPeak));
xyzStr = mat2str(xyz');
resTbl{curRow,3} = xyzStr(2:end-1);
curRow = curRow + 1;
end
end
A table with cluster sizes (k), peak values and coordinates is now stored in resTbl.
I think I encountered a small error in your code. The row "tmp = XYZ(1:3,idxClust);" will index XYZ which has as many columns as there are over-threshold-voxels with idxClust which is the index of voxels within the current cluster. Thus the returned coordinates will be the wrong ones. I modified it and my final working function looks like this:
Thanks!
function [resTbl] = get_cluster_maxima(map, voxthresh, sizethresh, npeaks)
%function [resTbl] = get_cluster_maxima(map, voxthresh, sizethresh, npeaks)
% Returns table with up to npeaks maxima for clusters retrieved from map
% and defined by voxthresh and sizethresh.
%e.g.:
%map = 'my_file.nii';
%voxthresh = 2.5;
%sizethresh = 10;
%npeaks = 3;
disp(map);
% load image
vol = spm_vol(map);
[Y, XYZ] = spm_read_vols(vol);
% voxel space
dim = size(Y);
[R,C,P] = ndgrid(1:dim(1), 1:dim(2), 1:dim(3));
RCP = [R(:)';C(:)';P(:)'];
% thresholding
ind = Y > voxthresh;
Y = Y(ind);
RCP = RCP(:, ind);
XYZ = XYZ(:, ind);
% get maxima
[N,Z,M,A,~] = spm_max(Y(:), RCP);
% get MNI coordinates for all clusters
[~, idxCoord] = ismember(M.', RCP.', 'rows');
XYZpeak = XYZ(1:3, idxCoord);
% create results table
resTbl = {'cluster-level','voxel-level','coordinates (mm)',[],[]; ...
'k','val','x','y','z'};
curRow = 3;
for iClust = 1:max(A)
% get associated voxels
idxClust = A == iClust;
% get associated coordinates
xyzClust = XYZpeak(:,idxClust);
% sort by value
[B,I] = sort(Z(idxClust),'descend');
% cluster extent
n = N(find(idxClust,1));
% cluster size threshold for printing
if n > sizethresh
resTbl{curRow,1} = n;
% iterate over peaks
for iPeak = 1:min(length(I), npeaks)
% value
resTbl{curRow,2} = B(iPeak);
% coordinates
xyz = xyzClust(:, I(iPeak));
resTbl{curRow,3} = xyz(1);
resTbl{curRow,4} = xyz(2);
resTbl{curRow,5} = xyz(3);
curRow = curRow + 1;
end
end
end
% display table
disp(resTbl);
% save csv
[p,n,~] = fileparts(map);
writetable(table(resTbl), fullfile(p, [n '_peaks.csv']), 'WriteVariableNames', 0);
end
- Badges
- Science topic
- Similar topics
- Computer Science
- Images