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])