function [V, Ci] = area_mc(px,py,rx,ry,N)
% function [V, Ci] = area_mc(px,py,rx,ry,N)
% This function estimates the area inside a polygon (px,py) using 
% Monte Carlo simulation. It returns the estimate of the volume as well as 
% the lower and upper values of confidence interval.
%
% ---------------------------------------------------------------
% INPUT
% px:       x-coordinate values of the polygon (nx1 double)
% py:       y-coordinate values of the polygon (nx1 double)
% rx:       x-coordinate values of the surrounding rectangle (4x1 double)
% yx:       y-coordinate values of the surrounding rectangle (4x1 double)
% N:        sample size (1x1 double)
% ---------------------------------------------------------------
% OUTPUT
% V:        Estimated area  (1x1 double)
% Ci:       Confidence intervals (2x1 double) [lower, upper]
% ---------------------------------------------------------------
% Created:
% 2018-02-19 Heikki Puustinen

% Determine minimum and maximum x- and y- values of the sourrounding
% reactangle.
xMin = min(rx);
xMax = max(rx);
yMin = min(ry);
yMax = max(ry);

% Generate 2N random numbers (uniform distribution U(0,1)).
sx = rand(N,1);
sy = rand(N,1);

% Scale to fit the reactangle (This is needed if rectancle is not a unit rectangle).
sx = xMin + (xMax - xMin) * sx;
sy = yMin + (yMax - yMin) * sy;

% Check which points are inside polygon (px,py). 
% 1: inside, 0: outside
in = inpolygon(sx,sy,px,py);

% Number of points inside the polygon.
n_in = sum(in);

% Fraction of points inside the polygon. 
% (Fraction = number of points in / number of sample points)
lambda = n_in / N; 

% Area of the surrounding rectangle 
% (Area = a * b)
A = (xMax - xMin) * (yMax - yMin);

% Value of the integral 
% (Area of the rectangle * fraction of points in).
V = A * lambda;

% To calculate confidence interval, we need variance.
% Each point is either inside or outside 
% => use sample variance of the binomial distribution:
% p^hat = X/n, var(X/n) = var(X)/n^2 
% var(p^hat) =  p * (1 - p) / n, n = sample size, p = lambda. 
var_V = lambda * (1 - lambda) / N;

% Using central limit theorem, approximate distribution of the integral
% being normal => confidence interval is: 
V_upper = V + 1.96 * sqrt(var_V);
V_lower = V - 1.96 * sqrt(var_V); 
% 1.96 is the value at which the inverse of the enormal cumulative
% distribution is 0.975 (1 - 0.05/2). Confidence level 95% -> alpha = 0.05

% Combine lower and upper values of the interval.
Ci = [V_lower, V_upper];
end