acv.m

function DEM = acv(DEM)

%ACV Anisotropic coefficient of variation (ACV) 
%
% Syntax
%
%     C = acv(DEM)
%
% Description
% 
%     The anisotropic coefficient of variation describes the general
%     geometry of the local land surface and can be used to distinguish
%     elongated from oval landforms.
%
% Input
%
%     DEM       digital elevation model (GRIDobj)
%
% Output
%
%     C         Anisotropic coefficient of variation (ACV) (GRIDobj)
%
% Example
%
%     DEM = GRIDobj('srtm_bigtujunga30m_utm11.tif');
%     A = acv(DEM);
%     imageschs(DEM,A)
%
% References
%
%     Olaya, V. 2009: Basic land-surface parameters. In: Geomorphometry. 
%     Concepts, Software, Applications, Hengl, T. & Reuter, H. I. (Eds.),
%     Elsevier, 33, 141-169.
% 
% See also: CONV2, FILTER2, BWDIST
%
% Author: Wolfgang Schwanghart (w.schwanghart[at]geo.uni-potsdam.de)
% Date: 17. August, 2017

DEM.Z = acvfun(DEM.Z);
DEM.name = 'ACV';



%% %%%%%%%%%%%%%%%%%%%%%%%
% subfunction
function C = acvfun(dem)
siz = size(dem);

dem     = padarray(dem,[2 2],nan);
inanpad = isnan(dem);
[~,L]   = bwdist(~inanpad);
dem     = dem(L);


k       = [ 1  0  1  0  1; ...
            0  0  0  0  0; ...
            1  0  0  0 -1; ...
            0  0  0  0  0; ...
           -1  0 -1  0 -1];
       
dz_AVG  = conv2(dem,k,'valid')/4;

F = {     [ 0  0  0  0  0; ...
            0  0  0  0  0; ...
            1  0  0  0 -1; ...
            0  0  0  0  0; ...
            0  0  0  0  0];...
            ...
          [ 1  0  0  0  0; ...
            0  0  0  0  0; ...
            0  0  0  0  0; ...
            0  0  0  0  0; ...
            0  0  0  0 -1];
            ...
          [ 0  0 -1  0  0; ...
            0  0  0  0  0; ...
            0  0  0  0  0; ...
            0  0  0  0  0; ...
            0  0 -1  0  0];
            ...
          [ 0  0  0  0 -1; ...
            0  0  0  0  0; ...
            0  0  0  0  0; ...
            0  0  0  0  0; ...
            1  0  0  0  0];...
     };
            
ACV = zeros(siz);
 
for r = 1:4;
    ACV = ACV + (conv2(dem,F{r},'valid') - dz_AVG).^2;
end

dem = dem(2:end-1,2:end-1);

F = {      [ 0  0  0;...
             1  0 -1;...
             0  0  0];...
             ...
           [ 1  0  0;...
             0  0  0;...
             0  0 -1];...
             ...
           [ 0  1  0;...
             0  0  0;...
             0 -1  0];...
             ...
           [ 0  0  1;...
             0  0  0;...
            -1  0  0];...
     };

for r = 1:4;
    ACV = ACV + (conv2(dem,F{r},'valid') - dz_AVG).^2;
end

dz_AVG = max(abs(dz_AVG),0.001);

C = log(1 + sqrt(ACV./8)./dz_AVG);