Get the peak position with subpixel precision

    Parameters
    ----------
    im: 2d array
        The image containing a peak
    wrap: boolean, defaults False
        True if the image reoresents a torric world

    Returns
    -------
    [y,x]: 2 numbers
        The position of the highest peak with subpixel precision

    Notes
    -----
    This is a bit hacky and could be improved

def get_peak_pos(im, wrap=False):
    """Get the peak position with subpixel precision

    Parameters
    ----------
    im: 2d array
        The image containing a peak
    wrap: boolean, defaults False
        True if the image reoresents a torric world

    Returns
    -------
    [y,x]: 2 numbers
        The position of the highest peak with subpixel precision

    Notes
    -----
    This is a bit hacky and could be improved
    """
    im = np.asarray(im)
    # remove invalid values (assuming im>0)
    im[np.logical_not(np.isfinite(im))] = 0
    # remove mean
    im = im - im.mean()
    # get maximum value
    argmax = im.argmax()
    dsize = im.size
    # get cut value (30% biggest peak)
    # TODO: choose less random value
    cut = .3 * im[argmax]
    # isolate peak
    peak = im > cut
    peak, __ = label(peak)
    # wrap border
    if wrap and peak[0] != 0 and peak[-1] != 0 and peak[0] != peak[-1]:
        peak[peak == peak[-1]] = peak[0]
    # extract peak
    peak = peak == peak[argmax]

    # get values along X and Y
    X = np.arange(dsize)[peak]
    Y = im[peak]
    # wrap border
    if wrap:
        # wrap X values d
        X[X > dsize // 2] -= dsize

    # remove argmax as in X**4 X should be small
    offset = X[Y == Y.max()][0]
    X -= offset

    # We want to fit in a radius of 3 around the center
    Y = Y[abs(X) < 3]
    X = X[abs(X) < 3]

    # if>2, use fit_log
    if peak.sum() > 2:
        ret, __ = gauss_fit_log(X, Y)
        # if fails, use center_of_mass
        if ret is np.nan:
            ret = center_of_mass(X, Y)
    elif peak.sum() > 1:
        # If only 2 pixel, gauss fitting is imposible, use center_of_mass
        ret = center_of_mass(X, Y)
    else:
        # 1 px peak is easy
        ret = X[0]

    """
    import matplotlib.pyplot as plt
    plt.figure()
    plt.plot(X,Y,'x',label='im')
    plt.plot([ret,ret],[1,Y.max()],label='logfit')
    plt.plot([X.min(),X.max()],[cut,cut])
    plt.plot([X.min(),X.max()],[im.std(),im.std()])
    #"""

    return ret + offset