%pylab inline
import hyperspy.api as hs
from statsmodels.nonparametric.smoothers_lowess import lowess
def DespikeEELS(spec, m=3., interp_outliers=True):
# Make a copy of the data
E = spec.axes_manager[0].axis.copy()
I = spec.data.copy()
# Lowess smooth.
smooth = lowess(I, E, frac=0.02)
# Using median, find inliers and outliers as # sigma away.
dist = np.abs(I-smooth[:,1])
d = np.abs(dist - np.median(dist))
mdev = np.median(d)
s = d/mdev if mdev else 0.
# Which points are inliers and which are out?
inliers = np.where(s<=m)[0]
outliers = np.where(s>m)[0]
if interp_outliers == False:
# Keep the inliers only (drop outliers).
E = E[inliers]
I = I[inliers]
spec2 = spec.deepcopy()
spec2.axes_manager[0].axis = E
spec2.data = I
return spec2
else:
# Interpolate the outliers to the lowess curve (retains same spectrum dimensions).
E[outliers] = smooth[outliers,0]
I[outliers] = smooth[outliers,1]
spec2 = spec.deepcopy()
spec2.axes_manager[0].axis = E
spec2.data = I
return spec2
def ProcessEELS(CoreLoss=None, LowLoss=None, BkgRange=None, DespikeThreshhold=None, NormRange=None):
# Load the core and low loss spectra.
cl = hs.load(CoreLoss)
if len(cl) > 1:
cl = cl[0]
if LowLoss is not None:
ll = hs.load(LowLoss)
# Align to ZLP.
ll.align_zero_loss_peak(also_align=[cl], subpixel=True)
# Estimate the thickness.
ll.estimate_thickness(threshold=5.).data
if DespikeThreshhold is not None:
cl = DespikeEELS(cl, m=DespikeThreshhold)
if BkgRange is not None:
x = cl.remove_background(background_type='PowerLaw',signal_range=BkgRange, fast=False)
cl.data = x
if LowLoss is not None:
deconv = cl.fourier_ratio_deconvolution(ll)
if NormRange is not None:
cl.data -= np.min(cl.data)
cl.data /= np.mean(cl.data[(cl.axes_manager[0].axis > NormRange[0]) & (cl.axes_manager[0].axis < NormRange[1])])
deconv.data -= np.min(deconv.data)
deconv.data /= np.mean(deconv.data[(deconv.axes_manager[0].axis > NormRange[0]) & (deconv.axes_manager[0].axis < NormRange[1])])
return ll, cl, deconv