kneeliverse.convex_hull

  1# coding: utf-8
  2
  3'''
  4The following module provides optimized methods
  5for the computation of a convex hullfor 2D.
  6'''
  7
  8__author__ = 'Mário Antunes'
  9__version__ = '1.0'
 10__email__ = 'mario.antunes@ua.pt'
 11__status__ = 'Development'
 12__license__ = 'MIT'
 13__copyright__ = '''
 14Copyright (c) 2021-2023 Stony Brook University
 15Copyright (c) 2021-2023 The Research Foundation of SUNY
 16'''
 17
 18
 19import functools
 20import numpy as np
 21
 22
 23def _ccw(a:np.ndarray, b:np.ndarray, c:np.ndarray) -> float:
 24    """
 25    Check if three points make a counter-clockwise turn.
 26
 27    Args:
 28        a (np.ndarray): numpy array with a single point (x, y)
 29        b (np.ndarray): numpy array with a single point (x, y)
 30        c (np.ndarray): numpy array with a single point (x, y)
 31
 32    Returns:
 33        float: $ \\gt 0$ if counter-clockwise; $ \\lt 0$ if clockwise; $ = 0$ if collinear
 34    """
 35    return (b[0] - a[0]) * (c[1] - a[1]) - (c[0] - a[0]) * (b[1] - a[1])
 36
 37
 38def _dist_points(pi, pj) -> np.ndarray:
 39    """
 40    Compute the distance between two points ($p_i$, $p_j$)
 41    """
 42    return np.linalg.norm(pi - pj)
 43
 44
 45def _compare_points(p0, pi, pj) -> int:
 46    o = _ccw(p0, pi, pj)
 47    if o == 0:
 48        return -1 if _dist_points(p0, pj) >= _dist_points(p0, pi) else 1
 49    else:
 50        return 1 if o > 0 else -1
 51
 52
 53def _sort_points(points:np.ndarray) -> np.ndarray:
 54    """
 55    Sort the points based on the polar angle to the first point.
 56    
 57    The functions does not calculate the angle, instead, it calculate the relative orientation 
 58    of two points to find out which point makes the larger angle. 
 59
 60    Args:
 61        points (np.ndarray): numpy array with the points (x, y)
 62
 63    Returns:
 64        np.ndarray: the sorted points
 65    """
 66    
 67    # find the smaller point and respective index
 68    p0 = min(points, key=lambda p: (p[0], p[1]))
 69    p0_idx = np.where(np.all(points==p0,axis=1))[0][0]
 70    
 71    # create a mask to select all points except the smaller point
 72    mask = np.full(len(points), True)
 73    mask[p0_idx] = False
 74    temp_list = points[mask]
 75
 76    # return a copy of the points sorted
 77    return [p0] + sorted(temp_list, key=functools.cmp_to_key(lambda x, y: _compare_points(p0, x, y)))
 78
 79
 80def graham_scan(points:np.ndarray) -> np.ndarray:
 81    """
 82    Returns an array of indexes that make up the convex hull surrounding the points.
 83    
 84    Args:
 85        points (np.ndarray): numpy array with the points (x, y)
 86
 87    Returns:
 88        np.ndarray: the indexes of the hull points
 89    """
 90    # convex_hull is a stack of points beginning with the leftmost point.
 91    stack = []
 92    sorted_points = _sort_points(points)
 93
 94    if len(sorted_points) >= 3:
 95        stack.append(sorted_points[0])
 96        stack.append(sorted_points[1])
 97        stack.append(sorted_points[2])
 98
 99    for i in range(3, len(sorted_points)):
100        p = sorted_points[i]
101        
102        while _ccw(stack[-2], stack[-1], p) >= 0:
103            stack.pop()
104        stack.append(p)
105
106    # the stack is now a representation of the convex hull, return it.
107    # convert the points into indexes
108    hull = [np.where(np.all(points==p,axis=1))[0][0] for p in stack]
109    return np.array(hull)
110
111
112def graham_scan_lower(points:np.ndarray) -> np.ndarray:
113    """
114    Returns an array of indexes that make up the lower convex hull.
115    
116    Args:
117        points (np.ndarray): numpy array with the points (x, y)
118
119    Returns:
120        np.ndarray: the indexes of the hull points
121    """
122    #sorted_points = sort_points(points)
123    # Add p0 and p1 to the stack
124    stack = [0, 1]
125
126    for i in range(2, len(points)):
127        # if we turn clockwise to reach this point, pop the last point from the stack, else, append this point to it.
128        while len(stack) > 1 and ccw(points[stack[-2]], points[stack[-1]], points[i]) <= 0:
129            stack.pop()
130        stack.append(i)
131    # the stack is now a representation of the convex hull, return it.
132    # convert the points into indexes
133    #x = points[:, 0]
134    #hull = [np.where(x == p[0])[0][0] for p in stack]
135    return np.array(stack)
136
137
138def graham_scan_upper(points:np.ndarray) -> np.ndarray:
139    """
140    Returns an array of indexes that make up the upper convex hull.
141    
142    Args:
143        points (np.ndarray): numpy array with the points (x, y)
144
145    Returns:
146        np.ndarray: the indexes of the hull points
147    """
148    #sorted_points = sort_points(points)
149    # Add p0 and p1 to the stack
150    stack = [0, 1]
151
152    for i in range(2, len(points)):
153        # if we turn clockwise to reach this point, pop the last point from the stack, else, append this point to it.
154        while len(stack) > 1 and ccw(points[i], points[stack[-1]], points[stack[-2]]) <= 0:
155            stack.pop()
156        stack.append(i)
157    # the stack is now a representation of the convex hull, return it.
158    # convert the points into indexes
159    #x = points[:, 0]
160    #hull = [np.where(x == p[0])[0][0] for p in stack]
161    return np.array(stack)