当前位置：首页 > article >正文

最小有向包围盒——2D平面

article 2024/12/4 17:57:04

介绍

主要步骤

代码

__init__.py

min_bounding_rect.py

min_rect.py

qhull_2d.py

结果

介绍

最小有向包围盒算法广泛应用于多个领域，包括：

计算几何：用于分析点集的边界特征。
图形学：用于碰撞检测和物体包围。
数据科学：在聚类分析、异常值检测等场景中用于计算数据的紧致边界。

本文提出了一种基于凸包和最小矩形算法的二维最小有向包围盒计算方法，并展示了如何通过代码实现这一算法。

主要步骤

凸包计算：使用 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