diff --git a/src/PyHyperScattering/BeamCentering.py b/src/PyHyperScattering/BeamCentering.py new file mode 100644 index 00000000..56c9e316 --- /dev/null +++ b/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') diff --git a/src/PyHyperScattering/DrawMask.py b/src/PyHyperScattering/DrawMask.py new file mode 100644 index 00000000..b456d68d --- /dev/null +++ b/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 diff --git a/src/PyHyperScattering/IntegrationUtils.py b/src/PyHyperScattering/GeoCheck.py similarity index 69% rename from src/PyHyperScattering/IntegrationUtils.py rename to src/PyHyperScattering/GeoCheck.py index 1ad1c537..99779fe8 100644 --- a/src/PyHyperScattering/IntegrationUtils.py +++ b/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( diff --git a/src/PyHyperScattering/util.py b/src/PyHyperScattering/util.py index 427fce77..bd1ddab1 100644 --- a/src/PyHyperScattering/util.py +++ b/src/PyHyperScattering/util.py @@ -1,7 +1,9 @@ from PyHyperScattering import Fitting from PyHyperScattering import HDR from PyHyperScattering import RSoXS -from PyHyperScattering import IntegrationUtils +from PyHyperScattering import GeoCheck +from PyHyperScattering import DrawMask #from PyHyperScattering import Nexus empty module as of 0.0.6-dev69 from PyHyperScattering import FileIO -from PyHyperScattering import PlotTools \ No newline at end of file +from PyHyperScattering import PlotTools +from PyHyperScattering import BeamCentering \ No newline at end of file