def interpolate_bilinear(p, img, img_size, transpose=True):
    """
    Bilinear interpolation.

    Parameters:
    ----------
    p : NDArray
        Float coordinates.
    img : NDArray
        original image.
    img_size : tuple of two inst
        Image size.
    transpose : bool, default True
        Whether do transpose of the output against input.

    Returns
    -------
    zz : NDArray
       Interpolated values.
    """
    p0 = p.floor().astype(int)
    p1 = p0 + 1
    x = p.slice_axis(axis=1, begin=0, end=1).squeeze(axis=1)
    y = p.slice_axis(axis=1, begin=1, end=2).squeeze(axis=1)
    x0 = p0.slice_axis(axis=1, begin=0, end=1).squeeze(axis=1)
    x1 = p1.slice_axis(axis=1, begin=0, end=1).squeeze(axis=1)
    y0 = p0.slice_axis(axis=1, begin=1, end=2).squeeze(axis=1)
    y1 = p1.slice_axis(axis=1, begin=1, end=2).squeeze(axis=1)
    x0 = x0.clip(0, img_size[0] - 1)
    x1 = x1.clip(0, img_size[0] - 1)
    y0 = y0.clip(0, img_size[1] - 1)
    y1 = y1.clip(0, img_size[1] - 1)
    z00 = img[:, (x0), (y0)].T
    z01 = img[:, (x0), (y1)].T
    z10 = img[:, (x1), (y0)].T
    z11 = img[:, (x1), (y1)].T
    x0 = x0.astype(p.dtype)
    x1 = x1.astype(p.dtype)
    y0 = y0.astype(p.dtype)
    y1 = y1.astype(p.dtype)
    w00 = ((x - x0) * (y - y0)).expand_dims(axis=1)
    w01 = ((x - x0) * (y1 - y)).expand_dims(axis=1)
    w10 = ((x1 - x) * (y - y0)).expand_dims(axis=1)
    w11 = ((x1 - x) * (y1 - y)).expand_dims(axis=1)
    zz = z00 * w11 + z10 * w10 + z01 * w01 + z11 * w00
    if not transpose:
        zz = zz.T
    return zz
