def estimate_foreground_ml_cupy(
    input_image,
    input_alpha,
    regularization=1e-5,
    n_small_iterations=10,
    n_big_iterations=2,
    small_size=32,
    block_size=(32, 32),
    return_background=False,
    to_numpy=True,
):
    """See the :code:`estimate_foreground` method for documentation."""
    h0, w0, depth = input_image.shape

    assert depth == 3

    n = h0 * w0 * depth

    if isinstance(input_image, cp.ndarray):
        input_image = input_image.astype(cp.float32).ravel()
    else:
        input_image = cp.asarray(input_image.astype(np.float32).flatten())

    if isinstance(input_alpha, cp.ndarray):
        input_alpha = input_alpha.astype(cp.float32).ravel()
    else:
        input_alpha = cp.asarray(input_alpha.astype(np.float32).flatten())

    F_prev = cp.zeros(n, dtype=cp.float32)
    B_prev = cp.zeros(n, dtype=cp.float32)

    F = cp.zeros(n, dtype=cp.float32)
    B = cp.zeros(n, dtype=cp.float32)

    image_level = cp.zeros(n, dtype=cp.float32)
    alpha_level = cp.zeros(h0 * w0, dtype=cp.float32)

    n_levels = (max(w0, h0) - 1).bit_length()

    w_prev = 1
    h_prev = 1

    def resize_nearest(dst, src, w_src, h_src, w_dst, h_dst, depth):
        grid_size = (
            div_round_up(w_dst, block_size[0]),
            div_round_up(h_dst, block_size[1]),
        )

        _resize_nearest(
            grid_size,
            block_size,
            (
                dst,
                src,
                w_src,
                h_src,
                w_dst,
                h_dst,
                depth,
            ),
        )

    resize_nearest(F, input_image, w0, h0, w_prev, h_prev, depth)
    resize_nearest(B, input_image, w0, h0, w_prev, h_prev, depth)

    for i_level in range(n_levels + 1):
        w = round(w0 ** (i_level / n_levels))
        h = round(h0 ** (i_level / n_levels))

        resize_nearest(image_level, input_image, w0, h0, w, h, depth)
        resize_nearest(alpha_level, input_alpha, w0, h0, w, h, 1)

        resize_nearest(F_prev, F, w_prev, h_prev, w, h, depth)
        resize_nearest(B_prev, B, w_prev, h_prev, w, h, depth)

        # Do more iterations for low resolution.
        n_iter = n_big_iterations
        if min(w, h) <= small_size:
            n_iter = n_small_iterations

        grid_size = (div_round_up(w, block_size[0]), div_round_up(h, block_size[1]))

        for i_iter in range(n_iter):
            ml_iteration(
                grid_size,
                block_size,
                (
                    F,
                    B,
                    F_prev,
                    B_prev,
                    image_level,
                    alpha_level,
                    w,
                    h,
                    np.float32(regularization),
                ),
            )

            F_prev, F = F, F_prev
            B_prev, B = B, B_prev

        w_prev = w
        h_prev = h

    # Reshape back to original shape.
    F = F.reshape(h0, w0, depth)
    B = B.reshape(h0, w0, depth)

    # Convert to NumPy if requested
    if to_numpy:
        F_out = cp.asnumpy(F)
        B_out = cp.asnumpy(B)
    else:
        F_out = F
        B_out = B

    if return_background:
        return F_out, B_out

    return F_out
See the :code:`estimate_foreground` method for documentation.
pymatting/foreground/estimate_foreground_ml_cupy.py