diff --git a/docs/developer-guide/internals/description.md b/docs/developer-guide/internals/description.md index 264305d9..5896b4ba 100644 --- a/docs/developer-guide/internals/description.md +++ b/docs/developer-guide/internals/description.md @@ -154,20 +154,31 @@ echo flag is set. Pulserver's reconstruction proxy reads the sequence file a series was played from and sends its description after the MRD XML header and before the first -acquisition, as one `TEXT` message: a JSON object under the key -`pulserver_sequence_description`, holding one entry per file of the sequence -chain. {func}`~blochsim.read_mrd_description` reads the stream up to that -message, returns the XML header and the descriptions in play order, and leaves -the stream at the first acquisition. A message other than a filename, config, -header or text before the description is refused. - -Each entry holds `subsequence_index`, `tr_duration_us`, `rf_raster_time_s`, -the events as three base64 arrays -- `type` (int32), `timestamp_us` -(float64) and `params` (float32, seven per event, in the wire order above) -- -and `rf_definitions`, each with its bandwidth, bands, and Pulseq-compressed -magnitude, phase and time shapes. The stream carries no shim definitions, no -per-sample slice-select gradient, no crusher dephasing or voxel size, and no -event actions. +acquisition, as ISMRMRD waveforms, which every MRD writer keeps with the series: +per file of the sequence chain a header (`waveform_id` 1024) and an event list +(1025), and one waveform per RF definition (1026). Each carries one vector of +32-bit words led by its length, over as many channels as keep a channel within +65 535 samples; floats are bit-cast and float64 values take two words, low then +high. + +- Header: format version, subsequence index, number of subsequences, number of + events, number of RF definitions, `tr_duration_us`, `rf_raster_time_s`. +- Event list: subsequence index, number of events, then ten words an event: + type, `timestamp_us` (float64), and the seven float32 parameters in the wire + order above. +- RF definition: subsequence index, id, bandwidth, number of bands, eight band + offsets, band bandwidth, a word left 0, then for the magnitude, phase and time + shapes whether each is given and its uncompressed and stored sample counts, + then the samples of those given. + +{func}`~blochsim.read_mrd_description` reads the stream up to the last of +them, passing the physiological waveforms over, returns the XML header and the +descriptions in play order, and leaves the stream at the first acquisition; +{func}`~blochsim.description_from_waveforms` reads them off the waveforms of a +stored series. A message other than a filename, config, header, text or +waveform before the description is refused. The description +carries no shim definitions, no per-sample slice-select gradient, no crusher +dephasing or voxel size, and no event actions. ## Gradients diff --git a/docs/explanations/description.md b/docs/explanations/description.md index bebfc70d..d4c41803 100644 --- a/docs/explanations/description.md +++ b/docs/explanations/description.md @@ -56,7 +56,7 @@ wire format are in {doc}`../developer-guide/internals/description`. - **Decodes the description pulserver sends on an MRD stream, one per file of the sequence chain, in play order.** :::{dropdown} Show code and tests - Code: {func}`~blochsim.read_mrd_description`. Tests: *a written stream reads back into the rows it was built from*; *every file of a chain is read in play order*; *a stream whose data arrives first is refused* (`test_mrd_reader.py`). + Code: {func}`~blochsim.read_mrd_description`, {func}`~blochsim.description_from_waveforms`. Tests: *a written stream reads back into the rows it was built from*; *every file of a chain is read in play order*; *a stream whose data arrives first is refused*; *the physiological waveforms ahead of the description are passed over*; *a kept series reads its description from its waveforms*; *a pulse is integrated on the raster its description states* (`test_mrd_reader.py`). ::: - **Reproduces the spoiled gradient-echo steady state from a description read either way.** diff --git a/docs/explanations/epg.md b/docs/explanations/epg.md index 373feb98..8b133d20 100644 --- a/docs/explanations/epg.md +++ b/docs/explanations/epg.md @@ -32,10 +32,10 @@ kernels hold. The derivations and the operator-by-operator detail are in :::{dropdown} Show code and tests Code: `transmit_field`, {class}`~blochsim.ShimDefinition`. Tests: *channels in antiphase cancel*; *the field is the complex sum written out* (`test_transmit_array.py`). ::: -- **Integrates a shaped pulse into the rotation it performs across the slice, or per voxel when its channel weights vary while it plays.** +- **Integrates a shaped pulse into the rotation it performs across the slice, or per voxel when its channel weights vary while it plays; a pulse selecting more than `NONSELECTIVE_THICKNESS_M` (20 mm) is taken at its centre rotation, a slab the encoding resolves.** :::{dropdown} Show code and tests - Code: `exact_slice_profile`, `transition_table`, `dynamic_pair`. Tests: *the slice shows up when the pulse selects one* (`test_exact_profile.py`); *the pair is what exponentiating each sample gives* (`test_dynamic_transmit.py`). + Code: `exact_slice_profile`, `transition_table`, `dynamic_pair`. Tests: *the slice shows up when the pulse selects one*; *a slab thicker than the threshold is profiled as its centre* (`test_exact_profile.py`); *the pair is what exponentiating each sample gives* (`test_dynamic_transmit.py`). ::: - **Plays a waveform sample by sample, with relaxation, precession and exchange between samples.** diff --git a/src/blochsim/sequence/__init__.py b/src/blochsim/sequence/__init__.py index f0c4865b..c299234c 100644 --- a/src/blochsim/sequence/__init__.py +++ b/src/blochsim/sequence/__init__.py @@ -5,6 +5,7 @@ __all__ = [ "SequenceDescription", "ShimDefinition", + "description_from_waveforms", "read_mrd_description", "FSEReadout", "Operator", @@ -53,7 +54,7 @@ ideal_rf_definition, # noqa: F401 rf_definition, ) -from ._mrd import read_mrd_description +from ._mrd import description_from_waveforms, read_mrd_description from ._operators import ( Delay, Dephase, diff --git a/src/blochsim/sequence/_accelerators.py b/src/blochsim/sequence/_accelerators.py index e1fc8b4a..17f9933f 100644 --- a/src/blochsim/sequence/_accelerators.py +++ b/src/blochsim/sequence/_accelerators.py @@ -5,6 +5,7 @@ __all__: list[str] = [] import functools +import math import os import weakref from collections import OrderedDict @@ -288,6 +289,40 @@ def _pack_for( ) +#: Thickness in m above which a pulse is taken as non-selective: a slab this +#: thick is resolved by the encoding, so each position of it sees the pulse's +#: centre rotation, not a profile across it. +NONSELECTIVE_THICKNESS_M = 20e-3 + + +def _selected_thickness_m(description: SequenceDescription) -> dict[int, float]: + """The thinnest slice each RF definition is played over, in m; 0 where unknown. + + A definition's ``bandwidth_hz`` over the event's slice-select gradient. A + description that states no gradient or no bandwidth states no thickness, + and its pulse keeps the mode its waveform asks for. + """ + thickness: dict[int, float] = {} + for event in description.events: + if event.type is not EventType.RF: + continue + key = event.rf_definition_id + bandwidth = float(description.rf_definitions[key].bandwidth_hz) + gradient = abs(float(event.slice_select_gradient_hz_per_m)) + here = bandwidth / gradient if bandwidth > 0.0 and gradient > 0.0 else 0.0 + thickness[key] = min(thickness.get(key, math.inf), here) + return thickness + + +def _selective(description: SequenceDescription, key: int) -> bool: + """Whether a definition selects a slice thin enough to profile across, or states none.""" + cached = getattr(description, "_thickness_cache", None) + if cached is None: + cached = _selected_thickness_m(description) + object.__setattr__(description, "_thickness_cache", cached) + return cached.get(key, math.inf) <= NONSELECTIVE_THICKNESS_M + + def _wants_a_table(description: SequenceDescription, slice_profile: Any) -> bool: """Whether the pulses are to be integrated into a rotation over the slice. @@ -311,6 +346,7 @@ def _wants_a_table(description: SequenceDescription, slice_profile: Any) -> bool } cached = any( description.rf_definitions[key].rf_mode() is RfMode.PROFILED + and _selective(description, key) for key in driven ) object.__setattr__(description, "_profiled_cache", cached) @@ -327,7 +363,9 @@ def _pulse_shapes( different shapes get one each, and the event carries which it reads. Ordered by first use, so the row a definition takes does not depend on how - the description happened to number it. + the description happened to number it. A pulse selecting more than + :data:`NONSELECTIVE_THICKNESS_M` is tabulated with no bandwidth, so every + position reads its centre rotation. Raises ------ @@ -339,7 +377,10 @@ def _pulse_shapes( continue key = event.rf_definition_id if key not in used: - used[key] = description.rf_definitions[key] + definition = description.rf_definitions[key] + if not _selective(description, key): + definition = replace(definition, bandwidth_hz=0.0) + used[key] = definition if not used: raise ValueError( "a slice profile is the rotation a pulse performs across the " diff --git a/src/blochsim/sequence/_mrd.py b/src/blochsim/sequence/_mrd.py index fb4fed20..4fed70c1 100644 --- a/src/blochsim/sequence/_mrd.py +++ b/src/blochsim/sequence/_mrd.py @@ -2,19 +2,19 @@ Pulserver's reconstruction proxy infers the description from the sequence file a series was played from and sends it after the MRD XML header and before the -first acquisition, as one ``TEXT`` message: a JSON object under the key -``pulserver_sequence_description``, holding one entry per file of the sequence -chain. What is read here is mechanical; which signal model the rows drive is -chosen by whoever reads them. +first acquisition, as ISMRMRD waveforms that every MRD writer keeps with the +series: per file of the sequence chain a header (``waveform_id`` 1024) and an +event list (1025), and one waveform per RF definition (1026), each one vector +of 32-bit words led by its length. What is read here is mechanical; which +signal model the rows drive is chosen by whoever reads them. """ from __future__ import annotations -__all__ = ["read_mrd_description"] +__all__ = ["description_from_waveforms", "read_mrd_description"] -import base64 -import json import struct +from collections.abc import Iterable from typing import Any, BinaryIO import numpy as np @@ -30,12 +30,19 @@ _trapezoid, ) -#: The key the description's JSON object opens with. -MESSAGE_KEY = "pulserver_sequence_description" - # MRD streaming message identifiers, and the fixed size of a FILENAME message. _FILENAME, _CONFIG, _HEADER, _CLOSE, _TEXT = 1, 2, 3, 4, 5 +_WAVEFORM = 1026 _FILENAME_BYTES = 1024 +# An ISMRMRD waveform header: its size, and where it holds the sample count, +# the channel count and the waveform id. +_WAVEFORM_HEADER_BYTES = 40 +_SAMPLES_AT, _CHANNELS_AT, _ID_AT = 28, 30, 36 + +#: Waveform ids of a description: header, event list, RF definition. +HEADER_ID, EVENTS_ID, RF_DEFINITION_ID = 1024, 1025, 1026 +#: Words an event takes in the event list. +_EVENT_WORDS = 10 def read_mrd_description( @@ -44,8 +51,9 @@ def read_mrd_description( """Read an MRD stream up to its sequence description. The messages ahead of the description -- a config, a config file name, the - XML header, other texts -- are read and the stream is left at the message - after the description, which is where its first acquisition is. + XML header, texts, other waveforms such as the physiological traces -- are + read and the stream is left at the message after the description's last + waveform, which is where its first acquisition is. Parameters ---------- @@ -70,6 +78,7 @@ def read_mrd_description( description. """ header = None + collected: list[tuple[int, np.ndarray]] = [] while True: (identifier,) = struct.unpack(" bytes: - read = stream.read(count) - if len(read) != count: - raise ValueError("the MRD stream ended before the sequence description") - return read +def description_from_waveforms(waveforms: Iterable[Any]) -> list[SequenceDescription]: + """Read the descriptions a series' waveforms carry. + + Parameters + ---------- + waveforms : iterable + The series' waveforms: ``ismrmrd.Waveform`` objects, as an HDF5 file + read with ``ismrmrd.Dataset.read_waveform`` gives them, or + ``(waveform_id, data)`` pairs. Those of other ids, such as the + physiological traces, are passed over. + + Returns + ------- + list of SequenceDescription + One per file of the sequence chain, in play order; empty when the + waveforms carry none. + + Raises + ------ + ValueError + If a description's header names more RF definitions or events than + its waveforms carry. + """ + headers: dict[int, np.ndarray] = {} + events: dict[int, np.ndarray] = {} + definitions: dict[int, list[np.ndarray]] = {} + for item in waveforms: + identifier, data = ( + item if isinstance(item, tuple) else (item.waveform_id, item.data) + ) + words = _words(data) + if int(identifier) == HEADER_ID: + headers[int(words[1])] = words + elif int(identifier) == EVENTS_ID: + events[int(words[0])] = words + elif int(identifier) == RF_DEFINITION_ID: + definitions.setdefault(int(words[0]), []).append(words) + return [ + _from_words(headers[index], events.get(index), definitions.get(index, [])) + for index in sorted(headers) + ] -def _array(packed: str, dtype: str) -> np.ndarray: - return np.frombuffer(base64.b64decode(packed), dtype=dtype) +def _words(data: Any) -> np.ndarray: + """The word vector of a description waveform, its leading length dropped.""" + flat = np.ascontiguousarray(np.asarray(data), dtype=" SequenceDescription: - kinds = _array(fields["type"], " bool: + """Whether every subsequence's header, events and RF definitions have arrived.""" + headers = [_words(data) for kind, data in collected if kind == HEADER_ID] + if not headers or len(headers) < int(headers[0][2]): + return False + events = sum(kind == EVENTS_ID for kind, _ in collected) + definitions = sum(kind == RF_DEFINITION_ID for kind, _ in collected) + return events == len(headers) and definitions == sum(int(h[4]) for h in headers) + + +def _from_words( + header: np.ndarray, events: np.ndarray | None, definitions: list[np.ndarray] +) -> SequenceDescription: + index = int(header[1]) + if events is None or len(definitions) != int(header[4]): + raise ValueError( + f"subsequence {index}: the header names {int(header[4])} RF definitions " + f"and its waveforms carry {len(definitions)}, with " + f"{'no' if events is None else 'an'} event list" + ) + tr_duration_us, raster_s = (float(v) for v in header[5:9].view(" RfDefinition: + floats = words.view(" bytes: + read = stream.read(count) + if len(read) != count: + raise ValueError("the MRD stream ended before the sequence description") + return read + + def _event(kind: EventType, time_us: float, row: np.ndarray) -> SequenceEvent: if kind is EventType.RF: return SequenceEvent.rf( @@ -142,15 +250,18 @@ def _event(kind: EventType, time_us: float, row: np.ndarray) -> SequenceEvent: return SequenceEvent.wait(time_us) -def _shape(fields: dict[str, Any] | None) -> RfShape | None: - if fields is None: - return None - return RfShape(int(fields["num_uncompressed"]), _array(fields["samples"], " RfDefinition: - magnitude = _shape(fields["magnitude"]) - time = _shape(fields["time"]) +def _built_definition( + *, + identifier: int, + bandwidth_hz: float, + num_bands: int, + band_offsets: tuple[float, ...], + band_bandwidth_hz: float, + magnitude: RfShape, + phase: RfShape | None, + time: RfShape | None, + raster_s: float, +) -> RfDefinition: envelope = magnitude.decompress().astype(np.float64) time_s = ( (np.arange(envelope.size) + 0.5) * raster_s @@ -158,16 +269,14 @@ def _definition(fields: dict[str, Any], raster_s: float) -> RfDefinition: else time.decompress(scale=raster_s).astype(np.float64) ) return RfDefinition( - id=int(fields["id"]), - bandwidth_hz=float(fields["bandwidth_hz"]), - num_bands=int(fields["num_bands"]), - band_frequency_offsets_hz=tuple( - float(value) for value in fields["band_frequency_offsets_hz"] - ), - band_bandwidth_hz=float(fields["band_bandwidth_hz"]), + id=identifier, + bandwidth_hz=bandwidth_hz, + num_bands=num_bands, + band_frequency_offsets_hz=band_offsets, + band_bandwidth_hz=band_bandwidth_hz, # Not carried by the stream: the stored envelope's squared integral. total_b1sq_power=float(_trapezoid(envelope**2, time_s)), magnitude=magnitude, - phase=_shape(fields["phase"]), + phase=phase, time=time, ) diff --git a/src/blochsim/sequence/_simulation.py b/src/blochsim/sequence/_simulation.py index 2d35a106..6fb96bce 100644 --- a/src/blochsim/sequence/_simulation.py +++ b/src/blochsim/sequence/_simulation.py @@ -370,7 +370,7 @@ def simulate( record: RecordMode = "all", nstates: int | None = None, across_slice: ExactSliceProfile | None = None, - rf_raster_time_s: float = 1e-6, + rf_raster_time_s: float | None = None, device: torch.device | str | None = None, events: Any = None, _every: bool = False, @@ -388,11 +388,17 @@ def simulate( a caller that holds the structure fixed and rebuilds only what varies; see :mod:`blochsim.optim`. It must be the packing of ``description``. + ``rf_raster_time_s`` is the raster the RF sample times count in, the + description's own ``rf_raster_time_s`` by default: a pulse integrated + on another raster turns through another angle. + Raises ------ RuntimeError: if no fused kernel can take this tissue -- a device that has none, or a build without the extension. """ + if rf_raster_time_s is None: + rf_raster_time_s = description.rf_raster_time_s if repetitions == AUTO: return self._settled( description, @@ -498,7 +504,7 @@ def simulate_subspace( record: RecordMode = "all", nstates: int | None = None, across_slice: ExactSliceProfile | None = None, - rf_raster_time_s: float = 1e-6, + rf_raster_time_s: float | None = None, device: torch.device | str | None = None, ) -> Subspace: """Simulate a dictionary and return its leading temporal basis.""" diff --git a/tests/sequence/test_exact_profile.py b/tests/sequence/test_exact_profile.py index 1bad85f0..b408595d 100644 --- a/tests/sequence/test_exact_profile.py +++ b/tests/sequence/test_exact_profile.py @@ -394,3 +394,39 @@ def record(*args, **kwargs): assert largest > 0.0 for reference, result in zip(expected, fast, strict=True): assert float((reference - result).abs().max()) / largest < 1e-5 + + +def _under_gradient(description, gradient_hz_per_m): + """The description with every pulse played under ``gradient_hz_per_m``.""" + from dataclasses import replace + + events = tuple( + replace(event, params=tuple(event.params[:6]) + (gradient_hz_per_m,)) + if event.type.name == "RF" + else event + for event in description.events + ) + return replace(description, events=events) + + +def test_a_slab_thicker_than_the_threshold_is_profiled_as_its_centre() -> None: + """A slab past NONSELECTIVE_THICKNESS_M is resolved by the encoding, not averaged over.""" + from blochsim.sequence._accelerators import NONSELECTIVE_THICKNESS_M + + bandwidth = 4.0e3 + slab = _under_gradient( + _describe(definition=_sinc(bandwidth)), + bandwidth / (5 * NONSELECTIVE_THICKNESS_M), + ) + slice_ = _under_gradient( + _describe(definition=_sinc(bandwidth)), + bandwidth / (0.25 * NONSELECTIVE_THICKNESS_M), + ) + centre = _signal(slab) + + assert (_signal(slab, exact_slice_profile(15)) - centre).abs().max() < 1e-5 * ( + centre.abs().max() + ) + assert (_signal(slice_, exact_slice_profile(15)) - centre).abs().max() > 0.05 * ( + centre.abs().max() + ) diff --git a/tests/sequence/test_mrd_reader.py b/tests/sequence/test_mrd_reader.py index 0ea9ac5a..d9b1a1d9 100644 --- a/tests/sequence/test_mrd_reader.py +++ b/tests/sequence/test_mrd_reader.py @@ -1,7 +1,7 @@ """Reading the sequence description an MRD stream carries, as Pulserver's proxy sends it. -The stream is written here the way the proxy writes it -- the XML header, one -``TEXT`` message, then the first acquisition -- from a spoiled gradient echo +The stream is written here the way the proxy writes it -- the XML header, the +description's waveforms, then the first acquisition -- from a spoiled gradient echo built with pypulseq, and the description read back is checked against the rows it was written from and against the spoiled gradient-echo steady state, which has never seen a stream. @@ -11,7 +11,6 @@ import base64 import io -import json import math import struct @@ -21,7 +20,12 @@ import torch from blochsim import Excitation, Simulator, SPGRReadout -from blochsim.sequence import EventType, RfUse, read_mrd_description +from blochsim.sequence import ( + EventType, + RfUse, + description_from_waveforms, + read_mrd_description, +) from blochsim.simulators import SPGRSimulator FLIP_DEG = 30.0 @@ -117,11 +121,15 @@ def _message(identifier, text): return struct.pack(" 1e-3