Path: news.mathworks.com!not-for-mail
From: <HIDDEN>
Newsgroups: comp.soft-sys.matlab
Subject: Re: Extract points (on a sphere) inside polygon. HELP!!!
Date: Fri, 26 Aug 2011 19:42:11 +0000 (UTC)
Organization: The MathWorks, Inc.
Lines: 72
Message-ID: <j38suj$c2b$1@newscl01ah.mathworks.com>
References: <5525763.1178619749040.JavaMail.jakarta@nitrogen.mathforum.org> <ellieandrogerxyzzy-0805071153420001@dialup-4.232.228.18.dial1.losangeles1.level3.net> <j36pk6$i50$1@newscl01ah.mathworks.com> <j37dji$d8l$1@newscl01ah.mathworks.com> <j38jke$a1g$1@newscl01ah.mathworks.com>
Reply-To: <HIDDEN>
NNTP-Posting-Host: www-04-blr.mathworks.com
Content-Type: text/plain; charset=UTF-8; format=flowed
Content-Transfer-Encoding: 8bit
X-Trace: newscl01ah.mathworks.com 1314387731 12363 172.30.248.35 (26 Aug 2011 19:42:11 GMT)
X-Complaints-To: news@mathworks.com
NNTP-Posting-Date: Fri, 26 Aug 2011 19:42:11 +0000 (UTC)
X-Newsreader: MATLAB Central Newsreader 1187260
Xref: news.mathworks.com comp.soft-sys.matlab:741386

"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