function constrainedDistanceMat = computeeuclideanaccessibledistance( validPositionsLogMat, spatialBinWidth, returnAdjacency )
    % Compute the distance between every pair of positions supplied in a
    % matrix of logical indicators of where the 'valid' positions are, in a
    % grid. The resolution can be arbitrarily precise.
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Uses the Johnson or Floyd-Warshall algorithm, depending on sparsity
    % (as implemented in the MatlabBG1 package).
    %
    % constrainedDistanceMat is a sparse upper triangular matrix: the
    % distance from position j to position i should be obtained by
    % extracting the distance from i to j using symmetry. This is the
    % distance in continuous units of space, i.e. in pixels.
    %
    % validPositionsLogMat is a m by n matrix of logicals, for which the
    % (i, j) entry is true if the corresponding position is considered
    % 'valid', or false otherwise.
    %
    % spatialBinWidth - width of a discrete grid square, in pixels.
    %
    % returnAdjacency - logical scalar. If true, we return the weighted
    % adjacency matrix rather than the distance matrix.
    %
    % Note: positions are indexed by column appearance, i.e. following
    % linear indexing rules, and only including valid positions. Thus, if
    % here are 5 positions with an X coordinate of 1, the position (Y, X) =
    % (1, 2) will have an index of 6: the 6th valid position.
    %
    % Y always indexes rows of position matrices, as the vertical spatial
    % dimension label. X is the horizontal dimension label, so it indexes
    % columns of position matrices.
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % CHANGELOG
    % 30.09.2014 - copied from old version. Implemented distance in
    % continuous (pixel) units.
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    if nargin < 2,
        error('Insufficient arguments!')
    end
    if nargin < 3,
        returnAdjacency = false;
    end
    if ~all(all(islogical(validPositionsLogMat))),
        error('validPositionsLogMat must be a logical matrix!')
    end
    
    [yDim, xDim] = size(validPositionsLogMat);
    
    % First we must derive the weighted adjacency matrix for the set of
    % valid positions.

    linearInds = find(validPositionsLogMat);
    spatialDim = length(linearInds);
    linearIndsMat = zeros([yDim, xDim]);
    linearIndsMat(linearInds) = 1:spatialDim;
    
    tri = (spatialDim - 1) * spatialDim / 2;
    weightsList = zeros([tri, 1]);
    upperDistanceMatCoords = zeros([tri, 2]);

    % Loop over valid positions.
    ind = 0;
    for sourceInd = 1:spatialDim,
        % Get coordinates of this position.
        [i, j] = ind2sub([yDim, xDim], linearInds(sourceInd));
        % Look at 'later' positions around this position:
        % First, below.
        if (i < yDim) && (validPositionsLogMat(i + 1, j)),
            ind = ind + 1;
            upperDistanceMatCoords(ind, :) = [sourceInd, linearIndsMat(i + 1, j)];
            weightsList(ind) = 1;
        end
        % Second, up and right.
        if (i > 1) && (j < xDim) && (validPositionsLogMat(i - 1, j + 1)),
            ind = ind + 1;
            upperDistanceMatCoords(ind, :) = [sourceInd, linearIndsMat(i - 1, j + 1)];
            weightsList(ind) = sqrt(2);
        end
        % Third, right.
        if (j < xDim) && (validPositionsLogMat(i, j + 1)),
            ind = ind + 1;
            upperDistanceMatCoords(ind, :) = [sourceInd, linearIndsMat(i, j + 1)];
            weightsList(ind) = 1;
        end
        % Fourth, down and right.
        if (i < yDim) && (j < xDim) && (validPositionsLogMat(i + 1, j + 1)),
            ind = ind + 1;
            upperDistanceMatCoords(ind, :) = [sourceInd, linearIndsMat(i + 1, j + 1)];
            weightsList(ind) = sqrt(2);
        end
    end
    
    weightsList = weightsList * spatialBinWidth;

    % Get rid of zero entries - position pairs that aren't adjacent - and
    % copy the upper part into the lower part.
    weightsList = [weightsList(1:ind); weightsList(1:ind)];
    upperDistanceMatCoords = [upperDistanceMatCoords(1:ind, :); upperDistanceMatCoords(1:ind, [2, 1])];

    % Create sparse weighted adjacency matrix.
    weightedAdjacencyMat = sparse(upperDistanceMatCoords(:, 1), upperDistanceMatCoords(:, 2), weightsList, spatialDim, spatialDim);
    if returnAdjacency,
        constrainedDistanceMat = weightedAdjacencyMat;
        return;
    end

    % Now we can compute the distances using MatlabBG1's graph algorithms.
%     constrainedDistanceMat = sparse(triu(all_shortest_paths(weightedAdjacencyMat)));
%     constrainedDistanceMat = triu(all_shortest_paths(weightedAdjacencyMat));
    constrainedDistanceMat = all_shortest_paths(weightedAdjacencyMat);

    % Alternative by Dijkstra's algorithm (slightly slower).
%     constrainedDistanceMat = zeros(spatialDim, spatialDim);
%     for ind = 1:spatialDim,
%         constrainedDistanceMat(:, ind) = shortest_paths(weightedAdjacencyMat, ind);
%     end
%     constrainedDistanceMat = sparse(triu(constrainedDistanceMat));
end