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
please find below what I use for his purpose. I hope it is useful.
Best regards,
Frank
PointInPolygonQ::usage="PointInPolygonQ[pt, poly] uses the winding-
number algorithm (Godkin and Pulli, 1984) to check, if point pt is
inside the closed polygon poly, which is given as list of its
vertices."
(*
checks, if a point is inside a polygon
pt: point as {lat [deg], lon [deg]} to test
poly: list of polygon vertices coordinates
GODKIN,C.B. AND J.J.PULLI: APPLICATION OF THE "WINDING-NUMBER
ALGORITHM" TO THE SPATIAL SORTING OF CATALOGED EARTHQUAKE DATA.
Bull. Seismol. Soc. Am. 74, 5, PP. 1845-1848, OCTOBER 1984
RETURN VALUE: 0 IF POINT OUTSIDE
+/-1 IF POINT INSIDE
2 IF POINT IS ON AN EDGE OR VERTEX
*)
PointInPolygonQ[pt_,poly_] := Module[
{ i,n,isicr,inside,px,py,pxx,pyy,x0,y0 },
n = Length[poly];
(* ACCUMULATE SIGNED CROSSING NUMBERS WITH INSIDE *)
inside = 0;
{x0,y0}=pt;
For[i=1,i < n,i++,
(* PROCEED AROUND POLYGON CHECKING EACH SEGMENT TO SEE IF
NEGATIVE X-AXIS WAS CROSSED
TRANSLATE COORDINATES OF POLYGON TO PUT TEST POINT AT
ORIGIN *)
{px,py} = poly[[i]];
{pxx,pyy} = poly[[i+1]];
isicr = ksicr[px - x0, py - y0, pxx - x0, pyy - y0];
(* STOP COUNTING IF POINT IS ON EDGE *)
If[isicr == 4, Return[2]];
inside += isicr;
];
(* CHECK SEGMENT FROM LAST VERTEX TO FIRST VERTEX *)
{px,py} = poly[[n]];
{pxx,pyy} = poly[[1]];
isicr = ksicr[px - x0, py - y0, pxx - x0, pyy - y0];
If[isicr == 4, Return[2]];
inside = (inside + isicr)/2;
Return[inside];
];
(*
COMPUTE SIGNED CROSSING NUMBER
A "SIGNED CROSSING NUMBER" TELLS WETHER A SEGMENT
(IN THIS CASE THE SEGMENT FROM (X1,Y1) TO (X2,Y2))
CROSSES THE NEGATIVE X-AXIS OR GOES THROUGH THE ORIGIN
THE RETURN VALUES ARE:
+2 IF SEGMENT CROSSES FROM BELOW
+1 IF SEGMENT EITHER ENDS ON -X-AXIS FROM BELOW OR STARTS
UPWARDS FROM -X-AXIS ("HALF CROSSING")
0 IF THERE IS NO CROSSING
-1 IF SEGMENT EITHER ENDS ON -X-AXIS FROM ABOVE OR STARTS
DOWNWARDS FROM -X-AXIS ("HALF CROSSING")
-2 IF SEGMENT CROSSES FROM ABOVE
+4 IF SEGMENT CROSSES THROUGH THE ORIGIN
*)
ksicr[x1_,y1_,x2_,y2_] := Module[
{
},
(* IF BOTH POINTS ARE ON THE SAME SIDE OF X-AXIS, RETURN 0 *)
If[N[y1*y2 > 0.], Return[0] (* no crossing *)];
(* CHECK IF SEGMENT CROSSES THROUGH THE ORIGIN *)
If[x1*y2 != x2*y1 || x1*x2 > 0.,
If[y1 * y2 < 0,
(*
COMPLETE CROSSING OF -X-AXIS?
BREAK INTO CASES ACCORDING TO CROSSING DIRECTION
*)
If[y1 > 0,
(* CASE Y1 > 0 > Y2 *)
If[y1 * x2 >= x1 * y2, Return[0];, (* no crossing *)
Return[-2]; (* downward crossing *)
];
,
(* CASE Y1 < 0 < Y2 *)
If[x1 * y2 >= y1 * x2, Return[0];, (* no crossing *)
Return[2]; (* upward crossing *)
];
];
,
(*
HALF CROSSING?
ONE END OF SEGMENT TOUCHES X-AXIS! WHICH END?
*)
If[y2 == 0,
(* HERE Y2==0 CHECK IF SEGMENT TOUCHES +X-AXIS *)
If[y1 == 0 || x2 > 0,
Return[0]; (* no crossing *)
,
(* UPWARD OR DOWNWARD? *)
If[y1 > 0.,
Return[-1]; (* Downward half crossing *)
,
Return[1]; (* Upward half crossing *)
];
];
,
(* HERE Y1==0 CHECK IF SEGMENT TOUCHES +X-AXIS *)
If[x1 > 0,
Return[0];
,
(* UPWARD OR DOWNWARD? *)
If[y2 > 0,
Return[1]; (* Upward half crossing *)
,
Return[-1]; (* Downward half crossing *)
];
];
];
(* HERE Y1=0 CHECK IF SEGMENT TOUCHES +X-AXIS *)
If[x1 > 0,
Return[0]; (* no crossing *)
];
(* UPWARD OR DOWNWARD? *)
If[y2 > 0.,
Return[-1]; (* Downward half crossing *),
Return[1]; (* Upward half crossing *)
];
];
,
Return[4];
];
];
the web is full of this test as well as any graphics
book. Look at
http://local.wasp.uwa.edu.au/~pbourke/geometry/insidepoly/
http://tog.acm.org/editors/erich/ptinpoly/
Regards
Jens
Here is my attempt. The idea is to move around the polygon and add up the
angles made with the point. If the point is outside this will be zero. The
steps are basically:
1) Check if the point is one of the vertices of the polygon.
2) Add the first polygon point to the end of the list of points.
3) Subtract the test point from each of the polygon vertices.
4) Partition the points into pairs.
5) Rotate each pair so the first one lies along the x axis. (To prevent
crossing the branch line)
6) Find the ArcTan of the second point. This will have a sign.
7) Sum the angles and take the absolute value.
8) Test if this is zero or some multiple of 2\[Pi].
PointInsidePolygonQ::usage =
"PointInsidePolygonQ[point,polygon] will return True if the point \
is on the boundary or inside the polygon and False otherwise.";
PointInsidePolygonQ[point_, polygon_] :=
Module[{work},
work = If[Head[polygon] === Polygon, First[polygon], polygon];
If[ \[Not] FreeQ[work, point], Return[True]];
work = # - point & /@ Join[work, {First[work]}];
work = Partition[work, 2, 1];
work = RotationTransform[{First[#], {1, 0}}]@# & /@ work;
work = Abs@Total[ArcTan @@ Last[#] & /@ work] // Chop // Rationalize;
TrueQ[work != 0]
]
Here is a graphical test for a simple polygon:
testpoints = testpoints = RandomReal[{-9, 9}, {100, 2}];
polypoints = {{-1.856, 3.293}, {1.257,
2.635}, {2.395, -0.6587}, {-1.018, -2.455}, {-3.293, -0.05988}};
Graphics[
{Lighter[Green, .8], Polygon[polypoints],
AbsolutePointSize[4],
{If[PointInsidePolygonQ[#, polypoints], Black, Red], Point[#]} & /@
testpoints},
PlotRange -> 10,
Frame -> True]
Here is a test for a more complex polygon.
testpoints = testpoints = RandomReal[{-9, 9}, {100, 2}];
polypoints = {{-3.653, 5.329}, {0.2395, 6.168}, {-0.8982,
1.138}, {-0.6587, 1.138}, {5.569, 3.234}, {6.527, -2.036}, {1.677,
0.479}, {-6.407, -1.976}, {-5.21, 2.635}, {1.856, -3.713}};
Graphics[
{Lighter[Green, .8], Polygon[polypoints],
AbsolutePointSize[4],
{If[PointInsidePolygonQ[#, Polygon[polypoints]], Black, Red],
Point[#]} & /@ testpoints},
PlotRange -> 10,
Frame -> True]
There might be simpler methods.
David Park
djm...@comcast.net
http://home.comcast.net/~djmpark/
> is there a way to test whether a point is inside a polygon? ie.
>
> PointInsidePolygonQ[point_,polygon_] -> True or False
[snip]
You could have a look at the different algorithms presented in Paul
Bourke's paper [1], pick up one of your liking, then code it from the
provided C source code to Mathematica syntax. How the algorithms work,
their benefits and draw backs, as well as tested C code and diagrams,
are neatly presented.
Depending of your level of familiarity with Mathematica programming, it
might be easier to translate the C syntax into equivalent Mathematica
procedural constructs, test the code, and then convert it into fast
functional Mathematica code.
If you have performance issues or troubles converting procedural code
(fast in C, but slow in Mathematica) into functional code (fast in
Mathematica), people in MathGroup will help you to improve performances
or overcome coding issues.
Regards,
--Jean-Marc
[1] Paul Bourke, "Determining If A Point Lies On The Interior Of A
Polygon",
http://local.wasp.uwa.edu.au/~pbourke/geometry/insidepoly/
Hi Mitch,
To write a routines that checks if a point is inside a polygon you must
realize that if you walk around the polygon you see a inner point
always on the same side. This is false for a outside point. Therefore,
the test is to check if the point in question: p0 is always on the same
side of the polygon sides.
Assume {p1,p1,..} are the polygon points, arrange in a definite order
(clockwise or anti-clockwise). How can we check on which side of {pi+1,
pi} p0 is? We may get a vector perpendicular to the side by rotating by
Pi/2: {{0,1},{-1,0}}.(pi+1 - pi). If we take the dot product with (p0
-pi) the sign will tell us the side. Here is an example:
pts = {{0, 0}, {1, 0}, {1, 1}, {0, 1}};
p0 = {0.5, 0.5};
sides= Subtract @@ # & /@ Partition[Append[pts, pts[[1]]], 2, 1];
perp= {{0,1},{-1,0}}.#& /@ sides;
Equal @@ Sign @ MapThread[Dot, {p0 - # & /@ pts, perp}]
The third lines calculates (pi+1 - pi), note that we added to first
point to the end of the list to close the polygon. The fourth turns each
vector by Pi/2. The fifth calculates p0-pi = {p0 - # & /@ pts and takes
the Dot product of corresponding vectors. Finally it takes the sign and
checks if all signs are the same.
hope this helps, Daniel
PointInPoly[{x_, y_}, poly_List] :=
Module[{i, j, c = False, npol = Length[poly]},
For[i = 1; j = npol, i <= npol, j = i++,
If[((((poly[[i, 2]] <= y) && (y <
poly[[j, 2]])) || ((poly[[j, 2]] <= y) && (y <
poly[[i, 2]])))
&& (x < (poly[[j, 1]] -
poly[[i, 1]])*(y - poly[[i, 2]])/(poly[[j, 2]] -
poly[[i, 2]]) + poly[[i, 1]])), c = \[Not] c
];
];
Return [c];
]
Cheers -- Sjoerd
On Feb 9, 12:31 pm, Mitch Murphy <mi...@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=
Cheers -- Sjoerd
On Feb 10, 12:50 pm, "David Park" <djmp...@comcast.net> wrote:
> Mitch,
>
> Here is my attempt. The idea is to move around the polygon and add up the
> angles made with the point. If the point is outside this will be zero. Th=
e
> steps are basically:
>
> 1) Check if the point is one of the vertices of the polygon.
> 2) Add the first polygon point to the end of the list of points.
> 3) Subtract the test point from each of the polygon vertices.
> 4) Partition the points into pairs.
> 5) Rotate each pair so the first one lies along the x axis. (To prevent
> crossing the branch line)
> 6) Find the ArcTan of the second point. This will have a sign.
> 7) Sum the angles and take the absolute value.
> 8) Test if this is zero or some multiple of 2\[Pi].
>
> PointInsidePolygonQ::usage =
> "PointInsidePolygonQ[point,polygon] will return True if the point \
> is on the boundary or inside the polygon and False otherwise.";
>
> PointInsidePolygonQ[point_, polygon_] :=
> Module[{work},
> work = If[Head[polygon] === Polygon, First[polygon], polygon]=
;
> If[ \[Not] FreeQ[work, point], Return[True]];
> work = # - point & /@ Join[work, {First[work]}];
> work = Partition[work, 2, 1];
> work = RotationTransform[{First[#], {1, 0}}]@# & /@ work;
> work = Abs@Total[ArcTan @@ Last[#] & /@ work] // Chop // Rationaliz=
e;
> TrueQ[work != 0]
> ]
>
> Here is a graphical test for a simple polygon:
>
> testpoints = testpoints = RandomReal[{-9, 9}, {100, 2}];
> polypoints = {{-1.856, 3.293}, {1.257,
> 2.635}, {2.395, -0.6587}, {-1.018, -2.455}, {-3.293, -0.05988}};
> Graphics[
> {Lighter[Green, .8], Polygon[polypoints],
> AbsolutePointSize[4],
> {If[PointInsidePolygonQ[#, polypoints], Black, Red], Point[#]} & /@
> testpoints},
> PlotRange -> 10,
> Frame -> True]
>
> Here is a test for a more complex polygon.
>
> testpoints = testpoints = RandomReal[{-9, 9}, {100, 2}];
> polypoints = {{-3.653, 5.329}, {0.2395, 6.168}, {-0.8982,
> 1.138}, {-0.6587, 1.138}, {5.569, 3.234}, {6.527, -2.036}, {1.677=
,
> 0.479}, {-6.407, -1.976}, {-5.21, 2.635}, {1.856, -3.713}};
> Graphics[
> {Lighter[Green, .8], Polygon[polypoints],
> AbsolutePointSize[4],
> {If[PointInsidePolygonQ[#, Polygon[polypoints]], Black, Red],
> Point[#]} & /@ testpoints},
> PlotRange -> 10,
> Frame -> True]
>
> There might be simpler methods.
>
> David Park
> djmp...@comcast.nethttp://home.comcast.net/~djmpark/
the problem was discussed by the mathgroup in 2000 and Luc Barthelet has
collected some results at
http://www.mathematica-users.org/webMathematica/wiki/wiki.jsp?pageName=Notebook:PointInsidePolygon.nb
Adriano Pascoletti
2009/2/9 Mitch Murphy <mi...@lemma.ca>
The routine that you submitted was much faster than the one I submitted.
That was basically because I was using a RotationTransform to obtain the
proper signed rotations for each side of the polygon. However, this can be
done much faster using complex arithmetic. So here is a new routine:
angle[{p1_, p2_}] :=
Module[{c1, c2},
{c1, c2} = #.{1, I} & /@ {p1, p2};
Arg[c2/c1]]
PointInsidePolygonQ::usage =
"PointInsidePolygonQ[point,polygon] will return True if the point \
is on the boundary or inside the polygon and False otherwise.";
SyntaxInformation[
PointInsidePolygonQ] = {"ArgumentsPattern" -> {_, _}};
PointInsidePolygonQ[point_, polygon_] :=
Module[{work = Join[polygon, {First[polygon]}]},
If[ \[Not] FreeQ[work, point], Return[True]];
work = # - point & /@ work;
(Total[angle /@ Partition[work, 2, 1]] // Chop) != 0
]
Here are graphical test routines for a simple polygon and one that folds on
itself.
testpoints = RandomReal[{-9, 9}, {5000, 2}];
polypoints = {{-1.856, 3.293}, {1.257,
2.635}, {2.395, -0.6587}, {-1.018, -2.455}, {-3.293, -0.05988}};
Graphics[
{Lighter[Green, .8], Polygon[polypoints],
AbsolutePointSize[2],
{If[PointInsidePolygonQ[#, polypoints], Black, Red], Point[#]} & /@
testpoints},
PlotRange -> 10,
Frame -> True,
ImageSize -> 400] // Timing
testpoints = testpoints = RandomReal[{-9, 9}, {5000, 2}];
polypoints = {{-3.653, 5.329}, {0.2395, 6.168}, {-0.8982,
1.138}, {-0.6587, 1.138}, {5.569, 3.234}, {6.527, -2.036}, {1.677,
0.479}, {-6.407, -1.976}, {-5.21, 2.635}, {1.856, -3.713}};
Graphics[
{Lighter[Green, .8], Polygon[polypoints],
AbsolutePointSize[2],
{If[PointInsidePolygonQ[#, polypoints], Black, Red], Point[#]} & /@
testpoints},
PlotRange -> 10,
Frame -> True,
ImageSize -> 400] // Timing
The following are the comparable graphical routines for the Godkin, Pulli
algorithm below.
testpoints = testpoints = RandomReal[{-9, 9}, {5000, 2}];
polypoints = {{-1.856, 3.293}, {1.257,
2.635}, {2.395, -0.6587}, {-1.018, -2.455}, {-3.293, -0.05988}};
Graphics[
{Lighter[Green, .8], Polygon[polypoints],
AbsolutePointSize[2],
{If[Abs@PointInPolygonQ[#, polypoints] == 1, Black, Red],
Point[#]} & /@ testpoints},
PlotRange -> 10,
Frame -> True,
ImageSize -> 400] // Timing
testpoints = testpoints = RandomReal[{-9, 9}, {5000, 2}];
polypoints = {{-3.653, 5.329}, {0.2395, 6.168}, {-0.8982,
1.138}, {-0.6587, 1.138}, {5.569, 3.234}, {6.527, -2.036}, {1.677,
0.479}, {-6.407, -1.976}, {-5.21, 2.635}, {1.856, -3.713}};
Graphics[
{Lighter[Green, .8], Polygon[polypoints],
AbsolutePointSize[2],
{If[Abs@PointInPolygonQ[#, polypoints] == 1, Black, Red],
Point[#]} & /@ testpoints},
PlotRange -> 10,
Frame -> True,
ImageSize -> 400] // Timing
The timings appear to be comparable, although the Godkin, Pulli routine
might be slightly faster.
From: Frank Scherbaum [mailto:Frank.S...@geo.uni-potsdam.de]
Mitch,
please find below what I use for his purpose. I hope it is useful.
Best regards,
Frank
Am Feb 9, 2009 um 11:32 AM schrieb Mitch Murphy:
>
> 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
>
PointInPolygonQ::usage="PointInPolygonQ[pt, poly] uses the winding-
Needs["VectorAnalysis`"]
SetCoordinates[Spherical];
mean = {1, t0, p0};
pt = {1, t, p};
xy = CoordinatesToCartesian /@ {mean, pt};
rot = RotationMatrix@xy;
LeafCount@rot
32798
I don't know what to tell Simplify about this, but it seem there are MANY
unnecessary Conjugate mentions:
Cases[rot, Conjugate[z_], Infinity] // Length
1575
and MANY unnecessary cases like this, too:
Cases[rot, Power[Abs[z_], 2], Infinity] // Length
1980
Each of these could be simply z^2, I think.
Trying to eliminate these with rules doesn't seem to help much, if any.
The first thing I thought of was PowerExpand, but
Map[PowerExpand, rot, {2}] // LeafCount
32798
What am I missing? How can it be that complicated in the first place?
Bobby
Needs["VectorAnalysis`"]
SetCoordinates[Spherical];
center = {1, 0, 0};
pt = {1, t, p};
xy = CoordinatesToCartesian /@ {center, pt};
rot = RotationMatrix@xy;
LeafCount@rot
1059
rot //. {Conjugate -> Identity, Abs[z_]^2 :> z^2} // Simplify
% // LeafCount
{{Cos[p]^2 (Cos[t] + Tan[p]^2), -Sin[2 p] Sin[t/2]^2,
Cos[p] Sin[t]}, {-Sin[2 p] Sin[t/2]^2, Cos[p]^2 + Cos[t] Sin[p]^2,
Sin[p] Sin[t]}, {-Cos[p] Sin[t], -Sin[p] Sin[t], Cos[t]}}
80
Bobby
On Wed, 11 Feb 2009 04:23:42 -0600, DrMajorBob <btr...@austin.rr.com>
wrote:
{
{"Daniel", 0.828, 1.219},
{"Daniel/Sjoerd", 0.688, 1.078},
{"Sedgewick/Sjoerd", 0.688, 0.954},
{"David1", 41.766, 81.89},
{"David2", 1.578, 2.531},
{"Godkin/Frank", 1.485, 2.359}
}
In the "Daniel/Sjoerd" method I have used Daniel's code and taken out
every intermediate step. This saves an additional 20%. The code is now
completely gibberish, but it still works ;-)
PointInsidePolygonQ[pt_, poly_] :=
Module[{},
Equal @@
Sign@MapThread[
Dot, {pt - # & /@
poly, {{0, 1}, {-1, 0}}.# & /@ (Subtract @@ # & /@
Partition[Append[poly, poly[[1]]], 2, 1])}]
]
Cheers --Sjoerd
On Feb 9, 12:31 pm, Mitch Murphy <mi...@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=
1) None of these solutions will work for a country that straddles the
Greenwich meridian or a continent that straddles a pole. Antarctica does
both.
CountryData["Antarctica", "Shape"]
To handle this in general, we'd need to rotate the globe so that the
country/continent in question doesn't have that problem... and test points
must be on the same hemisphere, too, to be inside a country.
(This is why I was asking about RotationMatrix.)
2) The poster mentioned Canada, which is made up of disjoint polygons:
Dimensions /@ First@canada;
Length@%
%%[[All, 1]] // Total
25
30458
Twenty-five polygons, more than 30K points... and that's the BASIC
version. Here's the other:
canada = CountryData["Canada", "FullPolygon"];
Dimensions /@ First@canada;
Length@%
%%[[All, 1]] // Total
971
43262
971 polygons, 43 thousand points! For the OP's application, he may need
VERY fast code.
3) In the PointInsidePolygon code, I'd use MemberQ and Append, not FreeQ
and Join, and I NEVER, EVER use Return. (But I found no difference in
speed.)
pointInsidePolygonQ[point_, polygon_] /; MemberQ[polygon, point] = True;
pointInsidePolygonQ[point_, polygon_] :=
0 != Chop@
Total[angle /@
Partition[# - point & /@ Append[polygon, First@polygon], 2, 1]]
4) Points on edges of the polygon test as OUTSIDE:
polypoints = {{-1.856, 3.293}, {1.257,
2.635}, {2.395, -0.6587}, {-1.018, -2.455}, {-3.293, -0.05988}};
PointInsidePolygonQ[(
2 polypoints[[1]] + 8 polypoints[[2]])/10, polypoints]
False
One solution for this would be to stop calculating angles if one of them
happens to be Pi. For instance,
Clear[angle, pointInsidePolygonQ]
angle[{p1_, p2_}] :=
Module[{c1, c2, a}, {c1, c2} = #.{1, I} & /@ {p1, p2};
a = Arg[c2/c1];
Chop[Pi - a] == 0 && Throw[True];
a
]
pointInsidePolygonQ[point_, polygon_] /; MemberQ[polygon, point] =
True;
pointInsidePolygonQ[point_, polygon_] :=
Catch[0 !=
Chop@Total[
angle /@
Partition[# - point & /@ Append[polygon, First@polygon], 2, 1]]]
polypoints = {{-1.856, 3.293}, {1.257,
2.635}, {2.395, -0.6587}, {-1.018, -2.455}, {-3.293, -0.05988}};
pointInsidePolygonQ[(
2 polypoints[[1]] + 8 polypoints[[2]])/10, polypoints]
True
This is slower for the test data, but it will be FASTER if many test
points are on the polygon.
Bobby
>> 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
>>
>
Oops, I meant the meridian 180 degrees from Greenwich.
Bobby
On Wed, 11 Feb 2009 17:38:29 -0600, DrMajorBob <btr...@austin.rr.com>
wrote:
Of course, we can write our own program to do just about
anything we want. If you want to know whether there is
a built in function "InPolygon", the answer is no. If
you want to if there are freely available packages that
include this functionality, the answer is yes:
CompG:
http://library.wolfram.com/infocenter/MathSource/5712/
and IMTEK:
http://www.imtek.de/simulation//mathematica/IMSweb/
Curiously, the functionality is clearly "in Mathematica",
it's just not exposed. You can see that a point in
polygon is implemented by playing with the following:
Graphics[Tooltip[Polygon[
{{0, 0}, {2, 0}, {2, 1}, {1, 1}, {1, 2}, {0, 2}}],
"in polygon"]]
In fact, there must be a large computational geometry
library contained in the Mathematica kernel to support
the extensive graphics routines. Hopefully, that
functionality will be exposed in the not too distant
future.
Mark McClure
Andrzej Kozlowski
On 12 Feb 2009, at 11:38, Sjoerd C. de Vries wrote:
> Since timing is often important for functions like this one, I've
> compared the various algorithms posted here using David's two tests.
> The results are below:
>
> {
> {"Daniel", 0.828, 1.219},
> {"Daniel/Sjoerd", 0.688, 1.078},
> {"Sedgewick/Sjoerd", 0.688, 0.954},
> {"David1", 41.766, 81.89},
> {"David2", 1.578, 2.531},
> {"Godkin/Frank", 1.485, 2.359}
> }
>
> In the "Daniel/Sjoerd" method I have used Daniel's code and taken out
> every intermediate step. This saves an additional 20%. The code is now
> completely gibberish, but it still works ;-)
>
> PointInsidePolygonQ[pt_, poly_] :=
> Module[{},
> Equal @@
> Sign@MapThread[
> Dot, {pt - # & /@
> poly, {{0, 1}, {-1, 0}}.# & /@ (Subtract @@ # & /@
> Partition[Append[poly, poly[[1]]], 2, 1])}]
> ]
>
> Cheers --Sjoerd
>
> On Feb 9, 12:31 pm, Mitch Murphy <mi...@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=
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
I missed a couple of useful optimizations. We can use Compile both in
preprocessing and in the predicate query itself.
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[
getSegsC[j, minx, len, eps, flatsegments]
, {j, nbins}];
{{minx, maxx}, {miny, maxy}, segmentbins}
]
getSegsC = Compile[
{{j, _Integer}, {minx, _Real}, {len, _Real}, {eps, _Real}, {segs,
_Real, 3}},
Module[{lo, hi},
lo = minx + (j - 1)*len - eps;
hi = minx + j*len + eps;
Select[segs,
Module[{xlo, xhi}, {xlo, xhi} = Sort[{#[[1, 1]], #[[2, 1]]}];
lo <= xlo <= hi || lo <= xhi <= hi || (xlo <= lo && xhi >=
= hi)
] &]]];
With this we can preprocess the polygons for Canada, using 1000 bins,
in around a minute on my machine.
In[346]:= canpoly = First[CountryData["Canada", "Polygon"]];
In[347]:= nbins = 1000;
Timing[{{xmin, xmax}, {ymin, ymax}, segmentbins} =
polyToSegmentList[canpoly, nbins];]
Out[337]= {55.3256, Null}
To repeat from the last note, there are almost certainly smarter ways
to do the preprocessing, so as to make it faster. For most
applications I can think of that would be more trouble than it is
worth, so it's not something I have attempted.
For the predicate evaluation we can do:
pointInPolygon[{x_, y_}, bins_, xmin_, xmax_, ymin_, ymax_] :=
Catch[Module[
{nbins = Length[bins], bin},
If[x < xmin || x > xmax || y < ymin || y > ymax, Throw[False]];
bin = Ceiling[nbins*(x - xmin)/(xmax - xmin)];
If[EvenQ[countIntersectionsC[bins[[bin]], x, y, ymin - 1.]], False,
True]
]]
countIntersectionsC = Compile[
{{segs, _Real, 3}, {x, _Real}, {yhi, _Real}, {ylo, _Real}},
Module[{tally = 0, yval, xlo, xhi, y1, y2},
Do[
{{xlo, y1}, {xhi, y2}} = segs[[j]];
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
]];
With this we can now process around 10,000 points in a second. I
selected the region so that most would be inside; this was so that we
would not gain speed due to a high percentage failing the basic
rectangle test.
In[352]:= n = 10000;
pts = Transpose[{RandomReal[{-115, -55}, {n}],
RandomReal[{45, 75}, {n}]}];
In[354]:= Timing[
inout = Map[pointInPolygon[#, segmentbins, xmin, xmax, ymin, ymax] &,
pts];]
Out[354]= {1.04284, Null}
In[355]:= Take[inout, 20]
Out[355]= {True, True, True, True, False, True, True, True, True, \
True, False, True, True, False, False, False, False, True, True, True}
I should make a few remarks about the speed/complexity. If we split
the n polygon segments into m bins, for m somewhat smaller than n,
then for "good" geographies we expect something like O(n/m). Figure
most segments appear in no more than two bins, most appear in only one
bin, and no bin has more than, say, three times as many segments as
any other.
To achieve better algorithmic complexity I believe one would need to
use something more complicated, along the lines of a triangulation.
With such a data structure, efficiently implemented, it might be only O
(log(n)) to see whether a new point falls into a triangle, and, if so,
whether that triangle is inside or outside the boundary. I doubt such
a tactic would be applicable for a polygonal set of the size in
question here.
If the goal is, say, to color all pixels inside Canada differently
than those outside (in effect, to query every pixel in the bounding
rectangle), then it would be faster simply to find all boundary points
for a given vertical. This transforms from a problem of query on
points to one of splitting segments (a query on lines, so to speak).
Much more efficient, I would guess.
Daniel Lichtblau
Wolfram Research
> >> 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["Cana=
da=