function DVH = dvh(dose,mask,d_bins)
% A dose volume histogram is calculated for a region of interest (roi)
% using dose bins specified by the user.
%
% Inputs:  dose = dose distribution
%          mask = mask of zeros and ones that specifies the roi
%          d_bins = dose bin boundaries, must be a vector
% 
% Output:  d_hist = dose volume histogram defined as:
%          d_hist(j) = percentage of voxels in the region of interested 
%                      receiving a dose of d_bins(j) or greater.

% ensure that d_bins is a vector
if prod(size(d_bins)) > length(d_bins)
    error('Third input argument must be a vector.');
end

% ensure that the mask is just ones and zeros
if ~isempty(find(mask~=0 & mask~=1))
    error('Mask must consist of zeros and ones only.');
end

% ensure that the mask and the dose matrix are the same size
if size(dose,1) ~= size(mask,1) | size(dose,2) ~= size(mask,2) ...
        | size(dose,3) ~= size(mask,3)
    error('The dose matrix and the mask must be the same size.');
end

Nbins = length(d_bins);  % number of bins in the dvh
Nvoxels = sum(mask(:));  % number of voxels in the mask

d_hist = zeros(1,length(d_bins));
DVH = zeros(size(d_hist));

ind = find(mask == 1);  % indices that lie in the roi

if ~isempty(ind)   
    dose_vector = dose(ind);  % vector of dose
    d_hist = histc(dose(ind),d_bins);  % create the differential histogram
    d_hist = 100*d_hist/length(ind); % renomalize
    DVH = flipud(cumsum(flipud(d_hist)));
end