Tomo flyscan: correct time-based trigger spacing and average PandA angles per projection

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)

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.")

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,

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)

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])