Tutorial for the HyperSpy For Electron Microscopy Data Analysis Workshop at the European Microscopy Conference
Copenhagen, 25 August, 2024
This workshop closely follows the demo notebook designed by Jordi Ferrer-Orri and Jonas Lähnemann, who both have my thanks.
Table of Contents:
We import the public functions (api = application programming interface) of HyperSpy.
Some functionalities of LumiSpy are directly available if the package is installed, but we can separately load it to access extra utilities.
Finally, matplotlib.pyplot provides some additional plotting functions and numpy numerical operations on arrays that we will use:
%matplotlib widget
# Use '%matplotlib widget' in JupterLab and '%matplotlib notebook' in JupyterNotebook for interactive inline functionality (e.g. on binder)
# For pop-up window plots on your local computer, use '%matplotlib tk' or '%matplotlib qt' instead
import hyperspy.api as hs
import lumispy as lum
import matplotlib.pyplot as plt
import numpy as np
For saving analyses, HyperSpy has its own hdf5-based data format .hspy.
HyperSpy supports a wide range of microscopy (and spectroscopy) related data file types:
.dm3/.dm4.sur (For the moment does not parse original_metadata to metadata).hdf5 (Still needs a parsing script at the moment that you can obtain from Delmic, will be contributed to RosettaSciIO).tif format (Reading the itex .img format will be available in RosettaSciIO)Further spectroscopy-related formats will be available soon (through RosettaSciIO):
.xml.tvf.wdfWe will load three files that we will use during the demo. A preprocessed dataset saved in the hspy format and two maps in the dm4 format:
We assume the file location as in the demo repository, if you downloaded the notebook and the data files individually, you might need to adapt the path.
You can also leave the path empty. A pop-up window will appear to select.
cl1 = hs.load("demo-files/01/01_demo.hspy")
cl2 = hs.load("demo-files/load_from_GatanFiles/asymmetric-peak_map.dm4")
cl3 = hs.load("demo-files/load_from_GatanFiles/quantum-emitters_map.dm4")
To see parameters that the function takes, in Jupyter, you can display the docstring by using a ?:
hs.load?
Signature: hs.load( filenames=None, signal_type=None, stack=False, stack_axis=None, new_axis_name='stack_element', lazy=False, convert_units=False, escape_square_brackets=False, stack_metadata=True, load_original_metadata=True, show_progressbar=None, **kwds, ) Docstring: Load potentially multiple supported files into HyperSpy. Supported formats: hspy (HDF5), msa, Gatan dm3, Ripple (rpl+raw), Bruker bcf and spx, FEI ser and emi, SEMPER unf, EMD, EDAX spd/spc, CEOS prz tif, and a number of image formats. Depending on the number of datasets to load in the file, this function will return a HyperSpy signal instance or list of HyperSpy signal instances. Any extra keywords are passed to the corresponding reader. For available options, see their individual documentation. Parameters ---------- filenames : None, (list of) str or (list of) pathlib.Path, default None The filename to be loaded. If None, a window will open to select a file to load. If a valid filename is passed, that single file is loaded. If multiple file names are passed in a list, a list of objects or a single object containing multiple datasets, a list of signals or a stack of signals is returned. This behaviour is controlled by the `stack` parameter (see below). Multiple files can be loaded by using simple shell-style wildcards, e.g. 'my_file*.msa' loads all the files that start by 'my_file' and have the '.msa' extension. Alternatively, regular expression type character classes can be used (e.g. ``[a-z]`` matches lowercase letters). See also the `escape_square_brackets` parameter. signal_type : None, str, default None The acronym that identifies the signal type. May be any signal type provided by HyperSpy or by installed extensions as listed by `hs.print_known_signal_types()`. The value provided may determines the Signal subclass assigned to the data. If None (default), the value is read/guessed from the file. Any other value would override the value potentially stored in the file. For example, for electron energy-loss spectroscopy use 'EELS'. If '' (empty string) the value is not read from the file and is considered undefined. stack : bool, default False Default False. If True and multiple filenames are passed, stacking all the data into a single object is attempted. All files must match in shape. If each file contains multiple (N) signals, N stacks will be created, with the requirement that each file contains the same number of signals. stack_axis : None, int or str, default None If None (default), the signals are stacked over a new axis. The data must have the same dimensions. Otherwise, the signals are stacked over the axis given by its integer index or its name. The data must have the same shape, except in the dimension corresponding to `axis`. new_axis_name : str, optional The name of the new axis (default 'stack_element'), when `axis` is None. If an axis with this name already exists, it automatically appends '-i', where `i` are integers, until it finds a name that is not yet in use. lazy : bool, default False Open the data lazily - i.e. without actually reading the data from the disk until required. Allows opening arbitrary-sized datasets. convert_units : bool, default False If True, convert the units using the `convert_to_units` method of the `axes_manager`. If False, does nothing. escape_square_brackets : bool, default False If True, and ``filenames`` is a str containing square brackets, then square brackets are escaped before wildcard matching with ``glob.glob()``. If False, square brackets are used to represent character classes (e.g. ``[a-z]`` matches lowercase letters). stack_metadata : {bool, int} If integer, this value defines the index of the signal in the signal list, from which the ``metadata`` and ``original_metadata`` are taken. If ``True``, the ``original_metadata`` and ``metadata`` of each signals are stacked and saved in ``original_metadata.stack_elements`` of the returned signal. In this case, the ``metadata`` are copied from the first signal in the list. If False, the ``metadata`` and ``original_metadata`` are not copied. show_progressbar : None or bool If ``True``, display a progress bar. If ``None``, the default from the preferences settings is used. Only used with ``stack=True``. load_original_metadata : bool, default True If ``True``, all metadata contained in the input file will be added to ``original_metadata``. This does not affect parsing the metadata to ``metadata``. reader : None, str, module, optional Specify the file reader to use when loading the file(s). If None (default), will use the file extension to infer the file type and appropriate reader. If str, will select the appropriate file reader from the list of available readers in HyperSpy. If module, it must implement the ``file_reader`` function, which returns a dictionary containing the data and metadata for conversion to a HyperSpy signal. print_info: bool, optional For SEMPER unf- and EMD (Berkeley)-files. If True, additional information read during loading is printed for a quick overview. Default False. downsample : int (1–4095), optional For Bruker bcf files, if set to integer (>=2) (default 1), bcf is parsed into down-sampled size array by given integer factor, multiple values from original bcf pixels are summed forming downsampled pixel. This allows to improve signal and conserve the memory with the cost of lower resolution. cutoff_at_kV : None, int, float, optional For Bruker bcf files and Jeol, if set to numerical (default is None), hypermap is parsed into array with depth cutoff at set energy value. This allows to conserve the memory by cutting-off unused spectral tails, or force enlargement of the spectra size. Bruker bcf reader accepts additional values for semi-automatic cutoff. "zealous" value truncates to the last non zero channel (this option should not be used for stacks, as low beam current EDS can have different last non zero channel per slice). "auto" truncates channels to SEM/TEM acceleration voltage or energy at last channel, depending which is smaller. In case the hv info is not there or hv is off (0 kV) then it fallbacks to full channel range. select_type : 'spectrum_image', 'image', 'single_spectrum', None, optional If None (default), all data are loaded. For Bruker bcf and Velox emd files: if one of 'spectrum_image', 'image' or 'single_spectrum', the loader returns either only the spectrum image, only the images (including EDS map for Velox emd files), or only the single spectra (for Velox emd files). first_frame : int, optional Only for Velox emd files: load only the data acquired after the specified fname. Default 0. last_frame : None, int, optional Only for Velox emd files: load only the data acquired up to specified fname. If None (default), load the data up to the end. sum_frames : bool, optional Only for Velox emd files: if False, load each EDS frame individually. Default is True. sum_EDS_detectors : bool, optional Only for Velox emd files: if True (default), the signals from the different detectors are summed. If False, a distinct signal is returned for each EDS detectors. rebin_energy : int, optional Only for Velox emd files: rebin the energy axis by the integer provided during loading in order to save memory space. Needs to be a multiple of the length of the energy dimension (default 1). SI_dtype : numpy.dtype, None, optional Only for Velox emd files: set the dtype of the spectrum image data in order to save memory space. If None, the default dtype from the Velox emd file is used. load_SI_image_stack : bool, optional Only for Velox emd files: if True, load the stack of STEM images acquired simultaneously as the EDS spectrum image. Default is False. dataset_path : None, str, list of str, optional For filetypes which support several datasets in the same file, this will only load the specified dataset. Several datasets can be loaded by using a list of strings. Only for EMD (NCEM) and hdf5 (USID) files. stack_group : bool, optional Only for EMD NCEM. Stack datasets of groups with common name. Relevant for emd file version >= 0.5 where groups can be named 'group0000', 'group0001', etc. ignore_non_linear_dims : bool, optional Only for HDF5 USID files: if True (default), parameters that were varied non-linearly in the desired dataset will result in Exceptions. Else, all such non-linearly varied parameters will be treated as linearly varied parameters and a Signal object will be generated. only_valid_data : bool, optional Only for FEI emi/ser files in case of series or linescan with the acquisition stopped before the end: if True, load only the acquired data. If False, fill empty data with zeros. Default is False and this default value will change to True in version 2.0. Returns ------- (list of) :class:`~.api.signals.BaseSignal` or subclass Examples -------- Loading a single file providing the signal type: >>> d = hs.load('file.dm3', signal_type="EDS_TEM") # doctest: +SKIP Loading multiple files: >>> d = hs.load(['file1.hspy','file2.hspy']) # doctest: +SKIP Loading multiple files matching the pattern: >>> d = hs.load('file*.hspy') # doctest: +SKIP Loading multiple files containing square brackets in the filename: >>> d = hs.load('file[*].hspy', escape_square_brackets=True) # doctest: +SKIP Loading multiple files containing character classes (regular expression): >>> d = hs.load('file[0-9].hspy') # doctest: +SKIP Loading (potentially larger than the available memory) files lazily and stacking: >>> s = hs.load('file*.blo', lazy=True, stack=True) # doctest: +SKIP Specify the file reader to use >>> s = hs.load('a_nexus_file.h5', reader='nxs') # doctest: +SKIP Loading a file containing several datasets: >>> s = hs.load("spameggsandham.nxs") # doctest: +SKIP >>> s # doctest: +SKIP [<Signal1D, title: spam, dimensions: (32,32|1024)>, <Signal1D, title: eggs, dimensions: (32,32|1024)>, <Signal1D, title: ham, dimensions: (32,32|1024)>] Use list indexation to access single signal >>> s[0] # doctest: +SKIP <Signal1D, title: spam, dimensions: (32,32|1024)> File: c:\users\gkusc\anaconda3\lib\site-packages\hyperspy\io.py Type: function
Each HyperSpy signal object has certain attributes that contain the relevant data about the axes, data and metadata.
To understand the HyperSpy datastructure, lets have a look at the dataset cl2 (Gatan file).
As LumiSpy is installed, the dataset is directly recognized as CL data and the signal_type set to CLSpectrum. (The fallback would be the more generic Signal1D if LumiSpy is not installed).
The signal class provides certain specific routines, for example conversion to energy axis in the case of luminescence data.
Our sample dataset has two navigation dimensions and one signal (spectral) dimension:
cl2
<CLSpectrum, title: 880_30K_map2a, dimensions: (30, 16|334)>
The information about the axes is stored in the axes_manager. Thus, we can get more details about the different axes, by calling the axes manager:
cl2.axes_manager
< Axes manager, axes: (30, 16|334) >
| Navigation axis name | size | index | offset | scale | units |
|---|---|---|---|---|---|
| x | 30 | 0 | -0.0 | 0.04001351818442345 | µm |
| y | 16 | 0 | -0.0 | 0.04001351818442345 | µm |
| Signal axis name | size | offset | scale | units | |
|---|---|---|---|---|---|
| Wavelength | 334 | 530.7898101451792 | 0.819013774394989 | nm |
The actual data (signal intensity) is stored in a numpy array:
cl2.data
array([[[ 5.125 , 4.625 , 4.5 , ..., 2.875 , 3. , -3.375 ],
[ 7.125 , 7.1875, 7. , ..., 4.3125, 3.75 , -1.5625],
[ 7.4375, 7.625 , 8.25 , ..., 6.3125, 6.8125, -0.0625],
...,
[11.6875, 11.25 , 12.4375, ..., 10.8125, 10.25 , 4.625 ],
[12.125 , 12. , 11.5 , ..., 12.25 , 11.3125, 5.3125],
[11.9375, 11.9375, 11.8125, ..., 12.375 , 12.125 , 4. ]],
[[ 8.25 , 8.5 , 8.375 , ..., 7.3125, 6.0625, 0.25 ],
[11.125 , 9.125 , 10.0625, ..., 9.4375, 9.875 , 3.4375],
[10.5 , 11.4375, 10.4375, ..., 11.3125, 9.4375, 4.4375],
...,
[14. , 13.8125, 13.25 , ..., 13.75 , 11.875 , 5.5625],
[13.75 , 15.0625, 14.4375, ..., 14.0625, 13.0625, 5.375 ],
[14.625 , 14.1875, 14.5625, ..., 13.0625, 12.75 , 5.8125]],
[[ 8.5 , 9.4375, 9.0625, ..., 8.75 , 7.75 , 1.8125],
[11.3125, 11.3125, 10.4375, ..., 10.3125, 9.875 , 3.1875],
[12. , 11.625 , 12.5 , ..., 10.875 , 10.8125, 4.875 ],
...,
[13.4375, 14.3125, 14.5625, ..., 15.375 , 12.8125, 5.9375],
[15. , 15.375 , 14.5 , ..., 14.1875, 13.1875, 6.5 ],
[14.6875, 15.4375, 14.5 , ..., 14. , 13.25 , 5.8125]],
...,
[[12.125 , 11.0625, 11.6875, ..., 10.5625, 10.9375, 9.5625],
[13.0625, 13.3125, 12.6875, ..., 11.75 , 10.9375, 5.875 ],
[14.125 , 14.375 , 13.8125, ..., 12.6875, 12.625 , 5.9375],
...,
[17.1875, 16.6875, 17.125 , ..., 16.375 , 15.1875, 8.625 ],
[16.625 , 16.5625, 16.0625, ..., 16.625 , 15.3125, 9.125 ],
[16.5 , 16.4375, 16.75 , ..., 16.5625, 14.875 , 8. ]],
[[11.75 , 11.0625, 10.4375, ..., 10.0625, 10.0625, 8.75 ],
[13.1875, 13.125 , 12.8125, ..., 12.9375, 10.8125, 5.1875],
[15.0625, 13.8125, 13. , ..., 13. , 12.1875, 5.75 ],
...,
[16.3125, 16.625 , 16.9375, ..., 15.5625, 15.6875, 8.0625],
[17.1875, 17.3125, 16.375 , ..., 16.75 , 14.5625, 8.5625],
[15.5625, 15.5 , 16.625 , ..., 15.9375, 15.4375, 8.75 ]],
[[11.0625, 11.0625, 11.0625, ..., 10.4375, 10.8125, 5.75 ],
[13.4375, 13.5 , 11.9375, ..., 12.4375, 12.9375, 3.9375],
[14.1875, 13.875 , 13.4375, ..., 12.5 , 11.4375, 5.375 ],
...,
[16.0625, 15.9375, 17.1875, ..., 15.1875, 15.4375, 7.4375],
[16.1875, 16.5 , 15.6875, ..., 16.0625, 15.5625, 7.875 ],
[16.8125, 17.5 , 16.625 , ..., 16.1875, 15.875 , 7.8125]]],
dtype=float32)
For most supported file formats, the metadata is automatically parsed into HyperSpy's metadata tree. It contains information about the measurement, but potentially also about post-processing:
cl2.metadata
In a separate tree, the complete metadata from the vendor format is read in (which follows different conventions depending on the format):
cl2.original_metadata
We can easily plot and explore the hyperspectral data:
(In the following, we will use the preprocessed dataset cl1. The sample contains MethylammoniumLead Bromine (MAPbBr3) perovskite single crystals fabricated by Alice Dearle.)
cl1.plot()
HBox(children=(Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Ba…
Plot the average CL spectrum of the whole map:
cl1.mean().plot()
HyperSpy has a powerful numpy (Matlab) style indexing mechanism that distinguishes between navigation and signal axes:
.inav[x1:x2,y1:y2].isig[s1:s2]The index parameters can be either:
For example, we can either plot a subset of the map in navigation space (selected using pixels as index):
cl1.inav[2:23,0:20].plot()
Or, we can plot the mean spectrum in a certain spectral range (selected using wavelength units):
cl1.isig[440.:600.].mean().plot()
Indexing can also be used for color-filtered (chromatic) imaging.
First, lets plot the panchromatic image:
(the object is transposed, so that we plot the intensity over navigation instead of signal dimensions)
cl1.T.mean().plot(cmap='viridis')
Secondly we can generate a monochromatic image by changing the navigation and signal axis and moving the red slider to the desired wavelength.
cl1.T.plot(cmap='viridis')
Now, we can plot the intensity in a selected spectral window (color-filtered image) using indexing:
cl1.isig[480.:550.].T.mean().plot(cmap='viridis')
Alternatively, we can interactively select a spectral window (color-filtered image) using regions of interest:
im = cl1.T
im.plot()
roi1 = hs.roi.SpanROI(left=455, right=485) #sets a digitalbandfilter
im_roi1 = roi1.interactive(im, color="red")
im_roi1_mean = hs.interactive(im_roi1.mean,
event=roi1.events.changed,
recompute_out_event=None)
im_roi1_mean.plot(cmap='viridis')
Required versions: hyperspy>=1.7.0 and lumispy>=0.2.0
HyperSpy has different types of axes:
UniformDataAxis is defined through an offset and a scale (delta between pixels)FunctionalDataAxis is defined through a UniformDataAxis and a function to convert the valuesDataAxis is defined through an axis vector/arrayThe wavelength scale of our sample object is a UniformDataAxis:
cl1.axes_manager
LumiSpy provides easy conversions of the signal axis to the energy scale:
It can either replace the axis in the existing object (default) or create a copy of the signal object with the new axis (inplace=False):
cl1_eV = cl1.to_eV(inplace=False)
The axes manager now contains an energy axis, which is a non-uniform DataAxis:
cl1_eV.axes_manager
To explore the data in the energy domain, we again plot the signal:
cl1_eV.plot()
To preserve the integrated intensity per spectral window, a Jacobian transformation has to be applied to the signal intensity:
As we require $I(E)dE = I(\lambda)d\lambda$, the scale transformation $E=hc/\lambda$ implies
$$I(E) = I(\lambda)\frac{d\lambda}{dE} = I(\lambda)\frac{d}{dE} \frac{h c}{E} = - I(\lambda) \frac{h c}{E^2}$$(where the minus sign just reflects the different directions of integration in the wavelength and energy domains)
This transformation is the default in LumiSpy, but can be deactivated by setting jacobian=False.
To visualize the effect of the Jacobian transformation, we can plot a signal with constant intensity before and after the transformation:
# Create a model signal with linear intensity
axis = {'offset': 300, 'scale': 4, 'units': 'nm', 'size': 101, 'name': 'Wavelength'}
s = hs.signals.Signal1D(np.ones(101), axes=[axis])
s.set_signal_type("Luminescence")
s2 = s.to_eV(inplace=False)
# Some additional arrays to help with visualizing the spectral bins during plotting
x = np.arange(9)*50+300
x2 = lum.nm2eV(x)
y2 = hs.signals.Signal1D(np.ones(9), axes=[{'offset': 300, 'scale': 50, 'size': 9,}])
y2.set_signal_type("Luminescence")
y2.to_eV()
# Plot comparative figures
fig1 = plt.figure(figsize=(10,4))
ax0 = plt.subplot(121)
plt.ylim(0,1.3)
plt.xlabel('Wavelength (nm)')
ax0.plot(s.axes_manager[0].axis,s.data,color='orange')
ax0.vlines(x,0,1,color='orange')
ax0.fill_between(s.axes_manager[0].axis,0,s.data, facecolor='orange', alpha=0.3)
ax1 = plt.subplot(122)
plt.ylim(0,0.43)
plt.xlabel('Energy (eV)')
ax1.plot(s2.axes_manager[0].axis,s2.data,color='orange')
ax1.vlines(x2,0,y2.data[::-1],color='orange')
ax1.fill_between(s2.axes_manager[0].axis,0,s2.data, facecolor='orange', alpha=0.3)
We will introduce basic fitting functionality using our sample dataset cl1, for more details see the Fitting_tutorial in the HyperSpy demos repository.
Note that HyperSpy has a range of built-in functions covering most needs that can be added as components to a model. However, it also has an intuitive mechanism to define custom functions.
First, we need to initialize the model:
m = cl1_eV.create_model()
Check the components of the model (should be empty, but for some types of signals, an automatic background component is added):
m.components
Now we need to create some components and add them to the model.
We will use a constant Offset and a Gaussian (defined through height and FWHM, hence GaussianHF):
bkg = hs.model.components1D.Offset()
g1 = hs.model.components1D.GaussianHF()
m.extend([g1, bkg])
m.components
To see the parameters of our components and their default values, we can print all parameter values:
m.print_current_values()
Lets set some sensible starting values for our components, for a position in the map where we know that there should be signal (as it is not the case everywhere for our sample dataset):
Note that we could also use the estimate parameters function of the GaussianHF component (commented out line). However, this approach does not work as well as using manual starting values for this rather noisy dataset where many pixels do not contain spectra, unless additional boundaries are set to the parameters below.
# g1.estimate_parameters(cl1_eV,2.3,2.5)
cl1.axes_manager.indices = (7,7)
g1.centre.value = 2.4 # Gaussian centre
g1.fwhm.value = 0.1 # Gaussian width
g1.height.value = 5 # Gaussian height
bkg.offset.value = 0.1 # background offset
We can also set boundaries (bmin and bmax) for some of the parameters:
g1.centre.bmax = g1.centre.value + 0.2
g1.centre.bmin = g1.centre.value - 0.2
g1.fwhm.bmin = 0.01
We can now fit the model at the chosen position, copy the result as starting value to all positions, and plot the result:
m.fit(bounded=True)
m.assign_current_values_to_all()
m.plot()
Again, we can also print the updated parameters:
m.print_current_values()
The model now has the result from our chosen pixel everywhere. Using this as optimized starting paramters, we can now fit all pixels. When plotting, we activate additional plotting of the individual components:
m.multifit(bounded=True, show_progressbar=True)
m.plot(plot_components=True)
To plot maps of the parameters of the Gaussian, we create signal objects from these datasets:
m_centre = g1.centre.as_signal()
m_centre.plot(cmap='bwr_r', centre_colormap=False) # Otherwise, it would be centred around `0` and we would see little difference between pixels
m_intensity = g1.height.as_signal()
m_intensity.plot(cmap='viridis')
We can use the fit model as basis to do a particle segmentation by creating a mask for all pixels, where the intensity is below the mean value:
mask_treshold = m_intensity.data.mean()
mask = m_intensity.data > mask_treshold #Returns a boolean matrix mask
plt.figure()
plt.imshow(mask)
We can now plot the previous graph of the centre-parametre, after applying the mask:
(m_centre * mask).plot(cmap='bwr_r', vmin=2.3, centre_colormap=False)
Working on the unprocessed dataset cl2, we can introduce some basic functions for artefact correction:
HyperSpy has an interactive tool for background removal that supports various functions, let's start by removing a simple offset:
Applycl2.plot()
cl2.remove_background(background_type="Offset")
The signal beyond 800 nm goes to negative values, so lets remove the last three pixels from every spectrum:
cl2 = cl2.isig[:-3]
cl2.plot()
There is also a tool for interactive removal of cosmic rays (pixels with sharp spikes), see Help for instructions.
In brief:
Inspect the derivative histogram
Set a sensible threshold to catch the outliers in the histogram
Iterate through Find next / Remove spike to continue for wrong identifications / remove identified spikes
Close when finished
NOTE: The interactive version does not work well with inline plotting (e.g. on binder)
cl2.spikes_removal_tool(interactive=False)
The current dataset is quite noisy. As the peak is broad in comparison with the spectral resolution, one way to improve that is by rebinning the data along the signal axis:
cl2 = cl2.rebin(scale=[1,1,2])
cl2.plot()
Additionally, HyperSpy provides three different functions for data smoothing:
smooth_lowess (lowess smoothing)smooth_savitzky_golay (Savitzky Golay filter)smooth_tv (total variation data smoothing)These functions can be run interactively to choose the right parameters, but the parameters can also be passed to the function:
cl2.plot()
cl2.smooth_lowess(number_of_iterations = 2)
cl2.plot()
In particular for asymmetric peaks, fitting might not always be the best way to determine peak characteristics (though asymmetric functions, such as the skew normal distribution are provided). Therefore HyperSpy provides a number of additional routines.
Peaks can be identified and characterized using the peak finder routine find_peaks1D_ohaver that is based on the downwards zero crossing of the first derivative.
As we have some side-peaks, we operate on a subrange of the wavelength axis.
peaks = cl2.isig[600.:].find_peaks1D_ohaver(maxpeakn=1)
The function returns a structured array that contains position, height and width for every pixel (potentially each for multiple peaks).
peaks[0,0]
You can also just determine the width of a peak directly from the signal without fitting a model to the data.
The default is to determine the FWHM, i.e. a factor=0.5. This value can be set to any other fraction of the peak height.
width = cl2.isig[600.:].estimate_peak_width(return_interval=True)
width[0].plot(cmap='viridis')
Especially for broad, asymmetric emission bands, the position of the maximum intensity might be of limited value. Therefore, LumiSpy provides an additional centroid function that determines the centre of mass of a peak.
Required version: lumispy>=0.2.2
Note that, as with fitting, it might make more sense to run these routines in the energy domain after a Jacobian transformation than to convert the result - in particular for broad emission bands.
com = cl2.isig[600.:].centroid()
com.plot(cmap='viridis')
Finally, to showcase the power of fitting with a programming language, lets fit the more complicated dataset cl3:
Fit a spectrum with a skew normal component on a broad Gaussian as background.
Then add a variable number of sharp Gaussians depending on the number of peaks found by the peak finder routine.
The presented dataset is a CL map of luminescent centers in h-BN on SiO$_2$, similar as discussed in Hernández-Mínguez et al., Phys. Rev. Appl. 10, 2331 (2018). The analysis is simplified for instructive purposes.
For reasons of calculation time, we run the fit on a single spectrum, but it could of course be looped across a complete spectral map - even with a variable number of peaks per spectrum. To apply functions defined for individual spectra on a complete spectral image, HyperSpy provides the map function.
cl3 = cl3.inav[10,8]
m = cl3.create_model()
g1 = hs.model.components1D.SkewNormal()
g2 = hs.model.components1D.GaussianHF()
m.extend((g1,g2))
g1.x0.value=655
g1.scale.value=50
g1.shape.value=2
g1.A.value=220
g1.x0.bmin=635
g1.x0.bmax=675
g1.x0.bounded=True
g2.centre.value=580
g2.fwhm.value=180
g2.height.value=20
m.fit()
m.plot()
def getPeaks(S2):
S2=S2.rebin(new_shape=[334])
S2.smooth_savitzky_golay(window_length=15,polynomial_order=2)
peaks = S2.find_peaks1D_ohaver(amp_thresh=0.2*np.max(S2.data),maxpeakn=10)[0]
return peaks
peaks
ng=10
S2 = cl3.deepcopy()
g = list()
for i in np.arange(ng):
g.append(hs.model.components1D.GaussianHF())
m.extend(g)
for i in np.arange(ng):
g[i].active_is_multidimensional = True
# The following code would need to be in a loop to run for a whole spectral image
# m.axes_manager.indices = cl3.axes_manager.indices
peaks = getPeaks(S2.inav[cl3.axes_manager.indices])
for i in np.arange(np.size(peaks)):
g[i].centre.value=peaks['position'][i]
g[i].centre.bmin=peaks['position'][i]-3
g[i].centre.bmax=peaks['position'][i]+3
g[i].centre.bounded=True
g[i].fwhm.value=5
g[i].fwhm.bmax=10
g[i].fwhm.bmin=1
g[i].fwhm.bounded=True
g[i].height.value=20
g[i].height.bmin=1
g[i].height.bounded=True
if np.size(peaks)<ng:
for i in np.arange(np.size(peaks),ng):
g[i].active = False
m.fit(bounded=True)
m.plot()