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