General beamcenter calibration

src/PyHyperScattering/BeamCentering.py

@@ -0,0 +1,88 @@
+import warnings
+import numpy as np
+import xarray as xr
+from pyFAI import azimuthalIntegrator
+from scipy.optimize import minimize
+from tqdm.auto import tqdm
+
+
+@xr.register_dataarray_accessor('beamcenter')
+class CenteringAccessor:
+
+    def __init__(self, xr_obj):
+        self._obj = xr_obj
+        self._pyhyper_type = 'reduced'
+        try:
+            self._chi_min = np.min(xr_obj.chi)
+            self._chi_max = np.max(xr_obj.chi)
+            self._chi_range = [self._chi_min, self._chi_max]
+        except AttributeError:
+            self._pyhyper_type = 'raw'
+        self.integrator = None
+        self.centering_energy = None
+
+    def create_integrator(self, energy,poni1=None,poni2=None):
+        if poni1 == None:
+            poni1 = self._obj.poni1
+        if poni2 == None:
+            poni2 = self._obj.poni2
+        return azimuthalIntegrator.AzimuthalIntegrator(
+            self._obj.dist, self._obj.poni1, self._obj.poni2,
+            self._obj.rot1, self._obj.rot2, self._obj.rot3,
+            pixel1=self._obj.pixel1, pixel2=self._obj.pixel2,
+            wavelength=1.239842e-6/energy
+        )
+
+    def optimization_func(self, x, image, obj, energy,
+                          q_min, q_max, pbar=None,
+                          chi_min=-179, chi_max=179,
+                          num_points=150):
+        #print(f'evaluating objective at {x}')
+        integrator = obj.create_integrator(energy,poni1=x[0],poni2=x[1])
+        _, image_int = integrator.integrate_radial(image, num_points,
+                                                   radial_range=(q_min, q_max),
+                                                   azimuth_range=(chi_min, chi_max),
+                                                   radial_unit='q_A^-1')
+        if pbar is not None:
+            pbar.update(1)
+        return np.var(image_int[image_int != 0])
+
+    def refine_geometry(self, energy, q_min, q_max,
+                        chi_min=-179, chi_max=179,
+                        poni1_guess=None, poni2_guess=None,
+                        bounds=None, method='Nelder-Mead',
+                        num_points=150, max_iter = 150, 
+                        verbose=False, ):
+        if not poni1_guess:
+            poni1_guess = self._obj.poni1
+            poni2_guess = self._obj.poni2
+        if (not self.integrator) | (self.centering_energy != energy):
+            self.integrator = self.create_integrator(energy=energy)
+            self.centering_energy = energy
+        if not bounds:
+            bounds = [(poni1_guess*0.9, poni1_guess*1.1),
+                      (poni2_guess*0.9, poni2_guess*1.1)]
+        options = {}
+        options['maxiter'] = max_iter
+        if verbose:
+            options['disp'] = True
+        image = self._obj.sel(energy=energy).values
+        with tqdm(total=1,desc='Optimizing Beamcenter') as pbar:
+            res = minimize(self.optimization_func, (poni1_guess, poni2_guess),
+                           args=(image, self, energy,
+                                 q_min, q_max, pbar, chi_min, chi_max,
+                                 num_points),
+                           bounds=bounds, method=method, options=options)
+
+        if res.success:
+            if (res.x == (self._obj.poni1,self._obj.poni2)).all():
+                print(f'Optimization successful, already at optimum values.  Nothing changed.')
+
+            else:
+                print(f'Optimization successful. Updating beamcenter to ({res.x[0]},{res.x[1]}), old values were ({self._obj.poni1},{self._obj.poni2})')
+                self._obj.attrs['poni1'] = res.x[0]
+                self._obj.attrs['poni2'] = res.x[1]
+
+            return res
+        else:
+            warnings.warn('Optimization was unsuccessful. Try new guesses and start again')

src/PyHyperScattering/DrawMask.py

@@ -0,0 +1,141 @@
+import warnings
+import xarray as xr
+import numpy as np
+import math
+
+try:
+    import holoviews as hv
+    import hvplot.xarray
+
+    import skimage.draw
+    import matplotlib.pyplot as plt
+    from matplotlib.colors import LogNorm,Normalize
+except (ModuleNotFoundError,ImportError):
+    warnings.warn('Could not import package for interactive integration utils.  Install holoviews and scikit-image.',stacklevel=2)
+import pandas as pd
+
+import json
+
+@xr.register_dataarray_accessor('drawmask')
+class DrawMask:
+    '''
+    Utility class for interactively drawing a mask in a Jupyter notebook.
+
+
+    Usage: 
+
+        Instantiate a DrawMask object using a PyHyper single image frame.
+
+        Call DrawMask.ui() to generate the user interface
+
+        Call DrawMask.mask to access the underlying mask, or save/load the raw mask data with .save or .load
+
+
+    '''
+
+    def __init__(self,xr_obj):
+        '''
+        Construct a DrawMask object
+
+        Args:
+            frame (xarray): a single data frame with pix_x and pix_y axes
+
+
+        '''
+        self._obj = xr_obj
+        self._pyhyper_type = 'reduced'
+        try:
+            self._chi_min = np.min(xr_obj.chi)
+            self._chi_max = np.max(xr_obj.chi)
+            self._chi_range = [self._chi_min, self._chi_max]
+        except AttributeError:
+            self._pyhyper_type = 'raw'
+
+        self.path_annotator = hv.annotate.instance()
+        self.poly = hv.Polygons([])
+
+    def ui(self,energy):
+        '''
+        Draw the DrawMask UI in a Jupyter notebook.
+
+
+        Returns: the holoviews object
+        '''
+        if energy == None:
+            if len(self._obj.shape) > 2:
+                try:
+                    img = self._obj.isel(system=0)
+                except KeyError as e:
+                    raise KeyError('This tool needs a single frame, not a stack!  .sel down to a single frame before starting!') from e
+        else:
+            img = self._obj.sel(energy=energy)
+
+        if len(img.shape) > 2:
+                warnings.warn('This tool needs a single frame, not a stack!  .sel down to a single frame before starting!',stacklevel=2)
+
+        self.frame = img
+
+        self.fig = self.frame.hvplot(cmap='terrain',clim=(5,5000),logz=True,data_aspect=1)
+
+
+        print('Usage: click the "PolyAnnotator" tool at top right.  DOUBLE CLICK to start drawing a masked object, SINGLE CLICK to add a vertex, then DOUBLE CLICK to finish.  Click/drag individual vertex to adjust.')
+        return self.path_annotator(
+                self.fig * self.poly.opts(
+                            width=self.frame.shape[0], 
+                            height=self.frame.shape[1], 
+                            responsive=False), 
+                annotations=['Label'], 
+            vertex_annotations=['Value'])
+
+
+    def save(self,fname):
+        '''
+        Save a parametric mask description as a json dump file.
+
+        Args:
+            fname (str): name of the file to save
+
+        '''
+        dflist = []
+        for i in range(len(self.path_annotator.annotated)):
+            dflist.append(self.path_annotator.annotated.iloc[i].dframe(['x','y']).to_json())
+
+        with open(fname, 'w') as outfile:
+            json.dump(dflist, outfile)
+
+    def load(self,fname):
+        '''
+        Load a parametric mask description from a json dump file.
+
+        Args:
+            fname (str): name of the file to read from
+
+        '''
+        with open(fname,'r') as f:
+            strlist = json.load(f)
+        print(strlist)
+        dflist = []
+        for item in strlist:
+            dflist.append(pd.read_json(item))
+        print(dflist)
+        self.poly = hv.Polygons(dflist)
+
+        self.path_annotator(
+                self.fig * self.poly.opts(
+                            width=self.frame.shape[1], 
+                            height=self.frame.shape[0], 
+                            responsive=False), 
+                annotations=['Label'], 
+            vertex_annotations=['Value'])
+
+
+    @property
+    def mask(self):
+        '''
+        Render the mask as a numpy boolean array.
+        '''
+        mask = np.zeros(self.frame.shape).astype(bool)
+        for i in range(len(self.path_annotator.annotated)):
+            mask |= skimage.draw.polygon2mask(self.frame.shape,self.path_annotator.annotated.iloc[i].dframe(['x','y']))
+
+        return mask

src/PyHyperScattering/IntegrationUtils.py → src/PyHyperScattering/GeoCheck.py

@@ -16,12 +16,23 @@
 
 import json
 
-class Check:
+@xr.register_dataarray_accessor('geocheck')
+class GeoCheck:
     '''
     Quick Utility to display a mask next to an image, to sanity check the orientation of e.g. an imported mask
 
     '''
-    def checkMask(integrator,img,img_min=1,img_max=10000,img_scaling='log',alpha=1):
+    def __init__(self, xr_obj):
+        self._obj = xr_obj
+        self._pyhyper_type = 'reduced'
+        try:
+            self._chi_min = np.min(xr_obj.chi)
+            self._chi_max = np.max(xr_obj.chi)
+            self._chi_range = [self._chi_min, self._chi_max]
+        except AttributeError:
+            self._pyhyper_type = 'raw'
+
+    def checkMask(self,integrator,energy=None,img_min=1,img_max=10000,img_scaling='log',alpha=1):
         '''
             draw an overlay of the mask and an image
 
@@ -32,6 +43,15 @@ def checkMask(integrator,img,img_min=1,img_max=10000,img_scaling='log',alpha=1):
                 img_max: max value to display
                 img_scaling: 'lin' or 'log'
         '''
+        if energy == None:
+            if len(self._obj.shape) > 2:
+                try:
+                    img = self._obj.isel(system=0)
+                except KeyError as e:
+                    raise KeyError('This tool needs a single frame, not a stack!  .sel down to a single frame before starting!') from e
+        else:
+            img = self._obj.sel(energy=energy)
+
         if len(img.shape) > 2:
                 warnings.warn('This tool needs a single frame, not a stack!  .sel down to a single frame before starting!',stacklevel=2)
 
@@ -43,7 +63,7 @@ def checkMask(integrator,img,img_min=1,img_max=10000,img_scaling='log',alpha=1):
         img.plot(norm=norm,ax=ax)
         ax.set_aspect(1)
         ax.imshow(integrator.mask,origin='lower',alpha=alpha)
-    def checkCenter(integrator,img,img_min=1,img_max=10000,img_scaling='log'):
+    def checkCenter(self,integrator,energy=None,img_min=1,img_max=10000,img_scaling='log'):
         '''
             draw the beamcenter on an image
 
@@ -54,6 +74,15 @@ def checkCenter(integrator,img,img_min=1,img_max=10000,img_scaling='log'):
                 img_max: max value to display
                 img_scaling: 'lin' or 'log'
         '''
+        if energy == None:
+            if len(self._obj.shape) > 2:
+                try:
+                    img = self._obj.isel(system=0)
+                except KeyError as e:
+                    raise KeyError('This tool needs a single frame, not a stack!  .sel down to a single frame before starting!') from e
+        else:
+            img = self._obj.sel(energy=energy)
+
         if len(img.shape) > 2:
                 warnings.warn('This tool needs a single frame, not a stack!  .sel down to a single frame before starting!',stacklevel=2)
 
@@ -70,7 +99,7 @@ def checkCenter(integrator,img,img_min=1,img_max=10000,img_scaling='log'):
         ax.add_patch(beamcenter)
         ax.add_patch(guide1)
         ax.add_patch(guide2)
-    def checkAll(integrator,img,img_min=1,img_max=10000,img_scaling='log',alpha=1,d_inner=50,d_outer=150):
+    def checkAll(self,integrator,energy=None,img_min=1,img_max=10000,img_scaling='log',alpha=1,d_inner=50,d_outer=150):
         '''
             draw the beamcenter and overlay mask on an image
 
@@ -81,6 +110,15 @@ def checkAll(integrator,img,img_min=1,img_max=10000,img_scaling='log',alpha=1,d_
                 img_max: max value to display
                 img_scaling: 'lin' or 'log'
         '''
+        if energy == None:
+            if len(self._obj.shape) > 2:
+                try:
+                    img = self._obj.isel(system=0)
+                except KeyError as e:
+                    raise KeyError('This tool needs a single frame, not a stack!  .sel down to a single frame before starting!') from e
+        else:
+            img = self._obj.sel(energy=energy)
+
         if len(img.shape) > 2:
                 warnings.warn('This tool needs a single frame, not a stack!  .sel down to a single frame before starting!',stacklevel=2)
 
@@ -98,6 +136,8 @@ def checkAll(integrator,img,img_min=1,img_max=10000,img_scaling='log',alpha=1,d_
         ax.add_patch(guide1)
         ax.add_patch(guide2)
         ax.imshow(integrator.mask,origin='lower',alpha=alpha)
+
+@xr.register_dataarray_accessor('drawmask')
 class DrawMask:
     '''
     Utility class for interactively drawing a mask in a Jupyter notebook.
@@ -114,7 +154,7 @@ class DrawMask:
 
     '''
 
-    def __init__(self,frame):
+    def __init__(self,xr_obj):
         '''
         Construct a DrawMask object
 
@@ -123,26 +163,42 @@ def __init__(self,frame):
 
 
         '''
-        if len(frame.shape) > 2:
-            warnings.warn('This tool needs a single frame, not a stack!  .sel down to a single frame before starting!',stacklevel=2)
+        self._obj = xr_obj
+        self._pyhyper_type = 'reduced'
+        try:
+            self._chi_min = np.min(xr_obj.chi)
+            self._chi_max = np.max(xr_obj.chi)
+            self._chi_range = [self._chi_min, self._chi_max]
+        except AttributeError:
+            self._pyhyper_type = 'raw'
 
-        self.frame=frame
 
-        self.fig = frame.hvplot(cmap='terrain',clim=(5,5000),logz=True,data_aspect=1)
-
-        self.poly = hv.Polygons([])
-        self.path_annotator = hv.annotate.instance()
 
-    def ui(self):
+    def ui(self,energy):
         '''
         Draw the DrawMask UI in a Jupyter notebook.
 
 
         Returns: the holoviews object
+        '''
+        if energy == None:
+            if len(self._obj.shape) > 2:
+                try:
+                    img = self._obj.isel(system=0)
+                except KeyError as e:
+                    raise KeyError('This tool needs a single frame, not a stack!  .sel down to a single frame before starting!') from e
+        else:
+            img = self._obj.sel(energy=energy)
+
+        if len(img.shape) > 2:
+                warnings.warn('This tool needs a single frame, not a stack!  .sel down to a single frame before starting!',stacklevel=2)
 
+        self.frame = img
+
+        self.fig = self.frame.hvplot(cmap='terrain',clim=(5,5000),logz=True,data_aspect=1)
 
-
-        '''
+        self.poly = hv.Polygons([])
+        self.path_annotator = hv.annotate.instance()
         print('Usage: click the "PolyAnnotator" tool at top right.  DOUBLE CLICK to start drawing a masked object, SINGLE CLICK to add a vertex, then DOUBLE CLICK to finish.  Click/drag individual vertex to adjust.')
         return self.path_annotator(
                 self.fig * self.poly.opts(