cosmo sphere offsets hdr

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.           #