Simultaneous imaging of hard and soft biological tissues in a low-field dental MRI scanner.

Magnetic Resonance Imaging (MRI) of hard biological tissues is challenging due to the fleeting lifetime and low strength of their response to resonant stimuli, especially at low magnetic fields. Consequently, the impact of MRI on some medical applications, such as dentistry, continues to be limited. Here, we present three-dimensional reconstructions of ex-vivo human teeth, as well as a rabbit head and part of a cow femur, all obtained at a field strength of 260 mT. These images are the first featuring soft and hard tissues simultaneously at sub-Tesla fields, and they have been acquired in a home-made, special-purpose, pre-medical MRI scanner designed with the goal of demonstrating dental imaging at low field settings. We encode spatial information with two pulse sequences: Pointwise-Encoding Time reduction with Radial Acquisition and a new sequence we have called Double Radial Non-Stop Spin Echo, which we find to perform better than the former. For image reconstruction we employ Algebraic Reconstruction Techniques (ART) as well as standard Fourier methods. An analysis of the resulting images shows that ART reconstructions exhibit a higher signal-to-noise ratio with a more homogeneous noise distribution.

Introduction

Magnetic Resonance Imaging (MRI[1]) plays an indispensable role in healthcare. In particular, its performance is unrivaled for soft tissues, being the only known technique capable of in-vivo imaging of deep biological tissues with high spatial resolution and tissue contrast while avoiding harmful ionizing radiation[2]. Despite MRI’s unquestionable success, imaging hard tissues (such as bone, tendons, dentin or enamel) remains problematic[3]. This is due to the fleeting lifetime and low strength of signals emitted by solids, as opposed to the case of soft and non-solid tissues[4]. In the latter, the deleterious dipole-dipole magnetic interaction between neighboring spins averages out much faster than any significant timescale in imaging processes. This results in strong, long-lived signals which are routinely exploited for the high quality images typical for MRI. In hard tissues, however, spins are to a good approximation fixed with respect to one another, and the dipole-dipole pattern vanishes exclusively at the so-called “magic angle”[5]. This means that every spin is subject to a noisy environment created by the surrounding spins. Since this contribution does not average out, the MR signal coherence of hard biological tissues is typically lost in hundreds of microseconds[6,7].

The above background justifies the rather meager penetration of MRI in dentistry[8]. Instead, odontologists heavily use X-ray based technologies[9]. However, these come together with a number of detrimental aspects, such as the use of ionizing radiation, low soft tissue contrast, and unreliable revelation of pulp diseases[10] or tooth cracks and fractures[11,12]. All these inconveniences can be overcome with MRI[12,13], albeit with its own particular challenges. Aside from the aforementioned technical difficulties in detecting hard biological tissues, MRI scanners are typically expensive to acquire, site, operate and maintain, mostly due to the high magnetic fields involved[14]. Low-field systems constitute a promising inexpensive alternative to standard MRI setups[14,15], but the reduced signal-to-noise ratio (SNR) can easily lead to long acquisition times incompatible with clinical conditions.

In this article we present “DentMRI - Gen I” (Fig. 1), a home-made special-purpose MRI scanner designed with the goal of demonstrating dental imaging at low field settings. We demonstrate the capability of “DentMRI - Gen I” to produce high quality combined images of soft and hard biological tissues at  260 mT. We use two pulse sequences: Pointwise-Encoding Time reduction with Radial Acquisition (PETRA[16,17]); and Double Radial Non-Stop Spin Echo (DRaNSSE), which we have devised to address sampling and contrast limitations we encounter with PETRA. These constitute, to our best knowledge, the first images featuring soft and hard tissues simultaneously at sub-Tesla fields. We reconstruct images of a rabbit head, bare human teeth and a cow femur with PETRA and DRaNSSE, using both Algebraic Reconstruction Techniques (ART[18-20]) and traditional Fourier transformation. We find that DRaNSSE can perform significantly better than PETRA and that ART results in higher quality reconstructions than Fourier spectral analysis with our settings.

Figure 1

(a) (Left) Photograph of the main magnet installed on the support structure. (Right) Photographs of “DentMRI - Gen I”, showing a general overview of the magnet and the components used for building the gradient and radio-frequency systems. (b) Sketch of the “DentMRI - Gen I” scanner.

Results

“DentMRI - Gen I” (Fig. 1a) is our first-generation MRI dental scanner, designed and built to demonstrate hard tissue imaging techniques in low magnetic field settings. Operating at a field of provided by a permanent magnet, all parts are standard commodities, 3D-printable, or inexpensive to machine. The total cost of all components adds up to €, a small fraction of the price of high field scanners (based on super-conducting magnets) previously used for dental imaging[8,12,13,16,21].

“Gen I” is meant for technology-demonstration purposes rather than for immediate clinical validation. We have therefore restricted our use to ex-vivo image acquisitions of samples accommodated in a cylindrical field of view of height  100 mm and diameter  45 mm, even if the gap between magnet poles is significantly larger. We also enclose the samples inside a Faraday cage to isolate our single radio-frequency (rf) coil from the magnetic gradient coils and electromagnetic interference noise present in the laboratory, which is otherwise unshielded. Further details can be found in the “Methods” section.

In the following, we will first provide images in which hard and soft biological tissues are simultaneously visible, we will then show 3D reconstructions of bare human teeth, and we will end the section with quantitative comparisons between the performance of PETRA and DRaNSSE pulse sequences, and between Fourier and ART-based mathematical reconstructions (See “Methods” section for further details). Table 1 contains the parameters used for all images in this section.

Table 1

Image acquisition parameters. “NA” stands for “not applicable”.

Image	Sequence	Flip angle (^∘)	Pulse time (us)	FOV (mm^3)	Pixel size (mm)	Dead time (us) / Acquisition time (us)	Bandwidth (kHz)	TE (us)	TR (ms)	Radial spokes	Single points	Scans	Scan time (min)
I. Fig. 2	PETRA	58	10	46 × 54 × 42	0.5	90 / 710	67.5	NA	10	4098	848	75	61
II. Fig. 2	PETRA	58	10	46 × 54 × 42	0.5	1,000 / 4,000	10.8	NA	10	4098	3904	25	33
III. Fig. 3	PETRA	67	9	40 × 40 × 40	0.5	100 /650	53	NA	15	2546	912	70	65
IV. Fig. 4	PETRA	90	9.1	46 × 48 × 32	1	85 / 915	24	NA	50	4896	40	7	29
V. Fig. 4	PETRA	90	9.1	46 × 48 × 32	1	1,000 / 2,000	8	NA	50	4896	1528	7	37
VI. Fig. 4	DRaNSSE	90	9.1	46 × 48 × 32	1	NA / 1,000	24	60	50	4896	NA	16	65
VII. Fig. 4	DRaNSSE	90	9.1	46 × 48 × 32	1	NA / 1,000	24	10,000	50	4896	NA	16	65
VIII. Fig. 5	DRaNSSE	90	10	44 × 52 × 42	1	NA / 1,000	26	60	50	1426	NA	26	31
IX. Fig. 5	DRaNSSE	90	10	44 × 52 × 42	1	NA / 1,000	26	10,000	50	1426	NA	26	31
X. Fig. 5	PETRA	90	10	44 × 52 × 42	1	90 / 910	5.2	NA	50	1426	64	12	15
XI. Fig. 5	PETRA	90	10	44 × 52 × 42	1	1,000 / 4,000	26	NA	50	1426	496	9	15

Combined tissue imaging: rabbit head

For initial demonstration purposes, we first show ex-vivo images of a rabbit head, which was soft boiled in tap water to delay tissue deterioration (Fig. 2a). Figure 2c contains selected slices from the full 3D reconstruction employing a PETRA sequence. The field of view is  mm, with an isotropic resolution of 0.5 mm. We excite the sample with a hard rf pulse of 10 s for a flip angle of  58 degrees, and the repetition time is TR = 10 ms. These parameters were used in two independent acquisitions: one with a short dead time (90 s, limited by ring-down in our setup) to read in the combined signal from hard and soft tissues (Fig. 2c top); and one with a long dead time (1 ms) to remove the short-lived contribution of teeth and skull tissues (Fig. 2c middle). Each radial acquisition lasted 710 s (4 ms) in the combined (soft) tissue scan, with a total of 4098 (4098) radial spokes, corresponding to an undersampling factor of 7 (7), and 848 (3904) single points to fill the k-space gap, requiring 75 (25) averages and a total scan time of 61 (33) min for the images in Fig. 2c. Both images are reconstructed with ART using and 2 iterations in Eq. (6), and denoised with Block-Matched Filters[22]. Reconstruction times take in a low-end GPU (see “Methods” ). Finally, we subtract one of the above scans from the other to produce the bottom image in Fig. 2c, where only hard tissues are highlighted. Even this basic post-processing is enough to identify the upper (A), bottom (B), and inner (C) teeth, as well as the rabbit jaw (D).

Figure 2

(a) Picture of the scanned rabbit head. (b) Picture of a rabbit skull. (c) Top: Single slices for 90 s dead time with PETRA; middle: the same slices for 1 ms dead time; bottom: difference between the above images. Further details can be found in the main text.

Ultra-short T2 tissue imaging: bare human teeth

The results shown in Fig. 2 demonstrate the capabilities of “DentMRI - Gen I” for imaging of soft and hard tissues simultaneously at sub-Tesla fields. Since human dental structures differ significantly from those in rabbits, below we demonstrate our scanner’s performance for a sample consisting of four human teeth embedded into a piece of ham emulating the gum (Fig. 3b).

Figure 3

(a) 2 dimensional slices of four human teeth embedded in a piece of pork ham (PETRA). (b) Photograph of the sample.

Figure 3a shows 2 dimensional slices obtained from a PETRA scan and applying ART. We acquired the image with a field of view of  mm and a nominal isotropic resolution of 0.5 mm. To excite the sample we used a hard pulse of 9 to produce a flip angle of  67 degrees. The dead time was set to 100 s and the acquisition time was 650 s with a repetition time of 15 ms. We acquired a total of 2546 radial spokes, corresponding to an undersampling factor of 8, and 912 single points. To increase the SNR, we averaged over 70 scans with a total scan time of  65 min. This image is reconstructed with ART using 2 iterations and in Eq. (6), and denoised with Block-Matched Filters[22]. Reconstruction takes in a low-end GPU (see “Methods” ). In Fig. 3a we identify the tooth and ham, visible in the photograph in Fig. 3b. The absence of pulp and nerve means there is no bright structure inside the dentin.

All in all, Fig. 3 demonstrates that human teeth can be imaged with high resolution at low magnetic fields, even if our scan times so far are hardly compatible with clinical use. We examine possible approaches to overcome this limitation in the Discussion section.

PETRA and DRaNSSE

As shown in Figs. 2 and 3, tissue contrast with weighting with PETRA is possible by image subtraction. This is a lengthy procedure, however, which can be significantly improved with DRaNSSE (Fig. 7b). Here we compare one approach against the other.

Figure 7

Sequence diagram for a single repetition of (a) PETRA and (b) DRaNSSE.

For these tests we use a piece of a cow femur (Fig. 4a), which is mostly composed of only two tissues (bone, with  ms, and marrow, with  ms) and therefore facilitates discrimination and image analysis. In Fig. 4b, images labeled as “PETRA long (short) ” correspond to an individual PETRA scan with long (short) dead times, where marrow (and bone) appear visible. For DRaNSSE we run a single scan, which can be used to reconstruct marrow (and bone) from the second (first) acquisition in the sequence, corresponding to images with the label “DRaNSSE SE2 (SE1)”.

Figure 4

(a) Photograph of cow bone sample; (b) raw image slices from 3 dimensional acquisitions with PETRA and DRaNSSE, reconstructed with ART and FT; (c) Signal to noise ratio along the red dotted lines in (b). White dashed lines highlight differences between ART and FT reconstructions. White arrows point at a small lump of soft tissue attached to the external surface of the bone.

Importantly, the total scan times for this study are kept the same, i.e. the sum of both PETRA scan durations is very close to the single DRaNSSE scan time ( 65 min). Also common to both sequences are: a flip angle of  90 degrees, a repetition time of 50 ms, a field of view of  mm, an isotropic voxel resolution of 1 mm, a sampling rate of 24 kHz (8 kHz for long PETRA), and a total of 4896 radial spokes in k-space, which is fully sampled in this case. For the short (long) PETRA acquisition we set the dead time to 85 s (1 ms), the radial acquisition time to 915 s (2 ms), we fill the center of k-space with 40 (1528) single points, and we average over 7 (7) acquisitions to increase the SNR with a total scan time of 29 (37) min. For DRaNSSE, TE1 is set to 60 and TE2 to 10 ms, with a common acquisition time of 1 ms. The overall DRaNSSE scan time was 65 min for 16 averages. Images are reconstructed with ART using 7 iterations and in Eq. (6). Reconstructions takes in a low-end GPU (see “Methods” ).

It is apparent from a qualitative comparison between the raw (unfiltered) image sets IV and VI in Fig. 4b, that DRaNSSE reconstructions feature a higher SNR than with PETRA. Figure 4c shows the SNR along the red dotted lines in Fig. 4b for four different cases, quantitatively reinforcing this observation (further details on these calculations can be found in “Methods” ). Voxels corresponding to bone and marrow tissue both feature a stronger SNR for DRaNSSE when reconstructed with ART (unlike with Fourier methods, which we discuss below). This is consistent with the fact that signal acquisition in DRaNSSE starts when spins are in full phase coherence after the refocusing pulse, while in PETRA the finite dead time means that spins have already started to dephase before data acquisition. Consequently, for SE1 echo times comparable to the dead time in PETRA, a stronger signal is expected in DRaNSSE, which also makes it possible to use an echo time TE1 shorter than the shortest dead time in PETRA. Furthermore, the simultaneous acquisition of both images in DRaNSSE allows for more averaging in the same total scan time. In the scans in Fig. 4 we acquired 16 averages with DRaNSSE, versus 7 for both short and long dead times with PETRA. In addition to the number of averages, the sampling rate also influences the noise level. We used 24 kHz for images IV, VI and VII and 8 kHz for image V. All in all, the SNR is a factor higher with DRaNSSE than with PETRA for soft tissues, and for hard tissues.

A similar comparison between DRaNSSE and PETRA acquisitions of a rabbit head shows that tissue contrast is also enhanced with the former. Figure 5 shows ART reconstruction slices from DRaNSSE (top) and PETRA (bottom) acquisitions. The right (left) column reconstructions show soft (and hard) tissues. Common to both sequences are: a total scan time of  min, a flip angle of  degrees, a repetition time of 50 ms, a field of view of  mm, an isotropic voxel resolution of 1 mm, a sampling rate of 26 kHz (5.2 kHz for long PETRA), and a total of 1426 radial spokes in k-space, corresponding to an undersampling factor of 5. For the short (long) PETRA acquisition we set the dead time to 90 s (1 ms), the radial acquisition time to 910 s (4 ms), we fill the center of k-space with 64 (496) single points, and we average over 12 (9) acquisitions to increase the SNR with a total scan time of 15 (15) min. For DRaNSSE, TE1 is set to 60 and TE2 to 10 ms, with a common acquisition time of 1 ms. The overall DRaNSSE scan time was 31 min for 26 averages. All images are produced using and 7 ART iterations. Reconstructions take in a low-end GPU (see “Methods” ).

Figure 5

Rabbit image slices obtained with DRaNSSE (top) and PETRA (bottom) for short time parameters (left) and long time parameters (right).

The red arrows in Fig. 5 show the upper (A) and inner (C) teeth (see Fig. 2b). The dead time of 1 ms in PETRA (image XI) is not long enough to fully remove the hard tissue signal, making it difficult to identify hard tissues from a comparison with X. One possible solution is to increase the dead time, but this requires longer overall scan times. For instance, if we double the dead time to  ms (still significantly shorter than the 10 ms echo time of in IX), then we need to increase the number of single points at the center of k-space by a factor . This adds up to 3968 single points keeping the rest of the above settings, prolonging the total scan time to  min. A second possibility is to increase the readout window proportionally to the dead time, i.e. going to an acquisition time of 8 ms in our example. However, such a long acquisition time is affected by decay, which degrades the sharpness of the image contours. This effect is already visible in XI for an acquisition time of 4 ms. Previously explored alternatives to short-T imaging[23] are the acquisition of a second gradient echo by flipping the gradient polarity[17], or the reduction of the number of single points by means of a ramped hybrid encoding[24]. However, they both lead to artifacts in our setup due to eddy currents induced during gradient switching. Another strategy is to acquire different sets of radial spokes with different gradient strengths, a technique known as WASPI (water- and fat-suppressed solid-state proton projection imaging)[25]. Our attempts to image with WASPI also led to artifacts, in this case due to effects (k-space points at the radial boundary between acquisitions are brighter in the strong gradient configuration).

A further advantage of DRaNSSE with respect to PETRA is that the relevant decoherence time constant is rather than , as a result of the induced spin echoes. This means we can wait until hard tissue signals have been strongly suppressed. Consequently, a comparison between images VIII and IX, where SE2 = 10 ms, shows a starker tissue contrast in regions A and C than between images X and XI.

Finally, a concluding remark regarding DRaNSSE: while it clearly performs better than PETRA in our setup, we occasionally observe reconstruction artifacts originated by stimulated echoes. One possible strategy to overcome this is to tailor the timings so as to not capture them, although it might prove inconvenient. As usual in spin echo sequences, it is therefore relevant to carefully calibrate coherent manipulations in the setup, to minimize the detrimental effect of stimulated echoes.

Fourier transform vs algebraic reconstruction techniques

Along with the choice of pulse sequence, the mathematical method employed strongly impacts the reconstruction quality. We used ART in most of the reconstructions shown so far, since we find it performs better than traditional Fourier analysis with our radial k-space sampling. Here we look into this topic in more depth.

The results in Fig. 4b allow for a first comparison. One obvious conclusion is that images IV and VI differ significantly in the encircled regions, where Fourier techniques artificially darken the area, whereas the signal intensity with ART is more homogeneous. We only observe this for low SNR conditions in the system, and we do not have an unambiguous explanation to this effect (this is not due to the excitation profile of the hard rf pulse, we do not use regularization terms in the iterative process, and our Encoding Matrix does not account for gradient non-linearities, magnetic field inhomogeneities or dephasing/relaxation effects). A second visible effect is that the contrast between the small lump of soft tissue and the external surface of the bone where it is attached (white arrows), is higher with ART than with Fourier transformation.

Figure 6 shows noise and SNR maps for ART and Fourier reconstructions of the cow femur (Fig. 6a with DRaNSSE and Fig. 6b with PETRA) and rabbit images (Fig. 6c,d both with PETRA and 1 mm and 0.5 mm, respectively), following procedures described in the “Methods” . We observe significant differences between the noise pattern obtained from Fourier reconstructions respect to ART. The former present a highly non-stationary noise for all images (i.e. the variance of noise depends on position). This is reflected by the coefficient variation (CV) of noise, as well as the average noise level, both shown in Table 2. We also observe for Fourier reconstructions that the maximum noise level happens at the gradient isocenter position, which is coherent with the radial sampling scheme due to the high sampling density close to the center of k-space. For ART, the matrix size strongly affects the level of non-stationarity, where the noise variation coefficient is 8 times smaller when compared with Fourier reconstruction for Fig. 6d, with double matrix size. All in all, we observe that, while PETRA images reconstructed with ART present a more stationary noise than Fourier reconstructions, the noise obtained with DRaNSSE is highly non-stationary even when reconstructed with ART.

Figure 6

Noise and SNR maps for ART and Fourier transform for images VI (a), V (b), X (c) and I (d).

Table 2

Estimated noise and SNR parameters corresponding to images in Fig. 6. SNR in soft and hard tissue on Fig. 6c,d correspond to tongue and inner tooth, respectively. CV stands for coefficient of variation (the standard deviation divided by the mean) and is the average operator.

Data set	Figure	⟨σ(x)⟩ ART / FT	CV(σ(x)) ART / FT	⟨SNR(x)⟩ ART / FFT	SNR soft tissue ART / FT	SNR hard tissue ART / FT
Cow 1 mm DRaNSSE	Fig. 6a	0.81/0.60 (×10^-7)	0.35/0.23	3.2/2.8	5.7/7.2	5.0/3.9
Cow 1 mm PETRA	Fig. 6b	0.92/0.74 (×10^-7)	0.14/0.20	2.4/2.4	4.7/4.0	2.1/1.4
Rabbit 1 mm PETRA	Fig. 6c	1.48/1.43 (×10^-7)	0.16/0.37	2.9/2.8	4.9/3.8	2.3/2.5
Rabbit 0.5 mm PETRA	Fig. 6d	4.10/4.84 (×10^-7)	0.05/0.41	2.4/1.6	2.8/2.2	2.2/2.2

In addition to the noise performance, there are also important differences in the SNR. Figures 6a,b show that ART reconstructions for both DRaNSSE and PETRA produce higher SNR than Fourier reconstructions in the hard tissue. Generally, ART produces SNR values that, in the worst case, are similar to those from Fourier reconstructions (see Fig. 6c,d and Table 2). Moreover, note that ART yields better results than Fourier transforms both when k-space is fully sampled (Fig. 6a,b) and when it is undersampled (Fig. 6c,d). Finally, since the results with ART depend on the reconstruction parameters and , and there is therefore room for potential improvement with respect to the results we present here, we conclude that ART is better suited for image reconstruction than Fourier analysis with our settings.

Discussion

In summary: we have demonstrated the capability of our new low-cost “DentMRI - Gen I” scanner to simultaneously image hard and soft biological tissues; we have devised a new pulse sequence (DRaNSSE) that, compared to standard short sequences such as PETRA, yields higher SNR images and enhanced tissue contrast; and we have shown that iterative techniques (ART) outperform traditional Fourier methods in all quantified metrics, except for the computational time required for the reconstruction.

Dental imaging at low field strengths can be advantageous for reasons beyond cost reduction[14]. For example, magnetic susceptibility effects, which lead to artifacts in the presence of metallic fillings or implants, scale linearly with field strength. Besides, times are shorter for weaker fields (which allows for faster repetition rates), coherence times are longer (somewhat mitigating the SNR loss due to the lower field strength), and specific absorption rate (SAR) limitations scale with the square of the field strength (enabling long echo trains for soft tissue imaging). Also off-resonant effects due to chemical shifts in the Larmor frequencies of different tissues are suppressed (tens of Hertz at 260 mT). Despite these advantages, we are not aware of previous dental MRI systems operated at sub-Tesla fields.

Nevertheless, a prospective “Gen II” scanner must pose solutions to significant remaining challenges if it is to be compatible with clinical conditions. First and foremost, scan times must be reduced (see last column in Table 1). For example, the acquisition times for the human teeth images shown are intolerably long. We therefore plan three major upgrades for “Gen II”: i) we will make use of balanced steady-state free precession protocols, which yield optimal SNR values and can be safely used at low magnetic fields[15,26]; ii) we will use quantum dynamical decoupling techniques such as WAHUHA[27] or CHASE[28] sequences, which can prolong the lifetime of the magnetic resonance signal of hard tissues; and iii) we will perform dual species MRI on protons and P nuclei, since the latter are more abundant in hard biological tissues and they provide complementary information[29]. Additionally, we are starting to investigate the possibility of slice selection with zero-echo time sequences for fast 2D imaging[30]. Aside from shortening scan times, a reduction in the overall footprint and weight would ease scanner siting. The heaviest (and most expensive) component is the permanent magnet (940 kg and 40 k€), which on top imposes the use of a sizable mechanical support structure. For this reason, we will place the new magnet with the yoke at the bottom and the poles pointing upwards (forming a “U”) so patients can lie on a flat bed. This reduces greatly the complexity and footprint of the required mechanical structure. Another important aspect is the radio-frequency system, which will differ significantly from the “Gen I” configuration. In “Gen II”, we will adapt the size and geometry to the human physiognomy, where neck and shoulders place strong boundary conditions, and we will count on dedicated transmit and receive channels, where the latter can be miniaturized and placed in direct contact with the subject’s face or teeth increasing significantly the SNR[21].

Besides dental imaging, a low cost MRI scanner capable of detecting hard tissues and solids may find application in different scenarios, including but not restricted to: head and extremity imaging[14,31]; the food industry (e.g. in inspection and selection tasks, see[32]); or rubber degradation control, which is relevant in, for example, the industries of mining or transport[33].

Methods

“DentMRI - Gen I” Scanner

Main magnet

In MRI systems the quantization axis, spin resonance (Larmor) frequency and, to a large extent, the signal-to-noise ratio of the detected signals are typically determined by a strong homogeneous magnetic field[1]. Here we employ a “C”-shaped permanent NdFeB magnet (Sabr Enterprises LLC, Fig. 1a), providing a main evolution field of 260 mT. The field is shimmed down to an homogeneity of 20 parts-per-million (ppm) over a spherical region of 150 mm in diameter. The magnet is the heaviest component of the scanner with a mass of around 940 kg for a distance between magnet poles of 220 mm. In order to distribute its weight over a 5 m floor surface (compliant with typical architectural regulations), we designed and constructed a support structure consisting mostly of profiles manufactured out of Aluminum 6063. Although this was specified to be non-magnetic, it induces field-strength inhomogeneities at the 4000 ppm level. These are mostly linear and we regularly shim down to 4 ppm during operation with the magnetic gradient system.

Magnetic gradient system

Linear magnetic field inhomogeneities (magnetic gradients) spatially encode information in MRI setups[1]. Intense gradient fields lead to efficient encoding, i.e. strong resolving power. “Gen I” is equipped with a gradient system capable of reaching strengths  mT/m along any spatial direction, sufficient for sub-mm resolution. The gradient system consists of three pairs of planar coils. Each pair produces a magnetic field pointing in the same direction as the main field, and a gradient in magnitude along a Cartesian axis in the laboratory frame of reference. Each of the z (x, y) coils is formed by one (two) lobes (see Fig. 1a). We fabricated the coils with hollow-tube Oxygen-Free-High-Conductivity (OFHC) copper, with low ohmic losses ( per pair) and the possibility of heat removal by running cooling water through the inner conduct. We can stably operate at 100% duty cycle with a total power consumption  kW. In order to pulse the magnetic gradient fields in short times of 50-100 s, we wound each coil into only a few loops for a total pair inductance in all three cases. The size and disposition of the coils ensures deviations % from perfect linearity over a spherical region of 10 cm in diameter. The space between loops is filled with translux D150 epoxy to avoid short circuits and improve mechanical stability. We drive the coils with bipolar amplifiers from International Electric Co. (GPA-400-750), which can ramp from 0 to ±400 A in 100 s with our loads.

Radio-frequency system

In MRI, rf electronics are required to excite the sample spins and detect their response for amplification, digitization and, ultimately, data processing as required for image reconstruction[1]. We employ a single coil for both sample excitation and signal detection. The solenoid, of length 100 mm, 52 mm in diameter and with 20 windings with copper wire of 0.4 mm diameter, is shown in Fig. 1a. The solenoid includes a gap capacitor to homogenize the current along the coil wire, resonant at 11 MHz with a quality factor . We keep the Q intentionally low to achieve excitation and detection bandwidths  MHz, compatible with high spatial resolution images. The rf electronics for fine tuning and impedance matching of the resonant circuit are placed next to the solenoid (Fig. 1b). Both the coil and matching electronics are enclosed in conducting Faraday boxes to avoid the deleterious effect of electromagnetic noise in the laboratory and coupling to the gradient coils. We 3D-printed the housing structure out of polylactic acid for mechanical stability (Fig. 1a).

In transmission (Tx) mode, we drive the solenoid from an rf power amplifier (RFPA-4/11-2000 from Barthel HF-Technik GmbH) fed by a direct digital synthesizer on our field-programmable-gate-array (FPGA) board (RadioProcessor-G from SpinCore Technologies Inc.). The RFPA output is low-pass filtered and sent to the solenoid coil (probe) for sample excitation. This scheme is shown in Fig. 1b.

A TxRx switch (Barthel HF-Technik GmbH, dead time at 11 MHz) commutes the system operation from transmission to reception (Rx). A low noise amplifier (Barthel HF-Technik GmbH, gain 39 dB and noise factor 1.0 dB) amplifies the analog signal induced by the precessing protons on the coil. After low-pass filtering, we digitize and digitally down-convert and filter the signal directly on the FPGA board.

Control system and graphical user interface

The FPGA board is the main component in the experimental control system. This is plugged into a Peripheral Component Interconnect (PCI) slot on the motherboard of a control computer, and allows us to: i) generate low power rf pulses which are amplified for coherent spin excitation; ii) generate three independent low-frequency and low-power outputs which are amplified and fed to the gradient coils for spatial information encoding; iii) read in, digitize, down-convert and filter the MR signal emitted by the sample; iv) execute all of the above operations synchronously; and v) communicate bidirectionally with the control computer.

We have programmed our own Graphical User Interface (GUI) in Matlab, where we design pulse sequences, set the individual pulse parameters independently, configure data acquisition and filtering settings, and visualize the received data and reconstructed objects. Additionally, we have written an intermediate layer in C/C++ to interact from the GUI with the provided drivers for the RadioProcessor-G.

To perform image reconstruction we use MATLAB in a PC with an Intel Core i7-7700 CPU (6.65 GHz, 4 main processors, 8 logic processors), an NVIDIA GeForce GTX 750 Ti GPU (640 CUDA cores, 1.4 TFLOPS, 2GB memory) and 32 GB RAM.

Pulse sequences

MRI pulse sequences are designed to manipulate the sample magnetization: resonant rf pulses lead to coherent rotations, and magnetic-field-gradient (or simply gradient) pulses encode spatial information. The short lifetime of hard tissue signals imposes the use of special-purpose pulse sequences. These typically force a high-intensity constructive interference signal (echo) over very short timescales, or even simultaneously excite the sample and detect its response[3]. Ultra-short echo time (UTE[34]) and zero-echo time (ZTE[16]) pulse sequences are prominent examples of the former approach, and SWeep Imaging with Fourier Transformation (SWIFT[13]) of the latter.

ZTE and SWIFT are better suited than UTE for ultra-short tissues such as those present in teeth. If the scanner hardware allows for rf excitations with sufficient bandwidth, ZTE is the most convenient choice[35]. Since the “DentMRI - Gen I” scanner has a rather small field of view and we have a fast 2 kW rf power amplifier, we acquired all images above with two variations of standard ZTE pulse sequences (Fig. 7): Pointwise Encoding Time Reduction with Radial Acquisition (PETRA[17]); and a sequence we have devised to overcome specific k-space (spatial frequency space) coverage and contrast limitations, Double Radial Non-Stop Spin Echo (DRaNSSE).

PETRA

In PETRA (Fig. 7a), a hard rf pulse homogeneously excites the sample after the onset of the encoding gradient fields. Data acquisition starts next, after a dead time usually set by the TxRx switch and the rf coil ring-down. Every acquisition follows a radial direction (spoke) in a 3-dimensional k-space. The center of k-space is not sampled in this way due to the finite dead time, so it is filled in a pointwise manner on a Cartesian grid. Once all radial spokes have been sampled, object reconstruction can be carried out either by inverse discrete Fourier transformation (which requires regridding and interpolation operations), or with Algebraic Reconstruction Techniques (ART[18-20]).

PETRA suffers from two main constraints. On the one hand, scan times can become exceedingly long when the pointwise-filled gap is large. This is accentuated at low fields, due to the slow ring-down of the Tx coil, which happens with a time constant . On the other hand, it is not possible to have contrast with PETRA, since the echo time (TE) is defined to be zero. Our solution to both problems is our new sequence DRaNSSE (Fig. 7b). When we do acquire with PETRA, we can subtract two independent scans with different parameters for contrast, following a previously existing scheme with UTE pulse sequences[36]. In one acquisition dead times are short, so hard tissues still appear bright, while in the other we use intentionally long dead times, to ensure hard-tissue contributions fade away before readout.

DRaNSSE

In DRaNSSE, after the initial rf pulse, two refocusing -pulses produce two different spin echoes. The first one (SE1) is induced at an echo time (TE1) as short as possible to include contributions from both hard and soft tissues. The second echo (SE2) occurs at a later time (TE2), when the short-lived signal from hard tissues has already faded away. A first acquisition takes place between SE1 and the second refocusing pulse, corresponding to radial spokes in k-space from to the maximum sampled value (); the second acquisition starts after the second -pulse and ends at TE2, sampling from to . The second acquisition could go beyond SE2, allowing for higher quality imaging of soft tissues. However, the occurrence of multiple unwanted stimulated echoes due to small experimental imperfections (STE1-STE3 in Fig. 7b) can degrade image quality. STE1 appears a time after the second -pulse, while STE2 and STE3 occur at times and after , respectively. Since in our sequences is much shorter than (the time between the start of the first acquisition and the center of the second refocusing pulse), we find it convenient to acquire the signal before SE2 to avoid reconstruction artifacts. Note that DRaNSSE is similar to radial multi-echo spin-echo with an echo train length of two[37], but with half-echo acquisition to minimize TE1 and avoid contamination from stimulated echoes.

Fourier and algebraic reconstruction techniques

Spatial encoding in MRI relies on inhomogeneous magnetic fields, which provide a Larmor or spin-precession frequency () dependent on the position of the nuclei in the Region of Interest (RoI). Mathematically, this can be expressed aswhere is the gyromagnetic ratio ( MHz/T for protons) and is the magnetic field at position at time t. As the pulse sequence evolves, the phase acquired by the spins depends on their position:During their precession, spins induce a time-varying signal with the interference of all spins on a nearby detector:where is the spin density distribution in the RoI. This signal is then digitized during a readout or acquisition window, and we apply one of the following mathematical tools to reconstruct an image.

Discrete Fourier transform reconstruction

Typically, scanners make use of linear gradient fields, in the presence of which the Larmor frequency varies also linearly with position. In this scenario, an inverse Fourier Transformation (FT) of s(t) suffices to reconstruct , since the integral in Eq. (2) after down-mixing is trivial and Eq. (3) becomesHere we have assumed, without loss of generality, that a gradient of strength points along the z-axis. Since the detected signal is discretized in time by an analog-to-digital converter, fast Fourier transform protocols can be used to map a discrete k-space onto a discrete reconstruction in real space[1].

Discrete Fourier transforms are computationally efficient and applicable in many useful scenarios, but this simple approach may yield suboptimal results when the sampling data does not adjust to a Cartesian grid in k-space. This occurs, for instance, if gradients are time dependent, discretization times are not homogeneously distributed or, as in the present work, k-space is sampled radially or following curved trajectories. In such cases, the acquired data must be pre-processed (e.g. with regridding, density compensation or interpolation operations) prior to Fourier transformation[38]. To perform Fourier reconstruction we first interpolate the k-space data to a Cartesian grid and then apply a Fast Fourier Transform protocol to the interpolated data.

Reconstruction based on an encoding matrix

An alternative to pre-processing and Fourier operations in non-Cartesian k-space sampling is to build a linear forward model for the system response to the applied pulse sequence, define a cost function for the reconstruction, optionally add regularization terms to penalize unrealistic results, and solve a linear inversion problem[39].

Once the time-dependent signal resulting from the interference of the precessing nuclei has been recorded and discretized, s(t) becomes a vector of length equal to the number of time steps , becomes a vector of length equal to the number of voxels , and becomes the Encoding Matrix with rows and columns. Equation (3) thus changes toand can be obtained by direct inversion of the Encoding Matrix as , or by any other means of solving the system of linear equations, e.g. by iterative algorithms such as Algebraic Reconstruction Techniques (ART)[18-20]. ART estimates based on the recursive equationwhere is a control parameter, is the component in vector , is the row in the encoding matrix , and can be set to zero. The estimated solution is updated times through Eq. (6), where stands for the overall number of ART iterations.

Although iterative methods such as ART are computationally slow compared to discrete Fourier transforms, we find in this work that they can vastly outperform Fourier reconstruction in other relevant metrics.

Estimation of the SNR

We assume that the final reconstruction signals in image domain follow a non-stationary Rician distribution[40], i.e., a Rician distribution where the noise parameter becomes position dependent: . This assumption is based on two facts:

Gaussianity and Rician: ART and FT reconstructions both produce signals with additive Gaussian Noise in each of the points. The magnitude of that Gaussian data produces Rician data if the real and imaginary parts of the former are independent and have the same parameter. Although this need not be strictly the case with our reconstruction processes, we show below that the Rician approach yields valid results.

Non-stationarity: ART and FT both perform a linear reconstruction of complex Gaussian data where different weights are applied to the original values. As an effect of the weights and the correlations introduced by the reconstruction, the final reconstructed data can show different noise properties in each voxel, i.e., every voxel can show a different value of . In our images, as is usually the case, varies slowly with and can be considered a low frequency signal[40,41].

In order to test the first assumption (the Rician nature of the data), we consider Fig. 6d, where the background noise is to a good approximation stationary and a single value of characterizes the distribution. The expected signal strength in the background is zero, and the Rician Probability Distribution Function (PDF) simplifies to a Rayleigh distribution. From the fits in Fig. 8 we conclude that the signal in the background is Rayleigh-distributed, which is compatible with our Rician assumption.

Figure 8

Fit of the background data of image I in Fig. 2 to a Rayleigh distribution. A fit to a Normal distribution is shown for comparison.

Thus, in order to estimate the SNR of the final volume we assume a non-stationary Rician model. We define the estimated SNR aswhere and are, respectively, an estimation of the original signal (in absence of noise) and and estimation of the variance of the noise in each point of the image.

In order to estimate the noise, we use a 3D extension of the homomorphic scheme proposed in Ref.[41]. The rough estimation of the original signal is provided by the second order moment of a Rician distribution. If is a Rician signal, then[40]We estimate the expectation by a local mean, calculated as the convolution of the original signal with a average kernel, . Hence:and the estimated SNR becomes:Note that, due to the use of the convolution kernel, the SNR will show some blurring in the edges of the image. However, for the comparison purposes of this work, that will not affect results.

Ethical statement

All animal parts were obtained from a local butcher and research was conducted following the 3R principles. Experiments using human teeth were carried out following Spanish regulations and under the ethical consent from the Spanish National Research Council (CSIC), agreement number 058/2020. The Ethical committee from Spanish National Research Council (CSIC) has approved the protocol followed for the study. Tooth samples were provided by Clínica Cortell Llobell S.L. and informed consent was obtained from the participants prior to study commencement.