docs: def median_spectrum_from_rectified_image(self, slitlet2d_rect,
                                               sigma_gaussian_filtering=0,
                                               nwinwidth_initial=5,
                                               nwinwidth_refined=7,
                                               times_sigma_threshold=5,
                                               minimum_threshold=None,
                                               npix_avoid_border=0,
                                               nbrightlines=None):
          """Median spectrum and line peaks from rectified image.
  
          In order to avoid the line ghosts, the line peaks are identified
          independently in the upper and lower halves of the rectified
          image. The final peaks correspond to lines that appear in both
          spectra.
  
          Parameters
          ----------
          slitlet2d_rect : numpy array
              Rectified slitlet image.
          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.
          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.
          times_sigma_threshold : float
              Times (robust) sigma above the median of the image to set
              the minimum threshold when searching for line peaks.
          minimum_threshold : float or None
              Minimum value of the threshold.
          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 : list
              List with maximum number of brightest lines to be employed
              in the wavelength calibration. The length of the list
              indicates the number of equal-size subintervals in the
              wavelength calibration direction to be considered.
              If this value is [0] or None, all the detected lines will be
              employed.
  
          Returns
          -------
          sp0 : numpy array
              Median spectrum.
          fxpeaks : numpy array
              Refined location of arc lines (in array index scale).
  
          """
  
          # protections
          naxis2, naxis1 = slitlet2d_rect.shape
          if naxis1 != self.bb_nc2_orig - self.bb_nc1_orig + 1:
              raise ValueError("Unexpected slitlet2d_rect naxis1")
          if naxis2 != self.bb_ns2_orig - self.bb_ns1_orig + 1:
              raise ValueError("Unexpected slitlet2d_rect naxis2")
  
          # lower, middle and upper spectrum trails
          ylower_line = \
              self.list_spectrails[self.i_lower_spectrail].y_rectified
          ymiddle_line = \
              self.list_spectrails[self.i_middle_spectrail].y_rectified
          yupper_line = \
              self.list_spectrails[self.i_upper_spectrail].y_rectified
  
          ilower = int(ylower_line + 0.5) - self.bb_ns1_orig
          imiddle = int(ymiddle_line + 0.5) - self.bb_ns1_orig
          iupper = int(yupper_line + 0.5) - self.bb_ns1_orig
  
          # median spectra using different image regions
          sp0_ini = np.median(slitlet2d_rect[ilower:(iupper + 1), :], axis=0)
          sp1_ini = np.median(slitlet2d_rect[ilower:(imiddle + 1), :], axis=0)
          sp2_ini = np.median(slitlet2d_rect[imiddle:(iupper + 1), :], axis=0)
  
          # gaussian filtering when requested (to avoid line saturation)
          if sigma_gaussian_filtering > 0:
              sp0 = ndimage.filters.gaussian_filter(
                  sp0_ini,
                  sigma=sigma_gaussian_filtering
              )
              sp1 = ndimage.filters.gaussian_filter(
                  sp1_ini,
                  sigma=sigma_gaussian_filtering
              )
              sp2 = ndimage.filters.gaussian_filter(
                  sp2_ini,
                  sigma=sigma_gaussian_filtering
              )
          else:
              sp0 = np.copy(sp0_ini)
              sp1 = np.copy(sp1_ini)
              sp2 = np.copy(sp2_ini)
  
          # compute threshold
          q25, q50, q75 = np.percentile(sp0, q=[25.0, 50.0, 75.0])
          sigma_g = 0.7413 * (q75 - q25)  # robust standard deviation
          threshold = q50 + times_sigma_threshold * sigma_g
          if abs(self.debugplot) >= 10:
              print("median...........:", q50)
              print("robuts std.......:", sigma_g)
              print("threshold........:", threshold)
          if minimum_threshold is not None:
              if minimum_threshold > threshold:
                  threshold = minimum_threshold
          if abs(self.debugplot) >= 10:
              print("minimum threshold:", minimum_threshold)
              print("final threshold..:", threshold)
  
          # initial location of the peaks (integer values)
          ixpeaks0 = find_peaks_spectrum(sp0, nwinwidth=nwinwidth_initial,
                                         threshold=threshold,
                                         debugplot=self.debugplot)
  
          # peaks in the lower and upper regions
          ixpeaks1 = find_peaks_spectrum(sp1, nwinwidth=nwinwidth_initial,
                                         threshold=threshold,
                                         debugplot=self.debugplot)
          ixpeaks2 = find_peaks_spectrum(sp2, nwinwidth=nwinwidth_initial,
                                         threshold=threshold,
                                         debugplot=self.debugplot)
  
          # the peaks are valid if the are also found in the lower and
          # upper regions (with a tolerance of +1 or -1 pixel)
          ixpeaks = []
          for ixpeak in ixpeaks0:
              l1 = ixpeak in np.concatenate((ixpeaks1, ixpeaks1+1, ixpeaks1-1))
              l2 = ixpeak in np.concatenate((ixpeaks2, ixpeaks2+1, ixpeaks2-1))
              if l1 and l2:
                  ixpeaks.append(ixpeak)
          ixpeaks = np.array(ixpeaks)
          if abs(self.debugplot) >= 10:
              print("Merged initial list of peaks:\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 <= len(sp0) - 1 - npix_avoid_border
              ixpeaks = ixpeaks[lok_ini * lok_end]
  
          # select a maximum number of brightest lines in each region
          if nbrightlines is None:
              pass
          elif len(nbrightlines) == 1 and nbrightlines[0] == 0:
              pass
          else:
              if abs(self.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 = sp0[ixpeaks_region]
                          spos = peak_fluxes.argsort()
                          ixpeaks_tmp = ixpeaks_region[spos[-nlines_in_region:]]
                          ixpeaks_tmp.sort()  # in-place sort
                          if abs(self.debugplot) >= 10:
                              print('ixpeaks in region........:\n', ixpeaks_tmp)
                          ixpeaks_filtered = np.concatenate((ixpeaks_filtered,
                                                             ixpeaks_tmp))
              ixpeaks = ixpeaks_filtered
              if abs(self.debugplot) >= 10:
                  print('ixpeaks filtered.........:\n', ixpeaks)
  
          # refined location of the peaks (float values)
          fxpeaks, sxpeaks = refine_peaks_spectrum(sp0, ixpeaks,
                                                   nwinwidth=nwinwidth_refined,
                                                   method="gaussian")
  
          if abs(self.debugplot) % 10 != 0:
              x = np.arange(self.bb_nc1_orig, self.bb_nc2_orig + 1)
              title = "Slitlet#" + str(self.islitlet) + " (median spectrum)"
              ax = ximplotxy(x, sp1, show=False, title=title,
                             xlabel='pixel coordinate (from 1 to NAXIS1)',
                             ylabel='number of counts',
                             **{'marker': ' ', 'label': 'lower region'})
              ax.plot(x, sp2, label='upper region')
              ax.plot(x, sp0, label='whole region')
              # mark peak location
              ax.plot(ixpeaks + self.bb_nc1_orig,
                      sp0[ixpeaks], 'o', label="initial location")
              ax.plot(fxpeaks + self.bb_nc1_orig,
                      sp0[ixpeaks], 'o', label="refined location")
              ax.legend()
              pause_debugplot(self.debugplot, pltshow=True, tight_layout=False)
  
          # return median spectrum and refined peak location
          return sp0, fxpeaks