Joint multi-field T1 quantification for fast field-cycling MRI.

PURPOSE: Recent developments in hardware design enable the use of fast field-cycling (FFC) techniques in MRI to exploit the different relaxation rates at very low field strength, achieving novel contrast. The method opens new avenues for in vivo characterizations of pathologies but at the expense of longer acquisition times. To mitigate this, we propose a model-based reconstruction method that fully exploits the high information redundancy offered by FFC methods.
METHODS: The proposed model-based approach uses joint spatial information from all fields by means of a Frobenius - total generalized variation regularization. The algorithm was tested on brain stroke images, both simulated and acquired from FFC patients scans using an FFC spin echo sequences. The results are compared to three non-linear least squares fits with progressively increasing complexity.
RESULTS: The proposed method shows excellent abilities to remove noise while maintaining sharp image features with large signal-to-noise ratio gains at low-field images, clearly outperforming the reference approach. Especially patient data show huge improvements in visual appearance over all fields.
CONCLUSION: The proposed reconstruction technique largely improves FFC image quality, further pushing this new technology toward clinical standards.

INTRODUCTION

The magnetic field dependency of the longitudinal and transverse relaxation times, also referred to as nuclear magnetic relaxation dispersion (NMRD), provides insight into the underlying structural order and dynamics of a wide range of molecular systems., In recent years, the dispersion of protons in particular has experienced increased interest for the investigation of biomarkers related to various pathological processes., , , The field‐dependent properties of such biomarkers are invisible to traditional MRI scanners, which operate only at one fixed main magnetic field strength and are restricted to the measurement of relaxation times corresponding to the field used. However, new MRI‐derived technologies are emerging that allow exploring different magnetic fields within a single system. One such technology is fast field‐cycling MRI, also known as FFC‐MRI or FFC imaging, which enables a modulation of the main magnetic field during an imaging sequence giving access to field‐dependent relaxation properties as a novel contrast mechanism. FFC imaging derives from MRI but uses radically different technologies to generate the main magnetic field, and both types of scanner offer different views on biological tissues.

Indeed, varying the main magnetic field within a defined range requires dedicated hardware and various approaches exist to realize FFC imaging systems. In the clinical field range, FFC imaging is implemented by means of a insert coil together with the superconducting magnet provided by a commercial MRI system for 1.5T,, , or 3T. This approach, also referred to as delta relaxation enhanced MR (dreMR), has auspicious applications for the detection and quantification of contrast agents with increased specificity and sensitivity., , , Several systems were also developed to access the endogenous dispersion of tissues in the low‐field regime., , , , Recently, a whole‐body FFC scanner approved for clinical imaging studies was reported, capable of reaching any field from 50 µT to 0.2T. Controlled variations of the magnetic field with this single resistive magnet design allow for multi‐field quantification over a wide range of field strength while retaining image quality down to ultra‐low fields. Pilot studies show promising potential for innovations in the imaging of osteoarthritis, sarcoma, or brain stroke, with potentially important applications in medicine as an in vivo assessment method of multi‐field and dispersion information.

Compared to conventional MRI, implementation of fast field‐cycling poses additional demands on power supplies, control electronics, magnet design, pulse sequences and image quality. Signal‐to‐noise ratio (SNR), therefore, is an important issue for FFC scanners to satisfy the latter. High fields benefit from an inherently high SNR as they rely on stable and homogeneous acquisition fields provided by superconducting magnets. Although the SNR is not a limiting factor for the individual images, contrast in dreMR is obtained by image subtraction and strongly depends on the dispersion of the contrast agent in use. The magnitude of the dreMR signal is rather small in comparison to the individual images (eg, about 2.5% in Ref. 12) and retaining sufficient high SNR may become an issue. Similarly, low‐field FFC imaging systems operate with acquisition fields of 0.2 T or less, which limits the SNR compared to conventional clinical fields due to its dependency to . Moreover, image quality deteriorates because of poor magnetic field homogeneity, field instabilities during operation and delays in the field ramps between different phases in the pulse sequence especially for ultra‐low evolution fields.

For all these reasons, both high‐ and low‐field FFC scanners may strongly benefit from SNR‐enhancing methods and a vast number of these have been developed for high field MRI images in the recent years., , , , , These can be divided into denoising and reconstruction‐based approaches. The former takes a series of noisy input images and tries to find a denoised solution by making use of a priori knowledge in the form of regularization. Regularization can use either spatial information, information from the acquired series, or a combination of both. Denoising approaches are generally simpler to implement and computation time is lower compared to reconstruction approaches, but they can not recover structures that were missed by the k‐space to image space transformation due to poor SNR. To this end, recent approaches rely on a constrained reconstruction process, incorporating the a priori knowledge in the image generation process., , , In the case of quantitative MRI, this approach can be taken one step further by including the non‐linear MRI signal model into the reconstruction process, thus, directly acting on the parameter maps of interest., , , This kind of fitting approach is known as model‐based reconstruction in the high field MRI regime. Most regularization strategies rely on some sort of sparsifying transform to separate image content from noise and artifacts. Commonly used transforms include finite differences based approaches, , and applications of the wavelet transformation., , The regularization functional used highly influences the appearance of the final reconstruction and should be chosen based on a priori knowledge about the given parameter map. It was shown that total generalized variation (TGV) priors are a superb choice for both image reconstruction and quantitative MRI, leading to high quality reconstruction results without the stair‐casing artifacts of total variation (TV). Similar to TV, TGV uses information from the image gradient in combination with the assumption that images typically consist of a few, discrete edges and, thus, fits in the concept of compressed sensing. Opposed to TV, patches between edges are not constrained to have a fixed value but can rather be linearly varying in the case of second order TGV (TGV2). To this end, stair‐casing artifacts can be avoided using TGV. Higher order TGV functionals allow for even higher degrees of freedom within the patches but are typically only needed in stereo imaging, such as red/green/blue (RGB).

These constrained reconstruction and fitting methods apply well to the estimation of maps. Standard quantification with inversion recovery sequences requires the acquisition of an image series with different inversion times, leading to high redundancy in the information collected that can be exploited by the regularization algorithm. FFC imaging adds an extra dimension to MRI by varying the magnetic field during the relaxation phase of the pulse sequence, thus providing an additional field series. This multi‐field data offers new possibilities to exploit information redundancies to improve the quantification process. These redundancies could be used by a model‐based reconstruction approach, incorporating the data from all measurements at different field strengths into one large optimization problem. Each field leads to an individual map that shares common information with the other fields, for example, most edges in the map should coincide. This information can be exploited by joining the individual regularization functionals in parametric dimension via a Frobenius norm. The Frobenius norm is the matrix equivalent to the ‐norm for vectors and links edge information in parametric dimension. Such an approach was shown to further improve the quality of the resulting parameter maps in the context of mapping from highly subsampled data.

Herein, we formulate the multi‐field FFC imaging parameter quantification as a non‐linear model‐based reconstruction problem with Frobenius type TGV2 regularization. With this formulation as a single optimization problem it is possible to exploit all the joint spatial information of the additional field dimension to stabilize the quantification process and hence enhance the image quality. The proposed method is evaluated on simulated numerical FFC imaging data as well as on in vivo datasets from two stroke patients and compared to Tikhonov regularized fits from individual fields, all fields combined, and regularization using the squared ‐norm of the gradient (H1‐regularization). The results show improved stability of the parameter quantification with excellent noise suppression properties. In particular, the proposed method reveals remarkable contrast between the lesion and surrounding tissues in case of ultra‐low fields.

THEORY

Fixing notation

Throughout the course of this work we fix the following notations. The image dimensions in 2D are denoted as and , defining the image space with defining a point at location . expresses the space of unknowns ‐maps, with and the number of evolution fields. The unknowns consist of the scaling , the correction factor for each field , and the field‐dependent relaxation times . The measured data space is denoted as and is constituted of measurements of 2D k‐spaces . Each measurement is a combination of an evolution field with an associated number of measurement time points , and the number of time points acquired at a specific evolution field . To simplify notation we will drop the indices and refer to as .

FFC imaging signal model

In the most general case of a FFC imaging pulse sequence, the main magnetic field is rapidly cycled between three different levels: polarization field , evolution field , and signal detection field . A designated pre‐polarization of the sample magnetization is not necessarily required in the high SNR regime and the polarization field can be set to the detection field, that is, . For simplicity, we will assume that this is the case, and we will refer to these fields as and the corresponding equilibrium magnetization as , respectively, as this does not alter the validity of our approach in the case of low‐field systems since the effect of polarization at a different field can be compensated by a polarization efficiency term that blends into the inversion efficiency parameter. A schematic of a typical inversion recovery FFC imaging pulse sequence can be seen in Figure 1: following an inversion RF pulse, the main magnetic field is cycled to the desired strength and the spin system undergoes a relaxation associated with the applied evolution field during a given evolution time . The longitudinal magnetization at the end of this evolution period is given bywhere and are the equilibrium magnetizations for the detection and evolution field, respectively. The relaxation time, corresponding to the evolution field applied, is given by and corrects for field‐dependent effects from ramping the field combined with non‐ideal inversion efficiency of the RF pulse. The equilibrium magnetization is proportional to the strength of the applied magnetic field, that is, and for the detection and evolution field, respectively. This can also be written as

FIGURE 1

Pulse sequence diagram of an inversion recovery FFC imaging pulse sequence. A pre‐polarization of the sample magnetization can be applied by cycling the main magnetic field to (dashed line) or by setting the polarization field to the detection field , that is, (solid line). After an inversion pulse at , the longitudinal magnetization evolves at the desired evolution field for a given evolution time. The following magnetization can be detected by any MRI acquisition module. Note that the MRI signal is both inverted and detected at . and represent the equilibrium magnetization for and , respectively

After the evolution period the signal is acquired at to ensure that the Larmor frequency of the spins corresponds to the tuning frequency of the receive RF coil. With Equations (1) and (2), and the unknown parameters , the acquired signal can be modeled for a specific evolution field and evolution time by the non‐linear signal equation , given bywith representing the Fourier transformation and sampling of k‐space.

Multi‐field parameter fitting

Acquiring several time points for a specific evolution field allows to quantify , , and . Typically, each field yields a different value; thus, the fitting process must be repeated for each evolution field, omitting joint information in the unknowns such as structural information. By combining the separate fitting steps into a single optimization problem it is possible to use the shared information to stabilize the quantification process. Furthermore, joint information between different parameter maps can be exploited by means of a Frobenius type functional. Such a fitting and regularization strategy has been successfully applied in other multi‐channel fitting problems., , Thus, we propose to apply a similar approach to quantify , , and from multi‐field FFC imaging data, using shared information, especially between the different maps.

Model‐based reconstruction framework

Assuming Gaussian noise corrupts the measurement data , it is possible to quantify the unknown parameter via a regularized non‐linear, minimum‐least‐squares problemwhich origins from a maximum a posteriori approach using Bayes' theory. denotes some non‐linear forward operator and reflects a priori knowledge about the unknowns by means of a regularization term. can be used to weight between data and regularization term and may either be a scalar value or a vector for each unknown in .

For multi‐field FFC data, is linked to the measurement data via We denote all measured data as and the corresponding mapping from unknown space to data space as . Thus, the optimization problem is defined as

is chosen as TGV2 regularization with a joint Frobenius norm on all unknowns . denotes the finite symmetric derivative and the auxiliary variable balances between the first and second derivative of the TGV2 functional. This type of regularization was shown to have favorable properties for multi‐parameter model‐based reconstruction.33 The terms resemble the Frobenius type TGV2 functionals, joining common spatial information of the unknown parameter maps by combining gradient information of all maps via an ‐norm in parameter direction.

Using the TGV2 model parameters and it is possible to balance the approximated first and second derivatives, avoiding the stair‐casing artifacts of TV while maintaining its favorable edge‐preserving features. The ratio of TGV2 is fixed throughout this work, as it was shown to yield good results for MRI image reconstruction. The numerical solution is done in analogy to Ref. 33 via a Gauss‐Newton (GN) approach. Leading to an inner GN problem of the form

The linearization is done via a Taylor series expansion of w.r.t. each unknown in at position . Constant terms are fused into to keep the notation clean. amounts to the Jacobian of , evaluated at . Introducing an additional weighted ‐norm penalty on improves convexity of the function. The weighting matrix can be used to resemble a Levenberg‐Marquardt update if . The regularization parameters and balance between the three terms and are reduced after each linearization step. Reducing the weights was shown to be beneficial in the context of the IRGN algorithm.

Using Fenchel duality the problem of non‐differentiability can be overcome and Equation (6) can be cast into a saddle‐point form

constitutes the linear operators encountered in Equation (6) within data and TGV2 norm, reflects the quadratic penalty on , and denotes the dual norms of the data and TGV2 term. Problems in such a form can be solved via a primal‐dual algorithm using a line‐search to speed up convergence. Mathematical details for each step are given in Supporting Information Text A. The update scheme written as pseudo code is given in Supporting Information Text B.

METHODS

Numerical FFC imaging data

To evaluate the proposed model‐based reconstruction approach numerical FFC imaging data were simulated using parameters measured from FFC imaging scans of brain stroke patients at the University of Aberdeen as part of a separate study (PUFFINS study, see details in section 3.2). The numerical phantom followed a schematic geometry and dispersive characteristics of an axial head scan with four regions (see Figure 3), representing the subcutaneous fat (region of interest [ROI] 1), the tissues surrounding the brain (ROI 2), the brain (ROI 3) and a stroke‐like lesion (ROI 4). values were simulated by means of a power‐law dispersion with model parameters and , , coarsely in line with proton ‐NMRD profiles of fat (ROI 1), white (ROI 2) and gray matter (ROI 3), and stroke lesions (ROI 4) measured in vivo from the PUFFINS patients cohort (data to be published). The values retained for the different evolution fields and times are summarized in Table 1 and Table 2, respectively, and served as a ground truth for the validation of the quantification. The numerical FFC imaging phantom was first generated as a vector graphic to be subsequently converted to matrix data to allow for any desired sampling resolution. In this case, we used an image resolution with matrix size of pixel, which is typical for the original FFC imaging of stroke patients. Tissue reference values were assigned for each ROI, and a data series was generated using the signal equation in Equation (3). Additionally, a small constant phase offset was introduced for each .

FIGURE 3

Multi‐field maps obtained from simulated FFC imaging inversion recovery data. The reference maps for three different evolution fields (200 mT, 21.1 mT, and 2.2 mT) are shown at the top. The different reconstruction methods are presented in each row. Standard refers to single field pixel‐wise fitting, multi‐field to combined field, pixel‐wise fitting approach and H1 to the model‐based approach with regularization using the squared ‐norm of the gradient. The proposed method is shown in the last row. The columns show increasing noise from left to right

TABLE 1

Evolution times and fields used to generate the simulated images and corresponding timings for the in vivo acquisitions

Application	Variable	Field (mT)	Evolution time (ms)
Simulation	tnevo_1	200	455	242	129	68	36
	tnevo_2	21.1	282	150	80	42	23
	tnevo_3	2.2	136	73	39	21	11
Patient I	tnevo_1	200	455	242	129	68	36
	tnevo_2	21.1	282	150	80	42	23
	tnevo_3	2.2	136	73	39	21	11
Patient II	tnevo_1	200	455	196	84	36
	tnevo_2	37	338	145	63	27
	tnevo_3	6.9	196	84	36	16
	tnevo_4	1.3	114	49	21	9

TABLE 2

Parameter values selected to generate the simulated images

Parameter	ROI	1	2	3	4
power law	a	5.6	4.4	2.6	3.8
	b	−0.1	−0.15	−0.3	−0.08
T_1(ms)	at 200 mT	152	178.5	237.3	231.4
	at 21 mT	121.3	127.3	120.7	193.2
	at 2.2 mT	96.8	90.8	61.3	161.3
α abs (a.u.)/phase (rad)	at 200 mT	1/0.5236
	at 21 mT	0.75/0.6981
	at 2.2 mT	0.6/0.8727
C(a.u.)	at all fields	1	1/3	2/3	2.03/3

Simulated proton density values were normalized to 1, resulting in a theoretical maximum signal amplitude of 1 for the simulated series. Zero‐mean Gaussian noise was added on both real and imaginary parts of the images to simulate the noise arising from the patient tissues and acquisition system. The noise amplitude was selected as a percentage of the theoretical maximum signal in the overall ground truth image series ranging from 1% to 4%, reflecting a typical range from the FFC images acquired. SNR directly after inversion at ranged from 33.3 to 8.3 in the white matter ROI and 66.7 to 16.7 for the gray matter ROI, respectively. Note that low‐field FFC images tend to exhibit markedly lower signal strength than higher‐field ones because of losses when the magnetic field switches, so their SNR was proportionally more affected using this approach. Finally, the image series was transformed to k‐space via a 2D Fourier transformation, as input for the proposed method with TGV2 and H1 regularization, respectively. The pixel‐wise fitting methods were applied to the image series.

In vivo FFC imaging data

The performance of the proposed method was tested on in vivo FFC imaging patient data. Two data sets obtained from patients scanned for a brain stroke were selected, as part of the PUFFINS study currently taking place at the University of Aberdeen. This study has been approved by the North of Scotland Research Ethics Committee (study number 16/NS/0136), and all the participants agreed for the clinical and FFC imaging data to be used anonymously for research purposes. The scans selected both present a lesion in the ultra‐low field regime that could not be easily observed at 200 mT, as illustrated in Figure 2 for patient I. Both cases were assessed from computed tomography (CT) and diffusion‐weighted MRI scans as embolic stroke for patient I and multiple embolic events for patient II. FFC measurements were performed using a whole‐body FFC scanner using a FFC inversion‐recovery spin echo sequence with an echo time (TE) of 24 ms, 20 kHz bandwidth, 8.37 MHz acquisition frequency, 10 mm slide thickness, and single slice acquisition. The images had a field of view of 290 mm and a resolution of 128 x 128 pixel in‐plane with 80 phase encode acquisitions and partial Fourier acquisition (80 lines out of 128). The sample was pre‐polarized at 200 mT for 300 ms before each evolution periods with the timings as shown in Table 1, for an acquisition time of 40 min.

FIGURE 2

FFC images of a stroke patient from the PUFFINS study (patient I). Image obtained at 200 mT evolution field show good signal after inversion (B) and after 455 ms evolution time (A) but low contrast in the lesion, while low‐field images at 21 mT (C) show good lesion contrast but most other tissue shows little to no signal

Data processing and corrections

The original raw image was reconstructed using partial Fourier completion to recover the correct image ratio. Phase‐encode artifacts were removed using a method previously published, but the images had not been filtered or further modified. While the noisy images were used as input for the standard pixel‐based fitting, the corresponding noisy k‐space data were used as input for the proposed fitting process using H1 and TGV2 regularization, respectively.

As reference method lsqnonlin of Matlab (The MathWorks, Inc.) was used for fitting Equation (4) field‐by‐field and pixel‐by‐pixel, where was replaced with Tikhonov regularization on the unknowns to stabilize fitting. Prior to fitting, images were smoothed in k‐space using the following filter functionwith denoting the cutoff radius, the k‐space location, and as parameter for the slope of the filter.

As a second reference method, all fields were combined for the pixel‐wise fitting without pre‐smoothing, similar to the proposed method. As a third reference, an H1 regularization was used as , that is, penalizing the squared ‐norm of the gradient of . The latter approach was implemented in Python and optimized using the proposed IRGN algorithm with an accelerated gradient descent optimizer for the inner iterations. Again, no pre‐smoothing was applied.

The analyses with the proposed method were done by implementing the FFC signal model in PyQMRI. All fittings were performed on a desktop PC equipped with an Intel(R) Core(TM) i7‐6700K CPU @ 4.00GHz with 64 gigabyte of RAM and a NVIDIA GeForce GTX 1080 Ti GPU with 12 gigabyte of RAM.

Optimization

The regularization weights and were reduced after each linearization step, following the iterative regularized Gauss‐Newton scheme. and were used as initial values and were reduced by a factor of and , respectively. To account for the typical smooth appearance of , corresponding regularization weights were multiplied by a factor of 10. The reduction steps were repeated down to and . In total, 12 linearization steps were performed. The number of primal‐dual iterations for each sub‐problem was doubled starting at 10 iterations up to 2000 iterations, that is, . If the relative decrease in the primal problem or the decrease of the primal‐dual gap was less than , the inner iteration was terminated. The step sizes of the primal‐dual algorithm used were determined via a line‐search, described in Algorithm 2. The same approach and the same weights have been used for the H1‐regularized reference method. Weights for the Tikhonov based approaches have been selected as small as possible to still achieve a stable fitting ().

RESULTS

The simulated high noise level can be seen as residual noise in the reconstructed maps of the pixel‐wise fitting approaches (Figure 3). Simultaneously, a difference to the simulated reference is visually noticeable in the pixel‐wise fitting approach. The H1 approach is able to reduce these outliers but suffers from blurring at image edges. The proposed model‐based method is able to reduce outliers throughout all noise levels and is visually closer to the simulated reference values. Plots of and in Supporting Information Figure S3‐2 show similar results. The single‐field pixel‐wise fitting approach even fails to capture the correct phase of the simulated phantom. A pixel‐wise relative absolute difference plot (Figure 4) confirms this visual impression of reduced noise using the proposed approach. The proposed method shows an up to 18‐fold lower mean error in the phantom, computed over all pixels, than standard pixel‐wise fitting. The error increases with increased noise level, as can be expected. Difference plots also reveal a slight bias of the proposed method. The bias of the methods is further assessed in 2D joint histogram plots (Figure 5). For these plots, values of all fields are combined to form a single plot. The proposed method shows slight underestimation of high values, as reflected by points lying below the identity line. However, noise could be greatly reduced compared to the pixel‐wise fitting, and different ranges are clearly separated and show a similar distribution as the simulated values. Fitting with the standard method took approximately 100 s. The proposed method took roughly 120 seconds.

FIGURE 4

Pixel‐wise relative absolute difference to the ground truth values. Numbers next to the difference images show mean relative absolute error within the phantom. All values are given in percent

FIGURE 5

2D histogram evaluation for maps in Figure 3, which were obtained from synthetic FFC imaging data by pixel‐wise fitting (standard), combined field pixel‐wise fitting (multi‐field), multi‐field H1, and joint model‐based reconstruction (proposed). The dashed line represents identity. Shown are reference values on the ordinate versus results obtained with the different reconstruction methods on the abscissa. Points below the identity line correspond to under‐estimation, points above to over‐estimation, respectively. All values are given in ms

The improvements in estimations held true when processing real FFC imaging data from stroke patients. The maps of unfiltered FFC images obtained using standard fitting‐based processing methods could not resolve anatomical features inside the brain region, as seen in Figures 6 and 7. Spatial regularization in combination with multi‐field fitting could greatly improve image contrast. The proposed method offers clear distinguishable structures in maps at 200 mT and is even able to recover some structural details in lower fields. It also assessed sharp features around the lesion area appearing at 37 mT and below in both patients. Fitting took approximately 65 and 150 s with the standard method for patient I and II, respectively, whereas the proposed method took 100 and 240 s for each patient, respectively.

FIGURE 6

In vivo multi‐field maps of a transverse slice of the brain of stroke patient I. From top to bottom, maps were obtained at three different evolution fields by pixel‐wise fitting of the signal model for each separately, combined field pixel‐wise fitting, multi‐field model‐based reconstruction with H1 regularization and by the proposed multi‐field model‐based reconstruction approach using the joint information of all three evolution fields (bottom row)

FIGURE 7

In vivo multi‐field maps of a transverse slice of the brain of stroke patient II. maps were obtained at four different evolution fields mT. The different reconstruction methods are given in each row. From top to bottom the methods are pixel‐wise fitting of the signal model for each separately, combined field pixel‐wise fitting, multi‐field model‐based reconstruction with H1 regularization and by the proposed multi‐field model‐based reconstruction approach using the joint information of all three evolution fields

The quality of the maps obtained allowed estimating the dispersion curves for different ROIs, as shown on Figure 8 for subcutaneous fat selected under the scalp, the area of the lesion observed at the lowest field strength, and white and gray matter as seen at the highest field strength (the ROIs are shown in Supporting Information Figure S4). The dispersion profiles of fatty tissues show large SDs, which may be attributed to the presence of various types of tissues within these ROIs, due to the relatively low resolution of the image. Otherwise, the dispersion profiles of white matter, gray matter, and the areas of the lesions are similar between the two patients. This is encouraging given that the two lesion have a similar diagnosis of ischemic stroke.

FIGURE 8

dispersion profiles obtained from the patients I (solid lines) and II (dashed lines) in the regions of subcutaneous fat (SC fat, orange) measured between the scalp and the brain, gray matter (GM, gray) measured over two centimeters of the cortical region, white matter (WM, blue) measured over the inner region of the lobes and the lesion (red). The error bars stand for twice the standard deviation of the values measured across the ROIs. The positioning of the ROIs is depicted in Supporting Information Figure S4

DISCUSSION

The approach used here has high potential to serve as a new standard procedure for fast post‐processing of FFC MRI data. As the phantom simulations showed, the noise in the reconstructed maps could be reduced very efficiently while preserving important anatomical details to a high extent. The algorithm outperforms established methods based on pixel‐wise fitting of the relaxation profiles yielding lower deviations from the reference values and significantly less variance (Figures 3, 4, 5). The improved stability results from the combination of information from all acquired fields and exploiting the existence of similar topological structures in the different unknowns. The improved stability is also reflected by increased accuracy of recovered pseudo proton density and correction factor values (Supporting Information Figures S1‐S3). Higher deviations of larger values in the reference methods are due to the Tikhonov regularization used,which penalizes the larger values than lower. Also the maps show significantly reduced variance, although there remains some bias that may be, at least partially, due to residual errors in . Another remarkable feature of the multi‐field methods is their ability to accurately recover the phase information in and , making phase correction prior to fitting obsolete. This in turn can improve maps, as no normalization with a noise phase estimate is necessary.

The advantages of the improved fitting approach become immanent in the in vivo applications (Figures 6 and 7). The standard approaches based on pixel‐wise fitting fail to reconstruct image details in both patients. In current practice, k‐space windowing filters are applied to recover usable information, but this dramatically reduces image resolution by filtering out the high‐frequency components of the image, which are responsible for the sharp features. In contrast, the joint regularization approach can recover clearly distinguishable gray and white matter regions at 200 mT on the two patient datasets, previously hidden in noise. The values obtained for the different regions of interest agree well between the patients, given the estimation of the error provided by the variation of the values within each ROI (Figure 8). The values were systematically higher in patient I than in patient II, which may be attributed to patient variability and different RF receive coil sensitivity relative to the used ROIs. In addition, lesion localization agrees well with conventional MRI and CT based imaging, shown in Supporting Information Figure S5.

As expected from the raw images, the largest contrast for stroke appeared below 0.1T, where values were larger than that of the surrounding tissues. A cutoff appears between 30 and 100 mT (or equivalently 1.2 to 4.2 MHz) above which the contrast disappears. This is consistent with the fact that higher clinical fields do not show significant changes in acute ischemic stroke. Clearer interpretations may be provided from the analysis of the full data set but a tentative explanations of this phenomenon could be made by taking into account the biological effects of ischemia. During acute ischemia neuron cells swell and burst, and this process is likely to disorganise large structures that interact with water over timescales that correspond to the cutoff frequency observed, that is, between 0.2 and 0.7 . The degradation of these components of the brain structure could have the effect to reduce the efficiency of the relaxation pathways at low magnetic fields, as observed here. Another possibility could be that the reduction of water mobility through the membranes of neurons may decrease the contribution of the intracellular water relaxation to the overall signal, which may dominate at low field but could be less efficient at higher fields. These explanation would be consistent with the absence of contrast at higher fields.

As in the phantom images the model‐based and spatially regularized methods proved to preserve anatomical features with high spatial frequencies because of using the existence of sharp edges for regularization. These approaches are increasingly accurate with the number of views that can be compared showing the same object, either as a repetition of a recording or as different acquisition of the same field of view, as it is the case here. Hence, FFC imaging can benefit from the high information redundancy obtained from the typical acquisition method, which repeats the measurement of the field of view at different evolution times and fields.

As the proposed approach is model‐based, and can therefore provide directly from the raw images, it could be used to reduce the number of steps required to process the image and limits data losses. However model‐based approaches also limit the amount of information that is extracted from the image, and properties not covered by the signal equation, may be missed. For instance, brain tissues are known to follow bi‐exponential relaxation because of the presence of the myelin sheath around the axons. Hence in a subsequent step, the model will be adapted to the type of scan, or following a test for potential multi‐exponential behavior.

Using direct reconstruction from k‐space opens up the possibility of undersampled image acquisition while maintaining high quality in the reconstructed maps. The proposed method allows for different kinds of undersampling and is not limited to Cartesian sampling or single slice acquisitions. While a single receive coil hast been used for the current study, the extension to a multi‐coil setup is straight forward. The combination of multiple receive coils and the potential of undersampling k‐space could be used to reduce acquisition time in FFC imaging which shall be subject of a future study. The gained acquisition time might lead to a clinical acceptable scan time using the three or four fields shown in this work or could be spent to investigate a multitude of different field strength. However, such extensions would require modifications to the phase correction algorithm, which is based on images.

Another advantage of direct reconstruction from k‐space data is the validity of the Gaussian noise assumption in the real and complex parts of k‐space. In the typically used magnitude images, noise is non‐linearly transformed, resulting in a Rician or non‐central Chi distribution. This invalidates the basic assumptions used to derive the ‐norm data fidelity term and can lead to a bias in the final solution. Even though the data term can be modified to account for these variations, the modified version need not be convex or differentiable. Thus, optimization of the correct function might lead to suboptimal solutions or demanding optimization algorithms. In practice, the favorable properties of the ‐norm usually outweigh the drawback of the bias to the theoretically optimal solution; thus, it is widely used.

A potential limitation of the proposed approach is the risk of cross‐contamination of the information between images due to the joint regularization., , It is assumed that features share the same edge position. If this assumption is violated in one parameter map, artificial edges might be introduced. The likelihood strongly depends on the used norm for joining the information. As we use a relative weak coupling by means of a Frobenius norm, such cross‐contamination is unlikely. It was shown in previous work that Frobenius norm based joint regularization does not show cross‐contamination in practice., , It might only occur if way too strong regularization weights are used; however, such cases would be discarded in practice as images would look unnatural.

The proposed reconstruction and fitting approach is integrated into a recently published Python framework for quantitative MRI. This framework allows for an easy adaption to different signal models and, thus, a broad application of the proposed method. Adaptions to the signal model can be made by simply editing text files. In addition, 3D regularization strategies are possible, which were shown to further improve reconstruction quality.,

CONCLUSIONS

We have successfully introduced joint TGV2 regularization to multi‐field quantification from FFC imaging. The highly significant improvements in estimation makes it now possible to obtain clinically usable multi‐field maps, and to produce reliable and comparable results. This shows exciting potential for the exploration of low magnetic fields and dispersion effects as illustrated here on two stroke patients.

FIGURE S1 Absolute value of multi‐field maps obtained from simulated FFC imaging inversion recovery data. The reference maps for three different evolution fields (200 mT, 21.1 mT and 2.2 mT) are shown at the top. The different reconstruction methods are presented in each row. Standard refers to single field pixel‐wise fitting, multi‐field to combined field, pixel‐wise fitting approach and H1 to the model‐based approach with regularization using the squared ‐norm of the gradient. The proposed method is shown in the last row. The columns show increasing noise from left to right. Values next to each figure represent the mean value within the simulated phantom in a.u.

FIGURE S2 Phase of multi‐field maps obtained from simulated FFC imaging inversion recovery data. The reference phase maps for three different evolution fields (200 mT, 21.1 mT and 2.2 mT) are shown at the top. 0The different reconstruction methods are presented in each row. Standard refers to single field pixel‐wise fitting, multi‐field to combined field, pixel‐wise fitting approach and H1 to the model‐based approach with regularization using the squared ‐norm of the gradient. The proposed method is shown in the last row. The columns show increasing noise from left to right. Values next to each figure represent the mean value within the simulated phantom in radiant

FIGURE S3 Absolute value and phase of maps obtained from simulated FFC imaging inversion recovery data. The reference map is shown at the top. In the left column, the multi‐field . maps were obtained by pixel‐wise fitting of each evolution field separately (standard), and the right column results from joint model‐based reconstruction of all three evolution fields together (proposed). The noise level increases from top to bottom from 1% to 4%

FIGURE S4 Regions of interest selected to extract the dispersion profiles in Figure 8 in patients I (left) and II (right). The regions for white matter are delineated in light blue dashed lines, grey matter in solid dark blue lines, fat in yellow dotted lines and lesions in red dot‐dashed lines

FIGURE S5 Images obtained from clinical examinations of Patient I (images on the left) and II (images on the right) from CT (top) and conventional MRI using 3D ‐FLAIR (centre) and DWI (bottom), showing only the slice that corresponds with the single slice acquisition of the FFC acquisition. The ‐FLAIR (‐weighted‐Fluid‐Attenuated Inversion Recovery) sequence had an isotropic resolution of 0.625 mm, a SPIR fat suppression, an inversion delay of 1650 ms, an echo time of 340 ms and a repetition time of 4800 ms. The DWI (Diffusion‐Weighted Imaging) had an in‐plane resolution of 0.8 mm, a slice thickness of 4 mm, a STIR fat suppression, an echo time of 77 ms, a repetition time of 3478 ms and a b‐factor of 1000. The results indicate ischemic strokes and are clearly visible in both MRI FLAIR and DWI images, with patient II exhibiting multiple small strokes. CT scans are less informative for ischemic stroke, as illustrated here on patient II

Click here for additional data file.