Find all quadrangles in intersection set
I want to take all the intersections of a set of rows and find all the convex quads they create. I'm not sure if there is an algorithm that is perfect for this, or if I need to loop through and create my own.
I have an array of strings and all their intersections.
Lines and intersections:
Quadrilaterals example 1:
Quadrilaterals example 2
In this case, I would go out with 8 quadrilaterals.
How can I achieve this? If there is no algorithm for this, how can I check each intersection with other intersections to determine if they form a convex quadrilateral?
source to share
There is a simple, not fast, complete brute force algorithm to find these quadrangles. However, first you need to clarify some definitions, especially the word "quadrangle". Do you think of it as four-sided if it has zero area, for example when all vertices are collinear? Do you consider it four-sided if it self-intersects or intersects? Do you think this if it is not convex? Do you think it is if two adjacent sides are straight (which includes consecutive vertices that are the same)? How about the polygon "doubles back" by itself, so the result looks like a triangle on one side?
Here's the top-level algorithm: Consider all combinations of line segments taken four at a time. If there are segments n
, then there are combinations n*(n-1)*(n-2)*(n-3)/24
. For each combination, look at the intersections of pairs of these segments: there will be no more than 6 intersections. Now let's see if you can make a quad from these intersections and segments.
It's brute force, but at least it's a polynomial at runtime O (n ^ 4). For your 8 line segment example, which means considering 70 segment combinations is not too bad. This could be accelerated by pre-calculating the intersection points: in your case there are no more than n*(n-1)/2
28 in your example.
Does this general algorithm suit your needs? Your last question is, "How can I test each intersection with other intersections to determine if they make a quadrilateral?" asking how to implement my expression "see if you can make a quad from these intersections and segments"? Still need an answer? You will need to clarify your definition of a quadrilateral before I can answer this question.
I'll explain the definitions of "quadrangle" more. This diagram shows four line segments in general position, where each segment intersects all others, and three points do not intersect at the same point.
Here are (some of) "quadrangles" arising from these four line segments and six intersection points.
- 1 simple and convex (ABDE)
- 1 simple and non-convex (ACDF)
- 1 intersection (BCEF)
- 4 triangles with an additional vertex on the side of the triangle (ABCE, ACDE, ABDF, ABFE). Note that the first two define the same area with different vertices, and the same applies to the last two.
- 4 "double-backs" that look like a triangle on one side (ACEF, BCDF, BDEF, CDEF)
Depending on how you define "quadrilateral" and "equal", you can get anywhere from 1 to 11 of them in this diagram. Wikipedia's definition will only include the first, second and fourth on my list - I'm not sure how this accounts for "duplicates" in my fourth group. And I'm not even sure I found all the possibilities in the diagram, so there might be even more.
I see that we now define a quadrilateral, denoted by four different line segments, which are sub-segments of the line segments that form a polygon that is strictly convex - the vertex angles are less and less than the right angle. This still leaves ambiguity in several edge cases - what if two line segments overlap at more than one point - but let's leave that aside, other than specifying that two such line segments do not have an intersection point. Then this algorithm, Python based pseudocode, should work.
We need a function intersection_point(seg1, seg2)
that returns the point of intersection of two given line segments, or None
if they are absent or the segments overlap. We also need a function polygon_is_strictly_convex(tuple of points)
that returns True
or False
depending on whether the tuple of points is a strictly convex polygon, with the addition that if any of the points None
is returned False
, both of these functions are standard in computational geometry. Note that "combination" in the following means that for each combination returned, the elements are in sorted order, so from (seg1, seg2)
and (seg2, seg1)
we will get exactly one of them. Python itertools.combinations()
does this beautifully.
intersections = {} # empty dictionary/hash table
for each combination (seg1, seg2) of given line segments:
intersections[(seg1, seg2)] = intersection_point(seg1, seg2)
quadrilaterals = emptyset
for each combination (seg1, seg2, seg3, seg4) of given line segments:
for each tuple (sega, segb, segc, segc) in [
(seg1, seg2, seg3, seg4),
(seg1, seg2, seg4, seg3),
(seg1, seg3, seg2, seg4)]:
a_quadrilateral = (intersections[(sega, segb)],
intersections[(segb, segc)],
intersections[(segc, segd)],
intersections[(segd, sega)])
if polygon_is_strictly_convex(a_quadrilateral):
quadrilaterals.add(a_quadrilateral)
break # only one possible strictly convex quad per 4 segments
Here is my actual, tested, Python 3.6 code that gives your eight polygons for your segments. First, here's a utility, geometry, assembled into a module rdgeometry
.
def segments_intersection_point(segment1, segment2):
"""Return the intersection of two line segments. If none, return
None.
NOTES: 1. This version returns None if the segments are parallel,
even if they overlap or intersect only at endpoints.
2. This is optimized for assuming most segments intersect.
"""
try:
pt1seg1, pt2seg1 = segment1 # points defining the segment
pt1seg2, pt2seg2 = segment2
seg1_delta_x = pt2seg1[0] - pt1seg1[0]
seg1_delta_y = pt2seg1[1] - pt1seg1[1]
seg2_delta_x = pt2seg2[0] - pt1seg2[0]
seg2_delta_y = pt2seg2[1] - pt1seg2[1]
denom = seg2_delta_x * seg1_delta_y - seg1_delta_x * seg2_delta_y
if denom == 0.0: # lines containing segments are parallel or equal
return None
# solve for scalars t_seg1 and t_seg2 in the vector equation
# pt1seg1 + t_seg1 * (pt2seg1 - pt1seg1)
# = pt1seg2 + t_seg2(pt2seg2 - pt1seg2) and note the segments
# intersect iff 0 <= t_seg1 <= 1, 0 <= t_seg2 <= 1 .
pt1seg1pt1seg2_delta_x = pt1seg2[0] - pt1seg1[0]
pt1seg1pt1seg2_delta_y = pt1seg2[1] - pt1seg1[1]
t_seg1 = (seg2_delta_x * pt1seg1pt1seg2_delta_y
- pt1seg1pt1seg2_delta_x * seg2_delta_y) / denom
t_seg2 = (seg1_delta_x * pt1seg1pt1seg2_delta_y
- pt1seg1pt1seg2_delta_x * seg1_delta_y) / denom
if 0 <= t_seg1 <= 1 and 0 <= t_seg2 <= 1:
return (pt1seg1[0] + t_seg1 * seg1_delta_x,
pt1seg1[1] + t_seg1 * seg1_delta_y)
else:
return None
except ArithmeticError:
return None
def orientation3points(pt1, pt2, pt3):
"""Return the orientation of three 2D points in order.
Moving from Pt1 to Pt2 to Pt3 in cartesian coordinates:
1 means counterclockwise (as in standard trigonometry),
0 means straight, back, or stationary (collinear points),
-1 means counterclockwise,
"""
signed = ((pt2[0] - pt1[0]) * (pt3[1] - pt1[1])
- (pt2[1] - pt1[1]) * (pt3[0] - pt1[0]))
return 1 if signed > 0.0 else (-1 if signed < 0.0 else 0)
def is_convex_quadrilateral(pt1, pt2, pt3, pt4):
"""Return True if the quadrilateral defined by the four 2D points is
'strictly convex', not a triangle nor concave nor self-intersecting.
This version allows a 'point' to be None: if so, False is returned.
NOTES: 1. Algorithm: check that no points are None and that all
angles are clockwise or all counter-clockwise.
2. This does not generalize to polygons with > 4 sides
since it misses star polygons.
"""
if pt1 and pt2 and pt3 and pt4:
orientation = orientation3points(pt4, pt1, pt2)
if (orientation != 0 and orientation
== orientation3points(pt1, pt2, pt3)
== orientation3points(pt2, pt3, pt4)
== orientation3points(pt3, pt4, pt1)):
return True
return False
def polygon_in_canonical_order(point_seq):
"""Return a polygon, reordered so that two different
representations of the same geometric polygon get the same result.
The result is a tuple of the polygon points. `point_seq` must be
a sequence of 'points' (which can be anything).
NOTES: 1. This is intended for the points to be distinct. If two
points are equal and minimal or adjacent to the minimal
point, which result is returned is undefined.
"""
pts = tuple(point_seq)
length = len(pts)
ndx = min(range(length), key=pts.__getitem__) # index of minimum
if pts[(ndx + 1) % length] < pts[(ndx - 1) % length]:
return (pts[ndx],) + pts[ndx+1:] + pts[:ndx] # forward
else:
return (pts[ndx],) + pts[:ndx][::-1] + pts[ndx+1:][::-1] # back
def sorted_pair(val1, val2):
"""Return a 2-tuple in sorted order from two given values."""
if val1 <= val2:
return (val1, val2)
else:
return (val2, val1)
And here is the code for my algorithm. I added a bit of complexity to use only the "canonical" shape of a pair of line segments and for a polygon to reduce memory usage for intersection and polygon containers.
from itertools import combinations
from rdgeometry import segments_intersection_point, \
is_strictly_convex_quadrilateral, \
polygon_in_canonical_order, \
sorted_pair
segments = [(( 2, 16), (22, 10)),
(( 4, 4), (14, 14)),
(( 4, 6), (12.54, 0.44)),
(( 4, 14), (20, 6)),
(( 4, 18), (14, 2)),
(( 8, 2), (22, 16))]
intersections = dict()
for seg1, seg2 in combinations(segments, 2):
intersections[sorted_pair(seg1, seg2)] = (
segments_intersection_point(seg1, seg2))
quadrilaterals = set()
for seg1, seg2, seg3, seg4 in combinations(segments, 4):
for sega, segb, segc, segd in [(seg1, seg2, seg3, seg4),
(seg1, seg2, seg4, seg3),
(seg1, seg3, seg2, seg4)]:
a_quadrilateral = (intersections[sorted_pair(sega, segb)],
intersections[sorted_pair(segb, segc)],
intersections[sorted_pair(segc, segd)],
intersections[sorted_pair(segd, sega)])
if is_strictly_convex_quadrilateral(*a_quadrilateral):
quadrilaterals.add(polygon_in_canonical_order(a_quadrilateral))
break # only one possible strictly convex quadr per 4 segments
print('\nThere are {} strictly convex quadrilaterals, namely:'
.format(len(quadrilaterals)))
for p in sorted(quadrilaterals):
print(p)
And the printout:
There are 8 strictly convex quadrilaterals, namely:
((5.211347517730497, 5.211347517730497), (8.845390070921987, 2.8453900709219857), (11.692307692307693, 5.692307692307692), (9.384615384615383, 9.384615384615383))
((5.211347517730497, 5.211347517730497), (8.845390070921987, 2.8453900709219857), (14.666666666666666, 8.666666666666668), (10.666666666666666, 10.666666666666666))
((5.211347517730497, 5.211347517730497), (8.845390070921987, 2.8453900709219857), (17.384615384615387, 11.384615384615383), (12.769230769230768, 12.76923076923077))
((6.0, 14.8), (7.636363636363637, 12.181818181818182), (10.666666666666666, 10.666666666666666), (12.769230769230768, 12.76923076923077))
((6.0, 14.8), (7.636363636363637, 12.181818181818182), (14.666666666666666, 8.666666666666668), (17.384615384615387, 11.384615384615383))
((9.384615384615383, 9.384615384615383), (10.666666666666666, 10.666666666666666), (14.666666666666666, 8.666666666666668), (11.692307692307693, 5.692307692307692))
((9.384615384615383, 9.384615384615383), (11.692307692307693, 5.692307692307692), (17.384615384615387, 11.384615384615383), (12.769230769230768, 12.76923076923077))
((10.666666666666666, 10.666666666666666), (12.769230769230768, 12.76923076923077), (17.384615384615387, 11.384615384615383), (14.666666666666666, 8.666666666666668))
source to share
AO (intersection_count 2 ) algorithm looks like this:
For each intersection:
Add the the intersection point to
a hash table with the lines as the key.
Let int be a lookup function that returns
true iff the inputted lines intersect.
RectCount = 0
For each distinct pair of intersections a,b:
Let A be the list of lines that pass
through point a but not through b.
Let B '' '' '' through b but not a.
For each pair of lines c,d in A:
For each pair of lines e,f in B:
If (int(c,e) and int(d,f) and
!int(c,f) and !int(d,e)) or
(int(c,f) and int(d,e) and
!int(c,e) and !int(d,f)):
RectCount += 1
source to share