2_AdvancedExamples_HyperSpy-LumiSpy-eXSpy.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Advanced Examples for working with *HyperSpy*, *LumiSpy* and *eXSpy*\n",
    "\n",
    "Tutorial for the **eBEAM2024 school on nano-optics with free electrons**\n",
    "\n",
    "> Aussois, September 1-13, 2024\n",
    "\n",
    "**Table of Contents:**\n",
    "\n",
    "- [Creating a signal from a text file](#Creating-a-signal-from-a-text-file)\n",
    "- [Axes types / Convert to energy scale](#Axes-types-/-Convert-to-energy-scale)\n",
    "- [Model fitting](#Model-fitting)\n",
    "- [Particle segmentation](#Particle-segmentation)\n",
    "- [Extracting a line profile using a ROI](#Extracting-a-line-profile-using-a-ROI)\n",
    "- [Using the map function](#Using-the-map-function)\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": [
    "# 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",
    "%matplotlib widget\n",
    "\n",
    "import hyperspy.api as hs\n",
    "import lumispy as lum\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "# Plot multiple inline figures side-by-side horizontally \n",
    "hs.preferences.Plot.widget_plot_style = 'horizontal'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Loading files\n",
    "\n",
    "We will load two files that we will use during the demo. A preprocessed dataset saved in the `hspy` format and two maps in the `dm4` (Gatan) format:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "cl1 = hs.load(\"data/01_demo.hspy\")\n",
    "cl3 = hs.load(\"data/quantum-emitters_map.dm4\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Creating a signal from a text file\n",
    "\n",
    "Read tabular data from a txt file using `numpy.loadtxt` and create a `Signal1D` object. The signal axis and the EELS intensity values are\n",
    "given by the first and second columns, respectively.\n",
    "\n",
    "The tabular data are taken from https://eelsdb.eu/spectra/la2nio4-structure-of-k2nif4/"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Read tabular data from a text file:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x, y = np.loadtxt(\"data/La2NiO4_eels.txt\", unpack=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define the axes dictionary of the signal from the array `x` and then create the signal with data array `y`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "axes = [\n",
    "    # use values from first column to define non-uniform signal axis\n",
    "    dict(axis=x, name=\"Energy\", units=\"eV\"),\n",
    "    ]\n",
    "\n",
    "s = hs.signals.Signal1D(y, axes=axes)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Check the signal:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "s"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Check the `axes manager` (we have a non-uniform axis, as it is defined by an array and not by `offset` and `scale` parameters):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "s.axes_manager"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Convert the non-uniform axis to a uniform axis, because non-uniform axes do not support all functionalities of HyperSpy. In this case, the error introduced during conversion to uniform scale is negligeable (only related to the limited precission of the saved axis array)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "s.axes_manager.signal_axes[0].convert_to_uniform_axis()\n",
    "s.axes_manager"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Set title of the dataset and label for the data axis:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "s.metadata.set_item(\"General.title\", \"La2NiO4 EELS\")\n",
    "s.metadata.set_item(\"Signal.quantity\", \"Intensity (counts)\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To retrieve an item from the metadata, you can call the respective node in the metadata tree:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "s.metadata.General.title"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the dataset (the title and y-axis label are set based on the metadata we just added):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "s.plot()"
   ]
  },
  {
   "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",
    "We turn to dataset `cl1`, a hyperspectral CL map. The *wavelength* scale of our sample object is a `UniformDataAxis`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cl1.axes_manager"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cl1.plot()"
   ]
  },
  {
   "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 signal axis is now a *non-uniform axis*:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "cl1_eV.axes_manager"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This axis is defined through an axis vector (the signal axis being the last axis of our `axes_manager`):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cl1_eV.axes_manager[-1].axis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To explore the data in the energy domain, we plot the converted 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 been 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 constant intensity and a copy with Energy axis\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": [
    "# Define some additional arrays to visualize 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 some extended 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 again 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": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "m = cl1_eV.create_model()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    },
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "source": [
    "**Check the components** of the model:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "m.components"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "**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": {
    "collapsed": false,
    "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": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "# g1.estimate_parameters(cl1_eV,2.3,2.5)\n",
    "cl1_eV.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 current position**. Fitting for a single pixel allows to optimize the starting parameters and boundaries before we run the more time consuming fit on all pixels of the map. We set `bounded=True` to activate the boundaries we defined above.\n",
    "\n",
    "We then 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. We use the `cmap` parameter to set different color maps for the maps:"
   ]
  },
  {
   "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": [
    "The plots are matplotlib objects, so we can use the `savefig` function to save the last (currently selected) figure:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.savefig('mygraph.svg')"
   ]
  },
  {
   "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, to mask pixels where the intensity is below the mean value:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "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": [
    "## Extracting a line profile using a ROI\n",
    "\n",
    "1. Interactively extract a line profile (with a certain width) from the intensity map obtained from the fit using `roi.Line2DROI`.\n",
    "2. Use `plot.plot_spectra` to plot several line profiles on the same figure.\n",
    "3. Save a profile data as ``msa`` file.\n",
    "\n",
    "We take the signal containing the intensity map from our fit above:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m_intensity"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Intialize Line-ROI (region of interest) from position (4896, 3894) to position (6250, 4901) of width 10 in calibrated axes units:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "line_roi = hs.roi.Line2DROI(4896, 3894, 6250, 4901, 10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Extract data along the ROI as new signal `profile` by “slicing” the signal and plot the profile:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "profile = line_roi(m_intensity)\n",
    "profile.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So far, slicing of the signal is not interactive. If you want to modify the line along which the profile is extracted, you can plot the image and display the ROI interactively (creates a new signal object). You can even display the same ROI on a second image to make sure that a profile is well placed on both images (we use the `m_centre` map here)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m_intensity.plot()\n",
    "profile1 = line_roi.interactive(m_intensity, color='green')\n",
    "m_centre.plot()\n",
    "profile2 = line_roi.interactive(m_centre, color='green')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can then drag and drop the ends of the ROI to adjust the placement. You can also change the width of the extracted profile by moving the dashed lines marking this width.\n",
    "\n",
    "If you want to later update the ROI initialization above with the modified parameters, you can print `(start x, start y, end x, end y, width)` as tuple:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(tuple(line_roi))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can now directly access the data of the profile objects, e.g. to plot both profiles in a single plot (though this is of limited value here, as the centre and intensity have quite different scales). We use the `plot.plot_spectra` function, which takes a list of signals:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hs.plot.plot_spectra([profile1, profile2])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since the profile is a signal object, you can operate on it using any other functionality provided by hyperspy, e.g. to save a profile as `.msa` text file:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "profile1.save('extracted-line-profile.msa', format='XY', overwrite=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Using the map function\n",
    "\n",
    "To perform an operation on the data at each coordinate, HyperSpy provides the `map()` function.\n",
    "\n",
    "As a simple example, we will apply `np.max` to get the maximum intensity from each spectrum. However, the  [`map` function](https://hyperspy.org/hyperspy-doc/current/user_guide/signal.html#iterating-external-functions-with-the-map-method) can be used to apply any functions defined for individual spectra on a complete spectral image."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cl1max = cl1.map(np.max, inplace=False)\n",
    "cl1max.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Obviously, our dataset contains a cosmic spike (single bright pixel) that we can easily remove with the `spikes_removal_tool`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cl1.spikes_removal_tool(interactive=False, max_num_bins=6000)\n",
    "cl1max = cl1.map(np.max, inplace=False)\n",
    "cl1max.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "However, basic mathematic functions such as `max` are directly implemented in *HyperSpy* and we can get the same result without using `map`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cl1max = cl1.max(axis=-1)\n",
    "cl1max"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As the resulting signal has navigation, but no signal dimensions, we have to transpose it if we want to change the colormap, as the navigator plot does not support different colormaps:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cl1max.T.plot(cmap='viridis')"
   ]
  },
  {
   "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",
    "3. Subtract the two background components from the signal.\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.*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cl3 = cl3.inav[10,8]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1. Initiate and fit the background model with sensible starting parameters:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true,
    "tags": []
   },
   "outputs": [],
   "source": [
    "ma = cl3.create_model()\n",
    "ga1 = hs.model.components1D.SkewNormal()\n",
    "ga2 = hs.model.components1D.GaussianHF()\n",
    "ma.extend((ga1,ga2))\n",
    "ga1.x0.value=655\n",
    "ga1.scale.value=50\n",
    "ga1.shape.value=2\n",
    "ga1.A.value=220\n",
    "ga1.x0.bmin=635\n",
    "ga1.x0.bmax=675\n",
    "ga1.x0.bounded=True\n",
    "ga2.centre.value=580\n",
    "ga2.fwhm.value=180\n",
    "ga2.height.value=20\n",
    "ma.fit()\n",
    "ma.plot(plot_components=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "2. Identify sharp peaks in spectrum and for each add a Gaussian component to the model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Function to identify sharp peaks in spectrum using rebinning and smoothing due to the high noise level\n",
    "def getPeaks(s):\n",
    "    s=s.rebin(new_shape=[334])\n",
    "    s.smooth_savitzky_golay(window_length=15,polynomial_order=2)\n",
    "    peaks = s.find_peaks1D_ohaver(amp_thresh=0.2*np.max(s.data),maxpeakn=10)[0]\n",
    "    return peaks"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Add ten active (multidimensional) components to model\n",
    "ng=10\n",
    "g = list()\n",
    "for i in np.arange(ng):\n",
    "    g.append(hs.model.components1D.GaussianHF())\n",
    "ma.extend(g) \n",
    "for i in np.arange(ng):\n",
    "    g[i].active_is_multidimensional = True\n",
    "# Identify peaks and use to set fit parameters for the Gaussian components\n",
    "# The following code would need to be in a function that is looped or passed to map to run for a whole spectral image\n",
    "peaks = getPeaks(cl3.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",
    "# Deactivate additional, not-needed components\n",
    "if np.size(peaks)<ng:\n",
    "    for i in np.arange(np.size(peaks),ng):\n",
    "        g[i].active = False\n",
    "# Fit the amended model and plot\n",
    "ma.fit(bounded=True)\n",
    "ma.plot(plot_components=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "3. Subtract the two background components from the signal."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Deactivate components relating to the Gaussian peaks and plot improved background model\n",
    "for i in np.arange(np.size(peaks)):\n",
    "    g[i].active = False\n",
    "ma.plot(plot_components=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# convert model to signal that then can be subtracted from a copy of the data\n",
    "ms = ma.as_signal()\n",
    "sa = cl3.deepcopy()\n",
    "sa.data-=ms.data\n",
    "# Due to high noise level, we rebin the signal for plotting\n",
    "sa.rebin(scale=[3]).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": "hyperspy-stable",
   "language": "python",
   "name": "hyperspy-stable"
  },
  "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.12.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}