Updates from the Beamline

patches/chx_outlier_detection.py

@@ -0,0 +1,98 @@
+def is_outlier(points,thresh=3.5,verbose=False):
+    """MAD test
+    """     
+    points.tolist()
+    if len(points) ==1:
+        points=points[:,None]
+        if verbose:
+            print('input to is_outlier is a single point...')
+    median = np.median(points)*np.ones(np.shape(points))#, axis=0)
+
+    diff = (points-median)**2
+    diff=np.sqrt(diff)
+    med_abs_deviation= np.median(diff)
+    modified_z_score = .6745*diff/med_abs_deviation
+    return modified_z_score > thresh
+
+def outlier_mask(avg_img,mask,roi_mask,outlier_threshold = 7.5,maximum_outlier_fraction = .1,verbose=False,plot=False):
+    """
+    outlier_mask(avg_img,mask,roi_mask,outlier_threshold = 7.5,maximum_outlier_fraction = .1,verbose=False,plot=False)
+    avg_img: average image data (2D)
+    mask: 2D array, same size as avg_img with pixels that are already masked
+    roi_mask: 2D array, same size as avg_img, ROI labels 'encoded' as mask values (i.e. all pixels belonging to ROI 5 have the value 5)
+    outlier_threshold: threshold for MAD test
+    maximum_outlier_fraction: maximum fraction of pixels in an ROI that can be classifed as outliers. If the detected fraction is higher, no outliers will be masked for that ROI.
+    verbose: 'True' enables message output
+    plot: 'True' enables visualization of outliers
+    returns: mask (dtype=float): 0 for pixels that have been classified as outliers, 1 else
+    dependency: is_outlier()
+
+    function does outlier detection for each ROI separately based on pixel intensity in avg_img*mask and ROI specified by roi_mask, using the median-absolute-deviation (MAD) method
+
+    by LW 06/21/2023
+    """
+    hhmask = np.ones(np.shape(roi_mask))
+    pc=1
+
+    for rn in np.arange(1,np.max(roi_mask)+1,1):
+        rm=np.zeros(np.shape(roi_mask));rm=rm-1;rm[np.where( roi_mask == rn)]=1
+        pixel = roi.roi_pixel_values(avg_img*rm, roi_mask, [rn] )
+        out_l = is_outlier((avg_img*mask*rm)[rm>-1], thresh=outlier_threshold)
+        if np.nanmax(out_l)>0: # Did detect at least one outlier
+            ave_roi_int = np.nanmean((pixel[0][0])[out_l<1])
+            if verbose: print('ROI #%s\naverage ROI intensity: %s'%(rn,ave_roi_int))
+            try:
+                upper_outlier_threshold = np.nanmin((out_l*pixel[0][0])[out_l*pixel[0][0]>ave_roi_int])
+                if verbose: print('upper outlier threshold: %s'%upper_outlier_threshold)
+            except:
+                upper_outlier_threshold = False
+                if verbose: print('no upper outlier threshold found')
+            ind1 = (out_l*pixel[0][0])>0; ind2 =  (out_l*pixel[0][0])< ave_roi_int
+            try:
+                lower_outlier_threshold = np.nanmax((out_l*pixel[0][0])[ind1*ind2])
+            except:
+                lower_outlier_threshold = False
+                if verbose: print('no lower outlier threshold found')
+        else:
+            if verbose:  print('ROI #%s: no outliers detected'%rn)
+
+        ### MAKE SURE we don't REMOVE more than x percent of the pixels in the roi 
+        outlier_fraction =  np.sum(out_l)/len(pixel[0][0])
+        if verbose: print('fraction of pixel values detected as outliers: %s'%np.round(outlier_fraction,2))
+        if outlier_fraction > maximum_outlier_fraction:
+            if verbose: print('fraction of pixel values detected as outliers > than maximum fraction %s allowed -> NOT masking outliers...check threshold for MAD and maximum fraction of outliers allowed'%maximum_outlier_fraction)
+            upper_outlier_threshold = False; lower_outlier_threshold = False
+
+        if upper_outlier_threshold:
+            hhmask[avg_img*rm > upper_outlier_threshold] = 0
+        if lower_outlier_threshold:
+            hhmask[avg_img*rm < lower_outlier_threshold] = 0
+
+        if plot:
+            if pc == 1: fig,ax = plt.subplots(1,5,figsize=(24,4))
+            plt.subplot(1,5,pc);pc+=1;
+            if pc>5: pc=1
+            pixel = roi.roi_pixel_values(avg_img*rm*mask, roi_mask, [rn] )
+            plt.plot( pixel[0][0] ,'bo',markersize=1.5 )
+            if upper_outlier_threshold or lower_outlier_threshold:
+                x=np.arange(len(out_l))
+                plt.plot([x[0],x[-1]],[ave_roi_int,ave_roi_int],'g--',label='ROI average: %s'%np.round(ave_roi_int,4))
+            if upper_outlier_threshold:
+                ind=(out_l*pixel[0][0])> upper_outlier_threshold
+                plt.plot(x[ind],(out_l*pixel[0][0])[ind],'r+')
+                plt.plot([x[0],x[-1]],[upper_outlier_threshold,upper_outlier_threshold],'r--',label='upper thresh.: %s'%np.round(upper_outlier_threshold,4))
+            if lower_outlier_threshold:
+                ind=(out_l*pixel[0][0])< lower_outlier_threshold
+                plt.plot(x[ind],(out_l*pixel[0][0])[ind],'r+')
+                plt.plot([x[0],x[-1]],[lower_outlier_threshold,lower_outlier_threshold],'r--',label='lower thresh.: %s'%np.round(upper_outlier_threshold,4))
+            plt.ylabel('Intensity') ;plt.xlabel('pixel');plt.title('ROI #: %s'%rn);plt.legend(loc='best',fontsize=8)
+
+    if plot:
+        fig,ax = plt.subplots()
+        plt.imshow(hhmask)
+        hot_dark=np.nonzero(hhmask<1)
+        cmap = plt.cm.get_cmap('viridis')
+        plt.plot(hot_dark[1],hot_dark[0],'+',color=cmap(0))
+        plt.xlabel('pixel');plt.ylabel('pixel');plt.title('masked pixels with outlier threshold: %s'%outlier_threshold)
+
+    return hhmask

patches/fix_get_sid_filenames.py

@@ -0,0 +1,29 @@
+def get_sid_filenames(hdr,verbose=False):
+    import glob
+    from time import strftime, localtime
+    start_doc = hdr.start
+    stop_doc = hdr.stop
+    success = False
+
+    ret = (start_doc["scan_id"], start_doc["uid"], glob.glob(f"{start_doc['data path']}*_{start_doc['sequence id']}_master.h5")) # looking for (eiger) datafile at the path specified in metadata
+    if len(ret[2])==0:
+        if verbose: print('could not find detector filename from "data_path" in metadata: %s'%start_doc['data path'])
+    else:
+         if verbose: print('Found detector filename from "data_path" in metadata!');success=True
+
+    if not success: # looking at path in metadata, but taking the date from the run start document
+        data_path=start_doc['data path'][:-11]+strftime("%Y/%m/%d/",localtime(start_doc['time']))
+        ret = (start_doc["scan_id"], start_doc["uid"], glob.glob(f"{data_path}*_{start_doc['sequence id']}_master.h5"))
+        if len(ret[2])==0:
+            if verbose: print('could not find detector filename in %s'%data_path)
+        else:
+             if verbose: print('Found detector filename in %s'%data_path);success=True
+
+    if not success: # looking at path in metadata, but taking the date from the run stop document (in case the date rolled over between creating the start doc and staging the detector)
+        data_path=start_doc['data path'][:-11]+strftime("%Y/%m/%d/",localtime(stop_doc['time']))
+        ret = (start_doc["scan_id"], start_doc["uid"], glob.glob(f"{data_path}*_{start_doc['sequence id']}_master.h5"))
+        if len(ret[2])==0:
+            if verbose: print('Sorry, could not find detector filename....')
+        else:
+             if verbose: print('Found detector filename in %s'%data_path);success=True
+    return ret 

patches/polygonmask_fix.py

@@ -0,0 +1,99 @@
+# fix for newer environments on jupyter hub: polygon( y,x, shape = image.shape) -> need to be specific about shape
+# Not needed for older, local environments on srv1,2,3. Fix above will probably also work there, BUT there would be a problem with disk <-> circle that has been renamed in skimage at some point
+# -> older, local environments on srv1,2,3 work without this fix, only apply when running on jupyter hub
+def create_multi_rotated_rectangle_mask(  image, center=None, length=100, width=50, angles=[0] ):
+    ''' Developed at July 10, 2017 by Y.G.@CHX, NSLS2
+     Create multi rectangle-shaped mask by rotating a rectangle with a list of angles
+     The original rectangle is defined by four corners, i.e.,
+         [   (center[1] - width//2, center[0]),
+             (center[1] + width//2, center[0]),
+             (center[1] + width//2, center[0] + length),
+             (center[1] - width//2, center[0] + length)
+          ]
+
+     Parameters:
+         image: 2D numpy array, to give mask shape
+         center: integer list, if None, will be the center of the image
+         length: integer, the length of the non-ratoted rectangle
+         width: integer,  the width of the non-ratoted rectangle
+         angles: integer list, a list of rotated angles
+
+    Return:
+        mask: 2D bool-type numpy array
+    '''
+
+    from skimage.draw import  polygon
+    from skimage.transform import rotate
+    cx,cy = center
+    imy, imx = image.shape
+    mask = np.zeros( image.shape, dtype = bool)
+    wy =  length
+    wx =   width
+    x = np.array( [ max(0, cx - wx//2), min(imx, cx+wx//2), min(imx, cx+wx//2), max(0,cx-wx//2 ) ])
+    y = np.array( [ cy, cy, min( imy, cy + wy) , min(imy, cy + wy) ])
+    rr, cc = polygon( y,x, shape = image.shape)
+    mask[rr,cc] =1
+    mask_rot=  np.zeros( image.shape, dtype = bool)
+    for angle in angles:
+        mask_rot +=  np.array( rotate( mask, angle, center= center ), dtype=bool) #, preserve_range=True)
+    return  ~mask_rot
+
+def create_cross_mask(  image, center, wy_left=4, wy_right=4, wx_up=4, wx_down=4,
+                     center_disk = True, center_radius=10
+                     ):
+    '''
+    Give image and the beam center to create a cross-shaped mask
+    wy_left: the width of left h-line
+    wy_right: the width of rigth h-line
+    wx_up: the width of up v-line
+    wx_down: the width of down v-line
+    center_disk: if True, create a disk with center and center_radius
+
+    Return:
+    the cross mask
+    '''
+    from skimage.draw import line_aa, line, polygon, disk
+
+    imy, imx = image.shape
+    cx,cy = center
+    bst_mask = np.zeros_like( image , dtype = bool)
+    ###
+    #for right part
+    wy = wy_right
+    x = np.array( [ cx, imx, imx, cx  ])
+    y = np.array( [ cy-wy, cy-wy, cy + wy, cy + wy])
+    rr, cc = polygon( y,x,shape=image.shape)
+    bst_mask[rr,cc] =1
+
+    ###
+    #for left part
+    wy = wy_left
+    x = np.array( [0,  cx, cx,0  ])
+    y = np.array( [ cy-wy, cy-wy, cy + wy, cy + wy])
+    rr, cc = polygon( y,x,shape=image.shape)
+    bst_mask[rr,cc] =1
+
+    ###
+    #for up part
+    wx = wx_up
+    x = np.array( [ cx-wx, cx + wx, cx+wx, cx-wx  ])
+    y = np.array( [ cy, cy, imy, imy])
+    rr, cc = polygon( y,x,shape=image.shape)
+    bst_mask[rr,cc] =1
+
+    ###
+    #for low part
+    wx = wx_down
+    x = np.array( [ cx-wx, cx + wx, cx+wx, cx-wx  ])
+    y = np.array( [ 0,0, cy, cy])
+    rr, cc = polygon( y,x,shape=image.shape)
+    bst_mask[rr,cc] =1
+
+    if center_radius!=0:
+        rr, cc = disk((cy, cx), center_radius, shape = bst_mask.shape)
+        bst_mask[rr,cc] =1
+
+
+    full_mask= ~bst_mask
+
+    return full_mask

patches/roi_nr_2019_3_0_1.py

@@ -0,0 +1,49 @@
+def get_roi_nr(qdict,q,phi,q_nr=True,phi_nr=False,q_thresh=0, p_thresh=0, silent=True, qprecision=5):
+    """
+    function to return roi number from qval_dict, corresponding  Q and phi, lists (sets) of all available Qs and phis
+    [roi_nr,Q,phi,Q_list,phi_list]=get_roi_nr(..)
+    calling sequence: get_roi_nr(qdict,q,phi,q_nr=True,phi_nr=False, verbose=True)
+    qdict: qval_dict from analysis pipeline/hdf5 result file
+    q: q of interest, can be either value (q_nr=False) or q-number (q_nr=True)
+    q_thresh: threshold for comparing Q-values, set to 0 for exact comparison
+    phi: phi of interest, can be either value (phi_nr=False) or q-number (phi_nr=True)
+    p_thresh: threshold for comparing phi values, set to 0 for exact comparison
+    silent=True/False: Don't/Do print lists of available qs and phis, q and phi of interest
+    by LW 10/21/2017
+    update by LW 08/22/2018: introduced thresholds for comparison of Q and phi values (before: exact match required)
+    update 2019/09/28 add qprecision to get unique Q
+    update 2020/3/12 explicitly order input dictionary to fix problem with environment after 2019-3.0.1
+    """
+    import collections
+    from collections import OrderedDict
+    qdict = collections.OrderedDict(sorted(qdict.items()))
+    qs=[]
+    phis=[]
+    for i in qdict.keys():
+        qs.append(qdict[i][0])
+        phis.append(qdict[i][1])  
+    qslist=list(OrderedDict.fromkeys(qs))
+    qslist = np.unique( np.round(qslist, qprecision ) )
+    phislist=list(OrderedDict.fromkeys(phis))
+    qslist=list(np.sort(qslist))
+    phislist=list(np.sort(phislist))
+    if q_nr:
+        qinterest=qslist[q]
+        qindices = [i for i,x in enumerate(qs) if np.abs(x-qinterest) < q_thresh]
+    else:
+        qinterest=q
+        qindices = [i for i,x in enumerate(qs) if np.abs(x-qinterest) < q_thresh] # new
+    if phi_nr:
+        phiinterest=phislist[phi]
+        phiindices = [i for i,x in enumerate(phis) if x == phiinterest]
+    else:
+        phiinterest=phi
+        phiindices = [i for i,x in enumerate(phis) if np.abs(x-phiinterest) < p_thresh] # new
+    ret_list=[list(set(qindices).intersection(phiindices))[0],qinterest,phiinterest,qslist,phislist] #-> this is the original
+    if silent == False:
+        print('list of available Qs:')
+        print(qslist)
+        print('list of available phis:')
+        print(phislist)
+        print('Roi number for Q= '+str(ret_list[1])+' and phi= '+str(ret_list[2])+': '+str(ret_list[0]))
+    return ret_list

pyCHX/XPCS_GiSAXS.py

@@ -4,6 +4,8 @@
 This module is for the GiSAXS XPCS analysis
 """
 
+
+
 from skbeam.core.accumulators.binned_statistic import BinnedStatistic1D, BinnedStatistic2D
 
 from pyCHX.chx_compress import (
@@ -771,7 +773,8 @@ def show_label_array_on_image(
     """
     ax.set_aspect("equal")
     if log_img:
-        im = ax.imshow(image, cmap=imshow_cmap, interpolation="none", norm=LogNorm(norm), **kwargs)  # norm=norm,
+        norm=LogNorm(vmin=vmin, vmax=vmax)
+        im = ax.imshow(image, cmap=imshow_cmap, interpolation="none", norm=norm)  # norm=norm,
     else:
         im = ax.imshow(image, cmap=imshow_cmap, interpolation="none", norm=norm, **kwargs)  # norm=norm,