Workshops/230324_CLEBIC_HyperSpy-LumiSpy.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Introduction to *HyperSpy* and *LumiSpy*\n",
    "\n",
    "Tutorial for the **Workshop on Cathodoluminescence and Electron Beam Induced Current of Semiconductor Nanostructures**\n",
    "\n",
    "Grenoble, March 24, 2023\n",
    "\n",
    "**Table of Contents:**\n",
    "\n",
    "- [Import packages](#Import-packages)\n",
    "- [Loading files](#Loading-files)\n",
    "- [Data structure / Axes handling](#Data-structure-/-Axes-handling)\n",
    "- [Metadata](#Metadata)\n",
    "- [Plot / Explore](#Plot-/-Explore)\n",
    "- [Indexing](#Indexing)\n",
    "- [Axes types / Convert to energy scale](#Axes-types-/-Convert-to-energy-scale)\n",
    "- [Model fitting](#Model-fitting)\n",
    "- [Correction of spectral defects](#Correction-of-spectral-defects)\n",
    "- [Peak identification / Centroid / Peak width](#Peak-identification-/-Centroid-/-Peak-width)\n",
    "- [Advanced fitting example](#Advanced-fitting-example)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Import packages\n",
    "\n",
    "We import the public functions (api = application programming interface) of `HyperSpy`.\n",
    "\n",
    "Some functionalities of `LumiSpy` are directly available if the package is installed, but we can separately load it to access extra utilities.\n",
    "\n",
    "Finally, `matplotlib.pyplot` provides some additional plotting functions and `numpy` numerical operations on arrays that we will use:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib widget\n",
    "# Use '%matplotlib widget' in JupterLab and '%matplotlib notebook' in JupyterNotebook for interactive inline functionality (e.g. on binder)\n",
    "# For pop-up window plots on your local computer, use '%matplotlib tk' or '%matplotlib qt' instead\n",
    "import hyperspy.api as hs\n",
    "import lumispy as lum\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Loading files\n",
    "\n",
    "For saving analyses, HyperSpy has its own hdf5-based data format `.hspy`.\n",
    "\n",
    "HyperSpy supports a wide range of microscopy (and spectroscopy) related **data file types**:\n",
    "\n",
    "- Gatan `.dm3/.dm4`\n",
    "- Attolight `.sur` (For the moment does not parse `original_metadata` to `metadata`)\n",
    "- Delmic `.hdf5` (Still needs a parsing script at the moment that you can obtain from Delmic, will be contributed to RosettaSciIO)\n",
    "- Hamamatsu streak camera images in `.tif` format (Reading the itex `.img` format will be available in RosettaSciIO)\n",
    "- ...\n",
    "\n",
    "Further spectroscopy-related formats will be available soon (through **RosettaSciIO**):\n",
    "\n",
    "- Horiba JobinYvon XML `.xml`\n",
    "- TriVista XML `.tvf`\n",
    "- Renishaw wire format `.wdf`\n",
    "\n",
    "We 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:\n",
    "\n",
    "*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.*\n",
    "\n",
    "*You can also leave the path empty. A pop-up window will appear to select.*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cl1 = hs.load(\"../demo-files/01/01_demo.hspy\")\n",
    "cl2 = hs.load(\"../demo-files/load_from_GatanFiles/asymmetric-peak_map.dm4\")\n",
    "cl3 = hs.load(\"../demo-files/load_from_GatanFiles/quantum-emitters_map.dm4\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To see **parameters** that the function takes, in Jupyter, you can **display the docstring** by using a `?`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "hs.load?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data structure / Axes handling\n",
    "\n",
    "Each HyperSpy signal object has certain attributes that contain the relevant data about the axes, data and metadata.\n",
    "\n",
    "To understand the HyperSpy datastructure, lets have a look at the dataset `cl2` (Gatan file).\n",
    "\n",
    "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).\n",
    "\n",
    "The **signal class** provides certain specific routines, for example conversion to energy axis in the case of luminescence data.\n",
    "\n",
    "Our sample dataset has **two navigation dimensions** and **one signal (spectral) dimension**:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "cl2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "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**:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "cl2.axes_manager"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The **actual data** (signal intensity) is stored in a numpy array:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "cl2.data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Metadata\n",
    "\n",
    "For most supported file formats, the metadata is automatically parsed into **HyperSpy's metadata tree**.\n",
    "It contains information about the measurement, but potentially also about post-processing:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cl2.metadata"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In a separate tree, the **complete metadata from the vendor format** is read in (which follows different conventions depending on the format):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "cl2.original_metadata"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot / Explore\n",
    "\n",
    "We can easily plot and explore the hyperspectral data:\n",
    "\n",
    "*(In the following, we will use the preprocessed dataset `cl1`. The sample contains MethylammoniumLead Bromine (MAPbBr3) perovskite single crystals fabricated by Alice Dearle.)*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    },
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "cl1.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the **average CL spectrum** of the whole map:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cl1.mean().plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Indexing\n",
    "\n",
    "HyperSpy has a powerful numpy (Matlab) style indexing mechanism that distinguishes between navigation and signal axes:\n",
    "\n",
    "- `.inav[x1:x2,y1:y2]`\n",
    "- `.isig[s1:s2]`\n",
    "\n",
    "The index parameters can be either:\n",
    "- Integer: Index in the axis array\n",
    "- Float: Value in calibrated axis units\n",
    "\n",
    "For example, we can either plot a subset of the map in navigation space (selected using pixels as index):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "cl1.inav[2:23,0:20].plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Or, we can plot the mean spectrum in a certain spectral range (selected using wavelength units):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "cl1.isig[440.:600.].mean().plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Chromatic imaging:\n",
    "\n",
    "Indexing can also be used for color-filtered (chromatic) imaging.\n",
    "\n",
    "First, lets plot the **panchromatic image**:\n",
    "\n",
    "*(the object is transposed, so that we plot the intensity over navigation instead of signal dimensions)*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "cl1.T.mean().plot(cmap='viridis')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we can **plot the intensity in a selected spectral window** (color-filtered image) using indexing:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "cl1.isig[480.:550.].T.mean().plot(cmap='viridis')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "Alternatively, we can interactively select a spectral window (color-filtered image) using regions of interest:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    },
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "im = cl1.T\n",
    "im.plot()\n",
    "roi1 = hs.roi.SpanROI(left=455, right=485) #sets a digitalbandfilter\n",
    "im_roi1 = roi1.interactive(im, color=\"red\")\n",
    "im_roi1_mean = hs.interactive(im_roi1.mean,\n",
    "                          event=roi1.events.changed,\n",
    "                          recompute_out_event=None)\n",
    "im_roi1_mean.plot(cmap='viridis')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "## Axes types / Convert to energy scale\n",
    "\n",
    "Required versions: hyperspy>=1.7.0 and lumispy>=0.2.0\n",
    "\n",
    "HyperSpy has different types of axes:\n",
    "- The standard `UniformDataAxis` is defined through an `offset` and a `scale` (delta between pixels)\n",
    "- A `FunctionalDataAxis` is defined through a `UniformDataAxis` and a `function` to convert the values\n",
    "- A more general `DataAxis` is defined through an `axis` vector/array\n",
    "\n",
    "The *wavelength* scale of our sample object is a `UniformDataAxis`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cl1.axes_manager"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*LumiSpy* provides easy conversions of the signal axis to the **energy scale**:\n",
    "\n",
    "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`):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    },
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "cl1_eV = cl1.to_eV(inplace=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The axes manager now contains an `energy` axis, which is a non-uniform `DataAxis`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "cl1_eV.axes_manager"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To explore the data in the energy domain, we again plot the signal:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "cl1_eV.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Jacobian transformation\n",
    "\n",
    "To preserve the integrated intensity per spectral window, a *Jacobian* transformation has to be applied to the signal intensity:\n",
    "\n",
    "As we require $I(E)dE = I(\\lambda)d\\lambda$, the scale transformation $E=hc/\\lambda$ implies\n",
    "\n",
    "$$I(E) = I(\\lambda)\\frac{d\\lambda}{dE} = I(\\lambda)\\frac{d}{dE}\n",
    "\\frac{h c}{E} = - I(\\lambda) \\frac{h c}{E^2}$$\n",
    "\n",
    "(where the minus sign just reflects the different directions of integration in the wavelength and energy domains)\n",
    "\n",
    "This transformation is the default in LumiSpy, but can be deactivated by setting `jacobian=False`.\n",
    "\n",
    "To visualize the effect of the *Jacobian transformation*, we can plot a signal with constant intensity before and after the transformation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Create a model signal with linear intensity\n",
    "axis = {'offset': 300, 'scale': 4, 'units': 'nm', 'size': 101, 'name': 'Wavelength'}\n",
    "s = hs.signals.Signal1D(np.ones(101), axes=[axis])\n",
    "s.set_signal_type(\"Luminescence\")\n",
    "s2 = s.to_eV(inplace=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Some additional arrays to help with visualizing the spectral bins during plotting\n",
    "x = np.arange(9)*50+300\n",
    "x2 = lum.nm2eV(x)\n",
    "y2 = hs.signals.Signal1D(np.ones(9), axes=[{'offset': 300, 'scale': 50, 'size': 9,}])\n",
    "y2.set_signal_type(\"Luminescence\")\n",
    "y2.to_eV()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Plot comparative figures\n",
    "fig1 = plt.figure(figsize=(10,4))\n",
    "ax0 = plt.subplot(121)\n",
    "plt.ylim(0,1.3)\n",
    "plt.xlabel('Wavelength (nm)')\n",
    "ax0.plot(s.axes_manager[0].axis,s.data,color='orange')\n",
    "ax0.vlines(x,0,1,color='orange')\n",
    "ax0.fill_between(s.axes_manager[0].axis,0,s.data, facecolor='orange', alpha=0.3)\n",
    "ax1 = plt.subplot(122)\n",
    "plt.ylim(0,0.43)\n",
    "plt.xlabel('Energy (eV)')\n",
    "ax1.plot(s2.axes_manager[0].axis,s2.data,color='orange')\n",
    "ax1.vlines(x2,0,y2.data[::-1],color='orange')\n",
    "ax1.fill_between(s2.axes_manager[0].axis,0,s2.data, facecolor='orange', alpha=0.3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Model fitting\n",
    "\n",
    "We will introduce basic fitting functionality using our sample dataset `cl1`, for more details see the `Fitting_tutorial` in the [HyperSpy demos repository](https://github.com/hyperspy/hyperspy-demos).\n",
    "\n",
    "Note that HyperSpy has a range of [built-in functions](https://hyperspy.org/hyperspy-doc/current/user_guide/model.html#pre-defined-model-components) covering most needs that can be added as components to a model. However, it also has an intuitive mechanism to [define custom functions](https://hyperspy.org/hyperspy-doc/current/user_guide/model.html#define-components-from-a-mathematical-expression).\n",
    "\n",
    "First, we need to **initialize the model**:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    },
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "m = cl1_eV.create_model()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "jupyter": {
     "outputs_hidden": false
    },
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "source": [
    "**Check the components** of the model (should be empty, but for some types of signals, an automatic background component is added):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "m.components"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "Now we need to **create some components** and **add them to the model**.\n",
    "\n",
    "We will use a constant `Offset` and a Gaussian (defined through height and FWHM, hence `GaussianHF`):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    },
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "bkg = hs.model.components1D.Offset()\n",
    "g1 = hs.model.components1D.GaussianHF()\n",
    "m.extend([g1, bkg])\n",
    "m.components"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To see the parameters of our components and their default values, we can **print all parameter values**:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    },
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "m.print_current_values()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "jupyter": {
     "outputs_hidden": false
    },
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "source": [
    "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):\n",
    "\n",
    "*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.*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    },
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "# g1.estimate_parameters(cl1_eV,2.3,2.5)\n",
    "cl1.axes_manager.indices = (7,7)\n",
    "g1.centre.value = 2.4        # Gaussian centre\n",
    "g1.fwhm.value = 0.1      # Gaussian width\n",
    "g1.height.value = 5      # Gaussian height\n",
    "bkg.offset.value = 0.1   # background offset"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also **set boundaries** (`bmin` and `bmax`) for some of the parameters:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "g1.centre.bmax = g1.centre.value + 0.2\n",
    "g1.centre.bmin = g1.centre.value - 0.2\n",
    "g1.fwhm.bmin = 0.01"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "We can now **fit the model at the chosen position**, copy the result as starting value to all positions, and **plot** the result:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "m.fit(bounded=True)\n",
    "m.assign_current_values_to_all()\n",
    "m.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "Again, we can also **print the updated parameters**:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "m.print_current_values()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "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:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "m.multifit(bounded=True, show_progressbar=True)\n",
    "m.plot(plot_components=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To plot maps of the parameters of the Gaussian, we create signal objects from these datasets:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "m_centre = g1.centre.as_signal()\n",
    "m_centre.plot(cmap='bwr_r', centre_colormap=False) # Otherwise, it would be centred around `0` and we would see little difference between pixels\n",
    "m_intensity = g1.height.as_signal()\n",
    "m_intensity.plot(cmap='viridis')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Particle segmentation\n",
    "\n",
    "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:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    },
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "mask_treshold = m_intensity.data.mean()\n",
    "mask = m_intensity.data > mask_treshold #Returns a boolean matrix mask\n",
    "plt.figure()\n",
    "plt.imshow(mask)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now plot the previous graph of the centre-parametre, after applying the mask:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "(m_centre * mask).plot(cmap='bwr_r', vmin=2.3, centre_colormap=False) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Correction of spectral defects\n",
    "\n",
    "Working on the unprocessed dataset `cl2`, we can introduce some basic functions for artefact correction:\n",
    "\n",
    "### Remove background (interactive)\n",
    "\n",
    "HyperSpy has an interactive tool for **background removal** that supports various functions, let's start by removing a **simple offset**:\n",
    "1. Select a region to be used to determine the background: On the signal plot click, drag and release\n",
    "2. You can still move the region or its boundaries with the mouse and inspect the different spectra using the navigator to make sure the region is right\n",
    "3. Press `Apply`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "cl2.plot()\n",
    "cl2.remove_background(background_type=\"Offset\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Remove last pixels from the spectrum\n",
    "\n",
    "The signal beyond 800 nm goes to negative values, so lets remove the last three pixels from every spectrum:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "cl2 = cl2.isig[:-3]\n",
    "cl2.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Remove spikes (interactive)\n",
    "\n",
    "There is also a tool for interactive removal of cosmic rays (pixels with sharp spikes), see `Help` for instructions.\n",
    "\n",
    "In brief:\n",
    "- Inspect the derivative histogram\n",
    "- Set a sensible threshold to catch the outliers in the histogram\n",
    "- Iterate through `Find next` / `Remove spike` to continue for wrong identifications / remove identified spikes\n",
    "- `Close` when finished\n",
    "\n",
    "- *NOTE: The interactive version does not work well with inline plotting (e.g. on binder)*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "cl2.spikes_removal_tool(interactive=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "### Data smoothing\n",
    "\n",
    "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:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "cl2 = cl2.rebin(scale=[1,1,2])\n",
    "cl2.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "Additionally, HyperSpy provides three different functions for **data smoothing**:\n",
    "\n",
    "- `smooth_lowess` (lowess smoothing)\n",
    "- `smooth_savitzky_golay` (Savitzky Golay filter)\n",
    "- `smooth_tv` (total variation data smoothing)\n",
    "\n",
    "These functions can be run interactively to choose the right parameters, but the parameters can also be passed to the function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "cl2.plot()\n",
    "cl2.smooth_lowess(number_of_iterations = 2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "cl2.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Peak identification / Centroid / Peak width\n",
    "\n",
    "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.\n",
    "\n",
    "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.\n",
    "\n",
    "*As we have some side-peaks, we operate on a subrange of the wavelength axis.*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "peaks = cl2.isig[600.:].find_peaks1D_ohaver(maxpeakn=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The function **returns a structured array** that contains `position`, `height` and `width` for every pixel (potentially each for multiple peaks)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "peaks[0,0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "You can also just determine the **width of a peak** directly from the signal without fitting a model to the data.\n",
    "\n",
    "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."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "width = cl2.isig[600.:].estimate_peak_width(return_interval=True)\n",
    "width[0].plot(cmap='viridis')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "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.\n",
    "\n",
    "Required version: lumispy>=0.2.2\n",
    "\n",
    "*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.*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "com = cl2.isig[600.:].centroid()\n",
    "com.plot(cmap='viridis')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can **add markers** for the *FWHM interval* and the *centre of mass* to the signal object and plot them on the spectra:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "mrk = hs.plot.markers.vertical_line(com.data, color='black')\n",
    "mrkl = hs.plot.markers.vertical_line(width[1].data, color='grey')\n",
    "mrkr = hs.plot.markers.vertical_line(width[2].data, color='grey')\n",
    "cl2.add_marker([mrk,mrkl,mrkr], permanent=True)\n",
    "cl2.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Advanced fitting example\n",
    "\n",
    "Finally, to showcase the power of fitting with a programming language, lets fit the more complicated dataset `cl3`:\n",
    "\n",
    "1. Fit a spectrum with a skew normal component on a broad Gaussian as background.\n",
    "\n",
    "2. Then add a variable number of sharp Gaussians depending on the number of peaks found by the peak finder routine.\n",
    "\n",
    "*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](https://doi.org/10.1103/PhysRevApplied.10.044031) (2018). The analysis is simplified for instructive purposes.*\n",
    "\n",
    "*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](https://hyperspy.org/hyperspy-doc/current/user_guide/signal.html#iterating-external-functions-with-the-map-method).*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "cl3 = cl3.inav[10,8]\n",
    "m = cl3.create_model()\n",
    "g1 = hs.model.components1D.SkewNormal()\n",
    "g2 = hs.model.components1D.GaussianHF()\n",
    "m.extend((g1,g2))\n",
    "g1.x0.value=655\n",
    "g1.scale.value=50\n",
    "g1.shape.value=2\n",
    "g1.A.value=220\n",
    "g1.x0.bmin=635\n",
    "g1.x0.bmax=675\n",
    "g1.x0.bounded=True\n",
    "g2.centre.value=580\n",
    "g2.fwhm.value=180\n",
    "g2.height.value=20\n",
    "m.fit()\n",
    "m.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def getPeaks(S2):\n",
    "    S2=S2.rebin(new_shape=[334])\n",
    "    S2.smooth_savitzky_golay(window_length=15,polynomial_order=2)\n",
    "    peaks = S2.find_peaks1D_ohaver(amp_thresh=0.2*np.max(S2.data),maxpeakn=10)[0]\n",
    "    return peaks"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "ng=10\n",
    "S2 = cl3.deepcopy()\n",
    "g = list()\n",
    "for i in np.arange(ng):\n",
    "    g.append(hs.model.components1D.GaussianHF())\n",
    "m.extend(g)\n",
    "for i in np.arange(ng):\n",
    "    g[i].active_is_multidimensional = True\n",
    "# The following code would need to be in a loop to run for a whole spectral image\n",
    "# m.axes_manager.indices = cl3.axes_manager.indices\n",
    "peaks = getPeaks(S2.inav[cl3.axes_manager.indices])\n",
    "for i in np.arange(np.size(peaks)):\n",
    "    g[i].centre.value=peaks['position'][i]\n",
    "    g[i].centre.bmin=peaks['position'][i]-3\n",
    "    g[i].centre.bmax=peaks['position'][i]+3\n",
    "    g[i].centre.bounded=True\n",
    "    g[i].fwhm.value=5\n",
    "    g[i].fwhm.bmax=10\n",
    "    g[i].fwhm.bmin=1\n",
    "    g[i].fwhm.bounded=True\n",
    "    g[i].height.value=20\n",
    "    g[i].height.bmin=1\n",
    "    g[i].height.bounded=True\n",
    "if np.size(peaks)<ng:\n",
    "    for i in np.arange(np.size(peaks),ng):\n",
    "        g[i].active = False\n",
    "m.fit(bounded=True)\n",
    "m.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    },
    "tags": []
   },
   "source": [
    "## Now try with your own data!\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}