diff --git a/.github/workflows/deploy_docs.yml b/.github/workflows/deploy_docs.yml index 7c2bfda..786d67d 100644 --- a/.github/workflows/deploy_docs.yml +++ b/.github/workflows/deploy_docs.yml @@ -15,7 +15,6 @@ jobs: - uses: actions/checkout@v4 with: fetch-depth: 0 - submodules: true - name: Install latest mdbook run: | tag=$(curl 'https://api.github.com/repos/rust-lang/mdbook/releases/latest' | jq -r '.tag_name') diff --git a/.gitmodules b/.gitmodules deleted file mode 100644 index 05f63bd..0000000 --- a/.gitmodules +++ /dev/null @@ -1,3 +0,0 @@ -[submodule "documentation/src/external/pulseq-zero"] - path = documentation/src/external/pulseq-zero - url = https://github.com/MRsources/pulseq-zero.git diff --git a/CHANGELOG.md b/CHANGELOG.md index 0440345..d9fd427 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,3 +1,9 @@ +- 1.2.0 + - add `fid_graph`: builds a graph with only `z0` and one `+` state per repetition, as an alternative to `compute_graph` without running the prepass + - fix: `import_pulseq_rs` stores `pulse.freq_offset` in Hz, matching the pydissect importer (was rad/s) (#197) + - speed up `import_pulseq_rs` by integrating gradient moments in NumPy instead of once per event, ~4-13x faster on typical sequences (#194) + - fix: cast values to native Python floats before writing the MRD header, so the XML no longer contains `np.float64(...)` (#189) + - fix: update pulseq-rs for full rotation extension support and soft-delay default - 1.1.2 - fix: diffusion for z states was missing a (2pi)² in the b-value computation - 1.1.1 diff --git a/Cargo.lock b/Cargo.lock index e029cfb..4597482 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -67,7 +67,7 @@ dependencies = [ [[package]] name = "mrzero_core" -version = "1.1.2" +version = "1.2.0" dependencies = [ "num-complex", "pulseq-rs", @@ -138,8 +138,8 @@ dependencies = [ [[package]] name = "pulseq-rs" -version = "0.2.2" -source = "git+https://github.com/pulseq-frame/pulseq-rs.git#b93aa3a492036862734f3d05ee0e2ef51268f409" +version = "0.3.0" +source = "git+https://github.com/pulseq-frame/pulseq-rs.git#9e763092357981892ecfec9438a9cbc51db07144" dependencies = [ "num-complex", "thiserror", diff --git a/Cargo.toml b/Cargo.toml index bf08510..c5e18ff 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -1,6 +1,6 @@ [package] name = "mrzero_core" -version = "1.1.2" +version = "1.2.0" edition = "2021" # See more keys and their definitions at https://doc.rust-lang.org/cargo/reference/manifest.html diff --git a/documentation/src/SUMMARY.md b/documentation/src/SUMMARY.md index 43c5361..718e8a4 100644 --- a/documentation/src/SUMMARY.md +++ b/documentation/src/SUMMARY.md @@ -11,7 +11,6 @@ - [Pulseq Integration](integration-pulseq.md) - [Matlab Integration](integration-matlab.md) - [Pulseq-zero](integration-pulseqzero.md) - - [ISMRMD]() - [FAQ](faq.md) --- @@ -22,7 +21,6 @@ - [Loading and generating phantoms](nifti-generate.md) - [Specification](nifti-spec.md) - [Coordinate System](nifti-coords.md) -- [Motion simulation]() - [State selection](state_selection.md) - [API overview](api.md) - [Sequence building blocks](api-sequence.md) @@ -35,12 +33,8 @@ - [MR-zero FFT/sign conventions](fft_sign_conventions.md) -# 📊 Pulseq-zero - -- [Introduction](pulseqzero-intro.md) -- [Math methods for differentiability](pulseqzero-math.md) - - --- [Literature](literature.md) + +[Changelog](changelog.md) diff --git a/documentation/src/api-simulation.md b/documentation/src/api-simulation.md index f6f645d..d956213 100644 --- a/documentation/src/api-simulation.md +++ b/documentation/src/api-simulation.md @@ -83,6 +83,23 @@ Returns the computed [`Graph`](#graph). | `size` | used to scale gradients to phantom FOV if the [`Sequence`](api-sequence.md#sequence) sets `normalize_grads=True` | +### `fid_graph` + +Alternative to [`compute_graph`](#compute_graph) that skips the pre-pass. +Returns a [`Graph`](#graph) that only contains `z0` and a single `+` state per repetition, which [`execute_graph`](#execute_graph) simulates as FID signal. +Every excitation pulse (`PulseUsage.EXCIT` or `UNDEF`) creates a new `+` state from `z0` and discards the previous one; refocusing pulses (`REFOC`) conjugate it, all other pulses only scale it. +Stimulated echoes and other states are never simulated, so this is only exact for perfectly spoiled sequences. + +```python +graph = mr0.fid_graph(seq) +signal = mr0.execute_graph(graph, seq, data) +``` + +| Parameter | Description | +| --------- | ----------- | +| `seq` | the [`Sequence`](api-sequence.md#sequence) to compute the graph for | + + ### `execute_graph` Takes the [`Sequence`](api-sequence.md#sequence), the [`SimData`](api-phantom.md#simdata) and the [`Graph`](#graph) and calculates the measured ADC signal. diff --git a/documentation/src/api.md b/documentation/src/api.md index ae9b6c4..8ba37aa 100644 --- a/documentation/src/api.md +++ b/documentation/src/api.md @@ -31,6 +31,7 @@ import MRzeroCore as mr0 - [`mr0.isochromat_sim`](api-simulation.md#isochromat_sim) - [`mr0.compute_graph`](api-simulation.md#compute_graph) - [`mr0.compute_graph_ext`](api-simulation.md#compute_graph_ext) + - [`mr0.fid_graph`](api-simulation.md#fid_graph) - [`mr0.Graph`](api-simulation.md#graph) - [`mr0.execute_graph`](api-simulation.md#execute_graph) - [Reconstruction](api-reconstruction.md) diff --git a/documentation/src/changelog.md b/documentation/src/changelog.md new file mode 100644 index 0000000..cd20421 --- /dev/null +++ b/documentation/src/changelog.md @@ -0,0 +1,3 @@ +# Changelog + +{{#include ../../CHANGELOG.md}} diff --git a/documentation/src/external/pulseq-zero b/documentation/src/external/pulseq-zero deleted file mode 160000 index b96f112..0000000 --- a/documentation/src/external/pulseq-zero +++ /dev/null @@ -1 +0,0 @@ -Subproject commit b96f112808aaa0fd4d59dcf795c2b7f83a1359d4 diff --git a/documentation/src/faq.md b/documentation/src/faq.md index fa2966a..be01000 100644 --- a/documentation/src/faq.md +++ b/documentation/src/faq.md @@ -30,7 +30,7 @@ seq = pp.Sequence() signal, ktraj_adc = mr0.util.simulate(seq) # seq is PyPulseq object ``` Note, this is will still generate a temporary seq file. -For avoiding seq files completely, check out [Pulseq-zero](https://github.com/pulseq-frame/pulseq-zero) +For avoiding seq files completely, check out [Pulseq-zero](https://github.com/MRsources/pulseq-zero) ## How do I use custom phantoms with B0 field modifications? diff --git a/documentation/src/integration-pulseqzero.md b/documentation/src/integration-pulseqzero.md index 632cf9c..06046f8 100644 --- a/documentation/src/integration-pulseqzero.md +++ b/documentation/src/integration-pulseqzero.md @@ -1 +1,8 @@ # Pulseq-zero + +[Pulseq-zero](https://github.com/MRsources/pulseq-zero) lets you define sequences with [PyPulseq](https://github.com/imr-framework/pypulseq) and optimize them with MR-zero, without writing `.seq` files in between. + +- [README](https://github.com/MRsources/pulseq-zero/blob/main/README.md): installation, usage and examples +- [MATH.md](https://github.com/MRsources/pulseq-zero/blob/main/MATH.md): math methods used for differentiability + +Examples are listed in the [Playground](playground.md). diff --git a/documentation/src/playground.md b/documentation/src/playground.md index d0490ce..d6c55cc 100644 --- a/documentation/src/playground.md +++ b/documentation/src/playground.md @@ -40,7 +40,7 @@ These .seq files could also be measured on any MRI scanner using a Pulseq interp ## MR-zero optimization Gradient descent optimizations using automatic differentiation by backpropagation. -Some notebooks use [pulseq-zero](https://github.com/pulseq-frame/pulseq-zero) for optimizable sequence definitions with PyPulseq. +Some notebooks use [pulseq-zero](https://github.com/MRsources/pulseq-zero) for optimizable sequence definitions with PyPulseq. | Notebook | ![](https://colab.research.google.com/assets/colab-badge.svg) | | -------- | ------------------------------------------------------------- | diff --git a/documentation/src/pulseqzero-intro.md b/documentation/src/pulseqzero-intro.md deleted file mode 100644 index 877b689..0000000 --- a/documentation/src/pulseqzero-intro.md +++ /dev/null @@ -1 +0,0 @@ -{{#include external/pulseq-zero/README.md}} diff --git a/documentation/src/pulseqzero-math.md b/documentation/src/pulseqzero-math.md deleted file mode 100644 index 22ae0e6..0000000 --- a/documentation/src/pulseqzero-math.md +++ /dev/null @@ -1 +0,0 @@ -{{#include external/pulseq-zero/MATH.md}} diff --git a/pyproject.toml b/pyproject.toml index 6ea6eda..5191343 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -4,7 +4,7 @@ build-backend = "maturin" [project] name = "mrzerocore" -version = "1.1.2" +version = "1.2.0" description = "Core functionality of MRzero" authors = [ {name = "Jonathan Endres", email = "jonathan.endres@uk-erlangen.de"}, diff --git a/python/MRzeroCore/__init__.py b/python/MRzeroCore/__init__.py index e7518c6..6a5c97b 100644 --- a/python/MRzeroCore/__init__.py +++ b/python/MRzeroCore/__init__.py @@ -24,6 +24,7 @@ from .simulation.isochromat_sim import isochromat_sim from .simulation.pre_pass import compute_graph, compute_graph_ext, merge_graphs, Graph from .simulation.main_pass import execute_graph +from .simulation.fid_graph import fid_graph from .simulation.sig_to_mrd import sig_to_mrd from .reconstruction import reco_adjoint from .pulseq.exporter import pulseq_write_cartesian diff --git a/python/MRzeroCore/sequence.py b/python/MRzeroCore/sequence.py index a1d9217..615bb31 100644 --- a/python/MRzeroCore/sequence.py +++ b/python/MRzeroCore/sequence.py @@ -12,6 +12,74 @@ from .pulseq.pulseq_loader import intermediate, PulseqFile, Adc, Spoiler +def _gradient_cumulative(g, times=None): + """Moment from the start of ``g``'s shape to a shape-relative time. + + The returned function matches ``g.integrate(0, s)`` for an array of times, + so an interval moment is a difference of two evaluations. Piecewise-linear + trapezoids are integrated in NumPy instead of one PyO3 call per event. + """ + duration = float(g.shape_duration) + scale = float(g.amp) + if times is None: + times = g.shape_times() + t = np.asarray(times, dtype=np.float64) + a = np.asarray(g.shape_amp(), dtype=np.float64) + zeros = lambda s: np.zeros(np.shape(s), dtype=np.float64) + if t.size == 0 or duration <= 0.0 or scale == 0.0: + return zeros + if t.size == 1: + slope = scale * float(a[0]) + + def cum_const(s, slope=slope, duration=duration): + s = np.clip(np.asarray(s, dtype=np.float64), 0.0, duration) + return slope * s + + return cum_const + + # Hold the endpoint amplitude on [0, first] and [last, duration], matching + # Shape::interpolate outside the stored samples. + if t[0] > 0.0: + t = np.concatenate(([0.0], t)) + a = np.concatenate(([a[0]], a)) + if t[-1] < duration: + t = np.concatenate((t, [duration])) + a = np.concatenate((a, [a[-1]])) + if t.size > 1 and np.any(np.diff(t) <= 0.0): + keep = np.empty(t.size, dtype=bool) + keep[0] = True + keep[1:] = np.diff(t) > 0.0 + t = t[keep] + a = a[keep] + if t.size < 2: + slope = scale * float(a[0]) if a.size else 0.0 + + def cum_flat(s, slope=slope, duration=duration): + s = np.clip(np.asarray(s, dtype=np.float64), 0.0, duration) + return slope * s + + return cum_flat + + prefix = np.empty(t.size, dtype=np.float64) + prefix[0] = 0.0 + np.cumsum(0.5 * (a[:-1] + a[1:]) * np.diff(t), out=prefix[1:]) + prefix *= scale + a = a * scale + + def cum(s, t=t, a=a, prefix=prefix, duration=duration): + s = np.clip(np.asarray(s, dtype=np.float64), 0.0, duration) + idx = np.clip(np.searchsorted(t, s, side="right") - 1, 0, t.size - 2) + t0 = t[idx] + t1 = t[idx + 1] + a0 = a[idx] + width = t1 - t0 + frac = np.divide(s - t0, width, out=np.zeros_like(s), where=width > 0.0) + v = a0 + (a[idx + 1] - a0) * frac + return prefix[idx] + 0.5 * (a0 + v) * (s - t0) + + return cum + + class PulseUsage(Enum): """:class:`Enum` of all pulse usages, needed for reconstruction. @@ -843,36 +911,50 @@ def _import_pulseq_rs(cls, file_name: str, blocks = list(interp.blocks) duration = interp.duration - # Sorted list of times - can use bisect to find blocks faster - block_starts = [b.start for b in blocks] - block_ends = [b.start + b.duration for b in blocks] - - # Per-axis sorted breakpoint times for fast lookups + # Per-axis sorted breakpoint times for fast lookups, plus a cumulative + # moment function so each repetition is integrated in one NumPy pass. flat_grad_times: dict[str, list[float]] = {"x": [], "y": [], "z": []} + grad_cum: dict[str, list] = {"x": [], "y": [], "z": []} + grad_origin: dict[str, list[float]] = {"x": [], "y": [], "z": []} + grad_end: dict[str, list[float]] = {"x": [], "y": [], "z": []} for b in blocks: for axis in ("x", "y", "z"): g = getattr(b, "g" + axis) if g is None: continue + shape_t = g.shape_times() t_off = b.start + g.delay - flat_grad_times[axis].extend(t_off + t for t in g.shape_times()) - - # Precomputed ADC lists for faster lookup + flat_grad_times[axis].extend(t_off + t for t in shape_t) + grad_cum[axis].append(_gradient_cumulative(g, shape_t)) + grad_origin[axis].append(t_off) + grad_end[axis].append(t_off + float(g.shape_duration)) + + # Precomputed ADC lists for faster lookup. Label values are stored as + # parallel columns so a repetition can slice them instead of reading + # one dict per sample. flat_adc_times: list[float] = [] flat_adc_phases: list[float] = [] - flat_adc_labels: list[dict[str, int]] = [] label_names: list[str] = [] + label_cols: dict[str, list[int]] = {} for b in blocks: if b.adc is None: continue lbl = b.adc.labels() - if not label_names and lbl is not None: - label_names = list(lbl.keys()) ts = b.adc.sample_times() ph = b.adc.sample_phases() flat_adc_times.extend(b.start + t for t in ts) flat_adc_phases.extend(ph) - flat_adc_labels.extend(lbl for _ in ts) + n_samples = len(ts) + if lbl is not None and not label_names: + label_names = list(lbl.keys()) + label_cols = {name: [] for name in label_names} + if label_names: + for name in label_names: + value = int(lbl[name]) if lbl is not None else 0 + label_cols[name].extend([value] * n_samples) + label_arr = { + name: np.asarray(col, dtype=np.int32) for name, col in label_cols.items() + } def events_axis_in(t0: float, t1: float, axis: str) -> list[float]: times = flat_grad_times[axis] @@ -880,30 +962,21 @@ def events_axis_in(t0: float, t1: float, axis: str) -> list[float]: hi = bisect.bisect_right(times, t1) return times[lo:hi] - def adcs_in(t0: float, t1: float) -> tuple[ - list[float], list[float], list[dict[str, int]] - ]: + def adcs_in(t0: float, t1: float) -> tuple[list[float], list[float], int, int]: lo = bisect.bisect_left(flat_adc_times, t0) hi = bisect.bisect_right(flat_adc_times, t1) - return flat_adc_times[lo:hi], flat_adc_phases[lo:hi], flat_adc_labels[lo:hi] - - def integrate_axis(axis: str, t0: float, t1: float) -> float: - if t1 <= t0: - return 0.0 - moment = 0.0 - # First block with end > t0, up to (excluding) first block with - # start >= t1 - can use bisect to narrow down relevant range - lo_idx = bisect.bisect_right(block_ends, t0) - hi_idx = bisect.bisect_left(block_starts, t1) - for j in range(lo_idx, hi_idx): - g = getattr(blocks[j], "g" + axis) - if g is None: - continue - lo = max(t0, block_starts[j]) - block_starts[j] - hi = min(t1, block_ends[j]) - block_starts[j] - if hi > lo: - moment += g.integrate(lo, hi) - return moment + return flat_adc_times[lo:hi], flat_adc_phases[lo:hi], lo, hi + + def moments(axis: str, edges: np.ndarray) -> np.ndarray: + total = np.zeros(edges.shape[0], dtype=np.float64) + origins = grad_origin[axis] + ends = grad_end[axis] + lo = bisect.bisect_right(ends, float(edges[0])) + hi = bisect.bisect_left(origins, float(edges[-1])) + cums = grad_cum[axis] + for j in range(lo, hi): + total += cums[j](edges - origins[j]) + return np.diff(total) seq = cls(normalized_grads=False) @@ -972,7 +1045,7 @@ def integrate_axis(axis: str, t0: float, t1: float) -> float: else: shim_arr = torch.as_tensor(shims, dtype=torch.float32) - adcs, adc_phases, adc_label_snapshots = adcs_in(rep_start, rep_end) + adcs, adc_phases, adc_lo, adc_hi = adcs_in(rep_start, rep_end) if exact_trajectories: first = pe @@ -1020,22 +1093,17 @@ def integrate_axis(axis: str, t0: float, t1: float) -> float: f"{len(adcs)} are ADC (starting at {adc_start})" ) - mom_x = np.empty(event_count, dtype=np.float64) - mom_y = np.empty(event_count, dtype=np.float64) - mom_z = np.empty(event_count, dtype=np.float64) - for k in range(event_count): - t0 = abs_times[k] - t1 = abs_times[k + 1] - mom_x[k] = integrate_axis("x", t0, t1) - mom_y[k] = integrate_axis("y", t0, t1) - mom_z[k] = integrate_axis("z", t0, t1) + edges = np.asarray(abs_times, dtype=np.float64) + mom_x = moments("x", edges) + mom_y = moments("y", edges) + mom_z = moments("z", edges) rep = seq.new_rep(event_count) rep.pulse.angle = torch.as_tensor(angle) rep.pulse.phase = torch.as_tensor(phase) rep.pulse.usage = rf_usage rep.pulse.shim_array = shim_arr - rep.pulse.freq_offset = 2*torch.pi * rf.freq + rep.pulse.freq_offset = rf.freq rep.pulse.duration = rf.shape_duration rep.pulse.grad = rf_grad if rep.pulse.freq_offset != 0 or rf_hasgrad: @@ -1051,11 +1119,12 @@ def integrate_axis(axis: str, t0: float, t1: float) -> float: rep.adc_usage[adc_start:adc_start + len(adcs)] = 1 rep.adc_phase[adc_start:adc_start + len(adcs)] = phases_t - if label_names and adc_label_snapshots and adc_start is not None: + if label_names and adc_hi > adc_lo and adc_start is not None: for name in label_names: t = torch.zeros(event_count, dtype=torch.int32) - for k, snap in enumerate(adc_label_snapshots): - t[adc_start + k] = snap[name] + t[adc_start:adc_start + (adc_hi - adc_lo)] = torch.from_numpy( + label_arr[name][adc_lo:adc_hi].copy() + ) rep.adc_labels[name] = t if print_stats: diff --git a/python/MRzeroCore/simulation/fid_graph.py b/python/MRzeroCore/simulation/fid_graph.py new file mode 100644 index 0000000..0a0ce44 --- /dev/null +++ b/python/MRzeroCore/simulation/fid_graph.py @@ -0,0 +1,85 @@ +from __future__ import annotations +from typing import Optional + +from ..sequence import Sequence, PulseUsage +from .pre_pass import Graph + + +class FidDistribution: + """Python stand-in for the distributions returned by the prepass. + + Provides all attributes read by :func:`execute_graph` and + :meth:`Graph.plot`. There is no prepass simulation, so the prepass metrics + are placeholders: every state is marked as fully relevant for the signal. + """ + def __init__(self, dist_type: str, + ancestor: Optional[tuple[str, FidDistribution]] = None + ) -> None: + # Same layout as the prepass edges: (relation, ancestor, rot_factor). + # The main pass recomputes the rotation, so the factor is unused. + self.ancestors = [(*ancestor, 1.0)] if ancestor else [] + self.dist_type = dist_type + # Set by the main pass + self.mag = None + self.kt_vec = None + # Never skip a state because of the signal thresholds + self.emitted_signal = 1.0 if dist_type == '+' else 0.0 + self.latent_signal = 1.0 + self.latent_signal_unormalized = 1.0 + # Fake prepass data, only used for plotting + self.prepass_mag = 1.0 + self.signal = self.emitted_signal + self.prepass_kt_vec = [0.0, 0.0, 0.0, 0.0] + + def __repr__(self) -> str: + return ( + f"FidDist(type: {self.dist_type}, " + f"#ancestors: {len(self.ancestors)})" + ) + + +def fid_graph(seq: Sequence) -> Graph: + """Create a graph that only simulates z0 and one transverse state. + + Alternative to :func:`compute_graph` without running the prepass. The + returned graph can be executed by :func:`execute_graph`. Echoes from + stimulated or refocused states (except for the single tracked one) are + not simulated, so this is only exact for perfectly spoiled sequences. + + The ``+`` state is created by every excitation (or undefined) pulse, + which discards previous transverse magnetisation. Refocusing pulses + conjugate it, all other pulses only scale it. + + Parameters + ---------- + seq : Sequence + The sequence that produces the returned graph + + Returns + ------- + Graph + Graph with one ``z0`` and at most one ``+`` state per repetition + """ + # The first repetition contains only one element: A fully relaxed z0 + graph = [[FidDistribution('z0')]] + + for rep in seq: + pulse = rep.pulse + prev = graph[-1] + dists = [FidDistribution('z0', ('zz', prev[0]))] + + if pulse.usage in (PulseUsage.EXCIT, PulseUsage.UNDEF): + dists.append(FidDistribution('+', ('z+', prev[0]))) + elif len(prev) == 2: + # Can't refocus if nothing is excited + if pulse.usage == PulseUsage.REFOC: + dists.append(FidDistribution('+', ('-+', prev[1]))) + else: + dists.append(FidDistribution('+', ('++', prev[1]))) + + # +z and -z are never used because there is no pulse.usage that flips + # excited magnetisation back to z + + graph.append(dists) + + return Graph(graph) diff --git a/python/MRzeroCore/simulation/sig_to_mrd.py b/python/MRzeroCore/simulation/sig_to_mrd.py index 45aec2b..6343ca9 100644 --- a/python/MRzeroCore/simulation/sig_to_mrd.py +++ b/python/MRzeroCore/simulation/sig_to_mrd.py @@ -245,10 +245,11 @@ def _seq_write_mrd_head( """Writes the pulseq sequence definitions and labels to the ISMRMRD header.""" def m_to_mm(x_in_m): - return x_in_m * 1e3 - + # Cast to a native python float to prevent "NP.FLOAT64(1.0)" in xml + return float(x_in_m) * 1e3 + def s_to_ms(x_in_s): - return x_in_s * 1e3 + return float(x_in_s) * 1e3 mrd_seq_params = ismrmrd.xsd.sequenceParametersType()