From ecc19f40b92d162344b36d27dac70109ae16eab4 Mon Sep 17 00:00:00 2001 From: Peter Dudenas Date: Thu, 2 Feb 2023 11:44:41 -0800 Subject: [PATCH 1/5] First attempt at beamcentering util --- src/PyHyperScattering/BeamCentering.py | 70 ++++++++++++++++++++++++++ src/PyHyperScattering/util.py | 3 +- 2 files changed, 72 insertions(+), 1 deletion(-) create mode 100644 src/PyHyperScattering/BeamCentering.py diff --git a/src/PyHyperScattering/BeamCentering.py b/src/PyHyperScattering/BeamCentering.py new file mode 100644 index 00000000..c6c04c7b --- /dev/null +++ b/src/PyHyperScattering/BeamCentering.py @@ -0,0 +1,70 @@ +import warnings +import numpy as np +import xarray as xr +from pyFAI import azimuthalIntegrator +from scipy.optimize import minimize + + +@xr.register_dataarray_accessor('util') +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): + 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, integrator, + q_min, q_max, + chi_min=-179, chi_max=179, + num_points=150): + poni1, poni2 = x + integrator.poni1 = poni1 + integrator.poni2 = poni2 + _, 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') + 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): + 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)] + image = self._obj.unstack('system').sel(energy=energy).values + res = minimize(self.optimization_func, (poni1_guess, poni2_guess), + args=(image, self.integrator, + q_min, q_max, chi_min, chi_max, + num_points), + bounds=bounds, method=method) + + if res.success: + print('Optimization successful. Updating beamcenter') + self._obj.attrs['poni1'] = res.x[0] + self._obj.attrs['poni2'] = res.x[1] + else: + warnings.warn('Optimization was unsuccessful. Try new guesses and start again') diff --git a/src/PyHyperScattering/util.py b/src/PyHyperScattering/util.py index 427fce77..f792066c 100644 --- a/src/PyHyperScattering/util.py +++ b/src/PyHyperScattering/util.py @@ -4,4 +4,5 @@ from PyHyperScattering import IntegrationUtils #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 From e66881b27f7b70e17d8146e48e5204be69a4afaf Mon Sep 17 00:00:00 2001 From: Peter Beaucage Date: Fri, 7 Apr 2023 11:20:28 -0400 Subject: [PATCH 2/5] Performance optimization in beamcentering --- src/PyHyperScattering/BeamCentering.py | 42 +++++++++++++++++--------- 1 file changed, 27 insertions(+), 15 deletions(-) diff --git a/src/PyHyperScattering/BeamCentering.py b/src/PyHyperScattering/BeamCentering.py index c6c04c7b..3eefc27a 100644 --- a/src/PyHyperScattering/BeamCentering.py +++ b/src/PyHyperScattering/BeamCentering.py @@ -3,6 +3,7 @@ import xarray as xr from pyFAI import azimuthalIntegrator from scipy.optimize import minimize +from tqdm.auto import tqdm @xr.register_dataarray_accessor('util') @@ -20,7 +21,11 @@ def __init__(self, xr_obj): self.integrator = None self.centering_energy = None - def create_integrator(self, energy): + 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, @@ -28,24 +33,26 @@ def create_integrator(self, energy): wavelength=1.239842e-6/energy ) - def optimization_func(self, x, image, integrator, - q_min, q_max, + def optimization_func(self, x, image, obj, energy, + q_min, q_max, pbar=None, chi_min=-179, chi_max=179, num_points=150): - poni1, poni2 = x - integrator.poni1 = poni1 - integrator.poni2 = poni2 - _, image_int = integrator.integrate_radial(image, num_points, + #print(f'evaluating objective at {x}') + obj.create_integrator(energy,poni1=x[0],poni2=x[1]) + _, image_int = obj.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): + num_points=150, max_iter = 150, + verbose=False, ): if not poni1_guess: poni1_guess = self._obj.poni1 poni2_guess = self._obj.poni2 @@ -55,15 +62,20 @@ def refine_geometry(self, energy, q_min, q_max, if not bounds: bounds = [(poni1_guess*0.9, poni1_guess*1.1), (poni2_guess*0.9, poni2_guess*1.1)] - image = self._obj.unstack('system').sel(energy=energy).values - res = minimize(self.optimization_func, (poni1_guess, poni2_guess), - args=(image, self.integrator, - q_min, q_max, chi_min, chi_max, - num_points), - bounds=bounds, method=method) + 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: - print('Optimization successful. Updating beamcenter') + 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] else: From 7a9b3a5f70a012182cec641c69ccd4e952a00e52 Mon Sep 17 00:00:00 2001 From: Peter Beaucage Date: Fri, 7 Apr 2023 14:19:00 -0400 Subject: [PATCH 3/5] refactor: move DrawMask and Check into accessors --- src/PyHyperScattering/BeamCentering.py | 12 +- src/PyHyperScattering/DrawMask.py | 141 ++++++++++++++++++ .../{IntegrationUtils.py => GeoCheck.py} | 86 +++++++++-- src/PyHyperScattering/util.py | 3 +- 4 files changed, 222 insertions(+), 20 deletions(-) create mode 100644 src/PyHyperScattering/DrawMask.py rename src/PyHyperScattering/{IntegrationUtils.py => GeoCheck.py} (69%) diff --git a/src/PyHyperScattering/BeamCentering.py b/src/PyHyperScattering/BeamCentering.py index 3eefc27a..2215603f 100644 --- a/src/PyHyperScattering/BeamCentering.py +++ b/src/PyHyperScattering/BeamCentering.py @@ -6,7 +6,7 @@ from tqdm.auto import tqdm -@xr.register_dataarray_accessor('util') +@xr.register_dataarray_accessor('beamcenter') class CenteringAccessor: def __init__(self, xr_obj): @@ -75,8 +75,12 @@ def refine_geometry(self, energy, q_min, q_max, bounds=bounds, method=method, options=options) if res.success: - 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] + 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] 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 f792066c..bd1ddab1 100644 --- a/src/PyHyperScattering/util.py +++ b/src/PyHyperScattering/util.py @@ -1,7 +1,8 @@ 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 From 95ee4ade92a009943fbabd1b9fd2c40247b2ea6a Mon Sep 17 00:00:00 2001 From: Peter Beaucage Date: Mon, 10 Apr 2023 12:27:41 -0400 Subject: [PATCH 4/5] Bugfix: actually use refined integrator in objective function --- src/PyHyperScattering/BeamCentering.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/PyHyperScattering/BeamCentering.py b/src/PyHyperScattering/BeamCentering.py index 2215603f..806bada2 100644 --- a/src/PyHyperScattering/BeamCentering.py +++ b/src/PyHyperScattering/BeamCentering.py @@ -38,8 +38,8 @@ def optimization_func(self, x, image, obj, energy, chi_min=-179, chi_max=179, num_points=150): #print(f'evaluating objective at {x}') - obj.create_integrator(energy,poni1=x[0],poni2=x[1]) - _, image_int = obj.integrator.integrate_radial(image, num_points, + 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') From 8e8b6dbfa1c78d2197dd48034a7dc2e7462f3a2b Mon Sep 17 00:00:00 2001 From: Andrew Levin Date: Wed, 4 Sep 2024 16:19:49 -0400 Subject: [PATCH 5/5] return refined beamcenter result --- src/PyHyperScattering/BeamCentering.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/PyHyperScattering/BeamCentering.py b/src/PyHyperScattering/BeamCentering.py index 806bada2..56c9e316 100644 --- a/src/PyHyperScattering/BeamCentering.py +++ b/src/PyHyperScattering/BeamCentering.py @@ -82,5 +82,7 @@ def refine_geometry(self, energy, q_min, q_max, 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')