.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_examples/01-basics/02-from-kspace-to-image.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_auto_examples_01-basics_02-from-kspace-to-image.py: ===================== From k-space to image ===================== This lesson reconstructs an undersampled Cartesian brain acquisition from its multichannel k-space to a coil-combined image, and shows what each step of a parallel-imaging and compressed-sensing pipeline contributes. Scan time in Cartesian MRI is proportional to the number of phase-encoding lines; skipping lines shortens the scan by the acceleration factor :math:`R`, but violates the Nyquist criterion and folds the image onto itself. Recovering an unaliased image from such data is what the receive coil array, and prior knowledge of the image, are used for. The acquisition is simulated from a BrainWeb tissue segmentation and the eight channels of BART's head-coil model, with one line in three acquired along the phase-encoding direction. The pipeline consists of coil compression, sensitivity calibration by ESPIRiT, and a regularized least-squares fit of the SENSE forward model .. math:: y = P F S x + \varepsilon, with :math:`S` the coil sensitivities, :math:`F` the Fourier transform, :math:`P` the sampling operator that keeps the acquired phase encodes, and :math:`\varepsilon` complex Gaussian noise. :doc:`../../explanation/encoding` states the model and :doc:`../../explanation/inverse-problems` the estimator. Shapes are C order, so a Cartesian k-space is ``(coils, z, y, x)`` with the readout along ``x`` and the phase encoding along ``y``; see :doc:`/explanation/data-layout`. **Learning objectives** - Simulate a multichannel Cartesian acquisition from a tissue segmentation. - Undersample the phase-encoding direction with a variable-density pattern around a fully sampled autocalibration (ACS) region. - Compress the channels with :func:`bartorch.tools.cc` and estimate their sensitivities with :func:`bartorch.tools.ecalib`. - Reconstruct with :func:`bartorch.apps.pics`, with a Tikhonov and with a wavelet sparsity penalty, and compare the results by error maps, NRMSE and SSIM. It builds on the conventions of :doc:`01-tensors-and-commands`. The sections after it examine calibration, regularization and the operator form of each step in turn; the next lesson, :doc:`../02-parallel-imaging/01-coil-calibration`, compares sensitivity estimators. .. GENERATED FROM PYTHON SOURCE LINES 52-154 .. code-block:: Python import csv from pathlib import Path import brainweb_dl import numpy as np import torch from brainweb_dl import get_mri import bartorch import bartorch.tools as bt from bartorch import apps, priors .. GENERATED FROM PYTHON SOURCE LINES 155-164 Phantom ------- BrainWeb [#brainweb]_ publishes a segmentation rather than an image: one membership map per tissue class, from which a table of relaxation times and proton densities gives the signal of a chosen acquisition. The volume ``brainweb-dl`` returns is indexed ``(inferior-superior, posterior-anterior, left-right)``, so its first axis selects an axial slice, and an image is drawn from its first row down, so flipping it puts anterior at the top. .. GENERATED FROM PYTHON SOURCE LINES 165-179 .. code-block:: Python SIZE = 192 COILS = 8 SLICE = 90 # axial, through the lateral ventricles TISSUES = (1, 2, 3, 4, 5, 6, 8) # everything the table gives relaxation times table = Path(brainweb_dl.__file__).parent / "data" / "brainweb1_tissues.csv" entries = list(csv.DictReader(table.open())) tissue_t1 = np.array([float(row["T1 (ms)"]) for row in entries], dtype=np.float32)[list(TISSUES)] tissue_t2 = np.array([float(row["T2 (ms)"]) for row in entries], dtype=np.float32)[list(TISSUES)] tissue_pd = np.array([float(row["PD (ms)"]) for row in entries], dtype=np.float32)[list(TISSUES)] fractions = np.flipud(get_mri(sub_id=0, contrast="fuzzy")[SLICE])[..., list(TISSUES)].copy() .. GENERATED FROM PYTHON SOURCE LINES 180-184 The slice is cropped to a square field of view around the head and resampled to the matrix reconstructed here. The crop leaves a margin, as a real field of view does: the aliased copies of an undersampled acquisition then fall partly outside the head. .. GENERATED FROM PYTHON SOURCE LINES 185-207 .. code-block:: Python MARGIN = 0.25 .. GENERATED FROM PYTHON SOURCE LINES 208-218 A membership-weighted average of the table gives :math:`T_1`, :math:`T_2` and the proton density at every voxel, and the spin-echo signal .. math:: S = \rho \, \left(1 - e^{-T_R/T_1}\right) e^{-T_E/T_2} turns those into the image the experiment measures. At a short repetition time and a short echo time the contrast is :math:`T_1`-weighted: white matter bright, cerebrospinal fluid dark, subcutaneous fat brightest of all. .. GENERATED FROM PYTHON SOURCE LINES 219-241 .. code-block:: Python TR, TE = 600.0, 12.0 # ms weights = memberships * torch.as_tensor(tissue_pd)[:, None, None] share = weights.sum(0).clamp(min=1e-6) T1 = (weights * torch.as_tensor(tissue_t1)[:, None, None]).sum(0) / share T2 = (weights * torch.as_tensor(tissue_t2)[:, None, None]).sum(0) / share proton_density = weights.sum(0) / weights.sum(0).max() signal = ( proton_density * (1 - torch.exp(-TR / T1.clamp(min=1e-3))) * torch.exp(-TE / T2.clamp(min=1e-3)) ) signal = torch.where(T1 > 0, signal, torch.zeros(())) signal = signal / signal.max() # A smooth quadratic phase stands in for the transmit and off-resonance phase # of a real object, so that nothing below depends on the image being real. grid_y, grid_x = torch.meshgrid( torch.linspace(-1.0, 1.0, SIZE), torch.linspace(-1.0, 1.0, SIZE), indexing="ij" ) image = (signal * torch.exp(0.8j * (grid_x**2 - 0.5 * grid_y**2))).to(torch.complex64) .. GENERATED FROM PYTHON SOURCE LINES 242-253 Coils ----- Each receive channel measures the object weighted by its complex sensitivity profile, :math:`x_c = S_c x`. The sensitivities here are BART's analytical head coil, evaluated on the image grid that :func:`bartorch.tools.grid` describes. Dividing them by their root sum of squares over the channels normalizes :math:`\sum_c |S_c|^2` to one, so that the optimal coil combination of the coil images is the image itself and a reconstruction can be compared against it directly. Complex Gaussian noise is then added to every k-space sample, as thermal noise is in the receiver chain. .. GENERATED FROM PYTHON SOURCE LINES 254-262 .. code-block:: Python sensitivities = bt.coils(t=bt.grid(D=(SIZE, SIZE, 1)), n=COILS)[:, 0] sensitivities = sensitivities / bartorch.rss(sensitivities, axes=(0,), keepdim=True) coil_images = sensitivities * image kspace = bartorch.fft(coil_images, axes=(-2, -1), unitary=True) kspace = bt.noise(kspace, n=1e-4, s=42) .. GENERATED FROM PYTHON SOURCE LINES 263-300 .. rst-class:: sphx-glr-horizontal * .. image-sg:: /auto_examples/01-basics/images/sphx_glr_02-from-kspace-to-image_001.png :alt: $T_1$, $T_2$, proton density, $T_1$-weighted image :srcset: /auto_examples/01-basics/images/sphx_glr_02-from-kspace-to-image_001.png :class: sphx-glr-multi-img * .. image-sg:: /auto_examples/01-basics/images/sphx_glr_02-from-kspace-to-image_002.png :alt: channel 2, channel 4, channel 6 :srcset: /auto_examples/01-basics/images/sphx_glr_02-from-kspace-to-image_002.png :class: sphx-glr-multi-img .. GENERATED FROM PYTHON SOURCE LINES 301-322 The first figure is the ground truth: the relaxation maps, drawn with the perceptually uniform colormaps recommended for relaxometry [#fuderer]_ -- lipari for :math:`T_1`, navia for :math:`T_2` -- and with a window that stops short of cerebrospinal fluid, the proton density, and the :math:`T_1`-weighted image they give. The second shows three of the eight sensitivities, magnitude above and phase below, with the outline of the head. Each magnitude is highest near its coil element and falls off across the head; the phase varies smoothly. These spatial variations are the extra encoding that parallel imaging uses to separate aliased voxels. Sampling -------- The readout is fully sampled, since it costs no scan time, and a subset of the phase encodes is acquired. The lines are drawn at random from a variable density that is highest at the k-space centre, where most of the signal energy is, with a block of 24 central lines, the autocalibration signal (ACS) region, acquired in full. ESPIRiT reads its calibration matrix from the ACS region, so an acquisition without one would need a separate calibration scan. The pattern is a column vector along the phase-encoding direction: it broadcasts over the readout and over the channels. .. GENERATED FROM PYTHON SOURCE LINES 323-344 .. code-block:: Python ACCELERATION = 3 CALIBRATION = 24 encodes = torch.arange(SIZE) - SIZE // 2 profile = (1.0 + 2.0 * encodes.abs() / SIZE) ** -3.0 centre = (encodes.abs() < CALIBRATION // 2).to(torch.float32) drawn = torch.multinomial( profile * (1.0 - centre), SIZE // ACCELERATION - CALIBRATION, replacement=False, generator=torch.Generator().manual_seed(11), ) lines = centre.clone() lines[drawn] = 1.0 pattern = lines.reshape(SIZE, 1).to(torch.complex64) measured = kspace[:, None] * pattern print(f"{int(lines.sum())} of {SIZE} phase encodes acquired, R = {SIZE / lines.sum():.1f}") .. rst-class:: sphx-glr-script-out .. code-block:: none 63 of 192 phase encodes acquired, R = 3.0 .. GENERATED FROM PYTHON SOURCE LINES 345-366 .. image-sg:: /auto_examples/01-basics/images/sphx_glr_02-from-kspace-to-image_003.png :alt: sampling pattern, acquired k-space, channel 0 :srcset: /auto_examples/01-basics/images/sphx_glr_02-from-kspace-to-image_003.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 367-381 In the pattern (readout horizontal, phase encoding vertical) every acquired phase encode is a full line; the lines cluster towards the centre and the ACS band is dense. Channel compression ------------------- Eight channels carry less independent information than eight images: the sensitivities overlap, and the singular value spectrum of the calibration matrix falls off. :func:`bartorch.tools.cc` returns the matrix that projects the channels onto their leading singular vectors [#huangcc]_, the virtual coils, and :func:`bartorch.tools.ccapply` applies it. Calibration, the encoding operator and every iteration then cost six channels rather than eight, at a negligible loss of the encoding capacity of the array. .. GENERATED FROM PYTHON SOURCE LINES 382-388 .. code-block:: Python VIRTUAL = 6 matrix = bt.cc(measured, p=VIRTUAL, M=True, r=CALIBRATION) compressed = bt.ccapply(measured, matrix, p=VIRTUAL) .. GENERATED FROM PYTHON SOURCE LINES 389-399 Sensitivity calibration ----------------------- ESPIRiT [#espirit]_ estimates the sensitivities from the ACS region alone: it builds a calibration matrix from all k-space neighbourhoods (kernels) in the region, and obtains the sensitivities at each voxel as the eigenvector of an operator derived from that matrix whose eigenvalue is one. Outside the object no eigenvalue is close to one; ``crop`` sets the maps to zero where the eigenvalue falls below it, which keeps the background out of the reconstruction. .. GENERATED FROM PYTHON SOURCE LINES 400-403 .. code-block:: Python maps = bt.ecalib(compressed, maps=1, calib_size=CALIBRATION, crop=0.8) .. GENERATED FROM PYTHON SOURCE LINES 404-421 Reconstruction -------------- Three reconstructions of the same data are compared. - The **zero-filled** reconstruction sets the missing phase encodes to zero, inverse-transforms each channel and combines them by root sum of squares. It uses no model of the encoding, so every missing line leaves aliasing. - **SENSE** [#sense]_ solves :math:`\min_x \|PFSx - y\|_2^2 + \lambda\|x\|_2^2` by conjugate gradients. The sensitivities unfold the aliasing, but the inversion amplifies the noise by the g-factor, which is highest where the coils cannot distinguish aliased voxels. - **Compressed sensing** [#lustig]_ replaces the Tikhonov term by an :math:`\ell_1` penalty on the wavelet coefficients, solved by FISTA [#beck]_. The random undersampling makes the aliasing incoherent, i.e. noise-like in the wavelet domain, and the sparsity penalty removes it together with the amplified noise. .. GENERATED FROM PYTHON SOURCE LINES 422-435 .. code-block:: Python channel_images = bartorch.ifft(compressed[:, 0], axes=(-2, -1), unitary=True) zero_filled = bartorch.rss(channel_images, axes=(0,)) sense = apps.pics(compressed, maps, l2=0.001, maxiter=60) wavelet = apps.pics( compressed, maps, regularizers=priors.Wavelet((-1, -2), 0.004), solver="fista", maxiter=100, ) .. GENERATED FROM PYTHON SOURCE LINES 436-442 The sensitivities ESPIRiT estimates and the ones the acquisition was simulated with differ by a phase that varies from voxel to voxel, so the reconstructed image does too, and the comparison is between magnitudes. ``pics`` returns the image in the units of the data it scaled, so :func:`bartorch.tools.nrmse` is called with ``scaled=True``, which fits a global factor before comparing. .. GENERATED FROM PYTHON SOURCE LINES 443-450 .. code-block:: Python results = {"zero-filled": zero_filled, "SENSE": sense, "wavelet CS": wavelet} for name, estimate in results.items(): error = bt.nrmse(image.abs(), estimate.abs(), scaled=True) similarity = bt.ssim(image.abs(), scaled(estimate, image)) print(f"{name:>12} NRMSE {error:.3f} SSIM {similarity:.3f}") .. rst-class:: sphx-glr-script-out .. code-block:: none zero-filled NRMSE 0.129 SSIM 0.541 SENSE NRMSE 0.101 SSIM 0.737 wavelet CS NRMSE 0.044 SSIM 0.926 .. GENERATED FROM PYTHON SOURCE LINES 451-474 .. rst-class:: sphx-glr-horizontal * .. image-sg:: /auto_examples/01-basics/images/sphx_glr_02-from-kspace-to-image_004.png :alt: reference, zero-filled, SENSE, wavelet CS :srcset: /auto_examples/01-basics/images/sphx_glr_02-from-kspace-to-image_004.png :class: sphx-glr-multi-img * .. image-sg:: /auto_examples/01-basics/images/sphx_glr_02-from-kspace-to-image_005.png :alt: zero-filled error, SENSE error, wavelet CS error :srcset: /auto_examples/01-basics/images/sphx_glr_02-from-kspace-to-image_005.png :class: sphx-glr-multi-img * .. image-sg:: /auto_examples/01-basics/images/sphx_glr_02-from-kspace-to-image_006.png :alt: reference, enlarged, zero-filled, SENSE, wavelet CS :srcset: /auto_examples/01-basics/images/sphx_glr_02-from-kspace-to-image_006.png :class: sphx-glr-multi-img .. GENERATED FROM PYTHON SOURCE LINES 475-489 The zero-filled image carries the aliasing of the missing phase encodes as blurring and ghosting along the vertical, phase-encoding direction. SENSE removes the coherent aliasing, but its error map shows noise amplified in the centre of the head, where the coil sensitivities are least distinct, and incoherent residual artefacts of the random sampling. The wavelet penalty suppresses both; in the enlarged region the cortical folding and the ventricle boundaries are sharper and the background of the brain is smooth. The NRMSE and SSIM printed above quantify the same ordering. How much the penalty removes depends on its weight, which is chosen here and not estimated: a larger weight removes more noise and more fine texture with it. :doc:`../02-parallel-imaging/01-coil-calibration` compares sensitivity estimators, and :doc:`../03-regularization/01-regularized-reconstruction` varies the weight. .. GENERATED FROM PYTHON SOURCE LINES 492-525 References ---------- .. [#brainweb] Collins DL, Zijdenbos AP, Kollokian V, Sled JG, Kabani NJ, Holmes CJ, Evans AC. Design and construction of a realistic digital brain phantom. *IEEE Trans Med Imaging* 17(3):463-468 (1998). https://doi.org/10.1109/42.712135 .. [#fuderer] Fuderer M, Wichtmann B, Crameri F, de Souza NM, Baeßler B, Gulani V, et al. Color-map recommendation for MR relaxometry maps. *Magn Reson Med* 93(2):490-506 (2025). https://doi.org/10.1002/mrm.30290 .. [#huangcc] Huang F, Vijayakumar S, Li Y, Hertel S, Duensing GR. A software channel compression technique for faster reconstruction with many channels. *Magn Reson Imaging* 26(1):133-141 (2008). https://doi.org/10.1016/j.mri.2007.04.010 .. [#espirit] Uecker M, Lai P, Murphy MJ, Virtue P, Elad M, Pauly JM, Vasanawala SS, Lustig M. ESPIRiT -- an eigenvalue approach to autocalibrating parallel MRI: where SENSE meets GRAPPA. *Magn Reson Med* 71(3):990-1001 (2014). https://doi.org/10.1002/mrm.24751 .. [#sense] Pruessmann KP, Weiger M, Scheidegger MB, Boesiger P. SENSE: sensitivity encoding for fast MRI. *Magn Reson Med* 42(5):952-962 (1999). https://doi.org/10.1002/(SICI)1522-2594(199911)42:5%3C952::AID-MRM16%3E3.0.CO;2-S .. [#lustig] Lustig M, Donoho D, Pauly JM. Sparse MRI: the application of compressed sensing for rapid MR imaging. *Magn Reson Med* 58(6):1182-1195 (2007). https://doi.org/10.1002/mrm.21391 .. [#beck] Beck A, Teboulle M. A fast iterative shrinkage-thresholding algorithm for linear inverse problems. *SIAM J Imaging Sci* 2(1):183-202 (2009). https://doi.org/10.1137/080716542 .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 2.496 seconds) .. _sphx_glr_download_auto_examples_01-basics_02-from-kspace-to-image.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: 02-from-kspace-to-image.ipynb <02-from-kspace-to-image.ipynb>` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: 02-from-kspace-to-image.py <02-from-kspace-to-image.py>` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: 02-from-kspace-to-image.zip <02-from-kspace-to-image.zip>` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_