Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
76 changes: 72 additions & 4 deletions src/ess/imaging/tools/analysis.py
Original file line number Diff line number Diff line change
Expand Up @@ -37,6 +37,23 @@ def blockify(
return out


def _is_1d_sorted_bin_edges(coords: sc.Coords, coord_name: str) -> bool:
return (
(coord := coords[coord_name]).ndim == 1
and coords.is_edges(coord_name)
and (
bool(sc.issorted(coord, dim=coord.dim, order='ascending'))
or bool(sc.issorted(coord, dim=coord.dim, order='descending'))
)
)


def _is_non_bin_edges(coords: sc.Coords, coord_name: str) -> bool:
return not any(
coords.is_edges(coord_name, dim=dim) for dim in coords[coord_name].dims
)


def resample(
image: sc.Variable | sc.DataArray,
sizes: dict[str, int],
Expand All @@ -56,21 +73,54 @@ def resample(
sizes:
A dictionary specifying the block sizes for each dimension.
For example, ``{'x': 4, 'y': 4}`` will create blocks of size 4x4.
Any dimensions with size ``1`` will be ignored.
If all sizes are set to ``1``,
it will not apply ``method`` and return a copy of the input resampling image.
method:
The reduction method to apply to the blocks. This can be a string referring to
any valid Scipp reduction method, such as 'sum', 'mean', 'max', etc.
Alternatively, a custom reduction function can be provided. The function
signature should accept a ``scipp.Variable`` or ``scipp.DataArray`` as first
argument and a set of dimensions to reduce over as second argument. The
function should return a ``scipp.Variable`` or ``scipp.DataArray``.


Notes
-----
If the image is a ``scipp.DataArray``,
bin edges in the resampling dimensions
will be preserved if they are 1-dimensional and sorted (they are dropped otherwise).
New bin edges will be chosen according to the resampling sizes.
For midpoint coordinates, new coordinates will be average values
of each resampled block.

.. warning::
The coordinates in the output may not be correct
if they are not sorted or not linear.

"""
# Filter the resample sizes first.
sizes = {dim: size for dim, size in sizes.items() if size != 1}
if not sizes:
return image.copy(deep=False)

blocked = blockify(image, sizes=sizes)
_method = getattr(sc, method) if isinstance(method, str) else method
out = _method(blocked, set(blocked.dims) - set(image.dims))
if 'position' in blocked.coords:
out.coords['position'] = blocked.coords['position'].mean(
set(blocked.dims) - set(image.dims)
)

if isinstance(image, sc.DataArray):
# Restore the coordinates dropped by the `_method` if possible.
_dropped_cnames = set(image.coords.keys()) - set(out.coords.keys())

for name in _dropped_cnames:
coord = image.coords[name]
if _is_1d_sorted_bin_edges(image.coords, name):
out.coords[name] = coord[coord.dim, :: sizes[coord.dim]]
elif _is_non_bin_edges(image.coords, name):
folded_coord = blocked.coords[name]
reduced_dims = set(folded_coord.dims) - set(coord.dims)
out.coords[name] = folded_coord.mean(reduced_dims)

return out


Expand All @@ -95,13 +145,31 @@ def resize(
The original sizes should be divisible by the specified sizes.
For example, ``{'x': 128, 'y': 128}`` will create an output image of size
128x128.
Any dimensions with same sizes to the resizing image will be ignored.
If the output image sizes will be same as the input image sizes,
it will not apply the ``method`` and return a copy of the input image.
method:
The reduction method to apply to the blocks. This can be a string referring to
any valid Scipp reduction method, such as 'sum', 'mean', 'max', etc.
Alternatively, a custom reduction function can be provided. The function
signature should accept a ``scipp.Variable`` or ``scipp.DataArray`` as first
argument and a set of dimensions to reduce over as second argument. The
function should return a ``scipp.Variable`` or ``scipp.DataArray``.


Notes
-----
If the image is a ``scipp.DataArray``,
bin edges in the resizing dimensions
will be preserved if they are 1-dimensional and sorted.
New bin edges will be chosen according to the resizing sizes.
For other coordinates, new coordinates will be average values
of each resized blocks.

.. warning::
The coordinates in the output may not be correct
if they are not sorted or not linear.

"""
block_sizes = {}
for dim, size in sizes.items():
Expand Down
95 changes: 94 additions & 1 deletion tests/imaging/tools/analysis_test.py
Original file line number Diff line number Diff line change
Expand Up @@ -3,6 +3,7 @@
import numpy as np
import pytest
import scipp as sc
from scipp.testing import assert_identical
from scitiff.io import load_scitiff

from ess import imaging as img
Expand All @@ -23,7 +24,14 @@ def test_resample() -> None:
assert resampled.sizes['y'] == da.sizes['y'] // 2


def test_resample_with_position_coord() -> None:
def test_resample_with_sizes_1() -> None:
da = load_scitiff(siemens_star_path())["image"]
resampled = img.tools.resample(da, sizes=dict.fromkeys(da.dims, 1))
assert_identical(resampled, da)
assert da is not resampled


def test_resample_with_2d_position_coord() -> None:
da = load_scitiff(siemens_star_path())["image"]
vectors = np.random.randn(*da.shape[1:], 3)
da.coords['position'] = sc.vectors(dims=['x', 'y'], values=vectors)
Expand All @@ -35,6 +43,50 @@ def test_resample_with_position_coord() -> None:
)


def test_resample_keep_coordinate() -> None:
da = load_scitiff(siemens_star_path())["image"]
# Overwrite the coordinate 't' to be a bin-edge.
da.coords['t'] = sc.array(dims=['t'], values=[1, 2, 3])
resampled = img.tools.resample(da, sizes={'x': 2, 'y': 2, 't': 3})
expected_t = sc.array(dims=['t'], values=[2.0], unit='dimensionless')
assert_identical(expected_t, resampled.coords['t'])


def test_resample_keep_bin_edge_coordinate() -> None:
da = load_scitiff(siemens_star_path())["image"]
# Overwrite the coordinate 't' to be a bin-edge.
da.coords['t'] = sc.array(dims=['t'], values=[0, 1, 2, 3])
resampled = img.tools.resample(da, sizes={'x': 2, 'y': 2, 't': 3})
expected_t = sc.array(dims=['t'], values=[0, 3], unit='dimensionless')
assert_identical(expected_t, resampled.coords['t'])


def test_resample_unsorted_bin_edge_coordinate_ignored() -> None:
da = load_scitiff(siemens_star_path())["image"]
# Overwrite the coordinate 't' to be a bin-edge.
da.coords['t'] = sc.array(dims=['t'], values=[0, 1, 3, 2])
resampled = img.tools.resample(da, sizes={'x': 2, 'y': 2, 't': 3})
assert 't' not in resampled.coords


def test_resample_2d_partial_bin_edge_coordinate_ignored() -> None:
data = sc.array(dims=['x', 'y'], values=[[1, 2], [3, 4]])
pos = sc.array(dims=['x', 'y'], values=[[1, 2, 3], [3, 4, 5]])
da = sc.DataArray(data=data, coords={'pos': pos})
assert not da.coords.is_edges('pos', 'x')
assert da.coords.is_edges('pos', 'y')
resampled = img.tools.resample(da, sizes={'x': 2, 'y': 2})
assert 'pos' not in resampled.coords


def test_resample_2d_coordinate_not_dropped_if_not_changed() -> None:
data = sc.array(dims=['x', 'y'], values=[[1, 2], [3, 4]])
pos = sc.array(dims=['x', 'y'], values=[[1, 2, 3], [3, 4, 5]])
da = sc.DataArray(data=data, coords={'pos': pos})
resampled = img.tools.resample(da, sizes={'x': 1, 'y': 1})
assert_identical(pos, resampled.coords['pos'])


def test_resample_mean() -> None:
da = load_scitiff(siemens_star_path())["image"]
resampled = img.tools.resample(da, sizes={'x': 2, 'y': 2}, method='mean')
Expand Down Expand Up @@ -73,6 +125,47 @@ def test_resize_callable() -> None:
assert resized.sizes['y'] == 256


def test_resize_same_sizes() -> None:
da = load_scitiff(siemens_star_path())["image"]
resized = img.tools.resize(da, sizes=da.sizes)
assert_identical(resized, da)
assert da is not resized


def test_resize_keep_coordinate() -> None:
da = load_scitiff(siemens_star_path())["image"]
# Overwrite the coordinate 't' to be a bin-edge.
da.coords['t'] = sc.array(dims=['t'], values=[1, 2, 3])
resampled = img.tools.resize(da, sizes={'t': 1})
expected_t = sc.array(dims=['t'], values=[2.0], unit='dimensionless')
assert_identical(expected_t, resampled.coords['t'])


def test_resize_keep_bin_edge_coordinate() -> None:
da = load_scitiff(siemens_star_path())["image"]
# Overwrite the coordinate 't' to be a bin-edge.
da.coords['t'] = sc.array(dims=['t'], values=[0, 1, 2, 3])
resized = img.tools.resize(da, sizes={'x': 256, 'y': 256, 't': 1})
expected_t = sc.array(dims=['t'], values=[0, 3], unit='dimensionless')
assert_identical(expected_t, resized.coords['t'])


def test_resize_unsorted_bin_edge_ignored() -> None:
da = load_scitiff(siemens_star_path())["image"]
# Overwrite the coordinate 't' to be a bin-edge.
da.coords['t'] = sc.array(dims=['t'], values=[0, 2, 1, 3])
resized = img.tools.resize(da, sizes={'x': 256, 'y': 256, 't': 1})
assert 't' not in resized.coords


def test_resize_2d_coordinate_not_dropped_if_not_changed() -> None:
data = sc.array(dims=['x', 'y'], values=[[1, 2], [3, 4]])
pos = sc.array(dims=['x', 'y'], values=[[1, 2, 3], [3, 4, 5]])
da = sc.DataArray(data=data, coords={'pos': pos})
resampled = img.tools.resize(da, sizes={'x': 2, 'y': 2})
assert_identical(pos, resampled.coords['pos'])


def test_resize_bad_size_requested_raises():
da = load_scitiff(siemens_star_path())["image"]
with pytest.raises(ValueError, match="Size of dimension 'x' .* is not divisible"):
Expand Down
Loading