function [offsets,distances]=cosmo_sphere_offsets(radius)
% computes sub index offsets for voxels in a sphere
%
% [offsets,distances]=cosmo_sphere_offsets(radius)
%
% Input:
% radius radius of the sphere (in voxel units)
%
% Output
% offsets Px3 sub indices relative to the origin (0,0,0).
% offsets(p,:)=[i, j, k] means that the euclidian distance
% between points at (i,j,k) and the origin is less than or
% or equal to radius
% distances Px1 distances from the origin (in voxel units).
%
% Example:
% % compute offsets for voxels that share a side or edge, but not those
% % that only share a corner (because sqrt(2) < 1.5 < sqrt(3)).
% [o, d]=cosmo_sphere_offsets(1.5);
% cosmo_disp(o);
% %|| [ 0 0 0
% %|| -1 0 0
% %|| 0 -1 0
% %|| : : :
% %|| 1 0 -1
% %|| 1 0 1
% %|| 1 1 0 ]@19x3
% cosmo_disp(d)
% %|| [ 0
% %|| 1
% %|| 1
% %|| :
% %|| 1.41
% %|| 1.41
% %|| 1.41 ]@19x1
%
% Notes:
% - this function computes distances in voxel space, not in world space.
% If voxels have the same size sof mm along all dimensions (i.e. are
% of size s x s x s mm^2 for some value of s; this property is known as
% "isotropic") then using radius=d/s will select voxels within
% distance d mm. This function is less suitable for spherical offsets
% in world space when voxels are non-isotropic.
% - offsets and distances are sorted by distance.
%
% See also sub2ind, ind2sub, cosmo_spherical_neighborhood
%
% # For CoSMoMVPA's copyright information and license terms, #
% # see the COPYING file distributed with CoSMoMVPA. #
