Decision

In the space of all lines in the graph, the lines along the corner correspond exactly to where the number or intersections are about to decrease. In other words, each of them forms a local maximum.

And for every row that goes through at least one corner, there is a corresponding row going through two corners with the same number of intersections.

The conclusion is that we only need to check the lines formed by the two rectangular corners, since they form a set that fully represents the local maxima of our problem. Of these, we choose the one that has the greatest intersection.

Complexity of time

This solution must first reconstruct all the lines that go through the two corners. The number of such lines is O (n ^ 2).

Then we need to count the number of intersections between this line and the rectangle. This can obviously be done in O (n) by comparing to each rectangle.

Maybe a more efficient way to proceed, but we know that this algorithm is then at most O (n ^ 3).

Python3 implementation

Here is a Python implementation of this algorithm. I was focusing on readability rather than efficiency, but it does exactly what is stated above.

def get_best_line(rectangles):
    """
    Given a set of rectangles, return a line which intersects the most rectangles.
    """

    # Recover all corners from all rectangles
    corners = set()
    for rectangle in rectangles:
        corners |= set(rectangle.corners)

    corners = list(corners)

    # Recover all lines passing by two corners
    lines = get_all_lines(corners)

    # Return the one which has the highest number of intersections with rectangles
    return max(
        ((line, count_intersections(rectangles, line)) for line in lines),
        key=lambda x: x[1])

      

        
        
        
      

    

This implementation uses the following helpers.

def get_all_lines(points):
    """
    Return a generator providing all lines generated
    by a combination of two points out of 'points'
    """
    for i in range(len(points)):
        for j in range(i, len(points)):
            yield Line(points[i], points[j])

def count_intersections(rectangles, line):
    """
    Return the number of intersections with rectangles
    """
    count = 0

    for rectangle in rectangles:
        if line in rectangle:
           count += 1

    return count

      

        
        
        
      

    

And here is a class definition that serves as a data structure for rectangles and strings.

import itertools
from decimal import Decimal

class Rectangle:
    def __init__(self, x_range, y_range):
        """
        a rectangle is defined as a range in x and a range in y.
        By example, the rectangle (0, 0), (0, 1), (1, 0), (1, 1) is given by
        Rectangle((0, 1), (0, 1))
        """
        self.x_range = sorted(x_range)
        self.y_range = sorted(y_range)

    def __contains__(self, line):
        """
        Return whether 'line' intersects the rectangle.
        To do so we check if the line intersects one of the diagonals of the rectangle
        """
        c1, c2, c3, c4 = self.corners

        x1 = line.intersect(Line(c1, c4))
        x2 = line.intersect(Line(c2, c3))

        if x1 is True or x2 is True \
                or x1 is not None and self.x_range[0] <= x1 <= self.x_range[1] \
                or x2 is not None and self.x_range[0] <= x2 <= self.x_range[1]:

            return True

        else:
            return False

    @property
    def corners(self):
        """Return the corners of the rectangle sorted in dictionary order"""
        return sorted(itertools.product(self.x_range, self.y_range))


class Line:
    def __init__(self, point1, point2):
        """A line is defined by two points in the graph"""
        x1, y1 = Decimal(point1[0]), Decimal(point1[1])
        x2, y2 = Decimal(point2[0]), Decimal(point2[1])
        self.point1 = (x1, y1)
        self.point2 = (x2, y2)

    def __str__(self):
        """Allows to print the equation of the line"""
        if self.slope == float('inf'):
            return "y = {}".format(self.point1[0])

        else:
            return "y = {} * x + {}".format(round(self.slope, 2), round(self.origin, 2))

    @property
    def slope(self):
        """Return the slope of the line, returning inf if it is a vertical line"""
        x1, y1, x2, y2 = *self.point1, *self.point2

        return (y2 - y1) / (x2 - x1) if x1 != x2 else float('inf')

    @property
    def origin(self):
        """Return the origin of the line, returning None if it is a vertical line"""
        x, y = self.point1

        return y - x * self.slope if self.slope != float('inf') else None

    def intersect(self, other):
        """
        Checks if two lines intersect.
        Case where they intersect: return the x coordinate of the intersection
        Case where they do not intersect: return None
        Case where they are superposed: return True
        """

        if self.slope == other.slope:

            if self.origin != other.origin:
                return None

            else:
                return True

        elif self.slope == float('inf'):
            return self.point1[0]

        elif other.slope == float('inf'):
            return other.point1[0]

        elif self.slope == 0:
            return other.slope * self.origin + other.origin

        elif other.slope == 0:
            return self.slope * other.origin + self.origin

        else:
            return (other.origin - self.origin) / (self.slope - other.slope)

      

        
        
        
      

    

Example

Here is a working example of the above code.

rectangles = [
    Rectangle([0.5, 1], [0, 1]),
    Rectangle([0, 1], [1, 2]),
    Rectangle([0, 1], [2, 3]),
    Rectangle([2, 4], [2, 3]),
]

# Which represents the following rectangles (not quite to scale)
#
#  *
#  *   
#
# **     **
# **     **
#
# **
# **

      

        
        
        
      

    

We can clearly see that the optimal solution is to find a line that goes over three rectangles, and that is indeed what it outputs.

print('{} with {} intersections'.format(*get_best_line(rectangles)))
# prints: y = 0.50 * x + -5.00 with 3 intersections