def find_fxpeaks(sp,
                 times_sigma_threshold,
                 minimum_threshold,
                 nwinwidth_initial,
                 nwinwidth_refined,
                 npix_avoid_border,
                 nbrightlines,
                 sigma_gaussian_filtering,
                 minimum_gaussian_filtering,
                 plottitle=None,
                 geometry=None,
                 debugplot=0):
    """Locate line peaks in array coordinates (from 0 to naxis1-1).

    Parameters
    ----------
    sp : 1d numpy array
        Spectrum to be wavelength calibrated.
    times_sigma_threshold : float
        Times robust sigma above the median to detect line peaks.
    minimum_threshold : float or None
            Minimum value of the threshold.
    nwinwidth_initial : int
        Width of the window where each peak must be found using
        the initial method (approximate)
    nwinwidth_refined : int
        Width of the window where each peak location will be
        refined.
    npix_avoid_border : int
            Number of pixels at the borders of the spectrum where peaks
            are not considered. If zero, the actual number will be
            given by nwinwidth_initial.
    nbrightlines : int or list of integers
        Maximum number of brightest lines to be employed in the
        wavelength calibration. If this value is 0, all the detected
        lines will be employed.
    sigma_gaussian_filtering : float
        Sigma of the gaussian filter to be applied to the spectrum in
        order to avoid problems with saturated lines. This filtering is
        skipped when this parameter is <= 0.
    minimum_gaussian_filtering : float
        Minimum pixel value to employ gaussian filtering. This value is
        employed only when sigma_gaussian_filtering is > 0.
    plottile : string
        Plot title.
    geometry : tuple (4 integers) or None
        x, y, dx, dy values employed to set the Qt backend geometry.
    debugplot : int
        Determines whether intermediate computations and/or plots
        are displayed. The valid codes are defined in
        numina.array.display.pause_debugplot.

    Returns
    -------
    fxpeaks : 1d numpy array
        Refined location of peaks in array index scale, i.e, from 0
        to naxis1 - 1.
    sxpeaks : 1d numpy array
        Line peak widths.

    """

    # spectrum dimension
    naxis1 = sp.shape[0]

    # apply gaussian filtering when requested
    if sigma_gaussian_filtering > 0:
        spf = ndimage.filters.gaussian_filter(
            sp,
            sigma=sigma_gaussian_filtering
        )
        lpreserve = sp < minimum_gaussian_filtering
        spf[lpreserve] = sp[lpreserve]
    else:
        spf = np.copy(sp)

    # threshold to search for peaks
    q50 = np.percentile(spf, q=50)
    sigma_g = robust_std(spf)
    threshold = q50 + times_sigma_threshold * sigma_g
    if abs(debugplot) >= 10:
        print("median....:", q50)
        print("robuts std:", sigma_g)
        print("threshold.:", threshold)
    if minimum_threshold > threshold:
        threshold = minimum_threshold
    if abs(debugplot) >= 10:
        print("minimum threshold:", minimum_threshold)
        print("final threshold..:", threshold)

    # initial location of the peaks (integer values)
    ixpeaks = find_peaks_spectrum(spf,
                                  nwinwidth=nwinwidth_initial,
                                  threshold=threshold)

    # select a maximum number of brightest lines in each region
    if len(nbrightlines) == 1 and nbrightlines[0] == 0:
        pass
    else:
        if abs(debugplot) >= 10:
            print('nbrightlines =', nbrightlines)
            print('ixpeaks in whole spectrum:\n', ixpeaks)
        region_size = (naxis1-1)/len(nbrightlines)
        ixpeaks_filtered = np.array([], dtype=int)
        for iregion, nlines_in_region in enumerate(nbrightlines):
            if nlines_in_region > 0:
                imin = int(iregion * region_size)
                imax = int((iregion + 1) * region_size)
                if iregion > 0:
                    imin += 1
                ixpeaks_region = \
                    ixpeaks[np.logical_and(ixpeaks >= imin, ixpeaks <= imax)]
                if len(ixpeaks_region) > 0:
                    peak_fluxes = spf[ixpeaks_region]
                    spos = peak_fluxes.argsort()
                    ixpeaks_tmp = ixpeaks_region[spos[-nlines_in_region:]]
                    ixpeaks_tmp.sort()  # in-place sort
                    if abs(debugplot) >= 10:
                        print('ixpeaks in region........:\n', ixpeaks_tmp)
                    ixpeaks_filtered = np.concatenate((ixpeaks_filtered,
                                                       ixpeaks_tmp))
        ixpeaks = ixpeaks_filtered
        if abs(debugplot) >= 10:
            print('ixpeaks filtered.........:\n', ixpeaks)

    # remove peaks too close to any of the borders of the spectrum
    if npix_avoid_border > 0:
        lok_ini = ixpeaks >= npix_avoid_border
        lok_end = ixpeaks <= naxis1 - 1 - npix_avoid_border
        ixpeaks = ixpeaks[lok_ini * lok_end]

    # refined location of the peaks (float values)
    fxpeaks, sxpeaks = refine_peaks_spectrum(spf, ixpeaks,
                                             nwinwidth=nwinwidth_refined,
                                             method="gaussian")

    # print peak location and width of fitted lines
    if abs(debugplot) >= 10:
        print(">>> Number of lines found:", len(fxpeaks))
        print("# line_number, channel, width")
        for i, (fx, sx) in enumerate(zip(fxpeaks, sxpeaks)):
            print(i, fx+1, sx)

    # display median spectrum and peaks
    if abs(debugplot) % 10 != 0:
        xplot = np.arange(1, naxis1 + 1, dtype=float)
        ax = ximplotxy(xplot, sp, show=False, geometry=geometry,
                       **{'label': 'original spectrum'})
        ax.set_xlabel('pixel (from 1 to NAXIS1)')
        ax.set_ylabel('counts')
        if plottitle is not None:
            ax.set_title(plottitle)
        if sigma_gaussian_filtering > 0:
            ax.plot(xplot, spf, label="filtered spectrum")
        ymin = sp.min()
        ymax = sp.max()
        dy = ymax - ymin
        ymin -= dy/20.
        ymax += dy/20.
        ax.set_ylim([ymin, ymax])
        # display threshold
        ax.axhline(y=threshold, color="black", linestyle="dotted",
                   label="detection threshold")
        # mark peak location
        ax.plot(ixpeaks + 1, spf[ixpeaks], 'bo', label="initial location")
        ax.plot(fxpeaks + 1, spf[ixpeaks], 'go', label="refined location")
        # legend
        ax.legend(numpoints=1)
        # show plot
        pause_debugplot(debugplot, pltshow=True)

    return fxpeaks, sxpeaks