Finding all quadrilaterals in a set of intersections

Question

I want to take all the intersections of a set of lines and find all the convex quadrilaterals they create. I am not sure if there is an algorithm that works perfect for this, or if I need to loop through and create my own.

I have an array of lines, and all their intersections.

Lines and intersections:

Example Quadrilaterals 1:

Example Quadrilaterals 2

In this case, I would come out with 8 quadrilaterals.

How can I achieve this? If there isn't an algorithm I can implement for this, how can I check each intersection with other intersections to determine if they make a convex quadrilateral?

Answer 1

There is a simple, non-speedy, brute-force over-all algorithm to find those quadrilaterals. However, first you would need to clarify some definitions, especially that of a "quadrilateral." Do you count it as a quadrilateral if it has zero area, such as when all the vertices are collinear? Do you count it as a quadrilateral if it self-intersects or crosses? Do you count it if it is not convex? Do you count it if two adjacent sides are straight (which includes consecutive vertices identical)? What about if the polygon "doubles back" on itself so the result looks like a triangle with one side extended?

Here is a top-level algorithm: Consider all combinations of the line segments taken four at a time. If there are n line segments then there are n*(n-1)*(n-2)*(n-3)/24 combinations. For each combination, look at the intersections of pairs of these segments: there will be at most 6 intersections. Now see if you can make a quadrilateral from those intersections and segments.

This is brute-force, but at least it is polynomial in execution time, O(n^4). For your example of 8 line segments that means considering 70 combinations of segments--not too bad. This could be sped up somewhat by pre-calculating the intersection points: there are at most n*(n-1)/2 of them, 28 in your example.

Does this overall algorithm meet your needs? Is your last question "how can I check each intersection with other intersections to determine if they make a quadrilateral?" asking how to implement my statement "see if you can make a quadrilateral from those intersections and segments"? Do you still need an answer to that? You would need to clarify your definition of a quadrilateral before I could answer that question.

---

I'll explain the definitions of "quadrilateral" more. This diagram shows four line segments "in general position," where each segment intersects all the others and no three segments intersect in the same point.

Here are (some of) the "quadrilaterals" arising from those four line segments and six intersection points.

• 1 simple and convex (ABDE)
• 1 simple and not convex (ACDF)
• 1 crossing (BCEF)
• 4 triangles with an extra vertex on a triangle's side (ABCE, ACDE, ABDF, ABFE). Note that the first two define the same region with different vertices, and the same is true of the last two.
• 4 "double-backs" which looks like a triangle with one side extended (ACEF, BCDF, BDEF, CDEF)

Depending on how you define "quadrilateral" and "equal" you could get anywhere from 1 to 11 of them in that diagram. Wikipedia's definition would include only the first, second, and fourth in my list--I am not sure how that counts the "duplicates" in my fourth group. And I am not even sure that I found all the possibilities in my diagram, so there could be even more.

---

I see we are now defining a quadrilateral as outlined by four distinct line segments that are sub-segments of the given line segments that form a polygon that is strictly convex--the vertex angles are all less than a straight angle. This still leaves an ambiguity in a few edge cases--what if two line segments overlap more than at one point--but let's leave that aside other than defining that two such line segments have no intersection point. Then this algorithm, pseudo-code based on Python, should work.

We need a function intersection_point(seg1, seg2) that returns the intersection point of the two given line segments or None if there is none or the segments overlap. We also need a function polygon_is_strictly_convex(tuple of points) that returns True or False depending on if the tuple of points defines a strictly-convex polygon, with the addition that if any of the points is None then False is returned. Both those functions are standard in computational geometry. Note that "combination" in the following means that for each returned combination the items are in sorted order, so of (seg1, seg2) and (seg2, seg1) we will get exactly one of them. Python's itertools.combinations() does this nicely.

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, which for your segments gives your eight polygons. First, here are the utility, geometry routines, collected into 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's 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 little complexity to use only a "canonical form" of a pair of line segments and for a polygon, to reduce the memory usage of the intersections and polygons 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 from that is:

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))

Answer 2

A O(intersection_count²) algorithm is as follows:

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