The area of a polygon (or any plane shape) can be evaluated by Monte Carlo integration. The process involves 'random' sampling of (x,y) points in the xy plane, and testing whether they fall inside or outside the shape. That is illustrated schematically below for a circle.

^

Steps:

1. Define a 2D region within which to sample points. In the present problem you must choose the rectangular box (aligned to the axes) that bounds the specified polygon.
2. By a 'random' process generate coordinates of a point within the bounding box. In the present problem a basic scheme using the uniform pseudorandom distribution available in MATLAB must be employed. Other schemes in use include quasi-random sampling (for greater uniformity of coverage) and stratified sampling (with more attention near the edges of the polygon).
3. Determine whether the sampled point is inside or outside the polygon. Due to working in double precision, it is extremely unlikely to sample a point falling exactly on the edge of the polygon, and consequently if it does happen it doesn't really matter whether it's counted as inside or outside or null.
4. Repeat steps 2–3 N times.
5. Based upon the proportion of sampled points falling within the polygon, report the approximate area of the polygon.

Inputs to your function will be N, the number of points to sample, and polygonX & polygonY, which together constitute an ordered list of polygon vertices. N will always be a positive integer. polygonX & polygonY will always describe a single, continuous outline; however, the polygon may be self-intersecting.

For polygons that are self-intersecting, you must find the area of the enclosed point set. In the case of the cross-quadrilateral, it is treated as two simple triangles (cf. three triangles in the first figure below), and the clockwise/counterclockwise ordering of points is irrelevant. A self-intersecting pentagram (left & middle of the second figure below) will therefore have the same area as a concave decagon (right).

^

EXAMPLE

 % Input
 N = 1000
 polygonX = [1 0 -1  0]
 polygonY = [0 1  0 -1]
 % Output
 A = 2.036

Explanation: the above polygon is a square of side-length √2 centred on the origin, but rotated by 45°. The exact area is hence 2. As the value of N is moderate, the estimate by Monte Carlo integration is moderately accurate (an error of 0.036 in this example, corresponding to 509 of the 1000 sampled points falling within the polygon). Of course, if the code were run again then a slightly different value of A would probably be returned, such as A=1.992 (corresponding to 498 of the 1000 sampled points falling within the polygon), due to the random qualities of the technique.