Optimization of magnetic flux density for fast MREIT conductivity imaging using multi-echo interleaved partial fourier acquisitions

Abstract

Background

Magnetic resonance electrical impedance tomography (MREIT) has been introduced as a non-invasive method for visualizing the internal conductivity and/or current density of an electrically conductive object by externally injected currents. The injected current through a pair of surface electrodes induces a magnetic flux density distribution inside the imaging object, which results in additional magnetic flux density. To measure the magnetic flux density signal in MREIT, the phase difference approach in an interleaved encoding scheme cancels out the systematic artifacts accumulated in phase signals and also reduces the random noise effect by doubling the measured magnetic flux density signal. For practical applications of in vivo MREIT, it is essential to reduce the scan duration maintaining spatial-resolution and sufficient contrast. In this paper, we optimize the magnetic flux density by using a fast gradient multi-echo MR pulse sequence. To recover the one component of magnetic flux density B _z , we use a coupled partial Fourier acquisitions in the interleaved sense.

Methods

To prove the proposed algorithm, we performed numerical simulations using a two-dimensional finite-element model. For a real experiment, we designed a phantom filled with a calibrated saline solution and located a rubber balloon inside the phantom. The rubber balloon was inflated by injecting the same saline solution during the MREIT imaging. We used the multi-echo fast low angle shot (FLASH) MR pulse sequence for MRI scan, which allows the reduction of measuring time without a substantial loss in image quality.

Results

Under the assumption of a priori phase artifact map from a reference scan, we rigorously investigated the convergence ratio of the proposed method, which was closely related with the number of measured phase encode set and the frequency range of the background field inhomogeneity. In the phantom experiment with a partial Fourier acquisition, the total scan time was less than 6 seconds to measure the magnetic flux density B _z data with 128×128 spacial matrix size, where it required 10.24 seconds to fill the complete k-space region.

Conclusion

Numerical simulation and experimental results demonstrated that the proposed method reduces the scanning time and provides the recovered B _z data comparable to what we obtained by measuring complete k-space data.

Background

Magnetic resonance electrical impedance tomography (MREIT) utilizes a magnetic resonance imaging (MRI) scanner to measure magnetic flux density B _z data inside an imaging object induced by the externally injected current. The internal current density distribution has been studied in magnetic resonance current density imaging (MRCDI) by measuring the whole magnetic flux density data B=(B _x ,B _y ,B _z ) [1, 2]. Combining MRCDI and electrical impedance tomography (EIT) technique, MREIT provides the cross-sectional conductivity images of the object with high spatial resolution [3–10]. Since an MRI scanner measures only one component B _z of B without rotating the imaging object, most MREIT algorithms assumed that the internal conductivity is isotropic and focused on visualizing its distribution by using one component of the magnetic flux density data B _z of B[11–18].

Recent MREIT imaging techniques have been developed with respect to both the capacity of measurement techniques and the numerical reconstruction algorithms. Experimental results from in vivo animal and human have been reported [19, 20] in MREIT. As an innovation of current MREIT, a fast MREIT imaging technique referring to the continuous monitoring of objects includes various wide application areas [21]. Recently, one of challenging problem in MREIT is to implement a new imaging technique with a very short acquisition time for the imaging of neural activities of brain related to the conductivity change.

Since current MREIT experiments suffer from poor SNR of the measured B _z data under the typical data acquisition durations and a small amount of injected current, it is important to reduce the scan time, while maintaining the spatial-resolution and sufficient contrast, for practical implementations of in vivo MREIT. Recently, to reduce the scan time in MREIT, Hamamura et al[22] reconstructed the interior conductivity using a single-shot spin-echo echo planar imaging (SS-SEPI) pulse sequence and Muftuler et al[23] used a SENSE-accelerated imaging technique to acquire phase signal by the injected current. One of basic approaches for maintaining the spatial resolution is to reduce the number of phase encoding steps because each phase encoding step requires a certain amount of time for execution. Since the MREIT techniques use an interleaved phase encoding acquisition scheme to double the B _z signal, Park et al[24] reconstructed the phase signal B _z by filling the skipped k-space region using the interleaved measurement property.

To obtain the static conductivity image in MREIT [19, 20], a spin-echo based MREIT pulse sequence has been predominantly used to reduce the background artifact and to increase the imaging quality. In real situations, it is difficult to employ the fast conventional MR pulse sequences because the noise standard deviation of B _z is inversely proportional to the width of injection current and the intensity of MR magnitude, simultaneously. An MREIT pulse sequence should be devised to enhance changes in MR phase images for given current amplitudes.

In this paper, we used the multi-echo fast low angle shot (FLASH) MR pulse sequence which allows the reduction of imaging time without any substantial loss in image quality. In addition, the multi-echo FLASH sequence maximizes the width of injection current extending the duration of injection current until the end of a readout gradient in MREIT. To reconstruct the internal conductivity distribution, most algorithms require at least two independent injection currents in an interleaved sense, which require relatively a long scanning duration [7, 11]. To reduce the scanning time, we adopt a partial phase encoding acquisition scheme using the multi-echo FLASH MREIT pulse sequence and rigorously investigate the relationship between the convergence ratio of the algorithm and the background field inhomogeneity [24]. We consider the discrete ℓ²-norm to evaluate the convergence ratio.

To show the feasibility of the proposed algorithm, we performed numerical simulations and compared the performance to the simulated true B _z data. We designed a cylindrical acrylic phantom filled with a calibrated saline solution and located a rubber balloon inside the phantom. The rubber balloon was inflated by injecting the same saline solution during the scan. The phantom was designed to provide a homogeneous magnitude image, but distinguishable signals of measured B _z between the inside and outside the balloon. The phantom experiment demonstrated that the proposed method reduces the scanning time and recovers the reasonable resolution of B _z , which is comparable to the recovered B _z using the complete k-space data.

Method

k-space signal and B _z data

In a conventional spin echo MREIT pulse sequence, both positive and negative currents of the same amplitude and duration are injected with reverse polarity. These injection currents with the pulse width of T _c accumulate extra phases. Corresponding k-space MR signals can be described as

$$S^{\pm }(k_x,k_y)={\int }_{\Omega }\rho (x,y)e^{\mathrm{i\delta }(x,y)}{e}^{\pm \mathrm{i\gamma }{B}_z(x,y)T_c}{e}^{i2\pi (k_{x}x+k_{y}y)}\mathit{\text{dxdy}}$$

(1)

where ρ is the T₂ weighted spin density, δ is any systematic phase artifact, and Ω is a field-of-view (FOV). Here, the superscript of S^±(k _x ,k _y ) denotes a brief notation for S⁺(k _x ,k _y ) and S⁻(k _x ,k _y ). For the standard coverage of k-space, we set

$$\left\{\begin{array}{l}{k}_x=\frac{\gamma }{2\pi }{G}_x(\mathrm{n\Delta t}-T_E)\text{for}n=-N_x/2,\cdots ,N_x/2-1\\ k_y=\frac{\gamma }{2\pi }\mathrm{m\Delta }{G}_y{T}_{\mathit{\text{pe}}}\text{for}m=-N_y/2,\cdots ,N_y/2-1\end{array}\right.$$

(2)

where γ=26.75×10⁷rad/T·s is the gyromagnetic ratio of hydrogen, Δ t is the time between samplings, G _x is the frequency encoding gradient strength, T _E is the echo time, Δ G _y is the phase encoding step, and T _pe is the phase encoding time. The induced magnetic flux density ±B _z is generated by the positive and negative injection currents I^±. Applying the inverse Fourier transform to the measured k-space data sets in (1), we can compute the magnetic flux density B _z as

$$B_z\left(\mathbf{r}\right)=\frac{1}{2\gamma T_c}{tan}^{-1}\left(\frac{\alpha \left(\mathbf{r}\right)}{\beta \left(\mathbf{r}\right)}\right)$$

(3)

where α and β are the imaginary and real part of $\left(\rho e^{\mathrm{i\delta }}{e}^{\mathrm{i\gamma }{B}_z{T}_c}\right)/\left(\rho e^{\mathrm{i\delta }}{e}^{-\mathrm{i\gamma }{B}_z{T}_c}\right)$, respectively [2].

The current MREIT method is based on the electromagnetic information embedded in the measured B _z data in order to visualize the conductivity(or current density) based on Biot-Savart law:

$$B_z\left(\mathbf{r}\right)=\frac{{\mu }_0}{4\pi }{\int }_{\Omega }\frac{(y-y^{\prime })J_x\left({\mathbf{r}}^{\prime }\right)-(x-x^{\prime })J_y\left({\mathbf{r}}^{\prime }\right)}{|\mathbf{r}-{\mathbf{r}}^{\prime }{|}^3}d{\mathbf{r}}^{\prime },\mathbf{r}=(x,y,z),{\mathbf{r}}^{\prime }=(x^{\prime },y^{\prime },z^{\prime })$$

(4)

where μ₀=4π 10⁻⁷Tm/A is the magnetic permeability of the free space. The current density J in Ω is given for the isotropic conductivity σ in Ω

$$\mathbf{J}\left(\mathbf{r}\right)=-\sigma \nabla u\left(\mathbf{r}\right)$$

(5)

and satisfies the following elliptic equation

$$\begin{array}{l}\nabla ·\left(\sigma \nabla u\right)=0\text{in}\Omega \\ -\sigma \nabla u·\nu =g\text{on}\mathrm{\partial \Omega }\text{and}{\int }_{\mathrm{\partial \Omega }}u\mathit{\text{ds}}=0\end{array}$$

(6)

where ν is the outward unit normal vector on ∂ Ω and g is an applied current density on the surface.

Recovery of complex $T_2^{\ast }$ weighted spin density using interleaved partial Fourier acquisition

To simplify, we develop a theory using a conventional cartesian k-space that has three zones within the k _y (phase-encode) domain; the central region (${\mathcal{P}}_0$), the positive region (${\mathcal{P}}^{+}$) and the negative region (${\mathcal{P}}^{-}$):

$$\left\{\begin{array}{l}{\mathcal{P}}^0=\left\{\right(k_x,k_y\left)\left|k_y=\frac{\gamma }{2\pi }\mathrm{m\Delta }{G}_y{T}_{\mathit{\text{pe}}},-N_c\le m\le N_c\right.\right\}\\ {\mathcal{P}}^{+}=\left\{\right(k_x,k_y\left)\left|k_y=\frac{\gamma }{2\pi }\mathrm{m\Delta }{G}_y{T}_{\mathit{\text{pe}}},m>N_c\right.\right\}\hfill \\ {\mathcal{P}}^{-}=\left\{\right(k_x,k_y\left)\left|k_y=\frac{\gamma }{2\pi }\mathrm{m\Delta }{G}_y{T}_{\mathit{\text{pe}}},m<-N_c\right.\right\}\end{array}\right.$$

(7)

Here, N _c denotes the number of partial Fourier over-sampling phase-encodes. The interleaved k-space data $S_p^{+}(k_x,k_y)$ and $S_p^{-}(k_x,k_y)$ as partially acquired for the phase $k_y=\frac{\gamma }{2\pi }\mathrm{m\Delta }{G}_y{T}_{\mathit{\text{pe}}}$ can be expressed

$$\begin{array}{ccc}{S}_p^{\pm }(k_x,k_y)& =& \left\{\begin{array}{ll}{S}^{\pm }(k_x,k_y),& (k_x,k_y)\in {\mathcal{P}}^{+}\cup {\mathcal{P}}^0\\ 0,& (k_x,k_y)\in {\mathcal{P}}^{-}\end{array}\right.\end{array}$$

(8)

where $S_p^{\pm }$ represents $S_p^{+}$ and $S_p^{-}$ simultaneously. We propose an algorithm to determine the $T_2^{\ast }$ weighted spin density ρ :

$${\rho }^{\pm }\left(\mathbf{r}\right):=\rho \left(\mathbf{r}\right)e^{\mathrm{i\delta }\left(\mathbf{r}\right)}{e}^{\pm \mathrm{i\gamma }{B}_z\left(\mathbf{r}\right)T_c}=F{T}^{-1}\left(S^{\pm }\right(k_x,k_y\left)\right)\left(\mathbf{r}\right)$$

(9)

using the partially scanned k-space data $S_p^{\pm }$.

The systematic phase artifact e^iδ(x,y), unavoidable artifacts due to the main field inhomogeneity and the mismatch between the center of data acquisition interval and echo formation, arises from a low frequency field, which mainly belongs to the central region ${\mathcal{P}}^0$. Including a small perturbed phase artifact e^iδ(x,y), we start with the initial guess $S_p^{\pm }$ and design an alternating procedure by updating the skipped k-space regions.

The recovered ${\rho }_p^{\pm }=F{T}^{-1}\left(S_p^{\pm }\right)$ can be formally expressed by the equation

$$\begin{array}{l}{\rho }_p^{\pm }\left(\mathbf{r}\right)=F{T}^{-1}\left(S^{\pm }\right(k_x,k_y\left)\right)\left(\mathbf{r}\right)+F{T}_{{\mathcal{P}}^{-}}^{-1}\left(S_p^{\pm }\right(k_x,k_y\left)\right)\left(\mathbf{r}\right)-F{T}_{{\mathcal{P}}^{-}}^{-1}\left(S^{\pm }\right(k_x,k_y\left)\right)\left(\mathbf{r}\right)\\ ={\rho }^{\pm }\left(\mathbf{r}\right)+E_p^{\pm }\left(\mathbf{r}\right)\end{array}$$

(10)

where $F{T}_{{\mathcal{P}}^{-}}^{-1}\left(S_p^{\pm }\right(k_x,k_y\left)\right)$ denotes the inverse Fourier transform by zero-filling the k-space except in the region ${\mathcal{P}}^{-}$ and the remainder term $E_p^{\pm }$ is

$$E_p^{\pm }\left(\mathbf{r}\right):=F{T}_{{\mathcal{P}}^{-}}^{-1}\left(S_p^{\pm }(k_x,k_y)-S^{\pm }(k_x,k_y)\right)\left(\mathbf{r}\right).$$

(11)

Let us define a support region ${\mathcal{D}}_{\delta }$ for a low spatially varying magnetic field due to background field inhomogeneities:

$${\mathcal{D}}_{\delta }:=\left\{\right(k_x,k_y\left)\right|\mathit{\text{FT}}\left(e^{-\mathrm{i\delta }}\right)(k_x,k_y)=0,k_y=\frac{\gamma }{2\pi }\mathrm{m\Delta }{G}_y{T}_{\mathit{\text{pe}}},-N_{\delta }\le m\le N_{\delta }\}$$

(12)

Observation 1. If ${\mathcal{D}}_{\delta }\subset {\mathcal{P}}^0$, $\mathit{\text{FT}}\left({\rho }_p^{\mp }{e}^{-2\mathrm{i\delta }}\right)(k_x,k_y)=\mathit{\text{FT}}\left({\rho }^{\mp }{e}^{-2\mathrm{i\delta }}\right)(k_x,k_y)$ for $(k_x,k_y)\in {\mathcal{P}}^{+}$.

The proof of observation 1 is provided in the Appendix A. By using the observation 1, we fill the skipped k-space regions in $S_p^{\pm }$

$$\begin{array}{lll}{S}_u^{\pm }(k_x,k_y)& =& \left\{\begin{array}{ll}{S}^{\pm }(k_x,k_y),& (k_x,k_y)\in {\mathcal{P}}^{+}\cup {\mathcal{P}}^0\\ \overline{\mathit{\text{FT}}\left({\rho }_p^{\mp }{e}^{-2\mathrm{i\delta }}\right)}(-k_x,-k_y),& (k_x,k_y)\in {\mathcal{P}}^{-}\end{array}\right.\end{array}$$

(13)

Observation 2. If ${\mathcal{D}}_{\delta }\subset {\mathcal{P}}^0$, the k-space $S_u^{\pm }$ in recovers the T₂ (or $T_2^{\ast }$) weighted spin density ρ^± without loss of information.

$$F{T}^{-1}\left(S_u^{\pm }\right)={\rho }^{\pm }=\rho e^{\mathrm{i\delta }}{e}^{\pm \mathrm{i\gamma }{B}_z(x,y)T_c}$$

The proof of observation 2 is provided in the Appendix 2. The observation 2 shows that the skipped region in the measured k-space S⁺ can be recovered by using the interleaved acquired S⁻ and the estimated background sensitivity map.

Convergence characteristics

When the common measured ${\mathcal{P}}^0$ does not cover the support region ${\mathcal{D}}_{\delta }$, N _δ >N _c , for the phase-encode $(k_x,k_y)\in {\mathcal{P}}^{+}$, the Fourier transform of ${\rho }_p^{-}{e}^{-2\mathrm{i\delta }}$ can be written by following the observation 1:

$$\mathit{\text{FT}}\left({\rho }_p^{-}{e}^{-2\mathrm{i\delta }}\right)(k_x,k_y)=\sum _{(k_x,k_m)\in {\mathcal{P}}^{+}\cup {\mathcal{P}}^0}\mathit{\text{FT}}\left({\rho }^{-}\right)(k_x,k_m)\mathit{\text{FT}}\left(e^{-2\mathrm{i\delta }}\right)(k_x,k_y-k_m)$$

(14)

From the relation (14), for the phase-encode $(k_x,k_y)\in {\mathcal{P}}^{+}$, we have

$$\begin{array}{l}\mathit{\text{FT}}\left({\rho }^{-}{e}^{-2\mathrm{i\delta }}\right)(k_x,k_y)-\mathit{\text{FT}}\left({\rho }_p^{-}{e}^{-2\mathrm{i\delta }}\right)(k_x,k_y)\\ ={\sum }_{(k_x,k_m)\in {\mathcal{P}}^{-}}\mathit{\text{FT}}\left({\rho }^{-}\right)(k_x,k_m)\mathit{\text{FT}}\left(e^{-2\mathrm{i\delta }}\right)(k_x,k_y-k_m)\end{array}$$

(15)

We set

$${\left[\psi \right]}_{\mathcal{R}}:=F{T}^{-1}\left(\mathit{\text{FT}}\left(\psi \right){|}_{\mathcal{R}}\right)$$

(16)

where $\mathcal{R}$ is a subregion of the k-space and $\mathit{\text{FT}}\left(\psi \right){|}_{\mathcal{R}}$ denotes the restriction of F T(ψ) to the region $\mathcal{R}$. The discrete ℓ₂-norm of ${\left[\psi \right]}_{\mathcal{R}}$ is equivalent to that of $\mathit{\text{FT}}\left(\psi \right){|}_{\mathcal{R}}$, i.e., $\parallel {\left[\psi \right]}_{\mathcal{R}}{\parallel }_{{\ell }_2}=C\parallel \mathit{\text{FT}}\left(\psi \right){|}_{\mathcal{R}}{\parallel }_{{\ell }_2}$, where the constant C is independent to the function ψ and the region $\mathcal{R}$.

When the skipped k-space region, ${\mathcal{P}}^{-}$, are filled the previously updated as in (13), we have

$$\begin{array}{l}\parallel {({\rho }^{-}-{\rho }_n^{-})}_{{\ell }_2}=C\parallel \mathit{\text{FT}}({\rho }^{-}-{\rho }_n^{-}){\parallel }_{{\ell }_2}\\ =C\parallel \mathit{\text{FT}}({\rho }^{-}-{\rho }_n^{-}){|}_{{\mathcal{P}}^{-}}{\parallel }_{{\ell }_2}\\ =C{∥\left(\overline{\mathit{\text{FT}}\left({\rho }^{+}{e}^{-2\mathrm{i\delta }}\right)}-\overline{\mathit{\text{FT}}\left({\rho }_{n-1}^{+}{e}^{-2\mathrm{i\delta }}\right)}\right){|}_{{\mathcal{P}}^{-}}(-k_x,-k_y)∥}_{{\ell }_2}\\ =C{∥\left(\mathit{\text{FT}}\left({\rho }^{+}{e}^{-2\mathrm{i\delta }}\right)-\mathit{\text{FT}}\left({\rho }_{n-1}^{+}{e}^{-2\mathrm{i\delta }}\right)\right){|}_{{\mathcal{P}}^{+}}(k_x,k_y)∥}_{{\ell }_2}\\ =C{∥\left(\mathit{\text{FT}}\left(\right({\rho }^{+}-{\rho }_{n-1}^{+}\left)e^{-2\mathrm{i\delta }}\right)\right){|}_{{\mathcal{P}}^{+}}(k_x,k_y)∥}_{{\ell }_2}\end{array}$$

(17)

For the k-space region $(k_x,k_y)\in {\mathcal{P}}^{+}$, we have the following identity

$$\begin{array}{l}\mathit{\text{FT}}\left(\left({\rho }^{+}-{\rho }_{n-1}^{+}\left)e^{-2\mathrm{i\delta }}\right)\right){|}_{{\mathcal{P}}^{+}}\right(k_x,k_y)\\ ={\sum }_{(k_x,k_m)\in {\mathcal{P}}^{-}}\mathit{\text{FT}}({\rho }^{+}-{\rho }_{n-1}^{+})(k_x,k_m)\mathit{\text{FT}}\left(e^{-2\mathrm{i\delta }}\right)(k_x,k_y-k_m)\\ +\Theta ({\rho }^{+}-{\rho }_{n-1}^{+},{\mathcal{P}}^0\cup {\mathcal{P}}^{+})\end{array}$$

(18)

where

$$\Theta ({\rho }^{+}-{\rho }_{n-1}^{+},{\mathcal{P}}^0\cup {\mathcal{P}}^{+}):=\sum _{(k_x,k_m)\in {\mathcal{P}}^0\cup {\mathcal{P}}^{\mathrm{+}}}\mathit{\text{FT}}({\rho }^{+}-{\rho }_{n-1}^{+})(k_x,k_m)\mathit{\text{FT}}\left(e^{-2\mathrm{i\delta }}\right)(k_x,k_y-k_m).$$

Since the updated complex density $\mathit{\text{FT}}\left({\rho }_n^{\pm }\right){|}_{{\mathcal{P}}^{+}\cup {\mathcal{P}}^0}=\mathit{\text{FT}}\left({\rho }^{\pm }\right){|}_{{\mathcal{P}}^{+}\cup {\mathcal{P}}^0}$, the remainder term $\Theta ({\rho }^{+}-{\rho }_{n-1}^{+},{\mathcal{P}}^0\cup {\mathcal{P}}^{+})=0$. Thus, the discrete ℓ²-norm of the difference between the true and the iteratively updated $T_2^{\ast }$ weighted spin density can be estimated

$$\begin{array}{l}\parallel \left(\mathit{\text{FT}}\left(\right({\rho }^{+}-{\rho }_{n-1}^{+}\left)e^{-2\mathrm{i\delta }}\right)\right){|}_{{\mathcal{P}}^{+}}{\parallel }_{{\ell }_2}^2\\ \le \parallel \mathit{\text{FT}}\left(\right({\rho }^{+}-{\rho }_{n-1}^{+}){|}_{{\mathcal{P}}^{-}}{\parallel }_{{\ell }_2}^2{\sum }_{k_y\in {\mathcal{P}}^{+}}\left[{\sum }_{(k_x,k_m)\in {\mathcal{P}}^{-}}|{\mathit{\text{FT}}\left(e^{-2\mathrm{i\delta }}\right)(k_x,k_y-k_m)|}^2\right]\end{array}$$

(19)

Detailed estimates of ℓ²-norm calculation are presented in the Appendix C Estimation of ℓ²-norm.

Define an estimator for the convergence of the proposed algorithm

$${\mathcal{Z}}_{2\delta ,{\mathcal{P}}^{\pm }}:=\sum _{k_y\in {\mathcal{P}}^{+}}\left[\sum _{(k_x,k_m)\in {\mathcal{P}}^{-}}{\left|\mathit{\text{FT}}\left(e^{-2\mathrm{i\delta }}\right)(k_x,k_y-k_m)\right|}^2\right]$$

(20)

Using the same procedures of (17) and (19), we have

$$\begin{array}{c}\parallel \mathit{\text{FT}}\left(\right({\rho }^{+}-{\rho }_{n-1}^{+}){|}_{{\mathcal{P}}^{-}}{\parallel }_{{\ell }_2}^2=\parallel \left(\mathit{\text{FT}}\left(\right({\rho }^{-}-{\rho }_{n-2}^{-}\left)e^{-2\mathrm{i\delta }}\right)\right){|}_{{\mathcal{P}}^{+}}(k_x,k_y){\parallel }_{{\ell }_2}\\ \le {\mathcal{Z}}_{2\delta ,{\mathcal{P}}^{\pm }}\parallel \mathit{\text{FT}}\left(\right({\rho }^{-}-{\rho }_{n-2}^{-}){|}_{{\mathcal{P}}^{-}}{\parallel }_{{\ell }_2}^2\end{array}$$

(21)

The relations (17), (19) and (21) show that the convergence of the proposed method depends on the estimator ${\mathcal{Z}}_{2\delta ,{\mathcal{P}}^{\pm }}$:

$$\parallel ({\rho }^{-}-{\rho }_n^{-}){\parallel }_{{\ell }_2}\le {\mathcal{Z}}_{2\delta ,{\mathcal{P}}^{\pm }}\parallel ({\rho }^{-}-{\rho }_{n-2}^{-}){\parallel }_{{\ell }_2}$$

(22)

Optimization of B _z using gradient multi-echo data

Since the noise standard deviation $s_{B_z}$ of the measured B _z is inversely proportional to injection current duration T _c and the SNR of MR magnitude image Υ _M [25, 26] as

$$s_{B_z}\left(\mathbf{r}\right)=\frac{1}{\sqrt{2}{T}_c{{\rm Y}}_M\left(\mathbf{r}\right)}$$

(23)

To reduce the noise level of B _z and the imaging time, we applied the proposed method to the gradient multi-echo pulse sequence as a fast MR imaging technique. Subsequently, by using the gradient multi-echo, it is possible to inject the current for a long duration to maximize the off-resonance phase.

When $T_{E_1}$ denotes the first echo time and Δ T _E is the echo spacing, the m-th echo time is $T_{E_m}=T_{E_1}+(m-1)\Delta T_E,m=1,\cdots ,N_E$, where N _E is the echo number. The phase artifact ${\delta }_{T_{E_m}}$ depends on the echo time $T_{E_m}$. Using a priori estimation of the phase artifact map ${\delta }_{T_{E_m}}$ from a reference scan, for a time varying functional MREIT technique, the measured k-space region including ${\mathcal{P}}^0$ and ${\mathcal{P}}^{+}$ can be determined by taking account of the imaging multi-echo times $T_{E_m}$, the echo number N _E and the repetition time T _R .

The recovered multiple $T_2^{\ast }$-weighted complex densities ${\rho }_{T_{E_m}}^{\pm },m=1,\cdots ,N_E,$ using the proposed algorithm in each echo time $T_{E_m}$ can be optimized to generate a representative measured $B_z^m$ data

$$B_z^m\left(\mathbf{r}\right)=\frac{1}{2\gamma T_c}\stackrel{-1}{tan}\left(\frac{{\alpha }^m\left(\mathbf{r}\right)}{{\beta }^m\left(\mathbf{r}\right)}\right)$$

(24)

where α^m and β^m are the imaginary and real parts of ${\rho }_{T_{E_m}}^{+}/{\rho }_{T_{E_m}}^{-}$, respectively. The recovered multiple $B_z^m$ data include pixel-by-pixel different noise level depending on the different imaging time $T_{E_m}$ and the width of injection current. The multiple measured $B_z^m$ data are optimally combined to reduce the noise level of B _z [27]:

$$B_z\left(\mathbf{r}\right)=\sum _{m=1}^{N_E}{\omega }_m\left(\mathbf{r}\right)B_z^m\left(\mathbf{r}\right)$$

(25)

where the point-wise weighting factor ω _m (r) is given as

$${\omega }_m\left(\mathbf{r}\right)=\frac{1/{\left(s_{B_z^m}\right(\mathbf{r}\left)\right)}^2}{\sum _{k=1}^{N_E}1/{\left(s_{B_z^k}\right(\mathbf{r}\left)\right)}^2}$$

(26)

Now, we setup an algorithm to reconstruct the magnetic flux density B _z data using the partially measured k-space data $S_p^{\pm }$ and involving the following steps:

• 1. Take an initial guess ${\rho }_0^{\pm }=F{T}^{-1}\left(S^{\pm }{|}_{{\mathcal{P}}^0\cup {\mathcal{P}}^{+}}\right)$.

• 2. Transform the n-th updated ${\rho }_n^{\pm }{e}^{-2\mathrm{i\delta }}$ to k-space by taking the Fourier transform.

• 3. Update the measured k-space data $S_{n+1}^{\pm }$ by filling the skipped region in $S_p^{\pm }$ with the transformed $\mathit{\text{FT}}\left({\rho }_n^{\pm }{e}^{-2\mathrm{i\delta }}\right)$ data.

• 4. Update ${\rho }_{n+1}^{\pm }$ by taking the two-dimensional inverse Fourier transform.

• 5. Stop if $\frac{\parallel {\rho }_{n+1}^{\pm }-{\rho }_n^{\pm }{\parallel }_{\Omega }}{\parallel {\rho }_n^{\pm }{\parallel }_{\Omega }}\le \epsilon$ where ε>0 is a given tolerance and ∥·∥ _Ω is a standard L ²-norm in Ω. Otherwise, repeat the process.

• 6. Reconstruct $B_z^m$ using (24) for m=1,2,⋯,N _E .

• 7. Determine a weighting factor map ω _m , m=1,2,⋯,N _E using (26) and reconstruct an optimally weighted B _z in (25).

Experimental setup

Numerical simulation setup

To validate the proposed algorithm, we performed numerical simulations with the two-dimensional finite-element model of a object 20×20 cm² with 256×256 rectangular elements and with the origin at its bottom-left, shown in Figure 1. We added the different complex field inhomogeneity artifacts to the simulated spin density image in Figure 1(a). The target magnetic flux density B _z in Figure 1(b) was generated by solving the elliptic equation (6) and by using the Biot-Savart law given by (4). Since the magnetic flux density B _z in Figure 1(b) is continuous and has no abrupt changes, we used |∇B _z | image displayed in Figure 1(c) to enhance image for the magnetic flux density. Figure 1(d) shows a partially measured k-space data corresponding to ρ⁺.

Figure 1

Simulation setup. a) T₂-weighted spin density, b) simulated magnetic flux density B _z image, c) intensity of |∇B _z image, d) magnitude image of the k-space data on the measured region.

The target conductivity distribution σ had different anomalies with different conductivity values and the amount of injection current was 10 mA. Set an applied current density g on the surface as

$$g(x,y)=\left\{\begin{array}{l}10,\text{if}|y-5|\le 2,\text{and}x=0\\ -10,\text{if}|y-5|\le 2,\text{and}x=10\\ 0,\text{otherwise.}\end{array}\right.$$

(27)

Phantom imaging experimental setup

For the practical application of the proposed method as a fast MREIT imaging, we designed a cylindrical phantom filled with the saline solution of conductivity 1 S/m (shown in Figure 2(a)), including a rubber balloon for the visualization of isotropic conductivity excluding other artifacts by any concentration gradient in the phantom. The inside of balloon was filled with the same saline solution and the volume of balloon was controlled by injecting saline solution during the imaging experiment. After positioning the phantom inside the bore of 3T MR scanner (Achieva TX, Philips Medical Systems, Best, The Netherlands) with 8 channel RF coil, we collected k-space data using the gradient multi-echo injection current nonlinear encoding (ICNE) pulse sequence which was originated from FALSH (Figure 2(b)). To obtain the MR magnitude and magnetic flux density (B _z ) images, it extends throughout the duration of injection current until the end of a readout gradient [28]. Since the multi-echo ICNE pulse sequence was synchronized to the injection currents with alternating polarity, it enabled to maximize the width of the injection currents and minimize the noise standard deviation of the measured B _z data. The maximum amplitude of injection current was 5 mA and the total imaging time was 10.24 second to fill the k-space in the interleaved sense. Since the total imaging time was corresponding to the whole k-space scan, the actual imaging time would be reduced to 5.52 and 5.92 seconds for the number of partial region $\mathcal{N}({\mathcal{P}}^0\cup {\mathcal{P}}^{+})=69$ and 74, respectively. The imaging parameters were followings: slice thickness 5 mm, number of imaging slices one, repetition time T _R =40 ms, echo spacing Δ T _E =6 ms, flip angle 40 degree, and multi-echo time $T_{E_m}=6+(m-1)\times 6$ ms for N _E =4. The FOV was 160 ×160 mm² with a matrix size of 128×128. The duration of current injection $T_{c_m}$ was almost same to the multi-echo time $T_{E_m}=6+(m-1)\times 6,m=1,2,3,4$, because the current was continuously injected until the end of the readout gradient.

Figure 2

Experimental setup. a) saline phantom with balloon, b) diagram of the ICNE-multi-echo MR pulse sequence based on a gradient echo.

Figure 3(a) shows the multiply acquired magnitude images $\left|{\rho }_{T_{E_m}}^{+}\right|,m=1,\cdots ,4$, where ${\rho }_{T_{E_m}}^{+}$ was the m-th measured $T_2^{\ast }$ weighted complex spin density, Figure 3(c) and (e) show the measured magnetic flux densities $B_z^m$ and the absolute of $\nabla B_z^m$ images at each echo m=1,⋯,4, respectively. The slope of $B_z^m$ reflecting the width of injected current linearly increased as the multi-echo time $T_{E_m}$ was increasing. Figure 3(b), (d), and (f) are the averaged images corresponding to Figure 3(a), (c), and (e), respectively.

Figure 3

Phantom imaging. a)$T_2^{\ast }$-weighted magnitude images $\left|{\rho }_{T_{E_m}}^{+}\right|,m=1,\cdots ,4$, b) averaged $T_2^{\ast }$ weighted magnitude image $\left|{\bar{\rho }}^{+}\right|:=\frac{1}{4}\sum _{m=1}^4\left|{\rho }_{T_{E_m}}^{+}\right|$, c) measured $B_z^m$ images at each echo m=1,⋯,4, d) averaged B _z image, $B_z:=\frac{1}{4}\sum _{m=1}^4{B}_z^m$, e)$|\nabla B_z^m|$ images at each echo m=1,⋯,4, f) averaged |∇B _z | image, $|\nabla B_z|:=\frac{1}{4}\sum _{m=1}^4|\nabla B_z^m|$.

Results

Simulation results

We compared the reconstructed magnetic flux density B _z using the complete k-space data to the B _z achieved using the partial k-space data. To evaluate the convergence characteristics of the proposed algorithm, we define the relative ℓ²-errors:

$$E\left({\psi }_n\right):=\frac{\parallel \psi -{\psi }_n\parallel }{\parallel \psi \parallel }$$

(28)

where ψ and ψⁿ are the recovered images with the complete k-space data and the n-th updated data using the partially measured k-space data, respectively, and ∥·∥ denotes the ℓ²-norm.

To investigate the estimator ${\mathcal{Z}}_{2\delta ,{\mathcal{P}}^{\pm }}$ for the convergence of the proposed iterative algorithm to fill the sipped k-space region ${\mathcal{P}}^{-}$, we fixed the number of partially measured k-space region as $\mathcal{N}({\mathcal{P}}^0\cup {\mathcal{P}}^{+})=138$, i.e., $\mathcal{N}\left({\mathcal{P}}^0\right)=20$ and $\mathcal{N}\left({\mathcal{P}}^{+}\right)=118$, and changed the frequency range of background field inhomogeneity. We generated several background field inhomogeneities changing the phase frequency range in k-space region by taking ${\mathcal{N}}_{\delta }=5,10,20,30$ where the background field inhomogeneity δ satisfies F T(e^iδ)(k _x ,k _y )=0 for $\left|k_y\right|>{\mathcal{N}}_{\delta }$.

Figure 4(a)-(d) show the background field inhomogeneities used in the reconstruction procedure and Table 1 shows the estimated ${\mathcal{Z}}_{2\delta ,{\mathcal{P}}^{\pm }}$ in each background field inhomogeneity for ${\mathcal{N}}_{\delta }=5,10,20,$ and 30, respectively.

Figure 4

Simulated background field inhomogeneity distributions. a-d) real part images of e^iδ for ${\mathcal{N}}_{\delta }=5,10,20,30$.

Table 1 Calculated estimator${\mathcal{Z}}_{2\delta ,{\mathcal{P}}^{\pm }}$ in (20) for${\mathcal{N}}_{\delta }=5,10,20,30$ and fixed measuredk-space data$\mathcal{N}({\mathcal{P}}^0\cup {\mathcal{P}}^{+})=138$

Figure 5(a) shows the reconstructed |∇B _z | images using the fixed background field inhomogeneity with${\mathcal{N}}_{\delta }=5$. From the top-left to the bottom-right, each image was corresponding to the j-th iterative updated |∇B _z |, j=0,1,⋯,10. The recovered B _z using the partially measured k-space data included a large amount of artifacts (the top-left image in Figure 5(a)). However, since the value of${\mathcal{Z}}_{2\delta ,{\mathcal{P}}^{\pm }}$ was small, the first updated magnetic flux density almost recovered the true B _z .

Figure 5

Reconstructed|∇B _z | images using the fixed background field inhomogeneity with${\mathcal{N}}_{\delta }=5$ and${\mathcal{N}}_{\delta }=30$.a) reconstructed |∇B _z | images using the fixed background field inhomogeneity with${\mathcal{N}}_{\delta }=5$. b) reconstructed |∇B _z | images using the fixed background field inhomogeneity with${\mathcal{N}}_{\delta }=30$. From top-left to bottom-right, each image corresponds to the j-th iterative updated |∇B _z |, j=0,1,⋯,9.

Figure 5(b) shows reconstructed |∇B _z | images using the fixed background field inhomogeneity with${\mathcal{N}}_{\delta }=30$ corresponding to Figure 5(a). Since the value of${\mathcal{Z}}_{2\delta ,{\mathcal{P}}^{\pm }}$ was 0.7738, the convergence ratio of ρ^± was relatively slow comparing to the field inhomogeneity with${\mathcal{N}}_{\delta }=5$. Table 2 shows the relative discrete ℓ²-errors of updated complex spin density for each iteration number and the convergence ratio of ρ^± was depending on the number of${\mathcal{N}}_{\delta }$.

Table 2 Relativeℓ²-errors of updated complex spin density for each iteration number(♯)

Table 3 shows the relative ℓ²-errors of reconstructed ∇B _z for each iteration number depending on the number of${\mathcal{N}}_{\delta }$. The decay rates of the relative ℓ²-error were very fast as the number of${\mathcal{N}}_{\delta }$ was small, but we needed relatively many iterations to approach the required accuracy as the number of${\mathcal{N}}_{\delta }$ was increase, even though the update procedure was rapidly computed by use of the fast Fourier transform.

Table 3 Relative discrete ℓ²-errors of∇B _z for each iteration number(♯)

Phantom experimental results

For the phantom experiment, we changed N _c =5,⋯,10 for the set${\mathcal{P}}^0$ to investigate the convergence behavior with respect to a given background field inhomogeneity. Using the collected k-space data with 8 channel RF coil and the gradient multi-echo by alternating readout gradient, we measured the$T_2^{\ast }$ weighted complex densities${\rho }_m^{\pm n},n=1,\cdots ,N_{\mathit{\text{CH}}},m=1,\cdots ,N_E$, where N _{C H} =8 denotes the coil number and N _E =4 is the echo number. Figure 6 shows the measured background field inhomogeneities by displaying the real part of$e^{2i{\delta }_{\mathit{\text{mn}}}}$ corresponding to the n-th coil and the m-th echo image. According to the increase of echo number, the accumulated background field inhomogeneity also increased.

Figure 6

Measured background field inhomogeneity distributions. Real part of$e^{2i{\delta }_{\mathit{\text{mn}}}},n=1,\cdots ,$N _{C H} , m=1,⋯,N _E , where N _{C H} =8 and N _E =4 denote the coil and echo numbers, respectively.

Figure 7(a) and (c) show the measured$T_2^{\ast }$ weighted magnitude and magnetic flux density B _z images at the time$T_{E_m},m=1,2,3,4$, using partially acquired k-space region${\mathcal{P}}^0\cup {\mathcal{P}}^{+}$ with N _c =5 by a transversally injected current. Although the amount of accumulated phase signal by the injected current increased as the echo time varied from$T_{E_1}$ to$T_{E_4}$, the magnitude image at the 4-th echo was more deteriorated comparing to the 1-st echo case. Figure 7(b) shows an averaged MR magnitude image at each echo time and Figure 7(d) is a weighted B _z image depending on the width of injected current using the phase signal in Figure 7(c). Figure 7(e)-(h) shows the measured magnitude and magnetic flux density B _z images corresponding to Figure 7(a)-(d) using partially acquired k-space region${\mathcal{P}}^0\cup {\mathcal{P}}^{+}$ with N _c =10.

Figure 7

Measured$T_2^{\ast }$ weighted magnitude and magnetic flux densityB _z images using${\mathcal{P}}^0\cup {\mathcal{P}}^{+}$ withN _c =5 and N _c =10.a) and e)$T_2^{\ast }$ weighted magnitude image at each echo time$T_{E_m},m=1,2,3,4$, with N _c =5 and N _c =10, respectively. b) and f) combined$T_2^{\ast }$ weighted magnitude image with N _c =5 and N _c =10, respectively. c) and g) recovered B _z image at each echo time$T_{E_m},m=1,2,3,4$ with N _c =5 and N _c =10, respectively. d) and h) weighted B _z image using multiple B _z image at each echo time with N _c =5 and N _c =10, respectively.

Comparing to the measured images in Figure 7(a)-(d), in contrast to the recovery of low phase frequency information corresponding to${\mathcal{P}}^0$, the increased background field inhomogeneity caused relatively high frequency artifacts.

Figure 8(a)-(d) shows iteratively updated$T_2^{\ast }$ weighted magnitude and magnetic flux density B _z images using${\mathcal{P}}^0\cup {\mathcal{P}}^{+}$ with N _c =5. We fixed the update iteration number as 20 for all experiments. When we fixed N _c =5, the 1-st and 2-nd recovered$T_2^{\ast }$ weighted complex densities in Figure 8(a) and (c) were relatively close to the recovered ones using the complete k-space data. However, as the phase artifact increased, the 3-rd and 4-th recovered$T_2^{\ast }$ weighted complex densities were deficient in reflecting full information of B _z signal. Especially, the 4-th updated magnetic flux density B _z image shows some defective region due to the insufficient recovery of$T_2^{\ast }$ weighted complex density.

Figure 8

Iteratively updated$T_2^{\ast }$ weighted magnitude and magnetic flux densityB _z images using${\mathcal{P}}^0\cup {\mathcal{P}}^{+}$ withN _c =5 andN _c =10.a) and e) recovered$T_2^{\ast }$ weighted magnitude image at each echo time$T_{E_m},m=1,2,3,4$, with N _c =5 and N _c =10, respectively. b) and f) combined$T_2^{\ast }$ weighted magnitude image using the recovered magnitude image at each echo time with N _c =5 and N _c =10, respectively. c) and g) recovered B _z image at each echo time$T_{E_m},m=1,2,3,4$,with N _c =5 and N _c =10, respectively. d) and h) weighted B _z image using multiple B _z image at each echo time with N _c =5 and N _c =10, respectively.

Figure 8(e)-(f) shows iteratively updated$T_2^{\ast }$ weighted magnitude and magnetic flux density B _z images corresponding to Figure 8(a)-(d) using${\mathcal{P}}^0\cup {\mathcal{P}}^{+}$ with N _c =10. When we used N _c =10, the updated$T_2^{\ast }$ weighted complex densities almost recovered the magnetic flux density B _z data comparing to those using the complete k-space data.

Discussion

We used the gradient multi-echo MREIT pulse sequence to reduce the imaging time and to maximize injection current duration. Since the MREIT techniques utilize accumulated phase signal by the injected current, it requires enough repetition time T _R to accumulate the phase signal. In this sense, the gradient multi-echo MREIT pulse sequence seems practical approach for the improvement of B _z quality as well as reducing the imaging time. In this paper, we used a partially acquired k-space data in the phantom experiment by filling the k-space as much as 74 line by line, results in 5.92 second to image the resolution of 128×128. Experimental results show that the proposed interleaved partial Fourier strategy for MREIT has a potential to reduce scan times and maintain the information of B z data comparable to what is obtained with complete k-space data.

The convergence ratio of the iteratively updated phase signal heavily depends on the frequency of the background filed inhomogeneity and the number of half-Fourier over-sampling phase-encodes${\mathcal{P}}^0$. Instead of the gradient multi-echo, if we use the spin multi-echo pulse sequence, the proposed iterative algorithm would rapidly recover T₂-weighted complex spin density due to a small amount of background field inhomogeneity. However, in spite of some advantages of the spin multi-echo MREIT pulse sequence, for a real-time MREIT imaging, MR pulse sequence should be carefully investigate by taking into account of the width of injection current, the scan duration and the low SNR of measured B_z signal.

In this paper, we assumed a priori background field inhomogeneity which is typically used in the sensitivity encoding (SENSE) as a fast MRI measurement technique. Since the MREIT techniques typically used interleaved acquisition by injecting alternative currents, it may be possible to extract background field inhomogeneity information under a low frequency range assumption and by cancelation of B_z information:

$$\begin{array}{lll}{\rho }^{+}(x,y){\rho }^{-}(x,y)& =& {\rho }^2(x,y)e^{2\mathrm{i\delta }(x,y)}\\ \frac{{\rho }^{+}(x,y)}{{\rho }^{-}(x,y)}& =& e^{2\mathrm{i\gamma }{B}_z(x,y)T_c}\end{array}$$

(29)

Several studies reported for the feasibility of MREIT to detect neural activities in the brain, directly [29, 30]. Functional MREIT technique is suggested to image brain activity via conductivity change related to neural activity through the fast MREIT pulse sequence. Our future study will focus on applying the proposed method to produce functional conductivity images of animal and/or human brain to pursue rapidly changing conductivity associated with neural activities.

Conclusion

In MREIT, the inherent challenges are to reduce the scan time and maintain current injection duration to make it feasible for the clinical applications. We developed an iterative method to optimize the measured magnetic flux density B_z using the multi-echo interleaved partial Fourier acquisitions for fast imaging in MREIT. The proposed method used a fast gradient multi-echo MR pulse sequence to reduce the scan time and to maximize the phase signal by injection current. Under the assumption of a priori background field inhomogeneity map, we rigorously investigated the convergence ratio of the proposed method using the discrete ℓ²-norm, which was closely related with the number of measured phase encode set and the frequency range of the background field inhomogeneity. To evaluate the proposed method, a specially designed conductivity phantom was used to provide a homogeneous magnitude, but it yielded distinguishable B_z signal inside and outside the anomaly. For the phantom experiment, total imaging time was 10.24 seconds to fill the complete k-space region in the interleaved sense and it was less than 6 seconds to fill the partial k-space region to implement the proposed method. The proposed interleaved partial Fourier strategy for the fast MREIT has a potential to reduce scan times and maintain the information of B_z data comparable to what is obtained with the complete k-space data.

Appendix

A Proof of Observation 1

For the phase-encode$(k_x,k_y)\in {\mathcal{P}}^{+}$, the Fourier transform of${\rho }_p^{-}{e}^{-2\mathrm{i\delta }}$ can be separated as

$$\begin{array}{ll}\hfill \mathit{\text{FT}}\left({\rho }_p^{-}{e}^{-2\mathrm{i\delta }}\right)(k_x,k_y)& =\left(\mathit{\text{FT}}\right({\rho }_p^{-})\ast \mathit{\text{FT}}(e^{-2\mathrm{i\delta }}\left)\right)(k_x,k_y)\\ ={\sum }_{(k_x,k_m)\in {\mathcal{P}}^{+}\cup {\mathcal{P}}^0\cup {\mathcal{P}}^{-}}\mathit{\text{FT}}\left({\rho }_p^{-}\right)(k_x,k_m)\mathit{\text{FT}}\left(e^{-2\mathrm{i\delta }}\right)(k_x,k_y-k_m)\\ ={\sum }_{(k_x,k_m)\in {\mathcal{P}}^{+}\cup {\mathcal{P}}^0}\mathit{\text{FT}}\left({\rho }_p^{-}\right)(k_x,k_m)\mathit{\text{FT}}\left(e^{-2\mathrm{i\delta }}\right)(k_x,k_y-k_m)\\ +{\sum }_{(k_x,k_m)\in {\mathcal{P}}^{-}}\mathit{\text{FT}}\left({\rho }_p^{-}\right)(k_x,k_m)\mathit{\text{FT}}\left(e^{-2\mathrm{i\delta }}\right)(k_x,k_y-k_m)\end{array}$$

(30)

where ∗ denotes the convolution with respect to k _y . Since the updated$S_p^{-}$ data conserve the measured data in${\mathcal{P}}^{+}\cup {\mathcal{P}}^0$,$\mathit{\text{FT}}\left({\rho }_p^{-}\right)(k_x,k_y)=\mathit{\text{FT}}\left({\rho }^{-}\right)(k_x,k_y)$ for$(k_x,k_y)\in {\mathcal{P}}^{+}\cup {\mathcal{P}}^0$. Thus, we have

$$\begin{array}{ll}\hfill \mathit{\text{FT}}\left({\rho }_p^{-}{e}^{-2\mathrm{i\delta }}\right)(k_x,k_y)=& {\sum }_{(k_x,k_m)\in {\mathcal{P}}^{+}\cup {\mathcal{P}}^0}\mathit{\text{FT}}\left({\rho }^{-}\right)(k_x,k_m)\mathit{\text{FT}}\left(e^{-2\mathrm{i\delta }}\right)(k_x,k_y-k_m)\\ +{\sum }_{(k_x,k_m)\in {\mathcal{P}}^{-}}\mathit{\text{FT}}\left({\rho }_p^{-}\right)(k_x,k_m)\mathit{\text{FT}}\left(e^{-2\mathrm{i\delta }}\right)(k_x,k_y-k_m)\\ ={\sum }_{(k_x,k_m)\in {\mathcal{P}}^{+}\cup {\mathcal{P}}^0}\mathit{\text{FT}}\left({\rho }^{-}\right)(k_x,k_m)\mathit{\text{FT}}\left(e^{-2\mathrm{i\delta }}\right)(k_x,k_y-k_m)\end{array}$$

(31)

Since the central phase-encode set${\mathcal{P}}^0$ includes all phase frequencies of the systematic phase artifact e^−iδ, the range of the phase frequency k _y −k _m for$(k_x,k_y)\in {\mathcal{P}}^{+}$ and$(k_x,k_m)\in {\mathcal{P}}^{-}$ is over 2N_δ. This means that F T(e^{−2i δ})(k _x ,k _y −k _m )=0. Thus, we have$\mathit{\text{FT}}\left({\rho }_p^{-}{e}^{-2\mathrm{i\delta }}\right)(k_x,k_y)=\mathit{\text{FT}}\left({\rho }^{-}{e}^{-2\mathrm{i\delta }}\right)(k_x,k_y)$ for$(k_x,k_y)\in {\mathcal{P}}^{+}$. The case for${\rho }_p^{+}{e}^{-2\mathrm{i\delta }}$ is similar.

B Proof of Observation 2

Since$\mathit{\text{FT}}\left({\rho }_p^{\mp }{e}^{-2\mathrm{i\delta }}\right)(k_x,k_y)=\mathit{\text{FT}}\left({\rho }^{\mp }{e}^{-2\mathrm{i\delta }}\right)(k_x,k_y)$ for$(k_x,k_y)\in {\mathcal{P}}^{+}$ due to the observation 1, we have

$$\mathit{\text{FT}}\left({\rho }^{+}{e}^{-2\mathrm{i\delta }}\right)(k_x,k_y)=\mathit{\text{FT}}\left(\rho e^{-\mathrm{i\delta }}{e}^{\mathrm{i\gamma }{B}_z{T}_c}\right)(k_x,k_y)\text{for}(k_x,k_y)\in {\mathcal{P}}^{+}$$

(32)

From the relation (32), by taking the complex conjugate, we recover the skipped k-space region${\mathcal{P}}^{-}$

$$\begin{array}{l}\mathit{\text{FT}}\left({\rho }^{-}\right)(k_x,k_y)=\mathit{\text{FT}}\left(\rho e^{\mathrm{i\delta }}{e}^{-\mathrm{i\gamma }{B}_z{T}_c}\right)(k_x,k_y)\\ =\overline{\mathit{\text{FT}}\left(\rho e^{-\mathrm{i\delta }}{e}^{\mathrm{i\gamma }{B}_z{T}_c}\right)}(-k_x,-k_y)\\ =\overline{\mathit{\text{FT}}\left({\rho }_p^{+}{e}^{-2\mathrm{i\delta }}\right)}(-k_x,-k_y)\end{array}$$

(33)

C Estimation of ℓ²-norm

The discrete ℓ²-norm of the difference between the true and the iteratively updated$T_2^{\ast }$ weighted spin density can be estimated as following:

$$\begin{array}{l}\parallel \left(\mathit{\text{FT}}\left(\right({\rho }^{+}-{\rho }_{n-1}^{+}\left)e^{-2\mathrm{i\delta }}\right)\right){|}_{{\mathcal{P}}^{+}}{\parallel }_{{\ell }_2}^2\\ ={\sum }_{k_y\in {\mathcal{P}}^{+}}\left|\mathit{\text{FT}}\left(({\rho }^{+}-{\rho }_{n-1}^{+})e^{-2\mathrm{i\delta }}\right)\right(k_x,k_y){|}^2\\ ={\sum }_{k_y\in {\mathcal{P}}^{+}}{\left|{\sum }_{(k_x,k_m)\in {\mathcal{P}}^{-}}\mathit{\text{FT}}({\rho }^{+}-{\rho }_{n-1}^{+})(k_x,k_m)\mathit{\text{FT}}\left(e^{-2\mathrm{i\delta }}\right)(k_x,k_y-k_m)\right|}^2\\ \le {\sum }_{k_y\in {\mathcal{P}}^{+}}{\left[{\sum }_{(k_x,k_m)\in {\mathcal{P}}^{-}}\left|\mathit{\text{FT}}\right({\rho }^{+}-{\rho }_{n-1}^{+}\left)\right(k_x,k_m\left)\mathit{\text{FT}}\right(e^{-2\mathrm{i\delta }}\left)\right(k_x,k_y-k_m\left)\right|\right]}^2\\ \le {\sum }_{k_y\in {\mathcal{P}}^{+}}\left[\parallel \mathit{\text{FT}}\left(\right({\rho }^{+}-{\rho }_{n-1}^{+}\left){|}_{{\mathcal{P}}^{-}}{\parallel }_{{\ell }_2}^2{\sum }_{(k_x,k_m)\in {\mathcal{P}}^{-}}\right|\mathit{\text{FT}}\left(e^{-2\mathrm{i\delta }}\right)(k_x,k_y-k_m){|}^2\right]\\ \le \parallel \mathit{\text{FT}}\left(\right({\rho }^{+}-{\rho }_{n-1}^{+}){|}_{{\mathcal{P}}^{-}}{\parallel }_{{\ell }_2}^2{\sum }_{k_y\in {\mathcal{P}}^{+}}\left[{\sum }_{(k_x,k_m)\in {\mathcal{P}}^{-}}\left|\mathit{\text{FT}}\right(e^{-2\mathrm{i\delta }}\left)\right(k_x,k_y-k_m){|}^2\right]\end{array}$$

(34)