Re: testing if a point is inside a polygon
- To: mathgroup at smc.vnet.net
- Subject: [mg96433] Re: testing if a point is inside a polygon
- From: Daniel Lichtblau <danl at wolfram.com>
- Date: Fri, 13 Feb 2009 03:45:46 -0500 (EST)
- References: <gmp0mt$btn$1@smc.vnet.net>
On Feb 9, 4:31 am, Mitch Murphy <mi... at lemma.ca> wrote:
> is there a way to test whether a point is inside a polygon? ie.
>
> PointInsidePolygonQ[point_,polygon_] -> True or False
>
> i'm trying to do something like ...
>
> ListContourPlot[table,RegionFunction->CountryData["Canada=
","Polygon"]]
>
> to create what would be called a "clipping mask" in photoshop.
>
> cheers,
> Mitch
The application seems to involve many such tests rather than a one-off
sort of query, so it is probably acceptable to spend time in
preprocessing, if it allows us to avoid O(n) work per query (n being
the number of segments). I'll show an implementation that bins the
segments according to x coordinates. Specifically, a bin will contain
a segment if either starting or terminating x coordinate of the
segment corresponds to x values in the bin, or else the entire bin's
range is between the starting and terminating values. Note that a
given segment might be in multiple bins. So long as we bin "sensibly",
all bins should have but a smallish fraction of the total number of
segments.
The preprocessing is a bit on the slow side because I do nothing
fancy. I guess one could use Sort and some smarts to make it faster
(which might be a good idea, if the plan is to do this for many
different objects, e.g. all countries).
Anyway, we create the bins of segments and also keep track of min and
max x and y values (these we use as a cheap way of ruling out points
not in an obvious bounding rectangle).
polyToSegmentList[poly_, nbins_] := Module[
{xvals, yvals, minx, maxx, miny, maxy, segments, flatsegments,
segmentbins, xrange, len, eps},
{xvals, yvals} = Transpose[Flatten[poly, 1]];
{minx, maxx} = {Min[xvals], Max[xvals]};
{miny, maxy} = {Min[yvals], Max[yvals]};
segments = Map[Partition[#, 2, 1, {1, 1}] &, poly];
flatsegments = Flatten[segments, 1];
xrange = maxx - minx;
eps = 1/nbins*len;
len = xrange/nbins;
segmentbins = Table[
lo = minx + (j - 1)*len - eps;
hi = minx + j*len + eps;
Select[flatsegments,
Module[{xlo, xhi}, {xlo, xhi} = Sort[{#[[1, 1]], #[[2, 1]]}];
lo <= xlo <= hi || lo <= xhi <= hi || (xlo <= lo && xhi >=
= hi)
] &],
{j, nbins}];
{{minx, maxx}, {miny, maxy}, segmentbins}
]
We preprocess for Canada.
In[125]:= canpoly = First[CountryData["Canada", "Polygon"]];
In[256]:= nbins = 100;
Timing[{{xmin, xmax}, {ymin, ymax}, segmentbins} =
polyToSegmentList[canpoly, nbins];]
Out[257]= {47.0648, Null}
A bit slow, but it's only done once. And it covers all 25 pieces, not
just the main body.
To test a point that is not outside the bounding rectangle, we will
drop a vertical from it to some point below the minimal y value, then
count intersections between that vertical and all segments in the bin
corresponding to the point's x coordinate.
pointInPolygon[{x_, y_}, bins_, xmin_, xmax_, ymin_, ymax_] :=
Catch[Module[
{nbins = Length[bins], bin, seglist, crosses},
If[x < xmin || x > xmax || y < ymin || y > ymax, Throw[False]];
bin = Ceiling[nbins*(x - xmin)/(xmax - xmin)];
seglist = bins[[bin]];
crosses = countIntersections[seglist, {x, y, ymin - 1}];
If[EvenQ[crosses], False, True]
]]
countIntersections[segs_, {x_, yhi_, ylo_}] := Module[
{tally = 0, seg, yval, xlo, xhi, y1, y2},
Do[
seg = segs[[j]];
{{xlo, y1}, {xhi, y2}} = seg;
If[(x < xlo && x < xhi) || (x > xlo && x > xhi), Continue[]];
yval = y1 + (x - xlo)/(xhi - xlo)*(y2 - y1);
If[ylo < yval < yhi, tally++];
, {j, Length[segs]}];
tally
]
Here are a few examples.
In[260]:= Timing[
pointInPolygon[{-86, 60}, segmentbins, xmin, xmax, ymin, ymax]]
Out[260]= {0.005, False}
In[261]:= Timing[
pointInPolygon[{-70, 55}, segmentbins, xmin, xmax, ymin, ymax]]
Out[261]= {0.003999, True}
In[262]:= Timing[
pointInPolygon[{-137, 64}, segmentbins, xmin, xmax, ymin, ymax]]
Out[262]= {0., True}
In[263]:= Timing[
pointInPolygon[{-122, 61}, segmentbins, xmin, xmax, ymin, ymax]]
Out[263]= {0.001, True}
In[266]:= Timing[
pointInPolygon[{-73, 70}, segmentbins, xmin, xmax, ymin, ymax]]
Out[266]= {0.002999, True}
Daniel Lichtblau
Wolfram Research