I have a set of non-intersecting lines, some of which are connected at vertices. I'm trying to find the smallest polygon, if one exists, that encloses a given point. So, in the image below, out of the list of all the line segments, given the point in red, I want to get the blue ones only. I'm using Python, but could probably adapt an algorithm from other languages; I have no idea what this problem is called.
First, remove all line segments that have at least one free endpoint, not coincident with any other segment. Do that repeatedly until no such segment remains.
Now you have a plane nicely subdivided into polygonal areas.
Find a segment closest to your point. Not the closest endpoint but the closest segment, mind you. Find out which direction along the segment you need (your point should be to the right of the directed segment). Go to the endpoint, turn right (that is, take the segment next to the one you came from, counting counterclockwise). Continue going to the next endpoint and turning right until you hit the closest segment again.
Then, check if the polygon encloses the given point. If it is not, then you have found an "island"; remove that polygon and the entire connected component it belongs to, and start over by selecting the nearest segment again. Connected components can be found with a simple DFS.
This gives you a clockwise-oriented polygon. If you want counterclockwise, which is often the default "positive" direction in both the software an the literature, start with the point to your left, and turn left at each intersection.
It surely helps if, given an endpoint, you can quickly find all segments incident with it.
This is really just an implementation of @n.m.'s answer. This is how far I got before the bounty expired; it's completely untested.
def smallestPolygon(point,segments):
"""
:param point: the point (x,y) to surrond
:param segments: the line segments ( ( (a,b),(c,d) ) , . . . )
that will make the polygon
(assume no duplicates and no intersections)
:returns: the segments forming the smallest polygon containing the point
"""
connected = list(segments)
def endPointMatches(segment1,segment2):
"""
Which endpoints of segment1 are in segment2 (with (F,F) if they're equal)
"""
if ( segment1 == segment2 or segment1 == segment2[::-1] ):
return ( False, False )
return ( segment1[0] in segment2 , segment1[1] in segment2 )
keepPruning = True
while ( keepPruning ):
keepPruning = False
for segment in connected:
from functors import partial
endPointMatcher = partial(endPointMatches,segment1=segment)
endPointMatchings = map(endPointMatcher,connected)
if ( not and(*map(any,zip(*endPointMatchings))) ):
connected.remove(segment)
keepPruning = True
def xOfIntersection(seg,y):
"""
:param seg: a line segment ( (x0,y0), (x1,y1) )
:param y: a y-coordinate
:returns: the x coordinate so that (x,y) is on the line through the segment
"""
return seg[0][0]+(y-seg[0][1])*(seg[1][0]-seg[0][0])/(seg[1][1]-seg[0][1])
def aboveAndBelow(segment):
return ( segment[0][1] <= point[1] ) != ( segment[1][1] <= point[1] )
# look for first segment to the right
closest = None
smallestDist = float("inf")
for segment in filter(aboveAndBelow,connected):
dist = xOfIntersection(segment,point[1])-point[0]
if ( dist >= 0 and dist < smallestDist ):
smallestDist = dist
closest = segment
# From the bottom of closest:
# Go to the other end, look at the starting point, turn right until
# we hit the next segment. Take that, and repeat until we're back at closest.
# If there are an even number of segments directly to the right of point[0],
# then point is not in the polygon we just found, and we need to delete that
# connected component and start over
# If there are an odd number, then point is in the polygon we found, and we
# return the polygon
Approach. I suggest to interpret the input as a PSLG, G, which consists of vertices and edges. Then your question reduces to finding the face of G that is hit by the point p. This is done by shooting a ray from p to some direction to find an edge of the boundary of the face and traverse the boundary of the face in some direction. However, the first edge hit could be a face that is not hit by p, but itself enclosed by the face hit by p. Hence, we may need to keep search along the ray emanated by p.
Implementational details. In the code below I shoot a ray to east direction and run around the face in clockwise direction, i.e., at each vertex I take the next counter-clockwise edge, until I end up at the first vertex again. The resulting face is returned as a sequence of vertices of G.
If you want to return a simple polygon then you have to clean-up the input graph G by pruning of trees in G such that only the simple faces remain.
def find_smallest_enclosing_polygon(G, p, simple=False):
"""Find face of PSLG G hit by point p. If simple is True
then the face forms a simple polygon, i.e., "trees"
attached to vertices are pruned."""
if simple:
# Make G comprise simple faces only, i.e., all vertices
# have degree >= 2.
done = False
while not done:
done = True
for v in [v in vertices if degree(v) <= 1]:
# Remove vertex v and all incident edges
G.remove(v)
done = False
# Shoot a ray from p to the east direction and get all edges.
ray = Ray(p, Vector(1, 0))
for e in G.iter_edges_hit(ray):
# There is no enclosing face; p is in the outer face of G
if e is None:
return None
# Let oriented edge (u, v) point clockwise around the face enclosing p
u, v = G.get_vertices(e)
if u.y < v.y
u, v = v, u
# Construct the enclosing face in clockwise direction
face = [u, v]
# While not closed
while face[-1] != face[0]:
# Get next counter-clockwise edge at last vertex at last edge.
# If face[-1] is of degree 1 then I assume to get e again.
e = G.get_next_ccw_edge(face[-2], face[-1])
face.append(G.get_opposite_vertex(e, face[-1]))
# Check whether face encloses p.
if contains(face, p):
return face
return None
Complexity. Let n denote the number of vertices. Note that in a PSLG the number of edges are in O(n). The pruning part may take O(n^2) time the way it is implemented above. But it could be one in O(n) time by identifying the degree-1 vertices and keep traversing from those.
The ray intersection routine can be implemented trivially in O(n) time. Constructing the face takes O(m) time, where m is the size of the polygon constructed. We may need to test multiple polygons but the sum of sizes of all polygons is still in O(n). The test contains(face, p) could be done by checking whether face contains an uneven number of edges returned by G.iter_edges_hit(ray), i.e., in O(m) time.
With some preprocessing you can build up a data structure that finds the face hit by p in O(log n) time by classical point location algorithms in computational geometry and you can store the resulting polygons upfront, for the simple and/or non-simple cases.
If you're doing this a number of times with the same lines and different points, it would be worthwhile preprocessing to figure out all the polygons. Then it's straightforward: draw a line from the point to infinity (conceptually speaking). Every time the you cross a line, increment the crossing count of each polygon the line is a part of. At the end, the first polygon with an odd crossing count is the smallest enclosing polygon. Since any arbitrary line will do just as well as any other (it doesn't even need to be straight), simplify the arithmetic by drawing a vertical or horizontal line, but watch out for crossing actual endpoints.
You could do this without preprocessing by creating the polygons as you cross each line. This basically reduces to n.m.'s algorithm but without all the special case checks.
Note that a line can belong to two polygons. Indeed, it could belong to more, but it's not clear to me how you would tell: consider the following:
+---------------------------+
| |
| +-------------------+ |
| | | |
| | +-----------+ | |
| | | | | |
| | | | | |
+---+---+-----------+---+---+
