Reputation: 11
I have some MultiPolygon of some cities of a country. For a given point, I want to check in which city it lies. (We can assume that every given point is located in this country thus it lies in exactly one city.)
One approach is to iterate each polygon and check if the point lies in that polygon by Point-in-Polygon algorithms. But because I may have many points, and Point-in-Polygon algorithms have at least O(n) complexity, it's not an appropriate way.
So I want an appropriate algorithm for this problem with these assumptions:
What algorithm can be used for this problem?
Upvotes: 1
Views: 852
Reputation: 2448
Point-in-Polygon seems necessary to me. But to speed up your algorithm, partition your space with a grid. You can have a quad that encopasses the city and a quad that is inside the city, to speed up your detection. Only if the point is inside the outer quad and not inside the inner quad use the point-in-poly algo. Quad Trees as suggested by Mehdi can be implemented in very fast ways. But the best solution is probably mapping to a grid as explained below, you will get a O(1) solution:
You can probably speed up the PIP algorithm tesselating the polygon too.
With an appropriate rounding function, you could use the coords of your point as indices for a bi-dimensional array and determine in O(1) if it's in a particular city, if it's not in a city or if it needs further investigation.
Example: point is (3.4560, 5.1547) Grid is 10x10. City1 is (3,4), (3,5) border of City1 crosses (3,3), (4,3), (4,4), (4,5), (4,6), (3,6), (2,6), (2,5), (2,4), (2,3). Rounding the point to integers you get 3,5 that can be used as indices for matrix that will have in Mat(3,5)=C1, where C1 stands for City1. You can use this trick with geographical coords with an arbitrary chosen resolution.
formula to map a coordinate to an index of the array:
Coordinate c in [a,b]
array index i in [0,n-1]
i = floor( (c-a) / ((b-a)/n) )
Upvotes: 2
Reputation:
You can get a significant speedup by working on a projection of the geometry on the vertical axis. Every polygon will degenerate to a segment, which you can insert in a segment tree https://en.wikipedia.org/wiki/Segment_tree. Then with query time O(log(N)+K) you will report the K polygons crossed by a straight line through the point.
If the polygons aren't too large, you can expect a reduction of the number of polygons to be tried from N to √N (roughly).
It may be advantageous to perform this process on the X axis as well and only process the polygons common to both results lists. This way, you will only test for inclusion the polygons that have their bounding box around the test point.
Upvotes: 1