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