ENH: Reduce memory use for convex_hull_image by checking hull inequal…

…ities in sequence (scikit-image#5020)

Co-authored-by: Juan Nunez-Iglesias <juan.nunez-iglesias@monash.edu>

skimage/morphology/convex_hull.py

@@ -18,6 +18,58 @@ def _offsets_diamond(ndim):
     return offsets
 
 
+def _check_coords_in_hull(gridcoords, hull_equations, tolerance):
+    r"""Checks all the coordinates for inclusiveness in the convex hull.
+
+    Parameters
+    ----------
+    gridcoords : (M, N) ndarray
+        Coordinates of ``N`` points in ``M`` dimensions.
+    hull_equations : (M, N) ndarray
+        Hyperplane equations of the facets of the convex hull.
+    tolerance : float
+        Tolerance when determining whether a point is inside the hull. Due
+        to numerical floating point errors, a tolerance of 0 can result in
+        some points erroneously being classified as being outside the hull.
+
+    Returns
+    -------
+    coords_in_hull : ndarray of bool
+        Binary 1D ndarray representing points in n-dimensional space
+        with value ``True`` set for points inside the convex hull.
+
+    Notes
+    -----
+    Checking the inclusiveness of coordinates in a convex hull requires
+    intermediate calculations of dot products which are memory-intensive.
+    Thus, the convex hull equations are checked individually with all
+    coordinates to keep within the memory limit.
+
+    References
+    ----------
+    .. [1] https://github.com/scikit-image/scikit-image/issues/5019
+
+    """
+    ndim, n_coords = gridcoords.shape
+    n_hull_equations = hull_equations.shape[0]
+    coords_in_hull = np.ones(n_coords, dtype=np.bool_)
+
+    # Pre-allocate arrays to cache intermediate results for reducing overheads
+    dot_array = np.empty(n_coords, dtype=np.float64)
+    test_ineq_temp = np.empty(n_coords, dtype=np.float64)
+    coords_single_ineq = np.empty(n_coords, dtype=np.bool_)
+
+    # A point is in the hull if it satisfies all of the hull's inequalities
+    for idx in range(n_hull_equations):
+        # Tests a hyperplane equation on all coordinates of volume
+        np.dot(hull_equations[idx, :ndim], gridcoords, out=dot_array)
+        np.add(dot_array, hull_equations[idx, ndim:], out=test_ineq_temp)
+        np.less(test_ineq_temp, tolerance, out=coords_single_ineq)
+        coords_in_hull *= coords_single_ineq
+
+    return coords_in_hull
+
+
 def convex_hull_image(image, offset_coordinates=True, tolerance=1e-10):
     """Compute the convex hull image of a binary image.
 
@@ -85,9 +137,9 @@ def convex_hull_image(image, offset_coordinates=True, tolerance=1e-10):
     else:
         gridcoords = np.reshape(np.mgrid[tuple(map(slice, image.shape))],
                                 (ndim, -1))
-        # A point is in the hull if it satisfies all of the hull's inequalities
-        coords_in_hull = np.all(hull.equations[:, :ndim].dot(gridcoords) +
-                                hull.equations[:, ndim:] < tolerance, axis=0)
+
+        coords_in_hull = _check_coords_in_hull(gridcoords,
+                                               hull.equations, tolerance)
         mask = np.reshape(coords_in_hull, image.shape)
 
     return mask