function [area_fraction,edges] = contribAreaHist(results,bins,units,normalized_by)
% Calculates and makes a histogram of the contributed area that a
% particular pore size makes, normalized by the image area (by default),
% or the total area of all pores, returning area fractions and edges. 
%
% [area_fraction,edges] = contribAreaHist(results,bins,units,normalized_by)
% 
% Inputs are 
% 1) results, a structure (user is prompted to choose a file if this
%   is absent or empty), 
% 2) bins, a variable which indicates the number of bins the histogram
%   should have (if bins is a scalar), or indicates the edges of the bins
%   (if bins is a vector). If empty or absent, bins defaults to 25.
% 3) units, a string, either 'pixels', 'nm', or empty.  If empty or absent,
%   defaults to 'nm'.  This means that you'll get an error if the results
%   are uncalibrated and you fail to specify 'pixels', but that's better
%   than thinking you're looking at calibrated data when you're not :)
% 4) normalized_by, a string, either 'image' or 'pores', which indicates
%   whether to normalize the contributed areas by the total image area or
%   the total pore area.  If empty or absent, defaults to 'image'.
%
% Outputs are
% 1) area_fraction, vector of histogram bar heights, the last of which will
%   always be zero (just so the lengths of area_fraction and edges match)
% 2) edges, vector of histogram bin edges
% 
% Example Usages:
%
% contribAreaHist;
%   Promts the user to select a results file, then plots a histogram with
%   25 bin of contributed area/image area vs pore diameter in nm. Reports
%   an error if the results chosen don't have a nm calibration.
%
% contribAreaHist(results);
%   Same as previous, except uses the supplied results structure intstead
%   of prompting for a file
%
% contribAreaHist(results,40);
%   Same as previous, except uses 40 bins instead of the default 25.
%
% contribAreaHist(results,40,'pixels');
%   Same as previous, except uses units of pixels rather than nm for the
%   diameters and bin edges
%
% contribAreaHist(results,40,'pixels','pores');
%   Same as previous, except normalized by the total pore area instead of
%   the total image area. 
%
% contribAreaHist(results,5:5:50);
%   Sets the histogram bin edges directly, every 5nm from 5 to 50 nm,
%   rather than just specifying the number of bins.  This allows for
%   non-equally spaced bins, but more importantly it would allow for the
%   same bin edges to be set for two different data sets, so that you can
%   really compare them side by side. 
%
% [area_fraction,edges] = contribAreaHist(...);
%   Returns a vector of area fractions (the heights of the histogram bars),
%   and the histogram edges. 
%
% [...] = contribAreaHist([],[],'pixels')
%   Same as with no arguments, except using 'pixels' instead of 'nm' for
%   units (i.e. prompt for results, default number of bins, and normalized
%   by total image area). 


if nargin==0 || isempty(results)
    %Load results from file
    [fname,pname] = uigetfile('*_r.mat','Load Results File');
    r = load([pname fname],'results');
    results = r.results;
end

% Set up defaults
if nargin<2 || isempty(bins)
    bins = (5:2:50);
end
if nargin<3 || isempty(units)
    units = 'nm';
end
if nargin<4 || isempty(normalized_by)
    normalized_by = 'image';
end

% Set units
switch units
    case {'nm','n'}
        ppn = results.pixpernm;
        if isempty(ppn)
            error('Results selected for histogram do not have a calibration, please request pixel data instead');
        end
        calibrated = 1; 
        units_label = '(nm)';
    case {'pixels','pix','p','pixel'}
        ppn = 1; 
        calibrated = 0;
        units_label = '(pixels)';
    otherwise
        error(['''' units ''' is not a valid value for specifying units']);
end

% Get diameters
diams = [results.pores.EquivDiameter]/ppn;

% Get areas
areas = [results.pores.Area]/ppn^2;

% Find bin edges for the diameters
if length(bins)==1
    % use hist to find bin edges automatically
    [dum,x] = hist(diams,bins);
    % x is bin centers, need edges
    wid = x(2)-x(1);
    edges = [x-wid/2, max(diams)+eps(max(diams))]; % have to add the final edge, and make sure it includes the maximum value
    edges(1) = min(diams);
else
    % user has specified bin edges
    edges = bins;
end
     
% Now convert bin edges to be area cutoffs instead
areaedges = pi*(edges/2).^2;

% Check for areas outside of the bounds
if any(areas<areaedges(1)) || any(areas>=areaedges(end))
    warning('Histogram bin edges do not cover all the data, data outside the bounds supplied has been ignored');
end

% Calculate contributed area for each bin
for i = 1:length(areaedges)-1
    inthisbin = areas(areas>=areaedges(i) & areas<areaedges(i+1)); 
    contrib_area(i) = sum(inthisbin);
end
contrib_area(length(areaedges)) = 0; % to match the style of histc.m.  (edges and contrib_area have to come out the same length)

% Normalize by appropriate quantity
switch normalized_by
    case 'image'
        % normalize by image area
        im_area = numel(results.L); % pixel area
        if calibrated
            im_area = im_area/ppn^2;
        end
        area_label = 'Total Image Area';
        area_fraction = contrib_area/im_area;
    case 'pores'
        % normalize by total pore area
        total_pore_area = sum(areas);
        area_fraction = contrib_area/total_pore_area;
        area_label = 'Total Pore Area';
    otherwise
        error('''normalized_by'' must be one of ''image'' or ''pores''');
end

    
% Make the actual histogram
figure; 
bar(edges,area_fraction,'histc');
xlabel(['Pore Diameter ' units_label],'FontSize',18);
ylabel(['Contributed Area / ' area_label],'FontSize',18);
title([fname(1:end-6) ' [' num2str(results.nmstats.mean_diam,3) ' nm avg ' num2str(results.stats.porosity*100,3) '% porosity]'],'FontSize', 18,'Interpreter','none');
set(gca,'FontSize',18);
xlim([0 45]);
ylim([0 6e-3]);