目录
介绍
最小有向包围盒算法广泛应用于多个领域，包括：
- 计算几何：用于分析点集的边界特征。
- 图形学：用于碰撞检测和物体包围。
- 数据科学：在聚类分析、异常值检测等场景中用于计算数据的紧致边界。
本文提出了一种基于凸包和最小矩形算法的二维最小有向包围盒计算方法，并展示了如何通过代码实现这一算法。
主要步骤
- 凸包计算：使用
qhull2D方法计算输入点集的凸包。这一步可以确保我们找到点集的外部边界，以此作为最小矩形的候选边界。 - 最小矩形计算：利用
minBoundingRect算法计算最小有向包围盒。算法返回矩形的宽度、高度、角点坐标和面积。
程序结构如下：
调用入口：
输入为numpy格式的n*2数组。
代码
__init__.py
from .min_rect import *
from .qhull_2d import *
from .min_bounding_rect import *
min_bounding_rect.py
import numpy as np
import sys
def minBoundingRect(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.arctan2(edges[i, 1], edges[i, 0])
# Check for angles in 1st quadrant
edge_angles = np.abs(edge_angles % (np.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
# print("Testing", len(edge_angles), "possible rotations for bounding box... \n")
for i in range(len(edge_angles)):
# Create rotation matrix to shift points to baseline
R = np.array([[np.cos(edge_angles[i]), np.cos(edge_angles[i]-(np.pi/2))],
[np.cos(edge_angles[i]+(np.pi/2)), np.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
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.cos(angle), np.cos(angle-(np.pi/2))],
[np.cos(angle+(np.pi/2)), np.cos(angle)]])
# Project convex hull points onto rotated frame
proj_points = np.dot(R, np.transpose(hull_points_2d)) # 2x2 * 2xn
# 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]
# 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)
# 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 min_bbox[2], min_bbox[3], corner_points,area
min_rect.py
import numpy as np
import matplotlib.pyplot as plt
from .qhull_2d import qhull2D
from .min_bounding_rect import minBoundingRect
def min_area_rect(xy_points):
# Find convex hull
hull_points = qhull2D(xy_points)
# Reverse order of points, to match output from other qhull implementations
hull_points = hull_points[::-1]
# print('Convex hull points: \n', hull_points, "\n")
# Find minimum area bounding rectangle
width, height, corner_points,area = minBoundingRect(hull_points)
# # Visualization
plt.figure()
plt.plot(xy_points[:, 0], xy_points[:, 1], 'o', label='Points')
plt.plot(hull_points[:, 0], hull_points[:, 1], 'r--', lw=2, label='Convex Hull')
# Plot the bounding box
box = np.vstack([corner_points, corner_points[0]]) # Close the box by repeating the first point
plt.plot(box[:, 0], box[:, 1], 'g-', lw=2, label='Min Bounding Box')
plt.legend()
plt.xlabel('X')
plt.ylabel('Y')
plt.title('Convex Hull and Minimum Area Bounding Box')
plt.grid(True)
plt.axis('equal')
plt.show()
return width, height,corner_points,area
qhull_2d.py
from __future__ import division
from numpy import *
def link(a, b):
return concatenate((a, b[1:]))
def edge(a, b):
return concatenate(([a], [b]))
def qhull2D(sample):
def dome(sample, base, depth=0, max_depth=1000):
h, t = base
dists = dot(sample - h, dot(((0, -1), (1, 0)), (t - h)))
if len(dists) == 0:
return base
# Handle cases where all distances are very small
if all(abs(dists) < 1e-10):
return base
outer = sample[dists > 0]
# print("dists:",dists,"outer:",outer)
# Handle empty outer case
if len(outer) == 0:
return base
pivot_index = argmax(dists)
pivot = sample[pivot_index]
# Ensure depth does not exceed maximum allowed
if depth > max_depth:
return base
# Recursive case
left_dome = dome(outer, edge(h, pivot), depth + 1, max_depth)
right_dome = dome(outer, edge(pivot, t), depth + 1, max_depth)
return link(left_dome, right_dome)
if len(sample) > 2:
axis = sample[:, 0]
base = take(sample, [argmin(axis), argmax(axis)], axis=0)
left_dome = dome(sample, base)
right_dome = dome(sample, base[::-1])
return link(left_dome, right_dome)
else:
return sample
结果
