areaintersection.m

function [area]=areaintersection(set1, set2, resolution)

% Function gives the approximate area of intersection of two polygons
%     whose vertices are given by set1 and set2. 
% 
% Note that implementation of this function requires the standard Matlab
%     function "inpolygon"
% 
% STANDARD CALL
% 
% area = areaintersection(set1, set2, resolution)
%     
%     
% INPUTS
% 
% set1 : an N x 2 matrix of the form (x1, y1; x1, y2; ...) giving the
%     vertices of a polygon
% set2 : an M x 2 matrix of the form (x1, y1; x1, y2; ...) giving the
%     vertices of a second polygon
% resolution : a positive scalar that determines the accuracy of the
%     approximation. The larger the value the greater the accuracy and the 
%     longer the evaluation time. The matlab function "inpolygon" is the
%     real limiting factor in terms of execution time
%     
% OUTPUT
% 
% area : a scalar value that represents the area of intersection of the two
%     polynomials in standard units squared
% 
% 
% 
%   Paul Koprowski 2007
%   paulkoprowski@hotmail.com
%   
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % 

% Checks the input arguments
if nargin~=3
    error('wrong number of inputs');
end
sset1=size(set1);
sset2=size(set2);
if sset1(1,1)<2
    error('set1 has too few rows');
end
if sset2(1,1)<2
    error('set2 has too few rows');
end
if sset1(1,2)<2
    error('set1 has too few columns');
end
if sset2(1,2)<2
    error('set2 has too few columns');
end
s=size(resolution);
if s(1,1)~=1
    error('resolution must be a scalar');
end
if s(1,2)~=1
    error('resolution must be a scalar');
end
% creates a hull of all point in both polygons
hset = [set1;set2];
    % Extreme boundaries of hull
maxx=max(hset(:,1));
minx=min(hset(:,1));
maxy=max(hset(:,2));
miny=min(hset(:,2));

% Generates matrix containing all points of both polygons and then some

stepwidthx = (maxx-minx)/resolution;
stepwidthy = (maxy-miny)/resolution;
hullx = ones(resolution+1,1)*[minx:stepwidthx:maxx];
hully = ([miny:stepwidthy:maxy]')*ones(1,resolution+1);

% Uses "inpolygon" function to determine points inside each polygon

[poly1] = inpolygon(hullx,hully,set1(:,1),set1(:,2));
[poly2] = inpolygon(hullx,hully,set2(:,1),set2(:,2));

% Finds the intersection points of the two polygons

int = and(poly1,poly2);

% Multiplies the area of each box in the hull by the number of boxes
% contained in each polygon to approximate the total area of intersection

area = (stepwidthx*stepwidthy*sum(sum(int)));