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.

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 (

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).

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

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.

Phantom used for the spin noise imaging experiments.

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.

Comparison of reconstructed

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 (

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

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

A one-dimensional (1D)

Basic acquisition sequences.

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:

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

Comparison of sliding window processing schemes for spin noise
data without

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

The resulting value for the sampling location is the weighted sum of the
four closest pixels (index

Scheme of our implementation of the SART algorithm (Andersen and
Kak, 1984). A ray is cast from every data point of the projection

The sum of all sampling points along the projection ray yields a new value
for each projection pixel

The algorithm can be initiated with a first guess of the image matrix

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

The projections are grouped by the values of the angle

3D SART reconstruction process. The 3D SART reconstruction
algorithm is realized as a two-step process via two successive
2D SART reconstructions.

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).

Qualitative comparison of resolution and signal-to-noise ratio for
different sliding window sizes using a 1D image of the phantom. Panel

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.

Iterative 3D SART reconstruction. The first 3D SART reconstruction uses no initial guess and reconstructs a low-resolution image with a high signal-to-noise ratio. This intermediate result is used as the initial guess for the second 3D SART reconstruction, yielding the final reconstructed image. Each 3D SART reconstruction is repeated multiple times with the same projections.

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

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

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

The supplement related to this article is available online at:

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.