Add utilities to compute image sharpness

docs/api-reference/index.md

@@ -29,5 +29,6 @@
 
    normalize
    io
+   tools
    workflow
 ```

docs/user-guide/tools/index.md

@@ -0,0 +1,9 @@
+# Tools
+
+```{toctree}
+---
+maxdepth: 1
+---
+
+sharpness-and-focus-point
+```

docs/user-guide/tools/sharpness-and-focus-point.ipynb

@@ -0,0 +1,218 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "0",
+   "metadata": {},
+   "source": [
+    "# Image sharpness and finding the focus point\n",
+    "\n",
+    "This notebook illustrates how to find the focus point of an imaging detector by computing the image sharpness as a function of camera position."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "1",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import scipp as sc\n",
+    "from ess import imaging as img\n",
+    "from ess.imaging import data\n",
+    "import scippnexus as sx\n",
+    "import numpy as np"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "2",
+   "metadata": {},
+   "source": [
+    "## Load the data\n",
+    "\n",
+    "We load the image data, and attach the camera stage position to the data array."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "3",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "dg = sx.load(data.get_tbl_orca_focussing_path())\n",
+    "da = dg[\"entry\"][\"instrument\"][\"orca_detector\"][\"data\"]\n",
+    "da.coords[\"position\"] = dg[\"entry\"][\"instrument\"][\"camera_stage\"][\"position_setpoint\"][\n",
+    "    \"value\"\n",
+    "].data\n",
+    "da = da.rename_dims(time=\"position\")\n",
+    "da"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "4",
+   "metadata": {},
+   "source": [
+    "The data contains 61 images, one for each camera position.\n",
+    "\n",
+    "We can plot the first image in the sequence to get an impression of what it looks like:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "5",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "da[\"position\", 0].plot(aspect=\"equal\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "6",
+   "metadata": {},
+   "source": [
+    "## Computing image sharpness\n",
+    "\n",
+    "The `ess.imaging.tools` module provides tools for image manipulation.\n",
+    "The `sharpness` function computes image sharpness over two spatial dimensions (that need to be specified).\n",
+    "The remaining `position` dimension will remain in the output."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "7",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "sharp = img.tools.sharpness(da, dims=[\"dim_1\", \"dim_2\"])\n",
+    "sharp"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "8",
+   "metadata": {},
+   "source": [
+    "A simple plot of the sharpness reveals a 'best position' where the sharpness peaks (around 180 mm):"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "9",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fig = sharp.plot()\n",
+    "fig"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "10",
+   "metadata": {},
+   "source": [
+    "## Find the sharpest image\n",
+    "\n",
+    "To accurately find the sharpest image, and the associated camera stage position with the best focus,\n",
+    "we filter outliers by applying a 75th percentile threshold to the data.\n",
+    "We then compute the mean position weighted by the sharpness on the remaining points."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "11",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "threshold = sc.scalar(np.percentile(sharp.values, 75))\n",
+    "sel = sharp.data >= threshold\n",
+    "subset = sharp[sel]\n",
+    "subset"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "12",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Compute weighted mean of selected points\n",
+    "pos = (subset.coords[\"position\"] * subset.data).sum() / subset.data.sum()\n",
+    "pos"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "13",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fig.ax.axhline(threshold.value, color=\"lightgray\", label=\"threshold\")\n",
+    "fig.ax.axvline(pos.value, color=\"k\", label=\"focus point\")\n",
+    "fig.ax.legend()\n",
+    "fig"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "14",
+   "metadata": {},
+   "source": [
+    "## Image comparison\n",
+    "\n",
+    "Finally, we can quickly compare the first image in the stack to the sharpest image determined above."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "15",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Show a smaller region to highlight details\n",
+    "dx = 100\n",
+    "xslice = (\"dim_2\", slice(363, 363 + dx))\n",
+    "yslice = (\"dim_1\", slice(90, 90 + dx))\n",
+    "\n",
+    "# Find the closest position to the focus point\n",
+    "closest = subset.coords[\"position\"][\n",
+    "    np.argmin(sc.abs(subset.coords[\"position\"] - pos).values)\n",
+    "]\n",
+    "\n",
+    "da[\"position\", 0][xslice][yslice].plot(aspect=\"equal\", title=\"First image\") + da[\n",
+    "    \"position\", closest\n",
+    "][xslice][yslice].plot(aspect=\"equal\", title=\"Sharpest image\")"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3 (ipykernel)",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}

pyproject.toml

@@ -48,6 +48,7 @@ test = [
     "pytest>=7.0",
     "pooch>=1.5",
     "tof>=24.7.0",
+    "scitiff>=24.6.0",
 ]
 
 [project.urls]

requirements/basetest.in

@@ -10,3 +10,4 @@
 pytest>=7.0
 pooch>=1.5
 tof>=24.7.0
+scitiff>=24.6.0

requirements/basetest.txt

@@ -1,10 +1,16 @@
-# SHA1:a2facda5876f52f2eaebf7967f39038ae196fb7c
+# SHA1:b13c3d6374fbf2f0a67f1c13edace05a1dc456ac
 #
 # This file was generated by pip-compile-multi.
 # To update, run:
 #
 #    requirements upgrade
 #
+annotated-types==0.7.0
+    # via pydantic
+attrs==25.3.0
+    # via
+    #   jsonschema
+    #   referencing
 certifi==2025.8.3
     # via requests
 charset-normalizer==3.4.2
@@ -19,6 +25,10 @@ idna==3.10
     # via requests
 iniconfig==2.1.0
     # via pytest
+jsonschema==4.25.0
+    # via scitiff
+jsonschema-specifications==2025.4.1
+    # via jsonschema
 kiwisolver==1.4.8
     # via matplotlib
 lazy-loader==0.4
@@ -33,6 +43,7 @@ numpy==2.3.2
     #   matplotlib
     #   scipp
     #   scipy
+    #   tifffile
 packaging==25.0
     # via
     #   lazy-loader
@@ -49,6 +60,10 @@ pluggy==1.6.0
     # via pytest
 pooch==1.8.2
     # via -r basetest.in
+pydantic==2.11.7
+    # via scitiff
+pydantic-core==2.33.2
+    # via pydantic
 pygments==2.19.2
     # via pytest
 pyparsing==3.2.3
@@ -57,15 +72,37 @@ pytest==8.4.1
     # via -r basetest.in
 python-dateutil==2.9.0.post0
     # via matplotlib
+referencing==0.36.2
+    # via
+    #   jsonschema
+    #   jsonschema-specifications
 requests==2.32.4
     # via pooch
+rpds-py==0.26.0
+    # via
+    #   jsonschema
+    #   referencing
 scipp==25.5.1
-    # via tof
+    # via
+    #   scitiff
+    #   tof
 scipy==1.16.1
     # via tof
+scitiff==25.7.0
+    # via -r basetest.in
 six==1.17.0
     # via python-dateutil
+tifffile==2025.6.11
+    # via scitiff
 tof==25.5.0
     # via -r basetest.in
+typing-extensions==4.14.1
+    # via
+    #   pydantic
+    #   pydantic-core
+    #   referencing
+    #   typing-inspection
+typing-inspection==0.4.1
+    # via pydantic
 urllib3==2.5.0
     # via requests

requirements/dev.txt

@@ -53,6 +53,7 @@ jsonschema[format-nongpl]==4.25.0
     #   jupyter-events
     #   jupyterlab-server
     #   nbformat
+    #   scitiff
 jupyter-events==0.12.0
     # via jupyter-server
 jupyter-lsp==2.2.6

requirements/nightly.in

@@ -6,6 +6,7 @@ graphviz
 tifffile>=2024.7.2
 pytest>=7.0
 pooch>=1.5
+scitiff>=24.6.0
 scipp
 --index-url=https://pypi.anaconda.org/scipp-nightly-wheels/simple/
 --extra-index-url=https://pypi.org/simple