Thread Subject:

Extract points (on a sphere) inside polygon. HELP!!!

In article
<5525763.1178619749040.JavaMail.jakarta@nitrogen.mathforum.org>, General
<mammasis82@hotmail.com> 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

ellieandrogerxyzzy@mindspring.com.invalid (Roger Stafford) wrote in message <ellieandrogerxyzzy-0805071153420001@dialup-4.232.228.18.dial1.losangeles1.level3.net>...
> In article
> <5525763.1178619749040.JavaMail.jakarta@nitrogen.mathforum.org>, General
> <mammasis82@hotmail.com> 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
http://www.mathworks.com.au/matlabcentral/newsreader/view_thread/147669

I computed


p1=cross_product(edge_sp_4,edge_sp_3)
p2=cross_product(edge_sp_3,edge_sp_2)
p3=cross_product(edge_sp_2,edge_sp_1)
p4=cross_product(edge_sp_1,edge_sp_4)

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

checks=.true.

endif


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.

"glanplon Bali" wrote in message <j30cfu$bu3$1@newscl01ah.mathworks.com>...
> 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
> http://www.mathworks.com.au/matlabcentral/newsreader/view_thread/147669
>
> I computed
>
> p1=cross_product(edge_sp_4,edge_sp_3)
> p2=cross_product(edge_sp_3,edge_sp_2)
> p3=cross_product(edge_sp_2,edge_sp_1)
> p4=cross_product(edge_sp_1,edge_sp_4)
>
> 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
>
> checks=.true.
>
> endif
>
> 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.
- - - - - - - - -
  I cannot tell what the difficulty is without some numerical data to work with. Perhaps you could give me a specific numerical example, hopefully a short one, that doesn't work for you.

  Your cross product expressions look backwards to me, but I can't really tell without the data. Remember, I said "... (if) the polygon interior lies to the left of the vector pointing from P1 to P2 .... then dot(P,cross(P1,P2)) >= 0". You seem to use cross(P2,P1) which points in the opposite direction.

  When you transformed longitude and latitude to cartesian coordinates I assume you used radian measure for them. Otherwise the transformation equations you describe would be very much in error, at least in matlab.

  Also remember I specified that all face angles and dihedral angles in the spherical polygon should be less than pi.

  I notice in that original 2007 article I got the parentheses wrong in the expression

 all(dot(P,cross(A,B,2),2)) >= 0

It should have been

 all(dot(P,cross(A,B,2),2) >= 0)

However, you have apparently rewritten it correctly in C so that shouldn't be a source of trouble.

Roger Stafford

Hi Roger,
As an example the four verticies ( V1 V2 V3 V4, in clockwise direction) of a quadilateral ( in longitude, latitude format ) are
V1 108.674866 -49.93645
V2 108.67069 -49.945087
V3 108.68399 -49.947758
V4 108.68817 -49.93913
Point (P) to be tested is 108.67943 longitude -49.9421 latitude

When I calculate the


p1=cross_product(edge_sp_4,edge_sp_3)
p2=cross_product(edge_sp_3,edge_sp_2)
p3=cross_product(edge_sp_2,edge_sp_1)
p4=cross_product(edge_sp_1,edge_sp_4)

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

inpolygon=.true.

endif

inpolygon is true if all the above conditions are satisfied.

edge_sp_1/2/3/4 are the cartisiean coordinates of the verticies v1/2/3/4. V1,V2,V3,V4 are arranged from lower left to lower right in clockwise direction.
 
I know this point is inside the polygon but my test is showing as outside. Can you suggest what is wrong that I am doing
GlanPlon

"glanplon Bali" wrote in message <j36pk6$i50$1@newscl01ah.mathworks.com>...
> Hi Roger,
> As an example the four verticies ( V1 V2 V3 V4, in clockwise direction) of a quadilateral ( in longitude, latitude format ) are
> V1 108.674866 -49.93645
> V2 108.67069 -49.945087
> V3 108.68399 -49.947758
> V4 108.68817 -49.93913
> Point (P) to be tested is 108.67943 longitude -49.9421 latitude
>
> When I calculate the
>
>
> p1=cross_product(edge_sp_4,edge_sp_3)
> p2=cross_product(edge_sp_3,edge_sp_2)
> p3=cross_product(edge_sp_2,edge_sp_1)
> p4=cross_product(edge_sp_1,edge_sp_4)
>
> 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
>
> inpolygon=.true.
>
> endif
>
> inpolygon is true if all the above conditions are satisfied.
>
> edge_sp_1/2/3/4 are the cartisiean coordinates of the verticies v1/2/3/4. V1,V2,V3,V4 are arranged from lower left to lower right in clockwise direction.
>
> I know this point is inside the polygon but my test is showing as outside. Can you suggest what is wrong that I am doing
> GlanPlon
- - - - - - - - -
Hello GlanPlon,

  Using your values for V1 through V4 and P and an earth radius of 6378 kilometers, I get

P1 = [x1,y1,z1] = [-1314.4477 3888.9851 -4881.2812] (kilometers)
P2 = [x2,y2,z2] = [-1313.9286 3888.3837 -4881.9000]
P3 = [x3,y3,z3] = [-1314.7583 3887.8630 -4882.0913]
P4 = [x4,y4,z4] = [-1315.2775 3888.4635 -4881.4733]
P = [xP,yP,zP] = [-1314.6033 3888.4243 -4881.6860]

for their corresponding earth-centered cartesian coordinates and this makes all four inequalities true:

 dot(P,cross(P1,P2)) >= 0
 dot(P,cross(P2,P3)) >= 0
 dot(P,cross(P3,P4)) >= 0
 dot(P,cross(P4,P1)) >= 0

This therefore places P in the interior of your polygon.

  As I mentioned before it appears that you are taking your cross products backwards. My original statement was:

"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."

  (Note that it is really unnecessary to multiply by the earth radius as far as these inequalities are concerned. You could just as well assume a radius of one.)

Roger Stafford

Hi Roger,
Thanks for the reply. I have got it working now.
The thing that I changed is
p1=cross_product(v1,v2)
p2=cross_product(v2,v3)
p3=cross_product(v3,v4)
p4=cross_product(v4,v1)

But get some additional points inside the polygon which I know are outside the polygon. For example

v1(1)= 108.96776
v1(2)= -49.994953

v2(1)= 108.96363
v2(2)=-50.003593

v3(1)=108.976974
v3(2)= -50.006226

v4(1) = 108.98111
v4(2) = -49.99759

The point to be tested is
Point: lon: 108.96824 lat -50.00923

Can you suggest here. Can you also give the value of the cross products from your code.

Just wonder if you do not forget to convert the lon/lat to radian?

Why don't you post your code?

Bruno

Hi
Yes I have done the conversion and have replicated the values P1, P2, P3,P4 in xyz that Roger has got.
I have also noticed that when I take Radius of earth of 6378 it works and not with 6371.

Another example which I am not getting correct decision is

Verticies
V1 108.72806 -49.94714
V2 108.72389 -49.95578
V3 108.73719 -49.95845
V4 108.74136 -49.949818

Point to be tested
 108.71934, -49.950127

This point is outside the polygon but i get it inside.

You might play with it:

v1(1)= 108.96776;
v1(2)= -49.994953;

v2(1)= 108.96363;
v2(2)=-50.003593;

v3(1)=108.976974;
v3(2)= -50.006226;

v4(1) = 108.98111;
v4(2) = -49.99759;

% Test grid
[wlon wlat] = ndgrid(108.96:1e-3:108.99,-50.01:1e-3:-49.99);

v = [v1; v2; v3; v4]';
w = [wlon(:) wlat(:)]';
 
lonlat2phiteta = @(lonlat) lonlat*pi/180;
spherical2cart = @(phitheta) [cos(phitheta(1,:)).*cos(phitheta(2,:));
                              sin(phitheta(1,:)).*cos(phitheta(2,:));
                              sin(phitheta(2,:))]; % radius 1, no need to bother with earth radius
lonlat2cart = @(lonlat) lonlat2phiteta(spherical2cart(lonlat));

% Convert to cartesian
P = lonlat2cart(v);
Q = lonlat2cart(w);

% Engine based on Roger's formula
N = cross(P(:,1:end),P(:,[2:end 1]));
np = size(N,2);
nq = size(Q,2);
Q = repmat(reshape(Q, [3 nq 1]), [1 1 np]);
N = repmat(reshape(N, [3 1 np]), [1 nq 1]);
inside = all(dot(N,Q)>=0,3);
inside = reshape(inside,size(wlon));

% Check
clf
plot(v(1,[1:end 1]), v(2,[1:end 1]), '-b');
hold on
h1 = plot(w(1,inside), w(2,inside), 'ob')
h2 = plot(w(1,~inside), w(2,~inside), '.c')
% Your point
h3 = plot(108.96824, -50.00923, 'rs');
axis equal
legend([h1 h2 h3], 'inside', 'outside', 'your point where is it?');

% Bruno

"glanplon" wrote in message <j38pel$c7$1@newscl01ah.mathworks.com>...

> Yes I have done the conversion and have replicated the values P1, P2, P3,P4 in xyz that Roger has got.
> I have also noticed that when I take Radius of earth of 6378 it works and not with 6371.
>

POST YOUR CODE PLEASE!!!

Bruno

"glanplon" wrote in message <j38jke$a1g$1@newscl01ah.mathworks.com>...
> Hi Roger,
> Thanks for the reply. I have got it working now.
> The thing that I changed is
> p1=cross_product(v1,v2)
> p2=cross_product(v2,v3)
> p3=cross_product(v3,v4)
> p4=cross_product(v4,v1)
>
> But get some additional points inside the polygon which I know are outside the polygon. For example
>
> v1(1)= 108.96776
> v1(2)= -49.994953
>
> v2(1)= 108.96363
> v2(2)=-50.003593
>
> v3(1)=108.976974
> v3(2)= -50.006226
>
> v4(1) = 108.98111
> v4(2) = -49.99759
>
> The point to be tested is
> Point: lon: 108.96824 lat -50.00923
>
> Can you suggest here. Can you also give the value of the cross products from your code.
- - - - - - - - -
clear
format long
Lon1 = 108.96776; Lat1 = -49.994953;
Lon2 = 108.96363; Lat2 = -50.003593;
Lon3 = 108.976974; Lat3 = -50.006226;
Lon4 = 108.98111; Lat4 = -49.99759;
LonP = 108.96824; LatP = -50.00923;
R = 6378; %km
f = pi/180;
P1 = R*[cos(Lat1*f)*cos(Lon1*f),cos(Lat1*f)*sin(Lon1*f),sin(Lat1*f)];
P2 = R*[cos(Lat2*f)*cos(Lon2*f),cos(Lat2*f)*sin(Lon2*f),sin(Lat2*f)];
P3 = R*[cos(Lat3*f)*cos(Lon3*f),cos(Lat3*f)*sin(Lon3*f),sin(Lat3*f)];
P4 = R*[cos(Lat4*f)*cos(Lon4*f),cos(Lat4*f)*sin(Lon4*f),sin(Lat4*f)];
P = R*[cos(LatP*f)*cos(LonP*f),cos(LatP*f)*sin(LonP*f),sin(LatP*f)];
C12 = cross(P1,P2);
C23 = cross(P2,P3);
C34 = cross(P3,P4);
C41 = cross(P4,P1);
D12 = dot(P,C12);
D23 = dot(P,C23);
D34 = dot(P,C34);
D41 = dot(P,C41);

% Results:
P1 = -1332.69003600969 3877.49933824585 -4885.47030999304
P2 = -1332.17111150384 3876.89862059990 -4886.08853944187
P3 = -1333.00097394469 3876.37591899445 -4886.27691991586
P4 = -1333.52032733296 3876.97608571394 -4885.65901055242
P = -1332.32678157670 3876.33679256683 -4886.49183133320

C12 = -5331.97250273079 -3359.09849285427 -1211.55900679063
C23 = -3284.29832367599 3803.82633614354 3913.62053094152
C34 = 5327.82976611704 3361.37825764623 1213.18614611495
C41 = 3288.02100339159 -3804.88425323460 -3916.78749214485

D12 = 3205.876172551885
D23 = -3204.160734776407
D34 = 3199.735891658813
D41 = 9618.821957089007

Therefore the point P does not lie within the polygon. It is
on the wrong side of the P2 to P3 edge.

Roger Stafford

Hi Roger
Thanks for sending the output of the cross product. Problem solved it works. The problem in my case ( fortran not matlab) was the precision due to which the cross product was not properly computed.
Coming back to your algorithm. I was wondering if conversion to Cartesian is really necessary. Can I take only lat,lon coordinate system if the sphere is Earth and directly apply your algorithm. Let me know if that can work as well.
GlanPlon

"glanplon" wrote in message <j3ggjd$nfu$1@newscl01ah.mathworks.com>...
> Coming back to your algorithm. I was wondering if conversion to Cartesian is really necessary. Can I take only lat,lon coordinate system if the sphere is Earth and directly apply your algorithm. Let me know if that can work as well.
- - - - - - - - -
  For regions many kilometers in size, using linear combinations of latitude and longitude directly rather than of cartesian coordinates becomes increasingly inaccurate. My advice to you would be to stick with Cartesian coordinates except on very small areas which can be considered very flat.

  Also remember that a degree difference in longitude at the equator is in the neighborhood of 111 kilometers, whereas at your latitude of -50 degrees a degree difference is only about 72 kilometers and at the poles this approaches zero.

  You should realize however that even the assumption of a spherical earth is not completely accurate. More accurately it is an oblate spheroid, though in fact it may be ever so slightly "pear" shaped. Read the following for a discussion of this:

 http://en.wikipedia.org/wiki/Figure_of_the_Earth

Roger Stafford