Invertibility of Multi-Energy X-ray Transform

The ME X-ray transform, as defined in Equation (LABEL:eq:mean_mapping), has physical meaning when $\bm{A}$ is in the positive subspace of ${\bm{R}}^{N}$, denoted as ${\bm{P}}^{N}$, where $A_{i}\geq 0$ for all $i$. However, the transform is mathematically valid over the domain ${\bm{R}}^{N}$. In order to construct a simply connected $M_{2}$, we will expand the domain of the mapping from ${\bm{P}}^{N}$.

Let $U_{0}={\bm{P}}^{N}$ and $V_{0}$ be the image of $U_{0}$ under the mapping $\mathcal{X}$. In $V_{0}$, the mean photon count $\bm{I}$ is bounded by $0<I_{m}\leq I_{0}d_{m}$. Furthermore, as $U_{0}$ is path connected and $\mathcal{X}$ is a continuous mapping, $V_{0}$ is path connected [26](page 150). From every point $\bm{I}_{i}$ in $V_{0}$, draw a straight line to the maximum-count point $I_{0}\bm{d}$ and define this line with end points as $V_{i}$. Every point in $V_{i}$ is bounded by $I_{im}\leq I_{m}\leq I_{0}d_{m}$. Define $M_{2}$ as the union of all $V_{i}$. As $V_{0}$ and $V_{i}$ are all path connected, $M_{2}$ is path connected. The space $M_{2}$ is simply connected if every closed curve in $M_{2}$ can be contracted to a point [25]. Define a closed curve $\phi(s):[0,1]\rightarrow M_{2}$. We can contract $\phi(s)$ to the maximum-count point $I_{0}\bm{d}$ through the following continuous function $H$: $[0,1]\times[0,1]\rightarrow M_{2}$,

As $M_{2}$ is path connected, the closed curve $\phi(s)$ can be contracted to any points in $M_{2}$ [26] (page 332). Therefore, $M_{2}$ is simply connected.

We assume that the Jacobian vanishes nowhere in ${\bm{R}}^{N}$. In other words, the mapping $\mathcal{X}$ is a local homeomorphism, which means that every point of $\bm{A}\in{\bm{R}}^{N}$ has a neighborhood that is homeomorphic to an open subset in the range. For every straight line $V_{i}(\bm{I})$, the corresponding preimage $u_{i}(\bm{A})$ in the domain ${\bm{R}}^{N}$ can be constructed by successive local inverses $\mathcal{X}^{-1}(\bm{I})$. The maximum-count point $I_{0}\bm{d}$ corresponds to only one point in $U_{0}$, and this point is the origin of the coefficient space. Therefore, the origin is one end point of all preimages. The point $\bm{I}_{i}$ may have multiple local inverses, which we can index with subscript $j$. The $j^{th}$ local inverse introduces an inverse curve $u_{ij}$, where the $j^{th}$ local inverse is the second end point of the corresponding preimage $u_{ij}$. Each $u_{ij}$ is connected and connected to $U_{0}$. We define the union of all $u_{ij}(\bm{A})$ as $U_{i}$. Furthermore, we define $M_{1}$ as the union of all $U_{i}$. $M_{1}$ is connected and a superset of ${\bm{P}}^{N}$. The expansion of the range and domain for $N=2$ is illustrated in Figure 3.

The interior points of $M_{1}$ and $M_{2}$ (excluding the boundaries) are both smooth N-dimensional manifolds, because they are open subsets of ${\bm{R}}^{N}$ [27] (page 19). We will limit our proof to the interior of $M_{1}$ and $M_{2}$ and discuss the boundary points in Section VI.

We define the support of the weighting functions $p_{m}(E)$ as $\Omega_{m}$. Using the Schwarz inequality we have

with equality if and only if $p_{m}(E)\propto\exp\left[-{\bm{A}}\cdot{\bm{f}}\left(E\right)\right]$. In many occasions, the equality condition is not attainable. For example, when the three basis functions given in Equation (2) are used, $\exp\left[-{\bm{A}}\cdot{\bm{f}}\left(E\right)\right]$ is not proportional to the $p_{m}(E)$ of the detectors illustrated in FIG. 2(b)-(d). If we define

then we have

Assume that the length $\left|\Omega_{m}\right|$ of each support set is finite. Replacing the integrand with its maximum possible value gives

Therefore we have another set of inequalities

Now we may choose an energy $E_{m}$ such that

The right-hand side satisfies $\ln\left(\sqrt{\left|\Omega_{m}\right|/\gamma_{m}}\right)\geq\ln(d_{m}I_{0}/I_{m})\geq 0$ through the Schwarz inequality as follows:

Each of the inequalities in Equation (16) forces the vector ${\bm{A}}$ to be on the same side of the corresponding hyperplane as the origin.

Therefore, for given mean photon count ${\bm{I}}$, where $d_{m}I_{0}\geq I_{m}>0$, the inequalities defined by Equation (9) and (16) force ${\bm{A}}$ to be in a bounded set defined by the first set of hyperplanes and the second set of hyperplanes. A typical picture of this scenario for $N=2$ is shown in Figure 4(a). Note that for a physical measurement, the corresponding coefficients $\bm{A}$ are further bounded by the coordinate planes. Here we focus on demonstrating that $\bm{A}$ is bounded for a given $\bm{I}$ even without the positivity constraints on $\bm{A}$.

Now we show that the mapping $\mathcal{X}:M_{1}\rightarrow M_{2}$ is a proper mapping. If we have a compact set $C$ in the data space $M_{2}$, where all of the data vectors are located, then there are maximum and minimum values for each $I_{m}$ over all ${\bm{I}}$ in $C$. The maximum value for $I_{m}$ determines the hyperplane ${\bm{A}}\cdot{\bm{n}}_{m}=\ln[(d_{m}I_{0})/I_{m}]$ that is closet to the origin. The minimum value for $I_{m}$ determines the hyperplane ${\bm{A}}\cdot{\bm{f}}\left(E_{m}\right)=\ln\left(\sqrt{\left|\Omega_{m}\right|/\gamma_{m}}\right)$ that is furthest away from the origin. Therefore, the set of ${\bm{A}}$ that are mapped into $C$ are contained in a region bounded by these two sets of hyperplanes. This bounded region together with its boundary form a closed and bounded set in $\mathbb{R}^{N}$, hence a compact set. As the map ${\mathcal{X}}\left({\bm{A}}\right)$ is continuous, the set of ${\bm{A}}$ that are mapped into the closed set $C$ is closed. The set of ${\bm{A}}$ that are mapped into $C$ is a closed subset of a compact set. This set is therefore also compact. As a result, the mapping ${\mathcal{X}}({\bm{A}})$ is proper.

The integrand in Equation (20) has several interesting symmetry properties. The factor $\prod_{m=1}^{M}e^{-{\bm{A}}\cdot{\bm{f}}(E_{m})}=e^{-{\bm{A}}\cdot{\sum_{m=1}^{M}\bm{f}}(E_{m})}$ has mirror symmetry about all hyperplanes $E_{i}=E_{j}$ for $i,j\in\{1,2,..,M\}$. The other factor, $\det[F(\bm{E})]$, has sign-switching mirror symmetry about the same hyperplanes, which can be described mathematically as:

A sign-switching mirror symmetry means that, when we switch the positions of two coordinates (odd permutation), the sign of the function changes but the absolute value of the function is preserved. For example, with two coordinates $E_{1}$ and $E_{2}$, $\det[F({E_{1},E_{2}})]=f_{1}(E_{1})f_{2}(E_{2})-f_{2}(E_{1})f_{1}(E_{2})=-\det[F({E_{2},E_{1}})]$. To illustrate the sign-switching mirror symmetry, we plotted $\det[F(E_{1},E_{2})]$ for the case when $f_{1}(E)$ and $f_{2}(E)$ are both Gaussian functions in Figure 4(b).

Now we can divide the space occupied by $\Omega$ into $M!$ subspaces with hyperplanes $E_{i}=E_{j}$ for $i,j\in\{1,2,..,M\}$. One of the subspace has property $E_{1}<E_{2}<...<E_{M}$ and we define this subspace as $\Omega_{1,2...M}$. For every point $(E_{1},E_{2},...,E_{M})$ in the subspace $\Omega_{1,2...M}$, there is a corresponding point $(E_{\sigma_{1}},E_{\sigma_{2}},...,E_{\sigma_{M}})$ in each of the remaining subspaces. Applying the sign-switching mirror symmetry of the determinant, we can further simplify the Jacobian to:

The symbol $M$ emphasizes that the integration is over all spectral measurements. The three factors of the integrand, $\prod_{m=1}^{M}e^{-{\bm{A}}\cdot{\bm{f}}(E_{m})}$, $\det[F({\bm{E}})]$ and $\det[P(\bm{E})]$ depend on the total attenuation, the basis functions and the energy-weighting functions, respectively.

When the Jacobian is non-vanishing everywhere in ${\bm{R}}^{N}$, we can construct the domain $M_{1}$ and the range $M_{2}$ and prove that the mapping $\mathcal{X}:M_{1}\rightarrow M_{2}$ is globally invertible. As a result, the ME X-ray transform defined on the domain $P^{N}$, which is a subset of $M_{1}$, is globally invertible. Therefore, we have proved the equivalence of global invertibility with local invertibility for an ME X-ray transform. This equivalence holds when K-edge basis functions are considered.

When the photoelectric/Compton/Rayleigh basis functions are used, the basis-function determinant $\det[F({\bm{E}})]$ is always negative in the subspace $\Omega_{123}$. This set of three basis functions is sufficient to describe an object when the materials of interest have no K-edges in the energy range used for imaging, e.g. soft tissue and bone. The proof of $\det[F({\bm{E}})]<0$ for all $\bm{E}$ that satisfy $E_{1}<E_{2}<E_{3}$ is provided in Appendix A. When K-edge materials are considered, the values of $\det[F({\bm{E}})]$ can be calculated numerically and the positive regions of $\det[F({\bm{E}})]$ can be avoided by adjusting the detector sensitivity or source spectrum in $p_{m}(E)$.

Now we apply the sufficient condition for invertibility to the DE scenario that has non-unique solutions discussed by Levine [7]. For DE X-ray imaging, a set of two basis functions, photoelectric and Compton, can be used. The basis-function determinant is always positive in the subspace $S_{12}$. Levine assumed the same detector response for the two measurements, hence the sign of $\det[P({\bm{E}})]$ is the same with the sign of $\det[S({\bm{E}})]$, which is the source-spectra determinant. Similar to the definition in Equation (22), the matrix $S({\bm{E}})$ can be written as

where $(E_{1},E_{2})$ can be any combinations of two energies. With these, $\det[S({\bm{E}})]=S_{1}(E_{1})S_{2}(E_{2})-S_{1}(E_{2})S_{2}(E_{1})$. Levine assumed that both source spectra $S_{1}(E)$ and $S_{2}(E)$ are not zero only at three energy points (30, 60, 100) keV. Hence, $\det[S({\bm{E}})]$ is not zero only when $(E_{1},E_{2})$ are combinations of the set $\{30,60,100\}$ keV. Within the subspace $S_{12}$, where $E_{1}$ is always less than $E_{2}$, $\det[S({\bm{E}})]$ is non-vanishing only at three points $(E_{1},E_{2})=(30,60)$, $(30,100)$, and $(60,100)$ keV. Given the two source spectra as $S_{1}(E)=(1,1,1)$ and $S_{2}(E)=(0.93,1.71,0.30)$ at $E=(30,60,100)$ keV, respectively. The values of $\det[S({\bm{E}})]$ at $(E_{1},E_{2})=(30,60)$, $(30,100)$, and $(60,100)$ keV are 0.78, -0.63 and -1.41, respectively. Therefore, $\det[S({\bm{E}})]$ does not have constant sign over the subspace $S_{12...N}$, hence the invertibility of the DE X-ray transform for the proposed scenario is not guaranteed. Note that this analysis only shows that the existence of non-unique solutions is possible, it does not prove their existence.

We then consider the sign of $\det[P({\bm{E}})]$ for the four types of detectors illustrated in Figure 2. Assuming the source spectrum $S(E)$ is same for different $m$, the weighting-function determinant $\det[P({\bm{E}})]$ has the same sign with the detector-response determinant $\det[D({\bm{E}})]$, where the $(i,j)$ element of matrix $D({\bm{E}})$ is the detector response of the $i^{th}$ measurement at $E_{j}$ denoted as $D_{i}(E_{j})$. The sign of $\det[D({\bm{E}})]$ for the four types of detectors are studied in Appendix A and the main results are presented as follows.

In conclusion, the mapping $\mathcal{X}:{\bm{A}}\rightarrow{\bm{I}}$ is globally invertible for these four types of detectors when measuring attenuation profiles without K-edges. For arbitrary detectors or systems with varying source spectra, the values of $\det[P({\bm{E}})]$ can be calculated numerically. One can alway maintain $\det[F({\bm{E}})]\det[P({\bm{E}})]\leq 0$ for any $E$ in $S_{12...N}$ by adjusting the energy-response function or the bin boundaries of the detectors.

We simulated a single X-ray path with different lengths of water as the attenuating media. The energy-bin boundaries of the detector were set at [30, 75, 100, 160] keV. The length of water ranged from 1 cm to 28 cm. For each length, 1000 sets of noisy data were generated, and the AM coefficients were estimated for each set of data. We calculated the mean and the variance of the estimated coefficients and compared the estimation uncertainty to the Crámer-Rao lower bound. The results are presented in Figure 5.

To check if the algorithm works for materials that are very different from water, we changed the material in the X-ray path to iron and the detector bin boundaries to [30, 110, 140, 160] keV. The bin boundaries were changed so that there are sufficient number of photons collected in all three bins. The length of iron ranged from 0.1 cm to 3 cm. The mean and variance of the 1000 repeated estimations at each length are presented in Figure 6.

In both scenarios, the mean of the estimates (black line) matches well with the true coefficient (red circle). The slight deviation between the true coefficient and the mean estimation at high attenuation region can be attributed to sampling error in Monte-Carlo simulation. The standard deviation of the estimates (purple area) is almost perfectly aligned with the Crámer-Rao lower bound. These results demonstrate that our estimation algorithm is unbiased and efficient. Furthermore, the AM coefficient corresponding to the Rayleigh scattering is estimable and the uncertainty in $\hat{A}_{rs}$ is comparable to the uncertainty in $\hat{A}_{pe}$, which corresponds to the photoelectric effect. However, the $\hat{A}_{pe}$ and $\hat{A}_{rs}$ are anti-correlated, which is demonstrated by a negative value of the element $[\mathrm{FIM}^{-1}]_{13}$. This correlation is probably due to the fact that the basis functions $f_{pe}(E)$ and $f_{rs}(E)$ have similar shapes. This correlation may partially explain why $\hat{A}_{pe}$ and $\hat{A}_{rs}$ have more variance than $\hat{A}_{cs}$, as shown in Figure 5. As a result of the correlation, the estimated total attenuation $\hat{\tau}(E)=\hat{A}_{pe}f_{pe}(E)+\hat{A}_{cs}f_{cs}(E)+\hat{A}_{rs}f_{rs}(E)$, which is the ultimate physical quantity we are interest in, is not very noisy, as will be shown in the reconstruction results of the phantom.

We simulate a two-dimensional fan-beam CT system (62^∘ fan-beam angle) with 360 views and 245 detectors. The field of view is 256x256 pixels with a pixel pitch of 1—1.5 mm. The same source, detector, and material database described in the previous section are used. X-ray attenuation is simulated according to Beer’s law, while scattering and detector imperfections are not considered. The detector energy bin boundaries are [30, 75, 100, 160] keV. The source photon budget for each beam path is $10^{7}$. The AM coefficients $\bm{A}$ are estimated for each beam path and the object represented by $\bm{a}(\bm{R})$ is reconstructed from $\hat{\bm{A}}$ using a filtered-back projection (FBP) algorithm.

The first phantom reconstructed is a circular water phantom of diameter about 30 cm with pixel-pitch 0.15 cm. This phantom and its reconstruction are plotted at $E=75.5$ keV in Figure 7. The reconstruction matches well with the object. As shown in the center-line plot in the right panel of Figure 7, the reconstruction has no cupping artifacts, which is typically associated with beam hardening.

A second phantom is a multi-material resolution phantom design inspired by Gong et.al. [33]. The length of the phantom is 20 cm with pixel pitch around 0.1 cm. The phantom, as shown in Figure 8 (left), is a Delrin block with 25 circular inserts in five rows and five columns. In each row, the inserts are made from the same material; in each column, the inserts have the same diameter. From top to bottom, the five materials for the inserts are water, polyvinyl chloride, cast magnesium, acrylic and methanol. The diameter of the inserts are from 0.6 to 1.8 cm with 0.3 cm step. The reconstruction of the resolution phantom are plotted at $E=75.5$ keV in Figure 8. In the plots, different shades of grey represent different materials. The reconstruction results show that the ML estimation algorithm works for a broad range of AM coefficients

Each reconstruction, which calls the ML estimation algorithm 88,200 times, takes approximately 130 seconds on a desktop with a quad-core central processing unit (CPU). The reconstruction can be further sped up using a graphic processing unit (GPU).