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("")
def compute_slitlet_boundaries ( filename , grism , spfilter , list_slitlets , size_x_medfilt , size_y_savgol , times_sigma_threshold , bounddict , debugplot = 0 ) : 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 ) 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 sfilename = os . path . basename ( filename ) instrument = image_header [ 'instrume' ] if instrument != 'EMIR' : raise ValueError ( "INSTRUME keyword is not 'EMIR'!" ) csu_config = CsuConfiguration . define_from_fits ( filename ) dtu_config = DtuConfiguration . define_from_fits ( filename ) grism_in_header = image_header [ 'grism' ] if grism != grism_in_header : raise ValueError ( "GRISM keyword=" + grism_in_header + " is not the expected value=" + grism ) spfilter_in_header = image_header [ 'filter' ] if spfilter != spfilter_in_header : raise ValueError ( "FILTER keyword=" + spfilter_in_header + " is not the expected value=" + spfilter ) rotang = image_header [ 'rotang' ] 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 ) 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 ) 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 ) slitlet2d_savgol = savgol_filter ( slitlet2d_smooth , window_length = size_y_savgol , polyorder = 2 , deriv = 1 , axis = 0 ) q25 , q50 , q75 = np . percentile ( slitlet2d_savgol , q = [ 25.0 , 50.0 , 75.0 ] ) sigmag = 0.7413 * ( q75 - q25 ) 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 ) labels2d_objects_pos , no_objects_pos = ndimage . label ( slitlet2d_savgol > q50 + times_sigma_threshold * sigmag ) labels2d_objects_neg , no_objects_neg = ndimage . label ( slitlet2d_savgol < q50 - times_sigma_threshold * sigmag ) 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 ) n_der_pos = 0 i_der_pos = 0 n_der_neg = 0 i_der_neg = 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 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 ] list_boundaries = [ ] for k in range ( 2 ) : xy_tmp = np . where ( labels2d_objects == list_slices_ok [ k ] ) xmin = xy_tmp [ 1 ] . min ( ) xmax = xy_tmp [ 1 ] . max ( ) xfit = [ ] yfit = [ ] 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 ) for xdum in range ( xmineff , xmaxeff + 1 ) : iok = np . where ( xy_tmp [ 1 ] == xdum ) y_tmp = xy_tmp [ 0 ] [ iok ] + sltlim . bb_ns1_orig 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 ) boundary = SpectrumTrail ( ) 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 ) pause_debugplot ( debugplot , pltshow = True ) 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 ( "" )
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.