% Plot probability density from a ray tracing. Based on density.m, but
% this version uses stratified sampling to avoid creating Moire patterns.
function density_image_ss = density_ss(xs, ys, L, delta)

N = numel(xs);

% Set up density density_image coordinates
Nx = 1+ceil((2+L)/delta);
Ny = 1+ceil(2/delta);
xcoord = (0:(Nx-1))*delta;
ycoord = (0:(Ny-1))*delta;
density_image_ss = zeros(Nx,Ny);

for i=1:N-1 % ith ray connects vertex i to vertex i+1.
   xlength = xs(i+1)-xs(i);
   ylength = ys(i+1)-ys(i);
   length = sqrt( xlength.^2 + ylength.^2 );
   Npoints = ceil(3*length/delta);
   s = ((0:Npoints-1)+rand(1,Npoints))/Npoints; % fractional arclength coordinate along ray.
      % The random displacement within each interval is an implementation of stratified sampling.
   xarr = xs(i) + s*xlength; % array of Npoints x-coordinates, equally spaced along ray
   yarr = ys(i) + s*ylength; % array of Npoints y-coordinates, equally spaced along ray
   xsub = ceil( (xarr + (1+L/2))/delta ); % array of x indices in density_image array
   ysub = ceil( (yarr + 1)/delta ); % array of y indices in density_image array
   density_image1Dsubs = Nx*(ysub-1) + xsub;
   density_image_ss(density_image1Dsubs) = density_image_ss(density_image1Dsubs) + 1;
end

hold off
imagesc(density_image_ss')
axis equal



