point_in_polygon test

Roderick Brown r.brown at earthsci.unimelb.edu.au
Wed Nov 10 23:19:03 EST 1999

>On Wed, 10 Nov 1999, Robb Hill wrote:
>> Can anyone offer advice?  Any good books on these types of algorithms or
>> code?
>  That's an excellent question! I, too, would like to have such a reference
>available. The last I heard -- and this goes back a dozen years when I first
>worked with ARC/Info -- is that the algorithms are rather closely held by
>those (like Jack Dagerman) who developed them. After all, the
>point-in-polygon and similar problems are why ESRI is able to charge the
>extremely high prices they do. It's also the reason why there are more
>raster-based analytical GISes than vector-based analytical offerings.
>  Please let me know what you find out. Have you posted to GIS-L?

This excercise seems to have confounded numerous folks in numerous areas of
endeavour it seems...I've appended two possible solutions (headers from
code) that I came across while searching for an answer some time ago. The
first is a simple minimalist solution, the second is a much more
sophisticated solution (can handle an array of polygons for the p_in_poly
test)...but will need more work to customise. If anybody else would like
the code (which can be freely and publicly distributed) let me know and
I'll forward same.

Regards, Roderick

 *                                                                         *
 *   INPOLY.C                                                              *
 *                                                                         *
 *   Copyright (c) 1995-1996 Galacticomm, Inc.  Freeware source code.      *
 *                                                                         *
 *   Please feel free to use this source code for any purpose, commercial  *
 *   or otherwise, as long as you don't restrict anyone else's use of      *
 *   this source code.  Please give credit where credit is due.            *
 *                                                                         *
 *   Point-in-polygon algorithm, created especially for World-Wide Web     *
 *   servers to process image maps with mouse-clickable regions.           *
 *                                                                         *
 *   Home for this file:  http://www.gcomm.com/develop/inpoly.c            *
 *                                                                         *
 *                                       6/19/95 - Bob Stein & Craig Yap   *
 *                                       stein at gcomm.com                   *
 *                                       craig at cse.fau.edu                 *
 *                                                                         *


Copyright (c) 1995  Raimund Seidel
All rights reserved.

This program is distributed WITHOUT ANY WARRANTY;
without even the implied warranty of MERCHANTABILITY or



	by R. Seidel

	last edited:  90 Feb 10

This file contains several routines for dealing with a general version
of the so-called planar point location problem, which can be described
as follows:  one is to preprocess a set of non-vertical straight line
segments into a data structure so that for any given query point p one
can determine quickly which segment is the first one directly below p.
It is assumed that no two segments have interior points in common.

The method employed is a randomized one that I have so far not gotten
around to write down formally in a paper.  It has the property that
a search requires fewer than 4*ln(n) comparisons in expectation.  It has
O(nlogn) expected preprocessing time and needs O(n) expected space.
All expectations are over "coin flips" in the algorithm and not over
assumed input distributions.

The Test Program allows for two usages:  either the input is a
set of segments as described above and for query points the next
segment below are returned, or the input is the sequence of corners
counter clockwise around a simple polygon and for query points it
is determined whether they lie inside this polygon or not.

Just compile this file and run it.  It should be fairly self-explanatory.

Some warnings:  There is only minimal checking whether the input
	is correct (whether there are improper intersections).
	Although the program should not crash, the answers will not
	be meaningful.
	There is also no overflow checking.  Currently everything
	is done in terms of integer coordinates.  Keeping their
	magnitudes smaller than 2^15 will definitely avoid overflow.
	Using a coordinate larger than `infinity' can cause a crash.

   The coordinate type can be changed from integer to something
   else quite easily by changing the first 5 definitions.  A change
   to floating point will of course open the usual Pandora's box
   of precision problems.  However, they are refined to the two
   macros point_above_segment and point_below_segment.

   These two macros along with segment_above_segment are also
   the only things that would need to be adjusted if one wanted to
   use this algorithm for segments that are not straight, but
   still x-monotone.

For any remarks, complaints, etc. send email to

	seidel at lyra.berkeley.edu


 Roderick Brown                           Tele: (Int+ 61 3) 9344-9868
 School of Earth Sciences                Dept.: (Int+ 61 3) 9344-7675
 The University of Melbourne         Facsimile: (Int+ 61 3) 9344-7761
 Parkville, 3052                Home Telephone: (Int+ 61 3) 9397-8848
 Australia                     Email: r.brown at earthsci.unimelb.edu.au


More information about the grass-user mailing list