the Creative Commons Attribution 4.0 License.

the Creative Commons Attribution 4.0 License.

# Nuclear spin noise tomography in three dimensions with iterative simultaneous algebraic reconstruction technique (SART) processing

### Stephan J. Ginthör

### Judith Schlagnitweit

### Matthias Bechmann

We report three-dimensional spin noise imaging (SNI) of nuclear spin density from spin noise data acquired by Faraday detection. Our approach substantially extends and improves the two-dimensional SNI method for excitation-less magnetic resonance tomography reported earlier (Müller and Jerschow, 2006). This proof of principle was achieved by taking advantage of the particular continuous nature of spin noise acquired in the presence of constant magnitude magnetic field gradients and recent advances in nuclear spin noise spectroscopy acquisition as well as novel processing techniques. In this type of projection–reconstruction-based spin noise imaging the trade-off between signal-to-noise ratio (or image contrast) and resolution can be adjusted a posteriori during processing of the original time-domain data by iterative image reconstruction in a unique way not possible in conventional rf-pulse-dependent magnetic resonance imaging (MRI). The 3D SNI is demonstrated as a proof of concept on a commercial 700 MHz high-resolution NMR spectrometer, using a 3D-printed polymeric phantom immersed in water.

- Article
(2029 KB) -
Supplement
(95 KB) - BibTeX
- EndNote

The phenomenon of nuclear spin noise, first predicted by Felix Bloch in 1946 (Bloch, 1946), can be ascribed to the incomplete cancellation of the fluctuating spin magnetic moments in a specimen. Owing to the extremely low amplitudes of these residual fluctuations of the bulk magnetic moment, very low-noise rf (radio-frequency) circuitry is required to separate nuclear spin noise signals from background noise (Müller et al., 2013; Ferrand et al., 2015; Pöschko et al., 2017). Therefore, experimental detection of nuclear spin noise succeeded only in 1985 (Sleator et al., 1985). Today readily available low-noise rf electronics enable one to observe nuclear spin noise in reasonable amounts of time, in particular if a cryogenically cooled probe circuit is used (Kovacs et al., 2005; Müller and Jerschow, 2006). In spite of low intrinsic sensitivity, nuclear spin noise-based spectroscopy and imaging techniques have an intriguing potential, mainly for three reasons: (1) the spin noise signal magnitude exceeds the thermal polarization-derived signal for very low numbers ($<\sim {\mathrm{10}}^{\mathrm{8}}$) of nuclear spins as it scales with the square root of the number of spins and not linearly. (2) The observation of undisturbed spin systems becomes possible. (3) Signals acquired in the absence of rf pulses are devoid of limitations imposed by pulse imperfections and bandwidth. The theoretical and technical aspects of nuclear spin noise detected by Faraday induction have been studied extensively in recent years (Marion and Desvaux, 2008; Nausner et al., 2009; Desvaux et al., 2009; Müller et al., 2013; Chandra et al., 2013; Ferrand et al., 2015; Pöschko et al., 2017; Ginthör et al., 2018). In the research we report here, a major sensitivity and image quality improvement in the case of spin noise imaging (SNI; Müller and Jerschow, 2006) is achieved by exploiting the tuning dependence of the spin noise line shape (Pöschko et al., 2014).

## 1.1 Magnetic resonance imaging

Spatial resolution in magnetic resonance imaging (MRI) relies on frequency encoding of the positions in magnetic field gradients (Kumar et al., 1975). Thus, spectra recorded in the presence of magnetic field gradients can be interpreted as the result of a forward radon transformation applied to the spin density function of the sample along the gradient direction (Deans, 2007).

Most of today's MRI approaches are based on the idea of Fourier imaging originally introduced by Richard R. Ernst (Kumar et al., 1975) and use rf pulses preceding excitation and evolution periods, in which phase encoding via field gradients is attained, followed by a two-dimensional (2D) Fourier transformation (FT; Wright, 1997; Brown et al., 2014). The resolution in the third dimension is commonly obtained by using slice selection, i.e., by applying rf pulses simultaneously with magnetic field gradients (Garroway et al., 1974). A multi-dimensional Fourier transform approach could, in principle, also work for noise data. However, it would suffer from even lower sensitivity due to transverse relaxation, in particular if no refocusing rf pulses are applied, which would counteract the goal of imaging without rf pulses. It should be noted here that when detecting noise from very low numbers of spins refocusing pulses are a viable option, in particular if indirect detection by optical means is used (Meriles et al., 2010).

## 1.2 Projection reconstruction

Alternatively, multi-dimensional MRI can be based on the principle of reconstruction from projected data in the direct (frequency or spatial) observation dimension (Chetih and Messali, 2015), which closely corresponds to the only approach available in radiation-based computer tomography (CT). In the case of MRI, projections are provided through acquisition of free induction decays, while different magnetic field gradients spanning the entire directional space are applied sequentially.

An intrinsic problem of the inverse radon transform is that it almost always produces artifacts if the projections contain substantial noise contributions or if the imaged distribution function is not smooth (Kabanikhin, 2008). In particular, the latter is a relevant restriction in most practical applications, as hard edges are very common features (e.g., bones in the human body). Therefore, alternatives to the inverse radon transform are required. The processing protocol presented here aims to minimize these artifacts while maintaining resolution limited by transverse relaxation and maximizing the spin noise to random noise contrast.

Different from spectroscopic applications of magnetic resonance, where the chemical shift is of main interest, for imaging purposes one often assumes that all spins inside the imaged object have indistinguishable chemical shifts. We are going to use that assumption here, being aware that methods to cope with imaging artifacts caused by non-uniform chemical shifts (e.g., water and fat) within a specimen exist (Dietrich et al., 2008) but may not be applicable straightforwardly in the case of spin noise detection.

In order to optimize the spin noise detection, we take advantage of the
progress in nuclear spin noise spectroscopy that has been achieved since the
introduction of 2D nuclear spin noise imaging (Müller and Jerschow,
2006). In order to obtain symmetrical line shapes and optimize receiver
sensitivity, the cryogenically cooled NMR probe is tuned to the SNTO (spin noise tuning optimum; Marion and Desvaux, 2008; Nausner et al., 2009; Pöschko et al., 2014). As residual static field gradients can lead to
spectral artifacts under conditions where radiation damping occurs (Pöschko et al., 2017), careful optimization of the basic magnetic field homogeneity (achieved by shimming) and sufficiently strong gradients, so
that ${T}_{\mathrm{2}}^{\ast -\mathrm{1}}$ exceeds the broadening caused by the radiation
damping rate ${\mathit{\lambda}}_{\mathrm{rd}}^{\mathrm{0}}$, are prerequisites for obtaining accurate
nuclear spin noise images. While the sensitivity of spin noise acquisition
depends in a highly non-linear way on the radiation damping rate *λ*_{rd} (McCoy and Ernst, 1989; Nausner et al., 2009; Bechmann and
Müller, 2017; Pöschko et al., 2017), it is proportional to the
radiation damping rate at equilibrium in this imaging regime.

Here, the only variables which depend on the instrumental setup are the
filling factor *η* and the probe quality factor *Q*, while the gyromagnetic ratio *γ* and *μ*_{0}, the permeability of the vacuum,
are immutable. Note that apart from maximizing the coupling between the
spins and the rf circuit by increasing *η* and *Q*, the SNR (signal-to-noise ratio) of spin-noise-detected experiments can be improved by reducing the noise from all other sources.

We demonstrate nuclear spin noise 3D tomography on the phantom shown and described in Fig. 1 using the acquisition and processing procedures outlined in Sect. 3.

The detected noise voltage was digitized quasi-continuously for each of 900 gradient directions uniformly distributed in 3D space. Each of the raw noise time-domain data blocks was divided into overlapping windows, which were Fourier transformed and the resulting power spectral data added up to yield individual projections for each gradient setting. Our newly developed iterative projection–reconstruction protocol combines projections obtained by Fourier transform of the time-domain noise data with different sliding window sizes (Desvaux et al., 2009) to obtain a 3D tomogram. This approach affords superior image quality with respect to resolution and contrast as compared to using a single fixed sliding window size only. A final reconstructed image of the phantom obtained from the projections of spin noise with different gradients by the optimized iterative reconstruction based on the simultaneous algebraic reconstruction technique (SART) introduced by Andersen and Kak (Andersen and Kak, 1984) can be seen in Fig. 2a and d.

In Fig. 2b and e we show the corresponding views, obtained using a fixed sliding window size of 1024 data points (corresponding to 103 ms).

This resolution is less than the maximum obtainable one based on the
experimentally determined water proton transverse relaxation time (*T*_{2}
= 380 ms) in a reference sample. But increasing the window size beyond 1024 data points reduces the SNR too much. The same raw data processed with
a small window size (128 data points) result in a correspondingly higher SNR but achieve an overly smoothed representation of the phantom as shown in
Fig. 2c and f. The window length of the smaller blocks was chosen
empirically by halving the length of the longer windows until the SNR was
acceptable while still being able to resolve courser details in the image.

In Fig. 2d–f images of cross sections of the phantom obtained by the different processing schemes are compared, demonstrating the flexibility of adjusting the contrast/resolution trade-off a posteriori from the same raw data. The high-resolution image (Fig. 2c and f) obtained by the conventional SART method might accurately represent the phantom, but the low SNR makes it difficult to draw a clear separation line between the phantom and the surrounding water in the cross section in Fig. 2f.

The low-resolution image in Fig. 2b and e, also obtained by SART, shows this boundary more clearly, but the resolution is too low to make out the correct shape of the phantom in the cross section in Fig. 2e. Only the image calculated with the iterative reconstruction method (Fig. 2a and d) shows a phantom with clear boundaries and a non-circular shape in the cross section in Fig. 2d, where the four “pockets” of water formed by the phantom can be resolved.

In Sect. 3 we describe the new processing procedure yielding the images in Fig. 2a and d, which, to our knowledge, is the most efficient way to obtain 3D images from spin noise data, currently. Most notably the continuous nature of spin noise time-domain data allows one to decide the resolution vs. contrast trade-off by reprocessing the raw data.

For completeness, we discuss the concept of indirect detection of the
spatial dimension in Fourier imaging with spin noise. Even though spin noise
has random phase, one can devise a way to encode spatial information in the
indirect dimension similarly to the way used for 2D spin-noise-detected spectroscopy (Chandra et al., 2013; Ginthör et al., 2018). This would
require a location-encoding gradient sandwiched between two acquisition
blocks and incremented in the usual way to yield the indirect *k*-space
dimension. This phase-encoding gradient modulates the relative phase of the
signal, which could be resolved by cross-correlation of two subsequent
acquisition blocks. However, this theoretical scheme suffers from excessive
relaxation losses occurring during this gradient and the two acquisition
periods involved. We have so far not succeeded in obtaining sufficient image
contrast using this scheme. Notably, with this indirect time-domain (*k*-space) approach one cannot take advantage of the sliding window
processing and a posteriori optimization as this acquisition scheme only affords discrete short data blocks.

The experiments were carried out on a high-resolution NMR spectrometer
equipped with a Bruker Avance III console connected to a Bruker Ascend 700
MHz magnet and a TCI cryoprobe (manufactured in 2011). The phantom is a
3D-printed helix made of PLA (polylactic acid) fitting tightly inside a
standard 5 mm NMR tube (Wilmad 535-PP). The tube is filled with a mixture of
H_{2}O:D_{2}O (9 : 1) and the filling height made equal to the height of
the phantom (50 mm). The rf probe is tuned to the SNTO (spin noise tuning optimum; Marion and Desvaux, 2008; Nausner et al., 2009; Pöschko et
al., 2014), and 3D shimming (using Bruker's TopShim) is performed on a sample (henceforth referred to as “shim sample”) which is prepared in the same
type of 5 mm tube (Wilmad 535-PP) as the phantom but filled with the same
mixture of H_{2}O:D_{2}O (9 : 1) to the same filling height without the
phantom. The radiation damping rate of this phantom in the probe used was
determined as ${\mathit{\lambda}}_{\mathrm{rd}}^{\mathrm{0}}=\mathrm{611.4}$ rad s^{−1} under the
conditions of the imaging experiment, while the transverse relaxation rate in the absence of radiation damping was $\mathrm{1}/{T}_{\mathrm{2}}=\mathrm{2.632}$ s^{−1}.

A one-dimensional (1D) ^{1}H spin noise spectrum was acquired to verify
the correctness of the setup and the shim. For this imaging experiment the
magnetic field gradients generated by the *X*, *Y*, and *Z* field correction
(shim) coils are used as the imaging gradients, after calibration for the
purpose of imaging. The gradient amplitude is measured by the broadening of
the peak in a spectrum of the test sample. For each projection direction the
magnitude of the applied magnetic field gradient must be the same,
independent of the direction set according to the *φ* and *θ*
values. (The definitions of the coordinate system and angles are given in
the Supplement.) This is achieved by calculating the
individual shim values via the trigonometric laws. A few compound gradient
settings (involving multiple individual gradients) are checked to verify
orthogonality of the Cartesian field gradient components that are used to calculate the gradient directions. After setup and calibration, the actual
phantom is inserted into the magnet without further shimming. A script
controls the setting of the gradients via the digital-to-analog converters of the shim system and initiates the acquisition. For each angle *φ*
and *θ* 30 projections were recorded, yielding 900 projections in
total at a gradient amplitude of 78 mT m^{−1}. The values are sampled uniformly
per angle in the range of 0 to 180^{∘}. The projections were
recorded with the following settings: for each projection angle two noise blocks were recorded (time-domain points: 1024*k*, spectral width: 5 kHz, spectral center: 4.670 ppm). The individual projections are recorded like standard 1D spin-noise-detected spectra (Müller and Jerschow, 2006; Nausner et al., 2009). The acquisition sequence is illustrated and compared
to conventional pulsed NMR in Fig. 3.

As in a conventional single-pulse NMR experiment (Fig. 3a), the direction of the projection is determined by the magnetic field gradient active during
acquisition. Gradients for different projection directions are set to the
same magnitude and thus only differ in their direction. The angles are laid
out in the following way: *V* angles for *φ* (angle between *X* axis and *X**Y* projection of the direction vector) are chosen uniformly. For each angle *φ* we again sample *U* angles for *θ* (angle between the *Z* axis and
direction vector) uniformly. A reference coordinate system is shown in the
Supplement. The *θ* angles have the same set of values
for all *φ* angles.

Due to the random phase of spin noise, direct accumulation of raw-phase sensitive data in the time domain (as is usually done for pulsed experiments) leads to signal cancellation. Instead, all acquired noise blocks (the spin noise equivalent of free induction decays are stored and Fourier transformed individually. After calculating the power or magnitude spectra, their addition yields the final projection (McCoy and Ernst, 1989; Müller and Jerschow, 2006; Nausner et al., 2009).

Longitudinal relaxation not being an issue in the absence of rf pulses, any recycling delay can be omitted. That allows one to record all blocks for one gradient in a single very long acquisition block, as indicated in Fig. 3d. During processing the data are split into adjustable smaller blocks which can be processed as described above for individually acquired ones. The acquisition of a single large noise block also allows one to take advantage of “sliding window processing”, which was introduced by Desvaux et al. (2009). By slicing a long acquisition block recorded into overlapping sub-blocks, one gains the option to compromise a posteriori (i.e., after acquisition of the raw data) between resolution and sensitivity, by adjusting the size of the smaller sub-blocks. The optimum overlap of the blocks or “selection windows” has been shown to be approximately one-seventh of the windows' size (Desvaux et al., 2009). Due to this overlap 7 times more blocks are used in the summation process, which improves the SNR by a factor of about √2, owing to the different accumulation behavior of correlated (spin) noise and uncorrelated (instrument and circuit) noise. After the sliding window splitting, the individual blocks are processed in the same way as individually acquired ones. The processing schemes from time-domain noise data to the final projection with and without the sliding window processing are compared in Fig. 4.

The processing was achieved by a custom Python script (available in the Supplement) using the numpy library (Oliphant, 2006) for general numerical calculations and for the FFT routine as well as scikit-image (Walt et al., 2014) for implementing the 2D SART reconstruction algorithm (see below). The iterative reconstruction method was used with two different resolutions. The lower-resolution images were processed with a sliding window length of 128 and an overlap of 110 data points, the higher-resolution ones with a length of 1024 and an overlap of 878 data points. The resulting three-dimensional (3D) images were visualized with ParaView (Ahrens et al., 2005).

We found the SART introduced by Andersen and Kak (Andersen and Kak, 1984) to be most suitable
as an alternative to the inverse radon transform for the spin noise imaging data. This algorithm sets up a set of linear equations (see Eq. 2) that
describe the dependency of the projected values on the distribution function
*D*(** r**) of the original sample (Andersen and Kak, 1984).

Integration occurs over the area *A* and ℘_{β,i} denotes the data point *i* in the projection ℘ at the angle *β*. The vector
** r** encodes the position of a point inside the sample. In
our case,

*A*is an area along the projection direction through the sample with the width of one data point. No exact solution of this system is possible because it is usually underdetermined (owing to a finite number of projection values ℘ compared to a continuous distribution function). Therefore

*D*(

**), the solution of this system of linear equations, is approximated by an appropriate method, e.g., the Kaczmarz method (Kaczmarz, 1937). For two dimensions the equations can be generated in the following way: given**

*r**V*data points for each projection, an image matrix ℑ of

*V*×

*V*pixels is created, representing the reconstructed image. The image matrix ℑ is a discrete approximation of the continuous distribution function

*D*(

**). There are multiple ways to formulate the forward projection process, but here a bilinear model is chosen in Eq. (3), as it is the basis for the SART reconstruction method (Andersen and Kak, 1984).**

*r*
℘_{β,i} denotes the data point *i* in the projection ℘ at the
angle *β*. The image matrix ℑ is sampled *j* times along the
projection direction for each point *i* in the projection ℘_{β}.

The resulting value for the sampling location is the weighted sum of the
four closest pixels (index *k*) in the image matrix, where *w* denotes the weighting factor for the individual pixel. This also shows the
underdetermination of the system. There are *N*^{2} variables (pixels in the image matrix) but only *x*×*V* equations (*x* represents the number of projection angles). The projections can be conceptually arranged around the matrix
ℑ at their respective projection angles. The reconstruction
works in the following way: from each projection ℘_{β} (corresponding to angle *β*) an updated image is calculated according
to Eq. (4) (Andersen and Kak, 1984).

Here *u* denotes the update function, ℘_{β} the projection at angle *β* and *r* the image reconstruction relaxation factor (with no relation to spin relaxation). *z* indicates the evolution of the image matrix
after an update has been applied. All projections ℘ are processed
sequentially in an order that maximizes the difference in angles between two
successive projections. The reconstruction is complete when all the
projections ℘ have been considered. The goal is to approximate *D*(** r**) with
ℑ. The update function

*u*performs three steps for each projection pixel ℘

_{β,i}. First, from each projection pixel ℘

_{β,i}a ray is cast into ℑ along the projection direction.

*S*intermediate values of the matrix ℑ are taken on that ray at equidistant locations ℒ

_{j}. At each location the values of the four closest matrix ℑ pixels, weighted by their distance to the sampling point, are summed to give the value assigned to the sampling location ℒ

_{j}. In Fig. 5 this process is illustrated.

The sum of all sampling points along the projection ray yields a new value
for each projection pixel ℘_{β,i}, named ℘_{calc β,i} according to Eq. (5).

The index variable *j* spans all sampling locations ℒ_{j} along the
ray; *k* is used to indicate the four surrounding pixels of a single sampling location ℒ_{j}; *w* is the weight used to calculate the
contribution of the neighboring pixel in question. This is repeated for all projection pixels ℘_{β,i} of a given projection ℘_{β}. The difference in value of the actual projection pixel ℘_{β,i} and the calculated projection pixel ℘_{calc β,i} is given
in Eq. (6).

This difference is then distributed back to the individual sampling
locations and in turn to their surrounding pixels, weighted by the
respective distance from the sampling location to the neighboring pixel, in a new previously empty image matrix ℑ_{update} as seen in
Eq. (7).

This matrix ℑ_{update} is the result of the update function
*u* and used to update the image matrix in Eq. (4).

The algorithm can be initiated with a first guess of the image matrix
ℑ_{z}=0. Starting with a reasonable guess as the initial
matrix instead of a zero matrix improves the reconstruction, provided the
guess is not too far off; in our case we use a low-resolution image.

The reconstruction is iterated a pre-determined number of times (in our case
two) with the same projections but always using the result from the previous step as the next initial guess. In each step the weight of narrow contributions increases (e.g., better-defined edges or increased contrast), while at the same time the noise level is rising. The number of iterations
depends on the quality of the original projections and is, for actual data,
best determined by visual inspection of the image improvement during the
iteration process. The relaxation factor *r* in Eq. (4) was set to 0.05 for all
reconstructions.

The projections are grouped by the values of the angle *φ*, yielding
*V* groups with *U* projections each. For each of these groups a SART image
reconstruction (Andersen and Kak, 1984) is performed, yielding *V* 2D images.
The rows (with respect to the *Z* axis) with the same index from all the images are grouped again: every first row in a group, every second in a group, and so on. Again, for each of these groups an image reconstruction is performed. This results in the final 3D image of the phantom. Figure 6
illustrates this procedure.

In conjunction with the SART image reconstruction the flexibility of large noise block acquisition and the sliding window processing turns out to be particularly advantageous. Shorter windows yield a higher number of blocks with fewer data points and hence lower resolution. The summation over a larger number of noise blocks, however, improves the SNR of the resulting projection. Larger windows have the opposite effect and improve the resolution while losing SNR. This fact is exemplified in Fig. 7. The initial image reconstruction (via SART; see above) is carried out with a zero seed matrix using projections computed from short sliding windows (e.g., 128 complex data points).

The reconstructed image from these low-resolution projections then serves as the starting image matrix for the second reconstruction using projections calculated with a longer sliding window. A schematic overview of the procedure can be seen in Fig. 8.

For additional clarity, we summarize the iterative reconstruction process:
the recorded noise blocks corresponding to different projection angles are processed with different increasing sliding window sizes. This procedure
yields multiple sets of projections that differ in their resolution and SNR. Then a first 2D SART image reconstruction step is done separately for
each set of 1D projections corresponding to a particular angle *φ*.
The first reconstruction uses no initial guess (i.e., a zero image matrix) and is reconstructed using the set of projections with the lowest resolution
and correspondingly highest SNR. The obtained 2D image is used as an initial guess for another reconstruction, which is improved by the next set of
projections with higher resolution but lower SNR. Continuing this series,
the obtained 2D images now form a sequence ranging from low-resolution, high-SNR images to the inverse high-resolution, low-SNR case. In order to use the lower-resolution images in the higher-resolution reconstructions, the image matrices need to be interpolated to match the required number of data
points. The sliding window sizes are doubled on each incrementation such that the resolution always increases by a factor of 2. This procedure is
used analogously in the second reconstruction step to obtain the final 3D
image. We note that the SART algorithm can also be applied to all three dimensions in a single step without intermediate 2D data (Mueller et al.
1999), but this has not been done here due to higher computational demands
and the easy availability of well-tested open-source 2D SART variants.

We have improved spin-noise-detected NMR imaging and extended it to three dimensions from the original 2D SNI technique published in 2006 (Müller and Jerschow, 2006). The unique properties of spin noise, in particular its average power being constant, while phase memory is lost with the usual time constant ${T}_{\mathrm{2}}^{\ast}$, and the absence of a defined starting point in time together with spin noise tuning optimization (Nausner, 2009; Pöschko, 2014) make it possible to use a quasi-continuous acquisition technique. Thus, one can obtain raw noise data of arbitrary duration while applying constant magnitude magnetic field gradients in different directions. The data can be processed in a unique way using dynamically adapted sliding windows to define fractionally overlapping data blocks, which are Fourier transformed and co-added after power spectrum computation.

Based on these fundamental steps, we have introduced a new kind of SART-based iterative image reconstruction technique, which yields 3D images that are superior in visual quality, improving SNR and resolution at the same time without introduction of artifacts. In this type of spin noise imaging the well-known trade-off between SNR (or image contrast) and resolution can be adjusted a posteriori during processing of the same original data by iterative image reconstruction, which is not applicable in conventional rf-pulse-dependent MRI.

The raw time domain spin noise data used to generate the 3D images in this paper have been uploaded to Zenodo and can be accessed by the https://doi.org/10.5281/zenodo.3967145 (Ginthör, 2020).

The supplement related to this article is available online at: https://doi.org/10.5194/mr-1-165-2020-supplement.

NM planned and supervised the research; SJG performed all experiments, programmed the processing algorithms and processed the data; JS devised the time-domain spin noise imaging experiment and contributed to the manuscript; MB and SJG programmed the spectrometer; NM and SJG wrote the manuscript.

The authors declare that they have no conflict of interest.

This work has been supported by the European Union through the ERDF INTERREG IV (RU2-EU-124/100–2010) program (ETC Austria-Czech Republic 2007–2013, project M00146, “RERI-uasb” for Norbert Müller).

This research has been supported by the European Regional Development Fund (grant no. RU2-EU-124/100–2010/M00146).

This paper was edited by Matthias Ernst and reviewed by two anonymous referees.

Ahrens, J., Geveci, B. and Law, C.: ParaView: An End-User Tool for Large Data Visualization, Elsevier, https://doi.org/10.1016/B978-012387582-2/50038-1, 2005.

Andersen, A. H. and Kak, A. C.: Simultaneous Algebraic Reconstruction Technique (SART): A superior implementation of the ART algorithm, Ultrasonic Imaging, 6, 81–94, https://doi.org/10.1016/0161-7346(84)90008-7, 1984.

Bechmann, M. and Müller, N.: Nonlinear Effects in NMR, Ann. Rep. NMR Spectrosc., 92, 199–226, https://doi.org/10.1016/bs.arnmr.2017.04.005, 2017.

Bloch, F.: Nuclear Induction, Phys. Rev., 70, 460–474, https://doi.org/10.1103/PhysRev.70.460, 1946.

Brown, R. W., Cheng, Y.-C. N., Haacke, E. M., Thompson, M. R., and Venkatesan, R., Eds.: Magnetic Resonance Imaging: Physical Principles and Sequence Design, John Wiley & Sons Ltd, Chichester, UK, 2014.

Chandra, K., Schlagnitweit, J., Wohlschlager, C., Jerschow, A., and Müller, N.: Spin-Noise-Detected Two-Dimensional Fourier-Transform NMR Spectroscopy, J. Phys. Chem. Lett., 4, 3853–3856, https://doi.org/10.1021/jz402100g, 2013.

Chetih, N. and Messali, Z.: Tomographic image reconstruction using filtered back projection (FBP) and algebraic reconstruction technique (ART), in: 2015 3rd International Conference on Control, Engineering Information Technology (CEIT), 25–27 May 2015, Tlemcen, Algeria, 1–6, 2015.

Deans, S. R.: The Radon Transform and Some of Its Applications, Courier Corporation, New York, USA, 2007.

Desvaux, H., Marion, D. J. Y., Huber, G., and Berthault, P.: Nuclear Spin-Noise Spectra of Hyperpolarized Systems, Angew. Chem. Int. Edit., 48, 4341–4343, https://doi.org/10.1002/anie.200901034, 2009.

Dietrich, O., Reiser, M. F., and Schoenberg, S. O.: Artifacts in 3-T MRI: Physical background and reduction strategies, Eur. J. Radiol., 65, 29–35, https://doi.org/10.1016/j.ejrad.2007.11.005, 2008.

Ferrand, G., Huber, G., Luong, M., and Desvaux, H.: Nuclear spin noise in NMR revisited, J. Chem. Phys., 143, 094201, https://doi.org/10.1063/1.4929783, 2015.

Garroway, A. N., Grannell, P. K., and Mansfield, P.: Image formation in NMR by a selective irradiative process, J. Phys. C: Solid State Phys., 7, L457–L462, https://doi.org/10.1088/0022-3719/7/24/006, 1974.

Ginthör, S. J.: Spin Noise Data, Zenodo, https://doi.org/10.5281/zenodo.3967145, 2020.

Ginthör, S. J., Chandra, K., Bechmann, M., Rodin, V. V., and Müller, N.: Spin-Noise-Detected Two-Dimensional Nuclear Magnetic Resonance at Triple Sensitivity, ChemPhysChem, 19, 907–912, https://doi.org/10.1002/cphc.201800008, 2018.

Kabanikhin, S. I.: Definitions and examples of inverse and ill-posed problems, J. Inverse. Ill-Pose. P., 16, 317–357, https://doi.org/10.1515/JIIP.2008.019, 2008.

Kaczmarz, S.: Angenäherte Auflösung von Systemen linearer Gleichungen, Bulletin International de l'Académie Polonaise des Sciences et des Lettres, Krakow, Poland, 355–357, 1937.

Kovacs, H., Moskau, D., and Spraul, M.: Cryogenically cooled probes—a leap in NMR technology, Prog. Nucl. Magn. Reson. Spectrosc., 46, 131–155, https://doi.org/10.1016/j.pnmrs.2005.03.001, 2005.

Kumar, A., Welti, D., and Ernst, R. R.: NMR Fourier zeugmatography, J. Magn. Reson., 18, 69–83, https://doi.org/10.1016/0022-2364(75)90224-3, 1975.

Marion, D. J.-Y. and Desvaux, H.: An alternative tuning approach to enhance NMR signals, J. Magn. Reson., 193, 153–157, https://doi.org/10.1016/j.jmr.2008.04.026, 2008.

McCoy, M. A. and Ernst, R. R.: Nuclear spin noise at room temperature, Chem. Phys. Lett., 159, 587–593, https://doi.org/10.1016/0009-2614(89)87537-2, 1989.

Meriles, C. A., Jiang, L., Goldstein, G., Hodges, J. S., Maze, J., Lukin, M. D., and Cappellaro, P.: Imaging mesoscopic nuclear spin noise with a diamond magnetometer, J. Chem. Phys., 133, 124105, https://doi.org/10.1063/1.3483676, 2010.

Mueller, K., Yagel, R., and Wheller, J. J.: Fast implementations of algebraic methods for three-dimensional reconstruction from cone-beam data, IEEE Transactions on Medical Imaging, 18, 538–548, https://doi.org/10.1109/42.781018, 1999.

Müller, N. and Jerschow, A.: Nuclear spin noise imaging, P. Natl. Acad. Sci. USA, 103, 6790–6792, https://doi.org/10.1073/pnas.0601743103, 2006.

Müller, N., Jerschow, A., and Schlagnitweit, J.: Nuclear Spin Noise, eMagRes, 2, 237–243, https://doi.org/10.1002/9780470034590.emrstm1314, 2013.

Nausner, M., Schlagnitweit, J., Smrečki, V., Yang, X., Jerschow, A., and Müller, N.: Non-linearity and frequency shifts of nuclear magnetic spin-noise, J. Magn. Reson., 198, 73–79, https://doi.org/10.1016/j.jmr.2009.01.019, 2009.

Oliphant, T. E.: A guide to NumPy, Trelgol Publishing, Austin, Texas, USA, 2006.

Pöschko, M. T., Schlagnitweit, J., Huber, G., Nausner, M., Horničáková, M., Desvaux, H., and Müller, N.: On the Tuning of High-Resolution NMR Probes, ChemPhysChem, 15, 3639–3645, https://doi.org/10.1002/cphc.201402236, 2014.

Pöschko, M. T., Rodin, V. V., Schlagnitweit, J., Müller, N., and Desvaux, H.: Nonlinear detection of secondary isotopic chemical shifts in NMR through spin noise, Nat. Commun., 8, 13914, https://doi.org/10.1038/ncomms13914, 2017.

Sleator, T., Hahn, E. L., Hilbert, C., and Clarke, J.: Nuclear-spin noise, Phys. Rev. Lett., 55, 1742–1745, https://doi.org/10.1103/PhysRevLett.55.1742, 1985.

Walt, S. van der, Schönberger, J. L., Nunez-Iglesias, J., Boulogne, F., Warner, J. D., Yager, N., Gouillart, E., and Yu, T.: scikit-image: image processing in Python, PeerJ, 2, e453, https://doi.org/10.7717/peerj.453, 2014.

Wright, G. A.: Magnetic resonance imaging, IEEE Signal Process. Mag., 14, 56–66, https://doi.org/10.1109/79.560324, 1997.