Fast algorithm to generate rectangles that contains a number of 2D points

Question

I have one problem which I'm struggling with. Given the following:

• an array all_points containing 2D points, each points is represented as a tuple (x, y).
• an array musthave_points containing the indices of points that are in all_points.
• an integer m, with m < len(all_points).

Return a list of rectangles, in which a rectangle is represented by a tuple containing its 4 vertices ((x0, y0), (x1, y1), (x2, y2), (x3, y3)), each rectangle must satisfy the conditions below:

1. Contains m points from all_points, these m points must lay completely inside the rectangle, i.e not on either 4 of the rectangle's edges.
2. Contains all points from musthave_points. If musthave_points is an empty list, the rectangles only need to satisfy the first condition.

If there's no such rectangle, return an empty list. Two rectangles are considered "identical" if they contain the same subset of points and there should not be "identical" rectangles in the output.

Note: One simple brute-force solution is to first generate all combinations of m points, each of them contains all points from musthave_points. For each combination, create one rectangle that covers all points in the combination. Then count the number of points that lays inside the rectangle, if the number of points is m, it's a valid rectangle. But that solution runs in factorial time complexity. Can you come up with something faster than that?

I already implemented the brute-force as shown below, but it's terribly slow.

import itertools
import numpy as np 
import cv2 
import copy 
import sys 

from shapely.geometry import Point
from shapely.geometry.polygon import Polygon

# Credit: https://github.com/dbworth/minimum-area-bounding-rectangle/blob/master/python/min_bounding_rect.py
def minBoundingRect(hull_points_2d):
    #print "Input convex hull points: "
    #print hull_points_2d

    # Compute edges (x2-x1,y2-y1)
    edges = np.zeros((len(hull_points_2d) - 1, 2)) # empty 2 column array
    for i in range(len(edges)):
        edge_x = hull_points_2d[i+1, 0] - hull_points_2d[i, 0]
        edge_y = hull_points_2d[i+1, 1] - hull_points_2d[i, 1]
        edges[i] = [edge_x,edge_y]

    # Calculate edge angles   atan2(y/x)
    edge_angles = np.zeros((len(edges))) # empty 1 column array
    for i in range(len(edge_angles)):
        edge_angles[i] = np.math.atan2(edges[i,1], edges[i,0])

    # Check for angles in 1st quadrant
    for i in range(len(edge_angles)):
        edge_angles[i] = np.abs(edge_angles[i] % (np.math.pi/2)) # want strictly positive answers

    # Remove duplicate angles
    edge_angles = np.unique(edge_angles)

    # Test each angle to find bounding box with smallest area
    min_bbox = (0, sys.maxsize, 0, 0, 0, 0, 0, 0) # rot_angle, area, width, height, min_x, max_x, min_y, max_y
    for i in range(len(edge_angles) ):
        R = np.array([[np.math.cos(edge_angles[i]), np.math.cos(edge_angles[i]-(np.math.pi/2))], [np.math.cos(edge_angles[i]+(np.math.pi/2)), np.math.cos(edge_angles[i])]])

        # Apply this rotation to convex hull points
        rot_points = np.dot(R, np.transpose(hull_points_2d)) # 2x2 * 2xn

        # Find min/max x,y points
        min_x = np.nanmin(rot_points[0], axis=0)
        max_x = np.nanmax(rot_points[0], axis=0)
        min_y = np.nanmin(rot_points[1], axis=0)
        max_y = np.nanmax(rot_points[1], axis=0)

        # Calculate height/width/area of this bounding rectangle
        width = max_x - min_x
        height = max_y - min_y
        area = width*height

        # Store the smallest rect found first (a simple convex hull might have 2 answers with same area)
        if (area < min_bbox[1]):
            min_bbox = (edge_angles[i], area, width, height, min_x, max_x, min_y, max_y)

    # Re-create rotation matrix for smallest rect
    angle = min_bbox[0]   
    R = np.array([[np.math.cos(angle), np.math.cos(angle-(np.math.pi/2))], [np.math.cos(angle+(np.math.pi/2)), np.math.cos(angle)]])

    # Project convex hull points onto rotated frame
    proj_points = np.dot(R, np.transpose(hull_points_2d)) # 2x2 * 2xn
    #print "Project hull points are \n", proj_points

    # min/max x,y points are against baseline
    min_x = min_bbox[4]
    max_x = min_bbox[5]
    min_y = min_bbox[6]
    max_y = min_bbox[7]
    #print "Min x:", min_x, " Max x: ", max_x, "   Min y:", min_y, " Max y: ", max_y

    # Calculate center point and project onto rotated frame
    center_x = (min_x + max_x)/2
    center_y = (min_y + max_y)/2
    center_point = np.dot([center_x, center_y], R)
    #print "Bounding box center point: \n", center_point

    # Calculate corner points and project onto rotated frame
    corner_points = np.zeros((4,2)) # empty 2 column array
    corner_points[0] = np.dot([max_x, min_y], R)
    corner_points[1] = np.dot([min_x, min_y], R)
    corner_points[2] = np.dot([min_x, max_y], R)
    corner_points[3] = np.dot([max_x, max_y], R)

    return (angle, min_bbox[1], min_bbox[2], min_bbox[3], center_point, corner_points) # rot_angle, area, width, height, center_point, corner_points

class PatchGenerator:
    def __init__(self, all_points, musthave_points, m):
        self.all_points = copy.deepcopy(all_points)
        self.n = len(all_points)
        self.musthave_points = copy.deepcopy(musthave_points)
        self.m = m

    @staticmethod
    def create_rectangle(points):
        rot_angle, area, width, height, center_point, corner_points = minBoundingRect(points)
        return corner_points 

    @staticmethod
    def is_point_inside_rectangle(rect, point):
        pts = Point(*point)
        polygon = Polygon(rect)

        return polygon.contains(pts)

    def check_valid_rectangle(self, rect, the_complement):
        # checking if the rectangle contains any other point from `the_complement`
        for point in the_complement:
            if self.is_point_inside_rectangle(rect, point):
                return False
        return True 

    def generate(self):

        rects = [] 
        # generate all combinations of m points, including points from musthave_points
        the_rest_indices = list(set(range(self.n)).difference(self.musthave_points))
        comb_indices = itertools.combinations(the_rest_indices, self.m - len(self.musthave_points))
        comb_indices = [self.musthave_points + list(inds) for inds in comb_indices]

        # for each combination
        for comb in comb_indices:
            comb_points = np.array(self.all_points)[comb]
            ## create the rectangle that covers all m points
            rect = self.create_rectangle(comb_points)

            ## check if the rectangle is valid 
            the_complement_indices = list(set(range(self.n)).difference(comb))
            the_complement_points = list(np.array(self.all_points)[the_complement_indices])

            if self.check_valid_rectangle(rect, the_complement_points):
                rects.append([comb, rect]) # indices of m points and 4 vertices of the valid rectangle

        return rects 

if __name__ == '__main__':
    all_points = [[47.43, 20.5 ], [47.76, 43.8 ], [47.56, 23.74], [46.61, 23.73], [47.49, 18.94], [46.95, 25.29], [54.31, 23.5], [48.07, 17.77],
                        [48.2 , 34.87], [47.24, 22.07], [47.32, 27.05], [45.56, 17.95], [41.29, 19.33], [45.48, 28.49], [42.94, 15.24], [42.05, 34.3 ],
                        [41.04, 26.3 ], [45.37, 21.17], [45.44, 24.78], [44.54, 43.89], [30.49, 26.79], [40.55, 22.81]]
    musthave_points =  [3, 5, 9]
    m = 17 
    patch_generator = PatchGenerator(all_points, musthave_points, 17)
    patches = patch_generator.generate()

Answer 1

Every such rectangle can be shrunk to the minimum size such that it still contains the same points. Thus you only need to check such minimal rectangles. Let n be the total number of points. Then there are at most n possible coordinates for the left side, and likewise for the other sides. For each possible pair of left and right side coordinates, you can do a linear sweep for the top and bottom coordinates. Final time complexity would be O(n^3).