function sandMats = formsandmatrices( resShifterMats, accessibleDistanceMat )
    % Compute the 'sand' matrices at each resolution, from the lowest (2 by
    % 2 grid) up to the highest.
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Distances are always measured from positions at highest resolution.
    % As we go down resolutions, the distance is measured to the nearest
    % high-resolution position within the 'chunked-up' regions.
    %
    % sandMats - cell array. Resolution indexes cells. Each cell is a
    % matrix with D rows and d columns, where D is the number of valid
    % positions at the highest resolution, and d is the number of valid
    % positions at the current resolution.
    %
    % resShifterMats - a cell array, with cells indexed by resolution of 2D
    % matrices that convert discrete positions at the previous resolution
    % (coarser, indexing rows) to discrete positions at the current
    % resolution (i.e., res, from 1 to maximum). Thus there are 4 columns
    % for each such matrix, and the first cell has only one row (zeroth
    % resolution is just one discretisation cell). Only valid positions are
    % listed in the columns. If there are less than 4, the remainder is
    % padded with NaNs.
    %
    % accessibleDistanceMat - the distance matrix for the highest
    % resolution. This is a sparse upper triangular matrix.
    %
    % There are 2^res discretisation grid cells covering the longest
    % spatial dimension; equivalently, the longest spatial dimension has
    % been divided into parts of size 2^(-res) times the whole length.
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % CHANGELOG
    % 31.10.2014 - copied from old version.
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    maxRes = length(resShifterMats);
    sandMats = cell([1, maxRes]);

    % Number of valid positions at highest resolution.
    nPositions = size(accessibleDistanceMat, 2);
    
%     accessibleDistanceMat = full(accessibleDistanceMat + accessibleDistanceMat');

    % Start at highest resolution, then work down.
    maxDist = max(max(accessibleDistanceMat));
    sandMats{maxRes} = log(1 + maxDist - accessibleDistanceMat);
    
    for res = (maxRes - 1):(-1):1,
        nPosThisRes = size(resShifterMats{res + 1}, 1);
        sandMats{res} = zeros([nPositions, nPosThisRes]);
        for iPos = 1:nPosThisRes,
            % We want to compute the minimum distance from every position
            % to a position somewhere in the regions at this resolution.
            % First get the highest-resolution positions corresponding to
            % this discrete position.
            highResPos = iPos;
            for iRes = (res + 1):maxRes,
                highResPos = nonzeros(resShifterMats{iRes}(highResPos, :));
            end
            sandMats{res}(:, iPos) = min(accessibleDistanceMat(:, highResPos), [], 2);
        end
        maxMinDist = max(max(sandMats{res}));
        sandMats{res} = log(1 + maxMinDist - sandMats{res});
    end
end