diff --git a/src/hextools/exporters/tomography.py b/src/hextools/exporters/tomography.py index 5495638..3238adb 100644 --- a/src/hextools/exporters/tomography.py +++ b/src/hextools/exporters/tomography.py @@ -8,6 +8,7 @@ from bluesky_tiled_plugins.clients.bluesky_run import BlueskyRunV3 from hextools.exporters.utils import ( + average_per_projection, confirm_required_start_fields, create_baseline_dataset, find_hdf_datasets, @@ -76,6 +77,24 @@ def export_tomo_flyscan_to_h5( for name, (relpath, dataset, _) in primary.items(): links["data" if single_det and name in det_names else name] = (relpath, dataset) + # When averaging, every PandA capture (e.g. rotation_angle) is written as the + # mean of each projection's images; the raw link is kept. + images_to_average = run_start.get("images_to_average", 1) + num_projections = run_start.get("num_points") + if not isinstance(images_to_average, int): + raise ValueError(f"images_to_average is not an int: {images_to_average!r}") + averaged = {} + for name, (_, _, ndim) in primary.items(): + if name in det_names or (images_to_average == 1 and ndim == 1): + continue + if not isinstance(num_projections, int): + raise ValueError(f"num_points is not an int: {num_projections!r}") + averaged[name] = average_per_projection( + run["primary"][name].read(), num_projections, images_to_average + ) + links[f"{name}_raw"] = links.pop(name) + logger.info(f"Averaging {name} over {images_to_average} images per projection") + if dark_flat_run is not None: for stream_name in ("dark", "flat"): found = find_hdf_datasets(dark_flat_run, stream_name, scan_id_dir_path) @@ -99,6 +118,9 @@ def export_tomo_flyscan_to_h5( for link_name, (relpath, dataset) in links.items(): data_grp[link_name] = h5py.ExternalLink(relpath, dataset) logger.info(f"Linked /entry/data/{link_name} -> {relpath}:{dataset}") + for name, values in averaged.items(): + data_grp.create_dataset(name, data=values) + logger.info(f"Wrote /entry/data/{name} ({len(values)} projections)") logger.info("Writing baseline data") create_baseline_dataset(run, h5file) diff --git a/src/hextools/exporters/utils.py b/src/hextools/exporters/utils.py index 6bc32c3..ef36203 100644 --- a/src/hextools/exporters/utils.py +++ b/src/hextools/exporters/utils.py @@ -51,6 +51,28 @@ def confirm_required_start_fields(run_start: dict[str, Any]) -> tuple[str, str, return tuple(field_values) # type: ignore +def average_per_projection( + values: np.ndarray, num_projections: int, images_to_average: int +) -> np.ndarray: + """Average a PandA capture to one value per projection. + + With ``images_to_average`` > 1 the PandA records every exposure. ophyd-async + describes each event as ``images_to_average`` values, so Tiled shapes the + capture (projections, images_to_average); any other shape is refused rather + than guessed. The averaged detector frame belongs at the mean of its + exposures' values. + """ + values = np.asarray(values) + if images_to_average == 1 and values.shape == (num_projections,): + return values + if values.shape == (num_projections, images_to_average): + return values.mean(axis=1) + raise ValueError( + f"Capture shape {values.shape} is not (num_projections, images_to_average)" + f" = ({num_projections}, {images_to_average})" + ) + + def create_baseline_dataset(run: BlueskyRunV3, h5file: h5py.File): if "baseline" not in run: logger.warning("No baseline data found in the run.") diff --git a/src/hextools/flyers.py b/src/hextools/flyers.py index adf2772..c8eac81 100644 --- a/src/hextools/flyers.py +++ b/src/hextools/flyers.py @@ -192,7 +192,8 @@ def construct_fly_info_models( pulse_step = travel_counts // (num_pulses - 1) else: pulse_width = max_exposure_time - pulse_step = move_time / num_pulses + # The first pulse fires at the start, so num_pulses - 1 steps reach the stop. + pulse_step = move_time / (num_pulses - 1) if num_pulses > 1 else move_time flyer_info = SingleAxisFlyscanInfo( start=start_in_counts, diff --git a/tests/test_export_averaging.py b/tests/test_export_averaging.py new file mode 100644 index 0000000..4e6c843 --- /dev/null +++ b/tests/test_export_averaging.py @@ -0,0 +1,43 @@ +"""Averaging PandA captures down to one value per projection for the .nxs export.""" + +import numpy as np +import pytest + +from hextools.exporters.utils import average_per_projection + +OFFSETS = np.array([-4, -2, 0, 2, 4]) + + +def test_averages_each_projections_captures(): + # 19 projections x 5 captures, as Tiled shapes a PandA dataset with + # collections_per_event=5 (events first). + centres = np.linspace(0, 180, 19) + captures = centres[:, None] + OFFSETS + np.testing.assert_allclose(average_per_projection(captures, 19, 5), centres) + + +def test_square_capture_averages_over_images(): + centres = np.linspace(0, 180, 5) + captures = centres[:, None] + OFFSETS + np.testing.assert_allclose(average_per_projection(captures, 5, 5), centres) + + +def test_refuses_captures_first(): + captures = (np.linspace(0, 180, 19)[:, None] + OFFSETS).T + with pytest.raises(ValueError, match="images_to_average"): + average_per_projection(captures, 19, 5) + + +def test_leaves_one_capture_per_projection_alone(): + angles = np.linspace(0, 180, 19) + np.testing.assert_array_equal(average_per_projection(angles, 19, 1), angles) + + +def test_refuses_a_shape_that_does_not_match(): + with pytest.raises(ValueError, match="images_to_average"): + average_per_projection(np.zeros((19, 4)), 19, 5) + + +def test_refuses_a_flattened_capture_when_averaging(): + with pytest.raises(ValueError, match="images_to_average"): + average_per_projection(np.zeros(95), 19, 5) diff --git a/tests/test_flyers.py b/tests/test_flyers.py new file mode 100644 index 0000000..67e13c8 --- /dev/null +++ b/tests/test_flyers.py @@ -0,0 +1,72 @@ +"""Where a single-axis flyscan's triggers land, for both trigger modes.""" + +import numpy as np +import pytest +from ophyd_async.fastcs.panda import PandaPcompDirection + +from hextools.exporters.utils import average_per_projection +from hextools.flyers import construct_fly_info_models + +RES = 0.001 # deg per encoder count + + +def _trigger_angles(flyer_info, motor_info, start): + """Angles the PandA records at each trigger, relative to ``start``.""" + k = np.arange(flyer_info.num_pulses) + if not flyer_info.time_based: + sign = 1 if flyer_info.direction == PandaPcompDirection.POSITIVE else -1 + counts = flyer_info.start + sign * k * flyer_info.pulse_step + return counts * flyer_info.position_scale + flyer_info.position_offset + # Constant velocity between start_position and end_position over time_for_move. + velocity = ( + motor_info.end_position - motor_info.start_position + ) / motor_info.time_for_move + positions = motor_info.start_position + velocity * k * flyer_info.pulse_step + return positions - start + + +def _models(num_projections, images_to_average, time_based, start=0.0, stop=180.0): + return construct_fly_info_models( + num_pulses=num_projections * images_to_average, + max_exposure_time=0.01, + start_position=start, + stop_position=stop, + encoder_resolution=RES, + max_motor_velocity=10.0, + current_position=0.0, + current_enc_position=0, + time_based=time_based, + ) + + +@pytest.mark.parametrize("time_based", [False, True]) +@pytest.mark.parametrize(("start", "stop"), [(0.0, 180.0), (180.0, 0.0)]) +def test_triggers_span_start_to_stop(time_based, start, stop): + flyer_info, motor_info = _models(19, 1, time_based, start, stop) + angles = _trigger_angles(flyer_info, motor_info, start) + np.testing.assert_allclose(angles, np.linspace(0, stop - start, 19), atol=1e-9) + + +# 19 x 5 triggers over 0-188 deg: 188000 counts split evenly into 94 steps. +@pytest.mark.parametrize("time_based", [False, True]) +def test_averaged_triggers_span_start_to_stop(time_based): + flyer_info, motor_info = _models(19, 5, time_based, stop=188.0) + angles = _trigger_angles(flyer_info, motor_info, 0.0) + np.testing.assert_allclose(angles, np.linspace(0, 188, 95), atol=1e-9) + + +def test_averaged_projection_angle_is_mean_of_its_images(): + flyer_info, motor_info = _models(19, 5, time_based=True) + angles = _trigger_angles(flyer_info, motor_info, 0.0) + groups = angles.reshape(19, 5) + np.testing.assert_allclose( + average_per_projection(groups, 19, 5), groups.sum(axis=1) / 5 + ) + assert angles[0] == pytest.approx(0.0) + assert angles[-1] == pytest.approx(180.0) + + +def test_single_time_based_pulse_fires_at_start(): + flyer_info, motor_info = _models(1, 1, time_based=True) + assert flyer_info.num_pulses == 1 + np.testing.assert_allclose(_trigger_angles(flyer_info, motor_info, 0.0), [0.0])