cosmo sphere offsets

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

    % how big a box should be to contain all voxels with the given radius
    side = ceil(radius) * 2 + 1;

    % make an array that is sufficienly large (actually too big)
    offsets = zeros(side^3, 3);

    % keep track of where we store indices in offsets:
    % for every location within the sphere (in the code below),
    % add one to row_pos and then store the locations
    row_pos = 0;

    % the grid positions (relative to the origin) to consider
    single_dimension_candidates = floor(-radius):ceil(radius);

    % Consider the candidates in all three spatial dimensions using
    % nested for loops and an if statement.
    % For each position:
    % - see if it is at most at distance 'radius' from the origin
    % - if that is the case, increase row_pos and store the position in
    %   'offsets'
    %
    for x = single_dimension_candidates
        for y = single_dimension_candidates
            for z = single_dimension_candidates
                if x^2 + y^2 + z^2 <= radius^2
                    row_pos = row_pos + 1;
                    offsets(row_pos, :) = [x y z];
                end
            end
        end
    end

    % cut off empty values at the end
    offsets = offsets(1:row_pos, :);

    % compute distances
    unsorted_distances = sqrt(sum(offsets.^2, 2));

    % sort distances and apply to offsets
    [distances, idxs] = sort(unsorted_distances);
    offsets = offsets(idxs, :);