function pltdvh(x, K, target, oar, dmax)
%
% pltdvh(x, K, target, oar, dmax)
%
% Function for plotting Dose-Volume Histograms (DVH)
%
% x      : the beamlets
% K      : the influence matrix
% target : structure for the target
% oar    : vector of structures for the organs
%          Brain stem, right parotid gland, left parotid gland
%          and spinal cord.
% dmax   : maximum dose for the plot (e.g. 80 Gy)
%
%

dose = K*x(:);
if (nargin < 4)
  error('Four or five input arguments are required.')
end

if nargin == 4
  dmax = 80;
end
Y = [];
% for target
x = x(:);
d = dose(target.ind);
y = dvh(d, dmax);
Y = [Y y];
L{1} = 'PTV';

% for Parotoid Left
d = dose(oar(3).ind);
y = dvh(d, dmax);
Y = [Y y];
L{2} = 'Parotoid L';

% for Parotoid Right
d = dose(oar(2).ind);
y = dvh(d, dmax);
Y = [Y y];
L{3} = 'Parotoid R';


% for Spinal cord
d = dose(oar(4).ind);
y = dvh(d, dmax);
Y = [Y y];
L{4} = 'Spinal cord';


% for Brainstem
d = dose(oar(1).ind);
y = dvh(d, dmax);
Y = [Y y];
L{5} = 'Brainstem';

figure;
t = 0:dmax;
plot(t,Y)
legend(L)
grid on
xlabel('Dose (Gy)')
ylabel('Fractional vol.')

%+++++++++++++++++++++++++++
function y = dvh(d,dmax)

m = length(d);
t = 0:dmax; t = t'; n = length(t);
y = zeros(n,1);
for i=n:-1:1
  I = find(d >= t(i));
  y(i) = y(i) + length(I);
  y(1:i-1)=y(1:i-1) + length(I);
  ind = setdiff(1:length(d),I);
  d = d(ind);
end
y = y/m;
