Source code for pixstem.fem_tools

import numpy as np
from dask import delayed
from tqdm import tqdm
from matplotlib.pyplot import subplots as pltsubplots
from hyperspy.signals import Signal1D
from hyperspy.io import load as hsload
import pixstem

'''
Plotting function used to display output of FEM calculations
'''


[docs]def fem_calc(s, centre_x, centre_y, show_progressbar=True): """ Parameters ---------- s : PixelatedSTEM signal centre_x, centre_y : scalar show_progressbar : bool, optional 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). Example ------- >>> s = ps.dummy_data.get_fem_signal() >>> fem_results = s.fem_analysis( ... centre_x=50, centre_y=50, ... show_progressbar=False) >>> fem_results['V-Omegak'].plot() """ offset = False if s.data.min() == 0: s.data += 1 # To avoid division by 0 offset = True results = {} 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) 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'] = Signal1D( Vrkdask.compute(progressbar=show_progressbar)) results['Vrek'] = 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'] = Signal1D(results['Vrk']) results['Vrek'] = 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 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)) fig, axes = pltsubplots(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') 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] = hsload(rootname + '_FEM_' + i + '.hdf5') return results