[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Re: determining if a point is "inside" or "outside" a shape

<dmarshall@ivory.trentu.ca> wrote in message
news:FJt6r8.GrL.A.ebony@news.trentu.ca...
> I have an interesting problem that I'm betting someone has solved already
> in some form or another.
>
> I have a set of x,y coordinates that describe a polygon, 2d in my case.
> How do I tell if another x,y point is inside the boundary of my polygon?

I wrote function POLY_INSIDE (below & attached) some time ago to do this. I
think it worked...

---
Mark Hadfield
m.hadfield@niwa.cri.nz  http://katipo.niwa.cri.nz/~hadfield/
National Institute for Water and Atmospheric Research
PO Box 14-901, Wellington, New Zealand


---------------------------------------------------------
;+
; NAME:
;   POLY_INSIDE
;
; PURPOSE:
;   This function determines if points are inside/outside a polygon.
;
; CATEGORY:
;   Geography, Geometry.
;
; CALLING SEQUENCE:
;   Result = POLY_INSIDE(X, Y, XP, YP)
;
; INPUTS:
;   X,Y:        X, Y position(s) defining the point(s) to be tested. Can
;               be vectors.
;
;   XP,YP:      Vectors of X, Y positions defining the polygon.
;
; INPUT KEYWORDS:
;   EDGE:       Set this keyword to accept edge (& vertex) points as
;               "inside". Default is to reject them. Note that this feature
;               is still not implemented properly.
;
; OUTPUTS:
;   The function returns 0 if the point is outside the polygon, 1 if
;   it is inside the polygon.
;
; PROCEDURE:
;   This routine uses the fact that the sum of angles between
;   lines from a point to successive vertices of a polygon is 0
;   for a point outside the polygon, and +/- 2*pi for a point
;   inside. A point on an edge will have one such successive angle
;   equal to +/- pi.
;
; EXAMPLE:
;
; MODIFICATION HISTORY:
;   Mark Hadfield, June 1995:
;       Written based on ideas in MATLAB routine INSIDE.M in the WHOI
;       Oceanography Toolbox  v1.4 (R. Pawlowicz, 14 Mar 94,
;       rich@boreas.whoi.edu).
;-
function POLY_INSIDE, X, Y, XP, YP, EDGE_INSIDE=edge

    ; Check geometry of the array of points. Note that this array
    ; is treated as a 1-D array internally to allow matrix operations,
    ; but the array structure and shape is restored on output.
    n = n_elements(X)  &  s = size(x)

    ; Make a local copy of polygon data
    ; If necessary, add a last point to close the polygon
    xpp = xp(*)  &  ypp = yp(*)  &  npp = n_elements(xpp)
    if xpp(npp-1) ne xpp(0) or ypp(npp-1) ne ypp(0) then begin
        xpp = [xpp,xpp(0)]  &  ypp = [ypp,ypp(0)]  &  npp = npp+1
    endif

    case npp of

        1:  if keyword_set(edge) then $
                return, fix(x eq xpp(0) and y eq ypp(0)) $
            else $
                return, reproduce(0,x)

     else:  begin

                ; The following matrix operations return (npp,n) arrays
containing
                ; distances from points to vertices
                dx = xpp#replicate(1D0,n) - replicate(1D0,npp)#x(*)
                dy = ypp#replicate(1D0,n) - replicate(1D0,npp)#y(*)

                angles = (atan(dy,dx))
                angles = angles(1:npp-1,*)-angles(0:npp-2,*)
                oor = where(angles le -!dpi,count)
                if count gt 0 then angles(oor) = angles(oor) + 2.*!dpi
                oor = where(angles gt !dpi,count)
                if count gt 0 then angles(oor) = angles(oor) - 2.*!dpi

                return,round(total(angles/!dpi,1)) ne 0


            end

    endcase

end




begin 666 poly_inside.pro
M.RL-"CL@3D%-13H-"CL@("!03TQ97TE.4TE$10T*.PT*.R!055)03U-%.@T*
M.R @(%1H:7,@9G5N8W1I;VX@9&5T97)M:6YE<R!I9B!P;VEN=',@87)E(&EN
M<VED92]O=71S:61E(&$@<&]L>6=O;BX-"CL-"CL@0T%414=/4EDZ#0H[(" @
M1V5O9W)A<&AY+"!'96]M971R>2X-"CL-"CL@0T%,3$E.1R!315%514Y#13H-
M"CL@("!297-U;'0@/2!03TQ97TE.4TE$12A8+"!9+"!84"P@65 I#0H[#0H[
M($E.4%544SH-"CL@("!8+%DZ(" @(" @("!8+"!9('!O<VET:6]N*',I(&1E
M9FEN:6YG('1H92!P;VEN="AS*2!T;R!B92!T97-T960N($-A;@T*.R @(" @
M(" @(" @(" @(&)E('9E8W1O<G,N#0H[#0H[(" @6% L65 Z(" @(" @5F5C
M=&]R<R!O9B!8+"!9('!O<VET:6]N<R!D969I;FEN9R!T:&4@<&]L>6=O;BX-
M"CL-"CL@24Y0550@2T595T]21%,Z#0H[(" @141'13H@(" @(" @4V5T('1H
M:7,@:V5Y=V]R9"!T;R!A8V-E<'0@961G92 H)B!V97)T97@I('!O:6YT<R!A
M<PT*.R @(" @(" @(" @(" @(")I;G-I9&4B+B!$969A=6QT(&ES('1O(')E
M:F5C="!T:&5M+B!.;W1E('1H870@=&AI<R!F96%T=7)E#0H[(" @(" @(" @
M(" @(" @:7,@<W1I;&P@;F]T(&EM<&QE;65N=&5D('!R;W!E<FQY+@T*.PT*
M.R!/55105513.@T*.R @(%1H92!F=6YC=&EO;B!R971U<FYS(# @:68@=&AE
M('!O:6YT(&ES(&]U='-I9&4@=&AE('!O;'EG;VXL(#$@:68-"CL@("!I="!I
M<R!I;G-I9&4@=&AE('!O;'EG;VXN#0H[#0H[(%!23T-%1%5213H-"CL@("!4
M:&ES(')O=71I;F4@=7-E<R!T:&4@9F%C="!T:&%T('1H92!S=6T@;V8@86YG
M;&5S(&)E='=E96X-"CL@("!L:6YE<R!F<F]M(&$@<&]I;G0@=&\@<W5C8V5S
M<VEV92!V97)T:6-E<R!O9B!A('!O;'EG;VX@:7,@, T*.R @(&9O<B!A('!O
M:6YT(&]U='-I9&4@=&AE('!O;'EG;VXL(&%N9" K+RT@,BIP:2!F;W(@82!P
M;VEN= T*.R @(&EN<VED92X@02!P;VEN="!O;B!A;B!E9&=E('=I;&P@:&%V
M92!O;F4@<W5C:"!S=6-C97-I=F4@86YG;&4-"CL@("!E<75A;"!T;R K+RT@
M<&DN#0H[#0H[($5804U03$4Z#0H[#0H[($U/1$E&24-!5$E/3B!(25-43U)9
M.@T*.R @($UA<FL@2&%D9FEE;&0L($IU;F4@,3DY-3H-"CL@(" @(" @5W)I
M='1E;B!B87-E9"!O;B!I9&5A<R!I;B!-051,04(@<F]U=&EN92!)3E-)1$4N
M32!I;B!T:&4@5TA/20T*.R @(" @("!/8V5A;F]G<F%P:'D@5&]O;&)O>" @
M=C$N-" H4BX@4&%W;&]W:6-Z+" Q-"!-87(@.30L#0H[(" @(" @(')I8VA 
M8F]R96%S+G=H;VDN961U*2X-"CLM#0IF=6YC=&EO;B!03TQ97TE.4TE$12P@
M6"P@62P@6% L(%E0+"!%1$=%7TE.4TE$13UE9&=E#0H-"B @(" [($-H96-K
M(&=E;VUE=')Y(&]F('1H92!A<G)A>2!O9B!P;VEN=',N($YO=&4@=&AA="!T
M:&ES(&%R<F%Y#0H@(" @.R!I<R!T<F5A=&5D(&%S(&$@,2U$(&%R<F%Y(&EN
M=&5R;F%L;'D@=&\@86QL;W<@;6%T<FEX(&]P97)A=&EO;G,L#0H@(" @.R!B
M=70@=&AE(&%R<F%Y('-T<G5C='5R92!A;F0@<VAA<&4@:7,@<F5S=&]R960@
M;VX@;W5T<'5T+@T*(" @(&X@/2!N7V5L96UE;G1S*%@I(" F("!S(#T@<VEZ
M92AX*0T*#0H@(" @.R!-86ME(&$@;&]C86P@8V]P>2!O9B!P;VQY9V]N(&1A
M=&$-"B @(" [($EF(&YE8V5S<V%R>2P@861D(&$@;&%S="!P;VEN="!T;R!C
M;&]S92!T:&4@<&]L>6=O;@T*(" @('AP<" ]('AP*"HI(" F("!Y<' @/2!Y
M<"@J*2 @)B @;G!P(#T@;E]E;&5M96YT<RAX<' I#0H@(" @:68@>'!P*&YP
M<"TQ*2!N92!X<' H,"D@;W(@>7!P*&YP<"TQ*2!N92!Y<' H,"D@=&AE;B!B
M96=I;@T*(" @(" @("!X<' @/2!;>'!P+'AP<"@P*5T@("8@('EP<" ](%MY
M<' L>7!P*# I72 @)B @;G!P(#T@;G!P*S$-"B @("!E;F1I9@T*#0H@(" @
M8V%S92!N<' @;V8-"@T*(" @(" @(" Q.B @:68@:V5Y=V]R9%]S970H961G
M92D@=&AE;B D#0H@(" @(" @(" @(" @(" @<F5T=7)N+"!F:7@H>"!E<2!X
M<' H,"D@86YD('D@97$@>7!P*# I*2 D#0H@(" @(" @(" @("!E;'-E("0-
M"B @(" @(" @(" @(" @("!R971U<FXL(')E<')O9'5C92@P+'@I#0H-"B @
M(" @96QS93H@(&)E9VEN#0H-"B @(" @(" @(" @(" @(" [(%1H92!F;VQL
M;W=I;F<@;6%T<FEX(&]P97)A=&EO;G,@<F5T=7)N("AN<' L;BD@87)R87ES
M(&-O;G1A:6YI;F<-"B @(" @(" @(" @(" @(" [(&1I<W1A;F-E<R!F<F]M
M('!O:6YT<R!T;R!V97)T:6-E<PT*(" @(" @(" @(" @(" @(&1X(#T@>'!P
M(W)E<&QI8V%T92@Q1# L;BD@+2!R97!L:6-A=&4H,40P+&YP<"DC>"@J*0T*
M(" @(" @(" @(" @(" @(&1Y(#T@>7!P(W)E<&QI8V%T92@Q1# L;BD@+2!R
M97!L:6-A=&4H,40P+&YP<"DC>2@J*0T*#0H@(" @(" @(" @(" @(" @86YG
M;&5S(#T@*&%T86XH9'DL9'@I*0T*(" @(" @(" @(" @(" @(&%N9VQE<R ]
M(&%N9VQE<R@Q.FYP<"TQ+"HI+6%N9VQE<R@P.FYP<"TR+"HI#0H@(" @(" @
M(" @(" @(" @;V]R(#T@=VAE<F4H86YG;&5S(&QE("TA9'!I+&-O=6YT*0T*
M(" @(" @(" @(" @(" @(&EF(&-O=6YT(&=T(# @=&AE;B!A;F=L97,H;V]R
M*2 ](&%N9VQE<RAO;W(I("L@,BXJ(61P:0T*(" @(" @(" @(" @(" @(&]O
M<B ]('=H97)E*&%N9VQE<R!G=" A9'!I+&-O=6YT*0T*(" @(" @(" @(" @
M(" @(&EF(&-O=6YT(&=T(# @=&AE;B!A;F=L97,H;V]R*2 ](&%N9VQE<RAO
M;W(I("T@,BXJ(61P:0T*#0H@(" @(" @(" @(" @(" @<F5T=7)N+')O=6YD
M*'1O=&%L*&%N9VQE<R\A9'!I+#$I*2!N92 P#0H-"@T*(" @(" @(" @(" @
=96YD#0H-"B @("!E;F1C87-E#0H-"F5N9 T*#0H`
`
end