.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_examples/02-parallel-imaging/03-noise-prewhitening.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_02-parallel-imaging_03-noise-prewhitening.py: ================== Noise prewhitening ================== This lesson measures how correlated noise between receive channels lowers the signal-to-noise ratio (SNR) of a SENSE reconstruction, and how much of it prewhitening with a noise-only acquisition recovers. The thermal noise of a receive array is correlated between channels, through mutual inductance between the coil elements and shared noise sources in the sample, and its level differs from one channel to another, through the elements' size, loading and preamplifier gain. Least squares is the maximum-likelihood estimator only for white noise, i.e. for a channel noise covariance proportional to the identity. A reconstruction that ignores the covariance :math:`\Psi` weights every channel equally and does not combine them with the optimal SNR [#roemer]_. Prewhitening transforms the data by a matrix :math:`W` with :math:`W \Psi W^H = I`, which makes the channel noise white; the sensitivities are then calibrated on, and the SENSE reconstruction run on, the whitened channels unchanged [#pruessmann]_. :math:`\Psi` is estimated from a noise scan, an acquisition with the RF transmitter off that most vendors run before every protocol. The SNR of the two reconstructions is measured by the pseudo-replica method [#robson]_: the same reconstruction is repeated on independent noise realizations added to one noise-free acquisition, and the standard deviation across repetitions is the noise of each voxel. This is how SNR and g-factor maps are obtained for iterative reconstructions, for which no closed-form noise propagation exists. The phantom and the coil sensitivities are built as in :doc:`../01-basics/02-from-kspace-to-image`; the cell that does it is hidden on this page and present in the script this page can be downloaded as. **Learning objectives** - Estimate the channel noise covariance from a noise scan and whiten the data with :func:`bartorch.tools.whiten`. - Measure an SNR map by the pseudo-replica method. - Quantify the SNR gain that prewhitening gives a SENSE reconstruction. It follows :doc:`02-nonlinear-inversion`. The next section starts with :doc:`../03-regularization/01-regularized-reconstruction`. .. GENERATED FROM PYTHON SOURCE LINES 47-233 .. 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 SIZE = 128 COILS = 8 .. GENERATED FROM PYTHON SOURCE LINES 234-243 Correlated channel noise ------------------------ The noise covariance of the simulation couples the channels by a correlation coefficient falling as :math:`0.5^{|i-j|}` with the distance between their indices, and gives the channels noise standard deviations between 0.6 and 1.6 times a common level. The noise scan measures the same channels without signal; its sample covariance is the estimate of :math:`\Psi` that :func:`bartorch.tools.whiten` inverts. .. GENERATED FROM PYTHON SOURCE LINES 244-264 .. code-block:: Python SIGMA = 0.01 # noise level, relative to the image's peak channels = torch.arange(COILS) levels = torch.linspace(0.6, 1.6, COILS) correlation = 0.5 ** (channels[:, None] - channels[None, :]).abs().float() covariance = levels[:, None] * correlation * levels[None, :] mixing = torch.linalg.cholesky(covariance).to(torch.complex64) generator = torch.Generator().manual_seed(2) def channel_noise(shape): """Complex Gaussian noise with the channel covariance above.""" white = torch.randn(COILS, *shape, generator=generator) white = (white + 1j * torch.randn(COILS, *shape, generator=generator)) / 2**0.5 return SIGMA * torch.einsum("ij,j...->i...", mixing, white.to(torch.complex64)) noise_scan = channel_noise((SIZE, SIZE))[:, None] .. GENERATED FROM PYTHON SOURCE LINES 265-269 The whitening matrix maps the measured covariance to the identity. The measure printed below is the mean magnitude of the off-diagonal covariance entries relative to the mean diagonal entry, which is zero for uncorrelated channels. .. GENERATED FROM PYTHON SOURCE LINES 270-289 .. code-block:: Python def channel_covariance(samples): """Sample covariance of the channels, over every other axis.""" flat = samples.reshape(COILS, -1) return (flat @ flat.conj().T) / flat.shape[1] def off_diagonal(matrix): diagonal = torch.diagonal(matrix).abs() return float((matrix.abs().sum() - diagonal.sum()) / (COILS * (COILS - 1)) / diagonal.mean()) measured = channel_covariance(noise_scan) whitened = channel_covariance(bt.whiten(noise_scan, noise_scan)) print(f"off-diagonal covariance: {off_diagonal(measured):.3f} measured") print(f" {off_diagonal(whitened):.3f} after whitening") .. rst-class:: sphx-glr-script-out .. code-block:: none off-diagonal covariance: 0.202 measured 0.000 after whitening .. GENERATED FROM PYTHON SOURCE LINES 290-309 .. image-sg:: /auto_examples/02-parallel-imaging/images/sphx_glr_03-noise-prewhitening_001.png :alt: measured, whitened :srcset: /auto_examples/02-parallel-imaging/images/sphx_glr_03-noise-prewhitening_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 310-322 The measured covariance has a strong diagonal whose entries grow with the channel index, the unequal noise levels, and off-diagonal bands, the correlation. After whitening it is the identity. Acquisition and calibration --------------------------- Every second phase encode is acquired (:math:`R = 2`), with 24 central lines kept as the ACS region. Each pseudo-replica adds a new noise realization to the same noise-free k-space. The sensitivities are calibrated once per pipeline, from the first replica: from the channels as measured, and from the whitened channels. .. GENERATED FROM PYTHON SOURCE LINES 323-344 .. code-block:: Python CALIBRATION = 24 REPLICAS = 32 lines = torch.zeros(SIZE) lines[::2] = 1.0 lines[SIZE // 2 - CALIBRATION // 2 : SIZE // 2 + CALIBRATION // 2] = 1.0 pattern = lines.reshape(SIZE, 1).to(torch.complex64) noiseless = bartorch.fft(sensitivities * image, axes=(-2, -1), unitary=True) def acquire(): """One replica: the noise-free k-space plus new noise, sampled.""" return ((noiseless + channel_noise((SIZE, SIZE))) * pattern)[:, None] first = acquire() maps_measured = bt.ecalib(first, maps=1, calib_size=CALIBRATION, crop=0.8) maps_whitened = bt.ecalib(bt.whiten(first, noise_scan), maps=1, calib_size=CALIBRATION, crop=0.8) .. GENERATED FROM PYTHON SOURCE LINES 345-351 Pseudo-replicas --------------- Both pipelines use the same reconstruction, conjugate-gradient SENSE with a small Tikhonov weight and a fixed number of iterations; they differ only in whether the data are whitened first. .. GENERATED FROM PYTHON SOURCE LINES 352-361 .. code-block:: Python plain, prewhitened = [], [] for _ in range(REPLICAS): data = acquire() plain.append(apps.pics(data, maps_measured, l2=1e-3, maxiter=30)) prewhitened.append(apps.pics(bt.whiten(data, noise_scan), maps_whitened, l2=1e-3, maxiter=30)) plain, prewhitened = torch.stack(plain), torch.stack(prewhitened) .. GENERATED FROM PYTHON SOURCE LINES 362-365 The SNR of a voxel is the magnitude of the mean reconstruction over the standard deviation across replicas. Both are reported over the white matter, where the phantom is homogeneous. .. GENERATED FROM PYTHON SOURCE LINES 366-382 .. code-block:: Python white_matter = memberships[CLASS["WM"]] > 0.8 def snr_map(replicas): return replicas.mean(0).abs() / replicas.std(0) for name, replicas in (("as measured", plain), ("prewhitened", prewhitened)): values = snr_map(replicas)[white_matter] mean, median = float(values.mean()), float(values.median()) print(f"{name:>12} white-matter SNR {mean:6.1f} (median {median:.1f})") gain = snr_map(prewhitened) / snr_map(plain) print(f"SNR ratio, prewhitened over as measured: {float(gain[white_matter].median()):.2f} (median)") .. rst-class:: sphx-glr-script-out .. code-block:: none as measured white-matter SNR 38.1 (median 37.7) prewhitened white-matter SNR 55.6 (median 54.3) SNR ratio, prewhitened over as measured: 1.42 (median) .. GENERATED FROM PYTHON SOURCE LINES 383-416 .. rst-class:: sphx-glr-horizontal * .. image-sg:: /auto_examples/02-parallel-imaging/images/sphx_glr_03-noise-prewhitening_002.png :alt: SNR, as measured, SNR, prewhitened :srcset: /auto_examples/02-parallel-imaging/images/sphx_glr_03-noise-prewhitening_002.png :class: sphx-glr-multi-img * .. image-sg:: /auto_examples/02-parallel-imaging/images/sphx_glr_03-noise-prewhitening_003.png :alt: SNR ratio :srcset: /auto_examples/02-parallel-imaging/images/sphx_glr_03-noise-prewhitening_003.png :class: sphx-glr-multi-img .. GENERATED FROM PYTHON SOURCE LINES 417-444 Prewhitening raises the SNR throughout the head, and the white-matter histogram shifts by the ratio printed above. The gain varies in space: it is largest in the posterior half of the head, where the channels with the lowest noise level (the first channels here) are most sensitive. Without whitening the least-squares fit weights those channels no more than the noisiest ones; after whitening each channel enters with the weight its noise level warrants. The size of the gain depends on the array's noise correlation and on the spread of its channel noise levels, and is measured here for one simulated covariance; with uncorrelated channels of equal noise level it is one. References ---------- .. [#roemer] Roemer PB, Edelstein WA, Hayes CE, Souza SP, Mueller OM. The NMR phased array. *Magn Reson Med* 16(2):192-225 (1990). https://doi.org/10.1002/mrm.1910160203 .. [#pruessmann] Pruessmann KP, Weiger M, Börnert P, Boesiger P. Advances in sensitivity encoding with arbitrary k-space trajectories. *Magn Reson Med* 46(4):638-651 (2001). https://doi.org/10.1002/mrm.1241 .. [#robson] Robson PM, Grant AK, Madhuranthakam AJ, Lattanzi R, Sodickson DK, McKenzie CA. Comprehensive quantification of signal-to-noise ratio and g-factor for image-based and k-space-based parallel imaging reconstructions. *Magn Reson Med* 60(4):895-907 (2008). https://doi.org/10.1002/mrm.21728 .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 3.179 seconds) .. _sphx_glr_download_auto_examples_02-parallel-imaging_03-noise-prewhitening.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: 03-noise-prewhitening.ipynb <03-noise-prewhitening.ipynb>` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: 03-noise-prewhitening.py <03-noise-prewhitening.py>` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: 03-noise-prewhitening.zip <03-noise-prewhitening.zip>` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_