<SYSTEM_TASK:>
Solve the following problem using Python, implementing the functions described below, one line at a time
<END_TASK>
<USER_TASK:>
Description:
 def find_xy_peak(img, center=None, sigma=3.0):
 """ Find the center of the peak of offsets """

    # find level of noise in histogram
    istats = imagestats.ImageStats(img.astype(np.float32), nclip=1,
                                   fields='stddev,mode,mean,max,min')
    if istats.stddev == 0.0:
        istats = imagestats.ImageStats(img.astype(np.float32),
                                       fields='stddev,mode,mean,max,min')
    imgsum = img.sum()

    # clip out all values below mean+3*sigma from histogram
    imgc = img[:, :].copy()
    imgc[imgc < istats.mode + istats.stddev * sigma] = 0.0
    # identify position of peak
    yp0, xp0 = np.where(imgc == imgc.max())

    # Perform bounds checking on slice from img
    ymin = max(0, int(yp0[0]) - 3)
    ymax = min(img.shape[0], int(yp0[0]) + 4)
    xmin = max(0, int(xp0[0]) - 3)
    xmax = min(img.shape[1], int(xp0[0]) + 4)
    # take sum of at most a 7x7 pixel box around peak
    xp_slice = (slice(ymin, ymax),
                slice(xmin, xmax))
    yp, xp = ndimage.measurements.center_of_mass(img[xp_slice])
    if np.isnan(xp) or np.isnan(yp):
        xp = 0.0
        yp = 0.0
        flux = 0.0
        zpqual = None
    else:
        xp += xp_slice[1].start
        yp += xp_slice[0].start

        # compute S/N criteria for this peak: flux/sqrt(mean of rest of array)
        flux = imgc[xp_slice].sum()
        delta_size = float(img.size - imgc[xp_slice].size)
        if delta_size == 0:
            delta_size = 1
        delta_flux = float(imgsum - flux)
        if flux > imgc[xp_slice].max():
            delta_flux = flux - imgc[xp_slice].max()
        else:
            delta_flux = flux
        zpqual = flux / np.sqrt(delta_flux / delta_size)
        if np.isnan(zpqual) or np.isinf(zpqual):
            zpqual = None

        if center is not None:
            xp -= center[0]
            yp -= center[1]
        flux = imgc[xp_slice].max()

    del imgc
    return xp, yp, flux, zpqual