% forest.m
% Generates a view of a forest from above as an instance of a 
% Poisson random field. The forest is contained in the unit square.
% The forest is concentrated towards the top-right corner of the square.

% Square has area 1 so generate Poisson(1 * MaxRate) r.v.
MaxRate = 1000;
N = random('Poisson', MaxRate);

% Generate the points and uniforms we will use to accept/reject.
% First two columns give XY locations. Third column is for accept/reject.
XY = rand(N, 3);

% Perform the acceptance/rejection
RateAtXY = MaxRate * (XY(:, 1) + XY(:, 2)).^2 / 4;
AcceptPoint = MaxRate * XY(:, 3) <= RateAtXY;
AcceptedPoints = XY(AcceptPoint > 0, 1:2);
plot(AcceptedPoints(:, 1), AcceptedPoints(:, 2), 'g*');
