CLIM mail archive

# [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.
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

```