Real space iterative reconstruction for vector tomography (RESIRE-V)

Tomography has had an important impact on the physical, biological, and medical sciences. To date, most tomographic applications have been focused on 3D scalar reconstructions. However, in some crucial applications, vector tomography is required to reconstruct 3D vector fields such as the electric and magnetic fields. Over the years, several vector tomography methods have been developed. Here, we present the mathematical foundation and algorithmic implementation of REal Space Iterative REconstruction for Vector tomography, termed RESIRE-V. RESIRE-V uses multiple tilt series of projections and iterates between the projections and a 3D reconstruction. Each iteration consists of a forward step using the Radon transform and a backward step using its transpose, then updates the object via gradient descent. Incorporating with a 3D support constraint, the algorithm iteratively minimizes an error metric, defined as the difference between the measured and calculated projections. The algorithm can also be used to refine the tilt angles and further improve the 3D reconstruction. To validate RESIRE-V, we first apply it to a simulated data set of the 3D magnetization vector field, consisting of two orthogonal tilt series, each with a missing wedge. Our quantitative analysis shows that the three components of the reconstructed magnetization vector field agree well with the ground-truth counterparts. We then use RESIRE-V to reconstruct the 3D magnetization vector field of a ferromagnetic meta-lattice consisting of three tilt series. Our 3D vector reconstruction reveals the existence of topological magnetic defects with positive and negative charges. We expect that RESIRE-V can be incorporated into different imaging modalities as a general vector tomography method. To make the algorithm accessible to a broad user community, we have made our RESIRE-V MATLAB source codes and the data freely available at https://github.com/minhpham0309/RESIRE-V.


Introduction
Tomography has had a radical impact on diverse fields ranging from medical diagnosis 1 to 3D structure determination of proteins 2 , crystal defects 3,4 and amorphous materials 5,6 , at the atomic resolution.Despite its very diverse applications, the central problem in tomography remains the same, that is, how to accurately reconstruct the 3D structure of an object from a number of projections with noise and incomplete data.The conventional reconstruction methods include filtered back projection (FBP) 1,2 , algebraic reconstruction technique (ART) 7 , simultaneous algebraic reconstruction technique (SART) 8 , and simultaneous iterative reconstruction technique (SIRT) 9,10 , which remain popular in tomographic applications.Recently, more advanced iterative algorithms have been developed for tomography, including equal slope tomography (EST) 11,12 , nonuniform fast Fourier transform (NUFFT) 13 , generalized Fourier iterative reconstruction (GENFIRE) 14,15 and real space iterative reconstruction (RESIRE) 5,16 .In particular, RESIRE, which uses the Radon transform as the forward projection and the Radon transpose as the back projection, is not only superior to other existing tomographic algorithms 16 , but also has been used to determine the 3D atomic structure of amorphous materials 5,6 and the chemical order and disorder in medium/high entropy alloys 17 .Despite all these applications, they only deal with scalar tomography, where each voxel in a 3D reconstruction has a magnitude but no direction.However, in some important applications, vector tomography is required, where each voxel has a magnitude and a direction such as the electric and magnetic field.Over the years, several vector tomography reconstruction methods have been developed, including vector electron tomography with Lorentz transmission electron microscopy and holography [18][19][20][21][22][23][24][25] , soft and hard x-ray vector tomography [26][27][28][29][30][31][32][33][34][35] .In particular, the combination of ptychography, a powerful coherent diffractive imaging method 36,37 , and vector tomography can in principle achieve the highest spatial resolution, which is only limited by the wavelength and the diffraction signal 27,31,34 .Very recently, we have merged soft-x-ray magnetic circular dichroism and ptychography with vector tomography to image the 3D topological magnetic monopoles and their interaction in a ferromagnetic meta-lattice with a spatial resolution of 10 nm 34 .Using RESIRE 16 , our vector tomography algorithm can accurately reconstruct the 3D magnetization vector field from multiple tilt series each with a limited number of experimental projections.Furthermore, due to the experimental error, the measured tilt angles may not always coincide with the true orientations of the projections.To tackle this problem, we implement an iterative angular refinement method to reduce the tilt angle error 16 , enabling us to obtain more accurate vector tomographic reconstruction.Here, we provide the mathematical foundation and implementation of our vector tomography algorithm, termed RESIRE-V.Both numerical simulations and experimental data have been used to demonstrate the effectiveness of this vector tomography algorithm.

Methods
We begin with some setup and conventions.First, we employ Euler angles to describe the orientation of a rigid body with respect to a fixed coordinate system.For example, the orientation representation ZYX used intensively in our research fits well with vector tomography experiments: samples are rotated about the Z-axis (in-plane rotation) before a set of tilt series (rotation about the Y-axis) are acquired.The last rotation about the X-axis is helpful in angular refinement.We use the notation Z φ Y θ X ψ to represent Euler angle rotations: the first rotation is about the Z-axis by an angle φ , followed by a rotation about the Y-axis by an angle θ , and ends with a rotation about the X-axis by an angle ψ, respectively.The corresponding rotation matrix ψ is defined to be the product of three single-axis rotation matrices about the Z, Y, and X axes by angles ψ, θ and φ respectively: For short notation, we write R θ instead of R Z φ Y θ X ψ where θ = {φ , θ , ψ} (no orientation is specified).In perfect experimental conditions where there is no X-axis rotation, ψ is zero.Otherwise, ψ can be non-zero and we use angular refinement to determine ψ.The convention finishes and we move to the formulation part.

Formulation
For an x-ray beam propagating along the z direction (standard unit vector ⃗ e z = [0, 0, 1] T ), only the z component of the magnetization contributes to the 2D signals.The contribution takes either positive or negative values depending on the left or right circular polarization.In the case of rotation, we need the inner product R θ M(R † θ ⃗ r), ⃗ k to count for the contribution.Here, M = [M x , M y , M z ] is the magnetization vector field, which is a function of the Cartesian coordinate vector⃗ r = (x, y, z), and R † θ is the adjoined, and also inverse and transpose, of R θ .Adding the non-magnetic term O and taking the integral along the z axis (projection), we obtain the 2D signal: linearity of the projection operator.For efficient implementation, the latter form is preferred over Eqn. 7.
The maximum likelihood function will appear different for other types of noise; however, the famous least-squares form can still handle other circumstances because of its simplicity and effectiveness.Eqn.7 is our final form of the vector tomography formulation, and the remaining part is designing a numerical scheme to solve this minimization.

RESIRE-V algorithm
We develop our algorithm based on the real space iterative technique.Noticing that the projection operator is linear, one can construct a matrix representation for each Π θ .We assume that the projections have size n × n, and the sample gets reconstructed with thickness n.In that case, each projection matrix Π θ has O(n3 ) non-zero elements.In the case of over-constraint, Eqn.7 can be solved using the normal equation.Otherwise, in the case of under-determined system, we need to add a regularizer to prevent overfitting.When adding a damping term as a regularizer, we have an overconstrainted system again and we can solve the equation using the normal equation as in the over-constraint case.In either case, storing projection matrices is tremendously expensive since the size expands in the cubic order of the projection size.Here, to save memory usage, we do not need to store the projection matrices but compute the forward projections at every iteration instead.This procedure will increase the number of computations; however, GPU parallel computing can help reduce the computational time significantly.
Our gradient descent algorithm incorporates two steps: forward projection and back projection.For the first step, we institute our 3D Radon transform from 2D Radon transform (Fig. 1), which can be found elsewhere 39,40 .The algorithm first divides pixels in a 3D image into four sub-pixels and projects each sub-pixel individually.Specifically, at a tilt angle, we compute the corresponding coordinate of each pixel and project it on the XY plane.The value of each sub-pixel is distributed proportionally to the four nearest neighbors, according to the distance between the projected location and the pixel centers.Supposing that the pixel projection lands on the center point of a bin, then the bin on the axis gets the entire value of the pixel.If the pixel projection hits the border between four bins, the pixel value is split evenly between these four bins.
Next, we establish the transpose of the Radon transform for the back-projection step.This process is similar to the forward projection but in reverse order.According to the distance between the projected location and the pixel centers, the four nearest neighbors to a projection sub-pixel proportionally contribute their values to the sub-pixel.If the pixel projection hits the border between 4 bins, the pixel takes a quarter value of each of these four bins.
After specifying the forward and back projection, we can now take the gradient of the error metric ε(M) in Eqn.7 with respect to each magnetization component: where Π T θ is the transpose operator of the Radon transform for Euler angles θ.As mentioned above, the second form of the gradient will be used for the C++/Cuda implementation.
Next, we show that the gradient is L-lipschitz and the algorithm will converge to the global minimum with an appropriate step size.Specifically, we want to find an L such that the following inequality is true: The Lipchitz constant L gets calculated as √ 3nN z where n and N Z are the number of projections and the thickness in pixels of the reconstruction, respectively.Hence, we can choose the step size to be 1/L for the convergence guarantee.Details of the proof can be found in the Supplementary, step-size analysis.The algorithm is finalized and described step by step in pseudocode 1 and Fig. 2.
For efficient implementation, the gradient w.r.t. each component ∂ M x will be accumulated.In addition, the step size is generalized to be t √ 3nN where t ≈ 1 is the normalized step size.By our analysis, t should be less than or equal to 1 for the convergence guarantee.The analysis uses triangle inequalities and takes into account the worst-case scenario.In practice where better scenarios are more popular, the algorithm can converge with t's values slightly larger than 1.The analysis is complete and we move to the discussion on conditions for vector tomography reconstruction.

Analysis: conditions for vector tomography reconstruction
Scalar tomography, a well-posed problem, only requires one dataset (single tilt axis) for reconstruction 42 , assuming the number of measurements is sufficient.The Fourier slice theorem can verify that a 2D image needs n sampling points on the Fourier domain for a unique reconstruction.This requirement implies that n 2 projections corresponding to these sampling points are sufficient for the scalar tomography reconstruction.Indeed, the actual number of measurements is much smaller than the theoretical value.
For a 3D vector field, all three magnetization components need to be recovered, giving a question about the particular requirement in the measured data.In his discussion in the 1980s, Norton showed that the reconstruction of a diverge-less 2D vector field appeared to be unique 43 .Prince gave a more generalized discussion of the reconstructions of arbitrary vector fields in the 1990s.He demonstrated that, for reconstructing an arbitrary n-dimensional vector field, n tomographic projection datasets in which the probe is sensitive to n different directions of the vector field needed to be acquired 44,45 .The idea of using more than one tilt rotation axes has been used successfully in scalar tomography to reduce the missing wedge artifacts 46,47 .That idea is also believed to be a key to solving the vector tomography problem.
In our research, we use the Fourier slice theorem to show specific experimental conditions for the reconstruction of arbitrary vector fields.The theorem states that the 2D Fourier transform of a 2D projection equals a 2D slice through the origin of the 3D Fourier transform of an object.The 2D slice is defined based on the corresponding rotation angle.In the case of noise-free, we apply the Fourier transform to both sides of Eqn.6.
Applying the Fourier slice theorem, we have a linear constraint involving the Fourier transforms mx , my and mz of the three magnetization components M x , M y and M z .This constraint applies to every Fourier point ⃗ ξ on a 2D Fourier slice through the origin.
Algorithm 1 RESIRE-V Input: a set of n projections with left and right polarization {P + θ i } n i=1 and {P − θ i } n i=1 and their corresponding tilt angles {θ i } n i=1 , the total number of iterations K, the step size t ≈ 1 Preprocessing: 1. Compute the mean images from the left and right polarization projections: b + θ i = 1 2 (P + θ i + P − θ i ). 2. Reconstruct the non-magnetic signal O from the mean images and use threshold to obtain its support.3. Compute the difference images from the left and right polarization projections: b for k = 0, . . ., K − 1 do for i = 1, . . ., n do Compute "forward projections" Compute "back projection" for each projection difference end for Update M k+1 : Apply support constraint and other regularizers if applicable.end for Output: M K Figure 2. RESIRE-V diagram: Inputs are the differences between the left and right polarization projection and the support from the scalar reconstruction.The algorithm uses a for loop to refine the magnetization vector field M. At each iteration, it calculates the forward projections and computes their differences with the measured ones.The residuals (or differences) are back-projected to yield gradients.The algorithm will use these gradients to update the magnetization and apply the support constraint.The step size t √ 3nN is replaced by s for simplification purposes.

5/12
The extra constraint ⟨ ⃗ ξ ,⃗ n θ ) = 0 is required by the requirement that a point belongs to a plan through the origin if the inner product between ⃗ ξ and the normal vector of the plan is zero.If sufficient measurements are provided, one can sample the values of all Fourier points on the frequency domain and we can extend Eqn.11 to every point ⃗ ξ on the 3D Fourier domain.
In order to separate mx ( ⃗ ξ ), my ( ⃗ ξ ) and mz ( ⃗ ξ ) for a given 3D frequency point − → ξ , we need to find three 2D Fourier slices whose normal vectors ⃗ n θ form a linear independent system in R 3 and that go through the origin and contains − → ξ .This is impossible since the set of normal vectors ⃗ n θ that satisfies the constraint ⟨ ⃗ ξ ,⃗ n θ ⟩ = 0 lies in a linear subspace of dimension two.
This analysis differs from Norton 43 and Phatak's theoretical development 18 , which analyze the reconstruction of the magnetic vector field instead.In that case, the authors can find a linearly independent system of three equations to separate the frequency signals of the magnetic vector field B. While the first two constraints are obtained from rotations, the last constraint is found by Gauss's law ∇ • B = 0 (since B is divergence-free).The magnetization vector field is not divergence-free but has another important property: the magnetization field can only exist in a magnetic material.Hence, one can utilize a support (defined as a 3D boundary of the magnetic material) as the necessary and complimentary constraint for the completeness of a magnetization reconstruction algorithm.
Furthermore, for the case of micro-magnetic and no external dynamics at the boundary, we can add in the boundary condition that the gradient of the magnetization is parallel to surface [48][49][50] , i.e. ∂ M ∂ n = 0.In practice, since the support and boundaries are difficult to get computed exactly, one should not enforce the constraint rigidly but relax it as a regularizer instead.We add this regularizer to the minimization (7): For the regularizer part, n ∂ Ω = (n 1 , n 2 , n 3 ) is the normal vector to the boundary surface ∂ Ω of the magnetic sample and ε is the regularizer parameter.The regularizer term ∥∇M • n∥ 2 ∂ Ω only takes places on the boundary and should not affect the magnetization within the magnetic structure.In further expansion, we can write the regularizer explicitly as ∥∇M ε is tunable and should be small for non-exact support.We can even ignore this regularizer (or set ε = 0) when the support cannot be computed accurately.In contrast, we can choose large ε for larger effect if the exact support is given.The final gradient will get computed with the extra term as below: Next, we discuss the robustness of each magnetization component reconstruction M x , M y , and M z concerning the in-plane rotation angles φ .The linear constraint for in-plane rotations as in Eqn. 12 reveals that the x and y components are coupled by linear factors sin θ cos φ and sin θ sin φ while the linear factor of z is only cos θ .So the x and y parts are coupled at a higher degree than the z component.As a result, the z component will get decoupling easier and yield a high-quality reconstruction than the other two.Assuming that two in-plane rotation φ 1 and φ 2 are chose, then these two angles should be chosen equally distanced on half of the unit circle to improve the robustness of the reconstruction.We can choose φ 1 = 0 o and φ 2 = 90 o as a simple option.
Side rotations can improve the robustness of the x and y components, but that approach is experimentally infeasible.The summary of our analysis is shown below: 1. In-plane rotations are necessary but not sufficient to decouple the Fourier coefficients of the three magnetization components.

6/12
2. Other constraints, such as support and boundary constraints, and regularizers, need invoking, if possible, for highly accurate reconstruction.
3. The z component gets reconstructed with higher quality than the x and y components in in-plane rotation systems.
With the help of a support constraint, we will show that highly accurate vector tomography reconstruction can be obtained numerically with in-plane rotations.

Vector tomography reconstruction of simulated data
In this simulation, the sample is a meta-lattice with a size of 100 × 100 × 100 pixels.The signals of the magnetization make up around 1.65% of the total signal.Two tilt series from two in-plane rotations where φ = 0 o and 90 o are inspected.For each tilt series, 45 projections of each left and right polarization P + θ and P − θ are generated in the range of [−66 o , 66 o ] with an increment of 3 degrees.So totally, we generate 180 projections with size 100 × 100 pixels.To make it realistic, we add Poisson noise to the projections by selecting a flux of 4e8 photons.This flux yields an SNR of 200 and less than 1% noise.Reconstructing the non-magnetic part is not the interest of this research.However, we need to assume that the support of the non-magnetic part is given since it plays an essential role in the reconstruction of the magnetization M.
Next, we take the left and right projection difference b − θ = 1 2 (P + θ − P − θ ).Since the magnetic part only makes up a fraction of the total signal, its SNR is much smaller than that of the non-magnetic part.The SNR of the projection difference is approximately 1.65% × 200 = 3.3, which is quite small.The high noise level in the projection differences causes the reconstruction of the magnetization M to be less robust than the scalar one.Assuming that the noise level in the non-magnetic part stays the same, then the robustness of the reconstruction will decline as the magnetization signals decrease relative to the non-magnetic signals.
Now we use our algorithm to reconstruct the three magnetization components M x , M y , and M z .The model and result are shown in Fig. 3. Since the support constraint is enforced, that is, the magnetization field only appears in the magnetic material.In addition, since we use two tilt series at φ = 0 o and φ = 90 o , the missing wedge artifact does not significantly affect the reconstruction.Fig. 3d-f, shows M x , M y , and M z components in the central slice along the z-axis, which are in good agreement with the model the qualities of reconstructions in all directions are comparable and as good as the model (Fig. 3a-c).To quantify the vector tomography reconstruction, we calculate the Fourier shell correlation of the three components between the model and the reconstruction (Fig. 3g).The large correlation coefficients indicate the good quality of the vector tomography reconstruction.Additionally, we observe that the reconstructed M z has higher quality than M x and M y , which is consistent with our analysis.Fig. 3h-i show the 3D magnetization vector field of the reconstruction and the central slice along the z-axis, respectively, which agree with the model (Fig. 3l-m).We also plot two topological defects with one positive topological charge (Fig. 3j) and the other negative charge (Fig. 3k), both of which are in accordance with the model (Fig. 3n-o).All these analyses confirm that RESIRE-V can reconstruct a high-quality 3D vector field from multiple tilt series each with a limited number of projections with a missing wedge.

Vector tomography reconstruction of an experimental data of a ferromagnetic meta-lattice
The ferromagnetic meta-lattice consists of 60 nm silica nanoparticles forming a face-centred cubic structure infiltered with nickle.The vector tomography experiment was conducted at the COSMIC beamline at the Advanced Light Source, Lawrence Berkeley National Lab.Circularly polarized x-rays of left-and right-helicity were used to achieve differential magnetic contrast based on x-ray magnetic circular dichrosim 26,51 .The x-ray energy was set as 856 eV, slightly above the nickel L 3 edge Three in-plane rotation angles (0°, 120°and 240°) were chosen, and tilt ranges is from -62°to +61°for each in-plane rotation angle.At each tilt angle, 2D diffraction patterns were reconstructed using the regularized ptychographic iterative engine, producing two projections with left and right polarization at each tilt angle (Fig. 4a-b).The scalar tomography reconstruction was performed from three sets of tilt series by summing each pair of the oppositely-polarized projections, from which a 3D support was obtained to separate the magnetic materials from the non-magnetic region.For the vector tomography reconstruction, the difference of the left-and right-circularly polarized projections was calculated (Fig. 4c), producing magnetic contrast projections of three independent tilt series.Using RESIRE-V with the support, we reconstructed the 3D magnetization vector field.Fig. 4d shows the 3D vector field of the magnified square region with dotted lines in Fig. 4a-c, where the colors represent the different directions of the vectors.A thinner slice of the magnified region and a topological defect with a positive charge are shown in Fig. 4e-f, respectively.A more detailed analysis of the 3D magnetization vector field and the topological defects in the ferromagnetic meta-lattice can be found elsewhere 34

Conclusion
We present the mathematical formulation and implementation of RESIRE-V, an iterative algorithm for the 3D reconstruction of the vector field.RESIRE-V requires the acquisition of multiple tilt series of projections and the algorithm iterates between these projections and a 3D structure by using a forward and a backward step.The forward and backward step consist of the Radon transform and a linear transformation, respectively.Our analysis indicates that incorporating a 3D support to separate the magnetic region from a non-magnetic region can help RESIRE-V achieve accurate and robust reconstruction of the 3D vector field.To validate RESIRE-V, we perform a numerical simulation of the 3D magnetization vector field in a meta-lattice.Using only two tilt series and a support, we reconstruct the 3D vector field with high accuracy.We also observe that the reconstructed z component has higher quality than the x and y components, which is consistent with our mathematical analysis.Finally, we apply RESIRE-V to an experimental data set of a ferromagnetic meta-lattice, consisting of three tilt series with different in-plane rotation angles.Each tilt series has two sets of projections with left and right polarization.By using a support constraint, we reconstruct the 3D magnetization vector field inside the ferromagnetic meta-lattice, showing topological defects with positive and negative charges.We expect that RESIRE-V can be a general vector tomography method for the 3D reconstruction of a wide range of vector fields

Figure 1 .
Figure 1.Illustration of the Radon transform in the 2D case from Matlab 41 :The algorithm first divides image pixels into four sub-pixels and projects them onto a 2D plane separately.The value of each sub-pixel is distributed proportionally to two nearest neighbors, according to the distance between the projected location and the pixel centers.The transpose of the Radon transform follows the same idealogy in the reverse order.According to the distance between the projected location and the pixel centers, the two nearest neighbors to a projection sub-pixel proportionally contribute their values to the sub-pixel. .

Figure 3 .Figure 4 .
Figure 3. Vector tomography reconstruction of simulated data.(a-f) Three magnetic components at the central slice in the z direction of the model (a-c) and vector tomography reconstruction(d-f) from the simulated data where the normalized cross correlations are 94.1%, 93.8% and 99.1% for M x , M y and M z respectively.(g) Fourier shell correlation of the three magnetic components, also confirming that the z component has higher quality than the x and y components.(h-k) 3D magnetization vector field of the model, including the overall vector field (h), the central slice along the z direction (i), two topological defect with positive charge (j) and negative charge (k), where the colors represent the different directions of the vectors.(l-o) Reconstructed 3D magnetization vector filed, including the overall vector field (l), the central slice along the z direction (m), two topological defects with positive charge (n) and negative charge (o), which are in good agreement with (h-k).