1. Introduction

Throughout the history of their use, X-rays have been used to make two-dimensional, and later, three-dimensional images encoding the distribution of elements or chemical species, which are expressed as scalar quantities. However, in many cases, one is interested either in magnetic moments (vector) or crystal orientation (tensor). The imaging of magnetic structures has been of major importance in the study of spintronic and other magnetic materials, and many methods have been developed for imaging magnetization in two dimensions. One important class of such methods makes use of the phenomenon of magnetic circular dichroism (MCD), in which the X-ray absorbance of a material differs between left- and right-circular polarization by an amount depending on the local magnetization [1,2]. This effect is readily used to create two-dimensional images in various sorts of X-ray microscopes such as ptychographs [3] or photoelectron emission microscopes [4]. Note that these two papers are only examples; many more exist in the literature. However, only recently has there been work in which the magnetic order, as a three-dimensional vector, been imaged in three dimensions [5]. It should be noted, however, that ultrasound tomographic reconstruction of flow fields in fluids is mathematically equivalent to MCD tomography and was theoretically investigated before the X-ray version was accomplished [6].

Another important class of materials is crystals, in which the orientation is an important part of the microstructure and plays an important role in understanding the properties and functions of crystalline materials. For instance, teeth are made of apatite, which is not a particularly strong or durable material, but the nanostructure of enamel makes them wear-resistant enough to last a lifetime [7]. While diffractive methods such as electron backscatter diffraction and X-ray diffraction tomography can provide very useful information in two and three dimensions, they have limitations in sample thickness (electron diffraction) or resolution on the nanoscale or fail if the material is not perfectly crystalline, as is often the case for biominerals. What is wanted is a method which can image the orientation of a crystalline material, over cubic microns of volume, with sub-micron or nanoscale resolution. Although methods such as PEEM-PIC mapping [8] and dichroic ptychography [9–11] provide the requisite resolution and sampling area in two dimensions, it can be very difficult to image orientations in three dimensions.

Dichroic tomography is an extension of the well-known technique of X-ray computed tomography, in which a sample is rotated in an X-ray beam and, at each of a set of rotation angles, a projection image (radiograph) is captured. This set of images is processed to produce a 3D reconstruction of the absorption density of the sample. In dichroic tomography, the incident beam is polarized, and so additional contrast is provided due to the crystal orientation (linear dichroism) or magnetization (circular). Magnetic linear dichroism, as observed in antiferromagnets, is equivalent to the linear dichroism seen in uniaxial crystals (see, for example [12]). In this case, orientation is replaced by antiferromagnetic order parameter with no change in the formalism and so is covered by discussion of X-ray linear dichroism (XLD).

While the object of dichroic tomography is to reconstruct a tensor at each voxel in a sample, this technique should not be confused with another, also called tensor tomography, in which scattering patterns from many projections are used to reconstruct 3D scattering patterns at each voxel. These patterns are often parameterized as tensors. [13–15]. This technique can provide information about the orientation of fibrils, but that is not equivalent to the orientation of crystals, as shown, for instance, in studies of teeth [7].

The purpose of this paper is to discuss the theory of how the properties of the sample map onto the projected images, and what measurements would be required in order to make it possible to reconstruct not only an overall absorption at each voxel, but the anisotropy. While the development is done in the context of reconstruction by filtered back-projection, the purpose is not to propose a reconstruction algorithm but to ask what would be required to make such an algorithm possible. This leads to considerations such as:

The next section goes over a number of definitions and discuss the framework and formalism used for the math. Then I’ll specify some assumptions being made and what they imply, including a review of filtered back-projection, after which comes the actual derivations.

This paper was originally written in the context of soft X-ray scanning transmission X-ray microscopy (STXM), but the results should apply to hard X-rays and to full-field imaging methods.

In the course of the development, some results will be derived which were previously known. These serve as examples and checks on the mathematical methods.

2. Preliminaries

This section is where definitions and formalism are laid out.

2.1 Assumptions

First, some assumptions:

The result of a normal tomographic reconstruction is the volume absorbance of the set of voxels defined by positions in the sample, $\mu (x,y,z)$ . Thus, the beam intensity after passing through a small thickness dl of material at this position is attenuated by a factor of $1 - \mu (x,y,z)dl$. For an isotropic material, this quantity is independent of beam direction or polarization, an assumption implicit in all standard tomographic reconstruction algorithms. The assumption of Beer’s Law, with an effective absorption at each voxel, which may depend on direction and polarization, is fundamental to most, if not all tomographic reconstruction algorithms.

An important way in which this assumption is violated is if the polarization of the beam changes as the beam propagates through the material. As an extreme example, consider a material which is so strongly dichroic that a small thickness of it acts like a polarizer. If the incident polarization is vertical but a grain has a vertical absorbing axis so it acts like a horizontal polarizer, the transmitted intensity is low. If, however, a grain oriented at 45° is placed before the original grain, the interposed grain turns the polarization, and now up to 1/4 of the incident beam goes through. If the order of the grains is switched, the transmission drops again. This shows that a simple line-integral treatment of absorption does not work in strong dichroism conditions. The transmission is not a scalar, but a tensor described not by a line integral but by a product of non-commuting matrices. This problem does not occur in the MCD case because the polarization remains left or right circular as the beam passes through the sample [16]. The propagation of X-rays through a linearly-dichroic material is described in a number of places including [9].

What can be done to make A1 true or nearly true? One way is to operate in a regime of weak dichroism and birefringence. This of course can cause trouble with noise and artifacts, but at least guarantees that the polarization remains nearly unchanged through the sample. Note that operating at an energy below or near the dichroic peak, even if the dichroism of absorption is small, may not be appropriate because the sample may still be birefringent enough to act like a stack of waveplates and scramble the polarization. In analysis, it may be possible to adopt some kind of tensor formalism, especially if real and imaginary parts of the transmission are known as would be the case if the projections are derived from ptychography. It may also be possible to come up with an algorithm incorporating a multislice model as the forward model in an iterative loop.

The next assumption refers to the use of ptychography:

In the presence of linear dichroism, this assumption should be applied carefully. The reconstruction of the transmitted field in principle needs to take into account the tensor nature of the index of refraction [9]. However, once the reconstruction is done, the squared magnitude of the transmitted field plays the role that transmission does in STXM. Similarly, the phase shift maps onto absorption in STXM (see for example [17]).

Now comes an assumption which may not be necessary but makes the math simpler:

The next assumption involves terms which will be defined below, so readers who don’t know what back-projection is should read to that section and come back. This assumption is:

This is actually a theorem, called in the literature the Fourier slice theorem [18], rather than just an assumption. The proof is as follows: The contribution to the back-projection from a single pixel of any of the projections is a line in real space, extending along the beam direction corresponding to that projection. Thus, its Fourier transform is a sheet in Fourier space, oriented perpendicular to the beam direction and passing through the origin. Now consider the Fourier transform of a single back-projection, which is the contribution to the back-projection from a single rotation. Since the lines from all points on the detector are parallel, they contribute to the same sheet in Fourier space, which now contains the Fourier transform of the projection. Since all the sheets are unique, it follows that the Fourier transform of the back-projection can be unfolded to give the Fourier transforms of each projection, which can be inverted to yield the set of projections. For a finite number of rotations, the sheets intersect each other at a set of measure zero compared with any one sheet, so a point on a sheet corresponds uniquely to a point on a projection’s Fourier transform.

Finally, an assumption which provides an important simplification for dichroic CT/CL:

Thus, we do not consider cone-beam imaging or imaging with a focused probe for which the thickness of the sample is significant compared with the depth of focus.

2.2 Definitions

In the following section we provide a number of definitions used throughout the paper. All unit vectors will be written with hatted symbols, e.g. $\hat{n}$ (beam direction), while non-normalized vectors such as a position $\vec{X}$ will be written with over-arrows.

The method consists of rotating the sample in the beam while taking radiographs. Each radiograph is processed assuming Beer’s Law to produce the projected absorption along the beam path. These are called projections, and are parameterized by the rotation angle γ and any other variables such as polarization. For phase contrast (e.g. from ptychography), the projected phase shift is the quantity given in the projections.

If the rotation axis is perpendicular to the beam, then we have computed tomography (CT) geometry and each projection may be regarded as a stack of projections of 2D slices. Otherwise, we have computed laminography (CL) geometry and no such slicing is possible – the whole volume must be reconstructed at once. While CL is more complex than CT, and as we shall see, has an intrinsic limitation due to missing information, it is often the only way to work with samples that are sheet-like in geometry, as CT will result in projections in which the sample is probed on-edge. This case results in ranges of rotation angles in which the data are useless because the sample is effectively too thick to see through. Also, for circular dichroism, CL offers the potential of detecting all three components of magnetization, while CT misses the component along the rotation axis.

Two coordinate systems will be used. One is the lab coordinate system (or reference frame), in which the beam comes from a fixed direction, the detector is in a fixed position, and the sample rotates. The other is the sample coordinate system (or reference frame) in which the sample is fixed and the source and detector pirouette around it. This is the coordinate system in which we want the eventual reconstructions to be defined.

It is assumed (A1) that the absorption through the sample (a projection) is described for each angle and voxel by a local effective absorption coefficient

The method of reconstruction we will consider involves the creation of back-projections, made by adding up single back-projections for each rotation angle. A single back-projection, parameterized by the rotation angle is made by projecting a projection back through the sample to create a 3D density:

The back-projection is the sum of the set of single back-projections, and looks like the object with a blurry halo around it. This halo, at least in the isotropic case, can be removed by a Fourier filtering procedure, resulting in a filtered back-projection as an approximation of the original sample.

To get from a set of projections or a back-projection to the unknown parameters of the sample, we must use some sort of mathematical method which essentially does a large linear-algebra problem. Reconstruction by filtered back-projection (FBP) is the core of one family of methods used for isotropic samples. The other major family of methods considers the set of local absorptions per pixel as a large vector of unknowns, which is related by a linear matrix operation to the set of projections, considered as a larger vector of observations. This is essentially a linear least-squares problem, but the response matrix is huge. The vector of unknowns has as many elements as there are voxels and the vector of observations has a length equal to the product of the number of angles and the number of pixels on the detector. Thus, for a small case of a 200 × 200 × 200-voxel sample volume, 360 rotations and a 200 × 200 pixel detector, we have 8 × 10⁶ unknowns and 1.44 × 10⁷ observations, so the matrix would have 1.16 × 10¹⁴ elements. Even on a supercomputer cluster, it would be impractical to store and invert such a matrix, even though it is sparse. Instead, iterative least-squares methods such as SIRT (Simultaneous Iterative Reconstruction Technique [19]) and LSQR [20] are used. Although these methods are better in some senses than FBP methods and have largely taken over the field, the argument above shows that back-projections contain all the information there is, so FBP methods can be used to investigate the conditions under which there is enough information to reconstruct anisotropic samples.

We will see that the local absorption coefficient varies with rotation angle according to a trigonometric series including sines and cosines of the rotation angle and twice the rotation angle. If one wants the coefficients in such an expansion, one way to do it is to average the function with a weight of $\sin \gamma ,\cos \gamma ,\sin 2\gamma ,or\;cos2\gamma $. Similarly, the projections $p(\vec{x},\gamma )$ will include trig-function variation, so one can get all the available information by taking averages over $\gamma $ such as $\left\langle {p(\vec{x},\gamma )\cos \gamma } \right\rangle $ with the angle brackets denoting an angular average over rotation angle $\gamma $ . These shall be denoted as weighted back-projections, and the trig-function coefficients of $\mu (\vec{x},\gamma )$ as moments. For some of the analysis, it will be more convenient to expand in terms of $\exp (im\gamma )$, with m ranging from -1 to 1 (MCD) and -2 to 2 (XLD). This is, of course, equivalent to an expansion in sines and cosines.

Much of the mathematical work was done using Mathematica [21]. In particular, the MatrixRank[] and Eigenvalues[] functions were used in order to evaluate the number of independent equations represented in some linear systems.

3. Dichroism framework

At this point, we can start to create the framework in which we work. First, I will define some geometry. This is where assumption A5 (beam direction is everywhere the same) comes in. Figure 1 shows the assumed geometry, which shows the general CL case and the notation for the relevant angles. This geometry was originally considered for a tomographic electron-microscope sample holder, but is general enough to be applied to other mechanical arrangements.

 

Fig. 1. Geometry of the measurement. The unit vectors $\hat{x},\hat{y},\hat{z}$ define the sample coordinate system with $\hat{z}$ the rotation axis. Lab coordinate axes are $\hat{X},\hat{Y},\hat{Z}$ with $\hat{Z}$ being the beam direction.

Download Full Size | PDF

The three angles provided by the mechanics are $\alpha $, the angle between the rotation axis and the beam, so $\alpha = 90^\circ $ is CT mode., $\beta $, a roll angle about the beam axis, and $\gamma $ the rotation angle, which is the only angle that varies between projections. These are essentially Euler angles, but there are other possible geometries, depending on the mechanics of the measurement, in which the rotations are not Euler angles. Further, the polarization in XLD measurements is defined in terms of $\chi $, the angle between the polarization vector and the horizontal ($\hat{X}$). The entire sample holder rotates about the X axis ($\alpha $), the plate holding the grid wags up and down by angle $\beta $, and the sample, idealized as a thin object such as a TEM grid, rotates in its own plane by angle $\gamma $. The rotation specified by $\beta $ is equivalent to rolling the sample about the beam axis. This action can be mimicked by rolling the polarization and detector about the beam axis, so $\beta $ has no non-trivial effect. Thus, we will take it to be 0 in the following. Now we can start to assemble the $(\hat{x}\hat{y}\hat{z})$ orthonormal triad in lab coordinates. Start with $\alpha = 90^\circ ;\;\beta = \gamma = 0$. Inspection of Fig. 1 shows that the $(\hat{x}\hat{y}\hat{z})$ triad in lab coordinates is:

Here and in the next few equations, thestack on the left-hand side is a stack of row vectors, so is a matrix. The determinant of this matrix is 1, showing that the parity of the sample coordinate system as illustrated in Fig. 1 is the same as that of the lab system. The first rotation to apply is $\alpha $, which differs from the other two in that it is done in the lab system, since the whole holder moves. This is a rotation about X, which mixes $\hat{Y}$ and $\hat{Z}$, so we have

Note that since we started with $\alpha = 90^\circ $, we have to adjust the right-hand rotation matrix to refer to a rotation of $90^\circ{-} \alpha $. Figure 1 shows the geometry with $\alpha$ between 0 and 90° and we see that Eq. (4) gets the sign right.

The next step is to apply the in-plane $\gamma $ rotation which mixes $\hat{x}$ and $\hat{y}$ :

Using Mathematica to do the algebra results in the final matrix:

This is the matrix which, when multiplying a vector in lab coordinates, yields its equivalent in sample coordinates. Even if the assumptions about the order of rotations is wrong, the use of this rotation matrix should still let us determine whether a given set of measurements is enough to perform reconstructions. The inverse of this matrix is its transpose.

Next, we will need the polarization and the beam direction, expressed in the sample frame of reference, in which $\hat{x},\hat{y},\hat{z}$ are the basis vectors. The beam direction, which I have been denoting by $\hat{n}$, is given by the third column of the rotation matrix:

Note that if we subtract the projections done at one polarization from corresponding projections done at the other, we remove ${\mu _{nd}}$, leaving the dichroic contribution. Decomposing this into single trig functions, we have

Here, ${m_{x,y,z}}$ are the components of $\vec{m}$ in sample coordinates. We see that in CT mode ($\alpha = 90^\circ $), we lose information about one of the components of the magnetization, ${m_z}.$ This is because in CT mode, the sample z axis is always perpendicular to the beam, so there is no dichroism due to magnetization in that direction. In CL mode, if we can get the three moments we can get the scaled magnetization because unless $\alpha $ is $0^\circ $ or 90$^\circ $, the set of linear Eq. (11) is non-singular.

For linear dichroism, the formalism is more complex. Many materials of interest are optically uniaxial due to their crystallinity (calcite, for instance) or nearly so (aragonite and vaterite). Similarly, antiferromagnetism induces a linear dichroism on resonance. If we assume uniaxiality with an orientation unit vector $\hat{p}$, then the absorption may be described according to Malus’ law as

Note that with Mathematica, the easiest way to extract these coefficients is to evaluate weighted integrals such as ${\mu _2} = (1/2\pi )\int_0^{2\pi } {\mu (\gamma )\exp ( - 2i\gamma )} $.

From the above, we see that from the ratio of, for instance, ${\mu _2}$ to ${\mu _1}$ we could extract ${p^{\prime}_z}$, and from ${\mu _{ {\pm} 1}}$, ${p^{\prime}_{x,y}}$ and then from ${\mu _0}$, the non-dichroic background ${\mu _{nd}}$. This only works in the case of laminography ($\alpha \ne 90^\circ $) and an oblique polarization. Thus, if we can somehow reconstruct the moments, we can get the information we want, in most cases without needing to switch polarizations (XLD) or tilt angles (MCD). We will see later on, however, that we can’t actually get all the moments without adding extra information.

There is another way to describe XLD which is mathematically more elegant than the above, but does not easily enforce the assumption of uniaxiality. This description also covers the more general case of biaxial materials such as orthorhombic or monoclinic crystals. The absorption is described by a symmetric tensor, which is proportional to the imaginary part of the dielectric tensor:

The advantage of this representation is that, unlike the uniaxial vector representation, the equations are linear in the unknowns so matrix algebra can be used, just as in the MCD case. Also, there’s no ${\pm} $ to keep track of. The moments are thus,

For the MCD case, if we can measure all three moments, we can get all components of the magnetization because we can solve for $\vec{m}$ in Eq. (21), provided $\alpha $ is not 0° or 90°. However, for XLD without some additional assumption or information such as perhaps another polarization, we can’t solve for $\vec{m}$ because there are 6 unknowns and only 5 observables. One possibility is to constrain to uniaxiality. Unfortunately, this constraint is a non-linear condition, requiring two of the eigenvalues of Q to be equal. Such a condition can’t be implemented simply by, for instance, adding another row to S. Another possibility is to take another scan at an energy at which there is no dichroism, then assume based on the composition of the sample a scaling relation between the observed absorption coefficient at that energy and the isotropic part at the energy at which the dichroic CT/CL is performed. That effectively gives the trace of Q, adding a sixth component to the moments vector and a sixth row to S. We can test this by adding another row to S consisting of (1,1,1,0,0,0) and asking Mathematica for the rank of the now 6 × 6 matrix. The result is 6, showing that adding the constraint on the trace of Q fixes the under-determination problem, given knowledge of the moments. However, this only works if the sample is homogeneous and its optical properties are very well-known, perhaps from measurements on a polycrystal.

If we use another polarization, then we have 10 observables and 6 unknowns, so one might think that there are plenty of constraints. Actually, the problem is just determined. To see why this is so, we can set up the equations for two polarizations:

At this point, it would appear that we can find all three components of the scaled magnetization (CL mode only) or all six components of the absorption tensor for XLD using two polarization (LCP/RCP for MCD, two angles for XLD). However, this is contingent on being able to measure all three (MCD) or five (XLD) moments. We will see that this is actually not possible; the set of projections does not yield enough information. This deficit is over and above the missing cone in Fourier space that’s inevitable in CL and which we will discuss later.

4. Back-projections and filtered back-projections

So far, it is now established that for MCD, reconstructing three moments of the absorption will provide all the information we need. For XLD, five moments would be sufficient if another constraint, such as uniaxiality, could be imposed. Now we need to figure out how to do this reconstruction, and whether and under what circumstances it’s possible. We have already been using assumption A1 (ability to describe transmission in terms of a local absorption) and A5 (collimated beam). From here on, we will use A3 (no missing angles) and A4 (back-projections map into projections). However, before getting into the details of reconstruction, I will digress a bit to discuss the missing information in laminography, even for isotropic materials.

4.1 Missing information in laminography

One way to see that there must be missing information in laminography is to consider the method applied to a sample which is layered normal to the rotation axis but is uniform in-plane. For this case, rotation of the sample produces no variation of absorption; all projections are the same. Therefore, the layering is completely unresolved and all one gets is the total absorption through the sample.

The effect may be described by considering the back-projection in Fourier space. Recall that the back-projection is the sum of many single back-projections. Each of these has the form of a bundle of lines of density extended along the beam direction for the given projection. The Fourier transform of such a bundle is a sheet oriented perpendicular to the beam direction. When the sample is rotated, the sheet, in the sample’s reference frame, rotates about the sample rotation axis, resulting in a double cone within which there is no intensity. Thus, Fourier components within this cone are not sampled and no information is available about them. The mathematics of this phenomenon are shown explicitly in [22] and [23]. Examples of the reconstruction artifacts produced by this missing cone are shown in [24] in their demonstration of cone-beam laminography on a laboratory apparatus. Note that although the missing cone was demonstrated mathematically in the context of filtered back-projection reconstruction, the effect is independent of reconstruction method because it’s a real gap in the available information. Theorem A4 shows that no reconstruction algorithm will dig up information that filtered back-projection can’t. As an example, the reconstructions shown in [24] were done using an iterative least-squares algorithm, not filtered back-projection, yet show missing-cone artifacts. However, with prior information about the unknown sample, it is possible to mitigate these effects. For example, it is shown in [25] that knowledge of the minimum and maximum values of the refractive index (they reconstructed phase contrast) can be used with an iterative modification of FBP to reduce the missing-cone artifacts.

4.2 Reconstruction of moments

We have seen that the unknown object may be described in terms of effective absorption coefficients which oscillate as a function of rotation angle $\gamma $. Suppose we want to derive a generalization of filtered back-projection which can reconstruct not just an angle-independent $\mu $ but all the moments $\{ {\mu _i}\} $. The back-projection corresponding to a single point of isotropic absorption at position $\vec{x} = {\vec{x}_0}$ looks like a double cone of non-zero values with the vertex at position ${\vec{x}_0}$, oriented along the rotation axis. The intensity is concentrated at this vertex because the lines from the different angles crowd together, thus producing a power-law singularity. This is why the unfiltered back-projection resembles the object but with a ‘haze’ around each point. Now consider the (unphysical) case in which only one of the non-zero moments appears, ${\mu _i}(\vec{x}) = \delta (\vec{x} - {\vec{x}_0}){\delta _{im}}$ . The back-projection produced by the “object” is still a double cone, but modulated with angle by a factor of ${e^{im\gamma }}$. Unlike the back-projection from an isotropic point absorber, the double cone with $m \ne 0$ does not have a concentration at the vertex. The appearance of the back-projections of point objects was verified by simulation using the ASTRA toolbox [26]. If we perform the back-projection of an object having intensity in all the moments ${\mu _m},m \in [ - 2,2]$ (XLD) or $m \in [ - 1,1]$ (MCD), we get a complex mixture of contributions from all the ${\mu _i}$ . The ${\mu _0}$ component will still be concentrated at points where ${\mu _0}$ is high, but simulations show that the other moments contribute significantly.

Now, if we perform the back-projection by first weighting each projection by ${e^{ - im\gamma }}$, we change the ${\mu _m}$ moment into what looks like ${\mu _0}$, somewhat the way a lock-in amplifier converts an oscillating signal into a DC one by mixing it with a reference of the same frequency. Thus, we could imagine performing a series of filtered back-projections on weighted projections, then solving for the moments. This process is illustrated graphically in Fig. 2. In this example, real-valued trig-function weightings are shown, rather than complex exponentials.

 

Fig. 2. A graphical representation of the effect of performing $\sin ,\cos (n\gamma )$-weighting on the contributions to weighted back-projections from the components of the absorption density. The double cones represent weighted projections of a point object with positive values shown as grey (brighter = higher), negative values shown as red and 0 as black. Note that the resulting angular components in the projections go up to n = 4. The weightings are shown as arrows connecting moments of the absorption (the circles on the left) with the resulting contributions to the back-projections.

Download Full Size | PDF

Here, the uniformly-shaded disk represents isotropic absorption, the single-headed arrows represent moments of order 1 and the lines with no arrows represent moments of order 2. The double cones represent back-projections of point objects with varying angular dependence. We see that the contributions of moments of order up to ${\pm} 2$, when weighted, yield contributions to order ${\pm} 4$. Thus, in order to go from observables (the back-projections) to unknowns (the moments), we must perform what is effectively a matrix inversion.

4.2.1 Generalized filtered back-projection

In ordinary filtered back-projection, the filtering is done by computing the FT of the double-conical back-projection and using its reciprocal (where the FT is not 0) as a deconvolution filter. Here, we generalize to a multicomponent problem. For this purpose, it is easier to work with complex components, as in Eq. (17), than sines and cosines as in Fig. 2. The weighted projections from an isotropic point object at $\vec{x} = 0$ are given by:

The second integral is one of the many integral forms of the J Bessel function, so:

The form for m = 0 is the one presented by [23]. From there, ${F_m}$ for non-0 m can be evaluated by integration by parts and the use of a recursion relation:

We need to be careful in evaluating F for ${k_z} = 0.$ At ${k_t} = {k_z} = 0$, the function is undefined because ${\lim _{{k_t} \to 0}}{\lim _{{k_z} \to 0}}{F_m} = 0$ but ${\lim _{{k_z} \to 0}}{\lim _{{k_t} \to 0}}{F_m}$ is infinite. It seems to be standard to take the value of the reciprocal transform ($F_m^{ - 1}$ in our case) to be 0 at the origin of k-space. For ${k_z} = 0,{k_t} \ne 0$ the integral is:

Similarly, in CT mode, ${k_z}$ is 0 everywhere, so ${F_m}$ is also given by Eq. (29). Now that we have the FTs. we find that the FT of the weighted back-projections is:

There is a special case which needs to be considered, which is that of $\vec{k} = 0$, which will be one of the points sampled in the discrete FT. The value of the FFT at the 0 point is the sum over the imaged volume of the operand. Consider a ‘sample’ which is uniform within some volume outside which it’s empty, and has only an m’th moment, which means that it has an unphysical $\mu = \exp (im\gamma ).$ Then, for the left-hand side of Eq. (32) we have ${\tilde{\mu }_n} = {\delta _{nm}}$ and ${\tilde{p}_n} = C{\delta _{nm}}$ where C is a constant. To evaluate this constant, it’s probably easiest to simulate the back-projections through a uniformly-filled volume and sum up all voxels of the back-projection. Since the proportionality of ${\tilde{\mu }_n}$ and ${\tilde{p}_n}$ holds for all values, it follows that the matrix in Eq. (32) must reduce to the identity matrix times C for the $\vec{k} = 0$ point.

All of the above is aimed at the problem of reconstructing the moments ${\mu _m}$ of the absorption coefficients. From that, we can go back to the formalism starting at Eq. (9) to get the material properties represented by $\vec{m}$ or ${Q_i}$. For instance, the development leading to Eq. (17) relates the moments to a solution for uniaxial linear dichroism with a single oblique polarization in laminography, while that leading to Eq. (11) shows that the scaled magnetization can be solved for with two polarizations in CL mode.

Summarizing, we find that if we can reconstruct the moments ${\mu _m}$ of the absorption coefficients, for instance by a weighted back-projection method, we could reconstruct the orientation and magnitude of the dichroism. With what set of measurements can we pull out all the moments uniquely?

4.2.2 Missing information

There is a major problem with the above approach if one tries to use it with a single polarization and tilt. It turns out that the $\tilde{G}$ matrix in Eq. (32) is singular. That means that it’s possible to add to the ${\tilde{\mu }_m}$ coefficients an extra contribution ${C_m}(\vec{k}){u_m}(z,t)$ where ${C_m}(\vec{k})$ is any function in k-space and ${u_m}(z,t)$ is one of the eigenvectors corresponding to a 0 eigenvalue (there is one for the 3 × 3 MCD case and three for XLD) of the matrix to be inverted. Several speculative workarounds can be proposed. A possible solution is to ask that outside of some support volume, all components of $\mu $ are zero. How would we get that support volume? Perhaps one could use a conventional method of reconstruction other than filtered back-projection, maybe at an energy at which dichroism is negligible. Some version of shrink-wrapping might also be possible. Thus, we can frame the problem as a minimization problem, with the norm of the absorption tensor outside the support volume as an additional cost function.

For XLD, one extra bit of information we have is that for many, if not most cases of interest, the optical properties are uniaxial, so the absorption coefficient tensor is given by Eq. (19). Counting up degrees of freedom yields four independent variables, i.e. two constraints. Since the matrix for XLD has three vanishing eigenvalues, it appears that uniaxiality does not provide enough constraint to remove the degeneracy. Also, the relationship between the weighted moments, in real space, is not simply linear, so linear methods don’t work, as shown in Eq. (17).

Suppose we add polarizations to get more information in XLD. As we did previously in discussing the connection between material parameters and moments of $\mu $, we can make a stack of the matrices connecting $\{ {\tilde{p}_m}\} $ and ${\mathbf Q}$ (Eqs. (20) and (32)) and test for rank. The resulting equation, in block-matrix form, is

It should be noted that since the matrix has more rows than columns, we’re now in the domain of a least-squares solution, though one of manageable size for each k-point, rather than the global one with hundreds of trillions of elements. We can thus solve the problem in one pass, rather than needing an iterative process.

For MCD, we do have a missing-information problem as was noted by [27], who showed that CL or CT done at a single orientation does not yield a unique solution for the magnetization. That result is reproduced in this work. While it’s true that the rank of ${S_{MCD}}$ is 3 for CL, the rank of the 3 × 3 matrix corresponding to the G-matrix in Eq. (32) is 2, so we can’t solve for all components of $\vec{m}$ . We can’t add polarizations because in MCD we only have two, and we’re using both to get rid of the non-dichroic term. This rank is preserved on multiplying on the right by ${S_{MCD}}$ . Thus, we need a constraint. On certain length scales, many materials have magnetizations either only in the plane of a film or perpendicular to it. If the film plane is the rotation plane, then the perpendicular case corresponds to $\vec{m} = (0,0,{m_z})$. In this case, tomography as we have been discussing it fails to deliver 3D information. If a CT geometry is used, the magnetization is perpendicular to the beam direction, so no dichroic signal will be observed. However, this case is somewhat artificial; the domain walls or other boundaries between areas of differing ${m_z}$ have local moments out of the normal direction, so the only information we can get comes from places where the constraint is violated. For the in-plane case, the nullspace of the $\tilde{G}S$ matrix contains one vector, which has components along $\hat{x},\hat{z}$ directions. This implies that anything one might add to $\vec{m}(\vec{x})$ that still satisfies the linear equation of the reconstruction will inevitably have a component in the $\hat{z}$-direction, which we’re not allowing. Unfortunately, many of the interesting structures we would want to look at such as skyrmions or hopfions obey no such constraints, as shown for instance in Fig. 3 of [27].

 

Fig. 3. A demonstration of an iterative reconstruction algorithm. A) Original (top) and reconstructed (bottom) orientations for the random polycrystal. The orientations are coded by the color map on the right-hand side, where $\theta $ is the polar angle of the optic axis relative to $\hat{z}$ and $\phi$ is the azimuthal angle of the optic axis relative to $\hat{x},\hat{y}$. The strength of the dichroism is coded as the lightness of the colors, which is why the isotropic material in the center shows as black. B) The same, but for the ‘spherulite’ phantom (orientations are radial). The four conditions used were $\{{(\chi ,\alpha )} \}= \{{(0,60^\circ ),(45^\circ ,60^\circ ),(90^\circ ,80^\circ ),(126^\circ ,80^\circ )} \}$.

Download Full Size | PDF

One way around the problem is to use multiple tilt angles, by analogy with the multiple polarizations and tilts used for XLD. The calculation must be done carefully because the matrix $\tilde{G}$ depends on the tilt angle for a given point in Fourier space so each block in the stack has its own $\tilde{G}$. On doing the same kind of calculation for MCD as above, but with two different values of $\alpha $, we find a matrix rank of 3, meaning that two tilts are necessary and sufficient. One group [28,29] varied both $\gamma $ and $\alpha $, but they inverted the role of the two angles, using only two $\gamma $ angles but scanning $\alpha $. Another group [27] performed magnetic tomography using a similar method to what’s discussed here, with two tilt angles and 180 rotations. For XLD, doing the same calculation for a single polarization shows that we need three tilts. The singularity is removed with three tilts only if the polarization is neither horizontal nor vertical.

5. Conclusions, demonstration and discussion

I have laid out a framework for reconstructing laminographic images of anisotropic samples under conditions of variable polarizations, as a generalization of filtered back-projection. Using this, it was possible to determine that for XLD, tomograms at at least three polarizations are required for full information to be retrieved. Further, the tilt angles for all the polarizations must not be equal. Alternatively, three tomograms at different tilt angles α measured at a common polarization other than horizontal or vertical provide enough information. This result is relevant for the use of dichroic tomography on beamlines whose sources are bending magnets or horizontally-polarized undulators, both of which have polarization fixed in the horizontal plane relative to gravity. Thus, to perform multi-tilt XLD tomography, the α axis must be tilted at an angle to the horizontal, i.e. β must be non-zero. For MCD, it is necessary to image at at least two tilt angles. These conditions hold even if algorithms other than filtered back-projection are used, because the back-projections hold all the information available.

We can summarize the above math by assembling a weighted FBP algorithm by analogy with the standard scalar algorithm. The simplest scalar back-projection algorithm may be summarized as:

This method has a number of drawbacks, including the inefficient use of a 3D FT-inverse FT pair, but it serves to illustrate the point. The corresponding algorithm sketch for a tensor order parameter is:

Filtered back-projection is a fast algorithm because it reconstructs in one pass, but it tends to do poorly with regard to missing angles and noise and artifacts from the missing cone in Fourier space. Therefore, an iterative reconstruction algorithm, similar in principle to what has been done for MCD [5,27], was developed. In this approach, the analysis is carried out entirely in real space, as a linear least-squares problem, using the LSMR algorithm [30] as implemented in Scipy [31]. One advantage of this family of algorithms is that regularization is available, which helps cut down on artifacts. In this approach, Eq.(1) is discretized:

In order to demonstrate that this algorithm works, two simple phantoms were simulated. One of these models a polycrystal. Ten points are randomly chosen in a cube, and the Voronoi neighborhoods of these are ‘colored’ with a uniaxial absorption tensor with a random orientation. The other phantom represents a spherulite. Here, a sphere is filled in with a uniaxial orientation with axis along the radial direction, except for a smaller sphere in the center which is filled with isotropic absorber. In both cases, the absorption along the optic axis is 2 (arbitrary units) and that perpendicular is 1. The absorption outside the sphere or cube is 0, representing empty space. The units can be arbitrary because the whole process is linear. When noise or other real-world effects are considered, this will not be the case. Figure 3 shows the comparison between the reconstructed and original phantoms in three slices. The orientation is mapped to color via the map shown, which is a modification of the one used in [7]. A similar exercise was carried out using the FBP algorithm sketched above.

The current algorithm does not include regularization, which will be helpful in the presence of noise. This could be remedied, at the cost of more computational effort, using, for instance, the normal-equations solver implemented in Pylops [32] or other methods as used in scalar or MCD tomographic reconstruction. The main difference between dichroic and scalar reconstruction is that the input data and unknowns are more complex, having indices other than just pixel/angle (projections) and voxel (density to be solved for). Note that a fixed support region can be implemented by modifying the code that maps between the single index $\{ qxyz\} $ and the four indices $q,x,y,z$ so that only voxels within the support volume are considered, with the rest of the voxels filled in with zeros.

The polarizations and tilts (listed in figure caption) were not chosen with any specific criterion in mind. We see that an accurate reconstruction is obtained.

The FBP algorithm was sketched out only to validate the use of FBP as a tool for examining the information content in a given set of projections. The iterative algorithm is more of a practical tool, though it has not been tried on actual data.

In laminography, there is a loss of information in a cone in Fourier space which is, again, independent of the reconstruction method used. It is possible that constraints such as shrink-wrapping, uniaxiality (XLD) or non-negativity could help with this problem. Another possible approach is to use machine-learning techniques to constrain reconstructions to what’s plausible, based on a training set of simulated tomograms constructed to have features similar to those found in real sample. This method was used to improve signal-to-noise and reduce required dosage in tomography by [33]. One may speculate that similar methods could be used to aid reconstruction in the case where an insufficient set of polarizations or tilts had been used for full recovery of the orientation. The training set must be carefully chosen so as not to impose constraints that the real sample doesn’t obey.

For XLD, all of these results are contingent on the assumption that the absorption through a line in the sample may be considered to be the line integral of the local absorption along the line. The local absorption must be independent of the absorption anywhere else, a condition which requires that the polarization must be preserved as the beam propagates through the sample, something that does not happen if dichroism is strong. This problem does not occur in MCD because polarization is preserved during propagation.