Repository: kmnhan/erlabpy (language: python)

Issue:
Add curvature and minimum gradient processing



Discussion / hints:

Fix patch:
diff --git a/.github/workflows/pr.yml b/.github/workflows/pr.yml
index ab43a80a..25557e91 100644
--- a/.github/workflows/pr.yml
+++ b/.github/workflows/pr.yml
@@ -28,8 +28,8 @@ jobs:
         run: |
           sudo apt update && sudo apt install -y libegl1-mesa-dev
           python -m pip install --upgrade pip
-          python -m pip install -v . pytest pytest-xdist
           python -m pip install -r requirements.txt
+          python -m pip install -v . pytest pytest-xdist
 
 
       - name: Test with pytest
diff --git a/.github/workflows/release.yml b/.github/workflows/release.yml
index 5f81fb6e..71ac74aa 100644
--- a/.github/workflows/release.yml
+++ b/.github/workflows/release.yml
@@ -26,8 +26,8 @@ jobs:
         run: |
           sudo apt update && sudo apt install -y libegl1-mesa-dev
           python -m pip install --upgrade pip
-          python -m pip install -v . pytest pytest-xdist
           python -m pip install -r requirements.txt
+          python -m pip install -v . pytest pytest-xdist
 
 
       - name: Test with pytest
diff --git a/docs/environment.yml b/docs/environment.yml
index f677e565..56674744 100644
--- a/docs/environment.yml
+++ b/docs/environment.yml
@@ -5,6 +5,7 @@ dependencies:
   - python>=3.11
   - pip>=19.3
   - ipykernel
+  - igor2>=0.5.6
   - numpy>=1.22.0
   - matplotlib>=3.8.0
   - hvplot
@@ -34,7 +35,5 @@ dependencies:
   - pip:
     - pyqt6>=6.2.2
     - csaps
-    - igor2 @ git+https://github.com/kmnhan/igor2@master
-    - arpes @ git+https://github.com/kmnhan/arpes@master
     - -e git+https://github.com/kmnhan/erlabpy.git#egg=erlab
     - sphinx-qt-documentation
diff --git a/docs/requirements.txt b/docs/requirements.txt
index 25ee242a..db7b14e9 100644
--- a/docs/requirements.txt
+++ b/docs/requirements.txt
@@ -30,6 +30,5 @@ qtawesome
 superqt
 csaps
 qtpy
-igor2 @ git+https://github.com/kmnhan/igor2@master
-arpes @ git+https://github.com/kmnhan/arpes@master
+igor2>=0.5.6
 -e git+https://github.com/kmnhan/erlabpy.git#egg=erlab
diff --git a/docs/source/conf.py b/docs/source/conf.py
index 6e55e2eb..4373965c 100644
--- a/docs/source/conf.py
+++ b/docs/source/conf.py
@@ -68,7 +68,7 @@
 autosummary_generate_overwrite = True
 
 autodoc_inherit_docstrings = False
-autodoc_typehints = "both"
+autodoc_typehints = "description"
 autodoc_typehints_format = "short"
 autodoc_class_signature = "mixed"
 autodoc_member_order = "bysource"
@@ -78,7 +78,7 @@
     # "exclude-members":("sigDataChanged",),
     "show-inheritance": False,
 }
-autodoc_typehints_description_target = "all"
+autodoc_typehints_description_target = "documented"
 
 # autodoc_type_aliases = {
 # "numpy.float64": "float",
@@ -127,7 +127,6 @@
     "xarray": ("https://docs.xarray.dev/en/stable/", None),
     "pandas": ("https://pandas.pydata.org/docs/", None),
     "pyqtgraph": ("https://pyqtgraph.readthedocs.io/en/latest/", None),
-    "arpes": ("https://arpes.readthedocs.io/en/latest/", None),
     "csaps": ("https://csaps.readthedocs.io/en/latest/", None),
     "iminuit": ("https://scikit-hep.org/iminuit/", None),
 }
@@ -150,22 +149,63 @@
 ]
 
 
-class ApsStyle(pybtex.style.formatting.unsrt.Style):
+class APSStyle(pybtex.style.formatting.unsrt.Style):
     """
     APS style for BibTeX formatting, adapted from the conf.py file of the `mitiq
     library<https://github.com/unitaryfund/mitiq>`_.
     """
 
+    def __init__(self, *args, **kwargs):
+        super().__init__(*args, **kwargs)
+        self.abbreviate_names = True
+
     def format_title(self, e, which_field, as_sentence=True):
         formatted_title = pybtex.style.template.field(
             which_field, apply_func=lambda text: text.capitalize()
         )
+
         formatted_title = pybtex.style.template.tag("em")[formatted_title]
         if as_sentence:
             return pybtex.style.template.sentence[formatted_title]
         else:
             return formatted_title
 
+    def format_editor(self, e, as_sentence=True):
+        editors = self.format_names("editor", as_sentence=False)
+        if "editor" not in e.persons:
+            return editors
+        result = pybtex.style.template.join(sep=" ")["edited by", editors]
+        if as_sentence:
+            return pybtex.style.template.sentence[result]
+        else:
+            return result
+
+    def format_address_organization_publisher_date(self, e, include_organization=True):
+        if include_organization:
+            organization = pybtex.style.template.optional_field("organization")
+        else:
+            organization = None
+        return pybtex.style.template.first_of[
+            pybtex.style.template.optional[
+                pybtex.style.template.join(sep=", ")[
+                    pybtex.style.template.sentence(add_period=False, sep=", ")[
+                        organization, pybtex.style.template.optional_field("publisher")
+                    ],
+                    pybtex.style.template.join(sep=", ")[
+                        pybtex.style.template.sentence(add_period=False)[
+                            pybtex.style.template.optional_field("address")
+                        ],
+                        pybtex.style.template.field("year"),
+                    ],
+                ]
+            ],
+            pybtex.style.template.join(sep=", ")[
+                organization,
+                pybtex.style.template.optional_field("publisher"),
+                pybtex.style.template.field("year"),
+            ],
+        ]
+
     def get_article_template(self, e):
         volume_and_pages = pybtex.style.template.first_of[
             # volume and pages
@@ -252,8 +292,38 @@ def get_incollection_template(self, e):
         ]
         return template
 
+    def get_inproceedings_template(self, e):
+        template = pybtex.style.formatting.toplevel[
+            pybtex.style.template.sentence[self.format_names("author")],
+            self.format_title(e, "title"),
+            pybtex.style.template.words[
+                "in",
+                pybtex.style.template.sentence[
+                    self.format_btitle(e, "booktitle", as_sentence=False),
+                    # self.format_volume_and_series(e, as_sentence=False),
+                    # self.format_chapter_and_pages(e),
+                    pybtex.style.template.field("year"),
+                    pybtex.style.template.sentence(add_period=False, sep=" ")[
+                        pybtex.style.template.optional[
+                            self.format_editor(e, as_sentence=False)
+                        ],
+                        pybtex.style.template.words(sep="")[
+                            "(", self.format_address_organization_publisher_date(e), ")"
+                        ],
+                    ],
+                    pybtex.style.template.sentence(add_period=True, sep=". ")[
+                        "p",
+                        pybtex.style.formatting.unsrt.pages,
+                        pybtex.style.template.optional_field("note"),
+                    ],
+                ],
+            ],
+            self.format_web_refs(e),
+        ]
+        return template
+
 
-pybtex.plugin.register_plugin("pybtex.style.formatting", "apsstyle", ApsStyle)
+pybtex.plugin.register_plugin("pybtex.style.formatting", "apsstyle", APSStyle)
 
 bibtex_bibfiles = ["refs.bib"]
 bibtex_default_style = "apsstyle"
@@ -296,8 +366,13 @@ def get_incollection_template(self, e):
 # https://www.sphinx-doc.org/en/master/usage/configuration.html#latex-options
 
 latex_engine = "lualatex"
-latex_table_style = ["booktabs"]
+latex_show_pagerefs = True
+latex_show_urls = "footnote"
+latex_table_style = ["booktabs", "colorrows"]
 latex_elements = {
     "fontpkg": r"""\usepackage{fontspec,unicode-math}
-""",
+\setmainfont{DejaVu Serif}
+\setsansfont{DejaVu Sans}
+\setmonofont{DejaVu Sans Mono}
+"""
 }
diff --git a/docs/source/notebooks/kconv.ipynb b/docs/source/notebooks/kconv.ipynb
index bb8a2a56..652164e4 100644
--- a/docs/source/notebooks/kconv.ipynb
+++ b/docs/source/notebooks/kconv.ipynb
@@ -21,7 +21,7 @@
     "Nomenclature\n",
     "------------\n",
     "\n",
-    "Momentum conversion in ERLabPy follows the nomenclature from :cite:t:`Ishida2018`. All experimental geometry are classified into 4 types. Definition of angles differ for each geometry.\n",
+    "Momentum conversion in ERLabPy follows the nomenclature from :cite:t:`ishida2018kconv`. All experimental geometry are classified into 4 types. Definition of angles differ for each geometry.\n",
     "\n",
     "For instance, image a typical Type 1 setup with a vertical slit that acquires maps by rotating about the `z` axis in the lab frame. In this case, the polar angle (rotation about `z`) is :math:`β`, and the tilt angle is :math:`ξ`.\n",
     "\n",
diff --git a/docs/source/refs.bib b/docs/source/refs.bib
index 237dc048..444ab586 100644
--- a/docs/source/refs.bib
+++ b/docs/source/refs.bib
@@ -1,17 +1,16 @@
-@article{Ishida2018,
-  author   = {Ishida, Y. and Shin, S.},
-  title    = {Functions to map photoelectron distributions in a variety of setups in angle-resolved photoemission spectroscopy},
-  journal  = {Rev. Sci. Instrum.},
-  volume   = {89},
-  number   = {4},
-  pages    = {043903},
-  keywords = {photoelectron microscopy},
-  doi      = {10.1063/1.5007226},
-  year     = {2018},
-  type     = {Journal Article}
+@article{ishida2018kconv,
+  author  = {Ishida, Y. and Shin, S.},
+  title   = {Functions to map photoelectron distributions in a variety of setups in angle-resolved photoemission spectroscopy},
+  journal = {Rev. Sci. Instrum.},
+  volume  = {89},
+  number  = {4},
+  pages   = {043903},
+  doi     = {10.1063/1.5007226},
+  year    = {2018},
+  type    = {Journal Article}
 }
 
-@article{Seah1979,
+@article{seah1979imfp,
   author  = {Seah, M. P. and Dench, W. A.},
   title   = {Quantitative electron spectroscopy of surfaces: A standard data base for electron inelastic mean free paths in solids},
   journal = {Surf. Interface Anal.},
@@ -24,13 +23,14 @@ @article{Seah1979
   type    = {Journal Article}
 }
 
-@inproceedings{Schirra2008,
+@inproceedings{schirra2008polygon,
   author    = {Schirra, Stefan},
   title     = {How Reliable Are Practical Point-in-Polygon Strategies?},
   booktitle = {Algorithms - ESA 2008},
   editor    = {Halperin, Dan and Mehlhorn, Kurt},
   series    = {Lect. Note. Comput. Sci},
   publisher = {Springer Berlin Heidelberg},
+  address   = {Berlin, Heidelberg},
   volume    = {5193},
   pages     = {744-755},
   isbn      = {978-3-540-87744-8},
@@ -38,3 +38,38 @@ @inproceedings{Schirra2008
   year      = {2008},
   type      = {Conference Proceedings}
 }
+
+@article{zhang2011curvature,
+  author  = {Zhang, P. and Richard, P. and Qian, T. and Xu, Y.-M. and Dai, X. and Ding, H.},
+  title   = {A precise method for visualizing dispersive features in image plots},
+  journal = {Rev. Sci. Instrum.},
+  volume  = {82},
+  number  = {4},
+  pages   = {043712},
+  doi     = {10.1063/1.3585113},
+  year    = {2011},
+  type    = {Journal Article}
+}
+
+@article{he2017mingrad,
+  author  = {He, Y. and Wang, Y. and Shen, Z.-X.},
+  title   = {Visualizing dispersive features in 2D image via minimum gradient method},
+  journal = {Rev. Sci. Instrum.},
+  volume  = {88},
+  number  = {7},
+  pages   = {073903},
+  doi     = {10.1063/1.4993919},
+  year    = {2017},
+  type    = {Journal Article}
+}
+
+@article{hoyer2017xarray,
+  title     = {xarray: {N-D} labeled arrays and datasets in {Python}},
+  author    = {Hoyer, S. and J. Hamman},
+  journal   = {Journal of Open Research Software},
+  volume    = {5},
+  number    = {1},
+  year      = {2017},
+  publisher = {Ubiquity Press},
+  doi       = {10.5334/jors.148}
+}
diff --git a/docs/source/userguide.rst b/docs/source/userguide.rst
index fa404639..e5a4ddbe 100644
--- a/docs/source/userguide.rst
+++ b/docs/source/userguide.rst
@@ -5,9 +5,11 @@ User Guide
 .. admonition:: Work in Progress
    :class: warning
 
-   The user guide is incomplete. For the full list of packages and modules provided by ERLabPy, see :doc:`reference`.
+   The user guide is incomplete. For the full list of packages and modules
+   provided by ERLabPy, see :doc:`reference`.
 
-This section contains some examples for getting started with ARPES data analysis and visualization.
+This section contains some examples for getting started with ARPES data analysis
+and visualization.
 
 
 Introduction
@@ -15,9 +17,14 @@ Introduction
 
 The following documentation assumes basic python programming experience.
 
-If you are new to scientific programming in python altogether, `Scientific Python Lectures <https://github.com/jrjohansson/scientific-python-lectures>`_ is a great resource to get started.
+If you are new to scientific programming in python altogether, `Scientific
+Python Lectures <https://github.com/jrjohansson/scientific-python-lectures>`_ is
+a great resource to get started.
 
-Data in ERLabPy are mostly represented by :mod:`xarray` objects. The user guide of the `xarray documentation <https://docs.xarray.dev/en/stable/index.html>`_ provides an excellent overview on data manipulation.
+Data in ERLabPy are mostly represented by :mod:`xarray` objects
+:cite:p:`hoyer2017xarray`. The `xarray user guide
+<https://docs.xarray.dev/en/stable/index.html>`_ and the `xarray tutorial
+<https://tutorial.xarray.dev/>`_ are great resources to get started with xarray.
 
 .. toctree::
    notebooks/io
@@ -28,5 +35,8 @@ Data in ERLabPy are mostly represented by :mod:`xarray` objects. The user guide
 Further Reading
 ===============
 
-- `Lectures on scientific computing with Python <https://github.com/jrjohansson/scientific-python-lectures>`_
-- `The beginner's guide to numpy <https://numpy.org/doc/stable/user/absolute_beginners.html>`_
+- `Lectures on scientific computing with Python
+  <https://github.com/jrjohansson/scientific-python-lectures>`_
+- `The beginner's guide to numpy
+  <https://numpy.org/doc/stable/user/absolute_beginners.html>`_
+- `Xarray tutorial <https://tutorial.xarray.dev/>`_
diff --git a/environment.yml b/environment.yml
index 7d17a044..1f5f3636 100644
--- a/environment.yml
+++ b/environment.yml
@@ -2,51 +2,28 @@ name: erlab
 channels:
   - conda-forge
 dependencies:
-  - python>=3.11
-  - ipython
-  - pip>=19.3
-  - setuptools
-  - h5netcdf
-  - varname
-  - numpy>=1.22.0
-  - scipy>=1.8.0
-  - matplotlib>=3.7.0
-  - h5py>=3.2.1
-  - netCDF4>=1.5.0
-  - astropy
-  - superqt>=0.3.1
-  - dill
-  - llvm-openmp>=13.0.0
-  - llvmlite
-  - asteval>=0.9.22
-  - xarray>=0.16.1
-  - numba>=0.56.2
-  - ase>=3.17.0
-  - bokeh>=2.0.0
-  - bottleneck>=1.3.2
-  - colorama
-  - imageio
-  - ipywidgets>=7.0.1
-  - joblib
+  - h5netcdf>=1.2.0
+  - igor2>=0.5.6
+  - iminuit>=2.25.2
+  - ipykernel
+  - joblib>=1.3.2
   - lmfit>=1.2.0
-  - pandas
-  - pint
-  - psutil
-  - pygments
+  - matplotlib>=3.8.0
+  - numba-progress>=1.0.0
+  - numba>=0.59.0
+  - numbagg>=0.8.1
+  - numpy>=1.26.0
+  - pip>=19.3
+  - pyperclip>=1.8.2
   - pyqtgraph>=0.13.1
-  - rx
-  - scikit-image
-  - scikit-learn>=0.24.0
-  - tbb
-  - titlecase
-  - tqdm
+  - python>=3.11
+  - qtawesome>=1.3.1
+  - qtpy>=2.4.1
+  - scipy>=1.12.0
+  - superqt>=0.6.2
+  - tqdm>=4.66.2
   - uncertainties>=3.0.1
-  - qtawesome
-  - numbagg
-  - numba-progress
-  - iminuit>=2.3
-  - pyperclip
+  - varname>=0.13.0
+  - xarray>=2024.02.0
   - pip:
     - pyqt6>=6.2.2
-    - igor2 @ git+https://github.com/kmnhan/igor2@master
-    - arpes @ git+https://github.com/kmnhan/arpes@master
diff --git a/environment_apple.yml b/environment_apple.yml
index d8cad1f0..270fa119 100644
--- a/environment_apple.yml
+++ b/environment_apple.yml
@@ -2,55 +2,32 @@ name: erlab
 channels:
   - conda-forge
 dependencies:
-  - python>=3.11
-  - ipython
-  - pip>=19.3
-  - setuptools
+  - h5netcdf>=1.2.0
+  - igor2>=0.5.6
+  - iminuit>=2.25.2
+  - ipykernel
+  - joblib>=1.3.2
   - libblas=*=*accelerate
   - libcblas=*=*accelerate
   - liblapack=*=*accelerate
   - liblapacke=*=*accelerate
-  - h5netcdf
-  - varname
-  - numpy>=1.22.0
-  - scipy>=1.8.0
-  - matplotlib>=3.7.0
-  - h5py>=3.2.1
-  - netCDF4>=1.5.0
-  - astropy
-  - superqt>=0.3.1
-  - dill
-  - llvm-openmp>=13.0.0
-  - llvmlite
-  - asteval>=0.9.22
-  - xarray>=0.16.1
-  - numba>=0.56.2
-  - ase>=3.17.0
-  - bokeh>=2.0.0
-  - bottleneck>=1.3.2
-  - colorama
-  - imageio
-  - ipywidgets>=7.0.1
-  - joblib
   - lmfit>=1.2.0
-  - pandas
-  - pint
-  - psutil
-  - pygments
+  - matplotlib>=3.8.0
+  - numba-progress>=1.0.0
+  - numba>=0.59.0
+  - numbagg>=0.8.1
+  - numpy>=1.26.0
+  - pip>=19.3
+  - pyperclip>=1.8.2
   - pyqtgraph>=0.13.1
-  - rx
-  - scikit-image
-  - scikit-learn>=0.24.0
-  - tbb
-  - titlecase
-  - tqdm
+  - python>=3.11
+  - qtawesome>=1.3.1
+  - qtpy>=2.4.1
+  - scipy>=1.12.0
+  - superqt>=0.6.2
+  - tqdm>=4.66.2
   - uncertainties>=3.0.1
-  - qtawesome
-  - numbagg
-  - numba-progress
-  - iminuit>=2.3
-  - pyperclip
+  - varname>=0.13.0
+  - xarray>=2024.02.0
   - pip:
     - pyqt6>=6.2.2
-    - igor2 @ git+https://github.com/kmnhan/igor2@master
-    - arpes @ git+https://github.com/kmnhan/arpes@master
diff --git a/environment_nogit.yml b/environment_nogit.yml
index 92298fa9..270fa119 100644
--- a/environment_nogit.yml
+++ b/environment_nogit.yml
@@ -2,53 +2,32 @@ name: erlab
 channels:
   - conda-forge
 dependencies:
-  - python>=3.11
-  - ipython
-  - pip>=19.3
-  - setuptools
+  - h5netcdf>=1.2.0
+  - igor2>=0.5.6
+  - iminuit>=2.25.2
+  - ipykernel
+  - joblib>=1.3.2
   - libblas=*=*accelerate
   - libcblas=*=*accelerate
   - liblapack=*=*accelerate
   - liblapacke=*=*accelerate
-  - h5netcdf
-  - varname
-  - numpy>=1.22.0
-  - scipy>=1.8.0
-  - matplotlib>=3.7.0
-  - h5py>=3.2.1
-  - netCDF4>=1.5.0
-  - astropy
-  - superqt>=0.3.1
-  - dill
-  - llvm-openmp>=13.0.0
-  - llvmlite
-  - asteval>=0.9.22
-  - xarray>=0.16.1
-  - numba>=0.56.2
-  - ase>=3.17.0
-  - bokeh>=2.0.0
-  - bottleneck>=1.3.2
-  - colorama
-  - imageio
-  - ipywidgets>=7.0.1
-  - joblib
   - lmfit>=1.2.0
-  - pandas
-  - pint
-  - psutil
-  - pygments
+  - matplotlib>=3.8.0
+  - numba-progress>=1.0.0
+  - numba>=0.59.0
+  - numbagg>=0.8.1
+  - numpy>=1.26.0
+  - pip>=19.3
+  - pyperclip>=1.8.2
   - pyqtgraph>=0.13.1
-  - rx
-  - scikit-image
-  - scikit-learn>=0.24.0
-  - tbb
-  - titlecase
-  - tqdm
+  - python>=3.11
+  - qtawesome>=1.3.1
+  - qtpy>=2.4.1
+  - scipy>=1.12.0
+  - superqt>=0.6.2
+  - tqdm>=4.66.2
   - uncertainties>=3.0.1
-  - qtawesome
-  - numbagg
-  - numba-progress
-  - iminuit>=2.3
-  - pyperclip
+  - varname>=0.13.0
+  - xarray>=2024.02.0
   - pip:
     - pyqt6>=6.2.2
diff --git a/environment_nogit_mkl.yml b/environment_nogit_mkl.yml
index cd533bf2..6c3debc7 100644
--- a/environment_nogit_mkl.yml
+++ b/environment_nogit_mkl.yml
@@ -2,53 +2,32 @@ name: erlab
 channels:
   - conda-forge
 dependencies:
-  - python>=3.11
-  - ipython
-  - pip>=19.3
-  - setuptools
+  - h5netcdf>=1.2.0
+  - igor2>=0.5.6
+  - iminuit>=2.25.2
+  - ipykernel
+  - joblib>=1.3.2
   - libblas=*=*mkl
   - libcblas=*=*mkl
   - liblapack=*=*mkl
   - liblapacke=*=*mkl
-  - h5netcdf
-  - varname
-  - numpy>=1.22.0
-  - scipy>=1.8.0
-  - matplotlib>=3.7.0
-  - h5py>=3.2.1
-  - netCDF4>=1.5.0
-  - astropy
-  - superqt>=0.3.1
-  - dill
-  - llvm-openmp>=13.0.0
-  - llvmlite
-  - asteval>=0.9.22
-  - xarray>=0.16.1
-  - numba>=0.56.2
-  - ase>=3.17.0
-  - bokeh>=2.0.0
-  - bottleneck>=1.3.2
-  - colorama
-  - imageio
-  - ipywidgets>=7.0.1
-  - joblib
   - lmfit>=1.2.0
-  - pandas
-  - pint
-  - psutil
-  - pygments
+  - matplotlib>=3.8.0
+  - numba-progress>=1.0.0
+  - numba>=0.59.0
+  - numbagg>=0.8.1
+  - numpy>=1.26.0
+  - pip>=19.3
+  - pyperclip>=1.8.2
   - pyqtgraph>=0.13.1
-  - rx
-  - scikit-image
-  - scikit-learn>=0.24.0
-  - tbb
-  - titlecase
-  - tqdm
+  - python>=3.11
+  - qtawesome>=1.3.1
+  - qtpy>=2.4.1
+  - scipy>=1.12.0
+  - superqt>=0.6.2
+  - tqdm>=4.66.2
   - uncertainties>=3.0.1
-  - qtawesome
-  - numbagg
-  - numba-progress
-  - iminuit>=2.3
-  - pyperclip
+  - varname>=0.13.0
+  - xarray>=2024.02.0
   - pip:
     - pyqt6>=6.2.2
diff --git a/pyproject.toml b/pyproject.toml
index 219f8798..3539fcba 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -21,25 +21,26 @@ classifiers = [
     "Operating System :: OS Independent",
 ]
 dependencies = [
-    "numpy>=1.22.0",
-    "numba>=0.56.2",
-    "numbagg",
-    "numba-progress",
-    "joblib",
-    "scipy>=1.8.0",
-    "lmfit>=1.2.0",
-    "uncertainties>=3.0.1",
+    "h5netcdf>=1.2.0",
+    "igor2>=0.5.6",
     "iminuit>=2.25.2",
+    "joblib>=1.3.2",
+    "lmfit>=1.2.0",
     "matplotlib>=3.8.0",
-    "xarray",
-    "h5netcdf",
-    "varname",
-    "tqdm",
-    "pyperclip",
-    "qtpy",
+    "numba-progress>=1.0.0",
+    "numba>=0.59.0",
+    "numbagg>=0.8.1",
+    "numpy>=1.26.0",
+    "pyperclip>=1.8.2",
     "pyqtgraph>=0.13.1",
-    "qtawesome",
-    "superqt",
+    "qtawesome>=1.3.1",
+    "qtpy>=2.4.1",
+    "scipy>=1.12.0",
+    "superqt>=0.6.2",
+    "tqdm>=4.66.2",
+    "uncertainties>=3.0.1",
+    "varname>=0.13.0",
+    "xarray>=2024.02.0",
 ]
 dynamic = ["version"]
 
@@ -58,6 +59,8 @@ package-dir = { "" = "src" }
 [tool.setuptools.dynamic]
 version = { attr = "erlab.__version__" }
 
+[tool.setuptools_scm]
+
 [tool.semantic_release]
 assets = []
 build_command = "pip install build && python -m build"
@@ -190,3 +193,6 @@ docstring-code-line-length = "dynamic"
 
 [tool.isort]
 profile = "black"
+
+[tool.pytest.ini_options]
+testpaths = "tests"
diff --git a/requirements.txt b/requirements.txt
index 24b06ba3..0e26a9ab 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -1,22 +1,21 @@
-numpy>=1.22.0
-numba>=0.56.2
-numbagg
-numba-progress
-joblib
-scipy>=1.8.0
-lmfit>=1.2.0
-uncertainties>=3.0.1
+h5netcdf>=1.2.0
+igor2>=0.5.6
 iminuit>=2.25.2
+joblib>=1.3.2
+lmfit>=1.2.0
 matplotlib>=3.8.0
-xarray
-h5netcdf
-varname
-tqdm
-pyperclip
+numba-progress>=1.0.0
+numba>=0.59.0
+numbagg>=0.8.1
+numpy>=1.26.0
+pyperclip>=1.8.2
 pyqt6>=6.2.2
-qtpy
 pyqtgraph>=0.13.1
-qtawesome
-superqt
-igor2 @ git+https://github.com/kmnhan/igor2@master
-arpes @ git+https://github.com/kmnhan/arpes@master
+qtawesome>=1.3.1
+qtpy>=2.4.1
+scipy>=1.12.0
+superqt>=0.6.2
+tqdm>=4.66.2
+uncertainties>=3.0.1
+varname>=0.13.0
+xarray>=2024.02.0
diff --git a/src/erlab/accessors.py b/src/erlab/accessors.py
index b524aca4..83979f17 100644
--- a/src/erlab/accessors.py
+++ b/src/erlab/accessors.py
@@ -798,7 +798,7 @@ def best_kp_resolution(self) -> float:
     def best_kz_resolution(self) -> float:
         r"""
         Estimates the minimum out-of-plane momentum resolution based on the mean free
-        path :cite:p:`Seah1979` and the kinetic energy.
+        path :cite:p:`seah1979imfp` and the kinetic energy.
 
         .. math:: \Delta k_z \sim 1/\lambda
 
diff --git a/src/erlab/analysis/__init__.py b/src/erlab/analysis/__init__.py
index 669b6e85..d971a3e9 100644
--- a/src/erlab/analysis/__init__.py
+++ b/src/erlab/analysis/__init__.py
@@ -13,6 +13,7 @@
    mask
    correlation
    gold
+   image
    interpolate
    kspace
    transform
diff --git a/src/erlab/analysis/fit/functions/__init__.py b/src/erlab/analysis/fit/functions/__init__.py
index f8e7ab04..2078df77 100644
--- a/src/erlab/analysis/fit/functions/__init__.py
+++ b/src/erlab/analysis/fit/functions/__init__.py
@@ -24,9 +24,9 @@
     "fermi_dirac_linbkg_broad",
     "step_linbkg_broad",
     "step_broad",
-    "PolyFunc",
+    "PolynomialFunction",
     "MultiPeakFunction",
-    "FermiEdge2dFunc",
+    "FermiEdge2dFunction",
 ]
 
 from erlab.analysis.fit.functions.general import (
@@ -43,7 +43,7 @@
 )
 
 from erlab.analysis.fit.functions.dynamic import (
-    PolyFunc,
+    PolynomialFunction,
     MultiPeakFunction,
-    FermiEdge2dFunc,
+    FermiEdge2dFunction,
 )
diff --git a/src/erlab/analysis/fit/functions/dynamic.py b/src/erlab/analysis/fit/functions/dynamic.py
index 18c2f071..c37b58b5 100644
--- a/src/erlab/analysis/fit/functions/dynamic.py
+++ b/src/erlab/analysis/fit/functions/dynamic.py
@@ -7,9 +7,9 @@
 __all__ = [
     "get_args_kwargs",
     "DynamicFunction",
-    "PolyFunc",
+    "PolynomialFunction",
     "MultiPeakFunction",
-    "FermiEdge2dFunc",
+    "FermiEdge2dFunction",
 ]
 import functools
 import inspect
@@ -95,7 +95,7 @@ def __call__(self, x: npt.NDArray[np.float64], **params) -> npt.NDArray[np.float
         raise NotImplementedError("Must be overloaded in child classes")
 
 
-class PolyFunc(DynamicFunction):
+class PolynomialFunction(DynamicFunction):
     """A callable class for a arbitrary degree polynomial.
 
     Parameters
@@ -119,7 +119,9 @@ def __call__(self, x, *coeffs, **params):
         if isinstance(x, np.ndarray):
             return np.polynomial.polynomial.polyval(x, coeffs)
         else:
-            coeffs = xr.DataArray(coeffs, coords={"degree": np.arange(self.degree + 1)})
+            coeffs = xr.DataArray(
+                np.asarray(coeffs), coords={"degree": np.arange(self.degree + 1)}
+            )
             return xr.polyval(x, coeffs)
 
 
@@ -261,16 +263,23 @@ def pre_call(
     def __call__(
         self, x: npt.NDArray[np.float64], **params: dict
     ) -> npt.NDArray[np.float64]:
+        if isinstance(x, xr.DataArray):
+            return x * 0.0 + self.__call__(x.values, **params)
+
         if self.convolve:
+            if "resolution" not in params:
+                raise TypeError(
+                    "Missing parameter `resolution` required for convolution"
+                )
             return do_convolve(x, self.pre_call, **params)
         else:
             return self.pre_call(x, **params)
 
 
-class FermiEdge2dFunc(DynamicFunction):
+class FermiEdge2dFunction(DynamicFunction):
     def __init__(self, degree=1) -> None:
         super().__init__()
-        self.poly = PolyFunc(degree)
+        self.poly = PolynomialFunction(degree)
 
     @property
     def argnames(self) -> list[str]:
@@ -296,4 +305,15 @@ def pre_call(self, eV, alpha, **params: dict):
         ) + params["offset"]
 
     def __call__(self, eV, alpha, **params: dict):
+        if isinstance(eV, xr.DataArray) and isinstance(alpha, xr.DataArray):
+            out = eV * alpha * 0.0
+            return out + self.__call__(eV.values, alpha.values, **params).reshape(
+                out.shape
+            )
+        if isinstance("eV", xr.DataArray):
+            eV = eV.values
+        if isinstance("alpha", xr.DataArray):
+            alpha = alpha.values
+        if "resolution" not in params:
+            raise TypeError("Missing parameter `resolution` required for convolution")
         return do_convolve_y(eV, alpha, self.pre_call, **params).flatten()
diff --git a/src/erlab/analysis/fit/models.py b/src/erlab/analysis/fit/models.py
index 9c666cdd..c71a0a4b 100644
--- a/src/erlab/analysis/fit/models.py
+++ b/src/erlab/analysis/fit/models.py
@@ -8,9 +8,6 @@
     "FermiEdge2dModel",
 ]
 
-
-import importlib
-
 import lmfit
 import numba
 import numpy as np
@@ -18,15 +15,10 @@
 import scipy.ndimage
 import xarray as xr
 
-if importlib.util.find_spec("arpes"):
-    from arpes.fits import XModelMixin
-else:
-    from lmfit.models import Model as XModelMixin
-
 from erlab.analysis.fit.functions import (
-    FermiEdge2dFunc,
+    FermiEdge2dFunction,
     MultiPeakFunction,
-    PolyFunc,
+    PolynomialFunction,
     fermi_dirac_linbkg_broad,
     step_linbkg_broad,
 )
@@ -108,58 +100,51 @@ def fit_edges_linear(x, data, len_fit) -> tuple[float, float, float, float]:
     return n0, m0, n1, m1
 
 
-class ExtendedAffineBroadenedFD(XModelMixin):
+class ExtendedAffineBroadenedFD(lmfit.Model):
     """
     Fermi-dirac function with linear background above and below the fermi level,
     convolved with a gaussian kernel.
     """
 
-    guess_dataarray = True
-
     @staticmethod
     def LinearBroadFermiDirac(
         x,
-        center=0,
-        temp=30,
+        center=0.0,
+        temp=30.0,
         resolution=0.02,
-        back0=1,
-        back1=0,
-        dos0=1,
-        dos1=0,
+        back0=0.0,
+        back1=0.0,
+        dos0=1.0,
+        dos1=0.0,
     ):
         return fermi_dirac_linbkg_broad(
             x, center, temp, resolution, back0, back1, dos0, dos1
         )
 
-    def __init__(
-        self, independent_vars=("x",), prefix="", missing="raise", name=None, **kwargs
-    ):
-        """Defer to lmfit for initialization."""
-        kwargs.update(
-            {"prefix": prefix, "missing": missing, "independent_vars": independent_vars}
-        )
+    def __init__(self, **kwargs):
         super().__init__(self.LinearBroadFermiDirac, **kwargs)
         self.set_param_hint("temp", min=0.0)
         self.set_param_hint("resolution", min=0.0)
 
     def guess(self, data, x, **kwargs):
-        """Make some heuristic guesses."""
         pars = self.make_params()
 
         len_fit = max(round(len(x) * 0.05), 10)
         dos0, dos1, back0, back1 = fit_edges_linear(x, data, len_fit)
-        efermi = x[
+        efermi = np.asarray(x)[
             np.argmin(np.gradient(scipy.ndimage.gaussian_filter1d(data, 0.2 * len(x))))
         ]
 
-        temp = 30
+        temp = 30.0
         if isinstance(data, xr.DataArray):
             try:
-                temp = data.attrs["temp_sample"]
+                temp = float(data.attrs["temp_sample"])
             except KeyError:
                 pass
 
-        pars[f"{self.prefix}center"].set(value=efermi, min=x.min(), max=x.max())
+        pars[f"{self.prefix}center"].set(
+            value=efermi, min=np.asarray(x).min(), max=np.asarray(x).max()
+        )
         pars[f"{self.prefix}back0"].set(value=back0)
         pars[f"{self.prefix}back1"].set(value=back1)
         pars[f"{self.prefix}dos0"].set(value=dos0)
@@ -173,13 +158,8 @@ def guess(self, data, x, **kwargs):
     guess.__doc__ = lmfit.models.COMMON_GUESS_DOC
 
 
-class StepEdgeModel(XModelMixin):
-    guess_dataarray = True
-
-    def __init__(
-        self, independent_vars=("x",), prefix="", missing="raise", name=None, **kwargs
-    ):
-        """Defer to lmfit for initialization."""
+class StepEdgeModel(lmfit.Model):
+    def __init__(self, independent_vars=("x",), prefix="", missing="raise", **kwargs):
         kwargs.update(
             {"prefix": prefix, "missing": missing, "independent_vars": independent_vars}
         )
@@ -187,7 +167,6 @@ def __init__(
         self.set_param_hint("sigma", min=0.0)
 
     def guess(self, data, x, **kwargs):
-        """Make some heuristic guesses."""
         pars = self.make_params()
 
         len_fit = max(round(len(x) * 0.05), 10)
@@ -209,10 +188,10 @@ def guess(self, data, x, **kwargs):
     guess.__doc__ = lmfit.models.COMMON_GUESS_DOC
 
 
-class PolynomialModel(XModelMixin):
+class PolynomialModel(lmfit.Model):
     def __init__(self, degree=9, **kwargs):
         kwargs.setdefault("name", f"Poly{degree}")
-        super().__init__(PolyFunc(degree), **kwargs)
+        super().__init__(PolynomialFunction(degree), **kwargs)
 
     def guess(self, data, x=None, **kwargs):
         pars = self.make_params()
@@ -230,11 +209,11 @@ def guess(self, data, x=None, **kwargs):
     guess.__doc__ = lmfit.models.COMMON_GUESS_DOC
 
 
-class MultiPeakModel(XModelMixin):
+class MultiPeakModel(lmfit.Model):
     """Model for fitting multiple Gaussian or Lorentzian peaks.
 
-    Most input parameters are passed to the
-    `erlab.analysis.fit.functions.MultiPeakFunction` constructor.
+    Most input parameters are passed to the :class:`MultiPeakFunction
+    <erlab.analysis.fit.functions.MultiPeakFunction>` constructor.
     """
 
     def __init__(
@@ -243,23 +222,13 @@ def __init__(
         peak_shapes: list[str] | str | None = None,
         fd: bool = True,
         convolve: bool = True,
-        independent_vars=None,
-        param_names=None,
-        nan_policy="raise",
-        prefix="",
-        name=None,
-        **kws,
+        **kwargs,
     ):
         super().__init__(
             MultiPeakFunction(
                 npeaks, peak_shapes=peak_shapes, fd=fd, convolve=convolve
             ),
-            independent_vars,
-            param_names,
-            nan_policy,
-            prefix,
-            name,
-            **kws,
+            **kwargs,
         )
 
         if self.func.convolve:
@@ -290,7 +259,7 @@ def guess(self, data, x=None, **kwargs):
     guess.__doc__ = lmfit.models.COMMON_GUESS_DOC
 
 
-class FermiEdge2dModel(XModelMixin):
+class FermiEdge2dModel(lmfit.Model):
     r"""A 2D model for a polynomial Fermi edge with a linear density of states.
 
     The model function can be written as
@@ -307,30 +276,17 @@ class FermiEdge2dModel(XModelMixin):
 
     """
 
-    n_dims = 2
-    dimension_order = ["eV", "alpha"]
-    guess_dataarray = True
-    fit_flat = True
-
     def __init__(
         self,
-        degree=2,
+        degree: int = 2,
         independent_vars=("eV", "alpha"),
-        prefix="",
-        nan_policy="raise",
         **kwargs,
     ):
-        kwargs.update(
-            {
-                "prefix": prefix,
-                "nan_policy": nan_policy,
-                "independent_vars": independent_vars,
-            }
-        )
-        super().__init__(FermiEdge2dFunc(degree), **kwargs)
+        kwargs.update({"independent_vars": independent_vars})
+        super().__init__(FermiEdge2dFunction(degree), **kwargs)
         self.name = f"FermiEdge2dModel (deg {degree})"
 
-    def guess(self, data, eV, alpha, negative=False, **kwargs):
+    def guess(self, data, eV, alpha, **kwargs):
         pars = self.make_params()
         for i in range(self.func.poly.degree + 1):
             pars[f"{self.prefix}c{i}"].set(value=0.0)
@@ -348,8 +304,11 @@ def guess(self, data, eV, alpha, negative=False, **kwargs):
 
         return lmfit.models.update_param_vals(pars, self.prefix, **kwargs)
 
-    def guess_fit(self, *args, **kwargs):
-        return super().guess_fit(*args, **kwargs)
+    def fit(self, data, *args, **kwargs):
+        if isinstance(data, xr.DataArray):
+            data = data.transpose("eV", "alpha").values
+        # Ensure flat fit
+        return super().fit(data.ravel(), *args, **kwargs)
 
     __init__.__doc__ = lmfit.models.COMMON_INIT_DOC.replace("['x']", "['eV', 'alpha']")
     guess.__doc__ = lmfit.models.COMMON_GUESS_DOC
diff --git a/src/erlab/analysis/gold.py b/src/erlab/analysis/gold.py
index 64935af1..cd3b9df1 100644
--- a/src/erlab/analysis/gold.py
+++ b/src/erlab/analysis/gold.py
@@ -9,9 +9,9 @@
     "resolution_roi",
 ]
 
-import importlib
 from collections.abc import Sequence
 
+import joblib
 import lmfit.model
 import matplotlib
 import matplotlib.figure
@@ -29,11 +29,9 @@
     StepEdgeModel,
 )
 from erlab.analysis.utilities import correct_with_edge
+from erlab.parallel import joblib_progress
 from erlab.plotting.general import autoscale_to, figwh
 
-if importlib.util.find_spec("arpes"):
-    import arpes.fits
-
 
 def edge(
     gold: xr.DataArray | xr.Dataset,
@@ -46,6 +44,7 @@ def edge(
     method: str = "leastsq",
     progress: bool = True,
     parallel_kw: dict | None = None,
+    parallel_obj: joblib.Parallel | None = None,
     return_full: bool = False,
     fixed_center: float | None = None,
     scale_covar: bool = True,
@@ -78,6 +77,9 @@ def edge(
         Whether to display the fitting progress, by default `True`.
     parallel_kw
         Additional keyword arguments for parallel fitting, by default `None`.
+    parallel_obj
+        The `joblib.Parallel` object to use for fitting, by default `None`. If provided,
+        `parallel_kw` will be ignored.
     return_full
         Whether to return the full fit results, by default `False`.
     fixed_center
@@ -99,21 +101,21 @@ def edge(
 
     """
     if fast:
-        params = {}
+        params = lmfit.create_params()
         model_cls = StepEdgeModel
     else:
         if temp is None:
             temp = gold.attrs["temp_sample"]
-        params = {
-            "temp": {"value": temp, "vary": vary_temp},
-        }
+        params = lmfit.create_params(temp={"value": temp, "vary": vary_temp})
         model_cls = ExtendedAffineBroadenedFD
 
+    model = model_cls()
+
     if parallel_kw is None:
         parallel_kw = {}
 
     if fixed_center is not None:
-        params["center"] = {"value": fixed_center, "vary": False}
+        params["center"].set(value=fixed_center, vary=False)
 
     if any(b != 1 for b in bin_size):
         gold_binned = gold.coarsen(alpha=bin_size[0], eV=bin_size[1], boundary="trim")
@@ -122,34 +124,76 @@ def edge(
     gold_sel = gold.sel(alpha=slice(*angle_range), eV=slice(*eV_range))
 
     # Assuming Poisson noise, the weights are the square root of the counts.
-    weights = np.sqrt(gold.sum("eV"))
+    weights = 1 / np.sqrt(gold_sel.sum("eV").values)
+
+    n_fits = len(gold_sel.alpha)
+
+    if parallel_obj is None:
+        if n_fits > 20:
+            parallel_kw.setdefault("n_jobs", -1)
+        else:
+            parallel_kw.setdefault("n_jobs", 1)
+        parallel_obj = joblib.Parallel(**parallel_kw)
+
+    def _fit(data, w):
+        pars = model.guess(data, x=data["eV"]).update(params)
+
+        res = model.fit(
+            data,
+            x=data["eV"],
+            params=pars,
+            method=method,
+            scale_covar=scale_covar,
+            weights=w,
+        )
+        return res
+
+    if progress:
+        with joblib_progress(desc="Fitting", total=n_fits) as _:
+            fitresults = parallel_obj(
+                joblib.delayed(_fit)(gold_sel.isel(alpha=i), weights[i])
+                for i in range(n_fits)
+            )
+    else:
+        fitresults = parallel_obj(
+            joblib.delayed(_fit)(gold_sel.isel(alpha=i), weights[i])
+            for i in range(n_fits)
+        )
 
-    fitresults = arpes.fits.broadcast_model(
-        model_cls,
-        gold_sel,
-        "alpha",
-        weights=weights,
-        params=params,
-        method=method,
-        progress=progress,
-        parallel_kw=parallel_kw,
-        scale_covar=scale_covar,
-        **kwargs,
-    )
     if return_full:
         return fitresults
 
-    center_arr, center_stderr = fitresults.F.p("center"), fitresults.F.s("center")
-    center_arr = center_arr.where(~center_stderr.isnull(), drop=True)
-    center_stderr = center_stderr.where(~center_stderr.isnull(), drop=True)
-    return center_arr, center_stderr
+    xval = []
+    res_vals = []
+
+    for i, r in enumerate(fitresults):
+        center_ufloat = r.uvars["center"]
+
+        if not np.isnan(center_ufloat.std_dev):
+            xval.append(gold_sel.alpha.values[i])
+            res_vals.append([center_ufloat.nominal_value, center_ufloat.std_dev])
+
+    xval = np.asarray(xval)
+    yval, yerr = np.asarray(res_vals).T
+
+    return (
+        xr.DataArray(yval, coords={"alpha": xval}),
+        xr.DataArray(yerr, coords={"alpha": xval}),
+    )
 
 
 def poly_from_edge(
     center, weights=None, degree=4, method="leastsq", scale_covar=True
 ) -> lmfit.model.ModelResult:
-    modelresult = PolynomialModel(degree=degree).guess_fit(
-        center, weights=weights, method=method, scale_covar=scale_covar
+    model = PolynomialModel(degree=degree)
+    pars = model.guess(center.values, x=center[center.dims[0]].values)
+    modelresult = model.fit(
+        center,
+        x=center[center.dims[0]].values,
+        params=pars,
+        weights=weights,
+        method=method,
+        scale_covar=scale_covar,
     )
     return modelresult
 
@@ -391,12 +435,14 @@ def resolution(
     gold_roi = gold_corr.sel(alpha=slice(*angle_range))
     edc_avg = gold_roi.mean("alpha").sel(eV=slice(*eV_range_fit))
 
-    params = {
-        "temp": {"value": gold.attrs["temp_sample"], "vary": False},
-        "resolution": {"value": 0.1, "vary": True, "min": 0},
-    }
-    fit = ExtendedAffineBroadenedFD().guess_fit(
-        edc_avg, params=params, method=method, scale_covar=scale_covar
+    params = lmfit.create_params(
+        temp={"value": gold_roi.attrs["temp_sample"], "vary": False},
+        resolution={"value": 0.1, "vary": True, "min": 0},
+    )
+    model = ExtendedAffineBroadenedFD()
+    params = model.guess(edc_avg, x=edc_avg["eV"]).update(params)
+    fit = ExtendedAffineBroadenedFD().fit(
+        edc_avg, x=edc_avg["eV"], params=params, method=method, scale_covar=scale_covar
     )
     if plot:
         plt.show()
@@ -437,12 +483,15 @@ def resolution_roi(
 ) -> lmfit.model.ModelResult:
     edc_avg = gold_roi.mean("alpha").sel(eV=slice(*eV_range))
 
-    params = {
-        "temp": {"value": gold_roi.attrs["temp_sample"], "vary": not fix_temperature},
-        "resolution": {"value": 0.1, "vary": True, "min": 0},
-    }
-    fit = ExtendedAffineBroadenedFD().guess_fit(
+    params = lmfit.create_params(
+        temp={"value": gold_roi.attrs["temp_sample"], "vary": not fix_temperature},
+        resolution={"value": 0.1, "vary": True, "min": 0},
+    )
+    model = ExtendedAffineBroadenedFD()
+    params = model.guess(edc_avg, x=edc_avg["eV"]).update(params)
+    fit = model.fit(
         edc_avg,
+        x=edc_avg["eV"],
         params=params,
         method=method,
         scale_covar=scale_covar,  # weights=1 / edc_stderr
diff --git a/src/erlab/analysis/image.py b/src/erlab/analysis/image.py
new file mode 100644
index 00000000..af33f2d4
--- /dev/null
+++ b/src/erlab/analysis/image.py
@@ -0,0 +1,522 @@
+"""
+Various image processing functions including tools for visualizing dispersive features.
+
+Note
+----
+For scipy-based filter functions, the default value of the `mode` argument is 'nearest',
+unlike the scipy default of 'reflect'.
+"""
+
+from collections.abc import Sequence
+
+import numpy as np
+import numpy.typing as npt
+import scipy
+import scipy.ndimage
+import xarray as xr
+from numba import carray, cfunc, types
+
+
+def gaussian_filter(
+    darr: xr.DataArray,
+    sigma: float | dict[str, float] | Sequence[float],
+    order: int | Sequence[int] | dict[str, int] = 0,
+    mode: str | Sequence[str] | dict[str, str] = "nearest",
+    cval: float = 0.0,
+    truncate: float = 4.0,
+    *,
+    radius: None | float | Sequence[float] | dict[str, float] = None,
+) -> xr.DataArray:
+    """Coordinate-aware wrapper around `scipy.ndimage.gaussian_filter`.
+
+    Parameters
+    ----------
+    darr
+        The input DataArray.
+    sigma
+        The standard deviation(s) of the Gaussian filter in data dimensions. If a float,
+        the same value is used for all dimensions, each scaled by the data step. If a
+        dict, the value can be specified for each dimension using dimension names as
+        keys. The filter is only applied to the dimensions specified in the dict. If a
+        sequence, the values are used in the same order as the dimensions of the
+        DataArray.
+    order
+        The order of the filter along each dimension. If an int, the same order is used
+        for all dimensions. See Notes below for other options. Defaults to 0.
+    mode
+        The boundary mode used for the filter. If a str, the same mode is used for all
+        dimensions. See Notes below for other options. Defaults to 'nearest'.
+    cval
+        Value to fill past edges of input if mode is 'constant'. Defaults to 0.0.
+    truncate
+        The truncation value used for the Gaussian filter. Defaults to 4.0.
+    radius
+        The radius of the Gaussian filter in data units. See Notes below. Defaults to
+        None.
+
+    Returns
+    -------
+    gaussian_filter : xarray.DataArray
+        The filtered array with the same shape as the input DataArray.
+
+    Note
+    ----
+    - The `sigma` and `radius` values should be in data coordinates, not pixels.
+    - The input array is assumed to be regularly spaced.
+    - `order`, `mode`, and `radius` can be specified for each dimension using a dict or
+      a sequence. If a dict, the value can be specified for each dimension using
+      dimension names as keys. If a sequence and `sigma` is given as a dictionary, the
+      order is assumed to be the same as the keys in `sigma`. If `sigma` is not a
+      dictionary, the order is assumed to be the same as the dimensions of the
+      DataArray.
+
+    See Also
+    --------
+    :func:`scipy.ndimage.gaussian_filter` : The underlying function used to apply the
+        filter.
+
+    Example
+    -------
+    >>> import numpy as np
+    >>> import xarray as xr
+    >>> import erlab.analysis as era
+    >>> darr = xr.DataArray(np.arange(50, step=2).reshape((5, 5)), dims=["x", "y"])
+    >>> darr
+    <xarray.DataArray (x: 5, y: 5)> Size: 200B
+    array([[ 0,  2,  4,  6,  8],
+        [10, 12, 14, 16, 18],
+        [20, 22, 24, 26, 28],
+        [30, 32, 34, 36, 38],
+        [40, 42, 44, 46, 48]])
+    Dimensions without coordinates: x, y
+    >>> era.image.gaussian_filter(darr, sigma=dict(x=1.0, y=1.0))
+    <xarray.DataArray (x: 5, y: 5)> Size: 200B
+    array([[ 3,  5,  7,  8, 10],
+        [10, 12, 14, 15, 17],
+        [20, 22, 24, 25, 27],
+        [29, 31, 33, 34, 36],
+        [36, 38, 40, 41, 43]])
+    Dimensions without coordinates: x, y
+
+    """
+    if np.isscalar(sigma):
+        sigma = {d: sigma for d in darr.dims}
+    elif not isinstance(sigma, dict):
+        sigma = dict(zip(darr.dims, sigma))
+
+    # Get the axis indices to apply the filter
+    axes = tuple(darr.get_axis_num(d) for d in sigma.keys())
+
+    # Convert arguments to tuples acceptable by scipy
+    if isinstance(order, dict):
+        order = tuple(order.get(d, 0) for d in sigma.keys())
+    if isinstance(mode, dict):
+        mode = tuple(mode[d] for d in sigma.keys())
+    if radius is not None:
+        if len(radius) != len(sigma):
+            raise ValueError("`radius` does not match dimensions of `sigma`")
+
+        if np.isscalar(radius):
+            radius = {d: radius for d in sigma.keys()}
+        elif not isinstance(radius, dict):
+            radius = dict(zip(sigma.keys(), radius))
+
+        # Calculate radius in pixels
+        radius: tuple[int, ...] = tuple(
+            round(r / (darr[d].values[1] - darr[d].values[0]))
+            for d, r in radius.items()
+        )
+
+    # Calculate sigma in pixels
+    sigma: tuple[float, ...] = tuple(
+        val / (darr[d].values[1] - darr[d].values[0]) for d, val in sigma.items()
+    )
+
+    return darr.copy(
+        data=scipy.ndimage.gaussian_filter(
+            darr.values,
+            sigma=sigma,
+            order=order,
+            mode=mode,
+            cval=cval,
+            truncate=truncate,
+            radius=radius,
+            axes=axes,
+        )
+    )
+
+
+def gaussian_laplace(
+    darr: xr.DataArray,
+    sigma: float | dict[str, float] | Sequence[float],
+    mode: str | Sequence[str] | dict[str, str] = "nearest",
+    cval: float = 0.0,
+    **kwargs,
+) -> xr.DataArray:
+    """Coordinate-aware wrapper around `scipy.ndimage.gaussian_laplace`.
+
+    This function calculates the Laplacian of the given array using Gaussian second
+    derivatives.
+
+    Parameters
+    ----------
+    darr
+        The input DataArray.
+    sigma
+        The standard deviation(s) of the Gaussian filter in data dimensions. If a float,
+        the same value is used for all dimensions, each scaled by the data step. If a
+        dict, the value can be specified for each dimension using dimension names as
+        keys. If a sequence, the values are used in the same order as the dimensions of
+        the DataArray.
+    mode
+        The mode parameter determines how the input array is extended beyond its
+        boundaries. If a string, the same mode is used for all dimensions. If a
+        sequence, the values should be the modes for each dimension in the same order as
+        the dimensions in the DataArray. If a dictionary, the keys should be dimension
+        names and the values should be the corresponding modes, and every dimension in
+        the DataArray must be present. Default is "nearest".
+    cval
+        Value to fill past edges of input if mode is 'constant'. Defaults to 0.0.
+    **kwargs
+        Additional keyword arguments to `scipy.ndimage.gaussian_filter`.
+
+    Returns
+    -------
+    gaussian_laplace : xarray.DataArray
+        The filtered array with the same shape as the input DataArray.
+
+    Note
+    ----
+    - `sigma` should be in data coordinates, not pixels.
+    - The input array is assumed to be regularly spaced.
+
+    See Also
+    --------
+    :func:`scipy.ndimage.gaussian_laplace` : The underlying function used to apply the
+        filter.
+    """
+    if np.isscalar(sigma):
+        sigma = {d: sigma for d in darr.dims}
+    elif not isinstance(sigma, dict):
+        sigma = dict(zip(darr.dims, sigma))
+
+    if len(sigma) != darr.ndim:
+        raise ValueError(
+            "`sigma` must be provided for every dimension of the DataArray"
+        )
+
+    # Convert mode to tuple acceptable by scipy
+    if isinstance(mode, dict):
+        mode = tuple(mode[d] for d in sigma.keys())
+
+    # Calculate sigma in pixels
+    sigma: tuple[float, ...] = tuple(
+        val / (darr[d].values[1] - darr[d].values[0]) for d, val in sigma.items()
+    )
+
+    return darr.copy(
+        data=scipy.ndimage.gaussian_laplace(
+            darr.values, sigma=sigma, mode=mode, cval=cval, **kwargs
+        )
+    )
+
+
+def gradient_magnitude(
+    input: npt.NDArray[np.float64],
+    dx: np.float64,
+    dy: np.float64,
+    mode: str = "nearest",
+    cval: float = 0.0,
+) -> npt.NDArray[np.float64]:
+    r"""Calculate the gradient magnitude of an input array.
+
+    The gradient magnitude is calculated as defined in Ref. :cite:p:`he2017mingrad`,
+    using given :math:`\Delta x` and :math:`\Delta y` values.
+
+    Parameters
+    ----------
+    input
+        Input array.
+    dx
+        Step size in the x-direction.
+    dy
+        Step size in the y-direction.
+    mode
+        The mode parameter controls how the gradient is calculated at the boundaries.
+        Default is 'nearest'. See `scipy.ndimage.generic_filter` for more information.
+    cval
+        The value to use for points outside the boundaries when mode is 'constant'.
+        Default is 0.0. See `scipy.ndimage.generic_filter` for more information.
+
+    Returns
+    -------
+    gradient_magnitude : numpy.ndarray
+        Gradient magnitude of the input array. Has the same shape as :code:`input`.
+
+    Note
+    ----
+    This function calculates the gradient magnitude of an input array by applying a
+    filter to the input array using the given dx and dy values. The filter is defined by
+    a kernel function that computes the squared difference between each element of the
+    input array and the central element, divided by the corresponding distance value.
+    The gradient magnitude is then calculated as the square root of the sum of the
+    squared differences.
+    """
+
+    dxy = np.sqrt(dx**2 + dy**2)
+    dist = np.array([[dxy, dy, dxy], [dx, 1.0, dx], [dxy, dy, dxy]]).flatten()
+
+    @cfunc(
+        types.intc(
+            types.CPointer(types.float64),
+            types.intp,
+            types.CPointer(types.float64),
+            types.voidptr,
+        )
+    )
+    def _kernel(values_ptr, len_values, result, data):
+        values = carray(values_ptr, (len_values,), dtype=types.float64)
+        val = 0.0
+        for i in range(9):
+            if i != 4:
+                val += ((values[i] - values[4]) / dist[i]) ** 2
+        result[0] = np.sqrt(val)
+        return 1
+
+    # https://github.com/jni/llc-tools/issues/3#issuecomment-757134814
+    func = scipy.LowLevelCallable(
+        _kernel.ctypes, signature="int (double *, npy_intp, double *, void *)"
+    )
+
+    return scipy.ndimage.generic_filter(input, func, size=(3, 3), mode=mode, cval=cval)
+
+
+def laplace(
+    darr, mode: str | Sequence[str] | dict[str, str] = "nearest", cval: float = 0.0
+) -> xr.DataArray:
+    """Coordinate-aware wrapper around `scipy.ndimage.laplace`.
+
+    This function calculates the Laplacian of the given array using approximate second
+    derivatives.
+
+    Parameters
+    ----------
+    darr
+        The input DataArray.
+    mode
+        The mode parameter determines how the input array is extended beyond its
+        boundaries. If a dictionary, the keys should be dimension names and the values
+        should be the corresponding modes, and every dimension in the DataArray must be
+        present. Otherwise, it retains the same behavior as in `scipy.ndimage.laplace`.
+        Default is 'nearest'.
+    cval
+        Value to fill past edges of input if mode is 'constant'. Defaults to 0.0.
+
+    Returns
+    -------
+    laplace : xarray.DataArray
+        The filtered array with the same shape as the input DataArray.
+
+    See Also
+    --------
+    :func:`scipy.ndimage.laplace` : The underlying function used to apply the filter.
+    """
+    if isinstance(mode, dict):
+        mode = tuple(mode[d] for d in darr.dims)
+    return darr.copy(data=scipy.ndimage.laplace(darr.values, mode=mode, cval=cval))
+
+
+def minimum_gradient(
+    darr: xr.DataArray, mode: str = "nearest", cval: float = 0.0
+) -> xr.DataArray:
+    """Minimum gradient method for detecting dispersive features in 2D data.
+
+    The minimum gradient is calculated by dividing the input DataArray by the gradient
+    magnitude. See Ref. :cite:p:`he2017mingrad`.
+
+    Parameters
+    ----------
+    darr
+        The 2D DataArray for which to calculate the minimum gradient.
+    mode
+        The mode parameter controls how the gradient is calculated at the boundaries.
+        Default is 'nearest'. See `scipy.ndimage.generic_filter` for more information.
+    cval
+        The value to use for points outside the boundaries when mode is 'constant'.
+        Default is 0.0. See `scipy.ndimage.generic_filter` for more information.
+
+    Returns
+    -------
+    minimum_gradient : xarray.DataArray
+        The minimum gradient of the input DataArray. Has the same shape as
+        :code:`input`.
+
+    Raises
+    ------
+    ValueError
+        If the input DataArray is not 2D.
+
+    Note
+    ----
+    - The input array is assumed to be regularly spaced.
+    - Any zero gradient values are replaced with NaN.
+    """
+
+    if darr.ndim != 2:
+        raise ValueError("DataArray must be 2D")
+
+    xvals = darr[darr.dims[1]].values
+    yvals = darr[darr.dims[0]].values
+
+    dx = xvals[1] - xvals[0]
+    dy = yvals[1] - yvals[0]
+
+    grad = gradient_magnitude(
+        darr.values.astype(np.float64), dx, dy, mode=mode, cval=cval
+    )
+    grad[grad == 0] = np.nan
+    return darr / grad
+
+
+def scaled_laplace(
+    darr,
+    factor: float = 1.0,
+    mode: str | Sequence[str] | dict[str, str] = "nearest",
+    cval: float = 0.0,
+) -> xr.DataArray:
+    r"""Calculate the Laplacian of a 2D DataArray with different scaling for each axis.
+
+    This function calculates the Laplacian of the given array using approximate second
+    derivatives, taking the different scaling for each axis into account.
+
+    .. math::
+        \Delta f \sim \frac{\partial^2 f}{\partial x^2} \left(\frac{\Delta x}{\Delta y}\right)^{\!2} + \frac{\partial^2 f}{\partial y^2}
+
+    See Ref. :cite:p:`zhang2011curvature` for more information.
+
+    Parameters
+    ----------
+    darr
+        The 2D DataArray for which to calculate the scaled Laplacian.
+    factor
+        The factor by which to scale the x-axis derivative. Negative values will scale
+        the y-axis derivative instead. Default is 1.0.
+    mode
+        The mode parameter determines how the input array is extended beyond its
+        boundaries. If a dictionary, the keys should be dimension names and the values
+        should be the corresponding modes, and every dimension in the DataArray must be
+        present. Otherwise, it retains the same behavior as in
+        `scipy.ndimage.generic_laplace`. Default is 'nearest'.
+    cval
+        Value to fill past edges of input if mode is 'constant'. Defaults to 0.0.
+
+    Returns
+    -------
+    scaled_laplace : xarray.DataArray
+        The filtered array with the same shape as the input DataArray.
+
+    Raises
+    ------
+    ValueError
+        If the input DataArray is not 2D.
+
+    Note
+    ----
+    The input array is assumed to be regularly spaced.
+
+    See Also
+    --------
+    :func:`scipy.ndimage.generic_laplace` : The underlying function used to apply the
+        filter.
+    """
+    if darr.ndim != 2:
+        raise ValueError("DataArray must be 2D")
+
+    xvals = darr[darr.dims[1]].values
+    yvals = darr[darr.dims[0]].values
+
+    dx = xvals[1] - xvals[0]
+    dy = yvals[1] - yvals[0]
+    weight = (dx / dy) ** 2
+
+    if factor > 0:
+        weight *= factor
+    elif factor < 0:
+        weight /= abs(factor)
+
+    if isinstance(mode, dict):
+        mode = tuple(mode[d] for d in darr.dims)
+
+    def d2_scaled(input, axis, output, mode, cval):
+        out = scipy.ndimage.correlate1d(input, [1, -2, 1], axis, output, mode, cval, 0)
+        if axis == 1:
+            out *= weight
+        return out
+
+    return darr.copy(
+        data=scipy.ndimage.generic_laplace(darr.values, d2_scaled, mode=mode, cval=cval)
+    )
+
+
+def curvature(darr: xr.DataArray, a0: float = 1.0, factor: float = 1.0) -> xr.DataArray:
+    """2D curvature method for detecting dispersive features.
+
+    The curvature is calculated as defined by :cite:t:`zhang2011curvature`.
+
+    Parameters
+    ----------
+    darr
+        The 2D DataArray for which to calculate the curvature.
+    a0
+        The regularization constant. Reasonable values range from 0.001 to 10. Default
+        is 1.0.
+    factor
+        The factor by which to scale the x-axis curvature. Negative values will scale
+        the y-axis curvature instead. Default is 1.0.
+
+    Returns
+    -------
+    curvature : xarray.DataArray
+        The 2D curvature of the input DataArray. Has the same shape as :code:`input`.
+
+    Raises
+    ------
+    ValueError
+        If the input DataArray is not 2D.
+
+    Note
+    ----
+    The input array is assumed to be regularly spaced.
+    """
+
+    if darr.ndim != 2:
+        raise ValueError("DataArray must be 2D")
+
+    xvals = darr[darr.dims[1]].values
+    yvals = darr[darr.dims[0]].values
+
+    dx = xvals[1] - xvals[0]
+    dy = yvals[1] - yvals[0]
+    weight = (dx / dy) ** 2
+
+    if factor > 0:
+        weight *= factor
+    elif factor < 0:
+        weight /= abs(factor)
+
+    dfdx, dfdy = np.gradient(darr.values, axis=(1, 0))
+    d2fdx2, d2fdydx = np.gradient(dfdx, axis=(1, 0))
+    d2fdy2 = np.gradient(dfdy, axis=0)
+
+    max_abs_dfdx_sq: float = np.max(np.abs(dfdx)) ** 2
+    max_abs_dfdy_sq: float = np.max(np.abs(dfdy)) ** 2
+
+    c0 = a0 * max(max_abs_dfdy_sq, weight * max_abs_dfdx_sq)
+
+    curv = (
+        (c0 + weight * dfdx**2) * d2fdy2
+        - 2 * weight * dfdx * dfdy * d2fdydx
+        + weight * (c0 + dfdy**2) * d2fdx2
+    ) / (c0 + weight * dfdx**2 + dfdy**2) ** 1.5
+
+    return darr.copy(data=curv)
diff --git a/src/erlab/analysis/kspace.py b/src/erlab/analysis/kspace.py
index 2e966fa6..8c1de2ca 100644
--- a/src/erlab/analysis/kspace.py
+++ b/src/erlab/analysis/kspace.py
@@ -3,7 +3,7 @@
 Typically, the user will not have to call this module directly, but will instead use the
 :func:`erlab.accessors.MomentumAccessor.convert` method.
 
-Angle conventions and function forms are based on Ref. :cite:p:`Ishida2018`.
+Angle conventions and function forms are based on Ref. :cite:p:`ishida2018kconv`.
 
 """
 
@@ -21,7 +21,7 @@
 class AxesConfiguration(enum.IntEnum):
     """
     Enum class representing different types of axes configurations. See Ref.
-    :cite:p:`Ishida2018`.
+    :cite:p:`ishida2018kconv`.
 
     """
 
@@ -67,7 +67,7 @@ def get_kconv_func(
         Dictionary of required angle parameters. If the configuration has a DA, the
         parameters should be `delta`, `chi`, `chi0`, `xi`, and `xi0`. Otherwise, they
         should be `delta`, `xi`, `xi0`, and `beta0`, following the notation in Ref.
-        :cite:p:`Ishida2018`.
+        :cite:p:`ishida2018kconv`.
 
     Returns
     -------
diff --git a/src/erlab/analysis/mask/__init__.py b/src/erlab/analysis/mask/__init__.py
index 787aeef1..a25d802a 100644
--- a/src/erlab/analysis/mask/__init__.py
+++ b/src/erlab/analysis/mask/__init__.py
@@ -2,7 +2,7 @@
 
 Polygon masking is adapted from the `CGAL
 <https://doc.cgal.org/5.3.2/Polygon/index.html>`_ library. More information on
-point-in-polygon strategies can be found in Ref. :cite:p:`Schirra2008`.
+point-in-polygon strategies can be found in Ref. :cite:p:`schirra2008polygon`.
 
 .. currentmodule:: erlab.analysis.mask
 
diff --git a/src/erlab/analysis/utilities.py b/src/erlab/analysis/utilities.py
index 0a954953..b4076c62 100644
--- a/src/erlab/analysis/utilities.py
+++ b/src/erlab/analysis/utilities.py
@@ -29,12 +29,12 @@ def shift(
     darr
         The array to shift.
     shift
-        The amount of shift to be applied along the specified dimension. If `shift` is a
-        DataArray, different shifts can be applied to different coordinates. The
-        dimensions of `shift` must be a subset of the dimensions of `darr`. For more
-        information, see the note below. If `shift` is a `float`, the same shift is
-        applied to all values along dimension `along`. This is equivalent to providing a
-        0-dimensional DataArray.
+        The amount of shift to be applied along the specified dimension. If
+        :code:`shift` is a DataArray, different shifts can be applied to different
+        coordinates. The dimensions of :code:`shift` must be a subset of the dimensions
+        of `darr`. For more information, see the note below. If :code:`shift` is a
+        `float`, the same shift is applied to all values along dimension `along`. This
+        is equivalent to providing a 0-dimensional DataArray.
     along
         Name of the dimension along which the shift is applied.
     shift_coords
@@ -149,7 +149,7 @@ def correct_with_edge(
     **shift_kwargs,
 ):
     """
-    Corrects the given data array `darr` using the edge correction method.
+    Corrects the given data array `darr` with the given values or fit result.
 
     Parameters
     ----------
@@ -157,20 +157,24 @@ def correct_with_edge(
         The input data array to be corrected.
     modelresult
         The model result that contains the fermi edge information. It can be an instance
-        of `lmfit.model.ModelResult`, a numpy array, or a callable function that takes
-        an array of angles and returns the corresponding energy value.
+        of `lmfit.model.ModelResult`, a numpy array containing the edge position at each
+        angle, or a callable function that takes an array of angles and returns the
+        corresponding energy value.
     shift_coords
-        Whether to shift the coordinates of the data array. Defaults to True.
+        If `True`, the coordinates of the output data will be changed so that the output
+        contains all the values of the original data. If `False`, the coordinates and
+        shape of the original data will be retained, and only the data will be shifted.
+        Defaults to `False`.
     plot
-        Whether to plot the original and corrected data arrays. Defaults to False.
+        Whether to plot the original and corrected data arrays. Defaults to `False`.
     plot_kw
-        Additional keyword arguments for the plot. Defaults to None.
+        Additional keyword arguments for the plot. Defaults to `None`.
     **shift_kwargs
-        Additional keyword arguments to `shift`.
+        Additional keyword arguments to `erlab.analysis.utilities.shift`.
 
     Returns
     -------
-    xarray.DataArray
+    corrected : xarray.DataArray
         The edge corrected data.
     """
     if plot_kw is None:
diff --git a/src/erlab/interactive/__init__.py b/src/erlab/interactive/__init__.py
index 8a636429..9547b430 100644
--- a/src/erlab/interactive/__init__.py
+++ b/src/erlab/interactive/__init__.py
@@ -22,13 +22,9 @@
 
 """
 
-__all__ = ["goldtool", "itool", "ktool", "dtool"]
-
-import importlib
+__all__ = ["dtool", "goldtool", "itool", "ktool"]
 
+from erlab.interactive.derivative import dtool
 from erlab.interactive.fermiedge import goldtool
 from erlab.interactive.imagetool import itool
 from erlab.interactive.kspace import ktool
-
-if importlib.util.find_spec("arpes"):
-    from erlab.interactive.derivative import dtool
diff --git a/src/erlab/interactive/colors.py b/src/erlab/interactive/colors.py
index 30b4e32a..163b8ff2 100644
--- a/src/erlab/interactive/colors.py
+++ b/src/erlab/interactive/colors.py
@@ -274,7 +274,9 @@ def __init__(
         self.vb.setMouseEnabled(x=False, y=True)
 
         self._colorbar = pg.ImageItem(
-            np.linspace(0, 1, 4096).reshape((-1, 1)), axisOrder="row-major"
+            np.linspace(0, 1, 4096).reshape((-1, 1)),
+            axisOrder="row-major",
+            autoDownsample=False,
         )
         self.addItem(self._colorbar)
 
diff --git a/src/erlab/interactive/curvefittingtool.py b/src/erlab/interactive/curvefittingtool.py
index ebe4f184..41472368 100644
--- a/src/erlab/interactive/curvefittingtool.py
+++ b/src/erlab/interactive/curvefittingtool.py
@@ -1,6 +1,7 @@
 import copy
 import sys
 
+import lmfit
 import pyqtgraph as pg
 import xarray as xr
 from qtpy import QtCore, QtWidgets
@@ -161,7 +162,7 @@ def mouseDragEvent(self, ev):
 
         elif QtCore.Qt.MouseButton.RightButton in ev.buttons():
             if ev.isStart():
-                self._start_width = self.param_widget.widgets["width"].value()
+                self._start_width = self.param_widgÏet.widgets["width"].value()
             ev.accept()
 
             val = self.mapToParent(ev.buttonDownPos() - ev.pos()).y()
@@ -384,9 +385,13 @@ def _refresh_plot_peaks(self):
         self.modelplot.setData(x=self.xdata, y=model.eval(x=self.xdata, **params))
 
     def do_fit(self):
-        res = self.model.guess_fit(
-            self.data,
-            params=self.params_dict,
+        params = lmfit.create_params(**self.params_dict)
+        model = self.model
+        params = model.guess(self.data, self.data[self.data.dims[0]]).update(params)
+        res = self.model.fit(
+            self.ydata,
+            x=self.xdata,
+            params=params,
             method=self._params_init.values["Method"],
         )
         print(res.best_values)
@@ -614,9 +619,13 @@ def _refresh_plot_peaks(self):
         self.modelplot.setData(x=self.xdata, y=model.eval(x=self.xdata, **params))
 
     def do_fit(self):
-        res = self.model.guess_fit(
-            self.data,
-            params=self.params_dict,
+        params = lmfit.create_params(**self.params_dict)
+        model = self.model
+        params = model.guess(self.data, self.data[self.data.dims[0]]).update(params)
+        res = self.model.fit(
+            self.ydata,
+            x=self.xdata,
+            params=params,
             method=self._params_init.values["Method"],
         )
         print(res.best_values)
diff --git a/src/erlab/interactive/derivative.py b/src/erlab/interactive/derivative.py
index e7d3e18d..a3d8b3a4 100644
--- a/src/erlab/interactive/derivative.py
+++ b/src/erlab/interactive/derivative.py
@@ -1,177 +1,14 @@
+__all__ = ["dtool"]
+
 import sys
+
 import numpy as np
 import pyqtgraph as pg
-
-from arpes.analysis.derivative import curvature, minimum_gradient
-from erlab.interactive.utilities import parse_data, AnalysisWidgetBase
-
 from qtpy import QtCore, QtWidgets
 from scipy.ndimage import gaussian_filter, uniform_filter
 
-# def move_mean_centered(a, window, min_count=None, axis=-1):
-#     w = (window - 1) // 2
-#     shift = w + 1
-#     if min_count is None:
-#         min_count = w + 1
-#     pad_width = [(0, 0)] * a.ndim
-#     pad_width[axis] = (0, shift)
-#     a = np.pad(a, pad_width, constant_values=np.nan)
-#     val = bn.move_mean(a, window, min_count=min_count, axis=axis)
-#     return val[(slice(None),) * (axis % a.ndim) + (slice(w, -1),)]
-
-
-# def move_mean_centered_multiaxis(a, window_list, min_count_list=None):
-#     w_list = [(window - 1) // 2 for window in window_list]
-#     pad_width = [(0, 0)] * a.ndim
-#     slicer = [
-#         slice(None),
-#     ] * a.ndim
-#     if min_count_list is None:
-#         min_count_list = [w + 1 for w in w_list]
-#     for axis in range(a.ndim):
-#         pad_width[axis] = (0, w_list[axis] + 1)
-#         slicer[axis] = slice(w_list[axis], -1)
-#     a = np.pad(a, pad_width, constant_values=np.nan)
-#     val = move_mean(
-#         a,
-#         numba.typed.List(window_list),
-#         numba.typed.List(min_count_list),
-#     )
-#     return val[tuple(slicer)]
-
-
-# def move_mean(a, window, min_count):
-#     if a.ndim == 3:
-#         return move_mean3d(a, window, min_count)
-#     elif a.ndim == 2:
-#         return move_mean2d(a, window, min_count)
-#     else:
-#         raise NotImplementedError
-
-
-# pg.setConfigOption('useNumba', True)
-
-__all__ = ["dtool"]
-
-
-# def boxcarfilt2(array, window_list, min_count_list=None):
-#     return move_mean_centered_multiaxis(array, window_list, min_count_list)
-
-
-def laplacian_O3(f, *varargs):
-    """
-    Third order accurate, 5-point formula. Not really laplacian, but second derivative
-    along each axis. Sum the outputs to get laplacian.
-    """
-    N = len(f.shape)  # number of dimensions
-    n = len(varargs)
-    if n == 0:
-        dx = [1.0] * N
-    elif n == 1:
-        dx = [varargs[0]] * N
-    elif n == N:
-        dx = list(varargs)
-    else:
-        raise SyntaxError("invalid number of arguments")
-
-    # use central differences on interior and first differences on endpoints
-
-    # print dx
-    outvals = []
-
-    # create slice objects --- initially all are [:, :, ..., :]
-    slice0 = [slice(None)] * N
-    slice1 = [slice(None)] * N
-    slice2 = [slice(None)] * N
-    slice3 = [slice(None)] * N
-    slice4 = [slice(None)] * N
-
-    otype = f.dtype.char
-    if otype not in ["f", "d", "F", "D"]:
-        otype = "d"
-
-    for axis in range(N):
-        # select out appropriate parts for this dimension
-        out = np.zeros(f.shape, f.dtype.char)
-
-        # http://www.sitmo.com/eq/262 (3rd order accurate)
-        slice0[axis] = slice(2, -2)
-        slice1[axis] = slice(None, -4)
-        slice2[axis] = slice(1, -3)
-        slice3[axis] = slice(3, -1)
-        slice4[axis] = slice(4, None)
-        # 1D equivalent -- out[2:-2] = (-f[:-4] + 16*f[1:-3] + -30*f[2:-2] +
-        # 16*f[3:-1] - f[4:])/12.0
-        out[slice0] = (
-            -f[slice1]
-            + 16.0 * f[slice2]
-            - 30.0 * f[slice0]
-            + 16.0 * f[slice3]
-            - f[slice4]
-        ) / 12.0
-
-        # http://www.sitmo.com/eq/260 (2nd order accurate; there's also a 3rd
-        # order accurate that requires 5 points)
-        slice0[axis] = slice(None, 2)
-        slice1[axis] = slice(3, 5)
-        slice2[axis] = slice(2, 4)
-        slice3[axis] = slice(1, 3)
-        # 1D equivalent -- out[0:2] = 2*f[0:2] - 5*f[1:3] + 4*f[2:4] - f[3:5]
-        out[slice0] = 2.0 * f[slice0] - 5.0 * f[slice3] + 4.0 * f[slice2] - f[slice1]
-
-        slice0[axis] = slice(-2, None)
-        slice1[axis] = slice(-5, -3)
-        slice2[axis] = slice(-4, -2)
-        slice3[axis] = slice(-3, -1)
-        # 1D equivalent -- out[0:2] = 2*f[0:2] - 5*f[1:3] + 4*f[2:4] - f[3:5]
-        out[slice0] = 2.0 * f[slice0] - 5.0 * f[slice3] + 4.0 * f[slice2] - f[slice1]
-
-        # divide by step size
-        axis_dx = dx[axis]
-        outvals.append(out / (axis_dx * axis_dx))
-
-        # reset the slice object in this dimension to ":"
-        slice0[axis] = slice(None)
-        slice1[axis] = slice(None)
-        slice2[axis] = slice(None)
-        slice3[axis] = slice(None)
-        slice4[axis] = slice(None)
-
-    if N == 1:
-        return outvals[0]
-    else:
-        return outvals
-
-
-class NoiseToolNew(AnalysisWidgetBase):
-    def __init__(self, data, *args, **kwargs):
-        super().__init__(*args, **kwargs)
-        self.set_input(data)
-
-    # def get_smoothfunc(self, amount, num, method="boxcar"):
-    #     if method == "boxcar":
-    #         def func(data):
-    #             out = data.copy(deep=True)
-    #             vals = data.values
-    #             for _ in range(num):
-    #                 vals = uniform_filter(
-    #                     vals,
-    #                     amount,
-    #                     mode="constant",
-    #                     origin=(0, 0),
-    #                 )
-    #             out.values = vals
-    #             return out
-    #     elif method == "gaussian":
-    #         def func(data):
-    #             out = data.copy(deep=True)
-    #             vals = data.values
-    #             sigma = [(amount[0] - 1) / 3, (amount[1] - 1) / 3]
-    #             for _ in range(num):
-    #                 vals = gaussian_filter(vals, sigma=sigma)
-    #             out.values = vals
-    #             return out
-    #     return func
+from erlab.analysis.image import curvature, minimum_gradient
+from erlab.interactive.utilities import parse_data
 
 
 class DerivativeTool(QtWidgets.QMainWindow):
@@ -273,17 +110,26 @@ def __init__(self, data=None, *args, **kwargs):
         self._curv_group = QtWidgets.QGroupBox("Curvature")
         curv_layout = QtWidgets.QGridLayout(self._curv_group)
         self._curv_alpha_spin = QtWidgets.QDoubleSpinBox(self._curv_group)
-        self._curv_alpha_spin.setRange(-30, 30)
-        self._curv_alpha_spin.setValue(0)
-        self._curv_alpha_spin.setSingleStep(0.05)
+        self._curv_alpha_spin.setRange(0, 100)
+        self._curv_alpha_spin.setDecimals(4)
+        self._curv_alpha_spin.setValue(1.0)
+        self._curv_alpha_spin.setSingleStep(0.001)
+        self._curv_factor_spin = QtWidgets.QDoubleSpinBox(self._curv_group)
+        self._curv_factor_spin.setRange(-1, 1)
+        self._curv_factor_spin.setDecimals(4)
+        self._curv_factor_spin.setValue(1.0)
+        self._curv_factor_spin.setSingleStep(0.0001)
         curv_layout.addWidget(self._curv_alpha_spin, 0, 0)
+        curv_layout.addWidget(self._curv_factor_spin, 1, 0)
 
         self._mingrad_group = QtWidgets.QGroupBox("Minimum Gradient")
         mingrad_layout = QtWidgets.QGridLayout(self._mingrad_group)
-        self._mingrad_delta_spin = QtWidgets.QSpinBox(self._mingrad_group)
-        self._mingrad_delta_spin.setRange(1, 100)
-        self._mingrad_delta_spin.setValue(1)
-        mingrad_layout.addWidget(self._mingrad_delta_spin, 0, 0)
+        self._mingrad_btn = QtWidgets.QPushButton("Calculate")
+        # self._mingrad_delta_spin = QtWidgets.QSpinBox(self._mingrad_group)
+        # self._mingrad_delta_spin.setRange(1, 100)
+        # self._mingrad_delta_spin.setValue(1)
+        # mingrad_layout.addWidget(self._mingrad_delta_spin, 0, 0)
+        mingrad_layout.addWidget(self._mingrad_btn, 0, 0)
 
         self._deriv_tab = QtWidgets.QWidget()
         deriv_content = QtWidgets.QGridLayout(self._deriv_tab)
@@ -331,12 +177,12 @@ def __init__(self, data=None, *args, **kwargs):
         self._smooth2_y_spin.valueChanged.connect(lambda: self.deriv_data())
         self._smooth2_n_spin.valueChanged.connect(self.set_smooth2_num)
         self._deriv_axis_combo.currentTextChanged.connect(self.deriv_data)
-        self._curv_alpha_spin.valueChanged.connect(
-            lambda beta: self.curv_data(beta=beta)
-        )
-        self._mingrad_delta_spin.valueChanged.connect(
-            lambda delta: self.mingrad_data(delta=delta)
-        )
+        self._curv_alpha_spin.valueChanged.connect(self.curv_data)
+        self._curv_factor_spin.valueChanged.connect(self.curv_data)
+        self._mingrad_btn.clicked.connect(lambda: self.mingrad_data())
+        # self._mingrad_delta_spin.valueChanged.connect(
+        #     lambda delta: self.mingrad_data(delta=delta)
+        # )
         for s in self._color_range_spin:
             s.valueChanged.connect(self._result_update_lims)
 
@@ -385,23 +231,27 @@ def set_smooth2_num(self, n):
         self.smooth2_num = n
         self.deriv_data()
 
-    def mingrad_data(self, delta=1):
+    # def mingrad_data(self, delta=1):
+    def mingrad_data(self):
         self.use_curv = False
         self.use_mingrad = True
         self.use_deriv = False
         # self.data_s_d = self.data_s.copy(deep=True)
-        self.data_s_d = -minimum_gradient(self.data_s, delta=delta)
+        # self.data_s_d = -minimum_gradient(self.data_s, delta=delta)
+        self.data_s_d = -minimum_gradient(self.data_s)
         self._result_update_lims()
 
-    def curv_data(self, alpha=1, beta=None):
+    def curv_data(self):
         self.use_curv = True
         self.use_deriv = False
         self.use_mingrad = False
-        if beta is None:
-            beta = self._curv_alpha_spin.value()
-        self.data_s_d = self.data_s.copy(deep=True)
-        self.data_s_d.values = curvature(
-            self.data_s / self.data_s.max(), alpha=alpha, beta=beta, values=True
+
+        # self.data_s_d = self.data_s.copy(deep=True)
+        self.data_s_d = curvature(
+            # self.data_s / self.data_s.max(), alpha=alpha, beta=beta, values=True
+            self.data_s,
+            a0=self._curv_alpha_spin.value(),
+            factor=self._curv_factor_spin.value(),
         )
         # m = self.data_s_d.max()
         # self.data_s_d /= m
@@ -506,7 +356,7 @@ def dtool(data, *args, **kwargs):
 if __name__ == "__main__":
     import erlab.io
 
-    gkmk_cvs = erlab.io.load_als_bl4(
+    gkmk_cvs = erlab.io.merlin.load(
         "/Users/khan/Documents/ERLab/CsV3Sb5/2021_Dec_ALS_CV3Sb5/211217 ALS BL4/csvtisb1/f_003.pxt"
     )
     alpha_new = np.linspace(gkmk_cvs.alpha[0], gkmk_cvs.alpha[-1], 1000)
diff --git a/src/erlab/interactive/imagetool/core.py b/src/erlab/interactive/imagetool/core.py
index 3dadf17e..0703fa34 100644
--- a/src/erlab/interactive/imagetool/core.py
+++ b/src/erlab/interactive/imagetool/core.py
@@ -6,7 +6,6 @@
 
 import collections
 import functools
-import gc
 import inspect
 import os
 import time
@@ -25,7 +24,7 @@
     pg_colormap_powernorm,
 )
 from erlab.interactive.imagetool.slicer import ArraySlicer
-from erlab.interactive.utilities import BetterAxisItem
+from erlab.interactive.utilities import BetterAxisItem, copy_to_clipboard
 
 if TYPE_CHECKING:
     from collections.abc import Callable, Iterable, Sequence
@@ -452,7 +451,6 @@ def on_close(self):
         if hasattr(self, "_data"):
             self._data.close()
             del self._data
-        gc.collect()
 
     def connect_axes_signals(self):
         for ax in self.axes:
@@ -1152,6 +1150,9 @@ def __init__(
         save_action = self.vb.menu.addAction("Save data as HDF5")
         save_action.triggered.connect(lambda: self.save_current_data())
 
+        copy_code_action = self.vb.menu.addAction("Copy selection code")
+        copy_code_action.triggered.connect(self.copy_selection_code)
+
         for i in (0, 1):
             self.getViewBoxMenu().ctrl[i].linkCombo.setVisible(False)
             self.getViewBoxMenu().ctrl[i].label.setVisible(False)
@@ -1308,7 +1309,7 @@ def add_cursor(self, update=True):
             item = self.image_cls(
                 self,
                 cursor=new_cursor,
-                autoDownsample=True,
+                autoDownsample=False,
                 axisOrder="row-major",
                 **self._item_kw,
             )
@@ -1509,6 +1510,14 @@ def save_current_data(self, fileName=None):
             fileName,
         )
 
+    @QtCore.Slot()
+    def copy_selection_code(self):
+        copy_to_clipboard(
+            self.array_slicer.qsel_code(
+                self.slicer_area.current_cursor, self.display_axis
+            )
+        )
+
     @property
     def display_axis(self) -> tuple[int, ...]:
         return self._display_axis
diff --git a/src/erlab/interactive/imagetool/fastbinning.py b/src/erlab/interactive/imagetool/fastbinning.py
index 51e6e0eb..6fadf7ad 100644
--- a/src/erlab/interactive/imagetool/fastbinning.py
+++ b/src/erlab/interactive/imagetool/fastbinning.py
@@ -6,13 +6,51 @@
 
 __all__ = ["fast_nanmean"]
 
-import numpy as np
-import numpy.typing as npt
 from collections.abc import Iterable
+
 import numba
-import numba.typed
 import numba.core.registry
+import numba.typed
 import numbagg
+import numpy as np
+import numpy.typing as npt
+
+_SIG_2_1 = [
+    numba.types.Array(numba.float64, 1, "C")(numba.types.Array(numba.float32, 2, "C")),
+    numba.types.Array(numba.float64, 1, "C")(numba.types.Array(numba.float32, 2, "A")),
+    numba.types.Array(numba.float64, 1, "C")(numba.types.Array(numba.float64, 2, "C")),
+    numba.types.Array(numba.float64, 1, "C")(numba.types.Array(numba.float64, 2, "A")),
+]
+_SIG_3_1 = [
+    numba.types.Array(numba.float64, 1, "C")(numba.types.Array(numba.float32, 3, "C")),
+    numba.types.Array(numba.float64, 1, "C")(numba.types.Array(numba.float32, 3, "A")),
+    numba.types.Array(numba.float64, 1, "C")(numba.types.Array(numba.float64, 3, "C")),
+    numba.types.Array(numba.float64, 1, "C")(numba.types.Array(numba.float64, 3, "A")),
+]
+_SIG_3_2 = [
+    numba.types.Array(numba.float64, 2, "C")(numba.types.Array(numba.float32, 3, "C")),
+    numba.types.Array(numba.float64, 2, "C")(numba.types.Array(numba.float32, 3, "A")),
+    numba.types.Array(numba.float64, 2, "C")(numba.types.Array(numba.float64, 3, "C")),
+    numba.types.Array(numba.float64, 2, "C")(numba.types.Array(numba.float64, 3, "A")),
+]
+_SIG_4_1 = [
+    numba.types.Array(numba.float64, 1, "C")(numba.types.Array(numba.float32, 4, "C")),
+    numba.types.Array(numba.float64, 1, "C")(numba.types.Array(numba.float32, 4, "A")),
+    numba.types.Array(numba.float64, 1, "C")(numba.types.Array(numba.float64, 4, "C")),
+    numba.types.Array(numba.float64, 1, "C")(numba.types.Array(numba.float64, 4, "A")),
+]
+_SIG_4_2 = [
+    numba.types.Array(numba.float64, 2, "C")(numba.types.Array(numba.float32, 4, "C")),
+    numba.types.Array(numba.float64, 2, "C")(numba.types.Array(numba.float32, 4, "A")),
+    numba.types.Array(numba.float64, 2, "C")(numba.types.Array(numba.float64, 4, "C")),
+    numba.types.Array(numba.float64, 2, "C")(numba.types.Array(numba.float64, 4, "A")),
+]
+_SIG_4_3 = [
+    numba.types.Array(numba.float64, 3, "C")(numba.types.Array(numba.float32, 4, "C")),
+    numba.types.Array(numba.float64, 3, "C")(numba.types.Array(numba.float32, 4, "A")),
+    numba.types.Array(numba.float64, 3, "C")(numba.types.Array(numba.float64, 4, "C")),
+    numba.types.Array(numba.float64, 3, "C")(numba.types.Array(numba.float64, 4, "A")),
+]
 
 
 @numba.njit(cache=True)
@@ -20,7 +58,7 @@ def _nanmean_all(a: npt.NDArray[np.float32 | np.float64]) -> np.float64:
     return np.nanmean(a)
 
 
-@numba.njit(cache=True, parallel=True)
+@numba.njit(_SIG_2_1, cache=True, parallel=True)
 def _nanmean_2_0(a: npt.NDArray[np.float32 | np.float64]) -> npt.NDArray[np.float64]:
     n = a.shape[1]
     output = np.empty(n, dtype=np.float64)
@@ -29,7 +67,7 @@ def _nanmean_2_0(a: npt.NDArray[np.float32 | np.float64]) -> npt.NDArray[np.floa
     return output
 
 
-@numba.njit(cache=True, parallel=True)
+@numba.njit(_SIG_2_1, cache=True, parallel=True)
 def _nanmean_2_1(a: npt.NDArray[np.float32 | np.float64]) -> npt.NDArray[np.float64]:
     n = a.shape[0]
     output = np.empty(n, dtype=np.float64)
@@ -38,7 +76,7 @@ def _nanmean_2_1(a: npt.NDArray[np.float32 | np.float64]) -> npt.NDArray[np.floa
     return output
 
 
-@numba.njit(cache=True, parallel=True)
+@numba.njit(_SIG_3_2, cache=True, parallel=True)
 def _nanmean_3_0(a: npt.NDArray[np.float32 | np.float64]) -> npt.NDArray[np.float64]:
     _, n0, n1 = a.shape
     output = np.empty((n0, n1), dtype=np.float64)
@@ -48,7 +86,7 @@ def _nanmean_3_0(a: npt.NDArray[np.float32 | np.float64]) -> npt.NDArray[np.floa
     return output
 
 
-@numba.njit(cache=True, parallel=True)
+@numba.njit(_SIG_3_2, cache=True, parallel=True)
 def _nanmean_3_1(a: npt.NDArray[np.float32 | np.float64]) -> npt.NDArray[np.float64]:
     n0, _, n1 = a.shape
     output = np.empty((n0, n1), dtype=np.float64)
@@ -58,7 +96,7 @@ def _nanmean_3_1(a: npt.NDArray[np.float32 | np.float64]) -> npt.NDArray[np.floa
     return output
 
 
-@numba.njit(cache=True, parallel=True)
+@numba.njit(_SIG_3_2, cache=True, parallel=True)
 def _nanmean_3_2(a: npt.NDArray[np.float32 | np.float64]) -> npt.NDArray[np.float64]:
     n0, n1, _ = a.shape
     output = np.empty((n0, n1), dtype=np.float64)
@@ -68,7 +106,7 @@ def _nanmean_3_2(a: npt.NDArray[np.float32 | np.float64]) -> npt.NDArray[np.floa
     return output
 
 
-@numba.njit(cache=True, parallel=True)
+@numba.njit(_SIG_3_1, cache=True, parallel=True)
 def _nanmean_3_01(a: npt.NDArray[np.float32 | np.float64]) -> npt.NDArray[np.float64]:
     n = a.shape[2]
     output = np.empty(n, dtype=np.float64)
@@ -77,7 +115,7 @@ def _nanmean_3_01(a: npt.NDArray[np.float32 | np.float64]) -> npt.NDArray[np.flo
     return output
 
 
-@numba.njit(cache=True, parallel=True)
+@numba.njit(_SIG_3_1, cache=True, parallel=True)
 def _nanmean_3_02(a: npt.NDArray[np.float32 | np.float64]) -> npt.NDArray[np.float64]:
     n = a.shape[1]
     output = np.empty(n, dtype=np.float64)
@@ -86,7 +124,7 @@ def _nanmean_3_02(a: npt.NDArray[np.float32 | np.float64]) -> npt.NDArray[np.flo
     return output
 
 
-@numba.njit(cache=True, parallel=True)
+@numba.njit(_SIG_3_1, cache=True, parallel=True)
 def _nanmean_3_12(a: npt.NDArray[np.float32 | np.float64]) -> npt.NDArray[np.float64]:
     n = a.shape[0]
     output = np.empty(n, dtype=np.float64)
@@ -95,7 +133,7 @@ def _nanmean_3_12(a: npt.NDArray[np.float32 | np.float64]) -> npt.NDArray[np.flo
     return output
 
 
-@numba.njit(cache=True, parallel=True)
+@numba.njit(_SIG_4_3, cache=True, parallel=True)
 def _nanmean_4_0(a: npt.NDArray[np.float32 | np.float64]) -> npt.NDArray[np.float64]:
     _, n0, n1, n2 = a.shape
     output = np.empty((n0, n1, n2), dtype=np.float64)
@@ -106,7 +144,7 @@ def _nanmean_4_0(a: npt.NDArray[np.float32 | np.float64]) -> npt.NDArray[np.floa
     return output
 
 
-@numba.njit(cache=True, parallel=True)
+@numba.njit(_SIG_4_3, cache=True, parallel=True)
 def _nanmean_4_1(a: npt.NDArray[np.float32 | np.float64]) -> npt.NDArray[np.float64]:
     n0, _, n1, n2 = a.shape
     output = np.empty((n0, n1, n2), dtype=np.float64)
@@ -117,7 +155,7 @@ def _nanmean_4_1(a: npt.NDArray[np.float32 | np.float64]) -> npt.NDArray[np.floa
     return output
 
 
-@numba.njit(cache=True, parallel=True)
+@numba.njit(_SIG_4_3, cache=True, parallel=True)
 def _nanmean_4_2(a: npt.NDArray[np.float32 | np.float64]) -> npt.NDArray[np.float64]:
     n0, n1, _, n2 = a.shape
     output = np.empty((n0, n1, n2), dtype=np.float64)
@@ -128,7 +166,7 @@ def _nanmean_4_2(a: npt.NDArray[np.float32 | np.float64]) -> npt.NDArray[np.floa
     return output
 
 
-@numba.njit(cache=True, parallel=True)
+@numba.njit(_SIG_4_3, cache=True, parallel=True)
 def _nanmean_4_3(a: npt.NDArray[np.float32 | np.float64]) -> npt.NDArray[np.float64]:
     n0, n1, n2, _ = a.shape
     output = np.empty((n0, n1, n2), dtype=np.float64)
@@ -139,7 +177,7 @@ def _nanmean_4_3(a: npt.NDArray[np.float32 | np.float64]) -> npt.NDArray[np.floa
     return output
 
 
-@numba.njit(cache=True, parallel=True)
+@numba.njit(_SIG_4_2, cache=True, parallel=True)
 def _nanmean_4_01(a: npt.NDArray[np.float32 | np.float64]) -> npt.NDArray[np.float64]:
     _, _, n0, n1 = a.shape
     output = np.empty((n0, n1), dtype=np.float64)
@@ -149,7 +187,7 @@ def _nanmean_4_01(a: npt.NDArray[np.float32 | np.float64]) -> npt.NDArray[np.flo
     return output
 
 
-@numba.njit(cache=True, parallel=True)
+@numba.njit(_SIG_4_2, cache=True, parallel=True)
 def _nanmean_4_02(a: npt.NDArray[np.float32 | np.float64]) -> npt.NDArray[np.float64]:
     _, n0, _, n1 = a.shape
     output = np.empty((n0, n1), dtype=np.float64)
@@ -159,7 +197,7 @@ def _nanmean_4_02(a: npt.NDArray[np.float32 | np.float64]) -> npt.NDArray[np.flo
     return output
 
 
-@numba.njit(cache=True, parallel=True)
+@numba.njit(_SIG_4_2, cache=True, parallel=True)
 def _nanmean_4_03(a: npt.NDArray[np.float32 | np.float64]) -> npt.NDArray[np.float64]:
     _, n0, n1, _ = a.shape
     output = np.empty((n0, n1), dtype=np.float64)
@@ -169,7 +207,7 @@ def _nanmean_4_03(a: npt.NDArray[np.float32 | np.float64]) -> npt.NDArray[np.flo
     return output
 
 
-@numba.njit(cache=True, parallel=True)
+@numba.njit(_SIG_4_2, cache=True, parallel=True)
 def _nanmean_4_12(a: npt.NDArray[np.float32 | np.float64]) -> npt.NDArray[np.float64]:
     n0, _, _, n1 = a.shape
     output = np.empty((n0, n1), dtype=np.float64)
@@ -179,7 +217,7 @@ def _nanmean_4_12(a: npt.NDArray[np.float32 | np.float64]) -> npt.NDArray[np.flo
     return output
 
 
-@numba.njit(cache=True, parallel=True)
+@numba.njit(_SIG_4_2, cache=True, parallel=True)
 def _nanmean_4_13(a: npt.NDArray[np.float32 | np.float64]) -> npt.NDArray[np.float64]:
     n0, _, n1, _ = a.shape
     output = np.empty((n0, n1), dtype=np.float64)
@@ -189,7 +227,7 @@ def _nanmean_4_13(a: npt.NDArray[np.float32 | np.float64]) -> npt.NDArray[np.flo
     return output
 
 
-@numba.njit(cache=True, parallel=True)
+@numba.njit(_SIG_4_2, cache=True, parallel=True)
 def _nanmean_4_23(a: npt.NDArray[np.float32 | np.float64]) -> npt.NDArray[np.float64]:
     n0, n1, _, _ = a.shape
     output = np.empty((n0, n1), dtype=np.float64)
@@ -199,7 +237,7 @@ def _nanmean_4_23(a: npt.NDArray[np.float32 | np.float64]) -> npt.NDArray[np.flo
     return output
 
 
-@numba.njit(cache=True, parallel=True)
+@numba.njit(_SIG_4_1, cache=True, parallel=True)
 def _nanmean_4_012(a: npt.NDArray[np.float32 | np.float64]) -> npt.NDArray[np.float64]:
     n = a.shape[3]
     output = np.empty(n, dtype=np.float64)
@@ -208,7 +246,7 @@ def _nanmean_4_012(a: npt.NDArray[np.float32 | np.float64]) -> npt.NDArray[np.fl
     return output
 
 
-@numba.njit(cache=True, parallel=True)
+@numba.njit(_SIG_4_1, cache=True, parallel=True)
 def _nanmean_4_013(a: npt.NDArray[np.float32 | np.float64]) -> npt.NDArray[np.float64]:
     n = a.shape[2]
     output = np.empty(n, dtype=np.float64)
@@ -217,7 +255,7 @@ def _nanmean_4_013(a: npt.NDArray[np.float32 | np.float64]) -> npt.NDArray[np.fl
     return output
 
 
-@numba.njit(cache=True, parallel=True)
+@numba.njit(_SIG_4_1, cache=True, parallel=True)
 def _nanmean_4_023(a: npt.NDArray[np.float32 | np.float64]) -> npt.NDArray[np.float64]:
     n = a.shape[1]
     output = np.empty(n, dtype=np.float64)
@@ -226,7 +264,7 @@ def _nanmean_4_023(a: npt.NDArray[np.float32 | np.float64]) -> npt.NDArray[np.fl
     return output
 
 
-@numba.njit(cache=True, parallel=True)
+@numba.njit(_SIG_4_1, cache=True, parallel=True)
 def _nanmean_4_123(a: npt.NDArray[np.float32 | np.float64]) -> npt.NDArray[np.float64]:
     n = a.shape[0]
     output = np.empty(n, dtype=np.float64)
diff --git a/src/erlab/interactive/imagetool/slicer.py b/src/erlab/interactive/imagetool/slicer.py
index 808e6dcb..de3098d2 100644
--- a/src/erlab/interactive/imagetool/slicer.py
+++ b/src/erlab/interactive/imagetool/slicer.py
@@ -576,18 +576,81 @@ def index_of_value(
         else:
             return _index_of_value_nonuniform(self.coords[axis], value)
 
-    def gen_selection_code(self, cursor: int, disp: Sequence[int]) -> str:
+    def isel_args(
+        self, cursor: int, disp: Sequence[int], int_if_one: bool = False
+    ) -> dict[str, slice | int]:
         axis = sorted(set(range(self._obj.ndim)) - set(disp))
-        slice_dict = {
-            self._obj.dims[ax]: self._bin_slice(cursor, ax, int_if_one=True)
-            for ax in axis
+        return {
+            self._obj.dims[ax]: self._bin_slice(cursor, ax, int_if_one) for ax in axis
         }
-        return f".isel(**{slice_dict!s}).squeeze()"
+
+    def qsel_args(self, cursor: int, disp: Sequence[int]) -> dict:
+        out: dict[str, float] = {}
+        binned = self.get_binned(cursor)
+
+        for dim, selector in self.isel_args(cursor, disp, int_if_one=True).items():
+            inc = self.incs[self._obj.dims.index(dim)]
+            order = int(-np.floor(np.log10(inc)) + 1)
+
+            if binned[self._obj.dims.index(dim)]:
+                coord = self._obj[dim][selector].values
+
+                out[dim] = np.round(coord.mean(), order)
+                width = np.round(abs(coord[-1] - coord[0]) + inc, order)
+
+                if not np.allclose(
+                    self._obj[dim]
+                    .sel({dim: slice(out[dim] - width / 2, out[dim] + width / 2)})
+                    .values,
+                    coord,
+                ):
+                    raise ValueError(
+                        "Bin does not contain the same values as the original data."
+                    )
+
+                out[dim + "_width"] = width
+
+            else:
+                out[dim] = np.round(self._obj[dim].values[selector], order)
+
+        return out
+
+    def qsel_code(self, cursor: int, disp: Sequence[int]) -> str:
+        if self._nonuniform_axes:
+            # Has non-uniform axes, fallback to isel
+            return self.isel_code(cursor, disp)
+
+        try:
+            qsel_kw = self.qsel_args(cursor, disp)
+        except ValueError:
+            return self.isel_code(cursor, disp)
+
+        dict_repr: str = ""
+        for k, v in qsel_kw.items():
+            dict_repr += f"{k}={v!s}, "
+        dict_repr = dict_repr.rstrip(", ")
+        return f".qsel({dict_repr})"
+
+    def isel_code(self, cursor: int, disp: Sequence[int]) -> str:
+        dict_repr: str = ""
+        for k, v in self.isel_args(cursor, disp, int_if_one=True).items():
+            dict_repr += f"{k}={v!s}, "
+        dict_repr = dict_repr.rstrip(", ")
+        return f".isel({dict_repr})"
 
     def xslice(self, cursor: int, disp: Sequence[int]) -> xr.DataArray:
-        axis = sorted(set(range(self._obj.ndim)) - set(disp))
-        slices = {self._obj.dims[ax]: self._bin_slice(cursor, ax) for ax in axis}
-        return self._obj.isel(**slices).squeeze()
+        isel_kw: dict[str, slice] = self.isel_args(cursor, disp, int_if_one=False)
+        binned_coord_average: dict[str, xr.DataArray] = {
+            k: self._obj[k][isel_kw[k]].mean()
+            for k, v in zip(self._obj.dims, self.get_binned(cursor))
+            if v
+        }
+        return (
+            self._obj.isel(**isel_kw)
+            .squeeze()
+            .mean(binned_coord_average.keys())
+            .assign_coords(**binned_coord_average)
+        )
 
     @QtCore.Slot(int, tuple, result=np.ndarray)
     def slice_with_coord(
diff --git a/src/erlab/plotting/general.py b/src/erlab/plotting/general.py
index a3663ccd..244055a0 100644
--- a/src/erlab/plotting/general.py
+++ b/src/erlab/plotting/general.py
@@ -416,6 +416,7 @@ def plot_array_2d(
     cax: matplotlib.axes.Axes | None = None,
     colorbar_kw: dict | None = None,
     imshow_kw: dict | None = None,
+    N: int = 256,
     **indexers_kwargs: dict,
 ):
     if lnorm is None:
@@ -456,7 +457,13 @@ def plot_array_2d(
     carr = carr.sel(**sel_kw)
 
     cmap_img, img = gen_2d_colormap(
-        larr.values, carr.values, cmap, lnorm=lnorm, cnorm=cnorm, background=background
+        larr.values,
+        carr.values,
+        cmap,
+        lnorm=lnorm,
+        cnorm=cnorm,
+        background=background,
+        N=N,
     )
 
     if colorbar:


Test patch:
diff --git a/tests/test_fit_functions_dynamic.py b/tests/test_fit_functions_dynamic.py
new file mode 100644
index 00000000..1cc4aae1
--- /dev/null
+++ b/tests/test_fit_functions_dynamic.py
@@ -0,0 +1,136 @@
+import numpy as np
+import pytest
+import xarray as xr
+
+from erlab.analysis.fit.functions.dynamic import (
+    FermiEdge2dFunction,
+    MultiPeakFunction,
+    PolynomialFunction,
+)
+from erlab.analysis.fit.functions.general import gaussian_wh, lorentzian_wh
+
+RAND_STATE = np.random.RandomState(1)
+
+
+def test_poly_func_call():
+    # Test case 1: Evaluate polynomial function with numpy array input
+    x = np.arange(5, dtype=np.float64)
+    coeffs = RAND_STATE.randn(3)
+    expected_result = np.polyval(np.asarray(list(reversed(coeffs))), x)
+    assert np.allclose(
+        PolynomialFunction(degree=2)(x, *coeffs),
+        expected_result,
+    )
+
+    # Test case 2: Evaluate polynomial function with xarray input
+    x = xr.DataArray(np.arange(5, dtype=np.float64), dims="x")
+    coeffs = RAND_STATE.randn(3)
+    expected_result = np.polyval(np.asarray(list(reversed(coeffs))), x)
+    result = PolynomialFunction(degree=2)(x, *coeffs)
+    assert isinstance(result, xr.DataArray)
+    assert np.allclose(result, expected_result)
+
+    # Test case 3: Evaluate polynomial function with parameters
+    x = np.arange(5, dtype=np.float64)
+    coeffs = RAND_STATE.randn(3)
+    expected_result = np.polyval(np.asarray(list(reversed(coeffs))), x)
+    params = dict(zip([f"c{i}" for i in range(3)], coeffs))
+    result = PolynomialFunction(degree=2)(x, **params)
+    assert np.allclose(result, expected_result)
+
+
+def test_multi_peak_function_call():
+    # Test case 1: Evaluate multi-peak function with numpy array input
+    x = np.linspace(-5, 5, 20, dtype=np.float64)
+    npeaks = 2
+    peak_shapes = ["lorentzian", "gaussian"]
+    fd = False
+    convolve = False
+    params = {
+        "p0_center": 0.0,
+        "p0_height": 1.0,
+        "p0_width": 0.5,
+        "p1_center": 2.0,
+        "p1_height": 0.5,
+        "p1_width": 0.2,
+        "lin_bkg": 0.1,
+        "const_bkg": 0.2,
+    }
+    expected_result = (
+        lorentzian_wh(x, params["p0_center"], params["p0_width"], params["p0_height"])
+        + gaussian_wh(x, params["p1_center"], params["p1_width"], params["p1_height"])
+        + params["lin_bkg"] * x
+        + params["const_bkg"]
+    )
+    assert np.allclose(
+        MultiPeakFunction(npeaks, peak_shapes, fd, convolve)(x, **params),
+        expected_result,
+    )
+    # Test case 2: With different peak shape signature
+    assert np.allclose(
+        MultiPeakFunction(npeaks, "l g", fd, convolve)(x, **params),
+        expected_result,
+    )
+
+    # Test case 3: Evaluate multi-peak function with xarray input
+    x = xr.DataArray(np.linspace(-5, 5, 20, dtype=np.float64), dims="x")
+    result_xr = MultiPeakFunction(npeaks, peak_shapes, fd, convolve)(x, **params)
+    assert isinstance(result_xr, xr.DataArray)
+    assert np.allclose(result_xr, expected_result)
+
+    # Test case 4: Evaluate multi-peak function with different peak shape signatures
+    x = np.linspace(-5, 5, 20, dtype=np.float64)
+    npeaks = 3
+    peak_shapes = ["gaussian", "lorentzian", "gaussian"]
+    fd = True
+    convolve = True
+    params = {
+        "p0_center": 1.0,
+        "p0_height": 0.5,
+        "p0_width": 0.2,
+        "p1_center": 2.0,
+        "p1_height": 1.0,
+        "p1_width": 0.5,
+        "p2_center": 3.0,
+        "p2_height": 0.8,
+        "p2_width": 0.3,
+        "lin_bkg": 0.2,
+        "const_bkg": 0.1,
+        "efermi": 0.0,
+        "temp": 30,
+        "offset": 0.5,
+        "resolution": 0.1,
+    }
+
+    expected_result = MultiPeakFunction(npeaks, "lorentzian", fd, convolve)(x, **params)
+
+    for peak_shapes in ("l l l", ["lorentzian"] * 3):
+        assert np.allclose(
+            MultiPeakFunction(npeaks, peak_shapes, fd, convolve)(x, **params),
+            expected_result,
+        )
+
+
+@pytest.mark.filterwarnings("ignore:overflow encountered in exp.*:RuntimeWarning")
+def test_fermi_edge_2d_function_call():
+    # Difficult to test function values, just check if equal for numpy and xarray
+    eV = np.linspace(-5, 5, 25, dtype=np.float64)
+    alpha = RAND_STATE.randn(20)
+    alpha.sort()
+    params = {
+        "c0": 1.0,
+        "c1": 0.5,
+        "lin_bkg": 0.1,
+        "const_bkg": 0.2,
+        "offset": 0.0,
+        "temp": 30.0,
+        "resolution": 0.02,
+    }
+    out_np = FermiEdge2dFunction()(eV, alpha, **params)
+
+    eV = xr.DataArray(eV, dims="eV")
+    alpha = xr.DataArray(alpha, dims="alpha")
+    out_xr = FermiEdge2dFunction()(eV, alpha, **params)
+
+    assert isinstance(out_xr, xr.DataArray)
+    assert np.allclose(out_np.reshape((25, 20)), out_xr)
diff --git a/tests/test_fit_functions_general.py b/tests/test_fit_functions_general.py
new file mode 100644
index 00000000..23845795
--- /dev/null
+++ b/tests/test_fit_functions_general.py
@@ -0,0 +1,175 @@
+import numpy as np
+import scipy.special
+from erlab.analysis.fit.functions.general import (
+    _gen_kernel,
+    do_convolve,
+    fermi_dirac,
+    fermi_dirac_linbkg,
+    fermi_dirac_linbkg_broad,
+    gaussian_wh,
+    lorentzian_wh,
+    step_broad,
+    step_linbkg_broad,
+)
+import pytest
+
+KB_EV = 8.617333262145179e-5
+
+
+def test_gen_kernel():
+    # Define test input values
+    x = np.array([0, 1, 2, 3, 4])
+    resolution = 2.0
+    pad = 5
+
+    # Call the function to generate the kernel
+    extended, gauss = _gen_kernel(x, resolution, pad)
+
+    # Define the expected results
+    expected_extended = np.arange(-4, 9, 1, dtype=np.float64)
+    expected_gauss = np.array(
+        [
+            7.16733764e-06,
+            9.17419217e-04,
+            2.93574150e-02,
+            2.34859320e-01,
+            4.69718639e-01,
+            2.34859320e-01,
+            2.93574150e-02,
+            9.17419217e-04,
+            7.16733764e-06,
+        ]
+    )
+
+    # Check if the generated extended and gauss arrays match the expected results
+    assert np.allclose(extended, expected_extended)
+    assert np.allclose(gauss, expected_gauss)
+
+
+def test_do_convolve():
+    # Test case 1: Convolve a linear function with a Gaussian kernel
+    x = np.linspace(0, 10, 5)
+
+    def testfunc(x):
+        return x
+
+    resolution = 1.0
+    expected = np.array(
+        [-6.61744490e-24, 5.87148334e00, 1.17429667e01, 1.76144500e01, 2.34859334e01]
+    )
+    result = do_convolve(x, testfunc, resolution)
+    assert np.allclose(result, expected)
+
+    # Test case 2: Convolve a sine function with a Gaussian kernel
+    x = np.linspace(0, 2 * np.pi, 5)
+
+    def testfunc(x):
+        return np.sin(x)
+
+    resolution = 0.5
+    expected = np.array(
+        [0.00000000e00, 2.95132925e00, 3.61433592e-16, -2.95132925e00, -7.22867185e-16]
+    )
+    result = do_convolve(x, testfunc, resolution)
+    assert np.allclose(result, expected)
+
+
+def test_gaussian_wh():
+    x = np.linspace(0, 10, 100)
+    center = 5.0
+    width = 2.0
+    height = 1.0
+    expected_result = height * np.exp(
+        -16 * np.log(2) * (1.0 * x - center) ** 2 / max(1e-10, width**2)
+    )
+    assert np.allclose(gaussian_wh(x, center, width, height), expected_result)
+
+
+def test_lorentzian_wh():
+    x = np.linspace(0, 10, 100)
+    center = 5.0
+    width = 2.0
+    height = 1.0
+    expected_result = height / (1 + 4 * ((1.0 * x - center) / max(1.0e-15, width)) ** 2)
+    assert np.allclose(lorentzian_wh(x, center, width, height), expected_result)
+
+
+@pytest.mark.filterwarnings("ignore:overflow encountered in exp.*:RuntimeWarning")
+def test_fermi_dirac():
+    x = np.linspace(0, 10, 100)
+    center = 5.0
+    temp = 2.0
+    expected_result = 1 / (1 + np.exp((1.0 * x - center) / max(1.0e-15, temp * KB_EV)))
+    assert np.allclose(fermi_dirac(x, center, temp), expected_result)
+
+
+@pytest.mark.filterwarnings("ignore:overflow encountered in exp.*:RuntimeWarning")
+def test_fermi_dirac_linbkg():
+    x = np.linspace(0, 10, 100)
+    center = 5.0
+    temp = 2.0
+    back0 = 1.0
+    back1 = 0.5
+    dos0 = 0.8
+    dos1 = 0.3
+    expected_result = (back0 + back1 * x) + (dos0 - back0 + (dos1 - back1) * x) / (
+        1 + np.exp((1.0 * x - center) / max(1.0e-15, temp * KB_EV))
+    )
+    assert np.allclose(
+        fermi_dirac_linbkg(x, center, temp, back0, back1, dos0, dos1), expected_result
+    )
+
+
+def test_fermi_dirac_linbkg_broad():
+    x = np.linspace(0, 10, 100)
+    center = 5.0
+    temp = 300.0
+    resolution = 0.1
+    back0 = 0.1
+    back1 = 0.2
+    dos0 = 0.3
+    dos1 = 0.4
+    expected_result = do_convolve(
+        x,
+        fermi_dirac_linbkg,
+        resolution=resolution,
+        center=center,
+        temp=temp,
+        back0=back0,
+        back1=back1,
+        dos0=dos0,
+        dos1=dos1,
+    )
+    assert np.allclose(
+        fermi_dirac_linbkg_broad(x, center, temp, resolution, back0, back1, dos0, dos1),
+        expected_result,
+    )
+
+
+def test_step_broad():
+    x = np.linspace(-10, 10, 100)
+    center = 0.0
+    sigma = 1.0
+    amplitude = 1.0
+    expected_result = (
+        amplitude
+        * 0.5
+        * scipy.special.erfc((1.0 * x - center) / max(1.0e-15, np.sqrt(2) * sigma))
+    )
+    assert np.allclose(step_broad(x, center, sigma, amplitude), expected_result)
+
+
+def test_step_linbkg_broad():
+    x = np.linspace(0, 10, 100)
+    center = 5.0
+    sigma = 2.0
+    back0 = 1.0
+    back1 = 0.5
+    dos0 = 2.0
+    dos1 = 1.0
+    expected_result = (back0 + back1 * x) + (dos0 - back0 + (dos1 - back1) * x) * (
+        step_broad(x, center, sigma, 1.0)
+    )
+    assert np.allclose(
+        step_linbkg_broad(x, center, sigma, back0, back1, dos0, dos1), expected_result
+    )
diff --git a/tests/test_image.py b/tests/test_image.py
new file mode 100644
index 00000000..16019475
--- /dev/null
+++ b/tests/test_image.py
@@ -0,0 +1,176 @@
+import numpy as np
+import xarray as xr
+import erlab.analysis as era
+
+
+def test_gaussian_filter():
+    # Create a test input DataArray
+    darr = xr.DataArray(np.arange(50, step=2).reshape((5, 5)), dims=["x", "y"])
+
+    # Define the expected output
+    expected_output = xr.DataArray(
+        np.array(
+            [
+                [3, 5, 7, 8, 10],
+                [10, 12, 14, 15, 17],
+                [20, 22, 24, 25, 27],
+                [29, 31, 33, 34, 36],
+                [36, 38, 40, 41, 43],
+            ]
+        ),
+        dims=["x", "y"],
+    )
+
+    # Apply the gaussian_filter function
+    result = era.image.gaussian_filter(darr, sigma={"x": 1.0, "y": 1.0})
+
+    # Check if the result matches the expected output
+    assert np.allclose(result, expected_output)
+
+
+def test_gaussian_laplace():
+    # Create a test input DataArray
+    darr = xr.DataArray(np.arange(50, step=2).reshape((5, 5)), dims=["x", "y"])
+
+    # Define the expected output
+    expected_output = xr.DataArray(
+        np.array(
+            [
+                [4, 4, 4, 4, 4],
+                [2, 2, 2, 2, 2],
+                [0, 0, 0, 0, 0],
+                [-2, -2, -2, -2, -2],
+                [-4, -4, -4, -4, -4],
+            ]
+        ),
+        dims=["x", "y"],
+    )
+
+    # Apply the gaussian_laplace function
+    result = era.image.gaussian_laplace(darr, sigma={"x": 1.0, "y": 1.0})
+
+    # Check if the result matches the expected output
+    assert np.allclose(result, expected_output)
+
+    # Additional test case
+    # Define the expected output
+    expected_output2 = xr.DataArray(
+        np.array(
+            [
+                [1, 1, 1, 1, 1],
+                [1, 1, 1, 1, 1],
+                [0, 0, 0, 0, 0],
+                [-1, -1, -1, -1, -1],
+                [-1, -1, -1, -1, -1],
+            ]
+        ),
+        dims=["x", "y"],
+    )
+
+    # Apply the gaussian_laplace function
+    result2 = era.image.gaussian_laplace(darr, sigma=2.0)
+
+    # Check if the result matches the expected output
+    assert np.allclose(result2, expected_output2)
+
+
+def test_laplace():
+    # Create a test input DataArray
+    darr = xr.DataArray(np.arange(50, step=2).reshape((5, 5)), dims=["x", "y"])
+
+    # Define the expected output
+    expected_output = xr.DataArray(
+        np.array(
+            [
+                [12, 10, 10, 10, 8],
+                [2, 0, 0, 0, -2],
+                [2, 0, 0, 0, -2],
+                [2, 0, 0, 0, -2],
+                [-8, -10, -10, -10, -12],
+            ]
+        ),
+        dims=["x", "y"],
+    )
+
+    # Apply the laplace function
+    result = era.image.laplace(darr)
+
+    # Check if the result matches the expected output
+    assert np.allclose(result, expected_output)
+
+
+def test_minimum_gradient():
+    # Create a test input DataArray
+    darr = xr.DataArray(np.arange(50, step=2).reshape((5, 5)), dims=["x", "y"])
+
+    # Define the expected output
+    expected_output = xr.DataArray(
+        np.array(
+            [
+                [0.0, 0.13608276, 0.27216553, 0.40824829, 0.58345997],
+                [0.49507377, 0.58834841, 0.68640647, 0.78446454, 0.89113279],
+                [0.99014754, 1.07863874, 1.17669681, 1.27475488, 1.38620656],
+                [1.48522131, 1.56892908, 1.66698715, 1.76504522, 1.88128033],
+                [2.91729983, 2.85773803, 2.9938208, 3.12990356, 3.17887766],
+            ]
+        ),
+        dims=["x", "y"],
+    )
+
+    # Apply the minimum_gradient function
+    result = era.image.minimum_gradient(darr).astype(np.float32)
+
+    # Check if the result matches the expected output
+    assert np.allclose(result, expected_output)
+
+
+def test_scaled_laplace():
+    # Create a test input DataArray
+    darr = xr.DataArray(
+        np.arange(50, step=2).reshape((5, 5)).astype(float), dims=["x", "y"]
+    )
+
+    # Define the expected output
+    expected_output = xr.DataArray(
+        np.array(
+            [
+                [12.0, 10.0, 10.0, 10.0, 8.0],
+                [2.0, 0.0, 0.0, 0.0, -2.0],
+                [2.0, 0.0, 0.0, 0.0, -2.0],
+                [2.0, 0.0, 0.0, 0.0, -2.0],
+                [-8.0, -10.0, -10.0, -10.0, -12.0],
+            ]
+        ),
+        dims=["x", "y"],
+    )
+
+    # Apply the scaled_laplace function
+    result = era.image.scaled_laplace(darr)
+
+    # Check if the result matches the expected output
+    assert np.allclose(result, expected_output)
+
+
+def test_curvature():
+    # Create a test input DataArray
+    darr = xr.DataArray(np.arange(50, step=2).reshape((5, 5)), dims=["x", "y"]) ** 2
+
+    # Define the expected output
+    expected_output = xr.DataArray(
+        np.array(
+            [
+                [0.11852942, 0.11855069, 0.11778077, 0.11184719, 0.10558571],
+                [0.16448492, 0.16107288, 0.15683689, 0.14772279, 0.1385244],
+                [0.17403091, 0.16649966, 0.15876051, 0.14719689, 0.1361662],
+                [0.09486956, 0.09038468, 0.08598027, 0.0783563, 0.07130584],
+                [0.05139264, 0.04942051, 0.04746361, 0.04241876, 0.03781145],
+            ]
+        ),
+        dims=["x", "y"],
+    )
+
+    # Apply the curvature function
+    result = era.image.curvature(darr).astype(np.float32)
+
+    # Check if the result matches the expected output
+    assert np.allclose(result, expected_output)
diff --git a/tests/test_utilities.py b/tests/test_utilities.py
new file mode 100644
index 00000000..006a26be
--- /dev/null
+++ b/tests/test_utilities.py
@@ -0,0 +1,25 @@
+import numpy as np
+import xarray as xr
+from erlab.analysis.utilities import shift
+
+
+def test_shift():
+    # Create a test input DataArray
+    darr = xr.DataArray(
+        np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]]).astype(float), dims=["x", "y"]
+    )
+
+    # Create a test shift DataArray
+    shift_arr = xr.DataArray([1, 0, 2], dims=["x"])
+
+    # Perform the shift operation
+    shifted = shift(darr, shift_arr, along="y")
+
+    # Define the expected result
+    expected = xr.DataArray(
+        np.array([[np.nan, 1.0, 2.0], [4.0, 5.0, 6.0], [np.nan, np.nan, 7.0]]),
+        dims=["x", "y"],
+    )
+
+    # Check if the shifted array matches the expected result
+    assert np.allclose(shifted, expected, equal_nan=True)