function [H_modbox,D_modbox] = frac_modbox(p,dividers)

%% License:
% "This work is licensed under the Creative Commons Attribution 4.0
% International (CC BY 4.0) License. To view a copy of this license, visit:
% http://creativecommons.org/licenses/by/4.0/"

%% please cite 
% --> Kristofer Marsch (2020)
% --> DOI for this data available at https://orcid.org/0000-0002-1857-506X


%% calculating Hurst Exponent and Fractal Dimension using modified box counting
% based on Malinverno (1990)

%% --parameters-- %%
% p: array of profile, first column horizontal values, second column heights
% profile must be de-trended
% dividers: number of divisions of axes, e.g. dividers = [2 4 8 16 32 64];

%% calculate H and D
xpoints = p(:,1)-min(p(:,1)); %move to first quadrant
ypoints = p(:,2)-min(p(:,2));

F = griddedInterpolant(xpoints,ypoints); %interpolate, as polygons visited by points will be counted
xq = linspace(min(xpoints),max(xpoints),10000); 
yq = F(xq);

xrange = max(xpoints)-min(xpoints); %define boxranges
yrange = max(ypoints)-min(ypoints);

%count boxes
count_div = 1;
result = zeros(length(dividers),3); %preallocate
for boxdiv = dividers

    x = min(xpoints): xrange/boxdiv :max(xpoints);
    y = min(ypoints): yrange/boxdiv :max(ypoints);
    [X,Y] = meshgrid(x,y);

    boxvisit = zeros(boxdiv^2,2); %preallocate
    count_box = 1;
    
    for i1 = 1:length(x)-1 %rows of grid
        
        for i2 = 1:length(x)-1 %columns of grid
            
            p1 = [X(i2,i1) Y(i2,i1)];
            p2 = [X(i2,i1+1) Y(i2,i1+1)];
            p3 = [X(i2+1,i1+1) Y(i2+1,i1+1)];
            p4 = [X(i2+1,i1) Y(i2+1,i1)];
            p5 = [X(i2,i1) Y(i2,i1)];

            xps = [p1(:,1);p2(:,1);p3(:,1);p4(:,1);p5(:,1)];
            yps = [p1(:,2);p2(:,2);p3(:,2);p4(:,2);p5(:,2)];

            [~,col] = find(xq >= p1(:,1) & xq < p2(:,1));
            usexq = xq(1,col);
            useyq = yq(1,col);
                        
            [in,on] = inpolygon(usexq,useyq,xps,yps); %points in box?
            
            boxvisit(count_box,1) = numel(xq(in));
            boxvisit(count_box,2) = numel(yq(on));
            count_box = count_box+1;
            
            usexq=[];
            useyq=[];
        
        end
        
    end
    
    no_box = nnz(boxvisit(:,1));

    result(count_div,1) = boxdiv;
    result(count_div,2) = no_box;
    result(count_div,3) = xrange/boxdiv;
    count_div = count_div+1;

end

X2 = [ones(length(log10(result(:,1))),1) log10(result(:,1))];
c = X2\log10(result(:,2));

D_modbox = c(2);
H_modbox = 2-D_modbox;
