Non-linear fitting with joint spatial regularization in Arterial Spin Labeling

Multi-Delay single-shot arterial spin labeling (ASL) imaging provides accurate cerebral blood flow (CBF) and, in addition, arterial transit time (ATT) maps but the inherent low SNR can be challenging. Especially standard fitting using non-linear least squares often fails in regions with poor SNR, resulting in noisy estimates of the quantitative maps. State-of-the-art fitting techniques improve the SNR by incorporating prior knowledge in the estimation process which typically leads to spatial blurring. To this end, we propose a new estimation method with a joint spatial total generalized variation regularization on CBF and ATT. This joint regularization approach utilizes shared spatial features across maps to enhance sharpness and simultaneously improves noise suppression in the final estimates. The proposed method is evaluated at three levels, first on synthetic phantom data including pathologies, followed by in vivo acquisitions of healthy volunteers, and finally on patient data following an ischemic stroke. The quantitative estimates are compared to two reference methods, non-linear least squares fitting and a state-of-the-art ASL quantification algorithm based on Bayesian inference. The proposed joint regularization approach outperforms the reference implementations, substantially increasing the SNR in CBF and ATT while maintaining sharpness and quantitative accuracy in the estimates.


Introduction
Arterial spin labeling (ASL) is a non-invasive MRI technique for quantifying local tissue perfusion (Detre et al., 1992).
The method utilizes magnetically labeled blood water by invertkristian.bredies@uni-graz.at (Kristian Bredies), rudolf.stollberger@tugraz.at (Rudolf Stollberger) ing the blood water spins upstream the imaging region. After waiting a specific period of time, called the post labeling delay (PLD) which accounts for the time the magnetically labeled blood needs to flow into the region of interest, an image is acquired. This so called label image is subtracted from a second image, the control image, acquired without magnetization alterations of the inflowing blood. From this difference image, also known as perfusion weighted image (PWI), the cerebral blood flow (CBF) can be quantified using a general kinetic model (Buxton et al., 1998). The recommended clinical ASL protocol (Alsop et al., 2014;Telischak et al., 2014) consists of single-delay pseudo-continuous ASL (pCASL) (Dai et al., 2008) combined with segmented 3D data acquisitions such as gradient and spin echo (GRASE) (Feinberg et al., 2009;Günther et al., 2005) or turbo spin echo (TSE) stack of spirals (SoSP) (Ye et al., 2000;Dai et al., 2008) readout due to efficient background suppression (Ye et al., 2000) and SNR gains (Alsop et al., 2014) of these methods. For single-delay acquisitions the signal in the PWI depends on both, the CBF and arterial transit time (ATT). Hence, the accuracy of CBF estimation from single-PLD ASL data is dependent on both factors.
Another limitation is that for cases with prolonged ATT (ATT > PLD) some of the labeled blood may remain in larger vessels.
This leads to bright spots and an overestimation in CBF in bigger vessels and an underestimation of CBF in brain areas. The bright spots are known as vascular artifacts and can complicate clinical diagnosis in patients with stroke, steno-occlusions, or moya-mayo disease (Zaharchuk, 2012). A way to improve the clinical interpretation of single-delay ASL images is by applying additional coefficient of variation maps obtained from the multiple PWIs (Mutsaerts et al., 2017). Another way to reduce misquantification is by using a longer PLD, ensuring that the blood has sufficient time to reach the tissue. However, this leads to longer acquisitions and additionally to a lower SNR due to the T1-relaxation of the labeled blood. Alternatively, multi-PLDs can be used to sample the inflowing blood at several time points, from a short PLD to long PLD. By fitting the acquired signal to a kinetic model, the potential bias in CBF due to unknown ATT can be reduced (Buxton et al., 1998).
However, the recommended segmented acquisitions have the drawback of a low temporal resolution with increased sensitivity to inter-segment motion. Therefore, only a limited number of PLDs can be acquired in a clinically acceptable time. Recently, accelerated single-shot 3D acquisition strategies (Ivanov et al., 2017;Boland et al., 2018;Spann et al., 2019) were implemented to overcome this drawback, at the cost of reduced SNR.
This makes the estimation of reasonable quantitative ATT and CBF maps from this low SNR perfusion weighted time series challenging. The standard voxel-wise non-linear least squares (NLLS) fitting approach leads to outliers in low-SNR voxels.
To this end, a weighted delay approach  proposed to reduce outliers in the quantitative maps. Further improvements could be achieved by inclusion of spatial priors on the CBF map  in a Bayesian inference model (BASIL (Chappell et al., 2010)). This stabilizes the fitting approach and reduces noise, ultimately leading to improved CBF estimates but introduces spatial blurring. Exploiting all available spatial information by means of joining the individual regularization of each unknown into a single, joint regularization functional can further improve reconstruction quality. Such an approach has been successfully applied in the context of relaxometry (Knoll et al., 2017;Wang et al., 2017;Maier et al., 2018). Joint regularization utilizes information present in each map, such as tissue boundaries, by means of advanced spatial regularization functionals to avoid the loss of small features and promotes overall sharper parameter maps. In this study, we propose a new non-linear fitting algorithm with joint spatial constraints on the CBF and ATT map to stabilize the estimation procedure and hence enhance the image quality. To improve the motion robustness of the 3D acquisition, we combine the proposed method with a single shot CAIPIRINHA accelerated 3D GRASE readout. The method is evaluated on synthetic phan-tom datasets including simulated pathologies, on six healthy subjects, as well as on seven stroke patients and compared to NLLS and BASIL without regularization on ATT (BASIL w/o) and with regularization on CBF and ATT (BASIL w/ ).

Fixing notation
Throughout the work we fix the following notations. The image dimensions in 3D are denoted as N i , N j , N k , defining the image space U = R N i ×N j ×N k with x = (i, j, k) defining a point at location (i, j, k) ∈ N 3 . u ∈ U N u expresses the space of unknown CBF-and ATT-maps with N u = 2 in this case. The measured data space is denoted as D = C N i ×N j ×N k and consists of N d perfusion weighted images derived from Control/Label (C/L)-pairs, measured at time t = (t 1 , t 2 , . . . , t N d ) ∈ R N d + .

Parameter fitting
From a statistical point of view the problem of identifying the unknown parameters u = (u 1 , u 2 , . . . , u N u ) ∈ U N u given a series of noisy measurements d = (d 1 , d 2 , . . . , d N d ) ∈ D N d can be solved via maximum likelihood estimation. Assuming the measurements at time t n are generated by some function A φ,t n : u → d n with fixed parameters φ the likelihood function of measuring d is given by p(d|u, A φ,t n ). The realization of p depends on the noise distribution in the measurements d. Under the assumption that additive independent and identically distributed zero-mean Gaussian noise with variance σ 2 (AWGN) corrupts the measurements d, the multivariate likelihood function turns into a product of single-variate functions. It is common to minimize the negative logarithm of the likelihood function, which is equivalent to maximizing the likelihood, as it turns the product into a sum and improves the numerical stability. Omitting constant terms with respect to u yields u * ∈ arg min u∈U Nu (1) which resembles the well known minimum least squares problem with · 2 being the standard L 2 -norm.
Typically, several measurements with varying sequence parameters are necessary to quantify tissue parameters. Especially in cases with a non-linear relationship between acquired signal and parameters, fitting is performed in an iterative fashion.

The ASL signal model
The quantification of CBF and ATT is based on the standard model for pseudo-continuous ASL (pCASL) (Buxton et al., 1998) which reads as where u = ( f, ∆) and f amounts to CBF in ml/g/s, but is normally quoted in ml/100g/min, and ∆ to ATT in seconds. The a priori known parameters of equation 2 are combined into the variable φ = (M 0α , T 1 , T 1b , τ). It is assumed that T 1 , the apparent longitudinal relaxation decay constant of the tissue, amounts to 1.33 seconds at 3T. T 1b is the longitudinal relaxation decay constant of blood, assumed to amount to 1.65 seconds at 3T . τ corresponds to the labeling duration, α is the labeling efficiency and set to 0.7 (Dai et al., 2008) and t n is the acquisition time point, i.e. the sum of post labeling delay and labeling duration, for the n th measurement. Further, the bloodbrain partition coefficient λ is assumed to be 0.9 ml/g (Herscovitch and Raichle, 1985) thus 1/T 1 app ( f ) = 1/T 1 + f /λ, and M 0α = αM 0 /λ with M 0 being the acquired proton density weighted image.

Regularization
As the acquired PWI images suffer from poor SNR the problem of quantifying CBF and ATT typically suffers from numerically instabilities. A method to incorporate a priori knowledge of the parameters u into the maximum likelihood estimation problem 1 is known as maximum a posteriori estimation and leads to with γ > 0 balancing between the data fidelity term and the regularization R. R(u) includes known information about u such as its statistical distribution or spatial features, e.g. u should consist of piece-wise constant areas. As the variance σ 2 is in general unknown, we will not consider σ 2 fixed but as something that can be chosen in the reconstruction process. Thus, we combine it with the regularization parameter γ. The introduced prior can lead to a biased estimate of u with reduced uncertainties (Brinkmann et al., 2017). Thus a trade-off between faithfulness to acquired data and the prior needs to be determined according to the expected noise in the data. The most basic form consists of classical Tikhonov regularization which penalizes outliers in the parameter maps in an L 2 -norm sense (Tikhonov and Arsenin, 1977). An extension to this basic form consists of penalizing the gradient of the maps which is known as H 1 regularization (Tikhonov and Arsenin, 1977), leading to a smoother appearance but comes at the cost of blurred edges. To preserve edges and to obtain a better visual impression, a sparsity promoting functional is usually preferred which can be realized by posing an L 1 -norm based constraint on the sparse domain of the unknowns (Donoho, 2006;Lustig et al., 2007). As u is usually not sparse in its native domain, a sparsifying transform such as a finite differences operation or a wavelet transformation is used. The total variation (TV) functional of Rudin-Osher-Fatemi (ROF) (Rudin et al., 1992) is based on an L 1 -norm combined with a forward finite differences operator. This combination can be interpreted as a spatial piece-wise constant prior which is known to be prone to stair-casing artifacts in the final reconstruction results . In order to avoid these stair-casing artifacts but leverage the edge-preserving feature of TV a generalization termed total generalized variation (TGV) functional was proposed by Bredies et al. (2010). In the context of MRI, TGV 2 , which enforces piece-wise linear solutions by balancing between a first order and approximated second order derivative, was shown to yield excellent reconstruction results, preserving fine details and edges while maintaining the denoising properties of TV (Knoll et al., 2010). In the discretized form the TGV 2 regularization is realized via a minimization problem of the following form The favorable properties of TGV 2 can be further improved by sharing common feature information between the unknown parameter maps by joining the TGV 2 functionals utilizing a Frobenius norm in parametric dimension (Bredies, 2014). Recently, this combination was shown to yield improved reconstruction results compared to separate regularization on each map in the context of quantitative T 1 mapping (Maier et al., 2018) and multi modal image reconstruction (Knoll et al., 2017). The combination by means of a Frobenius norm is justified by the assumption that quantitative maps share the same features at the same spatial positions. To incorporate the Frobenius norm the following adaptations to the TGV 2 semi-norm with v = (v 1,l , v 2,l , v 3,l ) N u l=1 ∈ U 3×N u constituting the approximation of 3D spatial derivatives, and for the symmetrized gradient χ = (χ 1,l , χ 2,l , χ 3,l , χ 4,l , χ 5,l , χ 6,l ) N u l=1 ∈ U 6×N u 2.5. The non-linear, non-smooth optimization problem The combination of TGV 2 with equation 3 leads to which is a non-linear problem in the unknowns u and nonsmooth due to the L 1 -norms of the TGV 2 functional. Recall that for the ASL signal, the non-linear operator A φ,t n (u) is defined by equation 2 and u amounts to u = ( f, ∆). A similar problem arises in model-based quantification of T 1 and M 0 (Roeloffs et al., 2016;Wang et al., 2017;Maier et al., 2018). The problem is thus solved in analogy via a two-step procedure. First the data fidelity term is linearized in a Gauss-Newton (GN) fashion, second the linearized, non-smooth sub-problem is solved using a primal-dual splitting algorithm. The linearized sub-problem for each linearization step k is given by Constant terms stemming from the linearization at position u k are fused with the data byd k n = d n − A φ,t n (u k ) + DA φ,t n u k and the matrix DA φ,t n | u=u k = ∂A φ,tn ∂u (u k ), i.e. the derivative of the signal with respect to each unknown, can be precomputed in each linearization step. The additional weighted L 2 -norm penalty 2 improves convexity of the function and resembles a Levenberg-Marquadt update if the weight matrix M is chosen as M k = diag(DA φ,t n | T u=u k DA φ,t n | u=u k ) or a simpler Levenberg update if M k is chosen as identity matrix. It was shown by Salzo and Villa (2012) that the GN approach converges with linear rate to a critical point for non-convex problems with non-differential penalty functions if the initialization is sufficiently close. By exploiting the Fenchel duality it is possible to transform the problem in equation 8 into a saddle-point which overcomes the non-differentiability issue of the L 1 -terms.
Problems of the form in (9) can be efficiently solved using a first order primal-dual splitting algorithm (Chambolle and Pock, 2010) in combination with a line search (Malitsky and Pock, 2018) to improve the convergence speed. The detailed derivation is given in the Appendix A. Pseudo-code for the implementation can be found in Appendix B.

Reference Methods
For comparison of the proposed algorithm we used the nonlinear least squares (NLLS) as well as the Bayesian Inference for Arterial Spin Labeling MRI (BASIL) method Groves et al., 2009 (Smith et al., 2004;Woolrich et al., 2009;Jenkinson et al., 2012) and uses Bayesian inference to estimate the unknown parameter maps. It incorporates fixed non-spatial priors as well as adaptive non-local spatial smoothing priors for the parameters. The spatial smoothing prior is used for CBF and is directly based on evidence in the data. The smoothing strength is adjusted based on the local support in the specific area in the data. The arterial (macro-vascular) contribution flag was set to "OFF" in BASIL to facilitate comparability to the proposed method which currently implements the pCASL model omitting the local arterial contribution. In addition to this standard form of BASIL, termed BASIL w/o, a simple duplication of the line associated with spatial priors in the starting script enables priors on both, CBF and ATT, as described in Groves et al., 2009). This modification serves as second BASIL reference and is termed BASIL w/.

Phantom generation
To evaluate the proposed method, synthetic ASL data was generated from brain T 1 and PD maps supplied by MRiLab (Liu et al., 2017)

Error propagation and stability
To asses the error propagation and stability due to the nonlinear fitting procedure we performed a pseudo replica analysis for all three methods. To this end, 100 different noise realizations with a standard deviation of 0.65 were simulated for Case 3. Due to the non-linear fitting process a Gaussian noise assumption in the parameter maps could be violated, thus the median and inter-quartile range between the 25 th and 75 th quartile were used for evaluation. We assessed potential biases using the medians of differences to the ground truth of the 100 realisations and compared the differences in the IQRs between the methods. For the synthetic dataset GM and WM binary masks, generated on the ground truth phantom, are employed.
Based on the down sampled GM and WM mask, low resolution mask were generated by thresholding the corresponding GM/WM masks with 0.7. Mask for the simulated lesions were generated in analogy.
Additionally, we compared the estimated CBF and ATT maps of the proposed method with the results of BASIL and NLLS by means of a relative difference to the numerical ground truth parameter maps.
To assess if differences in median or IQR are statistically  Bonferroni (1936).
Each method and tissue was considered as parallel test case.
Median and IQR of CBF and ATT were considered as separate cases. Results were considered statistically significant for p-values less than 0.05 (p<0.05).

In vivo measurements
All measurements were performed on a 3T MAGNETOM Prisma (Siemens Healthcare, Erlangen, Germany) system using a 20-channel head coil. Written informed consent was obtained by all healthy volunteers as well as by all patients following the local ethics committee's regulations. In total, six healthy volunteers, consisting of five male and one female subject with an age of 29.5±2.6 years were analyzed.Additionally, seven patients with ischemic stroke due to middle cerebral artery occlusion who received successful re-canalization therapy (i.e. intravenous thrombolysis followed by mechanical thrombectomy), consisting of six male and one female subject with an age of 57.1±13 years, were considered. Patient data was acquired 24 hours after recanalization therapy. ASL images were acquired using a prototype 3D pCASL sequence with a 2D CAIPIRINHA accelerated single-shot 3D GRASE readout (Ivanov et al., 2017) and two background suppression pulses (Vidorreta et al., 2013). Labeling efficiency for this sequence was experimentally determined in Vidorreta et al. (2013)  Additionally, for each healthy subject a T 1 weighted image was acquired using a 3D-MPRAGE sequence with the following imaging parameters: 1 mm 3 isotropic resolution, 176 slices, TR = 1900 ms, TE = 2.7 ms, TI = 900 ms, flip angle = 9 • , acquisition time = 5 min 58 s.

ASL Data Processing
The accelerated ASL images were reconstructed directly on the scanner console by means of a prototype reconstruction pipeline provided by the vendor. The reconstructed ASL images were motion corrected using Statistical Parameter Mapping (SPM)12 1 (Wellcome Trust Centre for Neuroimaging, University College London, UK) (Friston et al., 2007) and the ASL-Toolbox (Wang et al., 2008;Wang, 2012). This rigid-body based motion correction process involved three sub-steps as described in Wang (2012). After reconstruction and motion correction the perfusion weighted time series were calculated.
From this perfusion weighted time-series the CBF and ATT maps were estimated using the proposed method as well as the two reference methods. The fixed parameters φ amount to the same values as in the synthetic data set except for T 1 =1330 ms, the approximate tissue T 1 relaxation constant.

Anatomical Image Processing
For each healthy subject brain masks and PV estimates for GM and WM were computed using FSL (FMRIB Software Library, Oxford, UK (Jenkinson et al., 2012)) and BASIL.
In a first step, non-brain tissue was removed from the high resolution structural (T 1 weighted) image using the FSL tool BET (Smith, 2002). In a second step, PV estimates for GM and WM were obtained from the T1w image using the FSL tool FAST (Zhang et al., 2001). Third, the structural image and brain mask were registered to the mean ASL image using the FSL tool FLIRT with 6 degrees of freedom (Jenkinson and Smith, 2001;Jenkinson et al., 2002). The obtained transformation matrix served as initial guess for the next registration refinement step, implemented in BASIL. This step used the epi reg tool of FSL for boundary based registration of the perfusion image with the segmented white matter mask (Greve and Fischl, 2009). In the last step, the PV estimates for GM and WM were transformed to the ASL-space by a process that integrates over the volume of the low resolution voxels as described by Chappell et al. (2011) and implemented in the FSL tool applywrap.
Finally GM and WM binary masks in ASL space were computed by thresholding the PV estimates at 70% in WM and GM respectively. For the patient data brain masks were generated from the M 0 image using the FSL tool BET due to the missing T 1 weighted image.

Method Comparison
Healthy subjects were compared based on visual inspection of ATT and CBF for 1 and 4 acquired averages. In addition, WM and GM masks were used to compute median and IQR which were visualized with box-plots. Statistically significant differences in median and IQR between methods were assessed using Mann-Whitney-U tests, similar to the ones in the numerical simulation. p-values were adjusted for multiple comparisons (Bonferroni, 1936). Results were considered statistically significant for p-values less than 0.05 (p<0.05).
Stroke patients are compared based on visual inspection only.

Parameter optimization
To identify a good set of model and regularization parameters a grid search was performed on the synthetic dataset and in vivo healthy subjects. The resulting regularization parameters amounted to γ init = 10 −3 and δ init = 1 which were reduced by 0.5 and 0.1 respectively after each Gauss-Newton step down to γ f inal = 6.5 · 10 −6 and δ f inal = 10 −2 . A reduction of regularization parameters was observed to be beneficial for overall convergence in IRGN methods (Bakushinsky and Kokurin, 2004;Kaltenbacher et al., 2008;Kaltenbacher and Hofmann, 2010). Relative tolerance for convergence was set to 10 −8 between consecutive evaluations of function value. Regarding the inner iterations, 50 were used in the initial Gauss-Newton step and the number was increased by a factor of two until the maximum allowed number of 1000 iterations is reached , i.e.

Implementation
The proposed method is implemented in Python 3.7

Results
We compare the fitting quality of our proposed joint spatial TGV 2 regularization strategy to established quantification on the CBF-and ATT-maps at three levels:  (table 1). However, for the estimated ATT-map both methods show similar relative difference. Using spatial priors on CBF and ATT (BASIL w/ ) reduces this variations but also seems to introduce a slightly lower value in ATT (  . Case 1 has no simulated pathologies, Case 2 shows hyperperfusion in CBF only and Case 3 hyperperfusion in CBF and increased ATT in the corresponding area. In the first column the numerical ground truth is shown and in the following columns the estimated CBF-and ATT-maps from NLLS, BASIL without regularization on ATT (BASIL w/o), and with regularization on ATT (BASIL w/), and the proposed method without and with regularization on ATT, respectively. The proposed method with regularization on both unknown maps shows improved noise removal in CBF and ATT compared to the other methods due to joint spatial constraints. Median and 25% -75% IQR for selected ROIs is given in table 1 Fig. 2. CBF and ATT maps of synthetic phantoms of three cases with pathologies are shown. Each case represents a large partly occlusion of the arteria media combined with a small partly occlusion in frontal gray matter. No variation in CBF is simulated but each case shows an increase in ATT, which gets more severe from Case 4 to Case 6. The order of reference and displayed reconstruction algorithms is the same as in figure 1. The proposed method shows the least influence of the highly increased ATT on the CBF estimates and is able to recover higher ATT values in the affected areas than the other methods, as can be seen in table 1. Fig. 3. Pixel-wise relative absolute difference between the ground truth numerical reference and the quantitative maps, estimated with the algorithms of figure 1. The NLLS shows the greatest deviation in ATT and CBF. BASIL w/o reduces the relative difference in the CBF due to the spatial prior and BASIL w/ is able to reduce deviations even further. The proposed w/o method shows similar results on CBF and ATT as BASIL w/o. The least relative difference is achieved with the proposed method due to joint spatial constraints on CBF and ATT simultaneously.    ception of BASIL w/, which shows lower IQR than the proposed method. In GM CBF, the proposed method shows statistically significantly lower IQR than NLLS but both BASIL approaches are able to further reduce IQR than the proposed method. IQR of GM ATT shows the least deviations in the proposed method.
In the WM lesion the proposed method reduces IQR compared to NLLS in CBF and ATT statistically significant but no statistically significant difference to both BASIL methods is observed.
The CBF GM lesion shows statistically significant reduction of IQR using the proposed method over BASIL but no difference to NLLS. In the corresponding ATT, no statistically significant differences of IQR are observable. All p-values and median IQRs are reported in table 2.

Discussion
In this study we present a novel joint spatial regularization technique for quantitative ASL imaging, combining non-linear fitting with a TGV 2 functional. The proposed method poses a joint spatial TGV 2 prior on both, CBF and ATT, to improve the robustness of the fitting procedure. Synthetic ASL datasets with different pathologies as well as in vivo data from healthy and stroke patients with different SNR levels were considered.  Imposing prior knowledge on the unknown parameter maps leads to the fundamental problem of bias/variance trade-off, as the solution will depend on the used prior information. This is also true for the used MAP-based approach shown in this work.
The amount of bias, however, can be controlled by the used prior and the weight between data and prior information in optimization, respectively (Brinkmann et al., 2017). To this end, the NLLS approach without any regularization can be consid-  However, an extension to the complex model with the proposed method could be easily obtained by an adaptation of the signal equation used for fitting, which will be done in a future step.
As the proposed approach has been implemented into a Python toolbox (Maier et al., 2020), addition of new models can be achieved in a straight forward manner. Extension to other ASL models, e.g. pulsed ASL (PASL), is simply done by replacing the forward model in equation 2 with the appropriate one.
Simple models, not consisting of composed functions, can be included using a plain text file. Complex models need to be im-plemented in Python by the user but templates exists to help in the implementation process. A detailed description of the employed software and how to include new models can be found in Maier et al. (2020). An exemplary PASL fit for phantom Case 0 is given in Supplementary Material Figure S8.
In The contrary can be observed, the joint regularization produces the most stable CBF estimates of all methods. A totally wrong choice of the regularization weight compared to the supporting data could nevertheless introduce such errors. However, such a strong weight for regularization would also lead to a severely hampered visual impression and such fits would likely be discarded.  ASL imaging is very sensitive to signal variations from motion or changes in blood velocity due to the cardiac cycle (Verbree and van Osch, 2017). Currently, the proposed method does not directly account for these variations. Motion related variations are corrected for in a preprocessing step but no correction for blood velocity changes is applied. As it is planned to extend the method to use raw k-space data, motion could be included in the forward model, e.g. based on determined motion fields prior to fitting, as it is done in MOCHA. As estimation of motion directly from highly undersampled k-space data can be challenging, a robust estimation and correction needs to be found. Another possibility would be to include a motion term into the fitting process but this poses a mathematically challenging problem, especially for forming forward and adjoint operation pairs.
If ASL data is evaluated over large ROIs, NLLS seems to be the favourable method as it shows the least bias in our simulations (figure 7) for most investigated ROIs, especially in GM. In addition to the 3D acquisition used in this work, ASL is often performed in a 2D slice-by-slice fashion. Such data can also be fitted with the presented reconstruction framework, either, by applying regularization in 2D only or by adapting the third gradient direction to account for non-isotropic voxels. This adaptation amounts to a simple scaling of the gradient with respect to the ratio between in-plane and the acquired slice resolution, taking into account slice thickness and inter slice gap. Setting this scaling to zero equals 2D regularization. However, 3D regularization outperforms 2D, as has been shown in previous work (Huber et al., 2019;Maier et al., 2018).
Reconstruction on the GPU required roughly 1 GB of memory which should be available on any recent GPU. Computation speed varies with hardware and further reduction could be expected with the recent increase in GPU performance. If memory requirements supersede the available GPU memory, a double buffering strategy is available, as introduced in Maier et al.
(2020). New scanner consoles already often include GPUs, thus the proposed method could be directly integrated into the scanner reconstruction process.