Getting the centre of mass for an activation

Occasionally it is useful to get a blob's centre of mass, for instance in a binary mask where peaks are not available. This function presents a simplistic solution.

% Returns the centre of mass for blob(s) in a volume. Input either
% handle from spm_vol, or path to volume. If you aren't happy with
% the way the function divides the volume into blobs, you can mask
% off the blobs you like manually and save each blob to its own
% volume.
% Use: CoM = getCentreOfMass(vol)
% 27/10/2009 J Carlin, with thanks to Danny Mitchell and Russell
% Thompson

function CoM = getCentreOfMass(vol)

if isstr(vol)
    vol = spm_vol(vol);
end

v = spm_read_vols(vol);

v(isnan(v)) = 0;

% Identify and number code blobs
blobs = bwlabeln(v);

for i = 1:max(blobs(:))
    % Blob indices
    [x,y,z] = ind2sub(size(blobs),find(blobs == i));
    XYZ = [x, y, z];
    CoM(i,:) = mean(XYZ,1);
end

(JohanCarlin)