Monte Carlo method for estimating plane figure area

by

 

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)


Contact us