A module for operations on 2-dim point sets (e.g., images)

imageproc.py

import numpy as np
import numpy.typing as np_typing


try:
    import cv2 as opencv
    has_opencv = True
except ModuleNotFoundError:
    has_opencv = False
    
    
__all__ = [
    "PI",
    "convex_hull",
    "is_contour_convex",
    "rect_to_bbox",
    "to_bounding_rect",
    "to_min_area_rect",
]


PI = np.pi


_Bbox = np_typing.NDArray[float]  # num_points x 2
_Contour = np_typing.NDArray[float]  # num_points x 2
_Hull = np_typing.NDArray[float]  # num_points x 1 x 2
_Point2d = tuple[float, float]  # x, y
_Rect = tuple[_Point2d, tuple[float, float], float]  # center x, center y, width, height, angle


def _rad_to_2d_rot_mat(*angle: float) -> np_typing.NDArray[float]:
    """It converts an angle in radians into a 2-dim rotation matrix"""
    batched = False
    
    if len(angle) > 1:
        batched = True

    angle = np.asarray(angle)
    
    rot_mat = np.column_stack((
        np.cos(angle), 
        np.cos(angle - PI / 2), 
        np.cos(angle + PI / 2), 
        np.cos(angle)
    )).reshape((-1, 2, 2))
    
    if not batched:
        rot_mat = rot_mat.squeeze()
        
    return rot_mat


def _is_contour_convex(contour: _Contour) -> bool:
    """It tests a contour convexity
    
    It tests whether a simple contour without self-intersections is convex.
    """
    num_points = len(contour)
    
    curr_point_idx = num_points - 1
    prev_point_idx = (curr_point_idx - 1 + num_points) % num_points
    
    curr_point = contour[curr_point_idx]
    prev_point = contour[prev_point_idx]
    
    diff = curr_point - prev_point
    orientation = 0

    for idx in range(num_points):
        prev_point = curr_point
        curr_point = contour[idx]
        
        diff_idx = curr_point - prev_point
        
        dxdy = diff_idx[0] * diff[1]
        dydx = diff_idx[1] * diff[0]

        orientation |= (1 if dydx > dxdy else (2 if dydx < dxdy else 3))
        
        if orientation == 3:
            return False

        diff = diff_idx

    return True
    
    
def _convex_hull(contour: _Contour, clockwise: bool = False) -> _Hull:
    """It finds the convex hull of a 2-dim point set"""
    num_points = len(contour)
    
    sorted_points = contour[np.lexsort((contour[:, 1], contour[:, 0]))]

    upper_hull = []
    
    for idx in range(num_points):
        while len(upper_hull) >= 2 and np.cross(
            upper_hull[-1] - upper_hull[-2], sorted_points[idx] - upper_hull[-2]
        ) <= 0:
            upper_hull.pop()

        upper_hull.append(sorted_points[idx])

    lower_hull = []
    
    for idx in range(num_points - 1, -1, -1):
        while len(lower_hull) >= 2 and np.cross(
            lower_hull[-1] - lower_hull[-2], sorted_points[idx] - lower_hull[-2]
        ) <= 0:
            lower_hull.pop()
    
        lower_hull.append(sorted_points[idx])
        
    convex_hull = upper_hull[:-1] + lower_hull[:-1]
    
    # TODO: The Graham scan algorithm for sorting by polar angle,
    # https://en.wikipedia.org/wiki/Graham_scan
    
    if not clockwise:
        convex_hull = convex_hull[::-1]
        
    convex_hull = np.asarray(convex_hull)
    convex_hull = convex_hull.reshape((-1, 1, 2))

    return convex_hull
    
    
def _to_bounding_rect(contour: _Contour) -> _Rect:
    """It calculates the up-right bounding rectangle of a 2-dim point set"""
    min_x = np.min(contour[:, 0]).astype(float)
    min_y = np.min(contour[:, 1]).astype(float)
    max_x = np.max(contour[:, 0]).astype(float)
    max_y = np.max(contour[:, 1]).astype(float)
    
    center_x = (min_x + max_x) / 2
    center_y = (min_y + max_y) / 2
    
    width = max_x - min_x; height = max_y - min_y
    
    center = (center_x, center_y)
    extent = (width, height)
    
    rect = (center, extent, 0.)
    return rect
    
def _opencv_to_bounding_rect(contour: _Contour) -> _Rect:
    """It calculates the up-right bounding rectangle of a 2-dim point set"""
    tl_x, tl_y, width, height = opencv.boundingRect(contour)
    
    center_x = tl_x + width / 2
    center_y = tl_y + height / 2
    
    center = (center_x, center_y)
    extent = (width, height)
    
    rect = (center, extent, 0.)
    return rect
    
    
def _to_min_area_rect(contour: _Contour) -> _Rect:
    """It finds a (possibly) rotated rectangle of the minimum area enclosing a 2-dim point set
    
    It finds a (possibly) rotated rectangle using the rotating calipers algorithm,
    https://wikipedia.org/wiki/Rotating_calipers
    """
    hull = _convex_hull(contour)
    hull = hull.reshape((-1, 2))

    hull_edges_x = np.diff(hull[:, 0])
    hull_edges_y = np.diff(hull[:, 1])

    candidates = np.arctan2(hull_edges_y, hull_edges_x)
    candidates = np.abs(candidates % (PI / 2))
    candidates = np.unique(candidates)
        
    candidates_rot_mat = _rad_to_2d_rot_mat(*candidates)
        
    rot_hull = np.dot(candidates_rot_mat, hull.T)
    
    min_x_hull = np.nanmin(rot_hull[:, 0], axis=1)
    min_y_hull = np.nanmin(rot_hull[:, 1], axis=1)
    max_x_hull = np.nanmax(rot_hull[:, 0], axis=1)
    max_y_hull = np.nanmax(rot_hull[:, 1], axis=1)

    candidates_area = (max_x_hull - min_x_hull) * (max_y_hull - min_y_hull)
    
    rect_idx = np.argmin(candidates_area)
    
    angle = candidates[rect_idx]
    angle = np.rad2deg(angle)
    
    min_x = min_x_hull[rect_idx]
    min_y = min_y_hull[rect_idx]
    max_x = max_x_hull[rect_idx]
    max_y = max_y_hull[rect_idx]
    
    rot_center_x = (min_x + max_x) / 2
    rot_center_y = (min_y + max_y) / 2

    rot_center = (rot_center_x, rot_center_y)

    center_x, \
    center_y = np.dot(rot_center, candidates_rot_mat[rect_idx])

    width = max_x - min_x; height = max_y - min_y
    
    center = (center_x, center_y)
    extent = (width, height)
    
    rect = (center, extent, angle)
    return rect

    
def _rect_to_bbox(rect: _Rect) -> _Bbox:
    """It finds the four clockwise vertices of a (possibly) rotated rectangle"""
    center, extent, angle = rect

    angle = np.deg2rad(angle)

    rot_mat = _rad_to_2d_rot_mat(angle)
    rot_center = np.dot(rot_mat, np.asarray(center))

    rot_min_x = rot_center[0] - extent[0] / 2
    rot_min_y = rot_center[1] - extent[1] / 2
    rot_max_x = rot_center[0] + extent[0] / 2
    rot_max_y = rot_center[1] + extent[1] / 2
    
    rot_bbox = np.asarray([
        [rot_min_x, rot_min_y],
        [rot_max_x, rot_min_y],
        [rot_max_x, rot_max_y],
        [rot_min_x, rot_max_y]
    ])

    bbox = np.dot(rot_bbox, rot_mat)
    
    # TODO: The Graham scan algorithm for sorting by polar angle,
    # https://en.wikipedia.org/wiki/Graham_scan

    return bbox

    
is_contour_convex = opencv.isContourConvex if has_opencv else _is_contour_convex
convex_hull = opencv.convexHull if has_opencv else _convex_hull
to_bounding_rect = _opencv_to_bounding_rect if has_opencv else _to_bounding_rect
to_min_area_rect = opencv.minAreaRect if has_opencv else _to_min_area_rect
rect_to_bbox = opencv.boxPoints if has_opencv else _rect_to_bbox