function [dvh, dosebins] = dvhist(varargin)
% DVHIST Plot a DVH and/or compute the statistics of the plan
%
% Usage:
%   [dvh, dosebins, stats] = dvhist (dosemap, contour)
%
% INPUT:
%   dosemap = dose distribution, 2D or 3D
%   contour = index vector or logical matrix specifying the ROI
%   mode = 'relative' (default) or 'absolute'
%
% OUTPUT:
%   [] = plot the DVH
%   dvh = vector of dvh, ordinate of dvh
%   dosebins = dose bin locations, abscissa of the dvh
%
% See also: N/A
%
% TODO: change the way handles displayname or easier to plot output
%
% Copyleft (c) Xiaohu Mo
% Version: 1.3c

if nargin == 2
    % (dosemap, contour)
    dosemap = varargin{1};
    contour = varargin{2};
elseif nargin == 4
    % (dosemap, Geometry, ROIindex)
    dosemap = varargin{1};
    Geometry = varargin{2};
    ROIind = varargin{3};
    nfraction = varargin{4};
    contour = Geometry.ROIS{ROIind}.ind;
else
    error('unknown inputs');
end


if isempty(contour)
    error('contour is empty');
else
    dosevec = dosemap(contour);
end


% calculate DVH
if isempty(find(dosevec, 1))
    % if no voxel has dose, hist return -Nbins/2 ~ Nbins/2, affects plotting
    dvh = [100 0];
    dosebins = [0 1E-3];
else
    % Use histogram to calculate DVH
    [pdf dosebins] = hist(dosevec, 999);
    % clip negative bins
    pdf = pdf(dosebins >= 0);
    dosebins = dosebins(dosebins >= 0);
    dvh = fliplr(cumsum(fliplr(pdf))) / numel(dosevec) * 100;
end
% append the last bin
dosebins = [dosebins dosebins(end)+0.1];
dvh = [dvh 0];

end