def subdivide(im, divs=2):
    r"""
    Returns slices into an image describing the specified number of sub-arrays.

    This function is useful for performing operations on smaller images for
    memory or speed.  Note that for most typical operations this will NOT work,
    since the image borders would cause artifacts (e.g. ``distance_transform``)

    Parameters
    ----------
    im : ND-array
        The image of the porous media

    divs : scalar or array_like
        The number of sub-divisions to create in each axis of the image.  If a
        scalar is given it is assumed this value applies in all dimensions.

    Returns
    -------
    slices : 1D-array
        A 1-D array containing slice objects for indexing into ``im`` that
        extract the sub-divided arrays.

    Notes
    -----
    This method uses the
    `array_split package <https://github.com/array-split/array_split>`_ which
    offers the same functionality as the ``split`` method of Numpy's ND-array,
    but supports the splitting multidimensional arrays in all dimensions.

    Examples
    --------
    >>> import porespy as ps
    >>> import matplotlib.pyplot as plt
    >>> im = ps.generators.blobs(shape=[200, 200])
    >>> s = ps.tools.subdivide(im, divs=[2, 2])

    ``s`` contains an array with the shape given by ``divs``.  To access the
    first and last quadrants of ``im`` use:
    >>> print(im[tuple(s[0, 0])].shape)
    (100, 100)
    >>> print(im[tuple(s[1, 1])].shape)
    (100, 100)

    It can be easier to index the array with the slices by applying ``flatten``
    first:
    >>> s_flat = s.flatten()
    >>> for i in s_flat:
    ...     print(im[i].shape)
    (100, 100)
    (100, 100)
    (100, 100)
    (100, 100)
    """
    # Expand scalar divs
    if isinstance(divs, int):
        divs = [divs for i in range(im.ndim)]
    s = shape_split(im.shape, axis=divs)
    return s
r"""
    Returns slices into an image describing the specified number of sub-arrays.

    This function is useful for performing operations on smaller images for
    memory or speed.  Note that for most typical operations this will NOT work,
    since the image borders would cause artifacts (e.g. ``distance_transform``)

    Parameters
    ----------
    im : ND-array
        The image of the porous media

    divs : scalar or array_like
        The number of sub-divisions to create in each axis of the image.  If a
        scalar is given it is assumed this value applies in all dimensions.

    Returns
    -------
    slices : 1D-array
        A 1-D array containing slice objects for indexing into ``im`` that
        extract the sub-divided arrays.

    Notes
    -----
    This method uses the
    `array_split package <https://github.com/array-split/array_split>`_ which
    offers the same functionality as the ``split`` method of Numpy's ND-array,
    but supports the splitting multidimensional arrays in all dimensions.

    Examples
    --------
    >>> import porespy as ps
    >>> import matplotlib.pyplot as plt
    >>> im = ps.generators.blobs(shape=[200, 200])
    >>> s = ps.tools.subdivide(im, divs=[2, 2])

    ``s`` contains an array with the shape given by ``divs``.  To access the
    first and last quadrants of ``im`` use:
    >>> print(im[tuple(s[0, 0])].shape)
    (100, 100)
    >>> print(im[tuple(s[1, 1])].shape)
    (100, 100)

    It can be easier to index the array with the slices by applying ``flatten``
    first:
    >>> s_flat = s.flatten()
    >>> for i in s_flat:
    ...     print(im[i].shape)
    (100, 100)
    (100, 100)
    (100, 100)
    (100, 100)