function measures = fit_distribution(r,show_histogram)
% FIT_DISTRIBUTION     Fits lognormal distribution to data
%   MEASURES = FIT_DISTRIBUTION(R,SHOW_HISTOGRAM)
%   R is a vector of data
%   SHOW_HISTOGRAM is a flag idicating whether plot should be produced
%   MEASURES are valuses extracted from fitted distribution: median 
%   (geometric mean), geometric standard deviation, mean, variance and mode
%
%   Copyright 2018 Vedrana Andersen Dahl, vand@dtu.dk 

measures = NaN(1,5);
if numel(r)>1 
    pd = fitdist(r(:),'lognormal'); % pd.mu and pd.sigma are mean and standard deviation of the associated normal distribution
    measures(1) = exp(pd.mu); % median of the fitted distribution (geometric mean)
    measures(2) = exp(pd.sigma); % geometric standard deviation of the fitted distribution
    measures(3) = exp(pd.mu+pd.sigma.^2/2); % mean of lognormal distribution
    measures(4) = (exp(pd.sigma.^2)-1)*measures(3).^2; % variance of the lognormal distribution
    measures(5) = exp(pd.mu-pd.sigma.^2); % mode of the lognormal distribution 
end

if nargin>1 && show_histogram && numel(r)>1
    h = histfit(r(:),ceil(max(r(:))),'lognormal'); hold on
        x = h(2).XData;
    y = h(2).YData;
    title({['mean = ',num2str(measures(3)),', std =  ',num2str(sqrt(measures(4)))],...
        ['median = ',num2str(measures(1)),', gstd =  ',num2str(measures(2))],...
        ['mode = ',num2str(measures(5))]})
    l(1) = plot(measures(5)*[1,1],[0,interp1(x,y,measures(5))],'c--');    
    l(2) = plot(measures(1)*[1,1],[0,interp1(x,y,measures(1))],'g--');    
    l(3) = plot(measures(3)*[1,1],[0,interp1(x,y,measures(3))],'m--');    
    xlim([0,max(x)])
    xlabel('radius (pixels)'), ylabel('count')
    legend(l,{'mode','median (geometric mean)','mean'},'Location','NE')    
end