def make_pixel_distance(shape, xpix=None, ypix=None):
    """Fill a 2D array with dimensions `shape` with the distance of each
    pixel from a reference direction (xpix,ypix) in pixel coordinates.
    Pixel coordinates are defined such that (0,0) is located at the
    center of the corner pixel.

    """
    if np.isscalar(shape):
        shape = [shape, shape]

    if xpix is None:
        xpix = (shape[1] - 1.0) / 2.

    if ypix is None:
        ypix = (shape[0] - 1.0) / 2.

    dx = np.linspace(0, shape[1] - 1, shape[1]) - xpix
    dy = np.linspace(0, shape[0] - 1, shape[0]) - ypix
    dxy = np.zeros(shape)
    dxy += np.sqrt(dx[np.newaxis, :] ** 2 + dy[:, np.newaxis] ** 2)

    return dxy