Source code for pixstem.fem_tools

import matplotlib.pylab as plt
import numpy as np
from dask import delayed
from tqdm import tqdm
import hyperspy.api as hs
import pixstem


[docs]def fem_calc(s, centre_x=None, centre_y=None, show_progressbar=True): """Perform analysis of fluctuation electron microscopy (FEM) data as outlined in: T. L. Daulton, et al., Ultramicroscopy 110 (2010) 1279-1289. doi:10.1016/j.ultramic.2010.05.010 Parameters ---------- s : PixelatedSTEM Signal on which FEM analysis was performed centre_x, centre_y : int, optional All the diffraction patterns assumed to have the same centre position. show_progressbar : bool Default True Returns ------- results : Python dictionary Results of FEM data analysis, including the normalized variance of the annular mean (V-Omegak), mean of normalized variances of rings (V-rk), normalized variance of ring ensemble (Vrek), the normalized variance image (Omega-Vi), and annular mean of the variance image (Omega-Vk). Examples -------- >>> import pixstem.dummy_data as dd >>> import pixstem.fem_tools as femt >>> s = dd.get_fem_signal() >>> fem_results = femt.fem_calc( ... s, ... centre_x=128, ... centre_y=128, ... show_progressbar=False) >>> fem_results['V-Omegak'].plot() """ offset = False if centre_x is None: centre_x = np.int(s.axes_manager.signal_shape[0]/2) if centre_y is None: centre_y = np.int(s.axes_manager.signal_shape[1]/2) if s.data.min() == 0: s.data += 1 # To avoid division by 0 offset = True results = dict() results['RadialInt'] = ( s.radial_integration(centre_x=centre_x, centre_y=centre_y, normalize=False, show_progressbar=show_progressbar)) radialavgs = s.radial_integration(centre_x=centre_x, centre_y=centre_y, normalize=True, show_progressbar=show_progressbar) if radialavgs.data.min() == 0: radialavgs.data += 1 results['V-Omegak'] = ((radialavgs ** 2).mean() / (radialavgs.mean()) ** 2) - 1 results['RadialAvg'] = radialavgs.mean() if s._lazy: results['Omega-Vi'] = ((s ** 2).mean() / (s.mean()) ** 2) - 1 results['Omega-Vi'].compute(progressbar=show_progressbar) results['Omega-Vi'] = pixstem.pixelated_stem_class.PixelatedSTEM( results['Omega-Vi']) results['Omega-Vk'] = results['Omega-Vi'].radial_integration( centre_x=centre_x, centre_y=centre_y, normalize=True, show_progressbar=show_progressbar) oldshape = None if len(s.data.shape) == 4: oldshape = s.data.shape s.data = s.data.reshape( s.data.shape[0] * s.data.shape[1], s.data.shape[2], s.data.shape[3]) y, x = np.indices(s.data.shape[-2:]) r = np.sqrt((x - centre_x) ** 2 + (y - centre_y) ** 2) r = r.astype(np.int) nr = np.bincount(r.ravel()) Vrklist = [] Vreklist = [] for k in tqdm(range(0, len(nr)), disable=(not show_progressbar)): locs = np.where(r == k) vals = s.data.vindex[:, locs[0], locs[1]].T Vrklist.append(np.mean((np.mean(vals ** 2, 1) / np.mean(vals, 1) ** 2) - 1)) Vreklist.append(np.mean(vals.ravel() ** 2) / np.mean(vals.ravel()) ** 2 - 1) Vrkdask = delayed(Vrklist) Vrekdask = delayed(Vreklist) results['Vrk'] = hs.signals.Signal1D( Vrkdask.compute(progressbar=show_progressbar)) results['Vrek'] = hs.signals.Signal1D( Vrekdask.compute(progressbar=show_progressbar)) else: results['Omega-Vi'] = ((s ** 2).mean() / (s.mean()) ** 2) - 1 results['Omega-Vk'] = results['Omega-Vi'].radial_integration( centre_x=centre_x, centre_y=centre_y, normalize=True, show_progressbar=show_progressbar) oldshape = None if len(s.data.shape) == 4: oldshape = s.data.shape s.data = s.data.reshape( s.data.shape[0] * s.data.shape[1], s.data.shape[2], s.data.shape[3]) y, x = np.indices(s.data.shape[-2:]) r = np.sqrt((x - centre_x) ** 2 + (y - centre_y) ** 2) r = r.astype(np.int) nr = np.bincount(r.ravel()) results['Vrk'] = np.zeros(len(nr)) results['Vrek'] = np.zeros(len(nr)) for k in tqdm(range(0, len(nr)), disable=(not show_progressbar)): locs = np.where(r == k) vals = s.data[:, locs[0], locs[1]] results['Vrk'][k] = np.mean( (np.mean(vals ** 2, 1) / np.mean(vals, 1) ** 2) - 1) results['Vrek'][k] = np.mean( vals.ravel() ** 2) / np.mean(vals.ravel()) ** 2 - 1 results['Vrk'] = hs.signals.Signal1D(results['Vrk']) results['Vrek'] = hs.signals.Signal1D(results['Vrek']) if oldshape: s.data = s.data.reshape(oldshape) if offset: s.data -= 1 # Undo previous addition of 1 to input data return results
[docs]def plot_fem(s, results, lowcutoff=10, highcutoff=120, k_cal=None): """Produce a summary plot of all calculated FEM parameters Parameters ---------- s : PixelatedSTEM Signal on which FEM analysis was performed results : Dictionary Series of HyperSpy Signals containing results of FEM analysis performed on s lowcutoff : integer Position of low-q cutoff for plots highcutoff : integer Position of high-q cutoff for plots k_cal : float or None Reciprocal space unit of length per pixel in inverse Angstroms Returns ------- fig : Matplotlib Figure Examples -------- >>> s = ps.dummy_data.get_fem_signal() >>> import pixstem.fem_tools as femt >>> fem_results = s.fem_analysis( ... centre_x=50, ... centre_y=50, ... show_progressbar=False) >>> fig = femt.plot_fem(s, fem_results, 10, 120) """ if k_cal: xaxis = 2 * np.pi * k_cal * \ np.arange(0, len(results['RadialAvg'].data)) else: xaxis = np.arange(0, len(results['RadialAvg'].data)) if highcutoff > len(results['RadialAvg'].data): highcutoff = len(results['RadialAvg'].data) - 1 fig, axes = plt.subplots(3, 2, figsize=(9, 12)) axes[0, 0].imshow(np.log(s.mean().data + 1), cmap='viridis') axes[0, 0].set_title('Mean Pattern', size=15) axes[0, 0].set_yticks([]) axes[0, 0].set_xticks([]) axes[0, 1].plot(xaxis[lowcutoff:highcutoff], results['RadialAvg'].data[lowcutoff:highcutoff], linestyle='', marker='o') axes[0, 1].set_title('Mean Radial Profile', size=15) axes[0, 1].set_ylabel('Integrated Intensity (counts)') axes[0, 1].set_xlabel(r'k ($\AA^{-1}$)') axes[1, 0].plot(xaxis[lowcutoff:highcutoff], results['V-Omegak'].data[lowcutoff:highcutoff], linestyle='', marker='o') axes[1, 0].set_ylabel(r'V$_\Omega$(k)', size=15) axes[1, 0].set_xlabel(r'k ($\AA^{-1}$)') axes[1, 1].plot(xaxis[lowcutoff:highcutoff], results['Vrk'].data[lowcutoff:highcutoff], linestyle='', marker='o', label=r'$\overline{V}_r$(k)') axes[1, 1].plot(xaxis[lowcutoff:highcutoff], results['Vrek'].data[lowcutoff:highcutoff], linestyle='', marker='o', label=r'V$_{re}$(k)') axes[1, 1].legend() axes[1, 1].set_ylabel(r'$\overline{V}_r$(k),V$_{re}$(k)', size=15) axes[1, 1].set_xlabel(r'k ($\AA^{-1}$)') axes[2, 0].imshow(results['Omega-Vi'], cmap='viridis') axes[2, 0].set_title('Variance Image', size=15) axes[2, 0].set_yticks([]) axes[2, 0].set_xticks([]) axes[2, 1].plot(xaxis[lowcutoff:highcutoff], results['Omega-Vk'].data[lowcutoff:highcutoff], linestyle='', marker='o') axes[2, 1].set_ylabel(r'$\Omega_V$(k)', size=15) axes[2, 1].set_xlabel(r'k ($\AA^{-1}$)') fig.tight_layout() return fig
[docs]def save_fem(results, rootname): """Save dictionary of FEM results to individual HDF5 files. Parameters ---------- results : Dictionary Results of FEM analysis performed on s rootname : str Root name for output files. '_FEM_Parameter.hdf5' will be appended Examples -------- >>> import pixstem.fem_tools as femt >>> s = ps.dummy_data.get_fem_signal() >>> fem_results = s.fem_analysis( ... centre_x=50, ... centre_y=50, ... show_progressbar=False) >>> femt.save_fem(fem_results,'TestData') """ for i in results: results[i].save(rootname + '_FEM_' + i + '.hdf5')
[docs]def load_fem(rootname): """Load individual HDF5 files containing previously saved FEM parameters to dictionary. Parameters ---------- rootname : str Root name for output files. '_FEM_Parameter.hdf5' will be appended Returns --------- results : Dictionary Series of HyperSpy signals previously saved using saveFEM Examples -------- >>> import pixstem.fem_tools as femt >>> s = ps.dummy_data.get_fem_signal() >>> fem_results = s.fem_analysis( ... centre_x=50, ... centre_y=50, ... show_progressbar=False) >>> femt.save_fem(fem_results,'TestData') >>> new_results = femt.load_fem('TestData') """ results = {} keys = ['Omega-Vi', 'Omega-Vk', 'RadialAvg', 'RadialInt', 'V-Omegak', 'Vrek', 'Vrk'] for i in keys: results[i] = hs.load(rootname + '_FEM_' + i + '.hdf5') return results