docs: def compute_slitlet_boundaries(
          filename, grism, spfilter, list_slitlets,
          size_x_medfilt, size_y_savgol,
          times_sigma_threshold,
          bounddict, debugplot=0):
      """Compute slitlet boundaries using continuum lamp images.
  
      Parameters
      ----------
      filename : string
          Input continumm lamp image.
      grism : string
          Grism name. It must be one in EMIR_VALID_GRISMS.
      spfilter : string
          Filter name. It must be one in EMIR_VALID_FILTERS.
      list_slitlets : list of integers
          Number of slitlets to be updated.
      size_x_medfilt : int
          Window in the X (spectral) direction, in pixels, to apply
          the 1d median filter in order to remove bad pixels.
      size_y_savgol : int
          Window in the Y (spatial) direction to be used when using the
          1d Savitzky-Golay filter.
      times_sigma_threshold : float
          Times sigma to detect peaks in derivatives.
      bounddict : dictionary of dictionaries
          Structure to store the boundaries.
      debugplot : int
          Determines whether intermediate computations and/or plots are
          displayed.
  
      """
  
      # read 2D image
      hdulist = fits.open(filename)
      image_header = hdulist[0].header
      image2d = hdulist[0].data
      naxis2, naxis1 = image2d.shape
      hdulist.close()
      if debugplot >= 10:
          print('>>> NAXIS1:', naxis1)
          print('>>> NAXIS2:', naxis2)
  
  
      # ToDo: replace this by application of cosmetic defect mask!
      for j in range(1024):
          image2d[1024, j] = (image2d[1023, j] + image2d[1025, j]) / 2
          image2d[1023, j + 1024] = (image2d[1022, j + 1024] +
                                     image2d[1024, j + 1024]) / 2
          
      # remove path from filename
      sfilename = os.path.basename(filename)
  
      # check that the FITS file has been obtained with EMIR
      instrument = image_header['instrume']
      if instrument != 'EMIR':
          raise ValueError("INSTRUME keyword is not 'EMIR'!")
  
      # read CSU configuration from FITS header
      csu_config = CsuConfiguration.define_from_fits(filename)
  
      # read DTU configuration from FITS header
      dtu_config = DtuConfiguration.define_from_fits(filename)
  
      # read grism
      grism_in_header = image_header['grism']
      if grism != grism_in_header:
          raise ValueError("GRISM keyword=" + grism_in_header +
                           " is not the expected value=" + grism)
  
      # read filter
      spfilter_in_header = image_header['filter']
      if spfilter != spfilter_in_header:
          raise ValueError("FILTER keyword=" + spfilter_in_header +
                           " is not the expected value=" + spfilter)
  
      # read rotator position angle
      rotang = image_header['rotang']
  
      # read date-obs
      date_obs = image_header['date-obs']
  
      for islitlet in list_slitlets:
          if debugplot < 10:
              sys.stdout.write('.')
              sys.stdout.flush()
  
          sltlim = SlitletLimits(grism, spfilter, islitlet)
          # extract slitlet2d
          slitlet2d = extract_slitlet2d(image2d, sltlim)
          if debugplot % 10 != 0:
              ximshow(slitlet2d,
                      title=sfilename + " [original]"
                            "\nslitlet=" + str(islitlet) +
                            ", grism=" + grism +
                            ", filter=" + spfilter +
                            ", rotang=" + str(round(rotang, 2)),
                      first_pixel=(sltlim.bb_nc1_orig, sltlim.bb_ns1_orig),
                      debugplot=debugplot)
  
          # apply 1d median filtering (along the spectral direction)
          # to remove bad pixels
          size_x = size_x_medfilt
          size_y = 1
          slitlet2d_smooth = ndimage.filters.median_filter(
              slitlet2d, size=(size_y, size_x))
  
          if debugplot % 10 != 0:
              ximshow(slitlet2d_smooth,
                      title=sfilename + " [smoothed]"
                            "\nslitlet=" + str(islitlet) +
                            ", grism=" + grism +
                            ", filter=" + spfilter +
                            ", rotang=" + str(round(rotang, 2)),
                      first_pixel=(sltlim.bb_nc1_orig, sltlim.bb_ns1_orig),
                      debugplot=debugplot)
  
          # apply 1d Savitzky-Golay filter (along the spatial direction)
          # to compute first derivative
          slitlet2d_savgol = savgol_filter(
              slitlet2d_smooth, window_length=size_y_savgol, polyorder=2,
              deriv=1, axis=0)
  
          # compute basic statistics
          q25, q50, q75 = np.percentile(slitlet2d_savgol, q=[25.0, 50.0, 75.0])
          sigmag = 0.7413 * (q75 - q25)  # robust standard deviation
          if debugplot >= 10:
              print("q50, sigmag:", q50, sigmag)
  
          if debugplot % 10 != 0:
              ximshow(slitlet2d_savgol,
                      title=sfilename + " [S.-G.filt.]"
                            "\nslitlet=" + str(islitlet) +
                            ", grism=" + grism +
                            ", filter=" + spfilter +
                            ", rotang=" + str(round(rotang, 2)),
                      first_pixel=(sltlim.bb_nc1_orig, sltlim.bb_ns1_orig),
                      z1z2=(q50-times_sigma_threshold*sigmag,
                            q50+times_sigma_threshold*sigmag),
                      debugplot=debugplot)
  
          # identify objects in slitlet2d_savgol: pixels with positive
          # derivatives are identify independently from pixels with
          # negative derivaties; then the two set of potential features
          # are merged; this approach avoids some problems when, in
          # nearby regions, there are pixels with positive and negative
          # derivatives (in those circumstances a single search as
          # np.logical_or(
          #     slitlet2d_savgol < q50 - times_sigma_threshold * sigmag,
          #     slitlet2d_savgol > q50 + times_sigma_threshold * sigmag)
          # led to erroneous detections!)
          #
          # search for positive derivatives
          labels2d_objects_pos, no_objects_pos = ndimage.label(
              slitlet2d_savgol > q50 + times_sigma_threshold * sigmag)
          # search for negative derivatives
          labels2d_objects_neg, no_objects_neg = ndimage.label(
              slitlet2d_savgol < q50 - times_sigma_threshold * sigmag)
          # merge both sets
          non_zero_neg = np.where(labels2d_objects_neg > 0)
          labels2d_objects = np.copy(labels2d_objects_pos)
          labels2d_objects[non_zero_neg] += \
              labels2d_objects_neg[non_zero_neg] + no_objects_pos
          no_objects = no_objects_pos + no_objects_neg
  
          if debugplot >= 10:
              print("Number of objects with positive derivative:",
                    no_objects_pos)
              print("Number of objects with negative derivative:",
                    no_objects_neg)
              print("Total number of objects initially found...:", no_objects)
  
          if debugplot % 10 != 0:
              ximshow(labels2d_objects,
                      z1z2=(0, no_objects),
                      title=sfilename + " [objects]"
                                       "\nslitlet=" + str(islitlet) +
                            ", grism=" + grism +
                            ", filter=" + spfilter +
                            ", rotang=" + str(round(rotang, 2)),
                      first_pixel=(sltlim.bb_nc1_orig, sltlim.bb_ns1_orig),
                      cbar_label="Object number",
                      debugplot=debugplot)
  
          # select boundaries as the largest objects found with
          # positive and negative derivatives
          n_der_pos = 0  # number of pixels covered by the object with deriv > 0
          i_der_pos = 0  # id of the object with deriv > 0
          n_der_neg = 0  # number of pixels covered by the object with deriv < 0
          i_der_neg = 0  # id of the object with deriv < 0
          for i in range(1, no_objects+1):
              xy_tmp = np.where(labels2d_objects == i)
              n_pix = len(xy_tmp[0])
              if i <= no_objects_pos:
                  if n_pix > n_der_pos:
                      i_der_pos = i
                      n_der_pos = n_pix
              else:
                  if n_pix > n_der_neg:
                      i_der_neg = i
                      n_der_neg = n_pix
  
          # determine which boundary is lower and which is upper
          y_center_mass_der_pos = ndimage.center_of_mass(
              slitlet2d_savgol, labels2d_objects, [i_der_pos])[0][0]
          y_center_mass_der_neg = ndimage.center_of_mass(
              slitlet2d_savgol, labels2d_objects, [i_der_neg])[0][0]
          if y_center_mass_der_pos < y_center_mass_der_neg:
              i_lower = i_der_pos
              i_upper = i_der_neg
              if debugplot >= 10:
                  print("-> lower boundary has positive derivatives")
          else:
              i_lower = i_der_neg
              i_upper = i_der_pos
              if debugplot >= 10:
                  print("-> lower boundary has negative derivatives")
          list_slices_ok = [i_lower, i_upper]
  
          # adjust individual boundaries passing the selection:
          # - select points in the image belonging to a given boundary
          # - compute weighted mean of the pixels of the boundary, column
          #   by column (this reduces dramatically the number of points
          #   to be fitted to determine the boundary)
          list_boundaries = []
          for k in range(2):  # k=0 lower boundary, k=1 upper boundary
              # select points to be fitted for a particular boundary
              # (note: be careful with array indices and pixel
              # coordinates)
              xy_tmp = np.where(labels2d_objects == list_slices_ok[k])
              xmin = xy_tmp[1].min()  # array indices (integers)
              xmax = xy_tmp[1].max()  # array indices (integers)
              xfit = []
              yfit = []
              # fix range for fit
              if k == 0:
                  xmineff = max(sltlim.xmin_lower_boundary_fit, xmin)
                  xmaxeff = min(sltlim.xmax_lower_boundary_fit, xmax)
              else:
                  xmineff = max(sltlim.xmin_upper_boundary_fit, xmin)
                  xmaxeff = min(sltlim.xmax_upper_boundary_fit, xmax)
              # loop in columns of the image belonging to the boundary
              for xdum in range(xmineff, xmaxeff + 1):  # array indices (integer)
                  iok = np.where(xy_tmp[1] == xdum)
                  y_tmp = xy_tmp[0][iok] + sltlim.bb_ns1_orig  # image pixel
                  weight = slitlet2d_savgol[xy_tmp[0][iok], xy_tmp[1][iok]]
                  y_wmean = sum(y_tmp * weight) / sum(weight)
                  xfit.append(xdum + sltlim.bb_nc1_orig)
                  yfit.append(y_wmean)
              xfit = np.array(xfit)
              yfit = np.array(yfit)
              # declare new SpectrumTrail instance
              boundary = SpectrumTrail()
              # define new boundary
              boundary.fit(x=xfit, y=yfit, deg=sltlim.deg_boundary,
                           times_sigma_reject=10,
                           title="slit:" + str(sltlim.islitlet) +
                                 ", deg=" + str(sltlim.deg_boundary),
                           debugplot=0)
              list_boundaries.append(boundary)
  
          if debugplot % 10 != 0:
              for tmp_img, tmp_label in zip(
                  [slitlet2d_savgol, slitlet2d],
                  [' [S.-G.filt.]', ' [original]']
              ):
                  ax = ximshow(tmp_img,
                               title=sfilename + tmp_label +
                                     "\nslitlet=" + str(islitlet) +
                                     ", grism=" + grism +
                                     ", filter=" + spfilter +
                                     ", rotang=" + str(round(rotang, 2)),
                               first_pixel=(sltlim.bb_nc1_orig,
                                            sltlim.bb_ns1_orig),
                               show=False,
                               debugplot=debugplot)
                  for k in range(2):
                      xpol, ypol = list_boundaries[k].linspace_pix(
                          start=1, stop=EMIR_NAXIS1)
                      ax.plot(xpol, ypol, 'b--', linewidth=1)
                  for k in range(2):
                      xpol, ypol = list_boundaries[k].linspace_pix()
                      ax.plot(xpol, ypol, 'g--', linewidth=4)
                  # show plot
                  pause_debugplot(debugplot, pltshow=True)
  
          # update bounddict
          tmp_dict = {
              'boundary_coef_lower':
                  list_boundaries[0].poly_funct.coef.tolist(),
              'boundary_xmin_lower': list_boundaries[0].xlower_line,
              'boundary_xmax_lower': list_boundaries[0].xupper_line,
              'boundary_coef_upper':
                  list_boundaries[1].poly_funct.coef.tolist(),
              'boundary_xmin_upper': list_boundaries[1].xlower_line,
              'boundary_xmax_upper': list_boundaries[1].xupper_line,
              'csu_bar_left': csu_config.csu_bar_left(islitlet),
              'csu_bar_right': csu_config.csu_bar_right(islitlet),
              'csu_bar_slit_center': csu_config.csu_bar_slit_center(islitlet),
              'csu_bar_slit_width': csu_config.csu_bar_slit_width(islitlet),
              'rotang': rotang,
              'xdtu': dtu_config.xdtu,
              'ydtu': dtu_config.ydtu,
              'zdtu': dtu_config.zdtu,
              'xdtu_0': dtu_config.xdtu_0,
              'ydtu_0': dtu_config.ydtu_0,
              'zdtu_0': dtu_config.zdtu_0,
              'zzz_info1': os.getlogin() + '@' + socket.gethostname(),
              'zzz_info2': datetime.now().isoformat()
          }
          slitlet_label = "slitlet" + str(islitlet).zfill(2)
          if slitlet_label not in bounddict['contents']:
              bounddict['contents'][slitlet_label] = {}
          bounddict['contents'][slitlet_label][date_obs] = tmp_dict
  
      if debugplot < 10:
          print("")