[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: point inside polygon
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