Merge pull request #113 from jlaehne/doc-energy

Documentation: axis conversion

README.md

@@ -18,17 +18,16 @@ aimed at helping with the analysis of luminescence spectroscopy data
 (cathodoluminescence, photoluminescence, electroluminescence, Raman, SNOM).
 
 If analysis using LumiSpy forms a part of published work, please consider 
-recognising the code development by citing the project using the [Zenodo-DOI](https://doi.org/10.5281/zenodo.4640445).
+recognising the code development by citing the project using the
+[Zenodo-DOI](https://doi.org/10.5281/zenodo.4640445).
 
-Go to the documentation for instructions on how to install LumiSpy and get started with analysis: [Read the docs](https://lumispy.readthedocs.io).
+Go to the documentation for instructions on how to install LumiSpy and start an
+analysis: [Read the docs](https://lumispy.readthedocs.io).
 
-[Tutorials and example workflows](https://github.com/lumispy/lumispy-demos)
+[Tutorials and exemplary workflows](https://github.com/lumispy/lumispy-demos)
 have been curated as a series of Jupyter notebooks that you can work through 
-and modify to perform many common analyses. Simply:
-
-1. Download the `lumispy_demos` repository in your desired folder
-2. Load LumiSpy (see [installation guide](https://lumispy.readthedocs.io/en/latest/user_guide/installation.html))
-3. In Jupyter lab, navigate to the folder and start using the notebook
+and modify to perform many common analyses. These can be either downloaded and
+run locally or tried out using interactive online sessions.
 
 Everyone is welcome to contribute. Please read our
 [contributing guidelines](https://github.com/LumiSpy/lumispy/blob/main/CONTRIBUTING.rst) and get started!

doc/source/conf.py

@@ -42,6 +42,16 @@
 }
 intersphinx_disabled_domains = ["std"]
 
+# imgmath: Sphinx allows use of LaTeX in the html documentation, but not directly. It is first rendered to an image.
+# You can add here whatever preamble you are used to adding to your LaTeX document.
+imgmath_latex_preamble = r"""
+    \usepackage{xcolor}
+    \definecolor{mathcolor}{rgb}{0.8,0.3,0.1}
+    \everymath{\color{mathcolor}}
+    %\everydisplay{\color{mathcolor}}
+"""
+
+
 templates_path = ["_templates"]
 
 # -- Options for HTML output

doc/source/index.rst

@@ -33,6 +33,10 @@ aimed at helping with the analysis of luminescence spectroscopy data
 Check out the :ref:`installation-label` section for further information, including
 how to start using this project.
 
+Complementing this documentation, the `LumiSpy Demos <https://github.com/LumiSpy/lumispy-demos>`_
+repository contains curated Jupyter notebooks to provide tutorials and exemplary
+workflows.
+
 .. note::
 
    This project is under active development. Everyone is welcome to contribute. Please read our (see :ref:`contributing_label`) guidelines and get started!
@@ -56,6 +60,12 @@ Contents
 
    api/modules.rst
 
+.. toctree::
+   :maxdepth: 1
+   :caption: Tutorials
+
+   Demo notebooks <https://github.com/LumiSpy/lumispy-demos>
+
 .. toctree::
    :maxdepth: 1
    :caption: Release Notes

doc/source/user_guide/fitting_luminescence.rst

@@ -1,11 +1,18 @@
-.. _fitting_luminescence:
+.. _fitting_luminescence-label:
 
 Fitting luminescence data
 *************************
 
-LumiSpy is compatible with model fitting as Hyperspy. It can fit both in linear and non-linear axes.
+LumiSpy is compatible with HyperSpy model fitting. It can fit both in linear and non-linear axes.
 
 TODO: Link to Hyperspy guide
 TODO: Note on advantages of fitting Gaussians in the ``eV`` axis.
 TODO: Show how to extract the *modeled signal* with all/one component.
 
+.. _fitting_variance-label:
+
+Signal variance (noise)
+=======================
+
+TODO: Documentation on variance handling in the context of fitting
+in particular using ``s.estimate_poissonian_noise_variance()``

doc/source/user_guide/signal_axis.rst

@@ -1,14 +1,158 @@
-.. _signal_axis:
+.. _signal_axis-label:
 
 Non-uniform signal axes
 ***********************
 
-LumiSpy enables the use of non-linear axis (e.g. ``eV``) in an easy way.
+LumiSpy facilitates the use of `non-uniform axes 
+<https://hyperspy.org/hyperspy-doc/current/user_guide/axes.html#non-uniform-data-axis>`_,
+where the points of the axis vector are not uniformly spaced. This situation
+occurs in particular when converting a wavelength scale to energy (eV) or
+wavenumbers (e.g. for Raman shifts).
 
-The function :py:func:`~.signals.luminescence_spectrum.LumiSpectrum.to_eV` simplifies the definition
+The conversion of the signal axis can be performed using the functions 
+:py:meth:`~.signals.luminescence_spectrum.LumiSpectrum.to_eV`,
+:py:meth:`~.signals.luminescence_spectrum.LumiSpectrum.to_invcm` and
+:py:meth:`~.signals.luminescence_spectrum.LumiSpectrum.to_raman_shift`
+(alias for :py:meth:`~.signals.luminescence_spectrum.LumiSpectrum.to_invcm_relative`).
+If the unit of the signal axis is set, the functions can handle wavelengths in
+either nm or µm.
 
-TODO: Explain the Jacobian transform.
-TODO: Show how the transform is done ``s.to_eV()``.
+Accepted parameters are ``inplace=True/False`` (default is True), which
+determines whether the current signal object is modified or a new one is
+created, and ``jacobian=True/False`` (default is True, see
+:ref:`jacobian-label`).
 
-TODO: Note on signal noise.
+.. Note::
 
+    The non-uniform axis functionality will be available from HyperSpy v.1.7.
+    If this version is not yet available, you need to use the `development
+    branch <https://github.com/hyperspy/hyperspy>`_.
+
+
+.. _energy_axis-label:
+
+The energy axis
+===============
+
+The transformation from wavelength :math:`\lambda` to energy :math:`E` is
+defined as :math:`E = h c/ \lambda`. Taking into account the permittivity of
+air and doing a conversion from nm to eV, we get:
+
+.. math::
+
+    E = \frac{10^9 h c}{e \epsilon_r \lambda},
+
+where :math:`h` is the Planck constant, :math:`c` is the speed of light,
+:math:`e` is the elementary charge and :math:`\epsilon_r` is the relative
+permittivity of air.
+
+.. code-block:: python
+
+    >>> s2 = s.to_eV(inplace=False)
+    >>> s.to_eV()
+
+.. Note::
+
+    The relative permittivity of air :math:`\epsilon_r` is wavelength
+    dependent. This dependence is taken into account by LumiSpy based on the
+    analytical formula given by [Peck]_ valid from 185-1700 nm
+    (outside of this range, the permittivity values at the edges of the range
+    are used and a warning is raised).
+
+
+.. _wavenumber_axis-label:
+
+The wavenumber axis/Raman shifts
+================================
+
+The transformation from wavelength :math:`\lambda` to wavenumber
+:math:`\tilde{\nu}` (spatial frequency of the wave) is defined as
+:math:`\tilde{\nu} = 1/ \lambda`. The wavenumber is usually given in units of
+:math:`\mathrm{cm}^{-1}`.
+
+When converting a signal to Raman shift, i.e. the shift in wavenumbers from
+the exciting laser wavelength, the laser wavelength has to be passed to the function using the parameter
+``laser`` using the same units as for the original axis (e.g. 325 for nm or
+0.325 for µm) unless it is contained in the :ref:`metadata_structure` under
+``Acquisition_instrument.Laser.wavelength``.
+
+TODO: Automatically read laser wavelength from metadata if given there.
+
+.. code-block:: python
+
+    >>> s2 = s.to_invcm(inplace=False)
+    >>> s.to_invcm()
+    >>> s2 = s.to_raman_shift(inplace=False, laser=325)
+    >>> s.to_raman_shift(laser=325)
+
+
+.. _jacobian-label:
+
+Jacobian transformation
+=======================
+
+When transforming the signal axis, the signal intensity is automatically
+rescaled (Jacobian transformation), unless the ``jacobian=False`` option is
+given. Only converting the signal axis, and leaving the signal intensity
+unchanged, implies that the integral of the signal over the same interval would
+lead to different results depending on the quantity on the axis (see e.g.
+[Mooney]_.).
+
+For the energy axis as example, if we require :math:`I(E)dE = I(\lambda)d\lambda`,
+then :math:`E=hc/\lambda` implies
+
+.. math ::
+
+    I(E) = I(\lambda)\frac{d\lambda}{dE} = I(\lambda)\frac{d}{dE}
+    \frac{h c}{\lambda} = - I(\lambda) \frac{h c}{E^2}
+
+The minus sign just reflects the different directions of integration in
+the wavelength and energy domains. The same argument holds for the conversion
+from wavelength to wavenumber (just without the additional prefactors in the
+equation). The renormalization in LumiSpy is defined such that the intensity is
+converted from counts/nm (or counts/µm) to counts/meV. The following
+figure illustrates the effect of the Jacobian transformation:
+
+.. image:: images/jacobian.png
+  :width: 700
+  :alt: Illustration of the Jacobian transformation from wavelength (nm) to energy (eV).
+
+
+.. _jacobian_variance-label:
+
+Transformation of the variance
+------------------------------
+
+Scaling the signal intensities implies that also the stored variance of the
+signal needs to be scaled accordingly. According to :math:`Var(aX) = a^2Var(X)`,
+the variance has to be multiplied with the square of the Jacobian. This squared
+renormalization is automatically performed by LumiSpy if ``jacobian=True``.
+In particular, homoscedastic (constant) noise will consequently become
+heteroscedastic (changing as a function of the signal axis vector). Therefore,
+if the ``metadata.Signal.Noise_properties.variance`` attribute is a constant,
+it is converted into a :external:py:class:`hyperspy.signal.BaseSignal` object
+before the transformation.
+
+See :ref:`fitting_variance-label` for more general information on data variance
+in the context of model fitting and the HyperSpy documentation on `setting
+the noise properties
+<https://hyperspy.org/hyperspy-doc/current/user_guide/signal.html?highlight=variance_linear_model#setting-the-noise-properties>`_.
+
+.. Note::
+
+    If the Jacobian transformation is performed, the values of
+    ``metadata.Signal.Noise_properties.Variance_linear_model`` are reset to
+    their default values (``gain_factor=1``, ``gain_offset=0`` and ``correlation_factor=1``).
+    Should these values deviate from the defaults, make sure to run
+    :external:py:meth:`hyperspy.signal.BaseSignal.estimate_poissonian_noise_variance`
+    prior to the transformation.
+
+
+.. rubric:: References
+
+.. [Peck] E.R. Peck and K. Reeder, J. Opt. Soc. Am. **62**, 958
+    (1972). `doi:10.1364/JOSA.62.000958 <https://doi.org/10.1364/JOSA.62.000958>`_
+
+.. [Mooney] J. Mooney and P. Kambhampati, The Journal of
+    Physical Chemistry Letters **4**, 3316 (2013).
+    `doi:10.1021/jz401508t <https://doi.org/10.1021/jz401508t>`_

lumispy/signals/luminescence_spectrum.py

@@ -103,9 +103,9 @@ def to_eV(self, inplace=True, jacobian=True):
         units of nm unless the axis units are specifically set to µm.
 
         The intensity is converted from counts/nm (counts/µm) to counts/meV by
-        doing a Jacobian transformation, see e.g. Wang and Townsend, J. Lumin.
-        142, 202 (2013), doi:10.1016/j.jlumin.2013.03.052, which ensures that
-        integrated signals are correct also in the energy domain. If the
+        doing a Jacobian transformation, see e.g. Mooney and Kambhampati, J.
+        Phys. Chem. Lett. 4, 3316 (2013), doi:10.1021/jz401508t, which ensures
+        that integrated signals are correct also in the energy domain. If the
         variance of the signal is known, i.e.
         `metadata.Signal.Noise_properties.variance` is a signal representing
         the variance, a squared renormalization of the variance is performed.
@@ -262,8 +262,8 @@ def to_eV(self, inplace=True, jacobian=True):
 
     TO_INVCM_DOCSTRING = """
         The intensity is converted from counts/nm (counts/µm) to counts/cm^-1
-        by doing a Jacobian transformation, see e.g. Wang and Townsend,
-        J. Lumin. 142, 202 (2013), doi:10.1016/j.jlumin.2013.03.052, which
+        by doing a Jacobian transformation, see e.g. Mooney and Kambhampati,
+        J. Phys. Chem. Lett. 4, 3316 (2013), doi:10.1021/jz401508t, which
         ensures that integrated signals are correct also in the wavenumber 
         domain. If the variance of the signal is known, i.e.
         `metadata.Signal.Noise_properties.variance` is a signal representing the
@@ -439,12 +439,16 @@ def to_invcm(self, inplace=True, jacobian=True):
         branch of HyperSpy.    
     """
 
-    def to_invcm_relative(self, laser, inplace=True, jacobian=True):
+    def to_invcm_relative(self, laser=None, inplace=True, jacobian=True):
         """Converts signal axis of 1D signal to non-linear wavenumber axis
         (cm^-1) relative to the exciting laser wavelength (Stokes/Anti-Stokes
         shift). Assumes wavelength in units of nm unless the axis units are
         specifically set to µm.
         %s
+        laser: float or None
+            Laser wavelength in the same units as the signal axis. If None
+            (default), checks if it is stored in
+            `metadata.Acquisition_instrument.Laser.wavelength`.
         %s
         """
 
@@ -455,6 +459,25 @@ def to_invcm_relative(self, laser, inplace=True, jacobian=True):
                 " if the RELEASE_next_minor branch of HyperSpy is used."
             )
 
+        # check if laser wavelength is available
+        if laser == None:
+            if not self.metadata.has_item("Acquisition_instrument.Laser.wavelength"):
+                raise AttributeError(
+                    "Laser wavelength is neither given in the metadata nor passed"
+                    " to the function."
+                )
+            else:
+                laser = self.metadata.get_item(
+                    "Acquisition_instrument.Laser.wavelength"
+                )
+        # check if laser units make sense in respect to signal units
+        if (self.axes_manager.signal_axes[0].units == "µm" and laser > 10) or (
+            self.axes_manager.signal_axes[0].units == "nm" and laser < 100
+        ):
+            raise AttributeError(
+                "Laser wavelength units do not seem to match the signal units."
+            )
+
         invcmaxis, factor = axis2invcm(self.axes_manager.signal_axes[0])
 
         # convert to relative wavenumber scale
@@ -496,6 +519,9 @@ def to_invcm_relative(self, laser, inplace=True, jacobian=True):
 
     to_invcm_relative.__doc__ %= (TO_INVCM_DOCSTRING, TO_INVCMREL_EXAMPLE)
 
+    # Alias Method Name
+    to_raman_shift = to_invcm_relative
+
     def remove_background_from_file(self, background=None, inplace=False, **kwargs):
         """Subtract the background to the signal in all navigation axes.If no
         background file is passed as argument, the `remove_background()` from