Image sharpness and finding the focus point#

This notebook illustrates how to find the focus point of an imaging detector by computing the image sharpness as a function of camera position.

[1]:

import scipp as sc
from ess import imaging as img
from ess.imaging import data
import scippnexus as sx
import numpy as np

Load the data#

We load the image data, and attach the camera stage position to the data array.

[2]:

dg = sx.load(data.get_tbl_orca_focussing_path())
da = dg["entry"]["instrument"]["orca_detector"]["data"]
da.coords["position"] = dg["entry"]["instrument"]["camera_stage"]["position_setpoint"][
    "value"
].data
da = da.rename_dims(time="position")
da

Downloading file 'tbl-orca-focussing.hdf.zip' from 'https://public.esss.dk/groups/scipp/ess/imaging/1/tbl-orca-focussing.hdf.zip' to '/home/runner/.cache/essimaging/1'.
Unzipping contents of '/home/runner/.cache/essimaging/1/tbl-orca-focussing.hdf.zip' to '/home/runner/.cache/essimaging/1/tbl-orca-focussing.hdf.zip.unzip'

[2]:
Show/Hide data repr Show/Hide attributes
scipp.DataArray (61.00 MB)
    • position: 61
    • dim_1: 512
    • dim_2: 512
    • position
      (position)
      float64
      mm
      140.0, 141.0, ..., 199.0, 200.0
      Values:
      array([140., 141., 142., 143., 144., 145., 146., 147., 148., 149., 150.,
       151., 152., 153., 154., 155., 156., 157., 158., 159., 160., 161.,
       162., 163., 164., 165., 166., 167., 168., 169., 170., 171., 172.,
       173., 174., 175., 176., 177., 178., 179., 180., 181., 182., 183.,
       184., 185., 186., 187., 188., 189., 190., 191., 192., 193., 194.,
       195., 196., 197., 198., 199., 200.])
    • time
      (position)
      datetime64
      ns
      2025-06-11T10:07:14.603685891, 2025-06-11T10:07:19.443681068, ..., 2025-06-11T10:12:15.939916263, 2025-06-11T10:12:20.779694055
      Values:
      array(['2025-06-11T10:07:14.603685891', '2025-06-11T10:07:19.443681068',
       '2025-06-11T10:07:24.291694836', '2025-06-11T10:07:29.135681029',
       '2025-06-11T10:07:36.387704881', '2025-06-11T10:07:41.239693346',
       '2025-06-11T10:07:46.111723518', '2025-06-11T10:07:50.687718363',
       '2025-06-11T10:07:55.555682998', '2025-06-11T10:08:00.435680492',
       '2025-06-11T10:08:05.283693959', '2025-06-11T10:08:10.023683844',
       '2025-06-11T10:08:14.879729086', '2025-06-11T10:08:19.735700903',
       '2025-06-11T10:08:24.587683636', '2025-06-11T10:08:29.439675066',
       '2025-06-11T10:08:33.991718005', '2025-06-11T10:08:38.851754030',
       '2025-06-11T10:08:43.699704502', '2025-06-11T10:08:48.555669860',
       '2025-06-11T10:08:53.403690390', '2025-06-11T10:08:58.247718146',
       '2025-06-11T10:09:03.099704735', '2025-06-11T10:09:07.951692452',
       '2025-06-11T10:09:12.787697950', '2025-06-11T10:09:20.067721706',
       '2025-06-11T10:09:24.919730518', '2025-06-11T10:09:29.763706276',
       '2025-06-11T10:09:34.619683746', '2025-06-11T10:09:39.467698542',
       '2025-06-11T10:09:46.723715712', '2025-06-11T10:09:51.387727129',
       '2025-06-11T10:09:56.155690987', '2025-06-11T10:10:00.983695240',
       '2025-06-11T10:10:05.835713180', '2025-06-11T10:10:10.679685428',
       '2025-06-11T10:10:17.947739370', '2025-06-11T10:10:22.803701747',
       '2025-06-11T10:10:27.671682040', '2025-06-11T10:10:34.923718875',
       '2025-06-11T10:10:39.475679906', '2025-06-11T10:10:44.327714300',
       '2025-06-11T10:10:49.179701241', '2025-06-11T10:10:54.047723775',
       '2025-06-11T10:10:58.907683420', '2025-06-11T10:11:03.775723131',
       '2025-06-11T10:11:08.639689393', '2025-06-11T10:11:13.515701420',
       '2025-06-11T10:11:18.375689596', '2025-06-11T10:11:23.247713489',
       '2025-06-11T10:11:28.123680303', '2025-06-11T10:11:32.959681624',
       '2025-06-11T10:11:37.719680831', '2025-06-11T10:11:42.255678954',
       '2025-06-11T10:11:49.539686860', '2025-06-11T10:11:54.379713851',
       '2025-06-11T10:11:59.231692792', '2025-06-11T10:12:06.499695867',
       '2025-06-11T10:12:11.347714105', '2025-06-11T10:12:15.939916263',
       '2025-06-11T10:12:20.779694055'], dtype='datetime64[ns]')
    • (position, dim_1, dim_2)
      int32
      𝟙
      30213, 31248, ..., 7737, 7818
      Values:
      array([[[30213, 31248, 31972, ..., 21162, 20654, 19905],
        [30118, 31154, 31709, ..., 21594, 20768, 20430],
        [30200, 31051, 31501, ..., 20960, 20528, 20365],
        ...,
        [ 9970, 10262, 10097, ...,  7853,  7926,  8061],
        [ 9821,  9835, 10324, ...,  7721,  7617,  7816],
        [ 9373,  9988,  9827, ...,  7676,  7560,  7771]],

       [[30301, 31097, 31860, ..., 21286, 20226, 20102],
        [30427, 30812, 31838, ..., 21106, 20534, 20082],
        [29882, 31460, 31381, ..., 21053, 20535, 19926],
        ...,
        [ 9895, 10018,  9911, ...,  7965,  7917,  8088],
        [ 9624,  9767, 10171, ...,  7541,  7406,  7804],
        [ 9838,  9857, 10023, ...,  7556,  7388,  7705]],

       [[30044, 30657, 32174, ..., 20952, 20613, 20032],
        [29852, 31119, 31492, ..., 20668, 20736, 19924],
        [29956, 31202, 31819, ..., 21287, 20653, 20176],
        ...,
        [ 9786, 10056,  9942, ...,  7826,  7724,  8095],
        [ 9652,  9738, 10020, ...,  7573,  7685,  7855],
        [ 9637,  9689,  9995, ...,  7800,  7656,  7662]],

       ...,

       [[16267, 17469, 19609, ..., 15784, 15275, 14631],
        [16093, 17714, 19456, ..., 16212, 15494, 14393],
        [15566, 17114, 19024, ..., 16089, 15430, 14834],
        ...,
        [ 7755,  8114,  8318, ...,  7903,  8146,  7829],
        [ 7980,  7764,  8159, ...,  7680,  7954,  8030],
        [ 7802,  7929,  7885, ...,  7509,  7756,  7825]],

       [[15887, 17342, 19103, ..., 15943, 14901, 14436],
        [15416, 17493, 19189, ..., 16118, 15433, 14608],
        [15578, 17143, 19417, ..., 16151, 15445, 14567],
        ...,
        [ 7976,  8048,  8008, ...,  8038,  7926,  7682],
        [ 8040,  7800,  7903, ...,  7727,  8146,  7988],
        [ 7648,  7621,  8036, ...,  7444,  7760,  7836]],

       [[15713, 16996, 19101, ..., 15787, 15245, 14246],
        [15462, 17378, 18813, ..., 15957, 15037, 14255],
        [15247, 17023, 18774, ..., 16178, 15039, 14420],
        ...,
        [ 7876,  8267,  7935, ...,  8040,  8065,  7620],
        [ 7759,  7896,  7952, ...,  7567,  7853,  8073],
        [ 7553,  7772,  7796, ...,  7646,  7737,  7818]]],
      shape=(61, 512, 512), dtype=int32)

The data contains 61 images, one for each camera position.

We can plot the first image in the sequence to get an impression of what it looks like:

[3]:

da["position", 0].plot(aspect="equal")

[3]:
../_images/tools_sharpness-and-focus-point_5_0.svg

Computing image sharpness#

The ess.imaging.tools module provides tools for image manipulation. The sharpness function computes image sharpness over two spatial dimensions (that need to be specified). The remaining position dimension will remain in the output.

[4]:

sharp = img.tools.sharpness(da, dims=["dim_1", "dim_2"])
sharp

[4]:
Show/Hide data repr Show/Hide attributes
scipp.DataArray (2.55 KB)
    • position: 61
    • position
      (position)
      float64
      mm
      140.0, 141.0, ..., 199.0, 200.0
      Values:
      array([140., 141., 142., 143., 144., 145., 146., 147., 148., 149., 150.,
       151., 152., 153., 154., 155., 156., 157., 158., 159., 160., 161.,
       162., 163., 164., 165., 166., 167., 168., 169., 170., 171., 172.,
       173., 174., 175., 176., 177., 178., 179., 180., 181., 182., 183.,
       184., 185., 186., 187., 188., 189., 190., 191., 192., 193., 194.,
       195., 196., 197., 198., 199., 200.])
    • time
      (position)
      datetime64
      ns
      2025-06-11T10:07:14.603685891, 2025-06-11T10:07:19.443681068, ..., 2025-06-11T10:12:15.939916263, 2025-06-11T10:12:20.779694055
      Values:
      array(['2025-06-11T10:07:14.603685891', '2025-06-11T10:07:19.443681068',
       '2025-06-11T10:07:24.291694836', '2025-06-11T10:07:29.135681029',
       '2025-06-11T10:07:36.387704881', '2025-06-11T10:07:41.239693346',
       '2025-06-11T10:07:46.111723518', '2025-06-11T10:07:50.687718363',
       '2025-06-11T10:07:55.555682998', '2025-06-11T10:08:00.435680492',
       '2025-06-11T10:08:05.283693959', '2025-06-11T10:08:10.023683844',
       '2025-06-11T10:08:14.879729086', '2025-06-11T10:08:19.735700903',
       '2025-06-11T10:08:24.587683636', '2025-06-11T10:08:29.439675066',
       '2025-06-11T10:08:33.991718005', '2025-06-11T10:08:38.851754030',
       '2025-06-11T10:08:43.699704502', '2025-06-11T10:08:48.555669860',
       '2025-06-11T10:08:53.403690390', '2025-06-11T10:08:58.247718146',
       '2025-06-11T10:09:03.099704735', '2025-06-11T10:09:07.951692452',
       '2025-06-11T10:09:12.787697950', '2025-06-11T10:09:20.067721706',
       '2025-06-11T10:09:24.919730518', '2025-06-11T10:09:29.763706276',
       '2025-06-11T10:09:34.619683746', '2025-06-11T10:09:39.467698542',
       '2025-06-11T10:09:46.723715712', '2025-06-11T10:09:51.387727129',
       '2025-06-11T10:09:56.155690987', '2025-06-11T10:10:00.983695240',
       '2025-06-11T10:10:05.835713180', '2025-06-11T10:10:10.679685428',
       '2025-06-11T10:10:17.947739370', '2025-06-11T10:10:22.803701747',
       '2025-06-11T10:10:27.671682040', '2025-06-11T10:10:34.923718875',
       '2025-06-11T10:10:39.475679906', '2025-06-11T10:10:44.327714300',
       '2025-06-11T10:10:49.179701241', '2025-06-11T10:10:54.047723775',
       '2025-06-11T10:10:58.907683420', '2025-06-11T10:11:03.775723131',
       '2025-06-11T10:11:08.639689393', '2025-06-11T10:11:13.515701420',
       '2025-06-11T10:11:18.375689596', '2025-06-11T10:11:23.247713489',
       '2025-06-11T10:11:28.123680303', '2025-06-11T10:11:32.959681624',
       '2025-06-11T10:11:37.719680831', '2025-06-11T10:11:42.255678954',
       '2025-06-11T10:11:49.539686860', '2025-06-11T10:11:54.379713851',
       '2025-06-11T10:11:59.231692792', '2025-06-11T10:12:06.499695867',
       '2025-06-11T10:12:11.347714105', '2025-06-11T10:12:15.939916263',
       '2025-06-11T10:12:20.779694055'], dtype='datetime64[ns]')
    • (position)
      float64
      𝟙
      1.664e+08, 1.737e+08, ..., 3.636e+08, 3.592e+08
      Values:
      array([1.66367877e+08, 1.73721491e+08, 1.80515730e+08, 1.87450301e+08,
       1.95839357e+08, 2.03322584e+08, 2.11381659e+08, 2.21143653e+08,
       2.29351388e+08, 2.37960144e+08, 2.46379092e+08, 2.54095879e+08,
       2.62174234e+08, 2.69654120e+08, 2.77490212e+08, 2.87548060e+08,
       2.97669093e+08, 3.04265665e+08, 3.12310314e+08, 3.24953189e+08,
       3.32765238e+08, 3.38134294e+08, 3.46568909e+08, 3.53999035e+08,
       3.61792635e+08, 3.70058742e+08, 3.74338346e+08, 3.79159907e+08,
       3.85679708e+08, 3.89697049e+08, 3.91578219e+08, 3.98688639e+08,
       4.02523118e+08, 4.03725026e+08, 4.08685857e+08, 4.11394311e+08,
       4.10592086e+08, 4.13450995e+08, 4.17550265e+08, 4.20062414e+08,
       4.18278889e+08, 4.18660099e+08, 4.19266356e+08, 4.18551455e+08,
       4.16822756e+08, 4.14089584e+08, 4.15115642e+08, 4.14020582e+08,
       4.10765746e+08, 4.08197241e+08, 4.04466194e+08, 4.00390526e+08,
       3.97507864e+08, 3.96295749e+08, 3.91828606e+08, 3.83856874e+08,
       3.78202144e+08, 3.74145507e+08, 3.69857673e+08, 3.63584695e+08,
       3.59170710e+08])

A simple plot of the sharpness reveals a ‘best position’ where the sharpness peaks (around 180 mm):

[5]:

fig = sharp.plot()
fig

[5]:
../_images/tools_sharpness-and-focus-point_9_0.svg

Find the sharpest image#

To accurately find the sharpest image, and the associated camera stage position with the best focus, we filter outliers by applying a 75th percentile threshold to the data. We then compute the mean position weighted by the sharpness on the remaining points.

[6]:

threshold = sc.scalar(np.percentile(sharp.values, 75))
sel = sharp.data >= threshold
subset = sharp[sel]
subset

[6]:
Show/Hide data repr Show/Hide attributes
scipp.DataArray (1.49 KB)
    • position: 16
    • position
      (position)
      float64
      mm
      174.0, 175.0, ..., 188.0, 189.0
      Values:
      array([174., 175., 176., 177., 178., 179., 180., 181., 182., 183., 184.,
       185., 186., 187., 188., 189.])
    • time
      (position)
      datetime64
      ns
      2025-06-11T10:10:05.835713180, 2025-06-11T10:10:10.679685428, ..., 2025-06-11T10:11:18.375689596, 2025-06-11T10:11:23.247713489
      Values:
      array(['2025-06-11T10:10:05.835713180', '2025-06-11T10:10:10.679685428',
       '2025-06-11T10:10:17.947739370', '2025-06-11T10:10:22.803701747',
       '2025-06-11T10:10:27.671682040', '2025-06-11T10:10:34.923718875',
       '2025-06-11T10:10:39.475679906', '2025-06-11T10:10:44.327714300',
       '2025-06-11T10:10:49.179701241', '2025-06-11T10:10:54.047723775',
       '2025-06-11T10:10:58.907683420', '2025-06-11T10:11:03.775723131',
       '2025-06-11T10:11:08.639689393', '2025-06-11T10:11:13.515701420',
       '2025-06-11T10:11:18.375689596', '2025-06-11T10:11:23.247713489'],
      dtype='datetime64[ns]')
    • (position)
      float64
      𝟙
      4.087e+08, 4.114e+08, ..., 4.108e+08, 4.082e+08
      Values:
      array([4.08685857e+08, 4.11394311e+08, 4.10592086e+08, 4.13450995e+08,
       4.17550265e+08, 4.20062414e+08, 4.18278889e+08, 4.18660099e+08,
       4.19266356e+08, 4.18551455e+08, 4.16822756e+08, 4.14089584e+08,
       4.15115642e+08, 4.14020582e+08, 4.10765746e+08, 4.08197241e+08])
[7]:

# Compute weighted mean of selected points
pos = (subset.coords["position"] * subset.data).sum() / subset.data.sum()
pos

[7]:
Show/Hide data repr Show/Hide attributes
scipp.Variable (264 Bytes)
    • ()
      float64
      mm
      181.49986403473855
      Values:
      array(181.49986403)
[8]:

fig.ax.axhline(threshold.value, color="lightgray", label="threshold")
fig.ax.axvline(pos.value, color="k", label="focus point")
fig.ax.legend()
fig

[8]:
../_images/tools_sharpness-and-focus-point_13_0.svg

Image comparison#

Finally, we can quickly compare the first image in the stack to the sharpest image determined above.

[9]:

# Show a smaller region to highlight details
dx = 100
xslice = ("dim_2", slice(363, 363 + dx))
yslice = ("dim_1", slice(90, 90 + dx))

# Find the closest position to the focus point
closest = subset.coords["position"][
    np.argmin(sc.abs(subset.coords["position"] - pos).values)
]

da["position", 0][xslice][yslice].plot(aspect="equal", title="First image") + da[
    "position", closest
][xslice][yslice].plot(aspect="equal", title="Sharpest image")

[9]:
../_images/tools_sharpness-and-focus-point_15_0.svg