function varargout = discretisepositiondata( positionData, maxDimSize, testPositions )
    % Form the discrete 2D position data from processed position data and
    % the coordinates of 
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % discreteCoordsTrajectory will be a n by 2 matrix of positive integers
    % - 2D discrete grid coordinates.
    %
    % yDim, xDim - the size of the discretisation grid we use in the Y- and
    % X-dimensions respectively; i.e. the number of grid squares that cover
    % the environment in each dimension.
    %
    % positionData must already have been processed to have time in seconds
    % in the first column, zeroed on the same scale as spike times. The
    % second and third columns are the Y and X pixel coordinates,
    % respectively.
    %
    % maxDimSize is the number of discrete squares that will tile the
    % longest spatial dimension in the desired resolution. Must be a
    % positive integer.
    %
    % obsTimeInterval - if supplied, fixes a constant time interval between
    % observations. Starting with the first observation, we use linear
    % interpolation to obtain position observations at the imposed time
    % points.
    %
    % 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
    % 29.09.2014 - copied from old version.
    % 30.09.2014 - removed references to time indexing - we do not
    % discretise time. Removed linear interpolation (in time). Added
    % spatialBinWidth to output, for computing the centre coordinates of
    % grid cells later. Removed requirement for position data to be
    % positive.
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    if nargin < 2,
        error('You must at least supply the matrix of observed positions and a discretisation size!')
    end
    if size(positionData, 2) ~= 3,
        error('positionData must have 3 columns: time (seconds), Y coord and X coord!')
    end
    if rem(maxDimSize, 1) ~= 0,
        error('maxDimSize must be a positive integer!')
    end
    if maxDimSize < 1,
        error('maxDimSize must be a positive integer!')
    end

    % Convert to the discrete grid according to the specified resolution.
    % Get the spatial bin size, in pixels.
    
    if nargin == 3,
        dimSize = abs(max([positionData(:, 2:3); testPositions], [], 1) - min([positionData(:, 2:3); testPositions], [], 1));
    else
        dimSize = abs(max(positionData(:, 2:3), [], 1) - min(positionData(:, 2:3), [], 1));
    end

    maxY = dimSize(1);
    maxX = dimSize(2);
    maxXY = max([maxY, maxX]);
    % In case we had one position: dimension size should still be nonzero.
    if maxXY == 0,
        maxXY = 1;
    end
    spatialBinWidth = maxXY / maxDimSize;
    
    % Notice that if the discretisation size is the same as the resolution
    % of the data, we don't have to do anything.
    discreteCoordsTrajectory = floor(positionData(:, 2:3) / spatialBinWidth);

    if nargin == 3,
        testPositions = floor(testPositions / spatialBinWidth);
        testPositions = bsxfun(@plus, testPositions, 1 - min([discreteCoordsTrajectory; testPositions], [], 1));
        
        discreteCoordsTrajectory = bsxfun(@plus, discreteCoordsTrajectory, 1 - min([discreteCoordsTrajectory; testPositions], [], 1));

        % Get discrete resolution dimensions.
        dims = max([discreteCoordsTrajectory; testPositions], [], 1);
        
        % Use this for updated asymptotic means in simulated data.
        testPositions
    else
        % May need to re-zero, depending on the resolution.
        discreteCoordsTrajectory = bsxfun(@plus, discreteCoordsTrajectory, 1 - min(discreteCoordsTrajectory, [], 1));
   
        % Get discrete resolution dimensions.
        dims = max(discreteCoordsTrajectory, [], 1);
    end
     
    yDim = dims(1);
    xDim = dims(2);
    
    if nargout == 4,
        varargout = {discreteCoordsTrajectory, yDim, xDim, spatialBinWidth};
    elseif nargout == 5,
        varargout = {discreteCoordsTrajectory, yDim, xDim, spatialBinWidth, testPositions};
    end
    
end