.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "generated/gallery/04-non-cartesian/01_radial.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_04-non-cartesian_01_radial.py: =============== Radial sampling =============== The Cartesian gradient echo of :doc:`/generated/gallery/01-pulseq-basics/03_gradient_echo` changes the acquired line with a phase encode. This lesson replaces the phase encode with a rotation of the readout gradient itself, so that every repetition acquires a spoke through the centre of k-space. It establishes how many spokes such an acquisition requires and what ordering them by the golden angle changes. The measured quantity is the azimuthal gap between neighbouring spokes at the edge of k-space, computed from the sampling locations of the sequence and compared with the sample spacing along a spoke that the prescription sets. Learning objectives ------------------- After this lesson, you should be able to: - rotate a readout gradient and its prewinder into the imaging plane; - build one repetition per spoke from the Cartesian gradient echo; - measure the azimuthal gap at the edge of k-space from the sampling locations, and relate it to the Nyquist spoke count :math:`P = \tfrac{\pi}{2} N`; - compare uniform and golden-angle orderings of truncated acquisitions. .. GENERATED FROM PYTHON SOURCE LINES 29-39 .. GENERATED FROM PYTHON SOURCE LINES 40-49 A rotated readout ----------------- A spoke runs from one edge of k-space through the centre to the other, so the prewinder has half the readout area as it does on a Cartesian line, and the echo is at the middle of the acquisition window. The pair is then rotated about the slice axis by the angle of the spoke, which :func:`~pypulseqpp.rotate` does by resolving each gradient onto the two in-plane axes. .. GENERATED FROM PYTHON SOURCE LINES 49-109 .. code-block:: Python import numpy as np import pypulseqpp as pp system = pp.Opts( max_grad=32.0, grad_unit="mT/m", max_slew=130.0, slew_unit="T/m/s", rf_dead_time=100e-6, rf_ringdown_time=20e-6, adc_dead_time=10e-6, ) FOV = 220e-3 MATRIX = 128 THICKNESS = 5e-3 FLIP_ANGLE_DEG = 12.0 REPETITION_TIME = 10e-3 DWELL = 10e-6 rf, gz, gz_reph = pp.make_sinc_pulse( flip_angle=np.deg2rad(FLIP_ANGLE_DEG), duration=2e-3, slice_thickness=THICKNESS, apodization=0.5, time_bw_product=4.0, delay=system.rf_dead_time, system=system, use="excitation", return_gz=True, ) acquisition = MATRIX * DWELL raster = system.grad_raster_time gx = pp.make_trapezoid( channel="x", amplitude=MATRIX / FOV / acquisition, flat_time=raster * np.ceil(acquisition / raster), system=system, ) adc = pp.make_adc(num_samples=MATRIX, dwell=DWELL, delay=gx.rise_time, system=system) gx_pre = pp.make_trapezoid( channel="x", area=-(gx.amplitude * gx.rise_time / 2 + (MATRIX / 2 + 0.5) / FOV), duration=5e-4, system=system, ) spoiler = pp.make_crusher(4.0, FOV / MATRIX, channel="z", system=system)[0] # Both in-plane axes have a gradient for every spoke, but the vector amplitude # of a spoke is the readout amplitude at every angle, and each axis plays its # projection of it, so the per-axis limit binds where a spoke lies along an # axis and nowhere else. print( f"readout {1e3 * gx.amplitude / 42.576e6:.2f} mT/m, " f"per-axis limit {1e3 * system.max_grad / 42.576e6:.1f} mT/m" ) .. rst-class:: sphx-glr-script-out .. code-block:: none readout 10.68 mT/m, per-axis limit 32.0 mT/m .. GENERATED FROM PYTHON SOURCE LINES 110-115 One repetition per spoke ------------------------ The repetition is the Cartesian one with the phase encode removed and the readout rotated. Both in-plane axes have a gradient for every spoke. .. GENERATED FROM PYTHON SOURCE LINES 115-155 .. code-block:: Python GOLDEN_ANGLE = np.pi * (3.0 - np.sqrt(5.0)) / 2.0 def radial(spokes, golden=False): """A radial acquisition of the given number of spokes.""" if golden: angles = (np.arange(spokes) * GOLDEN_ANGLE) % np.pi else: angles = np.arange(spokes) * np.pi / spokes played = ( pp.calc_duration(rf, gz) + pp.calc_duration(gx_pre, gz_reph) + pp.calc_duration(gx) + pp.calc_duration(spoiler) ) recovery = pp.make_delay( pp.round_to_raster(REPETITION_TIME - played, system.block_duration_raster) ) seq = pp.Sequence(system=system) for angle in angles: seq.add_block(rf, gz) seq.add_block(*pp.rotate(gx_pre, angle=angle, axis="z"), gz_reph) seq.add_block(*pp.rotate(gx, angle=angle, axis="z"), adc) seq.add_block(spoiler) seq.add_block(recovery) return seq, angles SPOKES = 128 seq, angles = radial(SPOKES) ok, errors = seq.check_timing() print( f"timing {ok}, {seq.num_blocks} blocks, {SPOKES} spokes, {seq.duration()[0]:.3f} s" ) seq.paper_plot(tr=1) .. image-sg:: /generated/gallery/04-non-cartesian/images/sphx_glr_01_radial_001.png :alt: 01 radial :srcset: /generated/gallery/04-non-cartesian/images/sphx_glr_01_radial_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none timing True, 640 blocks, 128 spokes, 1.280 s .. GENERATED FROM PYTHON SOURCE LINES 156-161 The trajectory -------------- Every spoke passes through the centre of k-space, so the centre is sampled once per repetition and the periphery only where a spoke reaches it. .. GENERATED FROM PYTHON SOURCE LINES 161-164 .. code-block:: Python pp.plot.plot_kspace(seq, color_by="shot", plane="xy", show_trajectory=False) .. image-sg:: /generated/gallery/04-non-cartesian/images/sphx_glr_01_radial_002.png :alt: 01 radial :srcset: /generated/gallery/04-non-cartesian/images/sphx_glr_01_radial_002.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 165-183 Spokes against azimuthal gap ---------------------------- Along a spoke the samples are :math:`1/\mathrm{FOV}` apart, as the field of view requires. Between spokes the spacing grows with the distance from the centre, and at the edge it is the azimuthal arc between neighbouring spokes, .. math:: \Delta k_\varphi = \frac{\pi}{P} \, k_\mathrm{max} , for :math:`P` spokes over half a turn. Requiring it to be no larger than the radial spacing gives the familiar :math:`P = \tfrac{\pi}{2} N`: a radial acquisition needs more spokes than a Cartesian acquisition of the same matrix needs lines. The gap is measured from the sampling locations the sequence produces. .. GENERATED FROM PYTHON SOURCE LINES 183-248 .. code-block:: Python COUNTS = (32, 64, 128, 201, 256) radial_spacing = 1.0 / FOV def spoke_angles(sequence, spokes): """The angle of every spoke, from the samples the analysis reports.""" sampled = sequence.calculate_kspacePP()[0][:2].reshape(2, spokes, MATRIX) outermost = sampled[:, :, -1] return np.arctan2(outermost[1], outermost[0]) def largest_gap(sampled_angles): """The arc between neighbouring spokes at the edge, in 1/m.""" sorted_angles = np.sort(sampled_angles % np.pi) steps = np.diff(np.concatenate([sorted_angles, [sorted_angles[0] + np.pi]])) return steps.max() * MATRIX / (2 * FOV) uniform = [] for count in COUNTS: sampled, _ = radial(count) uniform.append( { "spokes": count, "gap": largest_gap(spoke_angles(sampled, count)), "scan": count * REPETITION_TIME, } ) nyquist = int(np.ceil(np.pi / 2 * MATRIX)) .. image-sg:: /generated/gallery/04-non-cartesian/images/sphx_glr_01_radial_003.png :alt: 01 radial :srcset: /generated/gallery/04-non-cartesian/images/sphx_glr_01_radial_003.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none spokes gap at the edge of the radial spacing 32 28.56 1/m 6.28 64 14.28 1/m 3.14 128 7.14 1/m 1.57 201 4.55 1/m 1.00 256 3.57 1/m 0.79 radial sample spacing 4.55 1/m, Nyquist at 202 spokes .. GENERATED FROM PYTHON SOURCE LINES 249-254 The measured gap falls as the reciprocal of the spoke count and crosses the radial sample spacing at the predicted count. An acquisition below it is undersampled at the periphery and not at the centre, which is why the artefact it produces is a streak from the edge of the object rather than the fold-over a Cartesian acquisition produces. .. GENERATED FROM PYTHON SOURCE LINES 256-263 Golden-angle ordering --------------------- Advancing the angle by :math:`\pi` times the golden ratio conjugate instead of by :math:`\pi/P` gives an ordering whose every prefix is nearly uniform, so the acquisition can be stopped, or divided into frames, at any length. The gap of a prefix is, however, never exactly that of the uniform ordering. .. GENERATED FROM PYTHON SOURCE LINES 263-292 .. code-block:: Python PREFIXES = np.arange(8, 257, 8) golden, _ = radial(PREFIXES.max(), golden=True) golden_angles = spoke_angles(golden, PREFIXES.max()) golden_gap = np.array([largest_gap(golden_angles[:count]) for count in PREFIXES]) uniform_gap = np.array( [largest_gap(np.arange(count) * np.pi / count) for count in PREFIXES] ) .. image-sg:: /generated/gallery/04-non-cartesian/images/sphx_glr_01_radial_004.png :alt: 01 radial :srcset: /generated/gallery/04-non-cartesian/images/sphx_glr_01_radial_004.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none golden-angle gap, over prefixes of 8 to 256 spokes: 1.17 to 1.89 times the uniform gap .. GENERATED FROM PYTHON SOURCE LINES 293-299 Every prefix is within a small factor of the uniform ordering of the same length, and no prefix leaves a gap of the kind a truncated uniform ordering would: stopping a uniform acquisition after half its spokes leaves half the angular range unsampled, while stopping a golden-angle one leaves the same range covered at half the density. The complete golden-angle set is not exactly uniform. .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 3.399 seconds) .. _sphx_glr_download_generated_gallery_04-non-cartesian_01_radial.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: 01_radial.ipynb <01_radial.ipynb>` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: 01_radial.py <01_radial.py>` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: 01_radial.zip <01_radial.zip>` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_