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(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