Home | About | Journals | Submit | Contact Us | Français |

**|**HHS Author Manuscripts**|**PMC2911978

Formats

Article sections

- Abstract
- I. Introduction
- II. The K-Bayes Model
- III. Simulation Study
- IV. Spin-echo data
- V. PRESS box data with single voxel spectroscopy validation
- VI. Summary and Discussion
- References

Authors

Related links

IEEE Trans Med Imaging. Author manuscript; available in PMC 2010 July 29.

Published in final edited form as:

Published online 2010 March 18. doi: 10.1109/TMI.2009.2037956

PMCID: PMC2911978

NIHMSID: NIHMS173694

Address Correspondence to: John Kornak University of California, San Francisco, Department of Radiology and Biomedical Imaging, Center for Molecular and Functional Imaging, China Basin Landing, 185 Berry Street, Suite 350, San Francisco, CA 94107 Tel: 415-353-4740 Fax: 415-353-9423 ; Email: ude.fscu@kanrok.nhoj

The publisher's final edited version of this article is available at IEEE Trans Med Imaging

See other articles in PMC that cite the published article.

A *k*-space-time Bayesian statistical reconstruction method (K-Bayes) is proposed for the reconstruction of metabolite images of the brain from proton (^{1}H) Magnetic Resonance Spectroscopic Imaging (MRSI) data. K-Bayes performs full spectral fitting of the data while incorporating structural (anatomical) spatial information through the prior distribution. K-Bayes provides increased spatial resolution over conventional discrete Fourier transform (DFT) based methods by incorporating structural information from higher-resolution co-registered and segmented structural Magnetic Resonance Images. The structural information is incorporated via a Markov random field (MRF) model that allows for differential levels of expected smoothness in metabolite levels within homogeneous tissue regions and across tissue boundaries. By further combining the structural prior model with a *k*-space-time MRSI signal and noise model (for a specific set of metabolites and based on knowledge from prior spectral simulations of metabolite signals) the impact of artifacts generated by low-resolution sampling is also reduced. The posterior mode estimates are used to define the metabolite map reconstructions, obtained via a generalized Expectation-Maximization algorithm. K-Bayes was tested using simulated and real MRSI datasets consisting of sets of *k*-space-time-series (the recorded free induction decays). The results demonstrated that K-Bayes provided qualitative and quantitative improvement over DFT methods.

Magnetic resonance (MR) techniques provide a means for mapping a number of physical and chemical parameters in the human body. In particular, magnetic resonance spectroscopic imaging (MRSI) [1, 2] combines the spatial sampling methods of structural magnetic resonance imaging (structural MRI) with sampling of a spectroscopic data dimension, thereby enabling discrimination of multiple resonance frequencies and mapping of multiple tissue metabolites. MRSI information is of diagnostic interest because alterations of metabolite concentrations can provide a more sensitive indication of disease than observing tissue structure via structural MRI alone.

The spatial resolution of MR techniques differs greatly, from approximately 1 mm^{3} for structural MRI (which images anatomy through measurement of variation in water concentration) to 1000-2000 mm^{3} for MRSI. This is because concentrations of metabolites are approximately a factor of 10,000 times smaller than that of water, leading to lower signal-to-noise ratio (SNR) for MRSI than structural MRI. The lower SNR must be partially offset by obtaining data at much lower spatial resolution, i.e. by sampling at fewer spatial frequencies, usually the center region of *k*-space (or spatial frequency-space).

Proton (^{1}H) MRSI datasets consist of a set of (discrete) time series with typical lengths on the order of 512 timepoints. Each series is recorded at an individual point in *k*-space, (*k _{x}, k_{y}*), where the various

The most widely used methods for reconstructing MRSI datasets focus on 4D DFT and spectral fitting or 3D spatial DFT followed by time domain fitting [3-6]. Conventional use of DFT for reconstruction of low-resolution MRSI data produces low-resolution images and artifacts associated with the limited sampling of spatial frequencies, such as Gibbs ringing, aliasing and blurring. Moreover, DFT has no mechanism for regularization via the incorporation of prior information, e.g. noise is transformed directly from *k*-space to image space. While a number of techniques have been suggested to address these issues individually (see e.g., Chapter 8 of Liang and Lauterbur [7]), a unified solution would be of great utility.

A number of general methods for improving spatial resolution within MR modalities have been proposed; however none of these has the dual advantage of combining a low-resolution *k*-space-time data generation model with prior information obtained from high-resolution structural (anatomical) scanning methods. Reconstruction to higher resolution using Bayesian image analysis methods for multi-contrast structural MRI is presented in [8]; however, they do not incorporate a full signal model and hence are subject to all of the problems with artifacts that are incurred from use of the DFT for reconstruction of low-resolution modalities. A model that incorporates a raw signal model for structural MRI was presented for maximum likelihood estimation in [9] and Bayesian maximum a posteriori reconstruction in [10]; however, the first model's focus was on high-resolution structural MRI, and the only prior information considered in [10] was a notion of overall *a priori* smoothness. Denney and Reeves [11] have developed a Bayesian approach for MRSI reconstruction that models the data in *k*-space-time and utilizes an edge preserving prior. They do not directly model the raw *k*-space-time series data of MRSI but instead use spectrally integrated estimates of the metabolites at each *k*-space point; for the paper they only considered the relatively strong water and lipid signals. Important uncertainty information is discarded when ignoring the original data and working with spectrally integrated estimates directly. Furthermore, the edge-preserving prior does not utilize tissue class information that could be valuable when determining effects specific to certain tissue types.

A related body of literature exists for the Bayesian reconstruction of emission tomography data using a co-registered and segmented structural MRI as prior information [12-14]. The idea is to model the distribution of counts at the receptors from individual voxels and combine this with structural prior spatial models of varying types, including models that incorporate line site detection for determining the tissue boundaries [12]; models that classify into a number of tissue regions within the emission tomography map [13]; and tissue composition models similar to the prior spatial model developed herein for MRSI [14].

A number of papers [15-20] aim to use tissue classification information to improve the reconstruction (and associated resolution) of MRSI data outside the Bayesian framework. The idea is to utilize basis function in order to smoothly represent signals within homogeneous tissue areas. Because it needs to be kept sparse (either by limiting the size of the basis set or by penalizing selection from a larger basis set), the basis function representation restricts the space of potential reconstructions. The earliest basis sets used were piece-wise constant functions such as those used in Spectral Localization by Imaging (SLIM) [15] and in Spectral Localization with Optimal Pointspread function (SLOOP) techniques, [16]. The SLOOP method extends on SLIM by optimizing the set of phase encoding gradients to match the pointspread function of the volumes of interest. The SLIM piece-wise constant basis function approach was also generalized to a wide range of basis function sets through the generalized series formulation [17]. A related Bayesian method is proposed by Haldar et al [21] that combines the basis set approach with prior structural information. The idea is to obtain *k*-space data at high-resolution, accepting significant loss of SNR in the data, and then use a structural prior distribution with “soft” edge classifications to recover the lost SNR.

In this report, an alternative image reconstruction method based on a combination of MRSI process modeling and Bayesian image analysis is presented for the enhancement of metabolite maps (K-Bayes). First developed for (non-dynamic) perfusion MRI [22, 23], the K-Bayes approach is here expanded considerably to accommodate the additional spectral fitting problem inherent to MRSI. K-Bayes combines information from high-resolution tissue segmented structural MRIs with low-resolution *k*-space-time MRSI data. The acquired MRSI data are supplemented by additional sources of prior information that characterize known relationships between local tissue structural patterns and the spatial distribution of metabolites in the brain. Therefore, K-Bayes effectively performs a full spectral fitting of the data along with incorporation of spatial information. To the authors’ knowledge, there are no other existing reconstruction methods that use a true Bayesian *k*-space-time signal acquisition model to directly reconstruct metabolite maps from low-resolution MRSI data. K-Bayes thereby provides a unified solution for reducing artifacts and increasing image resolution and accuracy.

The main modeling objectives of K-Bayes for MRSI data are: 1) to incorporate a full signal (and noise) model to relate the raw *k*-space-time data to the metabolite maps to be reconstructed and 2) to incorporate anatomical information from high-resolution structural MRIs. This is achieved within the framework of Bayesian image analysis modeling [24-26]. The definition of a full *k*-space-time signal acquisition model alongside a noise model forms the likelihood part of K-Bayes. The prior distribution for K-Bayes relates the spatial distribution of metabolites to a high-resolution and coregistered structural MRI that is segmented into gray matter (GM), white matter (WM), and cerebro-spinal fluid (CSF). The prior model incorporates knowledge of the expected ‘smoother’ behavior of metabolites within homogeneous tissue regions as opposed to less smooth behavior across the boundaries of different tissues; differences between gray and white matter metabolite concentrations have been determined to reach as much as 30% for creatine (Cr) in the temporal lobe and choline (Cho) in the parietal lobe [27]. That is, neighboring locations within homogeneous tissue type are expected to have more similar values than either 1) neighboring voxels of different tissue type, or 2) voxels of homogeneous tissue type that are further apart. These expectations are incorporated into a prior distribution via the use of Markov random field (MRF) models [25, 26, 28] that describe (probabilistically) the spatial distribution of metabolites: characterizing the different spatial correlation structure of metabolites within homogeneous tissue regions compared with that across tissue boundaries. Also incorporated into the prior distribution is the knowledge that metabolite levels should be zero outside the head and can be assumed negligible within CSF. The information in the data (incorporated via the likelihood) and the structural information used in the prior distribution are combined via Bayes’ Theorem. The combination leads to a posterior distribution of metabolite maps that is defined at the higher resolution of the structural MRI. K-Bayes reconstructed metabolite maps are obtained by maximizing the posterior distribution. Maximization is achieved by applying an iterative (generalized) Expectation-Maximization (EM) algorithm [29] to an expanded set of model variables.

The 4D multi-slice MRSI signal model for K-Bayes both corrects and expands upon the model given in Young et al. [4] and Soher et al. [5]. Each signal time series, *s*(*k _{x}, k_{y}, w,t*), for a fixed

$$\begin{array}{cc}\hfill & s({k}_{x},{k}_{y},w,t)=\sum _{m=0}^{M-1}\sum _{n=0}^{N\left(m\right)-1}\sum _{r=wc}^{(w+1)c-1}\sum _{p=0}^{P-1}\sum _{q=0}^{Q-1}{\int}_{z\left(r\right)-\frac{1}{2}}^{z\left(r\right)+\frac{1}{2}}{\int}_{y\left(q\right)-\frac{1}{2}}^{y\left(q\right)+\frac{1}{2}}{\int}_{x\left(q\right)-\frac{1}{2}}^{x\left(p\right)+\frac{1}{2}}{A}_{m}(x,y,z){L}_{mn}\times \hfill \\ \hfill & \mathrm{exp}\left\{-i\left[({\omega}_{m}+{\omega}_{mn}){B}_{0}(x,y,z)t+{\varphi}_{mn}+{\varphi}_{A}(x,y,z)+2\pi \left({k}_{x}x\u2215P+{k}_{y}y\u2215Q\right)\right]-\left[t\u2215{T}_{a}(x,y,z)\right]+{\left[t\u2215{T}_{b}(x,y,z)\right]}^{2}\right\}dxdydz.\hfill \end{array}$$

[1]

The multiple integral of Eq. [1] is evaluated over each voxel at the higher resolution of the structural MRI. The signal *s*(*k _{x}, k_{y}, w,t*) at each

Therefore, the quintuple summation in Equation [1] is over all metabolites, resonance lines and voxels in the *w*-th *k*-space slice. The *A _{m}*(

Note that for 3D *k*-space acquisition sequences (as opposed to multi-slice acquisition), Eq. [1] requires slight modifications so that *s*(*k _{x}, k_{y}, w,t*) becomes

The *T _{a}*,

An important issue that arises in the implementation of this model is how to choose the metabolites and associated resonance lines to include in the model. This precise choice would be experiment and TE-specific. Obviously all metabolites of interest in the study should be included, as well as any other metabolites that are thought to contribute a significant signal (thereby reducing residual noise levels). The papers by Young et al. [3, 4] and Soher et al. [5] address the choice of metabolites and spectral lines for inclusion in MRSI signal models in detail.

1. After performing the integration the model becomes:

$$\begin{array}{cc}\hfill & s({k}_{x},{k}_{y},w,t)=\mathrm{sinc}(\pi {k}_{x}\u2215P)\mathrm{sinc}(\pi {k}_{y}\u2215Q)\sum _{m=0}^{M-1}\sum _{n=0}^{N\left(m\right)-1}\sum _{r=wc}^{(w+1)c-1}\sum _{p=0}^{P-1}\sum _{q=0}^{Q-1}{A}_{m}(p,q,r){L}_{mn}\times \hfill \\ \hfill & \mathrm{exp}\left\{-i\left[({\omega}_{m}+{\omega}_{mn}){B}_{0}(p,q,r)t+{\varphi}_{mn}+{\varphi}_{A}(p,q,r)+2\pi \left({k}_{x}p\u2215P+{k}_{y}q\u2215Q\right)\right]-\left[t\u2215{T}_{a}(p,q,r)\right]+{\left[t\u2215{T}_{b}(p,q,r)\right]}^{2}\right\}.\hfill \end{array}$$

[2]

where sinc(*x*) = sin(*x*) / *x* . An implicit assumption for this model is that *A _{m}*,

The process *s*(*k _{x}*,

$$d({k}_{x},{k}_{y},w,t)=s({k}_{x},{k}_{y},w,t)+\epsilon ({k}_{x},{k}_{y},w,t).$$

[3]

The error/noise component is modeled as isotropic complex Gaussian noise: *ε*(*k _{x}*,

This leads to the following likelihood for the full dataset:

$$\pi (d\mid s)\propto \mathrm{exp}\left\{-\frac{1}{{\sigma}^{2}}\sum _{t}\sum _{w}\sum _{({k}_{x},{k}_{y})}{\mid d({k}_{x},{k}_{y},w,t)-s({k}_{x},{k}_{y},w,t)\mid}^{2}\right\}$$

[4]

where *π*(·) is used to denote a generic probability distribution. Note that indices are omitted to indicate the full array, i.e. writing *d* instead of *d*(*k _{x}*,

A segmented and co-registered structural MRI (that could be acquired within the same experimental session as the MRSI) provides a readily available and detailed source of prior spatial information. The acquired structural MRI is segmented into GM, WM, and CSF, using one of the many available segmentation algorithms [30, 31]. The goal of the prior distribution is to characterize known behavior of metabolites within homogeneous regions of tissue. The segmentation information is incorporated into the model via a prior distribution that satisfies two criteria: metabolite levels should be encouraged to vary smoothly across space within each of the tissue types of GM and WM, and sharper changes in metabolite levels should be considered more feasible at the identified GM-WM boundaries and at edges of the brain or ventricles. Regions outside the brain or areas of CSF should not have a detectable metabolite signal. These criteria provide information to the K-Bayes model regarding the spatial distribution of metabolites within homogeneous and heterogeneous tissue regions and can be mathematically represented within a prior distribution based on MRF models [26, 28, 32]. The prior model used here is defined at the spatial resolution of the structural MRI as follows:

$$\pi \left({A}_{m}\right)\propto \mathrm{exp}\left\{-\frac{1}{2}\left[\sum _{\langle (p,q,r),({p}^{\prime},{q}^{\prime},{r}^{\prime})\rangle}\left(\left[\sum _{\Delta \in \{B,G,W\}}\frac{1}{{\tau}_{m\Delta}^{2}}{I}_{\Delta}[(p,q,r),({p}^{\prime},{q}^{\prime},{r}^{\prime})]\right]{[{A}_{m}(p,q,r)-{A}_{m}({p}^{\prime},{q}^{\prime},{r}^{\prime})]}^{2}\right)\right]\right\}$$

[5]

with the added constraint that *A _{m}* is assumed zero everywhere both outside the brain and within regions consisting of CSF. (This condition can be relaxed if required.) The sum in Eq. [5] over < (

In order to obtain an estimate (reconstruction) of the metabolite maps, the posterior distribution is maximized with respect to the sets of metabolite amplitudes at all voxels. This is equivalent to maximizing the product of the likelihood and the prior distribution and is only one of many possible estimates that can be obtained from the posterior distribution. However, it is a common choice for Bayesian image analysis because a) it is easy to understand, and b) there is known (and relatively computationally efficient) methodology that can be used to obtain it [26]. A generalized form of the Expectation-Maximization (EM) procedure [29] was used for optimization. The algorithm is generated in a similar fashion to the algorithms in Green [35] and Miller et al. [9], where a set of unobservable latent variables is created to simplify the computational algorithm. These variables are treated as missing data within the EM algorithm. The latent variables for the EM reconstruction approach described here are defined as:

$$\begin{array}{cc}\hfill & z({k}_{x},{k}_{y},p,q,r,w,t,m)=\mathrm{sinc}(\pi {k}_{x}\u2215P)\mathrm{sinc}(\pi {k}_{y}\u2215Q)\sum _{n=0}^{N\left(m\right)-1}{A}_{m}(p,q,r){L}_{mn}\hfill \\ \hfill & \mathrm{exp}\left\{-i\left[({\omega}_{m}+{\omega}_{mn}){B}_{0}(p,q,r)t+{\varphi}_{mn}+{\varphi}_{A}(p,q,r)+2\pi \left({k}_{x}p\u2215P+{k}_{y}q\u2215Q\right)\right]-[t\u2215{T}_{a}(p,q,r)]+{[t\u2215{T}_{b}(p,q,r)]}^{2}\right\}\hfill \\ \hfill & +{\epsilon}_{z}({k}_{x},{k}_{y},p,q,r,w,t,m).\hfill \end{array}$$

[6]

The term *z*(*k _{x}*,

Note that the individual *z* -variables are not actually defined in a rigorous sense because of the uncertainty principle, i.e. it is impossible to simultaneously specify a point precisely in both *k* -space and image-space [36]. However, because we only consider different linear combinations of the *z* -variables (which are well defined) within the M-step of the EM optimization procedure, our methodology remains valid at all times.

For a certain set of starting values an iterative procedure is taken where firstly the expectation is calculated for each *z* -variable given the current values of the metabolite amplitudes *A _{m}* and all other

$$\begin{array}{c}\hfill z{({k}_{x},{k}_{y},p,q,r,w,t,m)}^{j}=\frac{1}{\nu}\left[\begin{array}{c}d({k}_{x},{k}_{y},w,t)-\mathrm{sinc}(\pi {k}_{x}\u2215P)\mathrm{sinc}(\pi {k}_{y}\u2215Q)\hfill \\ \sum _{m=0}^{M-1}\sum _{n=0}^{N\left(m\right)-1}{L}_{mn}\sum _{{r}^{\prime}=wc}^{(w+1)c-1}\sum _{({p}^{\prime},{q}^{\prime},{r}^{\prime})\ne (p,q,r)}\mathrm{exp}\left\{\begin{array}{c}-i\left[\begin{array}{c}({\omega}_{m}+{\omega}_{mn}){B}_{0}({p}^{\prime},{q}^{\prime},{r}^{\prime})t+{\varphi}_{mn}\hfill \\ +{\varphi}_{A}({p}^{\prime},{q}^{\prime},{r}^{\prime})+2\pi ({k}_{x}{p}^{\prime}\u2215P+{k}_{y}{q}^{\prime}\u2215Q)\hfill \end{array}\right]\hfill \\ -[t\u2215{T}_{a}({p}^{\prime},{q}^{\prime},{r}^{\prime})]+{[t\u2215{T}_{b}({p}^{\prime},{q}^{\prime},{r}^{\prime})]}^{2}\hfill \end{array}\right\}\hfill \end{array}\right]\hfill \\ \hfill \phantom{\rule{1em}{0ex}}\phantom{\rule{1em}{0ex}}\phantom{\rule{1em}{0ex}}\phantom{\rule{1em}{0ex}}\phantom{\rule{1em}{0ex}}\phantom{\rule{1em}{0ex}}\phantom{\rule{1em}{0ex}}\phantom{\rule{1em}{0ex}}\phantom{\rule{1em}{0ex}}\phantom{\rule{1em}{0ex}}\phantom{\rule{1em}{0ex}}\phantom{\rule{1em}{0ex}}\phantom{\rule{1em}{0ex}}\phantom{\rule{1em}{0ex}}\phantom{\rule{1em}{0ex}}\phantom{\rule{1em}{0ex}}\phantom{\rule{1em}{0ex}}\phantom{\rule{1em}{0ex}}\phantom{\rule{1em}{0ex}}\phantom{\rule{1em}{0ex}}\phantom{\rule{1em}{0ex}}\phantom{\rule{1em}{0ex}}\phantom{\rule{1em}{0ex}}\phantom{\rule{1em}{0ex}}\phantom{\rule{1em}{0ex}}\phantom{\rule{1em}{0ex}}\phantom{\rule{1em}{0ex}}\phantom{\rule{1em}{0ex}}\phantom{\rule{1em}{0ex}}\phantom{\rule{1em}{0ex}}\phantom{\rule{1em}{0ex}}\phantom{\rule{1em}{0ex}}\phantom{\rule{1em}{0ex}}\phantom{\rule{1em}{0ex}}\phantom{\rule{1em}{0ex}}\phantom{\rule{1em}{0ex}}\phantom{\rule{1em}{0ex}}\phantom{\rule{1em}{0ex}}\phantom{\rule{1em}{0ex}}\phantom{\rule{1em}{0ex}}\phantom{\rule{1em}{0ex}}\phantom{\rule{1em}{0ex}}\phantom{\rule{1em}{0ex}}\phantom{\rule{1em}{0ex}}\phantom{\rule{1em}{0ex}}\phantom{\rule{1em}{0ex}}\phantom{\rule{1em}{0ex}}\phantom{\rule{1em}{0ex}}\phantom{\rule{1em}{0ex}}\phantom{\rule{1em}{0ex}}\phantom{\rule{1em}{0ex}}+\frac{\nu -1}{\nu}z{({k}_{x},{k}_{y},p,q,r,w,t,m)}^{j-1}\hfill \end{array}$$

[7]

where *ν* = *cPQ* the number of high-resolution voxels in each structural MRI slice.

The corresponding M-step is:

$$\begin{array}{cc}\hfill & {A}_{m}{(p,q,r)}^{j}=\hfill \\ \hfill & \phantom{\rule{1em}{0ex}}\frac{1}{{\sigma}^{2}}\sum _{{k}_{x},{k}_{y}}\mathrm{sinc}(\pi {k}_{x}\u2215P)\mathrm{sinc}(\pi {k}_{y}\u2215Q)\sum _{(t\u2215\delta )=0}^{T-1}\mathrm{exp}\left\{-\left[t\u2215{T}_{a}(p,q,r)\right]+{\left[t\u2215{T}_{b}(p,q,r)\right]}^{2}\right\}\hfill \\ \hfill & \phantom{\rule{1em}{0ex}}\sum _{n=0}^{N\left(m\right)-1}{L}_{mn}\left(\begin{array}{c}\mathrm{Re}\left[z{({k}_{x},{k}_{y},p,q,r,w,t,m)}^{j}\right]\mathrm{cos}\left(\begin{array}{c}({\omega}_{m}+{\omega}_{mn}){B}_{0}(p,q,r)t+{\varphi}_{mn}\hfill \\ +{\varphi}_{A}(p,q,r)+2\pi ({k}_{x}p\u2215P+{k}_{y}q\u2215Q)\hfill \end{array}\right)-\hfill \\ \mathrm{Im}\left[z{({k}_{x},{k}_{y},p,q,r,w,t,m)}^{j}\right]\mathrm{sin}\left(\begin{array}{c}({\omega}_{m}+{\omega}_{mn}){B}_{0}(p,q,r)t+{\varphi}_{mn}\hfill \\ +{\varphi}_{A}(p,q,r)+2\pi ({k}_{x}p\u2215P+{k}_{y}q\u2215Q)\hfill \end{array}\right)\hfill \end{array}\right)\hfill \\ \hfill & \phantom{\rule{1em}{0ex}}\frac{+\sum _{\langle ({p}^{\prime},{q}^{\prime},{r}^{\prime})\in \delta (p,q,r)\rangle}\left(\sum _{\Delta \in \{B,G,W\}}\frac{1}{{\tau}_{m\Delta}^{2}}{I}_{\Delta}[(p,q,r,),({p}^{\prime},{q}^{\prime},{r}^{\prime})]\right){A}_{m}{({p}^{\prime},{q}^{\prime},{r}^{\prime})}^{{j}^{\ast}}}{\frac{1}{{\sigma}^{2}}\sum _{{k}_{x},{k}_{y}}{\left(\mathrm{sinc}(\pi {k}_{x}\u2215P)\mathrm{sinc}(\pi {k}_{y}\u2215Q)\right)}^{2}+\sum _{\langle ({p}^{\prime},{q}^{\prime},{r}^{\prime})\in \delta (p,q,r)\rangle}\sum _{\Delta \in \{B,G,W\}}\frac{1}{{\tau}_{m\Delta}^{2}}{I}_{\Delta}[(p,q,r),({p}^{\prime},{q}^{\prime},{r}^{\prime})]}\hfill \end{array}$$

[8]

where *j** can be taken as *j* – 1, or for slightly improved performance, the latest available update (either *j* – 1 or *j*).

To implement the K-Bayes algorithm the following steps are required.

- Perform conventional DFT-based reconstruction of the MRSI data using e.g. freely available SITools-FITT
- Spatially zero-fill all fitted parameter maps
*A*,_{m}*T*,_{a}*T*,_{b}*ϕ*,*B*_{0}, effectively interpolating the maps to structural MRI resolution. - Coregister and re-slice high-resolution structural MRI data to match interpolated DFT reconstructed maps.
- Segment re-sliced structural MRI into GM, WM, CSF and non-brain for use as prior information in K-Bayes.
- Use original
*k*-space-time data and segmented re-sliced structural MRI as input to K-Bayes. The interpolated DFT fitted metabolite maps can be used as start values for the K-Bayes reconstructed maps. - Iterate between E- and M-steps of the K-Bayes algorithm until the reconstructed metabolite maps have all converged to within some specified acceptable tolerance.

The simulation study compared K-Bayes and conventional DFT reconstruction under controlled conditions. The study utilized simplified versions of the full K-Bayes model described in Eqs. [1] and [2]. The first general simplification was that all phase and frequency offset terms were set to zero and that the *T _{a}* and

In 2D, with single metabolite frequencies, Eq. [2] reduces to:

$$\begin{array}{cc}\hfill & s({k}_{x},{k}_{y},t)=\mathrm{sinc}(\pi {k}_{x}\u2215P)\mathrm{sinc}(\pi {k}_{y}\u2215Q)\sum _{m=0}^{M-1}\sum _{p=0}^{P-1}\sum _{q=0}^{Q-1}{A}_{m}(p,q)\times \hfill \\ \hfill & \mathrm{exp}\left\{-i\left[{\omega}_{m}{B}_{0}(p,q)t+{\varphi}_{A}(p,q)+2\pi \left({k}_{x}p\u2215P+{k}_{y}q\u2215Q\right)\right]-\left[t\u2215{T}_{a}(p,q)\right]+{\left[t\u2215{T}_{b}(p,q)\right]}^{2}\right\}\hfill \end{array}$$

[9]

and the likelihood and prior distribution are similarly reduced to only two spatial dimensions. The extensions to 3D space and short TE data are straightforward (in terms of coding), but have increased computational burdens.

Figure 1 outlines the data simulation and reconstruction procedure. The simulated MRSI dataset was generated based on the Montreal Neurological Institute (MNI) structural MRI of the human brain [37] that had been segmented into GM, WM, and CSF. The *k*-space-time data was created by first generating a high-resolution metabolite intensity map with concentrations of 1.0 in GM and 0.5 in WM from the segmented MNI brain (128×128). The ensuing map was used to generate individual high-resolution metabolite maps for each of N-acetylaspartate (NAA), Cr and Cho with relative amplitudes of 1.0, 0.25, and 0.5, respectively. The difference in overall intensity for the three metabolite maps seen in Figure 2 is indicative of the relative metabolite concentrations. ‘Hotspot’ high-intensity regions were added in the WM regions for NAA and Cho such that the WM intensity was doubled (indicated by arrows in Figure 2). The hotspots were included for the purpose of testing whether K-Bayes was robust to spatial changes not matching prior expectations; i.e., to what extent would the MRF prior distribution smooth away metabolite variation that did not correspond to structural boundaries. Finally, nearest-neighbor spatial smoothing was applied across the metabolite maps to simulate a small bleeding effect across tissue boundaries: each voxel's value was replaced by the mean of the combination of itself and its four nearest neighbors.

Simulation study procedure: Simulated high-resolution metabolite maps are generated based on a structural MRI segmented into GM, WM, and CSF. At each high-resolution voxel the time-series signal was generated as the sum of exponentially decaying metabolite **...**

Reconstructions of MRSI simulated data for three metabolites each of which has a single resonance frequency. Each row corresponds to a single metabolite and the columns correspond to the truth, zDFT and K-Bayes reconstructions, respectively. The true **...**

High-resolution image-space-time data (length 128 timepoints and sampling interval of 1.0 milliseconds) were generated from the set of metabolite intensity maps. At each voxel the signal was generated as the sum of three exponentially decaying metabolite frequencies (NAA, Cr, Cho). The metabolite frequencies were simulated at 2.0 parts per million (ppm) for NAA, 3.0 ppm for Cr and 3.2 ppm for Cho. The signal's magnitude at each frequency was defined by the intensity of the voxel in the associated metabolite map. The signal was simulated separately for each point in high-resolution-image-space in order to approximate a continuous spatial process and the data were then DFT'd spatially into *k*-space-time. The central 32×32 region of *k*-space was then extracted to yield low-spatial-resolution *k*-space-time data. Finally, independent complex Gaussian noise was added at each *k*-space-timepoint with standard deviation of 0.1. This low-resolution data along with the original segmentation (i.e. tissue type) data became the input for K-Bayes.

We tested different combinations of values for σ^{2}, ${\tau}_{mB}^{2}$, ${\tau}_{mG}^{2}$, and ${\tau}_{mW}^{2}$ within the model, and found that the resulting reconstructions were reasonably robust to specification of the smoothing parameters in the MRF model. The actual values used for the simulation results shown in Figure 2 were σ^{2} = 0.1, ${\tau}_{mB}^{2}=2.0$, ${\tau}_{mG}^{2}=0.001$, and ${\tau}_{mW}^{2}=0.004$, with the values kept the same for all metabolites, *m*. The increased strength of *a priori* smoothness in GM over WM (i.e. ${\tau}_{mG}^{2}=0.001$ compared with ${\tau}_{mW}^{2}=0.004$) is required to balance the reduced connectedness amongst GM neighbors as compared with WM (i.e., because GM consists of the relatively thin surface of the cortex). The absolute value of σ^{2} does not directly influence the posterior estimation of the metabolite maps; only the relative magnitude of σ^{2} compared with the prior smoothing parameters is important. This lack of identifiability for the combined set of σ^{2}, ${\tau}_{mB}^{2}$, ${\tau}_{mG}^{2}$, and ${\tau}_{mW}^{2}$ can be seen by examination of the M-step of the EM algorithm in Eq. [8].

The quality of the K-Bayes reconstruction after suitable convergence of the algorithm was compared to that for the (spatially) zero-filled discrete Fourier transform (zDFT). The zDFT reconstruction was obtained by first performing a spatial DFT for each timepoint's image, leading to a time-series of images in image space. Subsequently, a least-squares time-domain fitting procedure was applied at each voxel to provide the metabolite amplitude estimates. The model fitted by least squares was the same model used to simulate the data, i.e., the sum of a set of exponentially decaying sinusoids. The metabolite amplitudes were the unknown parameters to be estimated and all other parameters were fixed to those used in the simulation. Zero-filling in the spatial dimension increased the spatial resolution of the DFT-based reconstruction to that of the gold standard high-resolution truth (i.e. the true intensity maps used to simulate the data) and the K-Bayes reconstruction. The use of consistent resolution permitted direct comparison of the different reconstruction methods across voxels with respect to the gold standard. We further examined bicubic spline interpolated DFT maps (sDFT) to determine whether the choice of interpolation scheme might be important. Comparisons with the gold standard high-resolution truth were made in terms of i) visual examinations of the reconstructed maps; ii) the average bias in each of GM, WM, and the hotspots; iii) the root mean square error (RMSE) across all tissue and in only the hotspots. Bias in this context is defined as $\frac{1}{N}{\sum}_{i=1}^{N}({y}_{i}-{x}_{i})$ and RMSE as $\sqrt{\frac{1}{N}{\sum}_{i=1}^{N}{({y}_{i}-{x}_{i})}^{2}}$, where *y* is the gold standard, *x* is the reconstructed image, and *i* indexes over the set of voxels for which the metric is to be evaluated.

Figure 2 shows the reconstructions for the *k*-space-time simulations. The maps consist of: Left, the truth (gold standard), i.e. the three original 128×128 maps for each of three metabolites; Center, the corresponding zDFT reconstructions; and Right, the K-Bayes reconstructions. There is considerable visual improvement in each of the K-Bayes reconstructions compared with zDFT. The sDFT reconstruction is very similar to the zDFT reconstruction and is therefore not shown. The K-Bayes reconstructions have increased spatial definition compared with zDFT. Furthermore, the zDFT reconstructed maps contain Gibbs ringing not visible in the K-Bayes ones. K-Bayes appears to preserve the hotspots similarly to the zDFT reconstruction, apparently signaling that in this example K-Bayes is as robust as zDFT when detecting metabolite level changes that are not defined by tissue structure boundaries.

The qualitative visual improvements are verified numerically in Table 1, where K-Bayes consistently does better than zDFT and sDFT; the two DFT-based schemes led to very similar results, but sDFT did slightly worse across all metrics than zDFT. Across the three metabolites, GM and WM bias for K-Bayes was always less than 6% of that for the DFT-based reconstructions, and RMSE for K-Bayes was always less than 50% of that for either DFT approach. Perhaps surprisingly, even in the hotspot we see that K-Bayes has less than 35% of the bias of the DFT-based reconstructions and less than 50% of the RMSE. Although K-Bayes imposes smoothness on the high-resolution reconstruction it does so with guidance from the structural prior. In contrast, the interpolation imposed on the DFT-based reconstructions implies an arbitrary level smoothness that ignores both structural information and the noise in the original signal (which is also incorporated into the interpolation).

Statistical comparison of the reconstruction methods for the simulated MRSI data reconstructions in terms of bias and RMSE. The first column gives the metric and the second column subsets this by metabolite. The final three columns respectively give the **...**

Figure 3 displays maps of the absolute residuals from the truth for each of the three metabolites NAA, Cr and Cho. The first two columns display the absolute zDFT and absolute K-Bayes residuals, respectively. These absolute residual maps are all constrained to be on the same scale. It is clear that for all metabolites K-Bayes has (on average) reduced absolute residuals compared with zDFT. This reduction in average intensity simply reflects the overall reduction in RMSE expressed in Table 1. The third column shows the zDFT minus K-Bayes difference (Diff.) in absolute residuals. The greatest reductions in residual error are seen as the brightest voxels, and they primarily occur in GM cortical regions and the WM hotspot. zDFT does a little better than K-Bayes for a small proportion (~10%) of voxels in the brain. These are shown as the darkest voxels in the Diff. map, typically close to the GM boundary in WM. The likely cause of the reduced accuracy of K-Bayes in these voxels is the effect of the cross-tissue-type boundary smoothing/variance parameter ${\tau}_{mB}^{2}$; although zDFT also smoothes across the boundary, it is doing so based on underestimated GM intensities and therefore the nearby WM intensities are less elevated. (It is noted earlier that the original spatial smoothing of the simulated data was applied in a highly localized fashion.)

Maps displaying the absolute values of residuals (relative to the truth) for the three metabolites NAA, Cr and Cho. Each row corresponds to a single metabolite and the columns correspond to the absolute zDFT residuals, absolute K-Bayes residuals and the **...**

To determine the extent of K-Bayes’ capability for providing gains over zDFT for highly localized anomalies we repeated the reconstruction procedure and performed comparisons for a range of hotspot sizes (not surprisingly the metrics did not vary much outside of the hotspot). Figure 4 displays plots of Bias and RMSE for zDFT and K-Bayes across the varying hotspot sizes (in terms of the square root of the number of voxels contained in the hotspot) for both the NAA and Cho hotspots. K-Bayes produced highly consistent improvements (reduced magnitude Bias and RMSE) over zDFT across hotpot sizes.

To determine the sensitivity of the K-Bayes reconstruction to the parameters of the anatomic distribution we performed reconstructions of NAA maps for 16 sets of parameters (hotspot radius as in original simulation). These parameter sets covered the ranges $0.1\le {\tau}_{B}^{2}\le 40.0$, $0.001\le {\tau}_{G}^{2}\le 1.0$ and $0.002\le {\tau}_{W}^{2}\le 5.0$: thus covering two orders of magnitude in each parameter. Figure 5 graphically summarizes the statistical metrics for the associated K-Bayes and zDFT fits. The first thing to note is that no matter which set of anatomic parameters was chosen, K-Bayes performed much better than zDFT in all metrics. Furthermore all K-Bayes reconstructions are tightly clustered relative to the zDFT value for all metrics (with just a couple of exceptions for the hotspot metrics). That is, despite large shifts in the parameters, the reconstructed images had very similar properties in all cases. To show that this relatively stable solution with respect to the parameters provides an appropriate balance between the accurate reconstruction of global effects, small focal effects, and noise removal, we must examine the individual box plots more closely. We note that there is consistently very tight clustering of the K-Bayes reconstructions in the non-hotspot metrics. For the hotspot metrics we see two exceptions. One set of parameters that leads to an outlier (lowest *a priori* smoothness level) only increases the RMSE. In these cases the hotspot is flattened by too strong smoothing parameters causing increased bias and RMSE simultaneously (all the error is to one side of the truth). Nonetheless, even these most limited of K-Bayes reconstructions still have greatly improved hotspot bias and RMSE compared with zDFT.

Boxplots of statistical metrics for zDFT reconstruction and 16 K-Bayes NAA reconstructions with differing anatomical parameter sets. zDFT is always the single point to the right (i.e. the worst score).

This sensitivity analysis indicates that K-Bayes is robust to a wide range of anatomic parameter sets and hence that the optimal balance between the anatomic information and the *k*-space-time data is well defined by the two data sources and not substantially influenced by any particular choice of parameter set.

Overall, these simulation results generate a very encouraging assessment of the potential gains afforded by K-Bayes in the reconstruction of MRSI datasets. However, they do provide close to an upper bound on the gains because the datasets are generated based on the K-Bayes model. That is, the signal model used is from the likelihood, and the metabolite maps are generated from the structural information (albeit with added blurring and hotspots). Inaccuracies in the signal model will equally propagate into errors for DFT-based reconstructions as well as K-Bayes. The potential realization of the relative gains of K-Bayes for real data reconstruction is therefore bounded by the extent to which differentials in metabolite levels respect tissue boundaries. Despite this limitation, some gains are made by K-Bayes in the simulation study even when metabolite differentials do not match the tissue boundaries, i.e. in the hotspot.

A spin-echo MRSI dataset from a male healthy volunteer, age 42 years, was reconstructed using K-Bayes and zDFT. The dataset was acquired at 3T with TR/TE = 1400/70ms, a 36×36 *k*-space matrix and a 256×256mm field-of-view (FOV). A corresponding GM/WM/CSF segmented magnetization prepared rapid acquisition gradient echo (MPRAGE) scan with nominal voxel size of 1×1×1.5mm across a 256×256mm FOV served as the structural information for K-Bayes. Each k-space-time-series consisted of 512 points with a sampling interval of 0.649ms, i.e., spectral bandwidth of 1540Hz. Figure 6 displays the structural MPRAGE scan and corresponding segmented map.

Structural view of spin-echo data slice: (a) MPRAGE structural scan (b) GM/WM segmented version of MPRAGE.

The data were initially processed using the SITools spectral analysis software [38]. Using the signal model of Young et al. [4] and Soher et al. [5], metabolite, residual water, and lipid signals were identified in all voxels and removed from all spectra while maintaining metabolite signals through the use of a singular value decomposition (SVD) based fitting algorithm. Residual signals left in most spectra outside the range of the predominant metabolites (1.9 to 3.9 PPM) appeared close to white noise. This ‘metabolite-only’ dataset was 1) DFT'd back into *k*-space-time for use in K-Bayes and 2) fitted using the SITools-FITT module to create NAA, Cr, and Cho metabolite maps. As part of the SITools-FITT, small baseline signal contributions due to macromolecular peaks in the metabolite region were identified and fitted. Metabolite maps were sinc interpolated to 256×256mm to match the voxel size of the MPRAGE acquisition and K-Bayes reconstruction. The spatially zero-filled frequency, phase, Lorentzian, and Gaussian decay estimates plus fitted macromolecular baseline signals were input into K-Bayes as fixed parameters. The zDFT estimates of the three metabolite signals (NAA, Cr, Cho) were used as starting maps for K-Bayes. Parameters for the structural prior distribution were obtained by trial and error. However, we found through a robustness study (not shown) that reconstructions were quite robust to within an order of magnitude change in the prior parameters. This requires further verification on multiple datasets. If K-Bayes is to be applied in multi-subject studies, we advocate that a structured approach to parameter setting be taken, such as that described in Section III.

The reconstructed maps in Figure 7 show that K-Bayes produces additional detail (relative to zDFT reconstruction) unrelated to the underlying anatomy as well as showing variation clearly related to anatomy. These encouraging results suggest that K-Bayes is capable of finding a reasonable balance between the accurate reconstruction of global effects, small focal effects, and noise.

Spin-echo metabolite map reconstructions: zDFT (left column) and K-Bayes (right column) reconstructions for NAA, Cr, and Cho.

To determine the sensitivity of K-Bayes to the parameters of the anatomic distribution, we performed reconstructions for 8 sets of anatomic parameters. For the three parameters, these 8 sets covered the range $0.1\le {\tau}_{B}^{2}\le 50.0$, $0.0005\le {\tau}_{G}^{2}\le 0.01$, $0.0005\le {\tau}_{W}^{2}\le 0.01$. These ranges are all greater than an order of magnitude. Although the ranges appear different from those in the simulation, they are in fact quite similar in the relative magnitude of the parameters to the (square of the) range of the data.

The eight K-Bayes reconstructed NAA maps are shown in Figure 8, contrasted with the single zDFT reconstruction. (We show only figures since we lack a gold standard with which to assess statistical metrics.) Clearly, variation in K-Bayes reconstructions is small compared to the difference between them and the zDFT reconstruction, indicating that K-Bayes is robust to a wide range of anatomical distribution parameter sets. Thus the optimal way to balance anatomic and *k*-space-time data is well defined in the data sources themselves. More precisely, the location of the posterior maximum (the optimal reconstruction) does not vary greatly for a wide range of prior smoothing parameter values; despite large shifts in parameters, the image at the posterior maximum is quite similar in all cases.

A proton resolved spectroscopy (PRESS) chemical shift imaging (CSI) (“PRESS box”) dataset was acquired from a single male healthy volunteer, age 40. The PRESS box covered a region of 80×80×15mm, acquired as a 32×32 *k*-space matrix over a 240×240mm FOV. A corresponding GM/WM/CSF segmented structural MRI (MPRAGE) with voxel dimensions 1×1×1.5mm and 256×256mm FOV was used as the structural information for K-Bayes. Each k-space-time-series consisted of 1024 points with a sampling interval of 0.649ms, i.e., spectral bandwidth of 1540Hz. Figure 9 displays the structural scan with a white box outlining the PRESS box region.

Spatial location of PRESS-CSI data: Top left: MPRAGE structural scan with PRESS box outlined in white. Top right: PRESS box zoomed. Bottom left: segmented MPRAGE. Bottom right: zoomed version of segmented PRESS box.

SITools-4DFT and SITools-FITT were employed for conventional spectral DFT and fitting of the data [38] using the signal model of Young et al. [4] and Soher et al. [5]. The complex baseline signal at each voxel was fit using a spline model with knots placed 30 frequency points apart in the spectrum. There was considerable residual water signal in the dataset, 10-100 times the area of the NAA peak, despite applying water suppression during signal acquisition. The residual water signal was typically observed to be Lorentzian shaped, but in many voxels was significantly asymmetric or had an actual shoulder. We determined experimentally that the residual water signal could be fit with three extra frequencies located at 4.6, 4.7, and 4.8ppm sufficiently well to provide, along with the spline baseline, reasonable removal of the residual water signals. After water and baseline signal removal, the residuals left in most spectra (outside the range of the predominant metabolites) appeared close to white noise. The resulting integrated spectral maps for each metabolite and frequency components of water signal were subsequently zero-filled to 240×240mm in the *k*-space dimensions to match the voxel size of the MPRAGE acquisition.

ZDFT reconstructed maps of the B_{0}, phase, Lorentzian and Gaussian decay map estimates (plus the fitted macromolecular baseline signals) from SITools-FITT were used as fixed quantities in K-Bayes. The corresponding zDFT map estimates from SITools-FITT of the three metabolite signals (NAA, Cr, Cho), as well as the three frequency components of water signal, were used as starting maps for K-Bayes. Parameters for the structural prior distribution were obtained by trial and error. There appeared to be more sensitivity in the reconstruction to the specific choice of parameter values for the PRESS box reconstruction. This is possibly due to the limitations of using K-Bayes with this particular PRESS box dataset, as discussed below).

Figure 10 displays zDFT and K-Bayes reconstructed maps of NAA, Cr, and Cho. Qualitatively, K-Bayes produces improved definition and reduced noise in the reconstructed metabolite maps. Although these PRESS results demonstrate feasibility, there are limitations. The PRESS box dataset is far from the ideal form for utilizing K-Bayes; there was a sharp drop-off at the top of the PRESS box that led to bias in the neighboring homogeneous tissue regions that do contain signal affected by the MRF prior (i.e., biased towards zero in regions of high signal). Note that the conservative approach was taken by assuming that the prior PRESS box segmented tissue map corresponded to the a priori defined PRESS box. This approach is in contrast to trying to determine the PRESS box after looking at the data. Further evidence of the inaccuracy in spatial localization of the PRESS box can be seen in the edge-effects of the upper right corner of the K-Bayes reconstructed Cho map: any signal outside the demarked PRESS box is likely being “squeezed into the edge of the box”.

PRESS-CSI metabolite map reconstructions: zDFT (left column) and K-Bayes (right column) reconstructions for NAA, Cr, and Cho.

A further limitation of both real MRSI data studies was that we only considered a single-slice 2D segmented structural MRI as prior information (the middle slice). Since both datasets consisted of 15mm thick slices, there would be variation in tissue composition along the slice in the *z* -direction. Taking *z* -direction variability into account via 3D spatial modeling would likely further improve K-Bayes reconstruction quality. This is discussed further in the next section.

Although MRSI provides useful regional information on metabolite levels in the brain, the results of conventional reconstruction methods are limited by the low spatial resolution inherent in the data acquisition imposed by typical signal to noise ratio considerations. K-Bayes can help quantify regional metabolite levels with much greater precision than that attainable through DFT/spectral-integration-based reconstruction. The increased precision is due to the enhanced spatial specification obtained by incorporating segmented structural MRI prior information, and the use of a full signal model for the *k*-space-time data. The increased precision of estimated metabolite levels could considerably benefit clinical studies that quantify regional differences in metabolite concentrations between different subject populations (e.g. patients with Alzheimer's disease versus healthy controls). The extra accuracy obtained by K-Bayes could be crucial in detecting real differences between disease states. Furthermore, the use of K-Bayes could expand to non-brain applications of MRSI. All that is required to implement K-Bayes is a high-resolution structural segmentation of the anatomy along with a physiologically reasonable model of the distribution of metabolites concentrations within and between structurally homogeneous regions.

In the current implementation of K-Bayes the maps *T _{a}*(

We know that the terms *T _{a}*(

Currently, *T _{a}*(

In the long run we aim to incorporate estimation of *T _{a}*(

In our real data examples (spin-echo and PRESS) the slice thickness of the MRSI was much thicker than that of the corresponding structural MRI slice used as prior information. The 15mm slice will contain partial voluming effects since the tissue structure will not be perfectly homogeneous across the MRSI slice. This is a notable limitation of the real data examples. The complete solution is to extend application of K-Bayes to full 3D-space datasets. Considerable computational enhancement will be required in order to make this practical.

In our spin-echo and PRESS data examples, we selected the central segmented structural MRI slice to define the anatomical prior information in the simplified 2D space case rather than the many other reasonable choices: e.g. defining a prior over blurred boundaries or taking the voxel-wise tissue class mode across all structural MRI slices. We chose the central slice partly for simplicity, partly because it represents the tissue class at the point of highest excitation in the slice, and finally because the central anatomic image represents a true slice through the brain, which is helpful for visualization purposes.

Two methodological issues should be addressed before K-Bayes can achieve more general applicability. The first is that the method requires good co-registration of segmented structural MRI with MRSI data. Similarly, poor quality segmentation can reduce the benefits of K-Bayes. Poor segmentation again corresponds to poor prior information, though unless the segmentation has high mis-classification rates, the problems associated with low-resolution MRSI data will be worse than any error introduced by inaccurate segmentation. The second (and primary) issue encountered with using K-Bayes is computational time. The 2D reconstruction examples were performed on a Dell Precision 370 desktop computer running RedHat Enterprise Linux 4.0 on a single thread of an Intel Xeon 2.33GHz processor. These were programmed in C and compiled with a gcc compiler. Reconstructions reached reasonable convergence in a day. Although it is impractical to achieve full convergence to the level of machine tolerance, the statistical metrics of the reconstructions did not change significantly after this time. Parallelization similar to that used in [10] will help to significantly reduce computational time. Because the algorithm is highly parallelizable, a large reduction in computational time can be achieved by distributing the processing and spreading the memory load amongst multiple processors. The benefit of parallelization will potentially be greater than linear speed up with the number of processors, because the significant memory burden will be distributed as well. Computational improvement of this level would potentially allow for K-Bayes to be used as a routine approach for obtaining high-resolution and accurate metabolite maps. Given the decreasing prices for multiprocessor supercomputers and random access memory (RAM), there is great potential for K-Bayes to become widely applicable in clinical settings.

This work was supported by a Society for Imaging Informatics in Medicine (SIIM) research grant; NIH grants R01 NS41946, P41 RR023953, R01 EB008387. We acknowledge the San Francisco Veterans Affairs Medical Center where part of this work was performed.

The E-step update applied to each missing data point *z*(*k _{x}*,

Bayes’ Theorem gives

$$\pi \left(z({k}_{x},{k}_{y},p,q,r,w,t,m)\mid d({k}_{x},{k}_{y},w,t),A\right)=\frac{\pi \left(d({k}_{x},{k}_{y},w,t)\mid z({k}_{x},{k}_{y},p,q,r,w,t,m),A\right)\pi \left(z({k}_{x},{k}_{y},p,q,r,w,t,m)\mid A\right)}{\pi \left(d({k}_{x},{k}_{y},w,t)\mid A\right)}.$$

Now,

$$\begin{array}{cc}\hfill & d({k}_{x},{k}_{y},w,t)\mid z({k}_{x},{k}_{y},p,q,r,w,t,m),A\sim \hfill \\ \hfill & CN\left(\begin{array}{c}\mathrm{sinc}(\pi {k}_{x}\u2215P)\mathrm{sinc}(\pi {k}_{y}\u2215Q)\hfill \\ \sum _{(t\u2215\delta )=0}^{T-1}\sum _{m=0}^{M-1}\sum _{n=0}^{N\left(m\right)-1}\sum _{{r}^{\prime}=wc}^{(w+1)c-1}\sum _{({p}^{\prime},{q}^{\prime},{r}^{\prime})\ne (p,q,r)}{A}_{m}({p}^{\prime},{q}^{\prime},{r}^{\prime}){L}_{mn}\mathrm{exp}\left\{\begin{array}{c}-i\left[\begin{array}{c}({\omega}_{m}+{\omega}_{mn}){B}_{0}({p}^{\prime},{q}^{\prime},r)t+{\varphi}_{mn}\hfill \\ +{\varphi}_{A}({p}^{\prime},{q}^{\prime},r)+2\pi \left({k}_{x}{p}^{\prime}\u2215P+{k}_{y}{q}^{\prime}\u2215Q\right)\hfill \end{array}\right]\hfill \\ -\left[t\u2215{T}_{a}({p}^{\prime},{q}^{\prime},r)\right]+{\left[t\u2215{T}_{b}({p}^{\prime},{q}^{\prime},r)\right]}^{2}\hfill \end{array}\right\}\hfill \end{array}+z({k}_{x},{k}_{y},p,q,r,w,t,m),\frac{{\sigma}^{2}(\nu -1)}{\nu}{I}_{2}\right).\hfill \end{array}$$

where *ν* = *cPQ*

Therefore,

$$\begin{array}{cc}\hfill & \pi \left(z({k}_{x},{k}_{y},p,q,r,w,t,m)\mid d({k}_{x},{k}_{y},w,t),A\right)=\hfill \\ \hfill & \propto \mathrm{exp}\left\{-\left[\begin{array}{c}\frac{\nu}{2{\sigma}^{2}(\nu -1)}{\mid z({k}_{x},{k}_{y},p,q,r,w,t,m)\mid}^{2}-\hfill \\ \left[\begin{array}{c}z({k}_{x},{k}_{y},p,q,r,w,t,m)\left\{\stackrel{\u2012}{\left(\begin{array}{c}d({k}_{x},{k}_{y},w,t)-\mathrm{sinc}(\pi {k}_{x}\u2215P)\mathrm{sinc}(\pi {k}_{y}\u2215Q)\sum _{m=0}^{M-1}\sum _{n=0}^{N\left(m\right)-1}\sum _{{r}^{\prime}=wc}^{(w+1)c-1}\sum _{({p}^{\prime},{q}^{\prime},{r}^{\prime})\ne (p,q,r)}{A}_{m}({p}^{\prime},{q}^{\prime},{r}^{\prime}){L}_{mn}\hfill \\ \mathrm{exp}\left\{\begin{array}{c}-i\left[({\omega}_{m}+{\omega}_{mn}){B}_{0}({p}^{\prime},{q}^{\prime},{r}^{\prime})t+{\varphi}_{mn}+{\varphi}_{A}({p}^{\prime},{q}^{\prime},{r}^{\prime})+2\pi \left({k}_{x}{p}^{\prime}\u2215P+{k}_{y}{q}^{\prime}\u2215Q\right)\right]\hfill \\ -\left[t\u2215{T}_{a}({p}^{\prime},{q}^{\prime},{r}^{\prime})\right]+{\left[t\u2215{T}_{b}({p}^{\prime},{q}^{\prime},{r}^{\prime})\right]}^{2}\hfill \end{array}\right\}\hfill \end{array}\right)}\right\}\hfill \\ +\stackrel{\u2012}{z({k}_{x},{k}_{y},p,q,r,w,t,m)}\left\{\left(\begin{array}{c}d({k}_{x},{k}_{y},w,t)-\mathrm{sinc}(\pi {k}_{x}\u2215P)\mathrm{sinc}(\pi {k}_{y}\u2215Q)\sum _{m=0}^{M-1}\sum _{n=0}^{N\left(m\right)-1}\sum _{r=wc}^{(w+1)c-1}\sum _{{p}^{\prime}\ne p,{q}^{\prime}\ne q,{r}^{\prime}\ne r}{A}_{m}({p}^{\prime},{q}^{\prime},{r}^{\prime}){L}_{mn}\hfill \\ \mathrm{exp}\left\{\begin{array}{c}-i\left[({\omega}_{m}+{\omega}_{mn}){B}_{0}({p}^{\prime},{q}^{\prime},{r}^{\prime})t+{\varphi}_{mn}+{\varphi}_{A}({p}^{\prime},{q}^{\prime},{r}^{\prime})+2\pi \left({k}_{x}{p}^{\prime}\u2215P+{k}_{y}{q}^{\prime}\u2215Q\right)\right]\hfill \\ -\left[t\u2215{T}_{a}({p}^{\prime},{q}^{\prime},{r}^{\prime})\right]+{\left[t\u2215{T}_{b}({p}^{\prime},{q}^{\prime},{r}^{\prime})\right]}^{2}\hfill \end{array}\right\}\hfill \end{array}\right)\right\}\hfill \end{array}\right]\hfill \\ +\frac{\nu}{2{\sigma}^{2}}{\mid z({k}_{x},{k}_{y},p,q,r,w,t,m)\mid}^{2}-\hfill \\ \left[\begin{array}{c}z({k}_{x},{k}_{y},p,q,r,w,t,m)\left(\stackrel{\u2012}{\mathrm{sinc}(\pi {k}_{x}\u2215P)\mathrm{sinc}(\pi {k}_{y}\u2215Q)\sum _{n=0}^{N\left(m\right)-1}{A}_{m}(p,q,r){L}_{mn}\mathrm{exp}\left\{\begin{array}{c}-i\left[\begin{array}{c}({\omega}_{m}+{\omega}_{mn}){B}_{0}(p,q,r)t+{\varphi}_{mn}\hfill \\ +{\varphi}_{A}(p,q,r)+2\pi \left({k}_{x}p\u2215P+{k}_{y}q\u2215Q\right)\hfill \end{array}\right]\hfill \\ -\left[t\u2215{T}_{a}(p,q,r)\right]+{\left[t\u2215{T}_{b}(p,q,r)\right]}^{2}\hfill \end{array}\right\}}\right)\hfill \\ +\stackrel{\u2012}{z({k}_{x},{k}_{y},p,q,r,w,t,m)}\left(\mathrm{sinc}(\pi {k}_{x}\u2215P)\mathrm{sinc}(\pi {k}_{y}\u2215Q)\sum _{n=0}^{N\left(m\right)-1}{A}_{m}(p,q,r){L}_{mn}\mathrm{exp}\left\{\begin{array}{c}-i\left[\begin{array}{c}({\omega}_{m}+{\omega}_{mn}){B}_{0}(p,q,r)t+{\varphi}_{mn}\hfill \\ +{\varphi}_{A}(p,q,r)+2\pi \left({k}_{x}p\u2215P+{k}_{y}q\u2215Q\right)\hfill \end{array}\right]\hfill \\ -\left[t\u2215{T}_{a}(p,q,r)\right]+{\left[t\u2215{T}_{b}(p,q,r)\right]}^{2}\hfill \end{array}\right\}\right)\hfill \end{array}\right]\hfill \end{array}\right]\right\}\hfill \end{array}$$

Completing the square, this leads to

$$\begin{array}{cc}\hfill & \pi \left(z({k}_{x},{k}_{y},p,q,r,w,t,m)\mid d({k}_{x},{k}_{y},w,t),A\right)\propto \hfill \\ \hfill & \text{exp}\left\{-\frac{{\nu}^{2}}{2{\sigma}^{2}(\nu -1)}{\mid \begin{array}{c}z({k}_{x},{k}_{y},p,q,r,w,t,m)-\hfill \\ \left({\scriptstyle \begin{array}{c}\frac{1}{\nu}\left[{\scriptstyle \begin{array}{c}d({k}_{x},{k}_{y},w,t)-\text{sinc}(\pi {k}_{x}\u2215P)\text{sinc}(\pi {k}_{y}\u2215Q)\sum _{m=0}^{M-1}\sum _{n=0}^{N\left(m\right)-1}{L}_{mn}\sum _{{r}^{\prime}=wc}^{(w+1)c-1}\sum _{({p}^{\prime},{q}^{\prime},{r}^{\prime})\ne (p,q,r)}{A}_{m}({p}^{\prime},{q}^{\prime},{r}^{\prime})\hfill \\ \text{exp}\left\{{\scriptstyle \begin{array}{c}-i\left[({\omega}_{m}+{\omega}_{mn}){B}_{0}({p}^{\prime},{q}^{\prime},{r}^{\prime})t+{\varphi}_{mn}+{\varphi}_{A}({p}^{\prime},{q}^{\prime},{r}^{\prime})+2\pi \left({k}_{x}{p}^{\prime}\u2215P+{k}_{y}{q}^{\prime}\u2215Q\right)\right]\hfill \\ -\left[t\u2215{T}_{a}({p}^{\prime},{q}^{\prime},{r}^{\prime})\right]+{\left[t\u2215{T}_{b}({p}^{\prime},{q}^{\prime},{r}^{\prime})\right]}^{2}\hfill \end{array}}\right\}\hfill \end{array}}\right]\hfill \\ +\frac{\nu -1}{\nu}\text{sinc}(\pi {k}_{x}\u2215P)\text{sinc}(\pi {k}_{y}\u2215Q){A}_{m}(p,q,r)\sum _{n=0}^{N\left(m\right)-1}{L}_{mn}\text{exp}\left\{{\scriptstyle \begin{array}{c}-i\left[{\scriptstyle \begin{array}{c}({\omega}_{m}+{\omega}_{mn}){B}_{0}(p,q,r)t+{\varphi}_{mn}\hfill \\ +{\varphi}_{A}(p,q,r)+2\pi \left({k}_{x}p\u2215P+{k}_{y}q\u2215Q\right)\hfill \end{array}}\right]\hfill \\ -\left[t\u2215{T}_{a}(p,q,r)\right]+{\left[t\u2215{T}_{b}(p,q,r)\right]}^{2}\hfill \end{array}}\right\}\hfill \end{array}}\right)\hfill \end{array}\mid}^{2}\right\}.\hfill \end{array}$$

Therefore, *π*(*z*(*k _{x}*,

$$\begin{array}{cc}\hfill & \frac{1}{\nu}\left[\begin{array}{c}d({k}_{x},{k}_{y},w,t)-\mathrm{sinc}(\pi {k}_{x}\u2215P)\mathrm{sinc}(\pi {k}_{y}\u2215Q)\sum _{m=0}^{M-1}{L}_{mn}\sum _{n=0}^{N\left(m\right)-1}\sum _{{r}^{\prime}=wc}^{(w+1)c-1}\sum _{({p}^{\prime},{q}^{\prime},{r}^{\prime})\ne (p,q,r)}{A}_{m}({p}^{\prime},{q}^{\prime},{r}^{\prime})\hfill \\ \mathrm{exp}\left\{-i\left[({\omega}_{m}+{\omega}_{mn}){B}_{0}({p}^{\prime},{q}^{\prime},{r}^{\prime})t+{\varphi}_{mn}+{\varphi}_{A}({p}^{\prime},{q}^{\prime},{r}^{\prime})+2\pi \left({k}_{x}{p}^{\prime}\u2215P+{k}_{y}{q}^{\prime}\u2215Q\right)\right]-\left[t\u2215{T}_{a}({p}^{\prime},{q}^{\prime},{r}^{\prime})\right]+{\left[t\u2215{T}_{b}({p}^{\prime},{q}^{\prime},{r}^{\prime})\right]}^{2}\right\}\hfill \end{array}\right]\hfill \\ \hfill & +\frac{\nu -1}{\nu}\mathrm{sinc}(\pi {k}_{x}\u2215P)\mathrm{sinc}(\pi {k}_{y}\u2215Q){A}_{m}(p,q,r)\sum _{n=0}^{N\left(m\right)-1}{L}_{mn}\mathrm{exp}\left\{-i\left[\begin{array}{c}({\omega}_{m}+{\omega}_{mn}){B}_{0}(p,q,r)t+{\varphi}_{mn}\hfill \\ +{\varphi}_{A}(p,q,r)+2\pi \left({k}_{x}p\u2215P+{k}_{y}q\u2215Q\right)\hfill \end{array}\right]-\left[t\u2215{T}_{a}(p,q,r)\right]+{\left[t\u2215{T}_{b}(p,q,r)\right]}^{2}\right\}\hfill \end{array}$$

and the *j* th E-step update for *z*(*k _{x}*,

$$\begin{array}{cc}\hfill z{({k}_{x},{k}_{y},p,q,r,w,t,m)}^{j}& =\frac{1}{\nu}\left[\begin{array}{c}d({k}_{x},{k}_{y},w,t)-\mathrm{sinc}(\pi {k}_{x}\u2215P)\mathrm{sinc}(\pi {k}_{y}\u2215Q)\sum _{m=0}^{M-1}\sum _{n=0}^{N\left(m\right)-1}{L}_{mn}\sum _{{r}^{\prime}=wc}^{(w+1)c-1}\sum _{({p}^{\prime},{q}^{\prime},{r}^{\prime})\ne (p,q,r)}{A}_{m}({p}^{\prime},{q}^{\prime},{r}^{\prime})\hfill \\ \mathrm{exp}\left\{-i\left[({\omega}_{m}+{\omega}_{mn}){B}_{0}({p}^{\prime},{q}^{\prime},{r}^{\prime})t+{\varphi}_{mn}+{\varphi}_{A}({p}^{\prime},{q}^{\prime},{r}^{\prime})+2\pi \left({k}_{x}{p}^{\prime}\u2215P+{k}_{y}{q}^{\prime}\u2215Q\right)\right]-\left[t\u2215{T}_{a}({p}^{\prime},{q}^{\prime},{r}^{\prime})\right]+{\left[t\u2215{T}_{b}({p}^{\prime},{q}^{\prime},{r}^{\prime})\right]}^{2}\right\}\hfill \end{array}\right]\hfill \\ \hfill & +\frac{\nu -1}{\nu}\mathrm{sinc}(\pi {k}_{x}\u2215P)\mathrm{sinc}(\pi {k}_{y}\u2215Q){A}_{m}(p,q,r)\sum _{n=0}^{N\left(m\right)-1}{L}_{mn}\mathrm{exp}\left\{\begin{array}{c}-i\left[\begin{array}{c}({\omega}_{m}+{\omega}_{mn}){B}_{0}(p,q,r)t+{\varphi}_{mn}+\hfill \\ {\varphi}_{A}(p,q,r)+2\pi \left({k}_{x}p\u2215P+{k}_{y}q\u2215Q\right)\hfill \end{array}\right]\hfill \\ -\left[t\u2215{T}_{a}(p,q,r)\right]+{\left[t\u2215{T}_{b}(p,q,r)\right]}^{2}\hfill \end{array}\right\}+\frac{\nu -1}{\nu}z{({k}_{x},{k}_{y},p,q,r,w,t,m)}^{j-1}\hfill \end{array}$$

To generate the M-step, the log-posterior distribution, log (*π*(*A _{m}*(

Differentiating the log-posterior distribution with respect to *A _{m}*(

$$\begin{array}{cc}\hfill & \frac{1}{{\sigma}^{2}}\left\{\begin{array}{c}\sum _{(t\u2215\delta )=0}^{T\u2215\delta -1}\sum _{{k}_{x},{k}_{y}}\left(\mathrm{sinc}(\pi {k}_{x}\u2215P)\mathrm{sinc}(\pi {k}_{y}\u2215Q){A}_{m}(p,q,r)\sum _{n=0}^{N\left(m\right)-1}{L}_{mn}\mathrm{Re}\left[\mathrm{exp}\left\{\begin{array}{c}-i\left[\begin{array}{c}({\omega}_{m}+{\omega}_{mn}){B}_{0}(p,q,r)t+{\varphi}_{mn}\hfill \\ +{\varphi}_{A}(p,q,r)+2\pi \left({k}_{x}p\u2215P+{k}_{y}q\u2215Q\right)\hfill \end{array}\right]\hfill \\ -\left[t\u2215{T}_{a}(p,q,r)\right]+{\left[t\u2215{T}_{b}(p,q,r)\right]}^{2}\hfill \end{array}\right\}\right]-\mathrm{Re}\left[z({k}_{x},{k}_{y},p,q,r,w,t,m)\right]\right)\hfill \\ \mathrm{sinc}(\pi {k}_{x}\u2215P)\mathrm{sinc}(\pi {k}_{y}\u2215Q)\mathrm{Re}\left[\mathrm{exp}\left\{-i\left[({\omega}_{m}+{\omega}_{mn}){B}_{0}(p,q,r)t+{\varphi}_{mn}+{\varphi}_{A}(p,q,r)+2\pi \left({k}_{x}p\u2215P+{k}_{y}q\u2215Q\right)\right]-\left[t\u2215{T}_{a}(p,q,r)\right]+{\left[t\u2215{T}_{b}(p,q,r)\right]}^{2}\right\}\right]\hfill \\ +\sum _{(t\u2215\delta )=0}^{T\u2215\delta -1}\sum _{{k}_{x},{k}_{y}}\left(\mathrm{sinc}(\pi {k}_{x}\u2215P)\mathrm{sinc}(\pi {k}_{y}\u2215Q){A}_{m}(p,q,r)\sum _{n=0}^{N\left(m\right)-1}{L}_{mn}\mathrm{Im}\left[\mathrm{exp}\left\{\begin{array}{c}-i\left[\begin{array}{c}({\omega}_{m}+{\omega}_{mn}){B}_{0}(p,q,r)t+{\varphi}_{mn}\hfill \\ +{\varphi}_{A}(p,q,r)+2\pi \left({k}_{x}p\u2215P+{k}_{y}q\u2215Q\right)\hfill \end{array}\right]\hfill \\ -\left[t\u2215{T}_{a}(p,q,r)\right]+{\left[t\u2215{T}_{b}(p,q,r)\right]}^{2}\hfill \end{array}\right\}\right]-\mathrm{Im}\left[z({k}_{x},{k}_{y},p,q,r,w,t,m)\right]\right)\hfill \\ \mathrm{sinc}(\pi {k}_{x}\u2215P)\mathrm{sinc}(\pi {k}_{y}\u2215Q)\mathrm{Im}\left[\mathrm{exp}\left\{-i\left[({\omega}_{m}+{\omega}_{mn}){B}_{0}(p,q,r)t+{\varphi}_{mn}+{\varphi}_{A}(p,q,r)+2\pi \left({k}_{x}p\u2215P+{k}_{y}q\u2215Q\right)\right]-\left[t\u2215{T}_{a}(p,q,r)\right]+{\left[t\u2215{T}_{b}(p,q,r)\right]}^{2}\right\}\right]\hfill \end{array}\right\}\hfill \\ \hfill & +\sum _{\langle ({p}^{\prime},{q}^{\prime},{r}^{\prime})\in \delta (p,q,r)\rangle}\left\{\left(\sum _{\Delta \in \{B,G,W\}}\frac{1}{{\tau}_{m\Delta}^{2}}{I}_{\Delta}[(p,q,r),({p}^{\prime},{q}^{\prime},{r}^{\prime})]\right)\left[{A}_{m}(p,q,r)-{A}_{m}({p}^{\prime},{q}^{\prime},{r}^{\prime})\right]\right\}=0\hfill \end{array}$$

[10]

where *δ*(*p*,*q*,*r*) is the set of neighboring voxels to voxel (*p*,*q*,*r*).

Rearranging [10] in terms of *A _{m}*(

$$\begin{array}{cc}\hfill {A}_{m}(p,q,r)=& \hfill \\ \hfill & \frac{1}{{\sigma}^{2}}\sum _{{k}_{x},{k}_{y}}\mathrm{sinc}(\pi {k}_{x}\u2215P)\mathrm{sinc}(\pi {k}_{y}\u2215Q)\sum _{(t\u2215\delta )=0}^{T-1}\mathrm{exp}\left\{-\left[t\u2215{T}_{a}(p,q,r)\right]+{\left[t\u2215{T}_{b}(p,q,r)\right]}^{2}\right\}\hfill \\ \hfill & \sum _{n=0}^{N\left(m\right)-1}{L}_{mn}\left(\begin{array}{c}\mathrm{Re}\left[z({k}_{x},{k}_{y},p,q,r,w,t,m)\right]\mathrm{cos}\left(\begin{array}{c}({\omega}_{m}+{\omega}_{mn}){B}_{0}(p,q,r)t+{\varphi}_{mn}\hfill \\ +{\varphi}_{A}(p,q,r)+2\pi \left({k}_{x}p\u2215P+{k}_{y}q\u2215Q\right)\hfill \end{array}\right)-\hfill \\ \mathrm{Im}\left[z({k}_{x},{k}_{y},p,q,r,w,t,m)\right]\mathrm{sin}\left(\begin{array}{c}({\omega}_{m}+{\omega}_{mn}){B}_{0}(p,q,r)t+{\varphi}_{mn}\hfill \\ +{\varphi}_{A}(p,q,r)+2\pi \left({k}_{x}p\u2215P+{k}_{y}q\u2215Q\right)\hfill \end{array}\right)\hfill \end{array}\right)\hfill \\ \hfill & \frac{+\sum _{\langle ({p}^{\prime},{q}^{\prime},{r}^{\prime})\in \delta (p,q,r)\rangle}\left(\sum _{\Delta \in \{B,G,W\}}\frac{1}{{\tau}_{m\Delta}^{2}}{I}_{\Delta}[(p,q,r),({p}^{\prime},{q}^{\prime},{r}^{\prime})]\right){A}_{m}({p}^{\prime},{q}^{\prime},{r}^{\prime})}{\frac{1}{{\sigma}^{2}}\sum _{{k}_{x},{k}_{y}}{\left(\mathrm{sinc}(\pi {k}_{x}\u2215P)\mathrm{sinc}(\pi {k}_{y}\u2215Q)\right)}^{2}+\sum _{\langle ({p}^{\prime},{q}^{\prime},{r}^{\prime})\in \delta (p,q,r)\rangle}\sum _{\Delta \in \{B,G,W\}}\frac{1}{{\tau}_{m\Delta}^{2}}{I}_{\Delta}[(p,q,r,),({p}^{\prime},{q}^{\prime},{r}^{\prime})]}\hfill \end{array}$$

[11]

Differentiating [10] again with respect to *A _{m}*(

$$\begin{array}{cc}\hfill & \frac{1}{{\sigma}^{2}}\left\{\begin{array}{c}\sum _{{k}_{x},{k}_{y}}{\left(\mathrm{sinc}(\pi {k}_{x}\u2215P)\mathrm{sinc}(\pi {k}_{y}\u2215Q)\sum _{n=0}^{N\left(m\right)-1}{L}_{mn}\mathrm{Re}\left[\mathrm{exp}\left\{\begin{array}{c}-i\left[({\omega}_{m}+{\omega}_{mn}){B}_{0}(p,q,r)t+{\varphi}_{mn}+{\varphi}_{A}(p,q,r)+2\pi \left({k}_{x}p\u2215P+{k}_{y}q\u2215Q\right)\right]\hfill \\ -\left[t\u2215{T}_{a}(p,q,r)\right]+{\left[t\u2215{T}_{b}(p,q,r)\right]}^{2}\hfill \end{array}\right\}\right]\right)}^{2}\hfill \\ \sum _{{k}_{x},{k}_{y}}{\left(\mathrm{sinc}(\pi {k}_{x}\u2215P)\mathrm{sinc}(\pi {k}_{y}\u2215Q)\sum _{n=0}^{N\left(m\right)-1}{L}_{mn}\mathrm{Im}\left[\mathrm{exp}\left\{\begin{array}{c}-i\left[({\omega}_{m}+{\omega}_{mn}){B}_{0}(p,q,r)t+{\varphi}_{mn}+{\varphi}_{A}(p,q,r)+2\pi \left({k}_{x}p\u2215P+{k}_{y}q\u2215Q\right)\right]\hfill \\ -\left[t\u2215{T}_{a}(p,q,r)\right]+{\left[t\u2215{T}_{b}(p,q,r)\right]}^{2}\hfill \end{array}\right\}\right]\right)}^{2}\hfill \end{array}\right\}\hfill \\ \hfill & +\sum _{\langle ({p}^{\prime},{q}^{\prime},{r}^{\prime})\in \delta (p,q,r)\rangle}\sum _{\Delta \in \{B,G,W\}}\frac{1}{{\tau}_{m\Delta}^{2}}{I}_{\Delta}[(p,q,r),({p}^{\prime},{q}^{\prime},{r}^{\prime})]>0\hfill \end{array}$$

[12]

The inequality of Equation [12] confirms that the solution given by [11] is a local maximum and therefore the *j* th M-step update for *A _{m}*(

$$\begin{array}{cc}\hfill {A}_{m}{(p,q,r)}^{j}=& \hfill \\ \hfill & \frac{1}{{\sigma}^{2}}\sum _{{k}_{x},{k}_{y}}\mathrm{sinc}(\pi {k}_{x}\u2215P)\mathrm{sinc}(\pi {k}_{y}\u2215Q)\sum _{(t\u2215\delta )=0}^{T-1}\mathrm{exp}\left\{-\left[t\u2215{T}_{a}(p,q,r)\right]+{\left[t\u2215{T}_{b}(p,q,r)\right]}^{2}\right\}\hfill \\ \hfill & \sum _{n=0}^{N\left(m\right)-1}{L}_{mn}\left(\begin{array}{c}\mathrm{Re}\left[z{({k}_{x},{k}_{y},p,q,r,w,t,m)}^{j}\right]\mathrm{cos}\left(\begin{array}{c}({\omega}_{m}+{\omega}_{mn}){B}_{0}(p,q,r)t+{\varphi}_{mn}\hfill \\ +{\varphi}_{A}(p,q,r)+2\pi \left({k}_{x}p\u2215P+{k}_{y}q\u2215Q\right)\hfill \end{array}\right)-\hfill \\ \mathrm{Im}\left[z{({k}_{x},{k}_{y},p,q,r,w,t,m)}^{j}\right]\mathrm{sin}\left(\begin{array}{c}({\omega}_{m}+{\omega}_{mn}){B}_{0}(p,q,r)t+{\varphi}_{mn}\hfill \\ +{\varphi}_{A}(p,q,r)+2\pi \left({k}_{x}p\u2215P+{k}_{y}q\u2215Q\right)\hfill \end{array}\right)\hfill \end{array}\right)\hfill \\ \hfill & \frac{+\sum _{\langle ({p}^{\prime},{q}^{\prime},{r}^{\prime})\in \delta (p,q,r)\rangle}\left(\frac{1}{{\tau}_{A}^{2}}+\frac{1}{{\tau}_{G}^{2}}{I}_{G}[(p,q,r),({p}^{\prime},{q}^{\prime},{r}^{\prime})]+\frac{1}{{\tau}_{W}^{2}}{I}_{W}[(p,q,r),({p}^{\prime},{q}^{\prime},{r}^{\prime})]\right){A}_{m}{(p,q,r)}^{j\ast}}{\frac{1}{{\sigma}^{2}}\sum _{{k}_{x},{k}_{y}}{\left(\mathrm{sinc}(\pi {k}_{x}\u2215P)\mathrm{sinc}(\pi {k}_{y}\u2215Q)\right)}^{2}+\sum _{\langle ({p}^{\prime},{q}^{\prime},{r}^{\prime})\in \delta (p,q,r)\rangle}\sum _{\Delta \in \{B,G,W\}}\frac{1}{{\tau}_{m\Delta}^{2}}{I}_{\Delta}[(p,q,r),({p}^{\prime},{q}^{\prime},{r}^{\prime})]}\hfill \end{array}$$

[13]

1. Brown TR, Kincaid BM, Ugurbil K. NMR chemical shift imaging in three dimensions. Proceedings of the National Academy of Sciences. 1982;79:3523–3526. [PubMed]

2. Maudsley AA, Hilal SK, Perman WH, Simon HE. Spatially resolved high resolution spectroscopy by four-dimensional NMR. J Magn Reson. 1983;51:147–152.

3. Young K, Govindaraju V, Soher BJ, Maudsley AA. Automated spectral analysis I: formation of a priori information by spectral simulation. Magn Reson Med. 1998 Dec;40:812–5. [PubMed]

4. Young K, Soher BJ, Maudsley AA. Automated spectral analysis II: application of wavelet shrinkage for characterization of non-parameterized signals. Magn Reson Med. 1998 Dec;40:816–21. [PubMed]

5. Soher BJ, Young K, Govindaraju V, Maudsley AA. Automated spectral analysis III: application to in vivo proton MR spectroscopy and spectroscopic imaging. Magn Reson Med. 1998 Dec;40:822–31. [PubMed]

6. de Beer R, van den Boogaart A, van Ormondt D, Pijnappel WW, den Hollander JA, Marien AJ, Luyten PR. Application of time-domain fitting in the quantification of in vivo 1H spectroscopic imaging data sets. NMR Biomed. 1992;5:171–178. [PubMed]

7. Liang Z-P, Lauterbur PC, IEEE Engineering in Medicine and Biology Society . Principles of magnetic resonance imaging : a signal processing perspective. SPIE Optical Engineering Press; IEEE Press; Bellingham, Wash. New York: 2000.

8. Hurn MA, Mardia KV, Hainsworth TJ, Kirkbridge J, Berry E. Bayesian Fused Classification of Medical Images. IEEE Transactions on Medical Imaging. 1996;15:850–858. [PubMed]

9. Miller IM, Schaewe TJ, Cohen SB, Ackerman JJH. Model-Based Maximum-Likelihood Estimation for Phase- and Frequency-Encoded Magnetic-Resonance-Imaging Data. Journal of Magnetic Resonance, Series B. 1995;107:10–221. [PubMed]

10. Schaewe TJ, Miller IM. Parallel Algorithms for Maximum A Posteriori Estimation of Spin Density and Spin-spin Decay in Magnetic Resonance Imaging. IEEE Transactions on Medical Imaging. 1995:362–373. [PubMed]

11. Denney TS, Jr., Reeves SJ. Bayesian image reconstruction from Fourier-domain samples using prior edge information. Journal of Electronic Imaging. 2005;14:043009, 1–11. 2005.

12. Ouyang X, Wong WH, Johnson VE, Hu X, Chen C-T. Incorporation of Correlated Structural Images in PET Image Reconstruction. IEEE Transactions on Medical Imaging. 1994;13:627–640. [PubMed]

13. Bowsher JE, Johnson VE, Turkington TG, Jaszczak RJ, Floyd CE, Jr., Coleman RE. Bayesian Reconstruction and Use of Anatomical A Priori Information for Emission Tomography. IEEE Transactions on Medical Imaging. 1996;15:673–686. [PubMed]

14. Sastry S, Carson RE. Multimodality Bayesian Algorithm for Image Reconstruction in Positron Emission Tomography: A Tissue Composition Model. IEEE Transactions on Medical Imaging. 1997;16:750–761. [PubMed]

15. Hu X, Levin DN, Lauterbur PC, Spraggins T. SLIM: spectral localization by imaging. Magn Reson Med. 1988;8:314–22. [PubMed]

16. Von Kienlin M, Mejia R. Spectral localization with optimal pointspread function. Journal of magnetic resonance. 1991;94:268–287.

17. Liang ZP, Lauterbur PC. A generalized series approach to MR spectroscopic imaging. Medical Imaging, IEEE Transactions on. 1991;10:132–137. [PubMed]

18. Liang ZP, Lauterbur PC. An efficient method for dynamic magnetic resonance imaging. Medical Imaging, IEEE Transactions on. 1994;13:677–686. [PubMed]

19. Jacob M, Zhu X, Ebel A, Schuff N, Liang ZP. Improved Model-Based Magnetic Resonance Spectroscopic Imaging. Medical Imaging, IEEE Transactions on. 2007;26:1305–1318. [PubMed]

20. Bao Y, Maudsley AA. Improved Reconstruction for MR Spectroscopic Imaging. Medical Imaging, IEEE Transactions on. 2007;26:686–695. [PMC free article] [PubMed]

21. Haldar JP, Hernando D, Song SK, Liang ZP. Anatomically constrained reconstruction from noisy data. Magnetic Resonance in Medicine. 2008;59:810. [PubMed]

22. Kornak J, Young K, Schuff N, Du A, Maudsley AA, Weiner MW. K-Bayes Reconstruction for Perfusion MRI I: Concepts and Application. Journal of Digital Imaging. 2008 [PMC free article] [PubMed]

23. Kornak J, Young K. K-Bayes Reconstruction for Perfusion MRI II: Modeling and Technical Development. Journal of Digital Imaging. 2008 [PMC free article] [PubMed]

24. Geman S, German D. Stochastic Relaxation, Gibbs Distributions, and the Bayesian Restoration of Images. IEEE Transactions on Pattern Analysis and Machine Intelligence. 1984;6:721–741. [PubMed]

25. Besag JE. Towards Bayesian Image Analysis. Journal of Applied Statistics. 1989;16:395–407.

26. Winkler G. Image Analysis, Random Fields, and Markov Chain Monte Carlo Methods: A Mathematical Introduction. 2nd ed. Springer; 2003.

27. Maudsley A, Domenig C, Govind V, Darkazanli A, Studholme C, Arheart K, Bloomer C. Mapping of brain metabolite distributions by volumetric proton MR spectroscopic imaging (MRSI) Magnetic Resonance in Medicine. 2009;61 [PMC free article] [PubMed]

28. Besag JE. Spatial Interaction and The Statistical Analysis of Lattice Systems (with discussion) Journal of the Royal Statistical Society, Series B. 1974;36:192–236.

29. Dempster AP, Laird NM, Rubin DB. Maximum Likelihood From Incomplete Data Via the EM Algorithm (with discussion) Journal of the Royal Statistical Society, Series B. 1977;39:1–38.

30. Ashburner J, Friston KJ. Multimodal Image Coregistration and Partitioning - A Unified Framework. NeuroImage. 1997;6:209–217. [PubMed]

31. Zhang Y, Brady M, Smith S. Segmentation of brain MR images through a hidden Markov random field model and the expectation-maximization algorithm. IEEE Trans Med Imaging. 2001 Jan;20:45–57. [PubMed]

32. Besag JE. Statistical Analysis of Non-lattice Data. The Statistician. 1975;24:179–195.

33. Kunsch HR. Intrinsic Autoregressions and Related Models on the Two-dimensional Lattice. Biometrika. 1987;74:517–524.

34. Besag JE, Kooperberg C. On Conditional and Intrinsic Autoregressions. Biometrika. 1995;82:733–746.

35. Green PJ. Bayesian Reconstructions From Emission Tomography Data Using A Modified EM Algorithm. IEEE Transactions on Medical Imaging. 1990;9:84–93. [PubMed]

36. Cohen L. Time-frequency analysis: theory and applications. 1995.

37. Cocosco CA, Kollokian V, Kwan RK-S, Evans AC. BrainWeb: Online Interface to a 3D MRI Simulated Brain Database. Proceedings of 3rd International Conference on Functional Mapping of the Human Brain; 1997.

38. Maudsley AA, Lin E, Weiner MW. Spectroscopic imaging display and analysis. Magn Reson Imaging. 1992;10:471–85. [PubMed]

39. Kaiser LG, Young K, Matson GB. Elimination of spatial interference in PRESS-localized editing spectroscopy. Magnetic resonance in medicine: official journal of the Society of Magnetic Resonance in Medicine/Society of Magnetic Resonance in Medicine. 2007;58:813. [PubMed]

40. Kaiser LG, Young K, Matson GB. Numerical simulations of localized high field 1H MR spectroscopy. Journal of Magnetic Resonance. 2008;195:67–75. [PMC free article] [PubMed]

41. Provencher S. Automatic quantitation of localized in vivo 1H spectra with LCModel. NMR in Biomedicine. 2001;14 [PubMed]

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |