https://stackoverflow.com/questions/68902521/how-many-points-of-intersections-are-possible-rectangle

import itertools
import math

import numpy as np
import cv2

# ============================================================================

intersections = [
    (406, 254), (123, 304), (120, 304), (405, 254)
    , (410, 253), (124, 303), (222, 286), (262, 56)
    , (381, 220), (278, 78), (260, 54), (408, 257)
    , (21, 86), (57, 162), (22, 87), (81, 215)
    , (122, 302), (22, 85), (266, 55), (261, 56)
    , (125, 73), (18, 86), (23, 87), (68, 189)
    , (119, 303), (331, 148), (264, 53), (407, 257)
    , (260, 54), (412, 256), (19, 88), (124, 302)
]

# NB: Overriding the original point array with a set
# of simpler lines for testing and demonstration...
intersections = [(10,10), (160,10), (80,10), (270,10), (81,110), (160,111)]

# ============================================================================

def draw_intersections(intersections):
    output = np.zeros((320,480,3), np.uint8)

    for intersection in intersections:
        cv2.drawMarker(output, intersection, (191,255,191))

    return output

# ============================================================================

def distance(a, b):
    dx, dy = (b[0] - a[0]), (b[1] - a[1])
    return math.sqrt(dx * dx + dy * dy)

def is_triangle_clockwise(pA, pB, pC):
    # https://stackoverflow.com/a/1165943/3962537
    total = ((pB[0] - pA[0]) * (pB[1] + pA[1])
        + (pC[0] - pB[0]) * (pC[1] + pB[1])
        + (pA[0] - pC[0]) * (pA[1] + pC[1]))
    return total >= 0

def get_inner_angle(pC, pA, pB):
    vAB = [(pA[0] - pB[0]), (pA[1] - pB[1])]
    vAC = [(pA[0] - pC[0]), (pA[1] - pC[1])]
    dot_product = (vAB[0] * vAC[0]) + (vAB[1] * vAC[1])
    dAB = distance(pA, pB)
    dAC = distance(pA, pC)
    cos_angle = np.clip(dot_product / dAB / dAC, -1.0, 1.0)
    angle = math.acos(cos_angle)
    angle_deg = math.degrees(angle)
    return angle_deg

# ============================================================================

def detect_rectangle(points):
    MINIMUM_DISTANCE = 20
    MAX_ANGLE_SUM_ERROR = 1.0
    MAX_CORNER_ANGLE_ERROR = 10.0

    pA = points[0] # Select anchror point

    distances = [
        distance(pA, points[1])
        , distance(pA, points[2])
        , distance(pA, points[3])
    ]
    order = np.argsort(distances)
    if distances[order[0]] < MINIMUM_DISTANCE:
        return None

    pC = points[order[2]+1] # Diagonal, most distant point
    pB, pD = points[order[0]+1], points[order[1]+1] # The other two corners

    # Make sure the points are in clockwise order (assuming it's a rectangle)...
    if not is_triangle_clockwise(pA, pB, pC):
        pB, pD = pD, pB

    angleA = get_inner_angle(pD, pA, pB)
    angleB = get_inner_angle(pA, pB, pC)
    angleC = get_inner_angle(pB, pC, pD)
    angleD = get_inner_angle(pC, pD, pA)
    angle_sum = angleA + angleB + angleC + angleD
    if math.fabs(360.0 - angle_sum) > MAX_ANGLE_SUM_ERROR:
        return None

    corner_angle_errors = np.abs(90.0 - np.array([angleA, angleB, angleC, angleD]))
    if np.max(corner_angle_errors) > MAX_CORNER_ANGLE_ERROR:
        return None

    return ([pA, pB, pC, pD], [angleA, angleB, angleC, angleD], angle_sum)

# ============================================================================

# Eliminate duplicate points
intersections = list(set(intersections))

# Process all combinations of 4 intersection points
for points in itertools.combinations(intersections, 4):
    result = detect_rectangle(points)
    if result:
        print(result)
        polygon = np.array(result[0], np.int32).reshape((-1,1,2))
        output = draw_intersections(intersections)
        cv2.polylines(output, [polygon], True, (0,255,255))
        cv2.imshow("", output)
        cv2.waitKey()