.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "generated/gallery/02-tours/06_spiral.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_generated_gallery_02-tours_06_spiral.py: ============== Spiral readout ============== A spiral readout acquires k-space along a spiral arm instead of along straight lines or spokes. The waveform of an arm is not a trapezoid: it is solved numerically against the system limits by :class:`~pypulseqpp.sequences.SpiralReadout2D`, and in this Tour you use the module only for the arms it designs. Three limits bound the traversal of an arm. Two belong to the gradient system, the maximum amplitude and the maximum slew rate. The third follows from the receiver: with a dwell time :math:`\Delta t` the trajectory may not advance further than :math:`1/\mathrm{FOV}` between samples, which caps the gradient amplitude at .. math:: G_\mathrm{bw} = \frac{1}{\gamma\, \Delta t\, \mathrm{FOV}} whatever the gradient system could deliver. The solver applies the lowest of the three, so the readout duration depends on the slew rate over part of the design space and not over the rest. You read the binding limit of each design from the waveform it produces, over a range of slew limits and sampling rates, and then relate the interleaf count to the duration of an arm and to the duration of a full set of interleaves. **Prerequisites:** lessons 5 and 9 of the :doc:`course `. .. GENERATED FROM PYTHON SOURCE LINES 31-151 .. code-block:: Python #: Gyromagnetic ratio of the proton (Hz/T). GAMMA = 42.576e6 FOV = 220e-3 MATRIX = 128 MAX_GRAD_MT_M = 40.0 .. GENERATED FROM PYTHON SOURCE LINES 152-160 Arm design ---------- The readout module designs the arm from the prescription and the system limits it is given, and stores the resulting gradient waveform as an attribute. ``design_interleaves`` sets the pitch of the spiral, against which the readout duration is measured. It is not the number of arms a scan plays. .. GENERATED FROM PYTHON SOURCE LINES 160-191 .. code-block:: Python import numpy as np import pypulseqpp as pp from pypulseqpp import sequences def spiral_arm(max_slew, sampling_rate_hz, interleaves=16): """One spiral arm designed at the given slew limit and sampling rate.""" system = pp.Opts( max_grad=MAX_GRAD_MT_M, grad_unit="mT/m", max_slew=max_slew, slew_unit="T/m/s", rf_dead_time=100e-6, rf_ringdown_time=30e-6, adc_dead_time=10e-6, ) excitation = sequences.SpatialSelectiveExcitation(system, 15.0, 5e-3) # The module solves against the package's derated limits, so those are the # ceilings a measurement of the waveform has to be read against. return pp.apply_system_derates(system), sequences.SpiralReadout2D( system, excitation.rf, excitation.gz, fov=FOV, matrix=MATRIX, design_interleaves=interleaves, readout_bandwidth_hz=sampling_rate_hz, ) .. GENERATED FROM PYTHON SOURCE LINES 192-195 The amplitude and slew the design reached are measured on the waveform it wrote, on its own raster and along the vector rather than per axis, because the two in-plane axes play at once. .. GENERATED FROM PYTHON SOURCE LINES 195-231 .. code-block:: Python def measure(system, arm): """Duration, waveform and the ceilings the design was bounded by.""" x = np.asarray(arm.gx.waveform) y = np.asarray(arm.gy.waveform) times = np.asarray(arm.gx.tt) raster = float(times[1] - times[0]) return { "readout": arm.n_samples / arm.bandwidth_hz, "time": times - times[0], "magnitude": np.hypot(x, y) / GAMMA, "slew": np.hypot(np.diff(x), np.diff(y)) / GAMMA / raster, "ceiling_grad": system.max_grad / GAMMA, "ceiling_slew": system.max_slew / GAMMA, "ceiling_bw": arm.bandwidth_hz / (GAMMA * FOV), } def limiting_ceiling(measured, tolerance=0.02): """Which ceiling caps the arm's amplitude, or ``'slew'`` if none of them does. An arm is at its slew limit wherever it is turning, so reaching that limit says nothing on its own. What decides whether more slew rate would shorten the traversal is whether the amplitude has reached a ceiling. """ cap = min(measured["ceiling_grad"], measured["ceiling_bw"]) if measured["magnitude"].max() < (1.0 - tolerance) * cap: return "slew" return ( "bandwidth" if measured["ceiling_bw"] < measured["ceiling_grad"] else "amplitude" ) .. GENERATED FROM PYTHON SOURCE LINES 232-238 Design space ------------ The slew limit is swept over the range a body gradient system covers, and the sampling rate over a range whose receiver cap runs from well below the gradient amplitude limit to above it. .. GENERATED FROM PYTHON SOURCE LINES 238-257 .. code-block:: Python SLEWS = (40.0, 60.0, 80.0, 100.0, 120.0, 150.0, 180.0, 210.0) RATES = (100e3, 250e3, 600e3) grid = {} for rate in RATES: grid[rate] = [] for slew in SLEWS: system, arm = spiral_arm(slew, rate) measured = measure(system, arm) grid[rate].append( {"slew_limit": slew, "regime": limiting_ceiling(measured), **measured} ) .. image-sg:: /generated/gallery/02-tours/images/sphx_glr_06_spiral_001.png :alt: 06 spiral :srcset: /generated/gallery/02-tours/images/sphx_glr_06_spiral_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none slew 100 kHz 250 kHz 600 kHz 40 10.96 ms band 7.49 ms slew 7.50 ms slew 60 10.96 ms band 6.11 ms slew 6.12 ms slew 80 10.96 ms band 5.30 ms slew 5.30 ms slew 100 10.96 ms band 4.74 ms slew 4.74 ms slew 120 11.00 ms band 4.40 ms band 4.34 ms slew 150 11.00 ms band 4.40 ms band 3.88 ms slew 180 11.00 ms band 4.40 ms band 3.54 ms slew 210 10.92 ms band 4.35 ms band 3.26 ms ampl .. GENERATED FROM PYTHON SOURCE LINES 258-273 The three rates behave differently. At the lowest, the receiver's cap is so far below the gradient amplitude limit that the arm reaches it within the first turn, and the duration is the same across the whole slew range. At the middle rate the amplitude climbs at the slew limit until it meets the receiver's cap, after which raising the limit changes the duration by about a percent. At the highest rate the receiver's cap is above the gradient amplitude limit, so the amplitude the arm settles at is the hardware's, and it is only reached at the top of the slew range; below that the arm is still climbing when it ends. The flat part of the middle curve is not exactly flat, and the reason is that an arm at constant amplitude is still turning. Holding :math:`|G|` while the direction rotates requires slew rate of its own, and the tighter the turn the more of it, so the slew limit continues to govern the first turns of an arm whose amplitude has already stopped growing. .. GENERATED FROM PYTHON SOURCE LINES 275-279 Waveform regimes ---------------- One design from each regime, with the ceilings drawn on the axes they bound. .. GENERATED FROM PYTHON SOURCE LINES 279-290 .. image-sg:: /generated/gallery/02-tours/images/sphx_glr_06_spiral_002.png :alt: slew-limited, receiver cap, amplitude cap :srcset: /generated/gallery/02-tours/images/sphx_glr_06_spiral_002.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 291-302 In the slew-limited design the amplitude is still climbing when the arm ends. In the other two it reaches a ceiling part way out and stays there, and the slew falls away from its limit once it does: the remaining traversal is at constant speed, and the only turning left is the angular one. The slew rate is at its limit early in every one of them, so reaching the slew limit does not by itself distinguish the three. The ceilings are drawn at the system's derated limits rather than at the numbers passed in. A design whose two in-plane axes play together is solved against a per-axis limit reduced by :math:`\sqrt{2}`, so that the vector magnitude drawn here respects the scalar limit. .. GENERATED FROM PYTHON SOURCE LINES 304-310 Interleaves against arm duration -------------------------------- Within one regime the pitch is the remaining free parameter: more interleaves cover k-space with shorter arms, and the set of them takes correspondingly longer to play. .. GENERATED FROM PYTHON SOURCE LINES 310-333 .. code-block:: Python interleaves = [] for count in (4, 8, 16, 32, 48): system, arm = spiral_arm(150.0, 250e3, interleaves=count) interleaves.append( { "interleaves": count, "readout": arm.n_samples / arm.bandwidth_hz, "scan": count * arm.duration, } ) .. image-sg:: /generated/gallery/02-tours/images/sphx_glr_06_spiral_003.png :alt: 06 spiral :srcset: /generated/gallery/02-tours/images/sphx_glr_06_spiral_003.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none arms readout per arm full set 4 17.15 ms 20.96 ms 0.084 s 8 8.61 ms 12.42 ms 0.099 s 16 4.40 ms 8.20 ms 0.131 s 32 2.27 ms 6.08 ms 0.195 s 48 1.60 ms 5.60 ms 0.269 s .. GENERATED FROM PYTHON SOURCE LINES 334-357 The arm duration falls almost as the reciprocal of the interleaf count while the time for a full set rises less than proportionally, because each repetition contains an excitation and a rewind whose duration does not depend on the pitch. Off-resonance and :math:`T_2^*` act over the readout duration, so the interleaf count is the remaining way to shorten it once the slew rate no longer does. Specification ------------- What this Tour built, stated the way you would ask an agent for it: .. code-block:: text With pypulseqpp, design a 2D spiral arm for a 220 mm FOV and a 128 matrix on a 40 mT/m system: a 15 degree SpatialSelectiveExcitation of 5 ms followed by SpiralReadout2D with design_interleaves=16. Measure the vector amplitude and slew rate of the arm's waveform against the derated limits and the receiver cap 1 / (gamma * dwell * FOV). Sweep the slew limit from 40 to 210 T/m/s at 100, 250 and 600 kHz sampling and report which limit binds. Then sweep the design interleaves from 4 to 48 at 150 T/m/s and 250 kHz and report the readout duration and the time of a full set of arms. .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 1.195 seconds) .. _sphx_glr_download_generated_gallery_02-tours_06_spiral.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: 06_spiral.ipynb <06_spiral.ipynb>` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: 06_spiral.py <06_spiral.py>` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: 06_spiral.zip <06_spiral.zip>` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_