Code covered by the BSD License

Monte Carlo method for estimating plane figure area

Konstantin Ninidze (view profile)

An area S(A) of a figure A, bounded by a Jordan curve is calculated.

S=Area_M(N)
```function S=Area_M(N)
% by Konstantin Ninidze
%
% Monte Carlo method for estimating the area for a plane figure.
% Suppose we need to calculate an area S(A) of a figure A, bounded by a
% Jordan curve (which in our case is  constructed by a cubic spline
% approximation). We'll place A inside a square of known area S(sq) and
% then place a known number of points (N) at random locations inside the
% square. After this we'll count the number of the random points (M) that lie
% inside the contour A and taking in mind, that the area S(A) of the object
% A is proportional to the number of points (M) that lie inside square
% we'll obtain a formula for the area: S(A)=S(sq)*M/N.

% A slight modification of Matlab "getcurve" function is needed to obtaine
% the closed curve, see mygetcurve1.
close
X=mygetcurve1;
X=[X,X(:,1)];
x=2*rand(1,N)-1;
y=2*rand(1,N)-1;
in = inpolygon(x,y,X(1,:),X(2,:));
M=sum(in);
S=4*M/N;
hold on
plot(x(in),y(in),'r+',x(~in),y(~in),'bo')
title([' area estimation is ',num2str(S)])
xlabel([' N= ',num2str(N)])

function [xy, spcv] = mygetcurve1
% do not repeat points!

w = [-1 1 -1 1];
clf,
axis(w), hold on, grid on

title('Use mouse clicks to pick points INSIDE the gridded area.')
pts = line('Xdata',NaN,'Ydata',NaN,'marker','o','erase','none');

maxpnts = 100; xy = zeros(2,maxpnts);
while 1
for j=1:maxpnts
[x,y] = ginput(1);
if isempty(x)||x<w(1)||x>w(2)||y<w(3)||y>w(4), break, end
xy(:,j) = [x;y];
if j>1
set(pts,'Xdata',xy(1,1:j),'Ydata',xy(2,1:j))
else
set(pts,'Xdata',x,'Ydata',y)
xlabel('When you are done, click OUTSIDE the gridded area.')
end
end
if j>1, break, end
xlabel(' You need to click INSIDE the gridded area at least once')
end

xy(:,j:maxpnts)=[];
if norm(xy(:,1)-xy(:,j-1))<.05, xy(:,j-1)=xy(:,1); end
set(pts,'Xdata',xy(1,:),'Ydata',xy(2,:),'erase','xor','linestyle','none')
xy=[xy,xy(:,1)];
spcv = cscvn(xy); fnplt(spcv)

```