Weierstraß-Institut für Angewandte Analysis und Stochastik

A BSTRACT . We present an implementation of a recently developed noise reduction algorithm for dMRI data, called multi-shell position orientation adaptive smoothing (msPOAS), as a toolbox for SPM. The method intrinsically adapts to the structures of different size and shape in dMRI and hence avoids blurring typically observed in non-adaptive smoothing. We give examples for the usage of the toolbox and explain the determination of experiment-dependent parameters for an optimal performance of msPOAS.


INTRODUCTION
Diffusion-weighted magnetic resonance imaging (dMRI) has developed into an extremely versatile tool for the in-vivo structural analysis of tissue, for example in the human brain [Johansen-Berg and Behrens, 2009].One reason is, that the diffusion signal obtained with the pulsed gradient spin sequence echo [PGSE, Stejskal and Tanner, 1965] directly relates, via three-dimensional Fourier transform, to the diffusion propagator which is the probability density of the underlying Random Walk the spin particles experience [Mitra and Sen, 1992].Therefore, if we measured the diffusion signal for all possible diffusion gradient directions, times and strengths, i.e. cover the whole q-space, we would know the full propagator.
Its spatial and directional dependence would allow us to infer on boundaries for the diffusing particles and hence the underlying structure.However, in practice, only a limited coverage of the q-space is feasible.
Therefore, a number of models have been developed in the past, which reveal at least partial information contained in the diffusion propagator.Most require dMRI data measured on at least one q-shell, that is characterized by a single b-value [Basser et al., 1994a] subsuming diffusion gradient strength and diffusion time.
The most prominent example of a diffusion model gives rise to Diffusion Tensor Imaging [DTI, Basser et al., 1994a,b].Surprisingly, although this model actually describes free diffusion in anisotropic media it has proven to relate well to the underlying tissue geometry in the brain in general, and to the main fiber directions in the white matter in particular [Johansen-Berg and Behrens, 2009].More sophisticated descriptions of the diffusion signal have been examined to infer on more complicated sub-voxel structure like multiple fiber directions.These include HARDI [Frank, 2001], Q-ball imaging [Tuch et al., 2002], tensor mixture models [Behrens et al., 2003, Assaf and Basser, 2005, Tabelow et al., 2012], higher order tensor approximations [Özarslan and Mareci, 2003, Liu et al., 2003, Jensen et al., 2005] and methods to determine the full diffusion propagator [Özarslan et al., 2006[Özarslan et al., , Cheng et al., 2010]], see Assemlal et al. [2011] for a recent review.
In any case the measures of interest, like fiber directions or quantitative measures like the fractional anisotropy (FA) in DTI, are estimated based on the raw diffusion images.Hence, the accuracy of the estimates depends on the data quality, which in turn typically requires retrospective correction of artifacts due to eddy-currents, motion [Mohammadi et al., 2010], susceptibility-related distortions [Ruthotto et al., 2012[Ruthotto et al., , 2013]], instrumental [Mohammadi et al., 2012a,b] or physiological noise [Mohammadi et al., 2013b,a].Signal-to-noise ratio (SNR) is especially low in dMRI because of the additional exponential dependence in the diffusion-weighted signal, see e.g.Stejskal and Tanner [1965].As a result the SNR in dMRI decreases with increased diffusion weighting, i.e. higher q-shells have lower SNR.Furthermore, beyond-tensor models with a larger number of parameters are more sensitive to data noise, making sophisticated denoising strategies mandatory for cutting-edge dMRI.In order to reduce noise in dMRI data a number of different approaches have been developed in recent years starting from Gaussian filtering [Westin et al., 1999], smoothing procedures in tensor space for DTI [Fletcher, 2004], to denoising algorithms based on partial differential equations [Ding et al., 2005, Parker et al., 2000, Duits and Franken, 2011] to name only a few.
Recently, we developed a position-orientation adaptive smoothing algorithm [POAS, Becker et al., 2012] based on the propagation-separation (PS) approach [Polzehl andSpokoiny, 2006, Becker andMathé, 2013].The method is edge-preserving and avoids blurring of the fine anisotropic structures observed in dMRI.The method has been extended to be applicable to dMRI measured on multiple shells and named multi-shell POAS [msPOAS, Becker et al., 2013].It capitalizes on the additional information on the different shells.Furthermore, several improvements make msPOAS feasible from a computational point of view.Finally, msPOAS can also be applied to single-shell dMRI data.Thus, in this paper we will only consider the more general msPOAS approach.
Compared to previous diffusion-model-based adaptive smoothing methods [e.g., Tabelow et al., 2008] (ms)POAS has the advantage that it directly denoises dMRI data without requiring any diffusion-model assumptions.Thus, it avoids a bias towards any diffusion model such as DTI, HARDI, or tensor mixture models (see Becker et al. [2012] for details).
(Ms)POAS has been originally implemented in the R language and environment for statistical computing and graphics [R Development Core Team, 2013].While the corresponding package dti [Tabelow and Polzehl, 2013] is easy to install [Polzehl and Tabelow, 2011] and allows for applying msPOAS to dMRI data using very few code lines, the use of R in the neuroimaging community is still limited.We therefore implemented msPOAS as a toolbox for Statistical Parametric Mapping (SPM) [Friston et al., 2006], the most widely used neuroimaging analysis package, to make the method available to a broader audience, see www.diffusiontools.com.
In this paper we shortly review the method in a simplified way, present the new toolbox for SPM, describe the usage of the toolbox, and suggest methods to determine experiment-dependent measures and to set method parameters.We present some worked examples with single-and multi-shell data.For a more in-depth review of the theory of POAS and msPOAS we refer the reader to the original work in Becker et al. [2012] and Becker et al. [2013].

METHODS
MsPOAS is a noise reduction method for dMRI data that is measured on at least one q-shell, i.e., for at least one b-value and a sufficiently large number of diffusion gradient directions.Although msPOAS also works for very few diffusion gradient directions, it particularly benefits from the information from more gradient directions [Becker et al., 2013].MsPOAS is also suitable for single-shell dMRI data.
MsPOAS in a nutshell.MsPOAS is based on the propagation-separation approach for point-wise estimation of piecewise smooth functions [Polzehl and Spokoiny, 2006], implemented as estimation of a piecewise constant function, to accelerate and stabilize the procedure.MsPOAS is built upon the assumption that the measurement space of voxel position ( v) and diffusion gradient orientation ( g) can be partitioned into regions with almost constant image intensity values separated by discontinuities, which is motivated by the piecewise continuous nature of the neuroanatomy.MsPOAS describes the multi-shell dMRI signal as a vector-valued function S (see next section) defined on the measurement space.In msPOAS a weighted mean of the observations is iteratively optimized to infer on the expected signal at v and g.The weights are composed as a product of two factors: (a) The first factor, which is common to non-adaptive kernel estimators (e.g. the Gaussian filter), defines the region of interest covered by the smoothing kernel (defined by a bandwidth h), (b) The second factor is crucial for introducing and defining adaptivity, because it characterizes the differences of diffusion weighted signals at the measurement points.The criteria for the adaptive factor are based on statistical information about the signal and noise in the MR signal.To improve the statistical sensitivity, the adaptive factor uses information from all shells (thus multi-shell POAS).Furthermore, it makes use of the estimates from the previous step of the iteration, successively improving the discrimination of the homogeneity regions.During iteration the local neighborhood increases while the adaptation becomes more restrictive, avoiding blurring at the edge of tissue structures.The rigorous derivation of msPOAS requires a detailed mathematical treatment, see Becker and Mathé [2013], Becker et al. [2013].
Review: Multi-shell position-orientation adaptive smoothing (msPOAS).The design space in diffusion MRI for a single q-shell forms an orientation space [Duits and Franken, 2011].Points m = ( v, g) in this space, denoted by R 3 × S 2 , are given by its voxel position v and gradient orientation g, see Becker et al. [2012].Denote S b (m) the measured signal value at voxel v, gradient direction g, and b-value.The non-diffusion weighted b 0 -image is denoted by S 0 ( v, 0) = S 0 ( v) and does of course not depend on g.In case of several b 0 -images we consider their mean.For each v and g we arrange all data (including the b 0 -image) as a vector of length B + 1, where B is the number of q-shells.As the diffusion directions do not necessarily coincide on all shells, this vector representation of the data may require some interpolation of values on a shell.For details, see Becker et al. [2013].
The measured random values S b (m) are distributed according to some probability distribution that is parametrized by a "true" intensity parameter θ, the noise standard deviation σ and a number of degrees of freedom 2L where L denote the effective number of coils for parallel imaging [Aja-Fernández et al., 2009].Then S b (m)/σ is assumed to be non-central χ-distributed with 2L degrees of freedom.
In msPOAS we assume that similar signal values in R 3 × S 2 extend over sets of neighboring points n = ( v n , g n ).This can be used to obtain an improved estimate for the image value at any point m ∈ R 3 × S 2 .If the definition of the neighborhood is specific for the point (and based on the data) we call the method adaptive, or non-adaptive otherwise.The notion of neighborhood typically requires the definition of a metric in the considered design space, here R 3 × S 2 .
MsPOAS is derived from non-adaptive kernel estimators [Nadaraya, 1964, Watson, 1964], see also [Fan and Gijbels, 1996].Then, for a given distance δ(m, n) between design points m and n, a non-adaptive kernel estimate S b (m) for the expected value at some bandwidth h is given by a weighted mean where K loc is some kernel function.Instead of the commonly used Gaussian kernel function we employ (with (x) + denoting the maximum of x and 0), due to its higher efficiency and computational simplicity, see [Fan and Gijbels, 1996].
MsPOAS is based on the propagation-separation approach [Polzehl andSpokoiny, 2006, Becker andMathé, 2013] and makes the following important extensions to the non-adaptive estimator: 1 Re-define the weighting schemes given in Eq. ( 2) by an additional term that evaluates the distance of the signal in two measurement points m and n making the weights adaptive, see Eq. ( 4), cf. the construction of bilateral filters. 2 Repeat the estimation step in Eq. (1) using adaptive weights (4) for an increasing (typically geometric) sequence of bandwidths h k for k = 1, ..., k instead of a single bandwidth h, see Eq. ( 2), in an iterative procedure, see Becker et al. [2013].This approach refers to scale-space methods.
MsPOAS then calculates at each iteration step k new estimates with adaptive weighting schemes mn employing a second kernel function In msPOAS we use a simple approximate distance in the design space R 3 × S 2 : (5) where ||.|| denotes the L 2 -norm in R 3 and κ steers the influence of the geodesic distance on the sphere S 2 [Hagmann et al., 2006, Becker et al., 2013].Specifically we use a decreasing sequence κ k = κ 0 /h k to restrict the smoothing on the sphere at later iteration steps, see Becker et al. [2013].
The statistical penalty s mn evaluates the difference between the estimators in points m and n in the previous iteration step and is defined by where KL denotes the Kullback-Leibler divergence between two non-central χ-distribution with 2L degrees of freedom and its expectation given as argument.This term has no analytic expression, in msPOAS we use the approximation ( 6) characterizes the variance reduction at each shell achieved due to averaging.For all specific details of msPOAS which are not covered here necessary for the handling of the b 0 -images, the modification for interpolated signal values and the initialization of the method, we refer to Becker et al. [2013].
Toolbox implementation and installation.The toolbox POAS4SPM for the neuroimaging software SPM has been implemented using C and Matlab.It comes as open-source software with GPL2.
The toolbox is part of the ACID toolbox for "Artefact correction in diffusion MRI" and can be downloaded from its homepage at http://www.diffusiontools.com/.It is listed on the SPM extension homepage http://www.fil.ion.ucl.ac.uk/spm/ext/, too.Installation is done by extracting the toolbox into the toolbox folder of SPM, and compiling all mex-files in the Preprocessing/POAS subfolder.Running the make_ACID.mutility in the cfiles folder of the ACID toolbox will automatically compile all necessary c-files.
Usage of the toolbox.As part of the ACID-Toolbox, the msPOAS module runs in the batch editor of SPM.In the menu of the batch editor msPOAS can be found at SPM -> ACID Toolbox -> Pre-Processing -> Choose POAS options -> POAS.One may load/save a batch file to use standard toolbox settings as usual.The toolbox options are defined as follows: DTI images: Choose the N images including N g diffusion weighted and N 0 non-diffusion weighted data files.
Diffusion directions: Add a 3 × N -array consisting of the diffusion gradient directions with normalized vectors that appear in the same order as the DTI images were entered.Choose a vector with three zeros k star: This is the parameter k of msPOAS that defines the number of iterations and thus the maximal location bandwidth h k .kappa: This is the parameter κ 0 of msPOAS that defines the initial ratio of the spatial and spherical distance in Eq. ( 5).lambda: This is the adaptation parameter λ of msPOAS, see Becker et al. [2012] and Becker et al. [2013] for more details.
sigma: The value σ is the noise level in the data and must be obtained from the data, see option Estimate sigma in the toolbox.Although the value of σ may vary spatially due to effects of parallel imaging, the current implementation of msPOAS assumes a homogeneous variance.The effect of a misspecified σ can be partly compensated by the choice of λ [Becker et al., 2013].
ncoils: This parameter specifies the parallel imaging factor L, i.e., the number of different receiver coils that contributed to the measured signal value.It may also vary with spatial location, but the current implementation of msPOAS assumes a global value for L. msPOAS has been shown to be relatively robust against misspecifications of L [ Becker et al., 2013].
After running the batch script, the smoothed diffusion weighted volumes are written to disk using poas as a prefix.Only one b0-image, obtained as smoothed average of all original b0-images, is written to disk, see Becker et al. [2012].For further processing the corresponding gradient orientations and b-values are written to disk as a mat-file.The input directory is used as the target directory for the script's output.Choice of the method parameters k , κ 0 , λ.The number of iteration steps k defines the maximum variance reduction in homogenous image regions, but also the numerical complexity of the method and thus the computation time, see Fig. 1.Large values of k may also lead to a step function approximation with a small step size.We therefore suggest a value between 10 and 12 for k .
The value for κ 0 defines the amount of initial smoothing on the sphere of diffusion weighting directions in the first step k = 0.This spherical smoothing stabilizes the initial estimates of msPOAS for a better noise reduction especially at very low SNR.On the other hand spherical smoothing comes with a cost of potential bias and loss of spherical resolution.In successive steps msPOAS therefore restricts the amount of spherical smoothing by the specific choice of the sequence κ k .Due to the decreasing variance of the estimates in the iteration process the statistical penalty in Eq. ( 6) generally increases.This leads to more restrictive adaptive weights and thus to reduced smoothing on the sphere.We recommend a value for κ 0 such that in the initialization step k = 0 for msPOAS an average number of 5 to 10 signal values on the sphere for neighboring gradient directions are included in the smoothing.This is equivalent to a choice cos −1 (1 − 5/ N ) < κ 0 < cos −1 (1 − 10/ N ), where N = N g /B denotes the mean number of measured gradient directions per shell [Becker et al., 2013].Larger values of κ 0 , are appropriate in case of small SNR or large number N g of gradient directions.For a graphical visualization of this choice and the specific values for the datasets used in this paper, see Fig. 2. λ is the adaptation bandwidth of the procedure that steers the amount of adaptation of msPOAS.For λ = 0 the adaptive weights in Eq. ( 4) vanish for all m = n.As a consequence, all estimates at any iteration step will coincide with the original data.In contrast, for λ = ∞ the adaptive weights coincide with the non-adaptive weights, i.e. w(k) mn = w mn for all m, n, see Eq. ( 4), and msPOAS finally generates a non-adaptive kernel estimate in the space R 3 × S 2 with kernel K loc , bandwidth h k and κ k = κ 0 /h k .
Basically, λ can be chosen to satisfy a propagation condition using simulation independent from the processed data [Becker et al., 2012] but for the specific noise distribution.In case of dMRI data, the distribution can be assumed to be a non-central χ-distribution with 2L degrees of freedom and noncentrality parameter θ.We suggest a value of λ = 12 for all datasets.Depending on the quality of the estimates for the data-dependent parameters σ, L slight adjustments may be required, as discussed below.How to estimate σ and L? MsPOAS requires two data-dependent input parameters σ and L, that fix the properties of the noise distribution and enter the definition of the statistical penalty s mn , see Eq. ( 6).
A suitable estimate for the standard deviation σ of the noise can be obtained by any method available to the user, see Aja-Fernández et al. [2009] for a review.The methods typically determine σ from the Rayleigh L = 1 or central χ-distribution L > 1 in the image background.
We implemented the method described in Constantinides et al. [1997] in the toolbox.It can be accessed in the batch editor via SPM -> ACID Toolbox -> Pre-Processing -> -> Choose POAS options -> Sigma estimation and requires as input the diffusion weighted images from the data, a binary mask file defining background voxel only and the specification of L, see below.
Running the estimation returns a mean value over all diffusion weighted images.It also writes a .txt-fileinto the data directory with individual values for each diffusion weighted volume.b 0 -image typically lead to larger estimates for σ, than the diffusion weighted volumes.Note, that the method of Constantinides et al. [1997] does not account for noise correlation due to multi-channel receiver coils [Hutton et al., 2012].Furthermore, non-background structure like ghosts within the defined background mask also leads to overestimates of σ.For msPOAS it is advisable to use a conservative small estimate for σ and potentially correct a misspecification by adjusting λ.
The parameter L depends on the reconstruction algorithm, see e.g.Aja-Fernández et al. [2011].It is very difficult to estimate from the data.For some reconstruction methods it can be shown that L = 1 [Sotiropoulos et al., 2013].Generally L equals the total number of receiver coils for a Sum-of-Squares reconstruction and is spatially varying for GRAPPA [Aja-Fernández et al., 2011].Fortunately, msPOAS has been shown to be relatively robust against misspecification of L [ Becker et al., 2013].We therefore suggest to use a value of L = 1 consistently in the estimation of σ and for msPOAS, if no other estimate is available.
Experimental data.Two healthy volunteers (male) participated in the study approved by the local ethics committee after giving written informed consent.The example data used in this paper has been acquired as follows: Experiments were performed on a MAGNETOM Trio, a Tim System 3T scanner (Siemens AG, Healthcare Sector, Erlangen, Germany).Two high-resolution diffusion magnetic resonance imaging (dMRI) data sets were acquired using a reduced field-of-view (FoV) technique [Heidemann et al., 2010], one multi-shell, and one single-shell data set.
For the multi-shell data, the 161 × 58mm FoV was centered on the motor cortex.It had 1.2mm isotropic resolution, and 10 percent slice gap, resulting in an effective slice thickness of 1.3mm.The images were acquired at 3 different b-values: 21 at b = 20s/mm 2 , 100 at b = 800s/mm 2 , and 100 at b = 2000s/mm 2 using the directions suggested by Caruyer et al. [2011].The total scan time was about 22min.This dataset was also used in Becker et al. [2013].
The 156 × 56mm FoV of the single-shell data set was centered on the thalamus with 1mm isotropic resolution, and 10 percent slice gap, resulting in an effective slice thickness of 1.1mm.The images were acquired at 2 different b-values: 5 at b = 0s/mm 2 and 64 at b = 1000s/mm 2 using the directions provided by Siemens.This dataset was acquired three times, giving a total scan time of about 20min.
Prior to POAS the data were corrected for motion and eddy current artifacts using the method detailed in Mohammadi et al. [2010], which is implemented as part of the ACID toolbox pipeline.For the analysis in this paper we then estimated the diffusion tensor and FA.
Parameter choices for the multi-shell dataset.We repeatedly defined an arbitrary region within the background of the data and used the method implemented in the toolbox and described in the Methods section to estimate the noise standard deviation σ.We consistently found a value of σ = 30.We used L = 1 for all calculations as in Becker et al. [2013].
The number of iteration steps was fixed at k = 12 that provided a good balance between computational costs and achieved noise reduction.
Parameter choices for the single-shell dataset.Estimation of the noise level in the image background for the single-shell dataset consistently provided a value of σ = 45.We used L = 1, k = 12, and λ = 12 for all calculations.This dataset was measured with three repetitions.First, we smoothed the data for a one-repetition dataset (using the first repetition), here the κ 0 value was used to be κ 0 = 0.5.
Then, we smoothed the data for all three repetitions, with κ 0 = 0.3.

RESULTS
The calculations on the described hardware using the optimal parameters given above took approximately 1200sec = 20min for the multi-shell dataset, 700sec ≈ 12min for the one-repetition singleshell dataset, and 1900sec ≈ 30min for the three repetition single-shell dataset.
In Figure 3 we compare an axial slice of the multi-shell data for the three used b-values 0, 800, 2000s/mm 2 before and after application of msPOAS.The gray scale has been adjusted such that the b 0 -image uses the full range of the scale.For the diffusion weighted images the scale is comparable among them but augmented in comparison with the b 0 -image for better contrast.
In Figure 4 we show one (randomly selected) diffusion weighted image before and after applying msPOAS, see Figs 4(a) and 4(b).For comparison we show an average of all three repetitions of the In Figure 5 we show the corresponding results evaluated in the diffusion tensor model by means of color-coded FA maps.
In Figures 6, 7, and 8 we show the dependence of the resulting color-coded FA map on the number of iteration steps k (k = 4, 12, 20, 28 only), the adaptation bandwidth λ (λ = 1, 5, 12, 100, ∞ only), and κ 0 , respectively.The amount of noise reduction increases with the number of iteration steps.At some stage no more improvements are achieved by increasing k .λ controls the amount of adaptation from full adaptation for very small λ, i.e. no smoothing at all, to non-adaptive smoothing for very large values.κ 0 has an influence on the achieved noise reduction via construction, smaller values lead to less noise reduction.

DISCUSSION
We implemented a method for adaptive denoising diffusion weighted MRI called msPOAS as a toolbox for SPM.The program is part of a general toolbox for artifact correction in diffusion MRI data named ACID.Dependence of the results on the method parameters and its interaction.MsPOAS requires the specification of data-dependent quantities σ and L and method parameters h , λ, and κ 0 .In this paper, we analyzed the dependence of the outcome of msPOAS on the choice of the method parameters.As shown in Figure 6 the effect of noise reduction increases with the number of iteration steps k .
In homogeneous compartments of the image, the smoothness increases with k .On the other hand, the adaptivity of the procedure obviously avoids blurring at the borders.Thus in principle increasing k improves the noise reduction effect of msPOAS without blurring structural border.This can be interpreted as an intrinsic local stopping criterion for the method.However, the computational costs of an increased k are rather high, see Fig. 1, such that a suitable compromise between computation time and required noise reduction has to be made, e.g.k = 12.
The adaptation bandwidth λ controls the adaptivity of the method, ranging from complete adaptation at λ = 0, where the original is not changed at all by msPOAS, to a non-adaptive estimate at λ = ∞, where the adaptivity of msPOAS is turned off, see Fig. 7. Best results can be achieved at λ = 12.
Different choices of λ have a similar effect on the results of msPOAS as the estimates for σ.Generally, if σ is underestimated, msPOAS will be to restrictive and only little smoothing effect will appear in the result.If σ is overestimated, this may lead to blurring in the result.This also means that the effect a misspecification of σ can be compensated (to some extent) by adjusting λ accordingly.
The choice of κ 0 influences the amount of smoothing on the sphere.For relatively good SNR κ 0 can be chosen smaller then the recommended value to reduce an estimation bias due to the violation of the a local constant signal function on the sphere.For lower SNR the initial estimates benefit from larger values of κ 0 through stabilization.Our construction of the sequences of bandwidths h k and κ k , see Eq. ( 5) and its discussion, automatically lead to increased noise reduction by larger values of κ 0 (with higher computational costs).
The choice of L has only a minor effect on the result of msPOAS, see Becker et al. [2013].Nevertheless, msPOAS may benefit from precise specification of L, if available.If L is unknown, we recommend to use L = 1.
Suggestions for parameter choices.We suggest the following procedure for the parameter choices for msPOAS: 2 Determine σ by some suitable method (e.g. the method given in the toolbox) using L as previously chosen.
7 Adjust parameters: slightly decrease λ if oversmoothing at borders occurs, which looks like Figs. 7(d) or 7(e).Slightly increase λ if the noise reduction in homogeneous regions is less then expected for the utilized k : In this case increasing k does not increase noise reduction in homogeneous regions.
8 Adjust k if more or less noise reduction for homogeneous regions is required.
9 Re-run msPOAS if adjustments are necessary.
The evaluation of the msPOAS result can be done at the level of diffusion weighted images or for diffusion model parameters, like FA maps.

CONCLUSION
MsPOAS is a powerful method for adaptive noise reduction in diffusion MRI data that is now available as a toolbox for SPM.We demonstrated and discussed the effect of different method parameters and data-dependent quantities on the results of msPOAS and gave recommendations for their choice and determination, respectively.

INFORMATION SHARING STATEMENT
POAS4SPM is part of the ACID-toolbox available at http://www.diffusiontools.com.

,
with sd L (x) denoting the standard deviation of a non-central χ-distribution with 2L degrees of freedom and expectation x.The factor Ñ (k−1)m,b for each b 0 -image.b-values:Add a 1 × N -array with b-values.They should appear in the same order as the DTI images were entered.The b-value is given in units of s/mm 2 .b 0 -images should have b = 0.If the data contains images with a small b-value (b < 100s/mm 2 ), which serve as reference image without directional information, mark them by using b = 0 as well.The diffusion-weighted images corresponding to different shells will be identified by their b-value.

FIGURE 1 .
FIGURE 1. Computation time in hours as a function of the number of iteration steps k for the single-shell dataset used in this paper and a typical set of parameters κ 0 = 0.5, λ = 12.The highlighted k = 12 provides a good balance between achieved variance reduction and computation time.

FIGURE 2 .
FIGURE 2. The choice of κ 0 depends on the mean number N of diffusion gradient directions per shell and a balance between achieved noise reduction and bias on the sphere.Choices for the datasets in this paper are indicated.
FIGURE 3. One slice of a diffusion weighted image w/o msPOAS for multi-shell data at all measured b-values.
FIGURE 4. One slice of a diffusion weighted image w/o msPOAS for the one repetition single-shell data, compared with the mean of the three-repetition data.For the same slice a high-resolution quantitative R1-image was depicted in Figure4(e) for anatomical reference.The R1 image was acquired with a multi-parameter protocol[Dick et al., 2012, Lutti et al., 2013, Sereno et al., 2013].

FIGURE 8 .
FIGURE 6. Dependence of FA after msPOAS on k for the one repetition single-shell data.