Re: Area of Intersection of two planar polygons
- To: mathgroup at smc.vnet.net
- Subject: [mg112249] Re: Area of Intersection of two planar polygons
- From: Valeri Astanoff <astanoff at gmail.com>
- Date: Mon, 6 Sep 2010 06:38:32 -0400 (EDT)
- References: <i5su76$p8d$1@smc.vnet.net>
On 4 sep, 09:58, Thomas Melehan <tpmele... at gmail.com> wrote:
> Does anyone have an algorithm for the area of the intersection of two
> planar polygons?
>
> TIA
>
> Tom=
Good day,
This is my -slow- way to do it with NIntegrate :
In[1]:=
halfLineCutsSegment[{{x0_?NumericQ, y0_?NumericQ},
{x1_?NumericQ, y1_?NumericQ}},
{{x2_?NumericQ, y2_?NumericQ},
{x3_?NumericQ, y3_?NumericQ}}]:=
(x3 (y0-y1)+x0 (y1-y3)+x1 (-y0+y3))*
(x3 (y0-y2)+x0 (y2-y3)+x2 (-y0+y3))> 0 &&
(x2 (y0-y1)+x0 (y1-y2)+x1 (-y0+y2))*
(x3 (-y0+y2)+x2 (y0-y3)+x0 (-y2+y3))> 0;
insideQ[pt:{_?NumericQ, _?NumericQ},
pol:{{_?NumericQ, _?NumericQ}..}]:=
Module[{p,rp,tp},
p = If[First[pol]==Last[pol],Most[pol],pol];
rp = RotateLeft[p];
tp = Transpose[{p,rp}];
And@@(OddQ /@ (Count[tp,seg_ /;
halfLineCutsSegment[{pt,#},seg]]& /@ (Mean /@ tp)))];
In[4]:= poly1 = {{20, 20}, {0, 20}, {0, 0}, {20, 0}, {5, 10}, {20,
20}};
In[5]:= insideQ[{10, 1}, poly1]
Out[5]= True
In[6]:= insideQ[{10, 10}, poly1]
Out[6]= False
In[7]:= poly2 = {{5, 0}, {20, 0}, {20, 10}, {5, 10}, {5, 0}};
In[8]:= xx1 = poly1[[All, 1]];
xx2 = poly2[[All, 1]];
yy1 = poly1[[All, 2]];
yy2 = poly2[[All, 2]];
xmin = Max[Min@xx1, Min[xx2]];
xmax = Min[Max@xx1, Max[xx2]];
ymin = Max[Min@yy1, Min[yy2]];
ymax = Min[Max@yy1, Max[yy2]];
In[16]:= NIntegrate[Boole[insideQ[{x, y}, poly1] && insideQ[{x, y},
poly2]],
{x, xmin, xmax}, {y, ymin, ymax}, PrecisionGoal -> 3] // Timing
Out[16]= {26.141, 75.0074}
The exact area is of course 75
--
V.Astanoff