.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_examples/07-tours/08-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_07-tours_08-noise-prewhitening.py: ================== Noise prewhitening ================== In this tour you measure how much signal-to-noise ratio (SNR) a SENSE reconstruction loses when the noise of the receive channels is correlated, and how much of it you get back by prewhitening with a noise-only acquisition. The thermal noise of a receive array is correlated between channels, through mutual inductance between the coil elements and noise sources shared in the sample, and its level differs from channel to channel, with the elements' size, loading and preamplifier gain. Least squares is the maximum-likelihood estimate only when the noise is white, that is when the channel noise covariance :math:`\Psi` is proportional to the identity. A reconstruction that ignores :math:`\Psi` weights every channel equally, and does not combine them with the best SNR [#roemer]_. Prewhitening multiplies the data by a matrix :math:`W` with :math:`W \Psi W^H = I`, which makes the channel noise white; you then calibrate the sensitivities on the whitened channels and run SENSE on them unchanged [#pruessmann]_. You estimate :math:`\Psi` from a noise scan, an acquisition with the RF transmitter off that most vendors run before every protocol. You measure the SNR by the pseudo-replica method [#robson]_: repeat the same reconstruction on independent noise realisations added to one noise-free acquisition, and take the standard deviation across the repetitions as the noise of each voxel. That is how you get SNR and g-factor maps for an iterative reconstruction, for which there is no closed form of the noise propagation. Prerequisites: lessons 2 and 3, :doc:`../01-basics/02-from-kspace-to-image` and :doc:`../02-parallel-imaging/01-coil-calibration`. **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. .. GENERATED FROM PYTHON SOURCE LINES 46-136 .. GENERATED FROM PYTHON SOURCE LINES 137-146 Correlated channel noise ------------------------ Your object is the BrainWeb slice of lesson 2, seen by BART's analytical eight-channel head coil; how it is built is hidden on this page. The noise of your simulated array couples the channels by a correlation coefficient that falls as :math:`0.5^{|i-j|}` with the distance between their indices, and the channels' noise standard deviations range from 0.6 to 1.6 times a common level. .. GENERATED FROM PYTHON SOURCE LINES 146-178 .. code-block:: Python import torch import bartorch import bartorch.tools as bt from bartorch import apps COILS = 8 SIGMA = 0.01 # noise level, relative to the image's peak sensitivities = bt.coils(t=bt.grid(D=(SIZE, SIZE, 1)), n=COILS)[:, 0] sensitivities /= bartorch.rss(sensitivities, axes=(0,), keepdim=True) 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, :] .. GENERATED FROM PYTHON SOURCE LINES 179-183 The noise scan records the same channels with no signal. :func:`bartorch.tools.whiten` takes the data and the noise scan, estimates :math:`\Psi` as the sample covariance of the noise scan, and whitens the data with it. Applied to the noise scan itself, it returns white noise: .. GENERATED FROM PYTHON SOURCE LINES 183-209 .. code-block:: Python noise_scan = channel_noise((SIZE, SIZE))[:, None] # (coils, 1, y, x) whitened_scan = bt.whiten(noise_scan, noise_scan) .. image-sg:: /auto_examples/07-tours/images/sphx_glr_08-noise-prewhitening_001.png :alt: measured, whitened :srcset: /auto_examples/07-tours/images/sphx_glr_08-noise-prewhitening_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none off-diagonal covariance: 0.202 measured 0.000 after whitening .. GENERATED FROM PYTHON SOURCE LINES 210-225 The printed number is the mean magnitude of the off-diagonal covariance entries relative to the mean diagonal entry, which is zero for uncorrelated channels. In the figure, 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 --------------------------- You acquire every second phase encode (:math:`R = 2`), with 24 central lines kept as the ACS region. Each pseudo-replica adds a new noise realisation to the same noise-free k-space. You calibrate the sensitivities once per pipeline, from the first replica: from the channels as measured, and from the whitened channels. .. GENERATED FROM PYTHON SOURCE LINES 225-252 .. code-block:: Python CALIBRATION = 24 REPLICAS = 32 lines = torch.zeros(SIZE) lines[::2] = 1.0 centre = slice(SIZE // 2 - CALIBRATION // 2, SIZE // 2 + CALIBRATION // 2) lines[centre] = 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.""" noisy = noiseless + channel_noise((SIZE, SIZE)) return (noisy * 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 253-259 Pseudo-replicas --------------- Both pipelines run the same reconstruction, conjugate-gradient SENSE with a small Tikhonov weight and a fixed number of iterations. The only difference is whether you whiten the data first. .. GENERATED FROM PYTHON SOURCE LINES 259-268 .. code-block:: Python plain, prewhitened = [], [] for _ in range(REPLICAS): data = acquire() plain.append(apps.pics(data, maps_measured, l2=1e-3, maxiter=30)) data = bt.whiten(data, noise_scan) prewhitened.append(apps.pics(data, maps_whitened, l2=1e-3, maxiter=30)) plain, prewhitened = torch.stack(plain), torch.stack(prewhitened) .. GENERATED FROM PYTHON SOURCE LINES 269-272 The SNR of a voxel is the magnitude of the mean reconstruction over the standard deviation across replicas. You report both pipelines over the white matter, where the phantom is homogeneous. .. GENERATED FROM PYTHON SOURCE LINES 272-318 .. code-block:: Python def snr_map(replicas): return replicas.mean(0).abs() / replicas.std(0) gain = snr_map(prewhitened) / snr_map(plain) .. rst-class:: sphx-glr-horizontal * .. image-sg:: /auto_examples/07-tours/images/sphx_glr_08-noise-prewhitening_002.png :alt: SNR, as measured, SNR, prewhitened :srcset: /auto_examples/07-tours/images/sphx_glr_08-noise-prewhitening_002.png :class: sphx-glr-multi-img * .. image-sg:: /auto_examples/07-tours/images/sphx_glr_08-noise-prewhitening_003.png :alt: SNR ratio :srcset: /auto_examples/07-tours/images/sphx_glr_08-noise-prewhitening_003.png :class: sphx-glr-multi-img .. rst-class:: sphx-glr-script-out .. code-block:: none as measured white-matter SNR 38.1 (median 37.7) prewhitened white-matter SNR 55.8 (median 54.5) SNR ratio, prewhitened over as measured: 1.43 (median) .. GENERATED FROM PYTHON SOURCE LINES 319-347 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 (the bottom of the image). Both the quietest channels, 0 and 1, and the noisiest, 6 and 7, are most sensitive there. Without whitening, the least-squares fit gives them all the same weight; after whitening, each channel enters with the weight its noise level warrants, and that matters most where channels of very different noise levels overlap. How large the gain is depends on the array's noise correlation and on the spread of its channel noise levels, and what you measured here is for one simulated covariance. With uncorrelated channels of equal noise level it is one. As a spec --------- What this tour built, stated the way you would ask an agent for it: .. code-block:: text Measure the SNR gain of noise prewhitening for a Cartesian SENSE reconstruction with bartorch. Simulate an 8-channel acquisition with a known channel noise covariance and a noise-only scan; undersample by 2 along y with a 24-line ACS block. Calibrate one set of ESPIRiT maps (bt.ecalib, crop 0.8) from the data as measured and from the data whitened by bt.whiten(data, noise_scan). Reconstruct 32 pseudo-replicas with apps.pics (l2=1e-3, 30 iterations) in both pipelines, and report the SNR map (|mean| / std across replicas), its white-matter mean and median, and the median prewhitened / as-measured SNR ratio. .. GENERATED FROM PYTHON SOURCE LINES 349-365 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 4.922 seconds) .. _sphx_glr_download_auto_examples_07-tours_08-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: 08-noise-prewhitening.ipynb <08-noise-prewhitening.ipynb>` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: 08-noise-prewhitening.py <08-noise-prewhitening.py>` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: 08-noise-prewhitening.zip <08-noise-prewhitening.zip>` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_