MATLAB: Variograms on Images

digital image processingdirectional variogramsgeostatisticsImage Processing ToolboxMATLABvariogram

I am working on applying geostatistics on Landsat imagery. My focus is on applications of variograms. I will be using n x n windows for generation of variogram plots and texture layers. From various online forums I learnt that the process would be faster if functions like blkproc, nlfilter, colfilt etc. are used instead of normal for loop based moving windows.
I see from the MATLAB help that I can integrate other custom functions to above said functions. But as I have to calculate directional variograms – in EW, NS, NE-SW and NW-SE directions I will have to consider only few pixels(cells of matrix) in front of the central pixel and not all surrounding pixels as in the case of filters. Can anyone suggest how I can go about/are there functions(or toolboxes) available to do such operations?

Best Answer

If I've understood correctly what you need, I'd use a function something like the following (which doesn't use blkproc or any of that family, but which does use MATLAB's vectorisation capability effectively).
function v = directionalVariogram(img, xoffset, yoffset)
%directionalVariogram computes empirical direction variogram
%   v = directionalVariogram(img, xoffset, yoffset) takes a 2D image array
%   and offsets in the x and y directions. It returns the mean of the
%   squared differences between pairs of pixels in the image such that the
%   spatial offsets between the pixels are as specified.
if xoffset < 0  % difference is symmetric so can force xoffset positive
    xoffset = -xoffset;
    yoffset = -yoffset;
end
% make offset and trimmed copies of image
if yoffset > 0
    imga = img(1+yoffset:end, 1+xoffset:end);
    imgb = img(1:end-yoffset, 1:end-xoffset);
else
    imga = img(1:end+yoffset, 1+xoffset:end);
    imgb = img(1-yoffset:end, 1:end-xoffset);
end
d = imga - imgb;
v = mean(d(:).^2);
end
To get the estimate for a distance of sqrt(2) pixels in the NE-SW direction, you call it like this, for example
v = directionalVariogram(img, 1, -1)
It will work with offsets in any direction and larger offsets (up to the size of the image) but the offsets have to be integers.
If you can't see how the function works, and you want to understand it, try putting in some calls to imshow to display imga and imgb, and call it with large offsets. You should also check it on some synthetic test data, for which you know what the answer should be, to make sure it does exactly what you expect.
Related Question