From d06e8dfe0fb1b9dc97e7e1efd1540ad0f1ebd71c Mon Sep 17 00:00:00 2001 From: Johannes Kasimir Date: Mon, 13 Apr 2026 09:36:28 +0200 Subject: [PATCH 1/3] fix: handle binned data in max res finder --- .../src/ess/imaging/tools/resolution.py | 17 +++++--- .../tools/maximum_resolution_monitor_test.py | 41 ++++++++++++++++++- 2 files changed, 51 insertions(+), 7 deletions(-) diff --git a/packages/essimaging/src/ess/imaging/tools/resolution.py b/packages/essimaging/src/ess/imaging/tools/resolution.py index c344df499..cd94f8f6e 100644 --- a/packages/essimaging/src/ess/imaging/tools/resolution.py +++ b/packages/essimaging/src/ess/imaging/tools/resolution.py @@ -7,7 +7,6 @@ def maximum_resolution_achievable( events: sc.DataArray, coarse_x_bin_edges: sc.Variable, coarse_y_bin_edges: sc.Variable, - time_bin_edges: sc.Variable, max_tries: int = 10, max_pixels_x: int = 2048, max_pixels_y: int = 2048, @@ -31,8 +30,6 @@ def maximum_resolution_achievable( Minimum acceptable resolution in ``x``. coarse_y_bin_edges: Minimum acceptable resolution in ``y``. - time_bin_edges: - Desired resolution in ``t``. max_tries: The maximum number of iterations before giving up. max_pixels_x: @@ -59,7 +56,17 @@ def maximum_resolution_achievable( nx = int(2**0.5 * lower_nx) + 1 ny = int(2**0.5 * lower_ny) + 1 - events = events.bin({time_bin_edges.dim: time_bin_edges}) + + events = events.copy(deep=False) + + if events.bins is not None: + image_dims = (coarse_x_bin_edges.dim, coarse_y_bin_edges.dim) + for c in image_dims: + if c not in events.bins.coords: + events.bins.coords[c] = sc.bins_like( + events, sc.midpoints(events.coords[c]) + ) + events = events.bins.concat(image_dims) for _ in range(max_tries): xbins = sc.linspace( @@ -88,7 +95,7 @@ def maximum_resolution_achievable( upper_nx = nx upper_ny = ny nx = min(round((lower_nx * nx) ** 0.5), upper_nx - 1) - ny = min(round((lower_ny * ny) ** 0.5), upper_nx - 1) + ny = min(round((lower_ny * ny) ** 0.5), upper_ny - 1) if upper_nx - lower_nx < 2 and upper_ny - lower_ny < 2: break diff --git a/packages/essimaging/tests/imaging/tools/maximum_resolution_monitor_test.py b/packages/essimaging/tests/imaging/tools/maximum_resolution_monitor_test.py index 5c2dabb37..a1b1a9a59 100644 --- a/packages/essimaging/tests/imaging/tools/maximum_resolution_monitor_test.py +++ b/packages/essimaging/tests/imaging/tools/maximum_resolution_monitor_test.py @@ -25,7 +25,6 @@ def test_finds_maximum_resolution(): events, sc.linspace('x', 0, 1, 2), sc.linspace('x', 0, 1, 2), - sc.linspace('t', 0, 2, 2), ) assert len(x_be) == 11 assert len(y_be) == 11 @@ -51,7 +50,6 @@ def test_finds_maximum_resolution_random(seed): events, sc.linspace('x', 0, 1, 2), sc.linspace('y', 0, 1, 2), - sc.linspace('t', 0, 2, 2), # Need enough tries to be sure we find the optimum max_tries=100, ) @@ -70,3 +68,42 @@ def test_finds_maximum_resolution_random(seed): .value == 0 ) + + +def test_finds_maximum_resolution_binned_input(): + np.random.seed(0) + n = np.random.randint(1000, 100_000) + events = sc.DataArray( + sc.ones(dims=['events'], shape=(n,)), + coords={ + 'x': sc.array(dims=['events'], values=np.random.random(n)), + 'y': sc.array(dims=['events'], values=np.random.random(n)), + 't': sc.array(dims=['events'], values=np.random.random(n)), + }, + ) + events = events.bin(x=100, y=100, t=500) + del events.bins.coords['x'] + del events.bins.coords['y'] + x_be, y_be = maximum_resolution_achievable( + events, + sc.linspace('x', 0, 1, 2), + sc.linspace('y', 0, 1, 2), + # Need enough tries to be sure we find the optimum + max_tries=100, + ) + + events.bins.coords['x'] = sc.bins_like(events, sc.midpoints(events.coords['x'])) + events.bins.coords['y'] = sc.bins_like(events, sc.midpoints(events.coords['y'])) + events = events.bins.concat(['x', 'y']) + + assert events.bin(x=x_be, y=y_be).bins.size().min().value > 0 + assert ( + events.bin( + x=sc.linspace('x', 0, 1, len(x_be) + 1), + y=sc.linspace('y', 0, 1, len(y_be) + 1), + ) + .bins.size() + .min() + .value + == 0 + ) From 622a383f0acca774f8cb8faef8f84df4f72df496 Mon Sep 17 00:00:00 2001 From: Johannes Kasimir Date: Mon, 13 Apr 2026 14:58:19 +0200 Subject: [PATCH 2/3] refactor: simplify and only ever reduce number of bins by factor of 2 in each dimension --- .../src/ess/imaging/tools/resolution.py | 137 ++++++------------ .../tools/maximum_resolution_monitor_test.py | 115 +++++++++------ 2 files changed, 109 insertions(+), 143 deletions(-) diff --git a/packages/essimaging/src/ess/imaging/tools/resolution.py b/packages/essimaging/src/ess/imaging/tools/resolution.py index cd94f8f6e..5e711edb5 100644 --- a/packages/essimaging/src/ess/imaging/tools/resolution.py +++ b/packages/essimaging/src/ess/imaging/tools/resolution.py @@ -1,3 +1,5 @@ +from math import floor, log2 + import numpy as np import scipp as sc from numpy.typing import NDArray @@ -5,121 +7,66 @@ def maximum_resolution_achievable( events: sc.DataArray, - coarse_x_bin_edges: sc.Variable, - coarse_y_bin_edges: sc.Variable, - max_tries: int = 10, - max_pixels_x: int = 2048, - max_pixels_y: int = 2048, - raise_if_not_maximum: bool = False, + image_dims: tuple[str, str], ): """ Estimates the maximum resolution achievable - given a desired binning in time. + given a desired binning. The maximum achievable resolution is defined as the resolution in ``xy`` such that - there is at least one event in every ``xyt`` pixel. + there is at least one event in every ``xy...`` pixel, + where ``...`` represents the rest of the binning dimenensions + of the event data, typically wavelength or time-of-flight. Parameters ------------- events: - 1D DataArray containing events with associated x, y, and t coordinates. - The names of the coordinates must not be `x`, `y` and `t`, - the names of the coordinates are taken from the provided ``bin_edges`` - for each respective dimension. - coarse_x_bin_edges: - Minimum acceptable resolution in ``x``. - coarse_y_bin_edges: - Minimum acceptable resolution in ``y``. - max_tries: - The maximum number of iterations before giving up. - max_pixels_x: - The maximum number of pixels in ``x``. - max_pixels_y: - The maximum number of pixels in ``y``. - raise_if_not_maximum: - Often it is not important to find the exact maximum resolution. - Therefore this parameter is ``False`` by default, and the function - returns an estimate of the maximum resolution. - If you want the returned resolution to be exactly the maximum resolution, - set the value of this parameter to ``True``. + DataArray binned in at least the ``image_dims`` dimensions and optionally more. + The number of bins in each ``image_dims`` dimension ought to be power of ``2``. Returns ------------- - The bin edges in x respectively y that define the - maximum achievable resolution. + The event data array folded to the coarsest possible resolution. """ + if events.bins is None: + raise ValueError( + 'Input data must be binned to the finest possible pixel binning.' + ' For best result, number of bins in each image dimension should' + ' be a power of 2.' + ) - lower_nx = coarse_x_bin_edges.size - lower_ny = coarse_y_bin_edges.size - upper_nx = max_pixels_x - upper_ny = max_pixels_y - - nx = int(2**0.5 * lower_nx) + 1 - ny = int(2**0.5 * lower_ny) + 1 - - events = events.copy(deep=False) - - if events.bins is not None: - image_dims = (coarse_x_bin_edges.dim, coarse_y_bin_edges.dim) - for c in image_dims: - if c not in events.bins.coords: - events.bins.coords[c] = sc.bins_like( - events, sc.midpoints(events.coords[c]) - ) - events = events.bins.concat(image_dims) + x, y = image_dims - for _ in range(max_tries): - xbins = sc.linspace( - coarse_x_bin_edges.dim, coarse_x_bin_edges[0], coarse_x_bin_edges[-1], nx - ) - ybins = sc.linspace( - coarse_y_bin_edges.dim, coarse_y_bin_edges[0], coarse_y_bin_edges[-1], ny - ) - min_counts_per_pixel = ( - events.bin( - { - coarse_x_bin_edges.dim: xbins, - coarse_y_bin_edges.dim: ybins, - } + for d in image_dims: + if events.sizes[d] % 2 == 1: + raise ValueError( + 'Input data must have an even number of bins in each image dimension.' + ' For best result, number of bins in each image dimension should' + ' be a power of 2.' ) - .bins.size() - .min() - ) - if min_counts_per_pixel.value > 0: - lower_nx = nx - lower_ny = ny - nx = max(min(round((upper_nx * nx) ** 0.5), nx * 2), lower_nx + 1) - ny = max(min(round((upper_ny * ny) ** 0.5), ny * 2), lower_ny + 1) - else: - upper_nx = nx - upper_ny = ny - nx = min(round((lower_nx * nx) ** 0.5), upper_nx - 1) - ny = min(round((lower_ny * ny) ** 0.5), upper_ny - 1) - - if upper_nx - lower_nx < 2 and upper_ny - lower_ny < 2: - break + refinements = ( + min( + floor(log2(events.sizes[image_dims[0]])), + floor(log2(events.sizes[image_dims[1]])), + ) + + 1 + ) - if raise_if_not_maximum and upper_nx - lower_nx >= 2 and upper_ny - lower_ny >= 2: - raise RuntimeError( - 'Maximal resolution was not found. Increase `max_tries` to search longer. ' - 'Or set `raise_if_not_maximum=False` if it is not necessary to locate the ' - 'maximum exactly.' + for i in range(refinements): + tmp_events = ( + events.fold(x, sizes={x: -1, '_x_aux': 2**i}) + .fold(y, sizes={y: -1, '_y_aux': 2**i}) + .bins.concat(['_x_aux', '_y_aux']) ) + min_counts_per_pixel = tmp_events.bins.size().min() + if min_counts_per_pixel.value > 0: + return tmp_events - return ( - sc.linspace( - coarse_x_bin_edges.dim, - coarse_x_bin_edges[0], - coarse_x_bin_edges[-1], - lower_nx, - ), - sc.linspace( - coarse_y_bin_edges.dim, - coarse_y_bin_edges[0], - coarse_y_bin_edges[-1], - lower_ny, - ), + raise ValueError( + 'Even at the coarsest pixel binning there are some pixels that have no events.' + ' Probably because the wavelength/time-of-arrival binning is too fine,' + ' or because number of bins in each image dimension is not a power of 2.' ) diff --git a/packages/essimaging/tests/imaging/tools/maximum_resolution_monitor_test.py b/packages/essimaging/tests/imaging/tools/maximum_resolution_monitor_test.py index a1b1a9a59..568a7a03e 100644 --- a/packages/essimaging/tests/imaging/tools/maximum_resolution_monitor_test.py +++ b/packages/essimaging/tests/imaging/tools/maximum_resolution_monitor_test.py @@ -21,17 +21,12 @@ def test_finds_maximum_resolution(): }, ) - x_be, y_be = maximum_resolution_achievable( - events, - sc.linspace('x', 0, 1, 2), - sc.linspace('x', 0, 1, 2), + rebinned = maximum_resolution_achievable( + events.bin(x=1024, y=1024), + ('x', 'y'), ) - assert len(x_be) == 11 - assert len(y_be) == 11 - assert x_be[0] == 0 - assert x_be[-1] == 1 - assert y_be[0] == 0 - assert y_be[-1] == 1 + assert rebinned.sizes['x'] == 8 + assert rebinned.sizes['y'] == 8 @pytest.mark.parametrize('seed', [0, 1, 2]) @@ -46,33 +41,27 @@ def test_finds_maximum_resolution_random(seed): 't': sc.ones(dims=['events'], shape=(n,)), }, ) - x_be, y_be = maximum_resolution_achievable( + events = events.bin(x=2**10, y=2**10) + del events.bins.coords['x'] + del events.bins.coords['y'] + rebinned = maximum_resolution_achievable( events, - sc.linspace('x', 0, 1, 2), - sc.linspace('y', 0, 1, 2), - # Need enough tries to be sure we find the optimum - max_tries=100, - ) - assert ( - events.bin(x=x_be, y=y_be, t=sc.linspace('t', 0, 2, 2)).bins.size().min().value - > 0 + ('x', 'y'), ) + + assert rebinned.bins.size().min().value > 0 assert ( - events.bin( - x=sc.linspace('x', 0, 1, len(x_be) + 1), - y=sc.linspace('y', 0, 1, len(y_be) + 1), - t=sc.linspace('t', 0, 2, 2), - ) + events.fold('x', sizes={'x': 2 * rebinned.sizes['x'], '_x_aux': -1}) + .fold('y', sizes={'y': 2 * rebinned.sizes['y'], '_y_aux': -1}) + .bins.concat(['_x_aux', '_y_aux']) .bins.size() .min() .value - == 0 - ) + ) == 0 def test_finds_maximum_resolution_binned_input(): - np.random.seed(0) - n = np.random.randint(1000, 100_000) + n = 100_000 events = sc.DataArray( sc.ones(dims=['events'], shape=(n,)), coords={ @@ -81,29 +70,59 @@ def test_finds_maximum_resolution_binned_input(): 't': sc.array(dims=['events'], values=np.random.random(n)), }, ) - events = events.bin(x=100, y=100, t=500) - del events.bins.coords['x'] - del events.bins.coords['y'] - x_be, y_be = maximum_resolution_achievable( - events, - sc.linspace('x', 0, 1, 2), - sc.linspace('y', 0, 1, 2), - # Need enough tries to be sure we find the optimum - max_tries=100, - ) - - events.bins.coords['x'] = sc.bins_like(events, sc.midpoints(events.coords['x'])) - events.bins.coords['y'] = sc.bins_like(events, sc.midpoints(events.coords['y'])) - events = events.bins.concat(['x', 'y']) + events = events.bin(x=128, y=128, t=500) + rebinned = maximum_resolution_achievable(events, ('x', 'y')) - assert events.bin(x=x_be, y=y_be).bins.size().min().value > 0 + assert rebinned.bins.size().min().value > 0 assert ( - events.bin( - x=sc.linspace('x', 0, 1, len(x_be) + 1), - y=sc.linspace('y', 0, 1, len(y_be) + 1), - ) + events.fold('x', sizes={'x': 2 * rebinned.sizes['x'], '_x_aux': -1}) + .fold('y', sizes={'y': 2 * rebinned.sizes['y'], '_y_aux': -1}) + .bins.concat(['_x_aux', '_y_aux']) .bins.size() .min() .value - == 0 + ) == 0 + + +def test_raises_if_bins_not_even(): + n = 100_000 + events = sc.DataArray( + sc.ones(dims=['events'], shape=(n,)), + coords={ + 'x': sc.array(dims=['events'], values=np.random.random(n)), + 'y': sc.array(dims=['events'], values=np.random.random(n)), + 't': sc.array(dims=['events'], values=np.random.random(n)), + }, + ) + events = events.bin(x=127, y=127, t=101) + with pytest.raises(ValueError, match='Input data must have an even number of bins'): + maximum_resolution_achievable(events, ('x', 'y')) + + +def test_raises_if_not_binned(): + n = 100_000 + events = sc.DataArray( + sc.ones(dims=['events'], shape=(n,)), + coords={ + 'x': sc.array(dims=['events'], values=np.random.random(n)), + 'y': sc.array(dims=['events'], values=np.random.random(n)), + 't': sc.array(dims=['events'], values=np.random.random(n)), + }, + ) + with pytest.raises(ValueError, match='Input data must be binned'): + maximum_resolution_achievable(events, ('x', 'y')) + + +def test_raises_if_reaches_coarsest_grid_without_success(): + n = 100 + events = sc.DataArray( + sc.ones(dims=['events'], shape=(n,)), + coords={ + 'x': sc.array(dims=['events'], values=np.random.random(n)), + 'y': sc.array(dims=['events'], values=np.random.random(n)), + 't': sc.array(dims=['events'], values=np.random.random(n)), + }, ) + events = events.bin(x=128, y=128, t=101) + with pytest.raises(ValueError, match='Even at the coarsest'): + maximum_resolution_achievable(events, ('x', 'y')) From 3f02742bfa4be702b58476e50d07a5e5f0d4a334 Mon Sep 17 00:00:00 2001 From: Johannes Kasimir Date: Mon, 13 Apr 2026 16:01:44 +0200 Subject: [PATCH 3/3] fix --- .../src/ess/imaging/tools/resolution.py | 41 +++++++++---------- .../tools/maximum_resolution_monitor_test.py | 40 +++++++++++++++--- 2 files changed, 54 insertions(+), 27 deletions(-) diff --git a/packages/essimaging/src/ess/imaging/tools/resolution.py b/packages/essimaging/src/ess/imaging/tools/resolution.py index 5e711edb5..1af31d93c 100644 --- a/packages/essimaging/src/ess/imaging/tools/resolution.py +++ b/packages/essimaging/src/ess/imaging/tools/resolution.py @@ -1,10 +1,16 @@ -from math import floor, log2 +from itertools import zip_longest import numpy as np import scipp as sc from numpy.typing import NDArray +def _all_divisors(n): + for i in range(1, n + 1): + if n % i == 0: + yield i + + def maximum_resolution_achievable( events: sc.DataArray, image_dims: tuple[str, str], @@ -22,41 +28,32 @@ def maximum_resolution_achievable( ------------- events: DataArray binned in at least the ``image_dims`` dimensions and optionally more. - The number of bins in each ``image_dims`` dimension ought to be power of ``2``. + The method works best when the number of bins in each ``image_dims`` dimension + is a power of ``2``. Returns ------------- - The event data array folded to the coarsest possible resolution. + The event data array folded to the finest resolution where there is at least + one event in every pixel. """ if events.bins is None: raise ValueError( - 'Input data must be binned to the finest possible pixel binning.' + 'Input data must be binned.' ' For best result, number of bins in each image dimension should' ' be a power of 2.' ) x, y = image_dims - for d in image_dims: - if events.sizes[d] % 2 == 1: - raise ValueError( - 'Input data must have an even number of bins in each image dimension.' - ' For best result, number of bins in each image dimension should' - ' be a power of 2.' - ) - - refinements = ( - min( - floor(log2(events.sizes[image_dims[0]])), - floor(log2(events.sizes[image_dims[1]])), - ) - + 1 - ) + xdivisors = list(_all_divisors(events.sizes[x])) + ydivisors = list(_all_divisors(events.sizes[y])) - for i in range(refinements): + for i, j in zip_longest(xdivisors, ydivisors): + i = xdivisors[-1] if i is None else i + j = ydivisors[-1] if j is None else j tmp_events = ( - events.fold(x, sizes={x: -1, '_x_aux': 2**i}) - .fold(y, sizes={y: -1, '_y_aux': 2**i}) + events.fold(x, sizes={x: -1, '_x_aux': i}) + .fold(y, sizes={y: -1, '_y_aux': j}) .bins.concat(['_x_aux', '_y_aux']) ) min_counts_per_pixel = tmp_events.bins.size().min() diff --git a/packages/essimaging/tests/imaging/tools/maximum_resolution_monitor_test.py b/packages/essimaging/tests/imaging/tools/maximum_resolution_monitor_test.py index 568a7a03e..6f17233b1 100644 --- a/packages/essimaging/tests/imaging/tools/maximum_resolution_monitor_test.py +++ b/packages/essimaging/tests/imaging/tools/maximum_resolution_monitor_test.py @@ -60,7 +60,8 @@ def test_finds_maximum_resolution_random(seed): ) == 0 -def test_finds_maximum_resolution_binned_input(): +def test_finds_maximum_resolution_time_binned_input(): + np.random.seed(0) n = 100_000 events = sc.DataArray( sc.ones(dims=['events'], shape=(n,)), @@ -84,7 +85,8 @@ def test_finds_maximum_resolution_binned_input(): ) == 0 -def test_raises_if_bins_not_even(): +def test_finds_maximum_resolution_non_even_number_of_bins(): + np.random.seed(0) n = 100_000 events = sc.DataArray( sc.ones(dims=['events'], shape=(n,)), @@ -94,12 +96,22 @@ def test_raises_if_bins_not_even(): 't': sc.array(dims=['events'], values=np.random.random(n)), }, ) - events = events.bin(x=127, y=127, t=101) - with pytest.raises(ValueError, match='Input data must have an even number of bins'): - maximum_resolution_achievable(events, ('x', 'y')) + events = events.bin(x=3**4 * 5, y=3**3 * 5, t=101) + rebinned = maximum_resolution_achievable(events, ('x', 'y')) + + assert rebinned.bins.size().min().value > 0 + assert ( + events.fold('x', sizes={'x': 45, '_x_aux': -1}) + .fold('y', sizes={'y': 15, '_y_aux': -1}) + .bins.concat(['_x_aux', '_y_aux']) + .bins.size() + .min() + .value + ) == 0 def test_raises_if_not_binned(): + np.random.seed(0) n = 100_000 events = sc.DataArray( sc.ones(dims=['events'], shape=(n,)), @@ -114,6 +126,7 @@ def test_raises_if_not_binned(): def test_raises_if_reaches_coarsest_grid_without_success(): + np.random.seed(0) n = 100 events = sc.DataArray( sc.ones(dims=['events'], shape=(n,)), @@ -126,3 +139,20 @@ def test_raises_if_reaches_coarsest_grid_without_success(): events = events.bin(x=128, y=128, t=101) with pytest.raises(ValueError, match='Even at the coarsest'): maximum_resolution_achievable(events, ('x', 'y')) + + +def test_one_dimension_denser_than_the_other(): + np.random.seed(0) + n = 100 + events = sc.DataArray( + sc.ones(dims=['events'], shape=(n,)), + coords={ + 'x': sc.array(dims=['events'], values=np.random.random(n)), + 'y': sc.array(dims=['events'], values=np.random.random(n)), + 't': sc.array(dims=['events'], values=np.random.random(n)), + }, + ) + events = events.bin(x=2, y=128, t=2) + rebinned = maximum_resolution_achievable(events, ('x', 'y')) + assert rebinned.bins.size().min().value > 0 + assert rebinned.sizes['y'] > 4