An algorithm for constrained one-step inversion of spectral CT data

Abstract

We develop a primal-dual algorithm that allows for one-step inversion of spectral CT transmission photon counts data to a basis map decomposition. The algorithm allows for image constraints to be enforced on the basis maps during the inversion. The derivation of the algorithm makes use of a local upper bounding quadratic approximation to generate descent steps for non-convex spectral CT data discrepancy terms, combined with a new convex-concave optimization algorithm. Convergence of the algorithm is demonstrated on simulated spectral CT data. Simulations with noise and anthropomorphic phantoms show examples of how to employ the constrained one-step algorithm for spectral CT data.

Fig. 1

Normalized spectrum of a typical X-ray source for CT operating at a potential of 120kV.

Fig. 2

Bone and brain maps derived from the FORBILD head phantom. Both images are shown in the gray scale window [0.9, 1.1].

Fig. 3

Convergence metrics for LSQ-TV and TPL-TV and for different values of λ with ideal, noiseless data. First, second, third, and fourth rows show the conditional primal-dual (cPD) gap, data discrepancy objective function, difference between the TV of estimated bone map and that of the phantom bone map, and same for the brain map TV. Note that the expressions for the gap and data discrepancy are different for TPL and LSQ; thus those quantities are not directly comparable.

Fig. 4

Convergence of the material map estimates to the phantom material maps for LSQ-TV and TPL-TV and for different values of λ with ideal, noiseless data.

Fig. 5

Difference between estimated brain and bone maps after 5,000 iterations and the corresponding phantom map shown in a 1% gray scale window [-0.01, 0.01] for TPL-TV and a 0.1% window [-0.001, 0.001] for LSQ-TV and different values of λ with ideal, noiseless data. The difference images are displayed in a region of interest around the sinus bones.

Fig. 6

Reconstructed bone map by use of TPL-TV from simulated noisy projection spectral CT transmission data. The material map TV constraints are varied according to fractions of the corresponding phantom material map TV.

Fig. 7

Reconstructed brain map by use of TPL-TV from simulated noisy projection spectral CT transmission data. The material map TV constraints are varied according to fractions of the corresponding phantom material map TV.

Fig. 8

Reconstructed bone map by use of LSQ-TV from simulated noisy projection spectral CT transmission data. The material map TV constraints are varied according to fractions of the corresponding phantom material map TV.

Fig. 9

Reconstructed brain map by use of LSQ-TV from simulated noisy projection spectral CT transmission data. The material map TV constraints are varied according to fractions of the corresponding phantom material map TV.

Fig. 10

Same as Fig. 3 except that only one value of λ is shown and the results are for noisy data and the TV constraints for the bone and brain maps are set to 1.1 × TV_bone and 1.1 × TV_brain, respectively. The TV constraint settings correspond to the center images in Figs. 6-9.

Fig. 11

Convergence of the material map estimates to the phantom material maps for LSQ-TV and TPL-TV and for noisy data with two different settings of the TV constraints. The TV factor applies to both the bone and brain maps, so that a TV factors of 1.1 and 1.2 correspond to the center and bottom, left images of Figs. 6-9.

Fig. 12

(Left) Chest phantom displayed at 70 KeV in a gray scale window of [0, 0.5] cm⁻¹. (Right) Reconstruction by use of unregularized TPL. The estimated material maps are combined to form the shown monochromatic image estimate at 70 KeV (gray scale is also [0, 1.0] cm⁻¹). For reference the TV values of the phantom and unconstrained reconstructed image are 2,587 and 7,686, respectively.

Fig. 13

Estimated monochromatic images by use of TPL-monoTV. The left column shows the complete image in a gray scale window of [0, 0.5] cm⁻¹. The right column magnifies a region of interest (ROI) in the right lung, and the gray scale is narrowed to [0, 0.1] cm⁻¹ in order to see the soft tissue detail. The top set of images correspond to the phantom. The location of the ROI is indicated in the left phantom image inset by use of the narrow [0, 0.1] cm⁻¹ gray scale. The second, third, and fourth rows correspond to images obtained by different TV constraints of the monoenergetic image at 70 KeV.

Fig. 14

Basis material maps: water (left), bone (middle), and Gadolinium contrast agent solution (right), corresponding to the monoenergetic image with TV of 1000 shown in Fig. 13. The basis material maps are shown in a gray scale window of [-0.2, 1.2]. The resulting reconstructed material maps agree well with the phantom maps in terms of structure, but interestingly the reconstructed maps show a larger noise level than the corresponding monochromatic image in Fig. 13, which is a linear combination of the shown material maps.

Fig. 15

Schematic illustrating the solution of ${max}_{g_m^{\prime }}\{g_m^{┬}{g}_m^{\prime }-\delta ({\Vert g_m^{\prime }\Vert }_1\le {\gamma }_m)\}$. The input vector g_m and the maximizing vector $g_m^{\prime }$ are indicated on a 2D schematic, but the argument applies for the full N_k-D space of g_m. Because $g_m^{\prime }$ is a vector of magnitudes, each component is non-negative $g_{m,k}^{\prime }\ge 0$. The indicator function confines $g_m^{\prime }$ below the line (hyper-plane), ${\sum }_k{g}_{m,k}^{\prime }={\gamma }_m$. The combination of these constraints confines $g_m^{\prime }$ to the schematic, shaded triangle. The maximizer $g_m^{\prime }$ is the vector that maximizes the dot product, $g_m^{┬}{g}_m^{\prime }$ (or equivalently the projection $g_m^{\prime }$ of onto g_m as indicated by the dashed line from the head of $g_m^{\prime }$ to the arrow indicating g_m). Maximization of this dot product is achieved by choosing $g_m^{\prime }={\gamma }_m{\widehat{e}}_{k-max}$ such that it is aligned along the unit vector corresponding to the largest component of g_m. The largest component of g_m is also known as the “infinity-norm”, ‖g_m‖∞. Thus we have (γ_mê_k−max)^┬g_m = γ_m ‖g_m‖∞.