| point inside polygon [message #11519] |
Tue, 31 March 1998 00:00  |
Philippe Peeters
Messages: 14
Registered: April 1996
|
Junior Member |
|
|
Hello,
Does anybody knows of an IDL function to test whether a given point is
inside a polygon? The coordinates of the polygon are real values,
positive and negative.
--
Philippe Peeters
------------------------------------------------------------ ------------
Belgian Institute for Space Aeronomy Tel: +32-2-373.03.81
Institut d'Aeronomie Spatiale de Belgique Fax: +32-2-374.84.23
3 Avenue Circulaire Email: philp@oma.be
B-1180 Brussels, Belgium http://www.oma.be/BIRA-IASB/
|
|
|
|
| Re: point inside polygon [message #19177 is a reply to message #11519] |
Mon, 06 March 2000 00:00  |
davidf
Messages: 2866
Registered: September 1996
|
Senior Member |
|
|
Klaus Scipal (kscipal@ipf.tuwien.ac.at) writes:
> the algorithm works fine until the point is lying on the polygon. in such a
> case the algorithm returns sometimes in and sometimes out. what i am looking
> for is an algorithm which discerns if a point is lying on the polygon and
> give only one answer either in or out but not one time in and the next time
> out.
Humm. Well, here is another one by Julie Greenwood and used
with her permission. I haven't used this one myself, but it
might be worth a try. It uses a nice coding style, which
always suggests to me that it will probably work. :-)
Cheers,
David
function IsPointInPolygon, xPts, yPts, XTest, YTest
;+
; NAME:
; IsPointInPolygon
;
; PURPOSE:
; Determine whether point (XTest,YTest) is within the boundary of
; the polygon defined by vectors xPts & yPts.
; Function returns 1 if point is within polygon, else returns 0
; Use ArePointsInPolygon to test entire array at a time
;
; AUTHOR: Julie Greenwood (julieg@oceanweather.com)
;
; CATEGORY: Geometry
;
; CALLING SEQUENCE:
; bResult = IsPointInPolygon (xPts, yPts, XTest, YTest)
;
; INPUTS:
; xPts, yPts - vectors containing vertices of convex polygon in $
; counterclockwise order. See argument B to routine Triangulate.pro
; xTest, yTest - point to test for within-ness
;
; Optional Inputs: None
;
; OUTPUTS:
; Function returns -1 if point is within polygon, else returns 0
;
; OPTIONAL OUTPUTS: None
;
; KEYWORD Parameters: None
;
; COMMON BLOCKS: None
;
; SIDE EFFECTS: None
;
; RESTRICTIONS:
; Polygon must be closed (first point is same as last)
;
; PROCEDURE:
; See: http://www.swin.edu.au/astronomy/pbourke/geometry/insidepoly /
;
; Consider a polygon made up of N vertices (xi,yi) where I ranges from
; 0 to N-1. The last vertex (xN,yN) is assumed to be the same as the
; first vertex (x0,y0), that is, the polygon is closed. To determine
; the status of a point (xp,yp) consider a horizontal ray emanating
; from (xp,yp) and to the right. If the number of times this ray
; intersects the line segments making up the polygon is even then the
; point is outside the polygon. Whereas if the number of intersections
; is odd then the point (xp,yp) lies inside the polygon.
;
; MODIFICATION HISTORY:
; jgg - 2-Dec-1999 - Created
;-
npts = n_elements(xPts)
if (npts ne n_elements(yPts)) then begin
print,' Error in IsPointInPolygon: vertex arrays must have same size.'
return,PointIsIn
endif
if (npts lt 2) then begin
print,' Error in IsPointInPolygon: polygon has less than 2 points.'
return,PointIsIn
endif
PointIsIn = 0
for ia = 0,nPts-1 do begin
ib = (ia+1) mod nPts
sLine = xPts[ia] + $
(yTest-yPts[ia]) * (xPts[ib]-xPts[ia]) / (yPts[ib]-yPts[ia])
bInRange = $
(( (yPts[ia] le yTest) and (yTest lt yPts[ib]) ) or $
( (yPts[ib] le yTest) and (yTest lt yPts[ia]) ))
if ( bInRange and (xTest lt sLine) ) then begin
PointIsIn = not PointIsIn
endif
endfor
return, PointIsIn
end
function ArePointsInPolygon, xPts, yPts, XTest, YTest
;+
; NAME:
; ArePointsInPolygon
;
; PURPOSE:
; Determine whether points (XTest,YTest) are within the boundary of
; the polygon defined by vectors xPts & yPts.
; Function returns array of same size and shape as XTest, containing
; 1 if point is within polygon, else 0
; Use IsPointInPolygon to test one point at a time
;
; AUTHOR: Julie Greenwood (julieg@oceanweather.com)
;
; CATEGORY: Geometry
;
; CALLING SEQUENCE:
; bResult = ArePointsInPolygon (xPts, yPts, XTest, YTest)
;
; INPUTS:
; xPts, yPts - vectors containing vertices of convex polygon in $
; counterclockwise order. See argument B to routine Triangulate.pro
; xTest, yTest - arrays of points to test for within-ness
;
; Optional Inputs: None
;
; OUTPUTS:
; Function returns array of same size and shape as XTest, containing
; -1 if point is within polygon, else 0
;
; OPTIONAL OUTPUTS: None
;
; KEYWORD Parameters: None
;
; COMMON BLOCKS: None
;
; SIDE EFFECTS: None
;
; RESTRICTIONS:
; Polygon must be closed (first point is same as last)
;
; PROCEDURE:
; See: http://www.swin.edu.au/astronomy/pbourke/geometry/insidepoly /
;
; Consider a polygon made up of N vertices (xi,yi) where I ranges from
; 0 to N-1. The last vertex (xN,yN) is assumed to be the same as the
; first vertex (x0,y0), that is, the polygon is closed. To determine
; the status of a point (xp,yp) consider a horizontal ray emanating
; from (xp,yp) and to the right. If the number of times this ray
; intersects the line segments making up the polygon is even then the
; point is outside the polygon. Whereas if the number of intersections
; is odd then the point (xp,yp) lies inside the polygon.
;
; MODIFICATION HISTORY:
; jgg - 3-Dec-1999 - Created
;-
nsizex = size(xTest)
nsizey = size(xTest)
if (total(nsizex ne nsizey)) then begin
help, xTest, yTest
print,' Error in ArePointsInPolygon: grid arrays must have same size.'
return,PointsIn
endif
if (nsizex[0] ne 2) then begin
help, xTest, yTest
print,' Error in ArePointsInPolygon: grid arrays must have 2 dimensions'
return,PointsIn
endif
npts = n_elements(xPts)
if (npts ne n_elements(yPts)) then begin
help, xPts, yPts
print,' Error in ArePointsInPolygon: vertex arrays must have same size.'
return,PointsIn
endif
if (npts lt 2) then begin
help, xPts, yPts
print,' Error in ArePointsInPolygon: polygon has less than 2 points.'
return,PointsIn
endif
PointsIn = LonArr(nsizex[1],nsizex[2])
for ia = 0,nPts-1 do begin
ib = (ia+1) mod nPts
sLine = xPts[ia] + $
(yTest-yPts[ia]) * (xPts[ib]-xPts[ia]) / (yPts[ib]-yPts[ia])
bInRange = $
(( (yPts[ia] le yTest) and (yTest lt yPts[ib]) ) or $
( (yPts[ib] le yTest) and (yTest lt yPts[ia]) ))
GotIt = where ( bInRange and (xTest lt sLine) , count)
if (count ne 0) then PointsIn[GotIt] = not PointsIn[GotIt]
endfor
return, PointsIn
end
--
David Fanning, Ph.D.
Fanning Software Consulting
Phone: 970-221-0438 E-Mail: davidf@dfanning.com
Coyote's Guide to IDL Programming: http://www.dfanning.com/
Toll-Free IDL Book Orders: 1-888-461-0155
|
|
|
|
|
|
|
|
| Re: point inside polygon [message #19356 is a reply to message #19177] |
Wed, 08 March 2000 00:00 |
Randall Frank
Messages: 17
Registered: October 1999
|
Junior Member |
|
|
Hello,
In IDL 5.3, you can use the IDLanROI object:
o=obj_new('IDLanROI',[[0,0],[0,1],[1,1],[1,0]])
print,o->containspoints([[0,0],[0,0.5],[0.5,0.5],[5,0.5]])
obj_destroy,o
This checks to see if the points [0,0] [0,0.5] [0.5,0.5] and [5,0.5]
lie inside the unit square polygon from the origin. The output
is:
3 2 1 0
for the points being on a vertex, on an edge, inside, or outside
the polygon. This should cover all the cases you were asking
about and gives you access to a bunch of other stats (see
IDLanROI::ComputeGeometry()).
Hope it helps...
David Fanning wrote:
>
> Klaus Scipal (kscipal@ipf.tuwien.ac.at) writes:
>
>> the algorithm works fine until the point is lying on the polygon. in such a
>> case the algorithm returns sometimes in and sometimes out. what i am looking
>> for is an algorithm which discerns if a point is lying on the polygon and
>> give only one answer either in or out but not one time in and the next time
>> out.
>
> Humm. Well, here is another one by Julie Greenwood and used
> with her permission. I haven't used this one myself, but it
> might be worth a try. It uses a nice coding style, which
> always suggests to me that it will probably work. :-)
>
> Cheers,
>
> David
>
> function IsPointInPolygon, xPts, yPts, XTest, YTest
> ;+
> ; NAME:
> ; IsPointInPolygon
> ;
> ; PURPOSE:
> ; Determine whether point (XTest,YTest) is within the boundary of
> ; the polygon defined by vectors xPts & yPts.
> ; Function returns 1 if point is within polygon, else returns 0
> ; Use ArePointsInPolygon to test entire array at a time
> ;
> ; AUTHOR: Julie Greenwood (julieg@oceanweather.com)
> ;
> ; CATEGORY: Geometry
> ;
> ; CALLING SEQUENCE:
> ; bResult = IsPointInPolygon (xPts, yPts, XTest, YTest)
> ;
> ; INPUTS:
> ; xPts, yPts - vectors containing vertices of convex polygon in $
> ; counterclockwise order. See argument B to routine Triangulate.pro
> ; xTest, yTest - point to test for within-ness
> ;
> ; Optional Inputs: None
> ;
> ; OUTPUTS:
> ; Function returns -1 if point is within polygon, else returns 0
> ;
> ; OPTIONAL OUTPUTS: None
> ;
> ; KEYWORD Parameters: None
> ;
> ; COMMON BLOCKS: None
> ;
> ; SIDE EFFECTS: None
> ;
> ; RESTRICTIONS:
> ; Polygon must be closed (first point is same as last)
> ;
> ; PROCEDURE:
> ; See: http://www.swin.edu.au/astronomy/pbourke/geometry/insidepoly /
> ;
> ; Consider a polygon made up of N vertices (xi,yi) where I ranges from
> ; 0 to N-1. The last vertex (xN,yN) is assumed to be the same as the
> ; first vertex (x0,y0), that is, the polygon is closed. To determine
> ; the status of a point (xp,yp) consider a horizontal ray emanating
> ; from (xp,yp) and to the right. If the number of times this ray
> ; intersects the line segments making up the polygon is even then the
> ; point is outside the polygon. Whereas if the number of intersections
> ; is odd then the point (xp,yp) lies inside the polygon.
> ;
> ; MODIFICATION HISTORY:
> ; jgg - 2-Dec-1999 - Created
> ;-
>
> npts = n_elements(xPts)
> if (npts ne n_elements(yPts)) then begin
> print,' Error in IsPointInPolygon: vertex arrays must have same size.'
> return,PointIsIn
> endif
> if (npts lt 2) then begin
> print,' Error in IsPointInPolygon: polygon has less than 2 points.'
> return,PointIsIn
> endif
>
> PointIsIn = 0
> for ia = 0,nPts-1 do begin
> ib = (ia+1) mod nPts
>
> sLine = xPts[ia] + $
> (yTest-yPts[ia]) * (xPts[ib]-xPts[ia]) / (yPts[ib]-yPts[ia])
>
> bInRange = $
> (( (yPts[ia] le yTest) and (yTest lt yPts[ib]) ) or $
> ( (yPts[ib] le yTest) and (yTest lt yPts[ia]) ))
>
> if ( bInRange and (xTest lt sLine) ) then begin
>
> PointIsIn = not PointIsIn
>
> endif
> endfor
>
> return, PointIsIn
>
> end
>
> function ArePointsInPolygon, xPts, yPts, XTest, YTest
> ;+
> ; NAME:
> ; ArePointsInPolygon
> ;
> ; PURPOSE:
> ; Determine whether points (XTest,YTest) are within the boundary of
> ; the polygon defined by vectors xPts & yPts.
> ; Function returns array of same size and shape as XTest, containing
> ; 1 if point is within polygon, else 0
> ; Use IsPointInPolygon to test one point at a time
> ;
> ; AUTHOR: Julie Greenwood (julieg@oceanweather.com)
> ;
> ; CATEGORY: Geometry
> ;
> ; CALLING SEQUENCE:
> ; bResult = ArePointsInPolygon (xPts, yPts, XTest, YTest)
> ;
> ; INPUTS:
> ; xPts, yPts - vectors containing vertices of convex polygon in $
> ; counterclockwise order. See argument B to routine Triangulate.pro
> ; xTest, yTest - arrays of points to test for within-ness
> ;
> ; Optional Inputs: None
> ;
> ; OUTPUTS:
> ; Function returns array of same size and shape as XTest, containing
> ; -1 if point is within polygon, else 0
> ;
> ; OPTIONAL OUTPUTS: None
> ;
> ; KEYWORD Parameters: None
> ;
> ; COMMON BLOCKS: None
> ;
> ; SIDE EFFECTS: None
> ;
> ; RESTRICTIONS:
> ; Polygon must be closed (first point is same as last)
> ;
> ; PROCEDURE:
> ; See: http://www.swin.edu.au/astronomy/pbourke/geometry/insidepoly /
> ;
> ; Consider a polygon made up of N vertices (xi,yi) where I ranges from
> ; 0 to N-1. The last vertex (xN,yN) is assumed to be the same as the
> ; first vertex (x0,y0), that is, the polygon is closed. To determine
> ; the status of a point (xp,yp) consider a horizontal ray emanating
> ; from (xp,yp) and to the right. If the number of times this ray
> ; intersects the line segments making up the polygon is even then the
> ; point is outside the polygon. Whereas if the number of intersections
> ; is odd then the point (xp,yp) lies inside the polygon.
> ;
> ; MODIFICATION HISTORY:
> ; jgg - 3-Dec-1999 - Created
> ;-
>
> nsizex = size(xTest)
> nsizey = size(xTest)
> if (total(nsizex ne nsizey)) then begin
> help, xTest, yTest
> print,' Error in ArePointsInPolygon: grid arrays must have same size.'
> return,PointsIn
> endif
> if (nsizex[0] ne 2) then begin
> help, xTest, yTest
> print,' Error in ArePointsInPolygon: grid arrays must have 2 dimensions'
> return,PointsIn
> endif
>
> npts = n_elements(xPts)
> if (npts ne n_elements(yPts)) then begin
> help, xPts, yPts
> print,' Error in ArePointsInPolygon: vertex arrays must have same size.'
> return,PointsIn
> endif
> if (npts lt 2) then begin
> help, xPts, yPts
> print,' Error in ArePointsInPolygon: polygon has less than 2 points.'
> return,PointsIn
> endif
>
> PointsIn = LonArr(nsizex[1],nsizex[2])
> for ia = 0,nPts-1 do begin
> ib = (ia+1) mod nPts
>
> sLine = xPts[ia] + $
> (yTest-yPts[ia]) * (xPts[ib]-xPts[ia]) / (yPts[ib]-yPts[ia])
>
> bInRange = $
> (( (yPts[ia] le yTest) and (yTest lt yPts[ib]) ) or $
> ( (yPts[ib] le yTest) and (yTest lt yPts[ia]) ))
>
> GotIt = where ( bInRange and (xTest lt sLine) , count)
> if (count ne 0) then PointsIn[GotIt] = not PointsIn[GotIt]
>
> endfor
>
> return, PointsIn
>
> end
>
> --
> David Fanning, Ph.D.
> Fanning Software Consulting
> Phone: 970-221-0438 E-Mail: davidf@dfanning.com
> Coyote's Guide to IDL Programming: http://www.dfanning.com/
> Toll-Free IDL Book Orders: 1-888-461-0155
--
rjf.
Randy Frank | ASCI Visualization
Lawrence Livermore National Laboratory | rjfrank@llnl.gov
B451 Room 2039 L-561 | Voice: (925) 423-9399
Livermore, CA 94550 | Fax: (925) 423-8704
|
|
|
|
comp.lang.idl-pvwave archive
Members
Search
Help
Login
Home
