"""Helpers for deriving per-frame maximum peak radius from peakfinding outputs."""
from __future__ import annotations
import numpy as np
import h5py
from coseda.io import config_to_paths, read_config
from coseda.logging_utils import log_start
def _centers_from_config(config):
"""Return centre definition from INI config as ('h5','h5') or (x0,y0)."""
params = config["Parameters"] if "Parameters" in config else {}
x_raw = params.get("peakfinding_x0", params.get("pf9_x0", "h5"))
y_raw = params.get("peakfinding_y0", params.get("pf9_y0", "h5"))
x_text = str(x_raw).strip()
y_text = str(y_raw).strip()
if x_text.lower() == "h5" and y_text.lower() == "h5":
return "h5", "h5"
try:
return float(x_text), float(y_text)
except ValueError:
# Fall back to HDF5 centres if values are malformed.
return "h5", "h5"
def _compute_chunk_maxres(n_peaks_chunk, peak_x_chunk, peak_y_chunk, center_x_chunk, center_y_chunk):
"""Compute max radius for one chunk of frames."""
out = np.zeros(len(n_peaks_chunk), dtype=np.float32)
n_slots = peak_x_chunk.shape[1]
for idx in range(len(n_peaks_chunk)):
n_peaks = int(n_peaks_chunk[idx])
if n_peaks <= 0:
continue
n_peaks = min(n_peaks, n_slots)
cx = float(center_x_chunk[idx])
cy = float(center_y_chunk[idx])
if not np.isfinite(cx) or not np.isfinite(cy):
continue
x_vals = peak_x_chunk[idx, :n_peaks]
y_vals = peak_y_chunk[idx, :n_peaks]
valid = np.isfinite(x_vals) & np.isfinite(y_vals)
if not np.any(valid):
continue
dx = x_vals[valid] - cx
dy = y_vals[valid] - cy
dist_sq = dx * dx + dy * dy
if dist_sq.size:
out[idx] = float(np.sqrt(np.max(dist_sq)))
return out
[docs]
def write_maxres_dataset(
h5file_path,
center_x="h5",
center_y="h5",
logfile_path=None,
chunk_size=2048,
):
"""
Create/update ``/entry/data/maxres`` with furthest-peak distance from centre.
Distances are in pixels, one value per frame.
"""
with h5py.File(h5file_path, "r+") as h5file:
data_group = h5file.get("entry/data")
if data_group is None:
raise KeyError("Missing group '/entry/data' in HDF5 file.")
required = ("nPeaks", "peakXPosRaw", "peakYPosRaw")
missing = [name for name in required if name not in data_group]
if missing:
raise KeyError(f"Missing dataset(s) in /entry/data: {', '.join(missing)}")
n_peaks_ds = data_group["nPeaks"]
peak_x_ds = data_group["peakXPosRaw"]
peak_y_ds = data_group["peakYPosRaw"]
if n_peaks_ds.ndim != 1:
raise ValueError("Dataset '/entry/data/nPeaks' must be 1D.")
if peak_x_ds.ndim != 2 or peak_y_ds.ndim != 2:
raise ValueError("Datasets '/entry/data/peakXPosRaw' and '/entry/data/peakYPosRaw' must be 2D.")
num_frames = int(n_peaks_ds.shape[0])
if peak_x_ds.shape[0] < num_frames or peak_y_ds.shape[0] < num_frames:
raise ValueError("Peak position datasets are shorter than nPeaks.")
use_h5_centers = (
isinstance(center_x, str)
and isinstance(center_y, str)
and center_x.lower() == "h5"
and center_y.lower() == "h5"
)
center_x_ds = None
center_y_ds = None
if use_h5_centers:
if "center_x" not in data_group or "center_y" not in data_group:
raise KeyError("Requested HDF5 centres, but '/entry/data/center_x' or '/entry/data/center_y' is missing.")
center_x_ds = data_group["center_x"]
center_y_ds = data_group["center_y"]
if center_x_ds.shape[0] < num_frames or center_y_ds.shape[0] < num_frames:
raise ValueError("Center datasets are shorter than nPeaks.")
else:
center_x = float(center_x)
center_y = float(center_y)
if "maxres" in data_group and data_group["maxres"].shape != (num_frames,):
del data_group["maxres"]
if "maxres" not in data_group:
data_group.create_dataset("maxres", shape=(num_frames,), dtype=np.float32)
maxres_ds = data_group["maxres"]
maxres_ds.attrs["units"] = "pixel"
maxres_ds.attrs["long_name"] = "Maximum peak radius from center"
if num_frames == 0:
log_start(logfile_path, "Updated /entry/data/maxres (0 frames).")
return
for start in range(0, num_frames, chunk_size):
end = min(start + chunk_size, num_frames)
n_peaks_chunk = np.asarray(n_peaks_ds[start:end], dtype=np.int64)
peak_x_chunk = np.asarray(peak_x_ds[start:end, :], dtype=np.float64)
peak_y_chunk = np.asarray(peak_y_ds[start:end, :], dtype=np.float64)
if use_h5_centers:
center_x_chunk = np.asarray(center_x_ds[start:end], dtype=np.float64)
center_y_chunk = np.asarray(center_y_ds[start:end], dtype=np.float64)
else:
chunk_len = end - start
center_x_chunk = np.full(chunk_len, center_x, dtype=np.float64)
center_y_chunk = np.full(chunk_len, center_y, dtype=np.float64)
maxres_ds[start:end] = _compute_chunk_maxres(
n_peaks_chunk=n_peaks_chunk,
peak_x_chunk=peak_x_chunk,
peak_y_chunk=peak_y_chunk,
center_x_chunk=center_x_chunk,
center_y_chunk=center_y_chunk,
)
log_start(logfile_path, "Updated /entry/data/maxres.")
[docs]
def write_maxres_from_config(configfile, logfile_path=None, chunk_size=2048):
"""Resolve centres from an INI file and write ``/entry/data/maxres``."""
_, _, _, inferred_logfile_path, _, h5file_path = config_to_paths(configfile)
if h5file_path is None:
raise ValueError(f"Could not resolve HDF5 path for config: {configfile}")
if logfile_path is None:
logfile_path = inferred_logfile_path
config = read_config(configfile)
center_x, center_y = _centers_from_config(config)
write_maxres_dataset(
h5file_path=h5file_path,
center_x=center_x,
center_y=center_y,
logfile_path=logfile_path,
chunk_size=chunk_size,
)