求解矩形相交问题¶
发布于：2020-09-01 | 分类：mathematics
本文求解矩形相交问题（Rectangle Intersection Problem）：已知n个水平或者竖直放置的矩形，对于每个矩形，求所有与之相交的矩形。一开始采用了暴力求解的方式，O(n^2)的时间复杂度在n<2000时尚能接受（~0.6s），继续增大n则导致计算时间爆炸了。实际上，矩形求交是一个研究过的问题，本文参考一篇论文中的算法进行了实现。
从判断两个矩形是否相交开始¶
首先，使用左下角(x0,y0)和右上角(x1,y1)两个点来描述一个标准矩形：
^ y
|
| (x1, y1)
4 | +-----+ +-------------+
| |A | |C |
3 | | | +-------------+ |
| | | |B | | |
2 | +-----+ | +-------------+
|(x0, y0) | |
1 | +-------------+
| x
+-------------------------------------->
0 1 2 3 4 5有多种思路可以判断相交，例如，两个矩形中心距离投影与边长的关系。这里从下图所示四种不相交的情况入手，其反面即为相交：
- B完全在A的上方：
B.y0 > A.y1 - B完全在A的下方：
B.y1 < A.y0 - B完全在A的左边：
B.x1 < A.x0 - B完全在A的右边：
B.x0 > A.x1
+----+ +--------------------+
| | | B |
| | +--------------------+
| | +----+
| | +----------+ | |
| B | | | | |
| | | A | | |
| | | | | B |
| | | | | |
| | +----------+ | |
+----+ | |
+---------------------+ | |
| B | | |
+---------------------+ +----+综合以上，得到矩形类及其判断相交的方法如下：
class Rect:
def __init__(self, points):
assert len(points)==4, 'invalid input'
self.x0, self.y0, self.x1, self.y1 = points
def intersects(self, rect):
''' if intersects with `rect`.'''
non_intersection = rect.x1<self.x0 or \
rect.x0>self.x1 or \
rect.y1<self.y0 or \
rect.y0>self.y1
return not non_intersection暴力求解¶
既然已经有了判断两个矩形是否相交的方法，自然想到两两判断所有矩形的暴力求解方法。最后用一个列表表示各个矩形的相交矩形，graph = [set() for i in range(num)]，也就是邻接表表示的图结构。
class Rects:
'''Collection of Rect instances.'''
def __init__(self):
self._rects = [] # type: list[Rect]
def append(self, instance): self._rects.append(instance)
def build_graph(self):
'''calculate the connectivity between every two rects by brute force method.'''
num = len(self._rects)
graph = [set() for i in range(num)]
for i in range(num):
for j in range(i+1, num):
if self._rects[i].intersects(self._rects[j]):
graph[i].add(j)
graph[j].add(i)
return graph
if __name__ == '__main__':
inputs = [(0.5,2,1.5,4), (2,1,4,3), (3,2,5,4)]
# initialize Rects
R = Rects()
for points in inputs: R.append(Rect(points))
# solving rectangle intersection
graph = R.build_graph()
print(graph)
# output:
# [set(), {2}, {1}]以上计算示例即为上图中的矩形A、B和C，结果与预期相符：A没有相交矩形，B与C相交，C与B相交。
优化的求解算法¶
暴力方法简单易行，但随着输入数量增大，性能急剧下降。于是寻求一种优化的算法，基本方向有两个：
- 以O(n*logn)代价沿着一个维度对矩形排序，这样可以避免每两个矩形都要进行一次判断
- 在一个维度的投影有重叠的情况下，只需按区间判断另一个维度的相交情况
本文对如下论文提出的算法进行了实现：
A Rectangle-Intersection Algorithm with Limited Resource Requirements
算法大意：
- 取所有
n个矩形的竖直边，并按非降的顺序排序得到共计m=2n条边的集合V - 然后重复如下描述的子过程
detect(V, m)：- 如果
m < 2则终止递归，否则， - 定义如下集合：
- 令
V1为V前一半即⌊m/2⌋个元素的集合， - 令
V2为V后一半即剩下元素的集合 - 令
S11为至少一条竖直边在V1中，且不经过V2区域的矩形的集合 - 令
S22为至少一条竖直边在V2中，且不经过V1区域的矩形的集合 - 令
S12为左竖直边在V1中且右竖直边跨过（spanning）V2区域的矩形的集合 - 令
S21为右竖直边在V2中且左竖直边跨过（spanning）V1区域的矩形的集合
- 令
- 分别对集合(
S12,S22)、(S21,S11)、(S12,S21)进行y区间是否重叠的判断stab(S12, S22)、stab(S21, S11)和stab(S12, S21) - 重复执行子过程：
detect(V1, ⌊m/2⌋)和detect(V2, m−⌊m/2⌋)
- 如果
以第一幅图矩形为例，第一次执行detect时的中间变量，此时无需执行stab操作。注意矩形B，虽然左竖直边在V1中，但右竖直边x=4不满足**跨过**V2，即要求x>5，所以B并不属于S12，同理分析亦不属于S21。
^ y
| + V1 = [0.5, 1.5, 2]
| (x1, y1) | V2 = [3, 4, 5]
4 | +-----+ | +-------------+ S11 = [A]
| |A | | |C | S22 = [C]
3 | | | +-------------+ | S12 = []
| | | |B | | | | S21 = []
2 | +-----+ | | +-------------+
|(x0, y0) | | |
1 | +-------------+
| | x
+-----------------+-------------------->
0 1 2 3 4 5
分别对上图的V1和V2执行detect操作：
^ y
| + V1 = [0.5]
| | V2 = [1.5, 2]
4 | +-----+ S11 = []
| |A | | S22 = [B]
3 | | | | +--- S12 = []
| | | | |B S21 = []
2 | +-----+ |
| | |
1 | | +---
| | x
+-------+------------------------------>
0 1 2 3 4 5 ^ y
| V1 = [3] +
| V2 = [4,5] |
4 | S11 = [] +-------------+
| S22 = [] |C | |
3 | S12 = [] ---------+ |
| S21 = [B] | | | |
2 | +-------------+
| B | |
1 | ---------+
| | x
+------------------------+------------->
0 1 2 3 4 5上面第二幅子图S21=[B]给出了集合S21的一个例子：右竖直边在V2，且左竖直边x=2跨过了V1即满足x<3。
注意到以上两图的len(V1)=1满足了detect的终止条件，接下来继续对V2执行detect：
^ y
| + V1 = [1.5]
| | V2 = [2]
4 | -+ | S11 = [A]
| A | | S22 = [B]
3 | | | +--+ S12 = []
| | | |B S21 = []
2 | -+ | |
| | |
1 | | +--+
| | x
+-----------+-------------------------->
0 1 2 3 4 5 ^ y
| V1 = [4] +
| V2 = [5] |
4 | S11 = [B] ---------+
| S22 = [] | C|
3 | S12 = [] --+ | |
| S21 = [C] | | |
2 | ---------+
| B | |
1 | --+ |
| | x
+------------------------------+------->
0 1 2 3 4 5终于，上图中可以执行stab(S21,S11)操作即判断B、C在y区间上的相交情况，显然它们是相交的；于是，最终得到一组相交矩形(B, C)。这里略去了stab的算法描述，即一维区间的求交问题。
综合来看，这个算法将矩形相交的判断分解为两个维度的相交，x方向通过排序和分治提高了效率，剩下的y方向就是一维区间的判断。根据S集合的定义，唯有(S12, S22)、(S21, S11)、(S12, S21)这三种组合在x方向出现重叠，所以算法中仅对它们进行y区间的stab判断。
于是，以上算法的一个Python实现：
def detect(V:list, num:int, index_groups:list):
''' divide and conquer in x-direction.
---
Args:
- V: rectangle-related x-edges data, [(index, (x0,y0,x1,y1), x), (...), ...]
- num: count of V instances, equal to len(V)
- index_groups: target adjacent list for connectivity between rects
'''
if num==1: return
# start/end points of left/right intervals
center_pos = int(num/2.0)
X0, X, X1 = V[0][-1], V[center_pos-1][-1], V[-1][-1]
# split into two groups
left = V[0:center_pos]
right = V[center_pos:]
# filter rects according to their position to each intervals
S11 = list(filter( lambda item: item[1][2]<=X, left ))
S12 = list(filter( lambda item: item[1][2]>=X1, left ))
S22 = list(filter( lambda item: item[1][0]>X, right ))
S21 = list(filter( lambda item: item[1][0]<=X0, right ))
# intersection in x-direction is fulfilled, so check y-direction further
stab(S12, S22, index_groups)
stab(S21, S11, index_groups)
stab(S12, S21, index_groups)
# recursive process
detect(left, center_pos, index_groups)
detect(right, num-center_pos, index_groups)
def stab(S1:list, S2:list, index_groups:list):
'''Check interval intersection in y-direction.'''
if not S1 or not S2: return
# sort
S1.sort(key=lambda item: item[1][1])
S2.sort(key=lambda item: item[1][1])
i, j = 0, 0
while i<len(S1) and j<len(S2):
m, a, _ = S1[i]
n, b, _ = S2[j]
if a[1] < b[1]:
k = j
while k<len(S2) and S2[k][1][1] < a[3]:
report_pair(int(m/2), int(S2[k][0]/2), index_groups)
k += 1
i += 1
else:
k = i
while k<len(S1) and S1[k][1][1] < b[3]:
report_pair(int(S1[k][0]/2), int(n/2), index_groups)
k += 1
j += 1
def report_pair(i:int, j:int, index_groups:list):
'''add pair (i,j) to adjacent list.'''
index_groups[i].add(j)
index_groups[j].add(i)进而，考虑改进算法的Rects类：
class Rects:
'''Collection of Rect instances.'''
def __init__(self):
self._rects = [] # type: list[Rect]
def append(self, instance): self._rects.append(instance)
def build_graph(self):
'''calculate rects connectivity by solving rectangle intersection problems.'''
num = len(self._rects)
graph = [set() for i in range(num)]
i_rect_x, i = [], 0
for rect in self._rects:
points = (rect.x0, rect.y0, rect.x1, rect.y1)
i_rect_x.append((i, points, rect.x0))
i_rect_x.append((i+1, points, rect.x1))
i += 2
i_rect_x.sort(key=lambda item: item[-1])
detect(i_rect_x, 2*num, graph) # use the improved algorithm
return graph
if __name__ == '__main__':
import random
from time import perf_counter
# initialize random coordinates
N = 50000
inputs = []
for i in range(N):
x0, y0 = random.random()*10000, random.random()*10000
x1 = x0 + random.random()*5
y1 = y0 + random.random()*5
inputs.append((x0,y0,x1,y1))
# initialize Rects
R = Rects()
for points in inputs: R.append(Rect(points))
# solving rectangle intersection
t0 = perf_counter()
graph = R.build_graph()
print(f'time: {perf_counter()-t0}')
# output:
# time: 1.2851152这样，处理50000个矩形相交问题的时间在2s以内。