Form a rectangle from plane points

Question

I have a set of red lines from which I get a set of green intersection points (visible on the screen):

Then I want to find the four points that most likely describe the rectangle (if there are several options, then choose the largest area). I read similar questions about how to find points that EXACTLY form a rectangle:

find if 4 points on a plane form a rectangle?

https://softwareengineering.stackexchange.com/questions/176938/how-to-check-if-4-points-form-a-square

There is an option to iterate over all four points and calculate the probability that they form a rectangle (or some coefficient of similarity to a rectangle). Suppose at the moment we are considering four points A, B, C, D. I tried 2 similarity functions:

1. ,

where <> denotes dot product, and || - vector norm.

2. ,

where std is the standard deviation of the distances from the vertices to the center of mass of the assumed rectangle, and mean is the average distance.

Both functions did not perform well.

Is there a way to introduce a function that is close to 1 when the four points of the plane are close to the vertices of the rectangle and equal to 0 when they are at the position farthest from the rectangle (assuming they are on 1 line)?

Answer 1

I can't really speak to finding an appropriate cost function for scoring what a "good" rectangle is. From the comments it looks like there's a lot of discussion, but no consensus. So for now I'm going to just use a scoring function that penalizes four-point shapes for having angles that are further away from 90 degrees. Specifically, I'm summing the squared distance. If you want to have a different scoring metric you can replace the calculation in the scoreFunc function.

I set up an interactive window where you can click to add points. When you press 'q' it'll take those points, find all possible combinations (not permutations) of 4 points, and then run the scoring function on each and draws the best.

I'm using a recursive, brute-force search. To avoid having a ton of duplicates I came up with a hashing function that works regardless of order. I used prime numbers to ID each point and the hashing function just takes the product of the ID's of the points. This ensures that (1,3,5,7) is the same as (3,1,7,5). I used primes because the product of primes is unique in this situation (they can't be factorized and clumped because they're primes).

After the search I have to make sure that the points are ordered in such a way that the lines aren't intersecting. I'm taking advantage of OpenCV's contourArea to do that calculation for me. I can swap the first point with it's horizontal and vertical neighbor and compare the areas to the original. "Bowtie" shapes from intersecting lines will have less area (I'm pretty sure they actually get zero area because they don't count as closed shapes) than a non-intersection shape.

import cv2
import numpy as np
import math

# get mouse click
click_pos = None;
click = False;
def mouseClick(event, x, y, flags, param):
    # hook to globals
    global click_pos;
    global click;

    # check for left mouseclick 
    if event == cv2.EVENT_LBUTTONDOWN:
        click = True;
        click_pos = (x,y);

# prime hash function
def phash(points):
    total = 1;
    for point in points:
        total *= point[0];
    return total;

# checks if an id is already present in list
def isInList(point, curr_list):
    pid = point[0];
    for item in curr_list:
        if item[0] == pid:
            return True;
    return False;

# look for rectangles
def getAllRects(points, curr_list, rects, curr_point):
    # check if already in curr_list
    if isInList(curr_point, curr_list):
        return curr_list;

    # add self to list
    curr_list.append(curr_point);

    # check end condition
    if len(curr_list) == 4:
        # add to dictionary (no worry for duplicates)
        rects[phash(curr_list)] = curr_list[:];
        curr_list = curr_list[:-1];
        return curr_list;

    # continue search
    for point in points:
        curr_list = getAllRects(points, curr_list, rects, point);
    curr_list = curr_list[:-1];
    return curr_list;

# checks if a number is prime
def isPrime(num):
    bound = int(math.sqrt(num));
    curr = 3;
    while curr <= bound:
        if num % curr == 0:
            return False;
        # skip evens
        curr += 2;
    return True;

# generate prime number id's for each point
def genPrimes(num):
    primes = [];
    curr = 1;
    while len(primes) < num:
        if isPrime(curr):
            primes.append(curr);

        # +2 to skip evens
        curr += 2;
    return primes;

# swap sides (fix intersecting lines issue)
def swapH(box):
    new_box = np.copy(box);
    new_box[0] = box[1];
    new_box[1] = box[0];
    return new_box;
def swapV(box):
    new_box = np.copy(box);
    new_box[0] = box[3];
    new_box[3] = box[0];
    return new_box;

# removes intersections
def noNoodles(box):
    # get three variants
    hbox = swapH(box);
    vbox = swapV(box);

    # get areas and choose max
    sortable = [];
    sortable.append([cv2.contourArea(box), box]);
    sortable.append([cv2.contourArea(hbox), hbox]);
    sortable.append([cv2.contourArea(vbox), vbox]);
    sortable.sort(key = lambda a : a[0]);
    return sortable[-1][1];

# 2d distance
def dist2D(one, two):
    dx = one[0] - two[0];
    dy = one[1] - two[1];
    return math.sqrt(dx*dx + dy*dy);

# angle between three points (the last point is the middle)
# law of cosines
def angle3P(p1, p2, p3):
    # get distances
    a = dist2D(p3, p1);
    b = dist2D(p3, p2);
    c = dist2D(p1, p2);

    # calculate angle // assume a and b are nonzero
    numer = c**2 - a**2 - b**2;
    denom = -2 * a * b;
    if denom == 0:
        denom = 0.000001;
    rads = math.acos(numer / denom);
    degs = math.degrees(rads);
    return degs;


# calculates a score
def scoreFunc(box):
    # for each point, calculate angle
    angles = [];
    for a in range(len(box)):
        prev = box[a-2][0];
        curr = box[a-1][0];
        next = box[a][0];
        angles.append(angle3P(prev, next, curr));

    # for each angle, score on squared distance from 90
    score = 0;
    for angle in angles:
        score += (angle - 90)**2;
    return score;

# evaluates each box (assigns a score)
def evaluate(boxes):
    sortable = [];
    for box in boxes:
        # INSERT YOUR OWN SCORING FUNC HERE
        sortable.append([scoreFunc(box), box]);
    sortable.sort(key = lambda a : a[0]);
    return sortable;

# set up callback
cv2.namedWindow("Display");
cv2.setMouseCallback("Display", mouseClick);

# set up screen
res = (600,600,3);
bg = np.zeros(res, np.uint8);

# loop
done = False;
points = [];
while not done:
    # reset display
    display = np.copy(bg);

    # check for new click
    if click:
        click = False;
        points.append(click_pos);

    # draw points
    for point in points:
        cv2.circle(display, point, 4, (0,200,0), -1);

    # show
    cv2.imshow("Display", display);
    key = cv2.waitKey(1);

    # check keypresses
    done = key == ord('q');

# generate prime number id's for each point
# if you have a lot of points, it would be worth it
# to just have a .txt file with a bunch of pre-gen primes in it
primes = genPrimes(len(points));
print(primes);
withPrimes = [];
for a in range(len(points)):
    withPrimes.append([primes[a], points[a]]);

# run brute-force search over all points
rects = {};
for a in range(len(withPrimes)):
    getAllRects(withPrimes, [], rects, withPrimes[a]);
print(len(rects));

# extract just the points (don't need the prime id's anymore)
boxes = [];
for key in rects:
    box = [];
    for item in rects[key]:
        box.append([item[1]]);
    boxes.append(np.array(box));

# go through all of the boxes and un-intersect their sides
for a in range(len(boxes)):
    boxes[a] = noNoodles(boxes[a]);

# draw each one to check for noodles
# for box in boxes:
#   blank = np.zeros_like(bg, np.uint8);
#   cv2.drawContours(blank, [box], -1, (255,255,255), -1);
#   cv2.imshow("Box", blank);
#   cv2.waitKey(0);

# noodles have been squared get best box
sortedBoxes = evaluate(boxes);
bestBox = sortedBoxes[0][1];

# draw
blank = np.zeros_like(bg, np.uint8);
cv2.drawContours(blank, [bestBox], -1, (255,255,255), -1);
for point in points:
    cv2.circle(blank, point, 4, (0,200,0), -1);
cv2.imshow("Best", blank);
cv2.waitKey(0);