20x fast convex_hull_image for 2D

docs/index.md

@@ -48,6 +48,8 @@ pip install imops[numba]  # additionally install Numba backend
 
 ::: imops.morphology.distance_transform_edt
 
+::: imops.morphology.convex_hull_image
+
 ::: imops.measure.label
 
 ::: imops.measure.center_of_mass

imops/__version__.py

@@ -1 +1 @@
-__version__ = '0.9.0'
+__version__ = '0.9.1'

imops/morphology.py

@@ -1,22 +1,27 @@
+from itertools import product
 from typing import Callable, Tuple, Union
 from warnings import warn
 
 import numpy as np
 from edt import edt
 from scipy.ndimage import distance_transform_edt as scipy_distance_transform_edt, generate_binary_structure
 from scipy.ndimage._nd_image import euclidean_feature_transform
+from scipy.spatial import ConvexHull, QhullError
+from skimage._shared.utils import warn as skimage_warn
 from skimage.morphology import (
     binary_closing as scipy_binary_closing,
     binary_dilation as scipy_binary_dilation,
     binary_erosion as scipy_binary_erosion,
     binary_opening as scipy_binary_opening,
 )
+from skimage.util import unique_rows
 
 from .backend import BackendLike, Cython, Scipy, resolve_backend
 from .box import add_margin, box_to_shape, mask_to_box, shape_to_box
 from .compat import _ni_support
 from .crop import crop_to_box
 from .pad import restore_crop
+from .src._convex_hull import _grid_points_in_poly, _left_right_bounds, _offset_unique
 from .src._fast_morphology import (
     _binary_dilation as cython_fast_binary_dilation,
     _binary_erosion as cython_fast_binary_erosion,
@@ -517,3 +522,77 @@ def distance_transform_edt(
         return result[0]
 
     return None
+
+
+def convex_hull_image(image, offset_coordinates=True):
+    """
+    Fast convex hull of an image. Similar to skimage.morphology.convex_hull_image with include_borders=True
+
+    Parameters
+    ----------
+    image: np.ndarray
+        input image, must be 2D
+    offset_coordinates: bool
+        If True, a pixel at coordinate, e.g., (4, 7) will be represented by coordinates
+        (3.5, 7), (4.5, 7), (4, 6.5), and (4, 7.5).
+        This adds some “extent” to a pixel when computing the hull.
+
+    Returns
+    -------
+    output: np.ndarray
+        resulting convex hull of the input image
+
+    Examples
+    --------
+    ```python
+    chull = convex_hull_image(x)
+    ```
+    """
+
+    ndim = image.ndim
+    if ndim != 2:
+        raise RuntimeError(f'convex_hull_image is currently implemented only for 2D arrays, got {ndim}D array')
+
+    if np.count_nonzero(image) == 0:
+        warn(
+            "Input image is entirely zero, no valid convex hull. " "Returning empty image",
+            UserWarning,
+        )
+        return np.zeros(image.shape, dtype=bool)
+
+    # In 2D, we do an optimisation by choosing only pixels that are
+    # the starting or ending pixel of a row or column.  This vastly
+    # limits the number of coordinates to examine for the virtual hull.
+    image = np.ascontiguousarray(image, dtype=np.uint8)
+
+    coords = _left_right_bounds(image)
+
+    if offset_coordinates:
+        coords = _offset_unique(coords)
+    else:
+        coords = unique_rows(coords)
+
+    # Find the convex hull
+    try:
+        hull = ConvexHull(coords)
+    except QhullError as err:
+        skimage_warn(f"Failed to get convex hull image. " f"Returning empty image, see error message below:\n" f"{err}")
+        return np.zeros(image.shape, dtype=bool)
+
+    vertices = hull.points[hull.vertices]
+
+    # return vertices
+
+    # If 2D, use fast Cython function to locate convex hull pixels
+    labels = _grid_points_in_poly(
+        np.ascontiguousarray(vertices[:, 0], dtype=np.float32),
+        np.ascontiguousarray(vertices[:, 1], dtype=np.float32),
+        image.shape[0],
+        image.shape[1],
+        len(vertices),
+    )
+
+    # edge points and vertices are included
+    mask = labels.view(bool)
+
+    return mask

imops/src/_convex_hull.pyx

@@ -0,0 +1,261 @@
+#cython: cdivision=True
+#cython: boundscheck=False
+#cython: nonecheck=False
+#cython: wraparound=False
+import numpy as np
+
+cimport cython
+cimport numpy as cnp
+from libc.math cimport ceilf, floorf
+
+
+cnp.import_array()
+
+FP_BOUND_DTYPE = np.dtype([
+    ('lb', np.float32),
+    ('rb', np.float32),
+    ('assigned', np.uint8)
+])
+
+INT_BOUND_DTYPE = np.dtype([
+    ('lb', np.int32),
+    ('rb', np.int32),
+    ('assigned', np.uint8)
+])
+
+
+def _grid_points_in_poly(float[:] vx, float[:] vy, Py_ssize_t M, Py_ssize_t N, Py_ssize_t nr_verts):
+    cdef intBound tmp_int_bound
+    cdef Py_ssize_t m, n, i
+    cdef float prev_x, prev_y, curr_x, curr_ys
+    cdef float tmp_from_x, tmp_from_y, tmp_to_x, tmp_to_y
+    cdef float lerp_t, bound_y
+    cdef Py_ssize_t x_set, x_start, x_stop
+
+    cdef cnp.ndarray[dtype=cnp.uint8_t, ndim=2, mode="c"] out = \
+         np.zeros((M, N), dtype=np.uint8)
+
+    cdef fpBound[:] fpBounds = np.empty(M, dtype=FP_BOUND_DTYPE)
+    cdef intBound[:] intBounds = np.empty(M, dtype=INT_BOUND_DTYPE)
+
+    for i in range(M):
+        fpBounds[i].assigned = False
+        fpBounds[i].lb = float('inf')
+        fpBounds[i].rb = -1
+        intBounds[i].assigned = False
+
+    prev_x = vx[nr_verts - 1]
+    prev_y = vy[nr_verts - 1]
+
+    # algorithm relies on vertex validity and counterclockwise orientation of the vertices
+    for i in range(nr_verts):
+        curr_x = vx[i]
+        curr_y = vy[i]
+
+        if prev_x == curr_x:
+            x_set = <int>(floorf(prev_x) if prev_y < curr_y else ceilf(prev_x))
+
+            fpBounds[x_set].assigned = True
+            fpBounds[x_set].lb = min(fpBounds[x_set].lb, prev_y, curr_y)
+            fpBounds[x_set].rb = max(fpBounds[x_set].rb, prev_y, curr_y)
+        else:
+            if prev_x < curr_x:
+                tmp_from_x = prev_x
+                tmp_from_y = prev_y
+                tmp_to_x = curr_x
+                tmp_to_y = curr_y
+            else:
+                tmp_from_x = curr_x
+                tmp_from_y = curr_y
+                tmp_to_x = prev_x
+                tmp_to_y = prev_y
+
+            # vertices are treated as points on image, so include x_stop
+            x_start = <int>ceilf(tmp_from_x)
+            x_stop = <int>floorf(tmp_to_x + 1)
+
+            for x_set in range(x_start, x_stop):
+                lerp_t = (x_set - tmp_from_x) / (tmp_to_x - tmp_from_x)
+                bound_y = lerp(tmp_from_y, tmp_to_y, lerp_t)
+
+                fpBounds[x_set].assigned = True
+                fpBounds[x_set].lb = min(fpBounds[x_set].lb, bound_y)
+                fpBounds[x_set].rb = max(fpBounds[x_set].rb, bound_y)
+
+        prev_x = curr_x
+        prev_y = curr_y
+
+    # bounds are computed as point interpolation
+    # so bounds must be valid indices for out array 
+    for m in range(M):
+        intBounds[m] = intify(fpBounds[m], 0, N - 1)
+
+    for m in range(M):
+        tmp_int_bound = intBounds[m]
+
+        if tmp_int_bound.assigned:
+            # Do not forget to fill right bound
+            for n in range(tmp_int_bound.lb, tmp_int_bound.rb + 1):
+                out[m, n] = True
+
+    return out
+
+
+# TODO: maybe use round instead of floorf and ceilf?
+cdef inline intBound intify(fpBound bound, Py_ssize_t min_idx, Py_ssize_t max_idx):
+    if bound.assigned:
+        return intBound(lb = max(min_idx, <int>floorf(bound.lb)), rb = min(max_idx, <int>ceilf(bound.rb)), assigned=True)
+
+    return intBound(lb=0, rb=0, assigned=False)
+
+
+cdef packed struct fpBound:
+    float lb
+    float rb
+    unsigned char assigned
+
+
+cdef packed struct intBound:
+    int lb
+    int rb
+    unsigned char assigned
+
+
+cdef inline float lerp(float y0, float y1, float t):
+    return y0 * (1 - t) + y1 * t
+
+
+cpdef _left_right_bounds(cnp.uint8_t[:,:] image):
+    cdef Py_ssize_t i, j, M = image.shape[0], N = image.shape[1], curr_pos = 0, left, right
+    cdef cnp.ndarray[dtype=int, ndim=2, mode="c"] left_right_bounds = np.zeros((2 * M, 2), dtype=np.int32)
+    cdef unsigned char found = False
+
+    for i in range(M):
+        found = False
+
+        for j in range(N):
+            if image[i, j]:
+                left = j
+                found = True
+                break
+
+        for j in range(N):
+            if image[i, N - 1 - j]:
+                right = N - 1 - j
+                found = True
+                break
+
+        if found:
+            left_right_bounds[2 * curr_pos, 0] = i
+            left_right_bounds[2 * curr_pos, 1] = left
+            left_right_bounds[2 * curr_pos + 1, 0] = i
+            left_right_bounds[2 * curr_pos + 1, 1] = right
+
+            curr_pos += 1
+
+    return np.ascontiguousarray(left_right_bounds[: 2 * curr_pos, :])
+
+
+cdef inline int set_unique_curr(float* expanded_bounds, int x, int l, int r):
+    if l == r:
+        expanded_bounds[0] = x
+        expanded_bounds[1] = l - 0.5
+
+        expanded_bounds[2] = x - 0.5
+        expanded_bounds[3] = l
+
+        expanded_bounds[4] = x
+        expanded_bounds[5] = l + 0.5
+
+        return 3
+    elif r == l + 1:
+        expanded_bounds[0] = x
+        expanded_bounds[1] = l - 0.5
+
+        expanded_bounds[2] = x - 0.5
+        expanded_bounds[3] = l
+
+        expanded_bounds[4] = x
+        expanded_bounds[5] = l + 0.5
+
+        expanded_bounds[6] = x - 0.5
+        expanded_bounds[7] = r
+
+        expanded_bounds[8] = x
+        expanded_bounds[9] = r + 0.5
+
+        return 5
+
+    else:
+        expanded_bounds[0] = x
+        expanded_bounds[1] = l - 0.5
+
+        expanded_bounds[2] = x - 0.5
+        expanded_bounds[3] = l
+
+        expanded_bounds[4] = x
+        expanded_bounds[5] = l + 0.5
+
+        expanded_bounds[6] = x
+        expanded_bounds[7] = r - 0.5
+
+        expanded_bounds[8] = x - 0.5
+        expanded_bounds[9] = r
+
+        expanded_bounds[10] = x
+        expanded_bounds[11] = r + 0.5
+
+        return 6
+
+
+cpdef _offset_unique(int[:,:] left_right_bounds):
+    cdef Py_ssize_t N = left_right_bounds.shape[0], i, curr_pos = 0
+    cdef cnp.ndarray[dtype=float, ndim=2, mode="c"] expanded_bounds = np.zeros((4 * N, 2), dtype=np.float32)
+
+    cdef int x_l_prev, y_l_prev, x_r_prev, y_r_prev, x_l_curr, y_l_curr, x_r_curr, y_r_curr, shift
+
+    x_l_prev = left_right_bounds[0, 0]
+    y_l_prev = left_right_bounds[0, 1]
+    x_r_prev = left_right_bounds[1, 0]
+    y_r_prev = left_right_bounds[1, 1]
+
+    curr_pos += set_unique_curr(&expanded_bounds[0, 0], x_l_prev, y_l_prev, y_r_prev)
+
+    for i in range(1, N // 2):
+        x_l_curr = left_right_bounds[2 * i, 0]
+        y_l_curr = left_right_bounds[2 * i, 1]
+        x_r_curr = left_right_bounds[2 * i + 1, 0]
+        y_r_curr = left_right_bounds[2 * i + 1, 1]
+
+        curr_pos += set_unique_curr(&expanded_bounds[curr_pos, 0], x_l_curr, y_l_curr, y_r_curr)
+
+        if x_l_prev + 1 == x_l_curr and (y_l_prev == y_l_curr or y_l_prev == y_r_curr):
+            pass
+        else:
+            expanded_bounds[curr_pos, 0] = x_l_prev + 0.5
+            expanded_bounds[curr_pos, 1] = y_l_prev
+            curr_pos += 1
+
+        if x_l_prev + 1 == x_l_curr and (y_r_prev == y_l_curr or y_r_prev == y_r_curr) and (y_r_prev != y_l_prev):
+            pass
+        else:
+            expanded_bounds[curr_pos, 0] = x_l_prev + 0.5
+            expanded_bounds[curr_pos, 1] = y_r_prev
+            curr_pos += 1
+
+        x_l_prev = x_l_curr
+        y_l_prev = y_l_curr
+        x_r_prev = x_r_curr
+        y_r_prev = y_r_curr
+
+    expanded_bounds[curr_pos, 0] = x_l_prev + 0.5
+    expanded_bounds[curr_pos, 1] = y_l_prev
+    curr_pos += 1
+
+    if y_r_prev != y_l_prev:
+        expanded_bounds[curr_pos, 0] = x_l_prev + 0.5
+        expanded_bounds[curr_pos, 1] = y_r_prev
+        curr_pos += 1
+
+
+    return np.ascontiguousarray(expanded_bounds[:curr_pos, :])