function [ R ] = randelement(E, SZ, P)
%RANDELEMENT Random elements distributed according to a user-specified distribution.
%   R = RANDELEMENT(E, SZ) returns an array of size SZ which contains
%   random uniformly distributed elements from E.
%
%   R = RANDELEMENT(E, SZ, P) returns an array of size SZ which contains
%   random elements from E. The distribution of the random elements 
%   is according to the relative weights in P. P must be of the same
%   dimension as E.
%
%   Author: medical-image-processing.com

% assume uniform distribution of
if (nargin < 3)
   P = ones(numel(E)); 
end

% determine the discrete cumulative distribution function
CS = cumsum(P);
CS = CS ./ CS(end);

% allocate storage for the output
R = rand(SZ);

% for each output element perform the inverse transform sampling
for i=1:numel(R)
   r = R(i);
   idx = find(r<CS);
   R(i) = E(idx(1));
end

end