"""Backend-independent description of an MR sequence state machine."""
from __future__ import annotations
__all__ = [
"AdcRole",
"EventAction",
"EventType",
"RfDefinition",
"RfMode",
"RfShape",
"RfUse",
"SequenceDescription",
"SequenceEvent",
"ShimDefinition",
"decompress_shape",
"ideal_rf_definition",
]
from collections.abc import Mapping
from dataclasses import dataclass, field
from enum import IntEnum, IntFlag
from typing import Any
import numpy as np
# ``np.trapz`` was renamed to ``np.trapezoid`` in NumPy 2.0.
_trapezoid = getattr(np, "trapezoid", None) or np.trapz
class EventType(IntEnum):
"""State-machine event type used by Pulserver SEQDESC streams."""
WAIT = 0
RF = 1
ADC = 2
class EventAction(IntFlag):
"""What the sequence does around an event, beyond the event itself.
A Pulseq file has no gradient "use" field, so a crusher cannot be told from
a phase encode by looking at waveforms -- that would mean tracking the
k-space moment through the whole TR. The sequence says instead what each
event's role implies, and the description carries it, which is what lets an
autogenerated stream reach the kernels without a policy being named.
``CRUSH_BEFORE`` and ``CRUSH_AFTER`` are the unbalanced gradients a
refocusing pulse sits between. ``SHIFT_AFTER`` is one played after a
readout, and ``SPOIL_AFTER`` is ideal transverse spoiling there instead.
The values are the bits the packed action word carries, so a description
and a kernel never hold two numberings of the same fact.
"""
NONE = 0
CRUSH_BEFORE = 1
CRUSH_AFTER = 2
SPOIL_AFTER = 8
SHIFT_AFTER = 16
class RfMode(IntEnum):
"""How a pulse's rotation is reached.
The three are alternatives rather than combinable. ``INSTANT`` turns
through a flip angle and a phase. ``PROFILED`` reads the rotation from a
table over slice position and effective flip. ``DYNAMIC`` reads it per
voxel, which is what a pulse whose channel weights vary while it plays
needs -- and which subsumes a profile, since a rotation integrated at the
voxel's own position is one.
"""
INSTANT = 0
PROFILED = 1
DYNAMIC = 2
class RfUse(IntEnum):
"""Pulseq RF-use tag.
The whole tag set, not the subset the operators that ship here happen to
emit. A file written elsewhere may carry any of them, and an operator
written later may emit one nothing emits today.
"""
UNKNOWN = 0
EXCITATION = 1
REFOCUSING = 2
INVERSION = 3
SATURATION = 4
PREPARATION = 5
OTHER = 6
class AdcRole(IntEnum):
"""Relationship between ADC events in a canonical event stream."""
NON_ACQUIRED = 0
SINGLE = 1
ECHO_CENTER = 2
NON_CENTER = 3
@dataclass(frozen=True)
class SequenceEvent:
"""One WAIT, RF, or ADC state-machine event.
``params`` deliberately accepts tensors. This lets a hard-coded sequence
description carry differentiable flip angles, phases, and timings without
introducing a second representation for sequence optimization.
Attributes
----------
type : EventType
Whether this is a wait, an RF pulse or an ADC.
timestamp_us : float or torch.Tensor
When it happens, in microseconds from the start of the stream.
params : tuple
The payload its type defines, read through the named properties rather
than by position.
action : EventAction
What the sequence does around it -- the crushers a refocusing pulse
sits between, the gradient or the ideal spoiling after a readout.
"""
type: EventType
timestamp_us: Any
params: tuple[Any, ...] = ()
action: EventAction = EventAction.NONE
@classmethod
def wait(cls, timestamp_us: Any) -> SequenceEvent:
"""Create a timing-only event."""
return cls(EventType.WAIT, timestamp_us)
@classmethod
def rf(
cls,
timestamp_us: Any,
definition_id: int,
use: RfUse,
amplitude_hz: Any,
phase_rad: Any = 0.0,
frequency_hz: Any = 0.0,
shim_id: int = 0,
slice_select_gradient_hz_per_m: Any = 0.0,
action: EventAction = EventAction.NONE,
) -> SequenceEvent:
"""Create an RF event using the Pulserver SEQDESC field order."""
return cls(
EventType.RF,
timestamp_us,
(
definition_id,
use,
amplitude_hz,
phase_rad,
frequency_hz,
shim_id,
slice_select_gradient_hz_per_m,
),
action,
)
@classmethod
def adc(
cls,
timestamp_us: Any,
role: AdcRole = AdcRole.SINGLE,
phase_rad: Any = 0.0,
*,
is_echo: bool = True,
action: EventAction = EventAction.NONE,
) -> SequenceEvent:
"""Create an ADC event using the Pulserver SEQDESC field order."""
return cls(EventType.ADC, timestamp_us, (role, phase_rad, int(is_echo)), action)
@property
def rf_definition_id(self) -> int:
self._require(EventType.RF)
return int(self.params[0])
@property
def rf_use(self) -> RfUse:
self._require(EventType.RF)
return RfUse(int(self.params[1]))
@property
def rf_amplitude_hz(self) -> Any:
self._require(EventType.RF)
return self.params[2]
@property
def rf_phase_rad(self) -> Any:
self._require(EventType.RF)
return self.params[3]
@property
def rf_frequency_hz(self) -> Any:
self._require(EventType.RF)
return self.params[4]
@property
def rf_shim_id(self) -> int:
self._require(EventType.RF)
return int(self.params[5])
@property
def slice_select_gradient_hz_per_m(self) -> Any:
self._require(EventType.RF)
return self.params[6]
@property
def adc_role(self) -> AdcRole:
self._require(EventType.ADC)
return AdcRole(int(self.params[0]))
@property
def adc_phase_rad(self) -> Any:
self._require(EventType.ADC)
return self.params[1]
@property
def is_echo(self) -> bool:
"""Whether this ADC position reaches the k-space origin."""
self._require(EventType.ADC)
return bool(int(self.params[2]))
def _require(self, expected: EventType) -> None:
if self.type is not expected:
raise AttributeError(
f"{expected.name} payload requested from {self.type.name} event"
)
@dataclass(frozen=True)
class RfShape:
"""One Pulseq derivative/run-length encoded shape.
Attributes
----------
num_uncompressed : int
How many samples the shape decodes to.
samples : array-like
The encoded samples, decoded by :meth:`decompress`.
"""
num_uncompressed: int
samples: Any = field(repr=False, compare=False)
def decompress(self, scale: float = 1.0) -> np.ndarray:
"""Return the decoded shape."""
return decompress_shape(self.samples, self.num_uncompressed, scale=scale)
@dataclass(frozen=True)
class RfDefinition:
"""RF definition referenced by RF events.
Attributes
----------
id : int
What an RF event names this pulse by.
bandwidth_hz : float
The pulse's bandwidth, which is the gradient the slice select is
quoted against.
num_bands : int
How many slices the pulse excites at once.
band_frequency_offsets_hz : tuple of float
Where each band sits, in Hz off the centre frequency.
band_bandwidth_hz : float
The bandwidth of one band.
total_b1sq_power : float
The pulse's integrated squared amplitude, which is what a power limit
is read against.
magnitude : RfShape or tuple of RfShape
The envelope, or one envelope per transmit channel.
phase : RfShape or tuple of RfShape, optional
The phase, matching ``magnitude``. ``None`` for a pulse of constant
phase.
time : RfShape, optional
Sample times, for a pulse on a non-uniform raster. ``None`` for a
uniform one.
gradient : RfShape, optional
The slice-select gradient per sample, in units ``bandwidth_hz`` is
quoted at. ``None`` where the gradient holds.
"""
id: int
bandwidth_hz: float
num_bands: int
band_frequency_offsets_hz: tuple[float, ...]
band_bandwidth_hz: float
total_b1sq_power: float
magnitude: RfShape | tuple[RfShape, ...]
phase: RfShape | tuple[RfShape, ...] | None = None
time: RfShape | None = None
gradient: RfShape | None = None
def __post_init__(self) -> None:
if isinstance(self.magnitude, tuple):
if not self.magnitude:
raise ValueError(f"RF definition {self.id} drives no channel")
if isinstance(self.phase, tuple):
if len(self.phase) != len(self.magnitude):
raise ValueError(
f"RF definition {self.id}: {len(self.magnitude)} "
f"magnitude channels against {len(self.phase)} phase"
)
elif self.phase is not None:
raise ValueError(
f"RF definition {self.id}: a per-channel magnitude takes a "
f"phase per channel, or none"
)
elif isinstance(self.phase, tuple):
raise ValueError(
f"RF definition {self.id}: a phase per channel takes a "
f"magnitude per channel"
)
@property
def channel_count(self) -> int:
"""How many transmit channels drive their own waveform."""
return len(self.magnitude) if isinstance(self.magnitude, tuple) else 1
def channel_envelopes(self) -> np.ndarray:
"""Return the normalized complex envelope per channel, ``(C, N)``."""
magnitudes = (
self.magnitude if isinstance(self.magnitude, tuple) else (self.magnitude,)
)
if isinstance(self.phase, tuple):
phases: tuple[RfShape | None, ...] = self.phase
else:
phases = (self.phase,) * len(magnitudes)
channels = []
for magnitude_shape, phase_shape in zip(magnitudes, phases, strict=False):
magnitude = magnitude_shape.decompress()
if phase_shape is None:
channels.append(magnitude.astype(np.complex64))
continue
phase = phase_shape.decompress(scale=2.0 * np.pi)
if phase.size != magnitude.size:
raise ValueError(
f"RF definition {self.id}: magnitude/phase size mismatch"
)
channels.append(magnitude * np.exp(1j * phase))
widths = {channel.size for channel in channels}
if len(widths) > 1:
raise ValueError(
f"RF definition {self.id}: channels of {sorted(widths)} samples"
)
return np.stack(channels)
def complex_envelope(self) -> np.ndarray:
"""Return the normalized complex RF envelope.
``(samples,)`` for a pulse every channel plays alike, ``(channels,
samples)`` for one whose channels each carry their own waveform.
"""
channels = self.channel_envelopes()
return channels[0] if self.channel_count == 1 else channels
def combined_envelope(self) -> np.ndarray:
"""Return the envelope a voxel of unit sensitivity on every channel sees.
The sum over channels, ``(samples,)``. This is the shape the flip
angle, the integral and the saturation are all read against, so each
keeps the meaning it has for a pulse driving a single channel.
"""
return self.channel_envelopes().sum(axis=0)
def gradient_waveform(self) -> np.ndarray | None:
"""Return the slice-select gradient per sample, or ``None`` if it holds.
In units of the gradient ``bandwidth_hz`` is quoted at, so a shape of
ones is the gradient the scanner does not move and reproduces the
constant case exactly. What a spin at a normalized position accrues
while sample ``s`` lasts is ``bandwidth_hz * position * g[s]`` in Hz.
Raises
------
ValueError: if the gradient and the envelope disagree on how many
samples the pulse has.
"""
if self.gradient is None:
return None
waveform = self.gradient.decompress()
samples = self.channel_envelopes().shape[1]
if waveform.size != samples:
raise ValueError(
f"RF definition {self.id}: the gradient carries "
f"{waveform.size} samples and the pulse {samples}"
)
return waveform
def sample_durations(self, *, rf_raster_time_s: float = 1e-6) -> np.ndarray:
"""Return how long each sample holds its field, in seconds.
``time`` names where along the pulse each sample sits, so what a sample
lasts is the distance between the midpoints either side of it, and the
whole spacing to its one neighbour at the two ends. A definition
declaring no time is played one sample per raster step.
This is what the rotation a pulse performs is composed from: the field
is held constant across a sample, so a longer sample turns further and
lets a spin off the slice centre accrue more phase.
Raises
------
ValueError: if the declared times and the envelope disagree on how
many samples the pulse has, or the times do not advance.
"""
samples = int(self.channel_envelopes().shape[1])
if self.time is None or samples < 2:
return np.full(samples, float(rf_raster_time_s), dtype=np.float64)
time_s = self.time.decompress(scale=rf_raster_time_s).astype(np.float64)
if time_s.size != samples:
raise ValueError(
f"RF definition {self.id}: {time_s.size} sample times against "
f"{samples} samples"
)
if not np.all(np.diff(time_s) > 0.0):
raise ValueError(
f"RF definition {self.id}: the sample times do not advance, so "
f"they do not say how long a sample lasts"
)
return np.gradient(time_s)
def rf_mode(self) -> RfMode:
"""Return which rotation this pulse's waveform asks the kernels for.
Memoized on the (immutable) definition: every event naming this pulse
asks the same question, and answering it reads the whole envelope.
"""
cached = getattr(self, "_mode_cache", None)
if cached is not None:
return cached
if self.channel_count > 1:
mode = RfMode.DYNAMIC
elif self._turns_about_one_axis():
mode = RfMode.INSTANT
else:
mode = RfMode.PROFILED
object.__setattr__(self, "_mode_cache", mode)
return mode
def _turns_about_one_axis(self) -> bool:
"""Whether a flip angle and a phase already say what this pulse does.
A rectangle played without slice selection holds its effective field
fixed for its whole duration, so integrating it lands on the rotation a
hard pulse performs -- and lands there exactly, where reading a
tabulated one at a flip between knots would not.
"""
if self.bandwidth_hz != 0.0:
return False
envelope = self.combined_envelope()
if envelope.size < 2:
return True
return bool(np.all(envelope == envelope[0]))
def integral(
self,
*,
rf_raster_time_s: float = 1e-6,
) -> complex:
"""Return the complex time integral of the normalized envelope.
The envelope is normalized, so this depends only on the pulse shape and
the raster, never on the amplitude of a particular occurrence. A
sequence typically reuses one definition for every refocusing pulse and
re-evaluates it on every simulation, so the result is memoized on the
(immutable) definition.
"""
cache = getattr(self, "_integral_cache", None)
if cache is None:
cache = {}
object.__setattr__(self, "_integral_cache", cache)
elif rf_raster_time_s in cache:
return cache[rf_raster_time_s]
envelope = self.combined_envelope()
if envelope.size < 2:
cache[rf_raster_time_s] = 0.0j
return 0.0j
if self.time is None:
time_s = (
np.arange(envelope.size, dtype=np.float64) + 0.5
) * rf_raster_time_s
else:
time_s = self.time.decompress(scale=rf_raster_time_s).astype(np.float64)
if time_s.size != envelope.size:
raise ValueError(
f"RF definition {self.id}: envelope/time size mismatch"
)
area = complex(_trapezoid(envelope, time_s))
cache[rf_raster_time_s] = area
return area
def saturation(self, *, rf_raster_time_s: float = 1e-6) -> float:
"""Return what this shape deposits in a bound pool, per flip squared.
The saturation an RF pulse applies to a pool with no transverse
magnetization is ``-pi gamma**2 int B1**2 dt * G(df)``. Writing the
pulse as a normalized envelope ``w`` driven to a flip angle,
``B1 = w / (gamma int w dt)`` per radian, so::
-pi gamma**2 int B1**2 dt == -pi int |w|**2 dt / |int w dt|**2
and the gyromagnetic ratio cancels: what is left depends on the shape
alone. Multiplied by the square of the flip a voxel actually takes and
by the lineshape at the offset it is played at, this is the exponent of
the bound pool's saturation.
In 1/s. A pulse with no area saturates nothing that a flip angle can
describe, so it reports zero.
"""
cache = getattr(self, "_saturation_cache", None)
if cache is None:
cache = {}
object.__setattr__(self, "_saturation_cache", cache)
elif rf_raster_time_s in cache:
return cache[rf_raster_time_s]
area = abs(self.integral(rf_raster_time_s=rf_raster_time_s))
envelope = self.combined_envelope()
if area == 0.0 or envelope.size < 2:
cache[rf_raster_time_s] = 0.0
return 0.0
if self.time is None:
time_s = (
np.arange(envelope.size, dtype=np.float64) + 0.5
) * rf_raster_time_s
else:
time_s = self.time.decompress(scale=rf_raster_time_s).astype(np.float64)
power = float(_trapezoid(np.abs(envelope) ** 2, time_s))
value = -np.pi * power / (area * area)
cache[rf_raster_time_s] = value
return value
def flip_angle(
self,
amplitude_hz: Any,
*,
rf_raster_time_s: float = 1e-6,
) -> tuple[Any, Any]:
"""Return ``(flip_rad, integral_phase_rad)`` for one RF occurrence.
Tensor amplitudes remain tensors, so both reverse- and forward-mode
automatic differentiation pass through a description unchanged.
"""
area = self.integral(rf_raster_time_s=rf_raster_time_s)
try:
import torch
except ImportError: # pragma: no cover - Torch is a required dependency
torch = None
if torch is not None and isinstance(amplitude_hz, torch.Tensor):
# The area is a Python number, so it multiplies a tensor without
# being made one first. A pulse whose envelope integrates to a real
# number -- every shape that is symmetric about its centre, which
# is most of them -- keeps the product real, and a real product is
# both cheaper per event and the only kind the real-subspace
# kernels can be handed.
signed_area = amplitude_hz * (area.real if area.imag == 0.0 else area)
return 2.0 * torch.pi * torch.abs(signed_area), torch.angle(signed_area)
signed_area = amplitude_hz * area
return 2.0 * np.pi * abs(signed_area), np.angle(signed_area)
[docs]
@dataclass(frozen=True)
class ShimDefinition:
"""One transmit-shim definition: a complex weight per channel.
The weights hold for the whole pulse, and every channel plays the same
waveform -- static parallel transmit. :mod:`._transmit` combines them with
the per-channel sensitivities into the field a voxel sees.
Attributes
----------
id : int
What an RF event names this shim by.
magnitudes : tuple of float
One weight per transmit channel.
phases_rad : tuple of float
One phase per transmit channel, in radians.
"""
id: int
magnitudes: tuple[float, ...]
phases_rad: tuple[float, ...]
[docs]
@dataclass(frozen=True)
class SequenceDescription:
"""Event stream and RF resources for one sequence or subsequence.
``crusher_dephasing_rad`` is the dephasing the sequence's unbalanced
gradient winds across one voxel of ``voxel_size_m``. Their ratio sets the
b-factor diffusion is damped by and the phase each dephasing order turns
through under flow, so a dephasing of zero leaves both terms out however
large the coefficient behind them.
``voxel_size_m`` also stands on its own: spins cross a voxel of that size
at their own speed whether or not a gradient is playing, which is what
washes the voxel out and replaces it with unexcited magnetization. Leaving
it ``None`` declares no voxel, and a spin velocity then has nothing to
cross.
Attributes
----------
subsequence_index : int
Which subsequence of a multi-part acquisition this is.
tr_duration_us : float or torch.Tensor
How long one playing of this stream holds the timeline, in
microseconds.
events : tuple of SequenceEvent
The stream itself, in the order it is played.
rf_definitions : dict
``{id: RfDefinition}`` -- the pulses the RF events name.
shim_definitions : dict
``{id: ShimDefinition}`` -- the transmit shims the RF events drive.
crusher_dephasing_rad : float
Dephasing wound across one voxel by the sequence's unbalanced
gradient.
voxel_size_m : float, optional
The voxel a spin crosses, and the length the dephasing is wound
across. ``None`` declares no voxel.
rf_raster_time_s : float
The dwell the RF shapes are sampled on. A pulse is integrated against
it, so a stream carrying shapes on a 2 microsecond raster and read on a
1 microsecond one turns half the angle it should.
"""
subsequence_index: int
tr_duration_us: Any
events: tuple[SequenceEvent, ...]
rf_definitions: dict[int, RfDefinition]
shim_definitions: dict[int, ShimDefinition] = field(default_factory=dict)
crusher_dephasing_rad: float = 0.0
voxel_size_m: float | None = None
rf_raster_time_s: float = 1e-6
[docs]
@classmethod
def from_operators(
cls,
*parts: Any,
rf_definitions: Mapping[int, RfDefinition] | None = None,
shim_definitions: Mapping[int, ShimDefinition] | None = None,
subsequence_index: int = 0,
crusher_dephasing_rad: float = 0.0,
voxel_size_m: float | None = None,
start_s: Any = 0.0,
) -> SequenceDescription:
"""Lay operators out end to end and return the stream they make.
The other way in besides a scanner: what a sequence written here plays,
in the same object one that arrived would be read into.
Parameters
----------
parts : Operator, optional
The operators, in the order they play.
rf_definitions : mapping, optional
The pulses the events drive, by id. One ideal pulse under id zero
when not given, which is what an operator naming no shape asks for.
shim_definitions : mapping, optional
The transmit shims the pulses are driven on, by id.
subsequence_index : int, optional
Which subsequence of a scan this is.
crusher_dephasing_rad : float, optional
The turn one crusher puts across a voxel. Diffusion and flow are
read off it and are zero without it.
voxel_size_m : float, optional
The distance that turn is put across.
start_s : float or torch.Tensor, optional
When the first operator starts.
Returns
-------
The stream, lasting as long as the operators laid out end to end.
"""
from ._operators import compose
events, duration_s = compose(*parts, start_s=start_s)
return cls(
subsequence_index=subsequence_index,
tr_duration_us=1e6 * duration_s,
events=events,
rf_definitions=dict(
{0: ideal_rf_definition()} if rf_definitions is None else rf_definitions
),
shim_definitions=dict(shim_definitions or {}),
crusher_dephasing_rad=crusher_dephasing_rad,
voxel_size_m=voxel_size_m,
)
[docs]
@classmethod
def from_pulseq(cls, source: Any, **settings: Any) -> SequenceDescription:
"""Return what one repetition of a Pulseq sequence plays.
Parameters
----------
source : str, Path or sequence
A ``.seq`` file, or a sequence in memory with pypulseq's reading
interface -- pypulseq's own ``Sequence``, or pypulseqpp's.
settings : Any
``tr_index`` to name which repetition to read, and the
``subsequence_index``, ``crusher_dephasing_rad`` and
``voxel_size_m`` the description carries.
Returns
-------
The stream, in the same object a scanner's own would arrive in.
"""
from ._pulseq import read_pulseq_description
return read_pulseq_description(source, **settings)
[docs]
def plot(self, axis: Any = None, *, upto_s: float | None = None) -> Any:
"""Draw the stream: pulses as stems, samples as markers.
What a layout actually laid down, which is easier to check here than
in the signal it produces. A pulse's stem is as tall as the angle it
turns.
Parameters
----------
axis : matplotlib.axes.Axes, optional
Where to draw. A new figure when not given.
upto_s : float, optional
Stop after this much of the stream, for a repetition long enough
that all of it is a smear.
Returns
-------
The axis drawn on.
Raises
------
ImportError: if matplotlib is not installed, which is not a
dependency of this package.
"""
try:
import matplotlib.pyplot as plt
except ImportError as reason: # pragma: no cover - depends on the env
raise ImportError(
"plotting a description needs matplotlib, which blochsim does "
"not depend on; pip install matplotlib"
) from reason
axis = plt.subplots()[1] if axis is None else axis
limit = float("inf") if upto_s is None else upto_s * 1e6
for event in self.events:
when = float(event.timestamp_us)
if when > limit:
break
if event.type is EventType.RF:
axis.vlines(
when / 1000.0,
0.0,
float(event.rf_amplitude_hz) * 180.0 / np.pi,
color="crimson",
lw=2.5,
)
elif event.type is EventType.ADC:
axis.plot(when / 1000.0, 0.0, "v", color="tab:blue", ms=9)
axis.plot([], [], color="crimson", lw=2.5, label="RF, height is the flip")
axis.plot([], [], "v", color="tab:blue", ms=9, label="ADC")
axis.set(xlabel="time [ms]", ylabel="flip angle [deg]")
axis.grid(alpha=0.3)
return axis
@property
def adc_events(self) -> tuple[SequenceEvent, ...]:
"""Return ADC events in state-machine order."""
return tuple(event for event in self.events if event.type is EventType.ADC)
def decompress_shape(
packed: Any,
num_uncompressed: int,
*,
scale: float = 1.0,
) -> np.ndarray:
"""Decode Pulseq derivative/run-length shape compression."""
packed = np.asarray(packed, dtype=np.float32).reshape(-1)
if num_uncompressed < 0:
raise ValueError("negative uncompressed shape length")
if packed.size == num_uncompressed:
return (packed * scale).astype(np.float32, copy=True)
if num_uncompressed == 0 and packed.size == 0:
return np.empty(0, dtype=np.float32)
delta = np.empty(num_uncompressed, dtype=np.float32)
source = 0
target = 0
while source < packed.size:
value = packed[source]
if source + 1 < packed.size and packed[source + 1] == value:
if source + 2 >= packed.size:
raise ValueError("truncated Pulseq RLE triplet")
count_value = float(packed[source + 2]) + 2.0
count = round(count_value)
if abs(count_value - count) > 1e-6 or count < 2:
raise ValueError("invalid Pulseq RLE repeat count")
if target + count > num_uncompressed:
raise ValueError("Pulseq RLE shape expands past declared length")
delta[target : target + count] = value
target += count
source += 3
else:
if target >= num_uncompressed:
raise ValueError("Pulseq shape expands past declared length")
delta[target] = value
target += 1
source += 1
if target != num_uncompressed:
raise ValueError(
f"Pulseq shape expanded to {target}, expected {num_uncompressed}"
)
return (np.cumsum(delta, dtype=np.float32) * scale).astype(np.float32)
[docs]
def rf_definition(
waveform: Any,
*,
dwell_s: float,
bandwidth_hz: float = 0.0,
definition_id: int = 0,
band_frequency_offsets_hz: Any = None,
band_bandwidth_hz: float | None = None,
rf_raster_time_s: float = 1e-6,
) -> RfDefinition:
"""Return the definition of a pulse from the envelope a scanner plays.
The envelope is whatever a designer or a Pulseq block hands over -- a
complex waveform on a uniform raster, in any amplitude units, one row per
transmit channel where there is more than one. What it is scaled to does
not matter and must not be relied on: it is rescaled here so that its
integral is ``1 / (2 pi)`` seconds, which is what makes an event's
amplitude, written in radians, the flip that event turns. A pulse handed
over peak-normalized and left that way turns a fraction of a degree and
returns a signal of nothing, which is the one mistake this exists to stop.
The pulse's own phase is kept rather than divided out, so a linear-phase
pulse keeps the axis it turns about.
Parameters
----------
waveform : array-like
The complex envelope, ``(samples,)`` or ``(channels, samples)``.
dwell_s : float
How long one sample is held, in seconds.
bandwidth_hz : float, optional
The pulse's bandwidth, which the slice select is quoted against. Zero
-- the default -- is a pulse that selects nothing and stays on the
instant operator, and a non-zero one is what makes a pulse slice
selective. There is no separate flag for that, here or in the stream a
description arrives on.
definition_id : int, optional
The identifier events name this pulse by.
band_frequency_offsets_hz : sequence of float, optional
Where each band sits, for a multiband pulse. One band at zero when not
given.
band_bandwidth_hz : float, optional
The bandwidth of one band, defaulting to ``bandwidth_hz``.
rf_raster_time_s : float, optional
The raster the sample times are written against, which is the one the
run will read them back with. Sample times are carried in units of it
rather than in seconds, so a pulse whose dwell is not the raster still
says how long it lasts.
Returns
-------
RfDefinition
The pulse, scaled so an event amplitude in radians is a flip in
radians.
Raises
------
ValueError
If the envelope is empty, if ``dwell_s`` is not positive, or if the
waveform integrates to nothing and so has no flip angle to be read
against.
"""
samples = np.asarray(waveform)
if samples.ndim == 1:
samples = samples[None, :]
if samples.ndim != 2 or samples.size == 0:
raise ValueError(
"an RF envelope is (samples,) or (channels, samples), and not empty"
)
if not dwell_s > 0.0:
raise ValueError(f"dwell_s must be positive, not {dwell_s!r}")
if not rf_raster_time_s > 0.0:
raise ValueError(f"rf_raster_time_s must be positive, not {rf_raster_time_s!r}")
count = samples.shape[-1]
# In units of the raster, which is how a description carries sample times
# and how a run reads them back.
times = RfShape(
count,
(np.arange(count, dtype=np.float64) * dwell_s / rf_raster_time_s).astype(
np.float32
),
)
def shapes(values: np.ndarray) -> Any:
rows = tuple(RfShape(count, row.astype(np.float32)) for row in values)
return rows[0] if len(rows) == 1 else rows
def built(values: np.ndarray, **extra: Any) -> RfDefinition:
return RfDefinition(
id=definition_id,
bandwidth_hz=float(bandwidth_hz),
num_bands=len(offsets),
band_frequency_offsets_hz=offsets,
band_bandwidth_hz=float(
bandwidth_hz if band_bandwidth_hz is None else band_bandwidth_hz
),
magnitude=shapes(np.abs(values)),
phase=shapes(np.angle(values) / (2.0 * np.pi)),
time=times,
**extra,
)
offsets = (
(0.0,)
if band_frequency_offsets_hz is None
else tuple(float(value) for value in band_frequency_offsets_hz)
)
if not 1 <= len(offsets) <= 8:
raise ValueError(
f"a pulse carries one to eight bands, not {len(offsets)}: that is "
f"how many a description has room for"
)
# Scaled against the integral the definition itself reads, rather than one
# taken here: the two would have to agree about the quadrature and about
# what the last sample's duration is, and only one of them decides.
area = abs(
built(samples, total_b1sq_power=0.0).integral(rf_raster_time_s=rf_raster_time_s)
)
if area <= 0.0:
raise ValueError(
"this envelope integrates to nothing, so it has no flip angle to "
"be read against"
)
scaled = samples / (2.0 * np.pi * area)
return built(scaled, total_b1sq_power=float((np.abs(scaled) ** 2).sum() * dwell_s))
def ideal_rf_definition(definition_id: int = 0) -> RfDefinition:
"""Return the definition of a pulse that turns instantly and selects nothing.
Its envelope integrates to exactly ``1 / (2 pi)`` seconds, so an event
amplitude expressed in radians comes back unchanged as the flip that event
turns. It declares no bandwidth, which is what keeps a sequence built from
it on the instant operator rather than on a tabulated rotation.
Parameters
----------
definition_id : int, optional
The identifier events name this pulse by.
Returns
-------
RfDefinition
The ideal hard pulse.
"""
duration_us = 1e6 / (2.0 * np.pi)
return RfDefinition(
id=definition_id,
bandwidth_hz=0.0,
num_bands=1,
band_frequency_offsets_hz=(0.0,),
band_bandwidth_hz=0.0,
total_b1sq_power=0.0,
magnitude=RfShape(2, np.ones(2, dtype=np.float32)),
time=RfShape(2, np.asarray([0.0, duration_us], dtype=np.float32)),
)
