01_EnergyConversion_GaussianFitting.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 01: General `lumispy` workflow\n",
    "\n",
    "## Energy conversion and model fitting\n",
    "\n",
    "This notebook shows:\n",
    "\n",
    "- Plotting cathodoluminescence data in different ways (and interactively)\n",
    "- Converting wavelength signal axis to energy signal axis\n",
    "(**Required versions: hyperspy>=1.7.0 and lumispy>=0.2**)\n",
    "- Gaussian fitting"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Import packages:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "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 lumispy as lum\n",
    "import hyperspy.api as hs\n",
    "import os\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Loading the pre-processed CL-SEM file"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load the `.hspy` file, which has already been pre-processed (background subtracted, spike removed and spectral corrected). See microscope-specific tutorials on how to do that.\n",
    "\n",
    "*You can also leave the path empty. A pop-up window will appear to select the| file from the browser.*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Load file\n",
    "path = os.path.relpath(\"demo-files/01/01_demo.hspy\")\n",
    "cl_sem = hs.load(path, signal_type='CL_SEM')\n",
    "\n",
    "print(cl_sem)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "cl_sem.metadata"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plotting data\n",
    "\n",
    "Plot the hyperspectral data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    },
    "pycharm": {
     "name": "#%%\n"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "cl_sem.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the average CL spectrum:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "cl_sem.mean().plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Chromatic imaging:\n",
    "\n",
    "Plot the panchromatic image:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "cl_sem.T.mean().plot(cmap='viridis')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "Get the colour filtered panchromatic images.\n",
    "Select the energy region to plot as coloured image interactively."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    },
    "pycharm": {
     "name": "#%%\n"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "im = cl_sem.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": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    },
    "pycharm": {
     "name": "#%%\n"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "im_filtered = roi1(im).mean()\n",
    "\n",
    "roi_width = roi1.right - roi1.left\n",
    "roi_centre = roi1.left + 0.5* roi_width\n",
    "\n",
    "im_filtered.metadata.General.title = \"Colour filtered image of {:.0f} $\\pm$ {:.0f} nm\".format(roi_centre, roi_width)\n",
    "im_filtered.plot(cmap='viridis')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "## Axes types and wavelength to energy conversion\n",
    "\n",
    "Required versions: hyperspy>=1.7.0 and lumispy>=0.2\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": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "cl_sem.axes_manager"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "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"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "cl_sem_eV = cl_sem.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": [
    "cl_sem_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": [
    "cl_sem_eV.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "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",
    "import numpy as np\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": [
    "## Fitting Gaussian\n",
    "\n",
    "We will introduce basic fitting functionality, for more details see the `Fitting_tutorial` in the [HyperSpy demos repository](https://github.com/hyperspy/hyperspy-demos).\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 = cl_sem_eV.create_model()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "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": {},
   "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": {
    "tags": []
   },
   "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": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "m.print_current_values()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "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):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "cl_sem_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": {},
   "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()\n",
    "m.assign_current_values_to_all()\n",
    "m.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "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": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "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": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    },
    "pycharm": {
     "name": "#%%\n"
    },
    "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": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    },
    "pycharm": {
     "name": "#%%\n"
    },
    "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"
    },
    "tags": []
   },
   "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": {
    "tags": []
   },
   "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": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "## END\n"
   ]
  }
 ],
 "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
}