function p = task_glcm(img, dx, dy, G)
%GLCM Compute the GLCM of an image. Not optimised for speed.

narginchk(3, 4);
if nargin == 3
    G = 1 + max(img(:));
end
p = zeros(G);
m_start = max(     1     ,      1     -dx);
n_start = max(     1     ,      1     -dy);
m_end   = min(size(img,1), size(img,1)-dx);
n_end   = min(size(img,2), size(img,2)-dy);
for m = m_start:m_end
    for n = n_start:n_end
        % IMG is assumed to have intensities in [0,G-1],
        % "1+" due to MATLAB's one-indexing.
        p(1+img(m,n), 1+img(m+dx,n+dy)) = p(1+img(m,n), 1+img(m+dx,n+dy)) + 1;
    end
end
p = p / ((1+m_end-m_start) * (1+n_end-n_start));