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?
>
>Robb,
>
> 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?
>
>Rich
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
FITNESS FOR A PARTICULAR PURPOSE.
**************************************************************
A RANDOMIZED PLANAR POINT LOCATION STRUCTURE
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.
Remarks:
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