From: <HIDDEN>
Newsgroups: comp.soft-sys.matlab
Subject: Re: Extract points (on a sphere) inside polygon. HELP!!!
Date: Tue, 23 Aug 2011 14:12:14 +0000 (UTC)
Organization: The MathWorks, Inc.
Lines: 87
Message-ID: <j30cfu$bu3$>
References: <> <>
Reply-To: <HIDDEN>
Content-Type: text/plain; charset=UTF-8; format=flowed
Content-Transfer-Encoding: 8bit
X-Trace: 1314108734 12227 (23 Aug 2011 14:12:14 GMT)
NNTP-Posting-Date: Tue, 23 Aug 2011 14:12:14 +0000 (UTC)
X-Newsreader: MATLAB Central Newsreader 2996139
Xref: comp.soft-sys.matlab:740969 (Roger Stafford) wrote in message <>...
> In article
> <>, General
> <> wrote:
> > Hi everyone,
> > 
> > I have created a sphere in Matlab and drawn some points on it. I am
> trying to extract some points directly from the sphere by drawing a
> polygon around those points. Unfortunately I cannot understand how to
> this. Can anybody help me?
> > 
> > Many Thanks
> > Regards
> > Alex
> --------------------
>   You haven't stated what coordinate system you are using to locate your
> various points and the vertices of the "polygon".  Whatever it is,
> presumably you can convert it to Cartesian coordinates with origin at the
> center of the sphere.  I will assume that the polygon vertices are
> connected by great circle arcs (face angles) which are less than than pi
> in arc angle.  I will assume moreover that the polygon is "convex" in the
> sense that all interior spherical angles (dihedral angles) of the polygon
> are less than pi.
>   Let P1 = (x1,y1,z1) and P2 = (x2,y2,z2) be two successive vertices of
> the polygon and ordered so that the polygon interior lies to the left of
> the vector pointing from P1 to P2.  If P = (x,y,z) is an arbitrary point,
> then
>  dot(P,cross(P1,P2)) >= 0
> if P lies on the side of the plane through P1, P2, and the origin, which
> contains the interior of the polygon.  So this gives you the algorithm. 
> Test a point P = (x,y,z) to see if it satisfies the above inequality for
> each pair of adjacent vertices in the polygon assuming they are ordered in
> the above "right hand" sense.
>   That can be done in the following way.  Let A be an n x 3 array in which
> successive rows give the x, y, and z coordinates of successive vertices
> ordered as above.  Define a shifted version of A:
>  B = A([2:n,1],:);
> Then the row vector P = [x,y,z] represents a point (x,y,z) in the interior
> of the polygon if and only if:
>  all(dot(P,cross(A,B,2),2)) >= 0
> is true.
> Roger Stafford

Hi Roger
I wanted to find out if a point is inside a polygon on earth or not and tried to use your amazing technique 

I computed 


if(dot_product(P,p1) >= 0 .AND.  &
dot_product(P,p2) >= 0 .AND.  &
dot_product(P,p3) >= 0 .AND.  &
dot_product(P,p4) >= 0 ) then



Here edge_sp_1 ,edge_sp_2,edge_sp_3,edge_sp_4 are actually verticies of the polygon that have been converted  to x y z using . P1,P2,P3,P4 are points arranged in counter clockwise direction

p(1) = R*cos(p_lon)*cos(p_lat)
p(2) = R*sin(p_lon)*cos(p_lat)
p(3) = R*sin(p_lat)

P is the point to be tested. 

It does not work always. Am i missing something.

To be true I give it a set of points that I know are inside/outside the polygon. It shows that some of the points that are detected inside are actually outside the polygon.