function resShifterMats = formresolutionshifters( validPositionsLogMat )
    % Form a cell array of matrices that convert positions of different
    % resolutions. Thus the r th cell maps positions at index r - 1,
    % indexing rows, to the (usually) 4 positions at resolution r covered
    % by the coarser discrete position. Only valid positions are mapped to.
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % validPositionsLogMat is a m by n matrix of logicals, where m, n are
    % the Y and X dimensions at the highest resolution, for which the (i,
    % j) entry is true if the corresponding position is considered 'valid',
    % or false otherwise.
    %
    % Note: a discrete position - at any resolution - is considered valid
    % if it contains a discrete position at the highest resolution that is
    % known to be valid via validPositionsLogMat.
    %
    % Zeros can appear in a row when there are fewer than 4 valid positions
    % at the higher resolution.
    %
    % The number of rows of the r th resolution matrix is equal to the
    % number of valid positions at resolution r - 1.
    %
    % The algorithm works as follows: each resolution 'inherits' a logical
    % matrix over all discrete grid cells at that resolution. We start at
    % the highest resolution; the log mat inherited first is thus the
    % validPositionsLogMat supplied as a parameter. The linear indices for
    % all valid positions at this resolution are placed at their position
    % on the grid. We aim to loop over every possible linear index at the
    % NEXT (one lower) resolution.
    % The index labels for the next resolution are computed thus: on the
    % discrete grid we list the linear indices for every possible position,
    % then we halve these, take the ceiling, double the columns, then
    % reshape. Finally we use the log mat to 'turn off' the invalid
    % positions (made zeros).
    % In our loop these labels are used to identify the regions where the
    % linear indices are for positions at the current resolution.
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % CHANGELOG
    % 31.10.2014 - copied from old version.
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    if ~all(all(islogical(validPositionsLogMat))),
        error('validPositionsLogMat must be a logical matrix!')
    end

    [yDim, xDim] = size(validPositionsLogMat);
    maxDim = max([yDim, xDim]);
    maxRes = ceil(log2(maxDim));
    resShifterMats = cell([1, maxRes]);
    
    % Initialise.
    currentValidPositionsLogMat = validPositionsLogMat;
        
    for res = maxRes:(-1):1,
        nValid = sum(sum(currentValidPositionsLogMat));

        % If there are an odd number of rows at this discretisation, our
        % linear indexing in nextResIndsMat can go wrong.
        if rem(yDim, 2) == 1,
            linearIndsMat = reshape(1:((yDim + 1) * xDim), [yDim + 1, xDim]);
            linearIndsMat = linearIndsMat(1:yDim, :);
        else
            linearIndsMat = reshape(1:(yDim * xDim), [yDim, xDim]);
        end

        validPosMat = zeros([yDim, xDim]);
        validPosMat(currentValidPositionsLogMat) = 1:nValid;

        % This matrix will have the linear indices of all positions at the
        % next (lower) resolution on the current resolution grid, i.e. they
        % will be replicated up to 4 times. Zeros will be in the place of
        % invalid positions. Note that these linear indices count all
        % positions regardless of being valid or not, so some numbers may
        % be skipped.
        nextResIndsMat = reshape(repmat(ceil(linearIndsMat / 2), [2, 1]), [yDim, xDim * 2]);
        nextResIndsMat = nextResIndsMat(:, 1:xDim);
        nextResIndsMat(~currentValidPositionsLogMat) = 0;

        % Update some of the structures for the next loop iteration.
        yDim = ceil(yDim / 2);
        xDim = ceil(xDim / 2);
        maxPosAtLowerRes = yDim * xDim;
        resShifterMats{res} = zeros([maxPosAtLowerRes, 4]);
        currentValidPositionsLogMat = false([yDim, xDim]);
        % Loop over all possible lower resolution positions. For those
        % which are valid, identify the current resolution positions with
        % them.
        validInds = unique(nonzeros(nextResIndsMat));
        % Use this linear index to correct for the indices in
        % nextResIndsMat skipping over invalid positions.
        linearInd = 0;
        for iPos = validInds',
            cellLocationLog = nextResIndsMat == iPos;
            linearInd = linearInd + 1;
            currentValidPositionsLogMat(iPos) = true;
            thisResIndsHere = validPosMat(cellLocationLog);
            resShifterMats{res}(linearInd, 1:length(thisResIndsHere)) = thisResIndsHere;
        end
        % Get rid of the all-zero rows (corresponding to invalid positions
        % at the next (lower) resolution).
        resShifterMats{res}(all(resShifterMats{res} == 0, 2), :) = [];
    end
end