CLIM mail archive
[Prev][Next][Index][Thread]
[Jeff Morrill: Re: point-in-polygon ]
OK, here's my contribution to the point-in-polygon algorithm. The issue
comes up every decade or so, as you can see from the references. I had
lost track of my copy, but JMorrill was able to find it. Thanks, Jeff.
k
------- Forwarded Message
To: Ken Anderson <KAnderson@BBN.COM>
Subject: Re: point-in-polygon
In-reply-to: Ken Anderson's message of Thu, 13 Jan 94 14:58:24 -0500.
Reply-To: jmorrill@BBN.COM
Date: Thu, 13 Jan 94 15:45:41 -0500
Message-ID: <5752.758493941@bbn.com>
From: Jeff Morrill <jmorrill@BBN.COM>
To: JMorrill@BBN.COM
Cc: Kanderson@BBN.COM, love@alcoa.com, spindler@alcoa.com,
ABoulanger@BBN.COM
From: Ken Anderson <KAnderson@BBN.COM>
Subject: point-in-polygon
Date: Thu, 13 Jan 94 14:58:24 -0500
Sender: kanderso@BBN.COM
I know i wrote fortran code a long time ago to decide if a point is in a
polygon, and i thought i had this in lisp, do you guys recall seening
anything like this from me?
k
This is probably what you were looking for. It was in scigraph.
jeff
(defun point-in-polygon (x y map-polygon polygon)
"Returns T if point X,Y is inside a closed boundary. Points on the
boundary edge are ambiguously either inside or outside. The
algorithm counts the number of line segements that intersect the
semi-infinite line from (x,y) to (+infinity,y). If there is an odd
number of such segements, the point is inside the polygon.
The function MAP-POLYGON maps a function over consecutive line
segments (X1 Y1 X2 Y2) of the polygon, such that the
first point is the same as the last. (The polygon need not be convex.)
See:
Shimrat, M., 1962, Algorithm 112: position of point relative to
polygon, CACM,5,434.
Hall, J.H., 1975, PTLOC - A FORTRAN subroutine for determining the
position of a point relative to a closed boundary, Math. Geol., 7,
75-79.
Anderson, K.R., 1975, Letter to the editor of Math. Geol.
The final version of this algorithm followed from discussions with
David Fitterman of the US Geological Survey, Jan 1975."
(let ((point-in NIL))
(funcall map-polygon
#'(lambda (x1 y1 x2 y2)
(if (eq (<= y y1) (> y y2)) ; Segement crosses ray.
(when (and (not (= y1 y2)) ; Ignore horizontal segement.
(< (- x x1 (/ (* (- y y1) (- x2 x1))
(- y2 y1))) ; Point is to left.
0))
;; (print (list x1 y1 x2 y2))
(setq point-in (not point-in)))))
polygon)
point-in))
------- End of Forwarded Message
Main Index |
Thread Index