diff --git a/packages/essimaging/src/ess/imaging/tools/analysis.py b/packages/essimaging/src/ess/imaging/tools/analysis.py index f62dc3968..129767e38 100644 --- a/packages/essimaging/src/ess/imaging/tools/analysis.py +++ b/packages/essimaging/src/ess/imaging/tools/analysis.py @@ -57,7 +57,7 @@ def _is_non_bin_edges(coords: sc.Coords, coord_name: str) -> bool: def resample( image: sc.Variable | sc.DataArray, sizes: dict[str, int], - method: str | Callable = 'sum', + method: str | Callable | None = None, ) -> sc.Variable | sc.DataArray: """ Resample an image by folding it into blocks of specified sizes and applying a @@ -83,6 +83,8 @@ def resample( 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``. + By default, it will concatenate bins if the input data is binned, and sum + otherwise. Notes @@ -105,8 +107,15 @@ def resample( 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)) + dims_to_reduce = set(blocked.dims) - set(image.dims) + if method is None: + if image.is_binned: + out = blocked.bins.concat(dims_to_reduce) + else: + out = blocked.sum(dims_to_reduce) + else: + _method = getattr(sc, method) if isinstance(method, str) else method + out = _method(blocked, dims_to_reduce) if isinstance(image, sc.DataArray): # Restore the coordinates dropped by the `_method` if possible. @@ -127,7 +136,7 @@ def resample( def resize( image: sc.Variable | sc.DataArray, sizes: dict[str, int], - method: str | Callable = 'sum', + method: str | Callable | None = None, ) -> sc.Variable | sc.DataArray: """ Resize an image by folding it into blocks of specified sizes and applying a @@ -155,6 +164,8 @@ def resize( 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``. + By default, it will concatenate bins if the input data is binned, and sum + otherwise. Notes diff --git a/packages/essimaging/src/ess/imaging/types.py b/packages/essimaging/src/ess/imaging/types.py index 9dad68452..e17217a3b 100644 --- a/packages/essimaging/src/ess/imaging/types.py +++ b/packages/essimaging/src/ess/imaging/types.py @@ -27,6 +27,7 @@ NXsource = snx.NXsource DetectorLtotal = unwrap_t.DetectorLtotal +TofDetector = unwrap_t.TofDetector WavelengthDetector = unwrap_t.WavelengthDetector PulseStrideOffset = unwrap_t.PulseStrideOffset LookupTable = unwrap_t.LookupTable diff --git a/packages/essimaging/src/ess/odin/workflows.py b/packages/essimaging/src/ess/odin/workflows.py index 9ba3f823b..e52086d9c 100644 --- a/packages/essimaging/src/ess/odin/workflows.py +++ b/packages/essimaging/src/ess/odin/workflows.py @@ -5,6 +5,7 @@ """ import sciline +from scippneutron.conversion.tof import tof_from_wavelength from ess.reduce.unwrap import GenericUnwrapWorkflow, WavelengthLutMode @@ -13,12 +14,15 @@ BeamMonitor2, BeamMonitor3, BeamMonitor4, + CorrectedDetector, DarkBackgroundRun, LookupTableRelativeErrorThreshold, NeXusMonitorName, OpenBeamRun, PulseStrideOffset, + RunType, SampleRun, + TofDetector, ) from .masking import providers as masking_providers @@ -42,6 +46,15 @@ def default_parameters() -> dict: } +def compute_detector_tof(da: CorrectedDetector[RunType]) -> TofDetector[RunType]: + """ + Compute the time-of-flight of neutrons from their wavelength. + """ + return da.transform_coords( + "tof", graph={"tof": tof_from_wavelength}, keep_intermediate=False + ) + + def OdinWorkflow( wavelength_from: WavelengthLutMode = "analytical", **kwargs ) -> sciline.Pipeline: @@ -63,6 +76,7 @@ def OdinWorkflow( wavelength_from=wavelength_from, **kwargs, ) + workflow.insert(compute_detector_tof) for key, param in default_parameters().items(): workflow[key] = param return workflow diff --git a/packages/essimaging/tests/conftest.py b/packages/essimaging/tests/conftest.py new file mode 100644 index 000000000..2f25b0185 --- /dev/null +++ b/packages/essimaging/tests/conftest.py @@ -0,0 +1,19 @@ +# SPDX-License-Identifier: BSD-3-Clause +# Copyright (c) 2026 Scipp contributors (https://github.com/scipp) + +from pathlib import Path + +import pytest + + +def pytest_addoption(parser: pytest.Parser) -> None: + parser.addoption("--file-output", help='Output folder for reduced data') + + +@pytest.fixture +def output_folder(request: pytest.FixtureRequest) -> Path: + if (path := request.config.getoption("--file-output")) is not None: + out = Path(path) + out.mkdir(parents=True, exist_ok=True) + return out + return request.getfixturevalue("tmp_path") diff --git a/packages/essimaging/tests/imaging/tools/analysis_test.py b/packages/essimaging/tests/imaging/tools/analysis_test.py index ec488326b..579fa5378 100644 --- a/packages/essimaging/tests/imaging/tools/analysis_test.py +++ b/packages/essimaging/tests/imaging/tools/analysis_test.py @@ -4,7 +4,7 @@ import pytest import scipp as sc from ess import imaging as img -from scipp.testing import assert_identical +from scipp.testing import assert_allclose, assert_identical from scitiff.io import load_scitiff from ess.imaging.data import siemens_star_path @@ -14,12 +14,30 @@ def test_blockify() -> None: da = load_scitiff(siemens_star_path())["image"] blocks = img.tools.blockify(da, {'x': 4, 'y': 4}) assert len(blocks.dims) == len(da.dims) + 2 - assert {da.sizes['x'] // 4, da.sizes['y'] // 4, 4}.issubset(blocks.sizes.values()) + assert blocks.sizes['x'] == da.sizes['x'] // 4 + assert blocks.sizes['y'] == da.sizes['y'] // 4 + + +def test_blockify_binned_data() -> None: + da = sc.data.binned_xy(1000, 16, 16) + blocks = img.tools.blockify(da, {'x': 2, 'y': 2}) + assert len(blocks.dims) == len(da.dims) + 2 + assert blocks.sizes['x'] == da.sizes['x'] // 2 + assert blocks.sizes['y'] == da.sizes['y'] // 2 def test_resample() -> None: da = load_scitiff(siemens_star_path())["image"] resampled = img.tools.resample(da, sizes={'x': 2, 'y': 2}) + assert len(resampled.dims) == len(da.dims) + assert resampled.sizes['x'] == da.sizes['x'] // 2 + assert resampled.sizes['y'] == da.sizes['y'] // 2 + + +def test_resample_binned_data() -> None: + da = sc.data.binned_xy(1000, 16, 16) + resampled = img.tools.resample(da, sizes={'x': 2, 'y': 2}) + assert len(resampled.dims) == len(da.dims) assert resampled.sizes['x'] == da.sizes['x'] // 2 assert resampled.sizes['y'] == da.sizes['y'] // 2 @@ -107,7 +125,15 @@ def test_resize() -> None: resized = img.tools.resize(da, sizes={'x': 128, 'y': 128}) assert resized.sizes['x'] == 128 assert resized.sizes['y'] == 128 - assert sc.identical(resized.sum(), da.sum()) + assert_identical(resized.sum(), da.sum()) + + +def test_resize_binned_data() -> None: + da = sc.data.binned_xy(1000, 16, 16) + resized = img.tools.resize(da, sizes={'x': 8, 'y': 8}) + assert resized.sizes['x'] == 8 + assert resized.sizes['y'] == 8 + assert_allclose(resized.sum(), da.sum()) def test_resize_mean() -> None: diff --git a/packages/essimaging/tests/odin/data_reduction_test.py b/packages/essimaging/tests/odin/data_reduction_test.py index 9b5fc03b7..96aba1cb4 100644 --- a/packages/essimaging/tests/odin/data_reduction_test.py +++ b/packages/essimaging/tests/odin/data_reduction_test.py @@ -1,21 +1,27 @@ # SPDX-License-Identifier: BSD-3-Clause # Copyright (c) 2025 Scipp contributors (https://github.com/scipp) +from pathlib import Path + import ess.odin.data # noqa: F401 import pytest import sciline as sl +import scipp as sc from ess import odin from ess.odin.beamline import choppers as odin_choppers +from ess.imaging import tools from ess.imaging.types import ( Filename, LookupTableFilename, + MaskingRules, NeXusDetectorName, NXsource, OpenBeamRun, Position, RawDetector, SampleRun, + TofDetector, WavelengthDetector, ) from ess.reduce import unwrap @@ -64,3 +70,43 @@ def test_can_compute_wavelength(run_type, wavelength_mode): da = wf.compute(WavelengthDetector[run_type]) assert "wavelength" in da.bins.coords + + +@pytest.mark.parametrize("run_type", [SampleRun, OpenBeamRun]) +@pytest.mark.parametrize("wavelength_mode", ["file", "analytical"]) +def test_can_compute_tof(run_type, wavelength_mode): + wf = _make_workflow(wavelength_mode) + wf[MaskingRules] = {} + da = wf.compute(TofDetector[run_type]) + + assert "tof" in da.bins.coords + + +def test_publish_reduced_scitiff(output_folder: Path): + wf = _make_workflow("analytical") + wf[MaskingRules] = {} + new_sizes = {'dim_0': 64, 'dim_1': 64} + tbins = sc.linspace('tof', 1.3e4, 1.5e5, 257, unit='us') + + sample = wf.compute(TofDetector[SampleRun]).drop_coords('detector_number') + res_sample = tools.resample(sample, sizes=new_sizes) + num = res_sample.hist(tof=tbins) + + openbeam = wf.compute(TofDetector[OpenBeamRun]).drop_coords('detector_number') + res_openbeam = tools.resample(openbeam, sizes=new_sizes) + den = res_openbeam.hist(tof=tbins) + + normed = num / den + + to_scitiff = normed.rename_dims(dim_0='y', dim_1='x', tof='t').drop_coords( + 'position' + ) + + assert "tof" in to_scitiff.coords + assert "Ltotal" in to_scitiff.coords + + from scitiff.io import save_scitiff + + save_scitiff( + to_scitiff, output_folder / 'bragg_edge_iron_normalized_16x16x256.tiff' + ) diff --git a/packages/essnmx/src/ess/nmx/types.py b/packages/essnmx/src/ess/nmx/types.py index 713e55cef..077e4ed37 100644 --- a/packages/essnmx/src/ess/nmx/types.py +++ b/packages/essnmx/src/ess/nmx/types.py @@ -4,18 +4,18 @@ import h5py import numpy as np -import sciline as sl import scipp as sc import scippnexus as snx from scippneutron.metadata import RadiationProbe, SourceType from ess.reduce.nexus import types as nexus_types -from ess.reduce.unwrap.types import LookupTable +from ess.reduce.unwrap import types as unwrap_types from ._display_helper import to_datagroup RunType = nexus_types.RunType SampleRun = nexus_types.SampleRun +TofDetector = unwrap_types.TofDetector class Compression(enum.StrEnum): @@ -275,13 +275,9 @@ class NMXLauetof: definitions: Literal['NXlauetof'] = 'NXlauetof' instrument: NMXInstrument sample: NMXSampleMetadata - lookup_table: LookupTable | None = None + lookup_table: unwrap_types.LookupTable | None = None reducer: NMXProgram = field(default_factory=NMXProgram) "Information of the reduction software." def to_datagroup(self) -> sc.DataGroup: return to_datagroup(self) - - -class TofDetector(sl.Scope[RunType, sc.DataArray], sc.DataArray): - """Detector data with time-of-flight coordinate.""" diff --git a/packages/essreduce/src/ess/reduce/unwrap/__init__.py b/packages/essreduce/src/ess/reduce/unwrap/__init__.py index 997e47573..88df2fbb7 100644 --- a/packages/essreduce/src/ess/reduce/unwrap/__init__.py +++ b/packages/essreduce/src/ess/reduce/unwrap/__init__.py @@ -33,6 +33,8 @@ LookupTableRelativeErrorThreshold, MonitorLtotal, PulseStrideOffset, + TofDetector, + TofMonitor, WavelengthDetector, WavelengthLutMode, WavelengthMonitor, @@ -64,6 +66,8 @@ "SimulationSeed", "SourceBounds", "TimeResolution", + "TofDetector", + "TofMonitor", "WavelengthDetector", "WavelengthLutMode", "WavelengthMonitor", diff --git a/packages/essreduce/src/ess/reduce/unwrap/types.py b/packages/essreduce/src/ess/reduce/unwrap/types.py index d6a44f57c..0846523dc 100644 --- a/packages/essreduce/src/ess/reduce/unwrap/types.py +++ b/packages/essreduce/src/ess/reduce/unwrap/types.py @@ -92,6 +92,14 @@ class MonitorLtotal(sl.Scope[RunType, MonitorType, sc.Variable], sc.Variable): """Total path length of neutrons from source to monitor.""" +class TofDetector(sl.Scope[RunType, sc.DataArray], sc.DataArray): + """Detector data with time-of-flight coordinate.""" + + +class TofMonitor(sl.Scope[RunType, MonitorType, sc.DataArray], sc.DataArray): + """Monitor data with time-of-flight coordinate.""" + + class WavelengthDetector(sl.Scope[RunType, sc.DataArray], sc.DataArray): """Detector data with wavelength coordinate."""