Bilinearly interpolate points from an image using NumPy or CuPy

bilinear_interpolate.py

def bilinear_interpolate(img, points, clipped_nan=True):
    """Bilinearly interpolate points from an image using NumPy or CuPy

    Args:
        img: Image of shape (Y, X) to interpolate from.
        points: array of shape (2, N) of (y, x) coordinates
        clipped_nan: If True, the value of coordinates outside the image shape 
            are set to nan. Otherwise they are clipped to the image edge.

    Returns:
        I: Array of bilinearly interpolated intensities

    """    
    y, x = points
    x0 = np.floor(x).astype(np.int32)
    x1 = x0 + 1
    y0 = np.floor(y).astype(np.int32)
    y1 = y0 + 1
    
    x0_clip = np.clip(x0, 0, img.shape[1])
    x1_clip = np.clip(x1, 0, img.shape[1])
    y0_clip = np.clip(y0, 0, img.shape[0])
    y1_clip = np.clip(y1, 0, img.shape[0])

    Ia = img[y0_clip, x0_clip]
    Ib = img[y1_clip, x0_clip]
    Ic = img[y0_clip, x1_clip]
    Id = img[y1_clip, x1_clip]

    wa = (x1-x) * (y1-y)
    wb = (x1-x) * (y-y0)
    wc = (x-x0) * (y1-y)
    wd = (x-x0) * (y-y0)
    
    I = (wa*Ia + wb*Ib + wc*Ic + wd*Id)
    
    if clipped_nan:
        mask = (x0 != x0_clip) | (x1 != x1_clip) | (y0 != y0_clip) | (y1 != y1_clip)
        I[mask] = np.nan
    return I