Silent 3D MR sequence for quantitative and multicontrast T1 and proton density imaging

This study aims to develop a silent, fast and 3D method for T1 and proton density (PD) mapping, while generating time series of T1-weighted (T1w) images with bias-field correction. Undersampled T1w images at different effective inversion times (TIs) were acquired using the inversion recovery prepared RUFIS sequence with an interleaved k-space trajectory. Unaliased images were reconstructed by constraining the signal evolution to a temporal subspace which was learned from the signal model. Parameter maps were obtained by fitting the data to the signal model, and bias-field correction was conducted on T1w images. Accuracy and repeatability of the method was accessed in repeated experiments with phantom and volunteers. For the phantom study, T1 values obtained by the proposed method were highly consistent with values from the gold standard method, R2 = 0.9976. Coefficients of variation (CVs) ranged from 0.09% to 0.83%. For the volunteer study, T1 values from gray and white matter regions were consistent with literature values, and peaks of gray and white matter can be clearly delineated on whole-brain T1 histograms. CVs ranged from 0.01% to 2.30%. The acoustic noise measured at the scanner isocenter was 2.6 dBA higher compared to the in-bore background. Rapid and with low acoustic noise, the proposed method is shown to produce accurate T1 and PD maps with high repeatability by reconstructing sparsely sampled T1w images at different TIs using temporal subspace. Our approach can greatly enhance patient comfort during examination and therefore increase the acceptance of the procedure.


Introduction
High resolution, T1-weighted (T1w) imaging is widely used in neuroscientific and clinical applications, for example, with demyelinating diseases (e.g. multiple sclerosis) (Spies et al 2013), brain development and pediatric imaging (Roque et al 2014), brain aging (Maniega et al 2015). The most common way of generating T1w images is by magnetization-prepared rapid gradient-echo or fast spin echo (FSE) sequences. Quantitative T1 mapping provides objective tissue characterization for evaluating pathology, and facilitates statistical modeling, multi-center and longitudinal studies (Tsialios et al 2017). The gold standard of quantitative T1 mapping is performed by inverting the longitudinal magnetization and then sampling the MR signal at multiple inversion times (TIs) along the exponential recovery curve, known as inversion recovery (IR) T1 mapping (Stikov et al 2015, Tsialios et al 2017. The method provides robust T1 quantification, but lengthy scanning time limits its clinical application. Widely used alternative methods to speed up the scanning include look-locker and the variable flip angle (VFA) method (Stikov et al 2015). However, the VFA method is highly sensitive to B1 errors, and requires careful consideration of the spoiling regimes, pulse sequence parameters and magnetization transfer effect (Schabel andMorrell 2008, Stikov et al 2015), whereas the look-locker method presents a tradeoff between spatial resolution, number of images reconstructed to sample the T1 relaxation curve and scanning time especially for 3D imaging (Henderson et al 1999).
T1 imaging can be accomplished by various k-space acquisition strategies, i.e. with the most commonly used Cartesian trajectory, or radial, spiral and other non-Cartesian trajectories. Radial imaging has wide clinical applications because it is less sensitive to motion artifacts and allows for reduced field of view (FOV) imaging (Chandarana et al 2011). A study on pediatric abdominal imaging suggested that 3D radial outperforms its Cartesian counterpart due to motion insensitivity (Roque et al 2014). The disadvantages of radial imaging are complications in reconstruction and k-space trajectory design, sensitivity to gradient inaccuracy and off-resonance effect (Wright et al 2014). Various radial trajectories have been proposed to improve the uniformity of readout distribution and sampling efficiency (Wong andRoos 1994, Lingala et al 2013). With each radial spoke, an equal amount of low and high spatial frequency data is collected to provide homogeneous image updates. This feature makes radial imaging well-suited for time-resolved, dynamic techniques which need to capture motion kinetics or contrast variations and have wide applications in angiography, functional MRI, dynamic contrast-enhanced imaging, cardiac imaging etc (Winkelmann et al 2006, Tsao andKozerke 2012). These techniques require radial profiles designed to allow uniform k-space sampling in each time frame. Winkelmann et al (2006) proposed an optimal radial profile based on the Golden Ratio, which yields a high flexibility in choosing an appropriate temporal resolution and is later adapted by Ehses et al (2013) for quantitative T1, T2 and proton density (PD) mapping using TrueFISP.
Several methods have been developed to speed up dynamic MRI using reduced data acquisition and without significantly compromising image quality. Based on the assumption that the contrast information is mostly contained in the low frequencies, several studies have employed view sharing (Ehses et al 2013, Kecskemeti et al 2016 for high resolution dynamic imaging with radial acquisition, mostly by using a k-space filter to overcome the challenge in obtaining high spatial and temporal resolution simultaneously. Other methods exploit redundancy in dynamic MRI that the image series typically exhibit spatiotemporal correlations, and can be reconstructed from undersampled k-t space data (Jung et al 2007, Lingala et al 2013. Model-based approaches called k-t BLAST/SENSE (Jung et al 2007) resolve aliasing artifacts in Fourier reciprocal x-f space using signal covariance learned from low-resolution training data. The method is widely applicable especially for the imaging of objects exhibiting quasiperiodic motion, such as heart, lung and abdomen. Tamir et al (2017b) applied low-dimensional temporal subspace constraints learned from signal model on dynamic FSE data, and reduced T2 blurring while generating multicontrast T2 weighted images at virtual echo times (TEs). These acceleration techniques are especially important for 3D imaging to achieve decent spatiotemporal resolution. Compared to 2D imaging, 3D imaging is free from non-perfect slice profiles and provides isotropic resolution which has the advantage of being reformatted and viewed in different planes.
Previous studies explored the use of silent T1w sequences for intracranial tumor patients and other clinical population (Alibek et al 2014, Ida et al 2015, Holdsworth et al 2018. Those sequences were achieved by using the magnetization prepared rotating ultrafast imaging sequence (RUFIS) (Madio and Lowe 1995), a 3D radial technique with minimal switching of frequency encoding gradient, leading to nearly inaudible scanning. These studies have reported comparable efficiency of silent T1w scanning in terms of signal/contrast to noise ratio and lesion conspicuity rated by radiologists, suggesting they can be used as viable quiet alternatives for conventional T1w imaging (Alibek et al 2014, Ida et al 2015, Holdsworth et al 2018. The acoustic noise originating from the rapid gradient switching can be up to 130 dBA at 3 T, causing patients discomfort, anxiety and even hearing loss, limiting the use of MRI especially in pediatric imaging, patients with hyperacusis etc (Schmitter et al 2008). Studies have also revealed vulnerability of the human fetus to excessive acoustic noise while external noise protections (ear plugs) are not applicable in fetal examinations (Committee on Environmental Health 1997). In addition, even with external noise protections, exposure to acoustic noise during 3 T MRI scan can cause altered cochlear function and temporary increase of hearing threshold (Radomskij et al 2002, Jin et al 2018. It has been shown that improved patients comfort during MRI scanning, mainly by reducing noise and confined sensation, is essential to reduce mental stress and claustrophobia occurrence during examination (Flaherty andHoskinson 1989, Dewey et al 2007). Therefore, silence is a desired feature and healthcare providers endeavor to develop noise reduction techniques. Hardware-based solutions, including Lorentz force balancing, rotating DC gradient and active vibration control, may considerably increase the complexity and costs of the MRI system. Whereas sequence-based methods can be incorporated into any existing MRI system and achieve noise reduction of 10-20 dBA (Schmitter et al 2008), thus they are more commonly used by various manufacturers.
The aim of this study was to develop a volumetric, fast and silent approach for quantitative T1 and PD mapping, while producing multiple bias-field corrected (Van de Moortele et al 2009) T1 images at different TIs using IR prepared RUFIS. It was implemented by combing an interleaved radial trajectory design and reconstruction of unaliased images using subspace constraints. The approach was validated in phantom and healthy volunteers, the accuracy and repeatability of the method was demonstrated. After the adiabatic inversion preparation, a train of low flip angle α excitation was applied with short TR. The gradient was active during RF excitation and was gradually reoriented with each TR. A waiting time (TD) was placed between consecutive interleaves to allow for free magnetization recovery. (b) Demonstration of acquisition strategy. Data was segmented along the readout (as illustrated in the dashed boxes on the left), and segments acquired at the same TI were grouped together to generate undersampled images. On the right is a representation of the k-space trajectory at one effective TI, which was first generated according to equation (1), and then realigned into segments with the same effective TI on different interleaves (segments illustrated in dashed boxes, and colors correspond to radial spokes allocated to different interleaves).

IR-RUFIS sequence implementation and undersampling strategy
The IR prepared RUFIS sequence is depicted in figure 1(a). The magnetization is prepared by an adiabatic inversion pulse followed by RUFIS readout, which consists of repetitive small flip angle α excitations and short repetition time (TR). A waiting time (TD) is placed between consecutive interleaves to allow for magnetization recovery. The signal model for original IR-snapshot-FLASH imaging was designed for Cartesian trajectory (Deichmann and Haase 1992), which assumes that the contrast is defined at the time point where the k-space center is acquired. When using radial trajectories, the center of k-space is constantly updated by each projection, making it possible to sample signal evolution along the temporal dimension.
Different from standard sequences, RUFIS utilizes nonselective hard pulse excitations in presence of the encoding gradients, immediately followed by 3D radial center-out data acquisition, resulting in nominal zero TE. After data acquisition in each TR, signal spoiling is achieved by the combination of gradient ramping before next excitation and RF phase cycling. To achieve uniform excitation in the FOV, short and intense RF excitation pulses are used to cover the imaging bandwidth spanned by the gradient. The achievable flip angle is limited by the B1 amplitudes (typically~15 µT on clinical MR scanners). Small and smooth gradient update between successive repetitions leads to robustness against eddy currents and silent scanning. Because of the small flip angle and short TR, the original RUFIS is PD weighted.
To achieve high spatial and temporal resolution, the data acquisition is conducted in an interleaved manner. Data segments with similar TI (adjacent position in the readout) from different interleaves are grouped together to generate undersampled images. The effective TI for each time frame is defined at the center of the acquisition window. In the work of Wong and Roos (1994), an interleaved radial trajectory was designed with the end point of each radial spoke uniformly distributed on a surface of sphere, and each interleave was rotated from one another. The expression for x, y, z components of the interleaved trajectory are as follows: where M is the number of interleaves, N is the number of spokes per interleave, n = 1, 2, ….,N and m = 1, 2,…,M. The current work needs a trajectory where: (i) each time frame contains data sample roughly uniform across 3D k-space and (ii) the stepwise gradient update should be small enough to maintain the silent feature of the sequence. We adapted the trajectory from equation (1) to meet the above requirements by realigning the trajectory interleaves into each time frame. The acquisition strategy and trajectory design are illustrated in figure 1(b).

MRI data acquisition
All experiments were conducted on a GE 3T MR750w scanner with a 12-channel GEM head array coil (GE Healthcare, Waukesha, WI). The general acquisition parameters for phantom and in vivo experiments were as follows: flip angle = 2 • , readout bandwidth = ± 15.6 kHz; The data acquisition began 40 ms after the IR pulse and lasted for approximately 5000 ms to acquire 2048 radial spokes for each interleave. Four dummy cycles were performed to reach steady state, followed by 32 interleaves of actual data acquisition. A waiting time of 1000 ms was applied between consecutive interleaves to allow signal recovery.
A phantom consisting of samples with different T1 values (Diagnostic Sonar, Livingston, UK) was used in this study. The T1 values of the phantom ranged from 300 to 1565 ms, covering the range of physiological T1 values expected in brain white matter (WM) and gray matter (GM). FOV = 20.1 cm, isotropic resolution = 1.5 mm, and TR = 2.46 ms. For comparison, the center slice of the phantom was also measured with a Cartesian inversion-recovery spin-echo (IR-SE) sequence at 25 different TIs ranging from 50 ms to 4000 ms, and TR = 13 s which is sufficient long to allow for a near-complete recovery. The reference T1 value was assessed by a three-parameter fit of the acquired images.
Six healthy volunteers, four males and two females, ages ranged from 28 to -33 years, participated in the data acquisition. The current study was approved by the local ethnic board and prior informed consents were obtained. FOV = 18.5-21 cm, isotropic resolution = 1.5 mm, TR = 2.30-2.56 ms, and scanning time 3:30 min. Two scans of each subject were acquired to test the repeatability of the proposed method.

Temporal subspace and reconstruction
The signal curve observed in the IR-RUFIS experiment is a function of the tissue parameters (T1, PD) and acquisition parameters (IR pulse and readout flip angle). Figure 2(c) shows the example of the signal evolution for different T1 values. The signal curves of different T1 values follow similar trends, which implies that they can be represented by low-rank subspace. Bloch simulation can be used to generate the dictionary M for T1 values of interest, where each column represents the signal evolution for a particular T1 value. By conducting singular value decomposition (SVD) of the dictionary, the subspace Φ K can be defined as the first K right singular vectors, and used to approximate the real signal: Where the modeled signal M can be described by the temporal basis coefficients α. Figure 2(c) illustrates the use of subspace to approximate signal evolution. The size of the subspace is chosen as a tradeoff between noise amplification and modeling error, as the noise variance in a subspace-constrained reconstruction increases linearly with the subspace size K (Tamir et al 2017b). The relative modeling errors for different subspace sizes are defined using the Frobenius norm and shown in figure 2(a). The radial k-space data was first segmented into time frames and then interpolated into a Cartesian grid by convolving the spokes by the commonly used Kaiser Bessel kernel. A dictionary was generated by Bloch simulation using T1 values from 100 ms to 5000 ms with 1 ms spacing, from which the temporal basis was calculated. The subspace size was chosen based on the largest allowed modeling error described in equation (3). In this study, we chose subspace size K = 4, which makes the relative modeling error less than 1% for all T1 values used in the simulation. The k-space data was projected onto the subspace before inverse Fourier transform to obtain K temporal coefficient images, which were used to estimate the coil sensitivity map. After Figure 3. Overview of the reconstruction process. 1. Undersampled radial k-space data in each time frame was interpolated into Cartesian grid and non-uniform sampling density was compensated. 2. K-space data time series were projected onto the temporal subspace, which was calculated from Bloch simulation. In the current study, number of time frames is 32 and subspace size K is 4. 3. Inverse fast Fourier transform, and coil combination were conducted in the temporal subspace to generate K coefficient images. 4. Temporal coefficient images were projected back to the original time series. 5. Quantitative T1 and PD maps were obtained by fitting the reconstructed T1w images. coil combination, the temporal coefficient images were projected back to the original time series. The overview of reconstruction process is illustrated in figure 3.

Parameter fitting and bias-field correction
The quantitative T1 map was computed by a pattern recognition algorithm (McGivney et al 2014) previously used in magnetic resonance fingerprinting, which compares the observed signal to each of the entries in the dictionary and finds the one with the maximum inner product as the best match. Then the T1 value of that entry is retrieved using the lookup table. The PD of each voxel was computed as the scaling factor between the observed signal and the dictionary entry. We employed a bias-field correction method based on image ratio (Van de Moortele et al 2009) to improve T1 contrast. The T1w images were divided by the PD map to eliminate receive RF coil sensitivity B − 1 and PD contrast. The effectiveness of the method was evaluated by comparing whole-brain intensity histograms of T1w images before and after the correction.

Acoustic noise measurement
Acoustic noise was measured using a Bruel&Kjaer sound level meter (Type 2250) equipped with MR compatible microphone (type 4189). Calibration of sound level meter was conducted before the measurement using a 94 dB and 114 dB sound source (Calibrator Type 4231). The acoustic noise was measured in 2 min duration for the in-bore ambient background, the proposed silent sequence (acquisition parameter: TR = 2.46 ms, TI = 40 ms, 1.5 mm isotropic resolution) and a IR prepared fast spoiled gradient echo sequence (acquisition parameter: TR/TE = 7.6 ms/2.8 ms, TI = 900 ms, in-plane resolution = 1.5 × 1.5 mm, slice thickness = 2 mm) as comparison. The A-weighted average (LAeq, dBA) and the C-weighted peak (LCpeak, dBC) sound pressure levels (SPL) were recorded. The microphone was placed in-bore at scanner isocenter inside the head coil. RF transmit was disabled during the measurement to avoid damage to the microphone.

Data analysis
The accuracy of the proposed method was examined in the phantom study by comparing to the gold-standard IR-SE sequence. The mean and standard deviations of T1 values were extracted from manually defined regions of interest (ROIs) in each sample. The Pearson correlation coefficient was calculated against the values obtained from IR-SE. The coefficients of variation (CVs) of the two scans were calculated to test the repeatability of the proposed method.
Subject-specific templates were created by registering the second scan to the first scan of each subject by rigid transformation using the MATLAB image processing toolbox. Brain extraction was conducted on composite T1w image using BET (FMRIB Software Library: www.fmrib.ox.ac.uk/fsl) (Smith 2002) and were visually inspected to ensure the quality. The mask was applied to T1 and PD maps. To visualize the distribution of T1 values, T1 histograms of a single subject and pooled histogram from all the subjects were

Phantom study
The results from the phantom experiment and repeatability test are shown in table 1 and figure 4. The effect of enforcing subspace constraint on signal evolution is shown in figure 4(b). T1 values measured with the proposed method are plotted against reference T1 values derived from the IR-SE experiment in figure 4(c). The proposed method showed high accuracy and only small deviation from the results of IR-SE, with correlation coefficient R 2 = 0.9976. The sample with the smallest T1 value showed the largest percentage difference compared to IR-SE. CVs ranged from 0.09% to 0.83%. T2 values of sample no. 2 and no. 3 were 130.45 ms and 41.90 ms, respectively, calibrated using spin echo. Note that for these two samples, the accuracy of the T1 measurement was not influenced by vastly different T2 values.

In vivo parameter mapping
The axial, sagittal and coronal T1 and PD maps from in vivo imaging are shown in figures 5(a) and (b). The volumetric T1 and PD maps with isotropic resolution can be reformatted into different orientations. The placements of GM and WM ROIs are shown in figure 5(c). Table 2 summarizes ROI analysis results for the in vivo experiment. Mean T1 values from WM and GM ROIs of all subjects ranged from 702 ms to 767 ms and 1033 ms to 1276 ms, respectively, in accordance with literature values (Zhu andPenn 2005, Wang et al 2018).
The PD ratios between WM and GM were calculated by comparing the combined WM and GM ROIs and ranged from 0.81 to 0.86 for all subjects, which were also consistent with literature values (Farace et al 1997).
The quantitative in vivo results exhibit high repeatability. For ROI analysis, the mean and standard deviations of CVs calculated across all subjects and ROIs were 0.40% ± 0.37% (max/min = 1.55%/0.01%). Figure 6(a) shows the whole-brain T1 histogram from a single subject, and figure 6(b) shows the pooled histogram from all subjects. The T1 value peaks of GM and WM can be clearly delineated for all subjects, and the peak values are reported in table 3. For histogram analysis of WM and GM peak values, the mean and standard deviations of CVs across all subjects were 0.58% ± 0.68% (max/min = 2.30%/0.08%). Figure 7(a) shows the representative images with different T1 contrast after bias-field correction, reconstructed at effective TIs ranging from 118 ms to 5000 ms. With the increasing effective TI from left to right, sequential nulling of WM, GM and CSF can be observed. Whole-brain intensity histograms of T1w images before and after bias-field correction are shown in figure 7(b). The correction method significantly improved the separation between GM and WM peaks, indicating a better tissue contrast.

Acoustic noise measurement
The LAeq and LCpeak values of in-bore ambient background were 72.2 dBA and 95.3 dBC; for the IR-RUFIS sequence, they were 74.8 dBA, 95.7 dBC; for the IR prepared fast spoiled gradient echo sequence which is widely used for T1 imaging, the LAeq was 104.9 dBA. The A-weighting is the most common frequency weighting which best resembles the loudness perceived by the human ear. The C-weighting is used in peak SPL measurements for evaluating very high-level or low-frequency sounds. Compared to in-bore ambient noise, the IR-RUFIS sequence caused LAeq and LCpeak acoustic noise levels to increase only 2.6 dBA and 0.4 dBC, respectively. Compared to IR gradient echo sequence, the proposed sequence achieved 30.1 dBA acoustic noise reduction, which corresponds to 97% reduction in sound pressure, since the acoustic measurement is in logarithmic scale.

Discussion
In the current study, we proposed a quantitative T1 and PD mapping method based on interleaved IR-RUFIS and temporal subspace constraint. The method is efficient in generating volumetric T1 and PD maps in clinically relevant resolution and time (whole-brain 1.5 mm isotropic resolution in 3.5 min) while producing multiple time-resolved multicontrast T1w images. T1 values obtained by the proposed method were compared with IR-SE sequence in phantom studies where the correlation coefficient R 2 = 0.9976. CVs ranged from 0.09% to 0.83% in repeatability test. In addition, the sequence is literally silent with the averaged SPL only increasing by 2.6 dBA compared to in-bore background. It is also worth noting that the accuracy of the current T1 mapping method was not influenced by T2 values. As demonstrated in the phantom results of    For in vivo validation, T1 values extracted from WM and GM ROIs are consistent with literature values. The T1 values in the WM range from 702 ms to 767 ms, and in GM 1033 ms to 1276 ms. The whole-brain T1 histogram is characterized by a distinct and sharp peak representing the T1 values found in WM, and a broader peak on the right representing T1 values from GM, extending into larger T1 values from CSF, which is consistent with previous studies (Van Walderveen et al 2003, Stikov et al 2015. The scan-rescan tests of the six subjects showed repeatability of CVs ranged from 0.01% to 2.30%. . The WM histogram peaks facilitate comparison by providing a distinguishable reference, as previous studies confirmed it to be a biomarker for multiple sclerosis (MS). Compared to normal subjects, the WM T1 histogram peak in MS patients is broader, with lower peak amplitude and shifted towards larger T1 values (Van Walderveen et al 2003, Vrenken et al 2006) because of T1 prolongation in normal appearing WM. It can be noted that the CSF peak is not Comparison of whole-brain intensity histogram of T1w images before (blue) and after (red) bias-field correction. The separation of WM and GM peaks is significantly improved after correction. distinctly visible on the T1 histogram in figure 5. This is possibly caused by the partial volume averaging effect (Van Walderveen et al 2003) or underestimation of long T1 values due to reduced dynamic range. We expect the measurement of long T1 values samples, such as CSF, to be more accurate with a longer acquisition window or waiting time.
Quantitative MRI aims to provide absolute measures to overcome bias between scans and sites. However, most T1 mapping literatures are very consistent within studies but not across studies (Stikov et al 2015). Inter-study discrepancies could be caused by a variety of reasons: the fitting routines, different study populations (both gender and age can affect T1) as well as the choice of sequences and parameters. Stikov et al (2015) examined the three most common T1 mapping methods: IR, Look-Locker and VFA. They found excellent agreement between the three methods in phantom experiments. But for in vivo imaging, the Look-Locker underestimates and VFA overestimates T1 values, with the peak of WM T1 histogram varies by more than 30%. These variations could be results of incomplete spoiling and B1 inaccuracy. Even when using the gold standard IR protocol, there are large variations of reported WM T1 values due to different experiment settings. Despite these inter-study discrepancies which make standardization difficult, methods with high repeatability are still reliable and sensitive for detection of changes between groups (Van Walderveen et al 2003, Stikov et al 2015. Poor scan-rescan repeatability would undermine the comparability of data, as the measurement variations could obscure pathological changes (Leary et al 1999). In the current study, the proposed method exhibits high repeatability in both phantom and in vivo experiments, which is advantageous for clinical application and longitudinal studies.
There are several sources of error to consider. T1 values of CSF were underestimated in the volunteer study, as discussed above. In the phantom study, the sample with the smallest T1 value showed the largest percentage error compared to the gold-standard method. The signal recovery is more rapid with a small T1 value; thus, it is more sensitive to the averaging effect within the time frame. We expect to improve the accuracy of small T1 value measurements by using shorter time frames. In addition, when a perfect inversion and excitation RF pulse are not achieved, the signal deviates from the simulation and introduces bias into the T1 mapping. T1 values can be underestimated with incomplete inversion pulse and overestimated when the effective flip angle of excitation pulse is smaller than the nominal flip angle. To account for RF transmit field inhomogeneity, some studies employ B1 mapping methods (Van de Moortele et al 2009, Stikov et al 2015. In the gold standard IR-SE method, the inversion efficiency is considered in the 3-parameter fitting. In the current study, we did not incorporate B + 1 errors into the model and fitting, because it will greatly increase the dictionary size and result in unstable estimation of T1 values. Although the accuracy is high for T1 values of interest in the brain, B + 1 errors should be considered in future studies to reduce systematic error. Some studies (Kecskemeti et al 2016, Tamir et al 2017b have proposed methods to reconstruct multiple images with different contrast along the relaxation curve in a single exam. In comparison to conventional methods, which only acquire a single contrast at the selected time point, these additional images can provide tissue dynamic information and enable subject-specific retrospective selection of the optimized contrast for diagnosis or segmentation purposes. In the current study, 32 T1w images with TIs ranging from 118 ms to 5000 ms were reconstructed, and nulling points of GM, WM and CSF can be observed in the sequential images. However, the B − 1 inhomogeneity impairs the image quality and intrinsic PD weighting tend to reduce the contrast in T1w images. These factors were eliminated using the technique described in Van de Moortele et al (2009). From the whole-brain intensity histogram of the T1w images after correction, the clear separation between WM and GM peaks indicate that the method achieves a good outcome in bias-field cancellation. The resulting purely T1w images have better anatomical quality and brain tissue differentiation, which are more suitable for segmentation, voxel-based morphometry and diagnosis applications (Van de Moortele et al 2009).
As a future extension for the current work, more contrasts can be acquired by replacing the IR with other preparation modules, and the same trajectory can be used to facilitate acceleration. Some studies (Fan et al 2014, Tamir et al 2017a applied IR together with T2/diffusion preparation to reconstruct multi-contrast images. The proposed method has potential application in pediatric imaging, after the motion effects being investigated by experiment or simulation. In addition, the proposed method can be combined with the long-T2 suppression technique to characterize tissue with short-T2 values, such as tendon or myelin.

Conclusion
We proposed a method for volumetric T1 and PD imaging using interleaved IR-RUFIS and temporal subspace constraint. The sequence is literally silent and showed high accuracy and repeatability in phantom and volunteer studies. These features could greatly improve patient comfort during examination and therefore increase the acceptance of the procedure.