docs: def median_slitlets_rectified(
          input_image,
          mode=0,
          minimum_slitlet_width_mm=EMIR_MINIMUM_SLITLET_WIDTH_MM,
          maximum_slitlet_width_mm=EMIR_MAXIMUM_SLITLET_WIDTH_MM,
          debugplot=0
      ):
      """Compute median spectrum for each slitlet.
  
      Parameters
      ----------
      input_image : HDUList object
          Input 2D image.
      mode : int
          Indicate desired result:
          0 : image with the same size as the input image, with the
              median spectrum of each slitlet spanning all the spectra
              of the corresponding slitlet
          1 : image with 55 spectra, containing the median spectra of
              each slitlet
          2 : single collapsed median spectrum, using exclusively the
              useful slitlets from the input image
      minimum_slitlet_width_mm : float
          Minimum slitlet width (mm) for a valid slitlet.
      maximum_slitlet_width_mm : float
          Maximum slitlet width (mm) for a valid slitlet.
      debugplot : int
          Determines whether intermediate computations and/or plots
          are displayed. The valid codes are defined in
          numina.array.display.pause_debugplot.
  
      Returns
      -------
      image_median : HDUList object
          Output image.
  
      """
  
      image_header = input_image[0].header
      image2d = input_image[0].data
  
      # check image dimensions
      naxis2_expected = EMIR_NBARS * EMIR_NPIXPERSLIT_RECTIFIED
  
      naxis2, naxis1 = image2d.shape
      if naxis2 != naxis2_expected:
          raise ValueError("NAXIS2={0} should be {1}".format(
              naxis2, naxis2_expected
          ))
  
      # check that the FITS file has been obtained with EMIR
      instrument = image_header['instrume']
      if instrument != 'EMIR':
          raise ValueError("INSTRUME keyword is not 'EMIR'!")
  
      # initialize output image
      if mode == 0:
          image2d_median = np.zeros((naxis2, naxis1))
      else:
          image2d_median = np.zeros((EMIR_NBARS, naxis1))
  
      # main loop
      for i in range(EMIR_NBARS):
          ns1 = i * EMIR_NPIXPERSLIT_RECTIFIED + 1
          ns2 = ns1 + EMIR_NPIXPERSLIT_RECTIFIED - 1
          sp_median = np.median(image2d[(ns1-1):ns2, :], axis=0)
  
          if mode == 0:
              image2d_median[(ns1-1):ns2, :] = np.tile(
                  sp_median, (EMIR_NPIXPERSLIT_RECTIFIED, 1)
              )
          else:
              image2d_median[i] = np.copy(sp_median)
  
      if mode == 2:
          # get CSU configuration from FITS header
          csu_config = CsuConfiguration.define_from_header(image_header)
  
          # define wavelength calibration parameters
          crpix1 = image_header['crpix1']
          crval1 = image_header['crval1']
          cdelt1 = image_header['cdelt1']
  
          # segregate slitlets
          list_useful_slitlets = csu_config.widths_in_range_mm(
              minwidth=minimum_slitlet_width_mm,
              maxwidth=maximum_slitlet_width_mm
          )
          list_not_useful_slitlets = [i for i in list(range(1, EMIR_NBARS + 1))
                                      if i not in list_useful_slitlets]
          if abs(debugplot) != 0:
              print('>>> list_useful_slitlets....:', list_useful_slitlets)
              print('>>> list_not_useful_slitlets:', list_not_useful_slitlets)
  
          # define mask from array data
          mask2d, borders = define_mask_borders(image2d_median, sought_value=0)
          if abs(debugplot) % 10 != 0:
              ximshow(mask2d.astype(int), z1z2=(-.2, 1.2), crpix1=crpix1,
                      crval1=crval1, cdelt1=cdelt1, debugplot=debugplot)
  
          # update mask with unused slitlets
          for islitlet in list_not_useful_slitlets:
              mask2d[islitlet - 1, :] = np.array([True] * naxis1)
          if abs(debugplot) % 10 != 0:
              ximshow(mask2d.astype(int), z1z2=(-.2, 1.2), crpix1=crpix1,
                      crval1=crval1, cdelt1=cdelt1, debugplot=debugplot)
  
          # useful image pixels
          image2d_masked = image2d_median * (1 - mask2d.astype(int))
          if abs(debugplot) % 10 != 0:
              ximshow(image2d_masked, crpix1=crpix1, crval1=crval1,
                      cdelt1=cdelt1, debugplot=debugplot)
  
          # masked image
          image2d_masked = np.ma.masked_array(image2d_median, mask=mask2d)
          # median spectrum
          image1d_median = np.ma.median(image2d_masked, axis=0).data
  
          image_median = fits.PrimaryHDU(data=image1d_median,
                                         header=image_header)
  
      else:
          image_median = fits.PrimaryHDU(data=image2d_median,
                                         header=image_header)
  
      return fits.HDUList([image_median])