3D Reconstruction of the Magnetic Vector Potential using Model Based Iterative Reconstruction

TL;DR

This paper introduces a model-based iterative reconstruction method for 3D magnetic vector potential imaging in TEM, improving accuracy over traditional filtered back projection methods, especially with incomplete data.

Contribution

The paper presents a novel MBIR algorithm that combines a forward TEM model with a prior to enhance 3D magnetic vector potential reconstructions from TEM data.

Findings

MBIR outperforms VFET in reconstruction quality

Quantitative improvements demonstrated on simulated data

Experimental data confirms enhanced reconstruction accuracy

Abstract

Lorentz Transmission Electron Microscopy (TEM) observations of magnetic nanoparticles contain information on the magnetic and electrostatic potentials. Vector Field Electron Tomography (VFET) can be used to reconstruct electromagnetic potentials of the nanoparticles from their corresponding LTEM images. The VFET approach is based on the conventional filtered back projection approach to tomographic reconstructions and the availability of an incomplete set of measurements due to experimental limitations means that the reconstructed vector fields exhibit significant artifacts. In this paper, we outline a model-based iterative reconstruction (MBIR) algorithm to reconstruct the magnetic vector potential of magnetic nanoparticles. We combine a forward model for image formation in TEM experiments with a prior model to formulate the tomographic problem as a maximum a-posteriori probability estimation problem (MAP). The MAP cost function is minimized iteratively to determine the vector potential. A comparative reconstruction study of simulated as well as experimental data sets show that the MBIR approach yields quantifiably better reconstructions than the VFET approach.

Figures

Tables

Table 1. Table 1: Root Mean Square Error (RMSE) analysis between the Shepp-Logan phantom and the reconstructions for different tomographic approaches.

Tomographic Method	RMSE
Back-Projection (BP)	$0.4635$
Filtered Back-Projection (FBP)	$0.1536$
Simultaneous Iterative Reconstruction Technique (SIRT)	$0.0793$
Model-Based Iterative Reconstruction (MBIR)	$0.0213$

Equations

$$\phi ({\mathbf{r}}_{\perp })={\phi }_e({\mathbf{r}}_{\perp })+{\phi }_m({\mathbf{r}}_{\perp })=\frac{\pi }{\lambda E_t}\underset{L}{\int }V({\mathbf{r}}_{\perp },z)dz-\frac{e}{\hslash }\underset{L}{\int }\mathbf{A}({\mathbf{r}}_{\perp },z)\cdot d\mathbf{r},$$

$${\tilde{\phi }}_{m,x}(k_u,k_v)=-\sin \theta {\tilde{A}}_y(k_u,k_v\cos \theta ,k_v\sin \theta )+\cos \theta {\tilde{A}}_z(k_u,k_v\cos \theta ,k_v\sin \theta );$$

$${\tilde{\phi }}_{m,y}(k_u,k_v)=-\sin \theta {\tilde{A}}_x(k_u\cos \theta ,k_v,k_u\sin \theta )+\cos \theta {\tilde{A}}_z(k_u\cos \theta ,k_v,k_u\sin \theta ).$$

$$\mathbf{k}.\tilde{\mathbf{A}}=k_x{\tilde{A}}_x+k_y{\tilde{A}}_y+k_z{\tilde{A}}_z=0.$$

$${\tilde{A}}_x$$

$${\tilde{A}}_y$$

$${\tilde{A}}_z$$

$$A_{\ast }(x,y,z)=\underset{0}{\stackrel{\pi }{\int }}\left[{\phi }_{m,x,\ast }(x,y\cos \theta +z\sin \theta )+{\phi }_{m,y,\ast }(x\cos \theta +z\sin \theta ,y)\right]d\theta ,$$

$${\phi }_{m,x,\ast }$$

$${\phi }_{m,y,\ast }$$

$$\mathbf{A}(\mathbf{k})=-\frac{\mathrm{i}{B}_0}{k^2}D(\mathbf{k})(\hat{\mathbf{m}}\times \mathbf{k}),$$

$$\text{NRMSE}=\frac{1}{{\varphi }_{\text{max}}-{\varphi }_{min}}\sqrt{\frac{1}{N}\underset{i=1}{\sum ^N}(\widehat{{\varphi }_i}-{\varphi }_i{)}^2}.$$

$${\hat{x}}_{\text{MAP}}=\underset{x}{\mathrm{argmin}}\left\{-\log \mathcal{P}(y|x)-\log \mathcal{P}(x)\right\}.$$

$$\log \mathcal{P}(y|x)\approx -\frac{1}{2}(y-Hx{)}^{T}W(y-Hx)+f(y),$$

$$\log \mathcal{P}(x)=-\sum _{\{i,j\}\in \mathcal{C}}{b}_{ij}\rho (x_i-x_j),$$

$$\rho (\Delta )=\frac{{\left|\Delta \right|}^p}{p{\sigma }_x^p}\left(\frac{{\left|\frac{\Delta }{T{\sigma }_x}\right|}^{q-p}}{1+{\left|\frac{\Delta }{T{\sigma }_x}\right|}^{q-p}}\right),$$

$${\hat{x}}_{\text{MAP}}=\underset{x}{\mathrm{argmin}}\left\{\frac{1}{2}{\Vert y-Hx\Vert }^2+\sum _{\{i,j\}\in \mathcal{C}}{b}_{ij}\rho (x_i-x_j)\right\},$$

$$c(x)=\frac{1}{2}{\Vert y-Hx\Vert }^2+\sum _{\{i,j\}\in \mathcal{C}}{b}_{ij}\rho (x_i-x_j).$$

$$\rho (\Delta ;{\Delta }^{\prime })=\frac{a_2}{2}{\Delta }^2,$$

$$a_2=\frac{{\rho }^{\prime }({\Delta }^{\prime })}{{\Delta }^{\prime }}.$$

$$\rho(\Delta;\Delta^{\prime})=\begin{cases}\begin{array}[]{c}\frac{\rho^{\prime}(\Delta^{\prime})}{2\Delta^{\prime}}\Delta^{2}\\ \frac{\rho^{\prime\prime}(0)}{2}\Delta^{2}\end{array}&\begin{array}[]{c}\text{if }\Delta\neq 0\\ \text{if }\Delta=0\end{array}\end{cases},$$

$$c(x;x^{\prime })=\frac{1}{2}{\Vert y-Hx\Vert }^2+\sum _{\{i,j\}\in \mathcal{C}}{\tilde{b}}_{ij}(x_i-x_j{)}^2,$$

$${\tilde{b}}_{ij}$$

$$x_i\leftarrow \frac{H^{T}y+2\sum _{\{i,j\}ϵ\mathcal{C}}{\tilde{b}}_{ij}{x}_j}{1+2\sum _{\{i,j\}ϵC}{\tilde{b}}_{ij}}.$$

$$\text{RMSE}=\sqrt{\frac{1}{N}\underset{i=1}{\sum ^N}(\widehat{{\varphi }_i}-{\varphi }_i{)}^2}.$$

$$x^{(\text{new})}\leftarrow x^{(\text{old})}+\lambda H^{T}W(y-H{x}^{(\text{old})})$$

$$\hat{\mathbf{x}}=\underset{\mathbf{x}}{\mathrm{argmin}}\left\{\frac{1}{2}{\Vert \mathbf{D}{\phi }_m-\mathbf{D}(\mathbf{Hx})\Vert }^2+\sum _{\{i,j\}\in \mathcal{C}}{\tilde{b}}_{ij}\mathbf{D}(\mathbf{x}(i)-\mathbf{x}(j){)}^2\right\}$$

$${\hat{x}}_x$$

$${\hat{x}}_y$$

$${\hat{x}}_z$$

$$c(\mathbf{x},x_x;x_x^{\prime })=\frac{1}{2}{\Vert D_x{\phi }_m-D_x(\mathbf{Hx})\Vert }^2+\sum _{\{i,j\}\in \mathcal{C}}{\tilde{b}}_{ij}{D}_x(\mathbf{x}(i)-\mathbf{x}(j){)}^2.$$

$$x_x(i)\leftarrow \frac{\left(H_x^T{\phi }_{m,x,x}+H_y^T{\phi }_{m,y,x}\right)+2\sum _{\{i,j\}\in \mathcal{C}}{\tilde{b}}_{ij}{x}_x(j)}{1+2\sum _{\{i,j\}\in \mathcal{C}}{\tilde{b}}_{ij}}.$$

$$x_y(i)$$

Full text

\DeclareCaptionType

algorithm

1 Introduction

In a Lorentz Transmission Electron Microscopy (LTEM) experiment, an electron propagating through a thin specimen experiences a Lorentz Force $\mathbf{F}_{L}=-e(\mathbf{E}+\mathbf{v}\times\mathbf{B})$ [1] due to the sample’s electrostatic field, $\mathbf{E}$, and magnetic field, $\mathbf{B}$; $-e$ is the electron’s charge and $\mathbf{v}$ its velocity. This (classical) force generates a deflection of the electron trajectory, which can be used to explain the Fresnel and Foucault observation modes [2]. A more robust explanation of the nature of the electron-specimen interaction involves quantum mechanics, in which the electron is described by a wave function $\psi(\mathbf{r}_{\perp})=a(\mathbf{r}_{\perp})e^{i\varphi(\mathbf{r_{\perp}})}$ [2]. Elastic scattering in the sample produces variations of the amplitude $a(\mathbf{r}_{\perp})$, whereas the electromagnetic potentials affect the phase $\varphi(\mathbf{r}_{\perp})$ of the wave; $\mathbf{r}_{\perp}$ is a vector normal to the propagation direction. Aharonov and Bohm [3] showed, in 1959, that the phase of the exit wave function encodes information on the sample’s electrostatic potential, $V(\mathbf{r}_{\perp},z)$, and magnetic vector potential, $\mathbf{A}(\mathbf{r}_{\perp},z)$, as follows:

[TABLE]

where $\hslash$ is the reduced Planck’s constant, $E_{t}$ is the total beam energy, and the integrals are carried out along the beam direction, $L$. The total phase shift, $\varphi$, consists of an electrostatic contribution, $\varphi_{e}$, and a magnetic contribution, $\varphi_{m}$. The phases are not directly observable, but their effect on the image contrast can be determined by considering the point spread function, $\mathcal{T}_{L}(\mathbf{r}_{\perp})$, of the Lorentz lens. The image intensity is then given by the modulus-squared of the convolution product $\psi(\mathbf{r}_{\perp})\varotimes\mathcal{T}_{L}(\mathbf{r_{\perp}})$ [2]. Hence, characterization of the electromagnetic potentials begins with phase shift retrieval from the image intensities, using either electron holography [4] or the transport-of-intensity equation (TIE) formalism [5], which is based on a through-focus series of Fresnel images. We use the linearity of the TIE [6] and time reversal symmetry to retrieve the individual phases $\varphi_{e}$ and $\varphi_{m}$.

Characterization of the electromagnetic fields is then achieved by performing scalar field and vector field tomographic reconstructions to determine $V(\mathbf{r})$ and $\mathbf{A}(\mathbf{r})$, respectively. We refer to the use of vector field tomography to reconstruct the electromagnetic potentials as Vector Field Electron Tomography (VFET) [7]. A schematic of the operations performed to complete an electromagnetic characterization task is shown in Fig. 1.

In this contribution, we primarily focus on vector field tomography to reconstruct $\mathbf{A}(\mathbf{r})$. In recent years, the VFET approach employed the filtered back projection (FBP) approach to perform the reconstructions [8, 9]. Although FBP yields a good estimate with a complete set of measurements, the typical missing wedge of TEM data significantly diminishes the quality of the reconstructions [10]. In addition, typical tilt series are obtained using an angular step size of $2^{\circ}$–$5^{\circ}$ to minimize the necessary pre-processing steps (image alignments) and reduce beam damage. These limitations, collectively, yield a reconstruction result that can exhibit substantial artifacts. To alleviate these problems, we resort to a more robust and statistically based tomographic reconstruction framework known as Model-Based Iterative Reconstruction (MBIR) to determine $\mathbf{A}(\mathbf{r})$. This approach has had considerable success in improving reconstruction quality in scalar tomography [11, 12].

In section 2, we briefly outline conventional VFET and show how we can reconstruct all three component of $\mathbf{A}(\mathbf{r})$ from just two tilt series; we also perform an error analysis of the quality of VFET reconstructions in the presence of a missing wedge. Next, in section 3, we provide an overview of the MBIR framework, and in section 4, we incorporate MBIR into the reconstruction of $\mathbf{A}(\mathbf{r})$ and compare the results with those from conventional VFET reconstructions.

2 Vector Field Electron Tomography

Vector field tomography is relatively new; it was not until 1988 when Norton [13], for the first time, outlined a mathematical model to determine the 2D fluid field from acoustic time travel measurements. Subsequent years saw extensions of 2D vector tomography to 3D cases; in particular, Juhlin [14] resolved the solenoidal part of a divergence free flow field using ultrasound Doppler measurements in 1992. In 2005, Lade et al. [8] presented the VFET model to reconstruct 3D vector fields from longitudinal and transverse measurements. In 2008, Phatak et al. [15] used the VFET approach to reconstruct the magnetic vector potential and induction of magnetic nanoparticles. Since this VFET model is still relevant for our new MBIR method, we devote this section to a brief review of the VFET framework.

Since tomographic reconstructions require a forward model to project the object being reconstructed, we begin by considering the computation of the magnetic phase shift. The relation $\varphi_{m}(\mathbf{r}_{\perp})=-\frac{e}{\hslash}\underset{L}{\int}\mathbf{A}(\mathbf{r}_{\perp},z)\cdot\mathrm{d}\mathbf{r}$ describes the magnetic phase shift obtained at $0^{\circ}$ tilt. To obtain the phase shift for a tilted sample we consider a tilt series around the $x$ axis (counterclockwise); the new coordinate vectors, $\mathbf{t},$ can be expressed in terms of the original ones, $\mathbf{r},$ by $\mathbf{r}=R_{\theta,x}\mathbf{t}$ where $\mathbf{r}=[x\ y\ z]$, $\mathbf{t}=[u\ v\ w]$, and $R_{\theta,x}$ is the counter-clockwise rotation matrix (Fig. 2(a)). From Fig. 2(b), the vectorial line element, $\mathrm{d}\mathbf{l}$, of the projection line $L(v,\theta)$, can be written as $\mathrm{d}\mathbf{l}=[\hat{y}\sin(\theta)-\hat{z}\cos(\theta)]\mathrm{d}l$.

Writing $\mathbf{A}(\mathbf{r})=A_{x}(x,y,z)\hat{x}+A_{y}(x,y,z)\hat{y}+A_{z}(x,y,z)\hat{z}$, a generic projection equation for the $x$ tilt series in Fourier space, $\tilde{\varphi}_{m,x}$, can be obtained as:

[TABLE]

a similar analysis for the $y$ tilt series produces

[TABLE]

Eqs. 2 and 3 represent the Fourier slice theorem for 3D vector fields for $x$ tilt series and $y$ tilt series respectively.

The formulation of the reconstruction procedure by means of the VFET approach begins by imposing a gauge constraint on the magnetic vector potential, i.e., $\nabla\cdot\mathbf{A}=0$. This constraint is written in Fourier space as:

[TABLE]

Combining eqs. 2, 3, and 4 and solving for the components of $\mathbf{A}(\mathbf{k})$ one obtains:

[TABLE]

In principle, one could interpolate the projection measurements from the polar grid to a Cartesian grid and perform an inverse Fourier transform (FT) to solve for $\mathbf{A}(\mathbf{r})$. However, this is a numerically unstable approach and the filtered back projection formulae, as outlined in ref. [16], facilitate a more stable reconstruction. In the present case, the FBP formulae become (the asterisk subscript represents $x$, $y$, or $z$):

[TABLE]

where

[TABLE]

In eqs. 15 and 22, $|k_{u}|$ and $|k_{v}|$ represent the filtering operation and $\mathcal{F}_{2}^{-1}$ the 2D inverse FT. The real-space vector quantities $\varphi_{m,x,*}$ and $\varphi_{m,y,*}$, are then used to determine the three components of $\mathbf{A}(\mathbf{r})$ through an inverse Hough transform. The gauge condition on $\mathbf{A}$ will be incorporated into the MBIR approach as well. However, we will replace the FBP formula of eq. 8 with a Bayesian inference-based objective function which will be described in more detail in the next section.

Before we outline the MBIR framework, we illustrate the numerical implementation of the VFET approach using uniformly and non-uniformly magnetized simulated magnetic nanoparticles (MNPs); in particular, we use a uniformly magnetized prismatic particle and a circular disk with a vortex state. The prismatic nanoparticle (NP) has dimensions of $\left[50\times 50\times 30\right]$ nm and magnetization direction, $\mathbf{m}=\left[\cos\frac{\pi}{6},\sin\frac{\pi}{6},0\right]$. The circular disk has a diameter of $60$ nm, a height of $30$ nm, and exhibits a counterclock-wise magnetization state with a sharp vortex [17, 18]. Both particles have a saturation magnetization of $B_{0}=1$ T. The two geometries and their magnetization states are depicted in Fig. 3.

Starting from the magnetization states, the Fourier transform of the magnetic vector potential can be expressed in terms of the shape function formalism of Ref. [17] as:

[TABLE]

where $D(\mathbf{k})$ is the shape amplitude. This expression allows for the computation of the magnetic phase shifts for two orthogonal tilt series around the $x$ and $y$ axes. In one series we use projections over the full tilt range $[-90^{\circ},+90^{\circ}]$ with a $2^{\circ}$ step size while in a second series, with a missing wedge, we use the sub-range of $[-70^{\circ},+70^{\circ}]$. After application of the VFET approach described above we compare the reconstructed potentials with their corresponding known values (i.e., ground truths), in terms of a normalized root mean square error (NRMSE) between the reconstructed result, $\hat{\phi}$, and the ground truth, $\phi$, as ($\phi$ represents one of the components of the vector potential):

[TABLE]

Fig. 4 compares sections of reconstructed vector potentials with their ground truths, and lists the corresponding NRMSE values. Blurring at edges and ring artifacts are clearly more prominent in the reconstructions from the missing wedge data set, which has a higher NRMSE value than the reconstruction using the full tilt series. Fig. 5 shows the NRMSE values corresponding to the $y$ planes in the range $[-35,35]$. These plots indicate that the reconstruction deduced from the missing wedge data set exhibits a higher level of deviation from the ground truth than its full range counterpart by about $5$–$25$%. In the following section, we will replace the filter based analytical relation of eq. 8 with a Bayesian statistics-based numerical relation.

3 The MBIR framework

We introduce the MBIR approach by first formulating its framework for a generic tomography problem. Let $x\in\mathbb{R}^{N}$ be the discrete vector of an unknown image and $y\in\mathbb{R}^{M}$ be the discrete vector of projection measurements of $x$ at a variety of angles. Then the MBIR approach performs a reconstruction by maximizing the joint probability distribution resulting from the likelihood function, $\mathcal{P}(y|x)$, and the prior distribution, $\mathcal{P}(x)$. This methodology is more commonly known as the maximum-a-posteriori (MAP) estimation technique [19]. Mathematically, the MAP estimate is expressed as:

[TABLE]

We use the Poisson probability mass function to model the likelihood term. We use a second order Taylor series expansion of the logarithm of the Poisson distribution [20] to approximate the first term of eq. 25 as:

[TABLE]

where $H$ is the forward projection matrix, $W$ is a diagonal noise weighting matrix and $f(y)$ is a term independent of the optimization variable $x$. Additionally, $H$ is an orthogonal matrix such that $H^{T}=H^{-1}$; hence, $H^{T}$ denotes the back projection operator.

The prior term in eq. 25 is modeled using a Markov Random Field (MRF) in:

[TABLE]

where $\mathcal{C}$ denotes the set of neighboring pixels, $b_{ij}$ is a non-causal symmetric weighing filter that is normalized such that $\sum b_{ij}=1$, and $\rho(\cdot)$ is the potential function. The choice of the MRF as the prior model is due to its proven usefulness in image processing as well as in tomographic reconstructions [21, 22]. In the present paper, we resort to a class of MRFs called the $q$-Generalized Gaussian Markov Random Field ($q$-GGMRF) to model the prior. Accordingly, the potential function in eq. 27 is given by:

[TABLE]

where $\Delta=x_{i}-x_{j}$, and $p,\ q,\ \sigma_{x}\text{ and }T$ are the $q$-GGMRF parameters. Typically, $1\leq p\leq q\leq 2$ is used to ensure strict convexity of the potential function and, subsequently, of the MAP optimization [23]; in this paper we use $q=2$. When $p$ is set to $2$, the potential function is quadratic and the prior model facilitates a reconstruction with smooth edges. On the other hand, when $p$ is close to $1$, the prior model performs sharp edge preserving reconstructions. Similarly, $\sigma_{x}$ is the variance of the prior distribution and its value is set to achieve a balance between noise and resolution. Finally, the constant $T$ determines the approximate threshold of transition between low and high contrast regions.

Substituting the log-likelihood expression of eq. 26 and the prior model of eq. 27 into eq. 25, and considering $W$ as the identity matrix, one obtains the MAP reconstruction to be the solution of the following optimization problem:

[TABLE]

where the cost (objective) function being minimized is given by

[TABLE]

Before proceeding with the minimization of the cost function, it should be noted that differentiating the cost function with respect to $x_{i}$ and setting the derivative to zero does not yield a closed form solution. Accordingly, one needs to make use of numerical methods, e.g., the Newton-Raphson method or the line search method, to determine the MAP estimate [20, 24]. However, these numerical methods can be computationally expensive to implement. Thus, we introduce the substitute or surrogate function to derive a closed form solution of the MAP estimate and subsequently facilitate a faster implementation [25]. A visual depiction of the intuition behind the substitute function is shown in Fig. 6. In accordance to the figure, symmetric bond method is employed to construct a substitute function, $\rho(\Delta;\Delta^{\prime})$, as an upper bound to the original function, $\rho(\Delta)$, such that minimizing the substitute function results in a lower value of the original function [26]. Accordingly, the potential function in eq. 27 is upper bounded by a symmetric and quadratic function of $\Delta$. Then its substitute form can be expressed as:

[TABLE]

where $a_{2}$ is a function of $\Delta^{\prime}$. Here, $a_{2}$ is determined by matching the gradient of $\rho(\Delta)$ and $\rho(\Delta;\Delta^{\prime})$ at $\Delta=\Delta^{\prime}$ to yield:

[TABLE]

Substituting the value of $a_{2}$ in eq. 31 results in the following symmetric bound substitute function:

[TABLE]

where the limiting value of $a_{2}=\rho^{\prime\prime}(0)$ is used when $\Delta=0$. Subsequently, we obtain the new MAP cost function as:

[TABLE]

where

[TABLE]

The surrogate cost function in eq. 34 is quadratic and so its minimum can be expressed in closed form. We use the Iterative Coordinate Descent (ICD) method to numerically minimize the surrogate cost function [25]. ICD works by sequentially minimizing the cost function of eq. 34 with respect to each pixel. The convergence of the ICD technique to solve the optimization problem can be inferred from the fact that the surrogate cost function in eq. 34 is strictly convex so that any locally optimal point is also globally optimal. Thus, we proceed to differentiate the surrogate cost function and solve for the minimum $x_{i}$ to yield each pixel update as:

[TABLE]

{algorithm}

[t]

1:initialize $u\leftarrow\text{FBP}(y)$

2:$v\leftarrow H^{T}y$

3:while not converged do

4:     for $i=1$ to $N$ do

5:         \text{for }\{i,j\}\in\mathcal{C}\text{ of }u_{i}\text{ determine }\tilde{b}_{ij}\text{ using eq.\leavevmode\nobreak\ \ref{eq:tilde_bij_val}}

6:         $u_{i}\leftarrow(v_{i}+2\underset{\{i,j\}\in\mathcal{C}}{\sum}\tilde{b}_{ij}u_{j})/(1+2\underset{\{i,j\}\in C}{\sum}\tilde{b}_{ij})$

7:     end for

8:     $e\leftarrow y-Hu$

9:     $u\leftarrow u+H^{T}e$

10:end while

Fourier Transform-based MBIR method for scalar tomography

A pseudocode to minimize the surrogate cost function for a pixel wise update using the ICD technique is listed in Algorithm 3. Note that, in contrast to the widely observed MBIR practice of updating the error sinogram, $e$, after each pixel update, we resort, instead, to updating $e$ after all $N$ pixels have been modified. This is primarily due to the fact that we will be using a Fourier-based forward model instead of a Radon based forward model [12] to minimize the cost function. The Fourier-based forward model is chosen to make the framework compatible with imposing a Coulomb gauge constraint when we extend it for the reconstruction of $\mathbf{A}(\mathbf{r})$.

We conclude this section with the implementation of the MBIR methodology to a 2D scalar tomography problem. We make use of the standard Shepp-Logan head phantom [27] to illustrate our results. The phantom or ground truth is depicted in Fig. 7(a); it is $[256\times 256]$ pixels in size and has values ranging from [math] to $1$, which we map on a blue-red color scale. In addition to the result from the MBIR approach (Fig. 7(e)), we also show results from other reconstruction techniques, such as back projection (BP) in Fig. 7(b) and FBP in Fig. 7(c). We used a band-limited Ram-Lak ramp filter as our filtering function [28, 29]. The result from an implementation of a conventional iterative tomography method known as the Simultaneous Iterative Reconstruction Technique (SIRT) is shown in Fig. 7(d).

Table 1 lists the root mean square error (RMSE) between the reconstructed results, $\phi$, and the true phantom, $\hat{\phi}$, for several tomographic reconstruction methods. Mathematically, the RMSE is expressed as:

[TABLE]

The RMSE value of the BP result was $0.4635$; such a high RMSE value was expected as the BP method is the simplest form of reconstruction technique and does not incorporate any filter or iterative technique in its model. The RMSE value decreases to about $0.1586$ for the FBP approach. Reconstruction by the SIRT method shows further improvement as the ring artifact, seen in the FBP result, is entirely eliminated, resulting in an RMSE value of $0.0793$. The effectiveness of the SIRT approach over analytical method like the FBP approach stems from the fact that the SIRT algorithm iteratively minimize the error sinogram as [30]:

[TABLE]

where $\lambda>0$ is a relaxation parameter. However, the ill-posed nature of the tomography problem limits the SIRT strategy to exhibit only semi-convergence [30, 31]. Accordingly, improvements seen in the first few SIRT iterations start to deteriorate as the number of iterations increases due to noise propagation. In our analysis, we show a locally converged SIRT result obtained after $10$ iterations with $\lambda=0.25$ and $W$ set as the identity matrix.

The MBIR reconstruction is obtained by iteratively minimizing the surrogate cost function of eq. 34 which, in turn, is achieved by following Algorithm 3. The $q$-GGMRF parameters used in this implementation were $p=1.1,\ T=0.001\text{ and }\sigma_{x}=0.8$. We used a $3\times 3$ non-causal weighting matrix, $b_{ij}$, to incorporate the influence of two nearest neighbors for any given pixel. The first nearest neighbors are weighted by a factor of $1/6$ while the second nearest nearest neighbors are weighted by a factor of $1/12$. The decreasing trend of the surrogate cost function over several iterations is shown in Fig. 8. An important point to note here is that the decreasing trend of the surrogate cost is a direct consequence of its strictly convex nature. Convergence after the cost minimization means that the reconstruction is globally minimal and unique. Hence, the coarseness seen along the inner region of the locally converged SIRT result is no longer visible in the MBIR result. Also, the superiority of the MBIR approach over any of the other tomography methods can be discerned from the RMSE analysis which reveals an error of only $0.0213$.

4 MBIR-based vector field reconstruction

The MBIR framework for vector field reconstruction revolves around solving an optimization problem similar in form to the one expressed in eq. 25. Specifically, we seek to formulate a framework that will iteratively work to minimize the sum of the squared differences between the data and its estimated forward projection, in combination with a regularizing prior function, $\mathcal{P}(x)$. However, in contrast to scalar tomography, vector reconstruction entails measurements from more than one tilt series; each tilt series involves a contribution from different components of $\mathbf{A}(\mathbf{r})$. We also impose a gauge constraint to reconstruct all three components of vector potential from just two tilt series. Accordingly, our generic MBIR framework from section 3 needs to be expanded to account for input measurements from two tilt series and imposition of a gauge constraint.

We begin MBIR based vector reconstruction by representing the input data by $\varphi_{m}=\left[\varphi_{m,x},\varphi_{m,y}\right]$, and subsequent estimates of the magnetic vector potential by $\mathbf{x}=\left[x_{x},x_{y},x_{z}\right]$. Let $\mathbf{H}=\left[H_{x},H_{y}\right]$ where $H_{x}$ and $H_{y}$ are operators for the forward projection of the magnetic vector potential for the $x$ and $y$ tilt series (counter-clockwise), in accordance with eqs. 2 and 3, respectively. Then the forward projection of $\mathbf{x}$ yields, $\mathbf{Hx}=\left[H_{x},H_{y}\right]\mathbf{x}=\left[U_{x},U_{y}\right]$. We define a deconvolution operator, $\mathbf{\mathcal{\mathbf{D}}}=[D_{x},D_{y},D_{z}]$, which will serve to evaluate the contributions of the magnetic phase shift (eqs. 15 and 22) to each component of vector potential $\mathbf{A}(\mathbf{r})$ separately. Thus any of the components $D_{*}$ of the deconvolution operator, when applied to the phase shift, yields $D_{*}\varphi_{m}=\left[D_{*}\varphi_{m,x},D_{*}\varphi_{m,y}\right]=\left[\varphi_{m,x,*},\varphi_{m,x,*}\right]$; and its application to $\mathbf{A}(\mathbf{r})$ yields $D_{*}\mathbf{A}(\mathbf{r})=A_{*}$.

Next, we define the prior model analogous to the one defined in section 2. We use a $q$-GGMRF as the potential function that minimizes the cost function for a given voxel based on the difference with its neighboring voxels. However, contrary to the 2D MBIR approach, we incorporate the influence of three nearest neighbors for any given voxel while evaluating the potential function $\rho(\cdot)$. Accordingly, the non-causal weighting matrix, $b_{ij}$, now contains $3\times 3\times 3$ values with first, second and third weighing factors assigned as, $9/132$, $9/264$, and $9/396$, respectively. Finally, we again make use of the surrogate majorization technique such that the $q$-GGMRF potential function assumes a quadratic form for $p\neq 2$.

Having defined all the pertinent variables, we can now express MBIR-based vector reconstruction as the solution to the following optimization problem:

[TABLE]

where $\tilde{b}_{ij}$ is determined using eq. 35. In practice, we do not directly solve eq. 39. Instead, we make use of the deconvolution operator, $\mathbf{D}$, and de-convolve eq. 39 into three MAP estimation problems as follows:

[TABLE]

Next, we identify the cost function associated with each of the potential estimates in eqs. 40, 41, and 42. For instance, the surrogate cost associated with the MAP estimate of $x_{x}$ is expressed as:

[TABLE]

Then we use the ICD algorithm, as detailed in section 3, to minimize the surrogate cost function of eq. 43. Differentiating eq. 43 with respect to $\mathbf{x}(i)$ and setting the result equal to zero, we obtain:

[TABLE]

In a similar manner, we solve for minimum $\mathbf{x}(i)$ corresponding to surrogate the cost functions of $\hat{x}_{y}$ and $\hat{x}_{z}$ in eqs. 41 and 42 to obtain:

[TABLE]

{algorithm}

[t]

1:Evaluate $\mathbf{A}(\mathbf{r})\leftarrow\text{VFET}$

2:$\mathbf{F}(\mathbf{r})\leftarrow\text{VFET}$

3:while not converged do

4:     for $i=1$ to $N$ do

5:         Determine corresponding $\tilde{b}_{ij}$ for $\{i,j\}\in\mathcal{C}$ of $A_{x}(i)$ using eq. 35

6:         $A_{x}(i)\leftarrow\frac{F_{x}(i)+2\underset{\{i,j\}\in\mathcal{C}}{\sum}\tilde{b}_{ij}A_{x}(j)}{1+2\underset{\{i,j\}\in C}{\sum}\tilde{b}_{ij}}$

7:     end for

8:     Update $A_{y}\text{ and }A_{z}$ in a similar manner similar

9:     Forward project $\mathbf{A}(\mathbf{r})$ using eqs. 2 and 3 to determine $U_{x}$ and $U_{y}$, respectively

10:     Determine the error sinogram for the two series as:

[TABLE]

11:     Use the results of imposing gauge constraint from eqs. 15 and 22 on $e\_\varphi_{m,x}$ and $e\_\varphi_{m,y}$ to determine $e\_\varphi_{m,x,x},e\_\varphi_{m,x,y},e\_\varphi_{m,x,z},e\_\varphi_{m,y,x},e\_\varphi_{m,y,y}$ and $e\_\varphi_{m,y,z}$

12:     Update the components of the vector potential as

[TABLE]

13:end while

MBIR method to reconstruct magnetic vector potential

Note that the MAP estimate of each component of the vector potential, $\hat{x}_{*}$, is determined by performing element-wise minimization of its associated cost function from eqs. 40-42 with respect to the overall magnetic vector potential, $\mathbf{x}(i)$, instead of the individual component, $x_{*}(i)$. This is primarily due to the fact that each component is reconstructed with the aid of two tilt series, $\varphi_{m,x}\text{ and }\varphi_{m,y}$, which in turn have contributions from all three components of the vector potential, $\mathbf{A}(\mathbf{r})$. Hence, each ICD iteration consists of simultaneously evaluating eq. 44-46 and then making use of the forward model and the deconvolution operator to update the estimates of the vector component as:

[TABLE]

A summary of the MBIR-based $\mathbf{A}(\mathbf{r})$ reconstruction is illustrated in the form of pseudo code in Algorithm 4. In the next two subsections we make use of this MBIR algorithm to reconstruct the magnetic vector potential of synthetic as well as experimental data sets.

4.1 Synthetic Data Set

In this section, we apply the MBIR algorithm to the reconstruction of $\mathbf{A}(\mathbf{r})$ for the cylindrical and prismatic MNPs discussed in section 2. The reconstruction was performed making use of the missing wedge projection set since this set replicates the limited angular scenario typical of TEM based experimental measurements. The $q$-GGMRF parameters used for the reconstruction comprised of $p=1.001,\ T=0.01\text{ and }\sigma_{x}=0.8$. Selected results retrieved from the MBIR approach are depicted in Figs. 9 and 10. These figures also include corresponding plots from the ground truth and the VFET approach for the purpose of qualitative comparison of reconstruction accuracy between the two approaches. We also present a quantitative comparison between the two methods in terms of NRMSE plots in Fig. 11.

A review of the plots in Figs. 4, 9, and 10 reveals that the low spatial resolutions, protrusions, ring artifacts and edge artifacts evident in the VFET reconstructions are significantly suppressed in the MBIR-based reconstructions. In case of the prismatic MNP, we note that the MBIR-based reconstruction displays values in closer proximity to the ground truth than the ones obtained from the VFET reconstructions, as shown in Figs. 9(a-c). This is confirmed by the NRMSE plots in Figs. 11a and 11c which show that the MBIR approach determines $A_{x}$ and $A_{z}$ with $10$ to $20\%$ less planar error (in $y$ axis) than those derived from the VFET approach. On planes, such as $(x,18,z)$, shown in Figs. 4(a-c), where the reconstructed vector potential suffers from strong edge artifact as a direct consequence of the limited angular range of projection measurement, the MBIR reconstruction performs significantly better than the VFET approach. In Fig. 9f, we observe that, although MBIR suppresses the ring artifact, the blurred edges persist; the loss of information due to the missing wedge is such that even an advanced algorithm as MBIR is unable to completely remove the artifacts observed in the VFET based reconstruction. This conclusion is also substantiated by the NRMSE plot in Fig. 11b which reveals that the gain in reconstruction accuracy having adopted the MBIR approach is merely $1$ to $2\%$ compared to the VFET approach.

For the cylindrical MNP, we see that the MBIR approach preserves the spatial resolution and retrieves vector potential values in closer proximity to the ground truth than the VFET approach, as depicted in Figs. 10(d-f). Accordingly, the NRMSE for $y$ planes of all three components of $\mathbf{A}(\mathbf{r})$ are less by about $10\%$ for MBIR reconstructions as compared to VFET reconstructions, as illustrated in Figs. 11(d-f). Also, edge and ring artifacts observed in the VFET reconstruction are substantially diminished in the MBIR reconstruction, as shown in Fig. 10(a-c). However, both MBIR and VFET results exhibit protrusion effects due to the missing wedge.

4.2 Experimental Data Set

In this final section we make use of a real data set to demonstrate the gain in reconstruction quality from making use of the MBIR approach. A 2D lattice of elongated $\text{Ni}_{80}\text{Fe}_{20}$ (permalloy/py) islands (“stadia”) was used for the comparative study of the MBIR and VFET techniques. The sample was fabricated on a JEOL $9300$ electron beam lithography system. A single layer of ZEP resist of $100$ nm thickness was coated on a Si/SiN substrate, followed by patterning of a square lattice with element shape parameters of $2\text{Lx}=290$ nm, $2\text{Ly}=130$ nm, and a lattice spacing of $a=390$ nm. A Py film of $20$ nm thickness was deposited on a seed layer of Cr ($3$ nm) using dc magnetron sputtering at $3$ mTorr pressure and $50$ W power. The pattern was transferred by a lift-off process. This was followed by optical lithography and wet-etching of Si to create electron transparent windows on $3$ mm square grids, which could be loaded directly into the TEM for observation. The microscopy was performed using the JEOL 2100F TEM equipped with a dedicated Lorentz lens and a spherical aberration corrector [32].

The projection measurements acquired for tomographic reconstruction consist of $x$ and $y$ tilt series with angles ranging from $-50^{\circ}$ to $+50^{\circ}$ at a step size of $1^{\circ}$. Fig. 12 depicts one of the projections as a holographic plot, $\cos(100\varphi_{m})$, of the magnetic phase shift of the sample at $0^{\circ}$ tilt. Each phase map has a resolution of $256\times 256$ pixels ($6$ nm per pixel). We implemented our tomographic reconstruction on a 3D voxel grid with $256^{3}$ nodes. First, we employed the VFET approach (as discussed in section 2) to determine an estimate of $\mathbf{A}(\mathbf{r})$ of the Py sample. These results are then used to initialize the MBIR algorithm. The $q$-GGMFR parameters used in the MBIR algorithm are $q=2.0$, $p=1.001$, $T=0.01$, and $\sigma_{x}=0.8$. The cost function in eq. 39 was then monotonically decreased over $35$ iterations to determine a MAP estimate of $\mathbf{A}(\mathbf{r})$ of the Py sample. The forward model calculation for each iteration was distributed over $24$ parallel threads using OpenMP and required about $18$ minutes to complete the iteration. Some of the $\mathbf{A}(\mathbf{r})$ results obtained from the VFET and MBIR approaches are depicted in Fig. 13.

A comparison of the plots deduced from the two methods reveals that the MBIR results show a significant gain in spatial resolution. The coarseness of the VFET result is substantially reduced in the MBIR result. Moreover, we see appropriately segmented and smoothly transitioning $\mathbf{A}(\mathbf{r})$ values in the MBIR results. The blurring artifacts observed in some of the VFET reconstruction sections are considerably diminished in the MBIR results. These reconstruction gains in the MBIR-based $\mathbf{A}(\mathbf{r})$ results are observed throughout the 3D spatial region when compared to their VFET counterparts.

Since the reconstructed quantity is a 3D vector field, we conclude this section with a 3D rendering (Fig. 14) of the magnetic vector potential $\mathbf{A}(\mathbf{r})$ resulting from the MBIR reconstruction approach. Fig. 14(a) shows a color-coded integrated (along the beam direction) magnetic induction map, derived from the magnetic phase shift by a gradient operation. The white rectangle delineates the region used for the 3D rendering of the vector field in Fig. 14(b). To reduce the complexity of the 3D rendering, the field vectors (shown as small cones with the appropriate orientation) are only drawn for the horizontal center plane going through the Py islands and four vertical planes corresponding to the dashed lines in Fig. 14(a). The yellow arrow indicates the viewing direction for the 3D rendering in Fig. 14(b). The individual islands are colored according to the in-plane direction of the integrated magnetic induction following the color wheel in (a). Note that the induction directions, which are oriented along the local curl of the magnetic vector potential, are properly oriented; i.e., the opposite rotation of the vector field around the green and red islands is clearly visible, and application of the right-hand rule for the curl results in the correct direction of the integrated magnetic induction. Furthermore, the square of islands on the left front of the image, and the two squares in the back, are all in a vortex state, i.e., the magnetization of the islands circulates clockwise or counterclockwise. The corresponding section through the vector field shows vectors which are all blue ($\mathbf{A}$ pointing down) or all orange ($\mathbf{A}$ pointing up). For the square on the front right, however, two of the islands are blue, so that there is no closed vortex and the vector field is much more complex in this area, reflecting the local frustrated state of the magnetization pattern.

5 Conclusion

In this contribution, we have incorporated the Bayesian inference-based MBIR technique into the tomographic reconstruction process for the 3D magnetic vector potential. The MBIR model is formulated by combining a vector slice-based forward model and a $q$-GGMRF-based prior model. An objective function was constructed and minimized using the Iterative Coordinate Descent technique to determine a MAP estimate of $\mathbf{A}(\mathbf{r})$. The MAP estimate shows diminished edge, ring and blurring artifacts when directly compared to reconstructions carried out with the more traditional filtered back-projection reconstruction from the VFET methodology. This qualitative conclusion is confirmed by the NRMSE values for simulated magnetic nan-particles with uniform as well as non-uniform magnetization configurations and prismatic and cylindrical particle shapes. We have implemented the MBIR technique to reconstruct the magnetic vector potential $\mathbf{A}(\mathbf{r})$ of a 2D array of Py islands. The resulting vector field was compared to an integrated magnetic induction map derived from the magnetic phase shift and is in good agreement.

Acknowledgements

PK and MDG acknowledge financial support from a U.S. DOE grant (DE-FG02-01ER45893) which supported the development of the 2D MBIR approach, and an NSF grant (DMR #1306296), which supported the final development of the 3D MBIR reconstructions. AM and CB acknowledge an AFOSR MURI grant # FA9550-12-1-0458 for financial support. The experimental work by CP was supported by the U.S. Department of Energy, Office of Science, Basic Energy Sciences, Materials Sciences and Engineering Division. Use of Center for Nanoscale Materials was supported by the U.S. Department of Energy, Office of Science, Office of Basic Energy Sciences, under contract no. DE-AC02-06CH11357. The authors would also like to acknowledge the computational facilities of the Materials Characterization Facility at CMU under grant # MCF-677785.