An oblique projection modification technique (OPMT) for fast multispectral CT reconstruction

Suppose that after $n$ iterations, we have the estimation $\left({{{\boldsymbol{f}}}_{1}}^{(n)},{{{\boldsymbol{f}}}_{2}}^{(n)},\cdots ,{{{\boldsymbol{f}}}_{N}}^{(n)}\right)$ of the basis material images, and make the first-order Taylor expansion of formula (6) at $\left({{{\boldsymbol{f}}}_{1}}^{(n)},{{{\boldsymbol{f}}}_{2}}^{(n)},\cdots ,{{{\boldsymbol{f}}}_{{\rm{N}}}}^{(n)}\right)$ to obtain the linearized polychromatic projection equation

$$\begin{eqnarray}&&\displaystyle \sum _{i=1}^{N}\displaystyle \frac{{{\rm{\Theta }}}_{i,k,l}^{(n)}}{{q}_{k,l}^{(n)}}{R}_{l}({{\boldsymbol{f}}}_{i}-{{{\boldsymbol{f}}}_{i}}^{(n)})={p}_{k,l}-{p}_{k,l}^{(n)},l\in {\xi }_{k},k=1,2,\ldots ,{\rm{N}}\end{eqnarray} \tag{ 7 }$$

where,

$$\begin{eqnarray}\begin{array}{lll}{p}_{k,l}^{(n)} & = & -\,\mathrm{ln}\,\displaystyle \sum _{m=1}^{{{\rm{M}}}_{k}}{S}_{k,m}\delta {e}^{-\displaystyle {\sum }_{i=1}^{{\rm{N}}}{\theta }_{i,m}{R}_{l}{{{\boldsymbol{f}}}_{i}}^{(n)}}\\ {q}_{k,l}^{(n)} & = & \displaystyle \sum _{m=1}^{{{\rm{M}}}_{k}}{S}_{k,m}\delta {e}^{-\displaystyle {\sum }_{i=1}^{{\rm{N}}}{\theta }_{i,m}{R}_{l}{{{\boldsymbol{f}}}_{i}}^{(n)}}\\ {{\rm{\Theta }}}_{i,k,l}^{(n)} & = & \displaystyle \sum _{m=1}^{{{\rm{M}}}_{k}}{S}_{k,m}\delta {\theta }_{i,m}{e}^{-\displaystyle {\sum }_{i=1}^{{\rm{N}}}{\theta }_{i,m}{R}_{l}{{{\boldsymbol{f}}}_{i}}^{(n)}}.\end{array}\end{eqnarray} \tag{ 8 }$$

In the case of geometric inconsistency, suppose that we only know one polychromatic projection corresponding to a specific spectrum along the given x-ray path. We assume that the known polychromatic projection data ${p}_{1,l}$ is obtained from the spectrum ${S}_{1}(E)$ along the x-ray path ${\xi }_{1}.$ And the values of ${p}_{k,l},k=2,3,\cdots ,{\rm{N}}$ are unknown from the spectrum ${S}_{k}(E),$ $k=2,3,\cdots ,{\rm{N}}$ along the same path. On the current x-ray path, the projection equation (7) is transformed into the following linear equations with N variables.

$$\begin{eqnarray}&&{a}_{11}{x}_{1}+{a}_{12}{x}_{2}+\ldots +{a}_{1{\rm{N}}}{x}_{{\rm{N}}}={b}_{1}^{{\rm{known}}},l\in {\xi }_{1}\end{eqnarray} \tag{ 9.1 }$$

$$\begin{eqnarray*}&&\vdots \end{eqnarray*}$$

$$\begin{eqnarray}&&{a}_{n1}{x}_{1}+{a}_{n2}{x}_{2}+\ldots +{a}_{n{\rm{N}}}{x}_{{\rm{N}}}={b}_{n}^{{\rm{unknown}}},l\in {\xi }_{1}\end{eqnarray} \tag{ 9.n }$$

$$\begin{eqnarray*}&&\vdots \end{eqnarray*}$$

$$\begin{eqnarray}&&{a}_{{\rm{N}}1}{x}_{1}+{a}_{{\rm{N}}2}{x}_{2}+\ldots +{a}_{{\rm{NN}}}{x}_{{\rm{N}}}={b}_{{\rm{N}}}^{{\rm{unknown}}},l\in {\xi }_{1}\end{eqnarray} \tag{ 9.N }$$

where,

$$\begin{eqnarray}\begin{array}{ccc}{a}_{ij} & = & \displaystyle \frac{{{\rm{\Theta }}}_{j,i,l}^{(n)}}{{q}_{i,l}^{(n)}}\\ {x}_{n} & = & {R}_{l}{{\boldsymbol{f}}}_{n}\\ {b}_{1}^{{\rm{known}}} & = & {p}_{1,l}-{p}_{1,l}^{(n)}+\displaystyle \sum _{k=1}^{N}\displaystyle \frac{{{\rm{\Theta }}}_{k,1,l}^{(n)}}{{q}_{1,l}^{(n)}}{R}_{l}{{{\boldsymbol{f}}}_{k}}^{(n)}\\ {b}_{n}^{{\rm{unknown}}} & = & {p}_{n,l}^{{\rm{unknown}}}-{p}_{n,l}^{(n)}+\displaystyle \sum _{k=1}^{N}\displaystyle \frac{{{\rm{\Theta }}}_{k,n,l}^{(n)}}{{q}_{n,l}^{(n)}}{R}_{l}{{{\boldsymbol{f}}}_{k}}^{(n)},n=2,3,\cdots ,{\rm{N}}\end{array}\end{eqnarray} \tag{ 10 }$$

where ${b}_{1}^{{\rm{known}}}$ is the value calculated from the measured polychromatic projection data ${p}_{1,l}$ on the current x-ray path. Here, ${p}_{n,l}^{{\rm{unknown}}},$ $n=2,3,\cdots ,N$ are the unknown polychromatic projection data along the current x-ray path, ${b}_{n}^{{\rm{unknown}}},$ $n=2,3,\cdots ,{\rm{N}}$ are unknown items, while ${a}_{ij},\left(i=1\cdots {\rm{N}}\right.,$ $\left.j=1\cdots {\rm{N}}\right)$ are known.

For better statement of the main idea of the proposed method, let us focus on the iterative process of the E-ART algorithm before getting back to this. As mentioned before about the E-ART algorithm, the polychromatic projection equation corresponding to only one spectrum along a certain x-ray path is used each time. Here, we mark this linearized polychromatic projection equation as ${{\rm{H}}}_{1},$ formed as equation (9.1). As shown in figure 1(a), the basis material image is updated from the current iteration point ${R}_{l}{{\boldsymbol{f}}}^{(n)}=\left({R}_{l}{{{\boldsymbol{f}}}_{1}}^{(n)},{R}_{l}{{{\boldsymbol{f}}}_{2}}^{(n)},\cdots ,{R}_{l}{{{\boldsymbol{f}}}_{N}}^{(n)}\right)$ to ${{\rm{H}}}_{1}$ as an orthogonal projection, and the updated iteration point is ${R}_{l}{{\boldsymbol{f}}}^{(n+1)}=\left({R}_{l}{{{\boldsymbol{f}}}_{1}}^{(n+1)},{R}_{l}{{{\boldsymbol{f}}}_{2}}^{(n+1)},\cdots ,{R}_{l}{{{\boldsymbol{f}}}_{N}}^{(n+1)}\right).$ Aimed at the MCT reconstruction problem, on one x-ray path, the condition number of projection equations corresponding to multiple spectra is high; thus this kind of update mode slows down the convergence speed.

Zoom In Zoom Out Reset image size

Returning to the subject, it can be observed from equations $(9.i),i=2,3,\cdots ,{\rm{N}}$ that the values of the right-hand side of equations $(9.i)$ are unknown, but the coefficients ${a}_{ij}$ are known; that is, the normal direction of the hyperplane represented by each equation in $(9.i)$ is determined. Consequently, we can calculate the normal direction of the hyperplane formed by the normal directions of the ${\rm{N}}-1$ equations in N-dimensional space (note that the normal direction of the hyperplane is H₂). H₂ contains the directional information of the linearized equation set corresponding to the unknown spectra along the current x-ray path.

Based on this observation, we can use H₁ and H₂ to construct the oblique projection modification direction of the iteration point ${R}_{l}{{\boldsymbol{f}}}^{(n)}$ to H₁. Here, H₁ represents the equation of polychromatic projection calculated by the known spectrum along the current x-ray path, and H₂ shows the directional information of the corresponding projection from unknown spectra along the same path. As shown in figure 1(b), $dir1$ represents the direction of H₂, and $dir2$ represents the normal direction of the hyperplane represented by H₁. In fact, $dir2$ is also the direction of the E-ART algorithm updating estimations as the orthogonal projection from the current iteration point ${R}_{l}{{\boldsymbol{f}}}^{(n)}$ to H₁. The direction of oblique projection modification from the current iteration point ${R}_{l}{{\boldsymbol{f}}}^{(n)}$ to H₁ is obtained by the weighted summation of direction $dir1$ and direction $dir2.$ Then, a one-dimensional search method is used to calculate the new iteration point. Because the oblique projection modification direction comprehensively refers to the directional information of the polychromatic projection equation corresponding to the known spectrum along the current x-ray path and the unknown spectra on the same path, the convergence speed of the proposed algorithm is significantly improved compared with the ASD-NC-POCS and E-ART algorithms, in which only the projection equation of the known spectrum on the current x-ray path is considered.

This section describes how to select direction $dir1$ and direction $dir2$ in the proposed algorithm, and gives the formula of the new iteration points.

The direction of orthogonal projection from the current iteration point $\left({R}_{l}{{{\boldsymbol{f}}}_{1}}^{(n)},{R}_{l}{{{\boldsymbol{f}}}_{2}}^{(n)},\cdots ,{R}_{l}{{{\boldsymbol{f}}}_{N}}^{(n)}\right)$ to the hyperplane (9.1) can be calculated by (11).

$$\begin{eqnarray}\begin{array}{ccc}dir21 & = & \displaystyle \frac{{b}_{1}^{{\rm{known}}}-\displaystyle \sum _{k=1}^{{\rm{N}}}{a}_{1,k}{x}_{k}}{\displaystyle \sum _{k=1}^{N}{\left({a}_{1,k}\right)}^{2}}\left({a}_{1,1},{a}_{1,2},\cdots ,{a}_{1,n},\cdots ,{a}_{1,{\rm{N}}}\right)\\ & = & \displaystyle \frac{{p}_{1,l}-{p}_{1,l}^{(n)}}{\displaystyle \sum _{k=1}^{{\rm{N}}}{\left(\tfrac{{{\rm{\Theta }}}_{k,1,l}^{(n)}}{{q}_{1,l}^{(n)}}\right)}^{2}}\left(\displaystyle \frac{{{\rm{\Theta }}}_{1,1,l}^{(n)}}{{q}_{1,l}^{(n)}},\displaystyle \frac{{{\rm{\Theta }}}_{2,1,l}^{(n)}}{{q}_{1,l}^{(n)}},\cdots ,\displaystyle \frac{{{\rm{\Theta }}}_{n,1,l}^{(n)}}{{q}_{1,l}^{(n)}},\cdots ,\displaystyle \frac{{{\rm{\Theta }}}_{{\rm{N}},1,l}^{(n)}}{{q}_{1,l}^{(n)}}\right).\end{array}\end{eqnarray} \tag{ 11 }$$

Take the unit vector of (11) as the normal direction of the hyperplane represented by ${{\rm{H}}}_{1},$ i.e.

$$\begin{eqnarray}&&dir2=\displaystyle \frac{dir21}{\parallel dir21\parallel }\end{eqnarray} \tag{ 12 }$$

where $\parallel \cdot \parallel$ is a vector modulus operator.

According to $(9.i),i=2,3,\cdots ,{\rm{N}},$ vectors $({a}_{i1},{a}_{i2},\cdots ,{a}_{i{\rm{N}}}),$ $i=2,3,\cdots ,{\rm{N}}$ are the normal directions of the hyperplanes represented by $(9.i),i=2,3,\cdots ,{\rm{N}}$, respectively. The normal direction of the hyperplane formed by these ${\rm{N}}-1$ vectors in N-dimensional space is expressed as follows:

$$\begin{eqnarray}{\rm{\det }}=\left|\begin{array}{rrrrr}{{\boldsymbol{e}}}_{1} & {{\boldsymbol{e}}}_{2} & {{\boldsymbol{e}}}_{3} & \cdots & {{\boldsymbol{e}}}_{{\rm{N}}}\\ {a}_{2,1} & {a}_{2,2} & {a}_{2,3} & \cdots & {a}_{2,{\rm{N}}}\\ \cdots & \cdots & \cdots & \cdots & \cdots \\ {a}_{n,1} & {a}_{n,2} & {a}_{n,3} & \cdots & {a}_{n,{\rm{N}}}\\ \cdots & \cdots & \cdots & \cdots & \cdots \\ {a}_{{\rm{N}},1} & {a}_{{\rm{N}},2} & {a}_{{\rm{N}},3} & \cdots & {a}_{N,N}\end{array}\right|=\displaystyle \sum _{i=1}^{{\rm{N}}}{D}_{i}{{\boldsymbol{e}}}_{i}\end{eqnarray} \tag{ 13 }$$

where ${{\boldsymbol{e}}}_{1},{{\boldsymbol{e}}}_{2},\cdots ,{{\boldsymbol{e}}}_{{\rm{N}}}$ is a set of standard orthogonal bases in N-dimensional space and $\left|\cdot \right|$ is a determinant operator.

${D}_{n}={(-1)}^{n+1}\left|\begin{array}{ccccccc}{a}_{2,1} & {a}_{2,2} & \cdots & {a}_{2,n-1} & {a}_{2,n+1} & \cdots & {a}_{2,{\rm{N}}}\\ {a}_{3,1} & {a}_{3,2} & \cdots & {a}_{3,n-1} & {a}_{3,n+1} & \cdots & {a}_{3,{\rm{N}}}\\ \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots \\ {a}_{{\rm{N}},1} & {a}_{{\rm{N}},2} & \cdots & {a}_{{\rm{N}},n-1} & {a}_{{\rm{N}},n+1} & \cdots & {a}_{N,N}\end{array}\right|.$ ${D}_{n}$ is the algebraic cofactor of element ${{\boldsymbol{e}}}_{n}$ of the determinant $\det .$ Based on the properties of the algebraic cofactor of the determinant, we can get (14)

$$\begin{eqnarray}&&\displaystyle \sum _{j=1}^{{\rm{N}}}{D}_{j}{a}_{i,j}=0,i=2,3,\cdots ,{\rm{N}}.\end{eqnarray} \tag{ 14 }$$

So, vector $\left({D}_{1},{D}_{2},\cdots ,{D}_{{\rm{N}}}\right)$ is orthogonal to $({a}_{i1},{a}_{i2},\cdots ,{a}_{i{\rm{N}}}),$ $i=2,3,\cdots ,{\rm{N}}.$ Then, the normal direction of the hyperplane formed by the normal directions of the ${\rm{N}}-1$ equations in N-dimensional space can be expressed as (15)

$$\begin{eqnarray}&&dir11=\left({D}_{1},{D}_{2},\cdots ,{D}_{{\rm{N}}}\right)\,or\,dir12=-dir11.\end{eqnarray} \tag{ 15 }$$

Take the direction with an acute angle with $dir2$ in the two directions of ${{\rm{H}}}_{2}$ represented by formula (15) as the direction $dir1$

$$\begin{eqnarray}dir1=\left\{\begin{array}{c}\begin{array}{cc}\displaystyle \frac{dir11}{\parallel dir11\parallel }, & if\left\langle dir11,dir2\right\rangle \gt 0\end{array}\\ \begin{array}{cc}\displaystyle \frac{dir12}{\parallel dir12\parallel }, & if\left\langle dir12,dir2\right\rangle \gt 0\end{array}\end{array}\right.\end{eqnarray} \tag{ 16 }$$

where $\left\langle ,\right\rangle$ is the inner product of two vectors.

Let

$$\begin{eqnarray}&&dir={\lambda }_{1}dir1+{\lambda }_{2}dir2\end{eqnarray} \tag{ 17 }$$

be the direction of oblique projection modification from the current iteration point $\left({R}_{l}{{{\boldsymbol{f}}}_{1}}^{(n)},{R}_{l}{{{\boldsymbol{f}}}_{2}}^{(n)},\cdots ,{R}_{l}{{{\boldsymbol{f}}}_{N}}^{(n)}\right)$ to the known polychromatic projection equation ${{\rm{H}}}_{1}$ along the current x-ray path, where ${\lambda }_{1}$ and ${\lambda }_{2}$ are the coefficients weighting the two directions. The updated line integrals in this direction are $\left({R}_{l}{{{\boldsymbol{f}}}_{1}}^{(n+1)},\cdots ,{R}_{l}{{{\boldsymbol{f}}}_{n}}^{(n+1)},\right.$ $\left.\cdots ,{R}_{l}{{{\boldsymbol{f}}}_{{\rm{N}}}}^{(n+1)}\right)$

$$\begin{eqnarray}&&\left(\begin{array}{c}{R}_{l}{{{\boldsymbol{f}}}_{1}}^{(n+1)}\\ \vdots \\ {R}_{l}{{{\boldsymbol{f}}}_{n}}^{(n+1)}\\ \vdots \\ {R}_{l}{{{\boldsymbol{f}}}_{{\rm{N}}}}^{(n+1)}\end{array}\right)=\left(\begin{array}{c}{R}_{l}{{{\boldsymbol{f}}}_{1}}^{(n)}\\ \vdots \\ {R}_{l}{{{\boldsymbol{f}}}_{n}}^{(n)}\\ \vdots \\ {R}_{l}{{{\boldsymbol{f}}}_{{\rm{N}}}}^{(n)}\end{array}\right)+\displaystyle \frac{\left({p}_{1,l}-{p}_{1,l}^{(n)}\right)* di{r}^{T}}{\left\langle \left(\tfrac{{{\rm{\Theta }}}_{1,1,l}^{(n)}}{{q}_{1,l}^{(n)}},\tfrac{{{\rm{\Theta }}}_{2,1,l}^{(n)}}{{q}_{1,l}^{(n)}},\cdots ,\tfrac{{{\rm{\Theta }}}_{n,1,l}^{(n)}}{{q}_{1,l}^{(n)}},\cdots ,\tfrac{{{\rm{\Theta }}}_{{\rm{N}},1,l}^{(n)}}{{q}_{1,l}^{(n)}}\right),dir\right\rangle }.\end{eqnarray} \tag{ 18 }$$

Then, we can get the updated basis material images by applying the linear reconstruction operator ${{R}_{l}}^{-1}$ to the line integrals, i.e.

$$\begin{eqnarray}&&\left(\begin{array}{c}{{{\boldsymbol{f}}}_{1}}^{(n+1)}\\ \vdots \\ {{{\boldsymbol{f}}}_{n}}^{(n+1)}\\ \vdots \\ {{{\boldsymbol{f}}}_{{\rm{N}}}}^{(n+1)}\end{array}\right)=\left(\begin{array}{c}{{{\boldsymbol{f}}}_{1}}^{(n)}\\ \vdots \\ {{{\boldsymbol{f}}}_{n}}^{(n)}\\ \vdots \\ {{{\boldsymbol{f}}}_{{\rm{N}}}}^{(n)}\end{array}\right)+{{R}_{l}}^{-1}\left(\displaystyle \frac{\left({p}_{1,l}-{p}_{1,l}^{(n)}\right)* di{r}^{T}}{\left\langle \left(\tfrac{{{\rm{\Theta }}}_{1,1,l}^{(n)}}{{q}_{1,l}^{(n)}},\tfrac{{{\rm{\Theta }}}_{2,1,l}^{(n)}}{{q}_{1,l}^{(n)}},\cdots ,\tfrac{{{\rm{\Theta }}}_{n,1,l}^{(n)}}{{q}_{1,l}^{(n)}},\cdots ,\tfrac{{{\rm{\Theta }}}_{{\rm{N}},1,l}^{(n)}}{{q}_{1,l}^{(n)}}\right),dir\right\rangle }\right).\end{eqnarray} \tag{ 19 }$$

In the case of geometric consistency, on the current x-ray path, the projection equation (7) is transformed into the following linear equations with N variables

$$\begin{eqnarray}&&{a}_{11}{x}_{1}+{a}_{12}{x}_{2}+\ldots +{a}_{1{\rm{N}}}{x}_{{\rm{N}}}={b}_{1}^{{\rm{known}}},l\in {\xi }_{1}\end{eqnarray} \tag{ 20.1 }$$

$$\begin{eqnarray*}&&\vdots \end{eqnarray*}$$

$$\begin{eqnarray}&&{a}_{n1}{x}_{1}+{a}_{n2}{x}_{2}+\ldots +{a}_{n{\rm{N}}}{x}_{{\rm{N}}}={b}_{n}^{{\rm{known}}},l\in {\xi }_{1}\end{eqnarray} \tag{ 20.n }$$

$$\begin{eqnarray*}&&\vdots \end{eqnarray*}$$

$$\begin{eqnarray}&&{a}_{{\rm{N}}1}{x}_{1}+{a}_{{\rm{N}}2}{x}_{2}+\ldots +{a}_{{\rm{NN}}}{x}_{{\rm{N}}}={b}_{{\rm{N}}}^{{\rm{known}}},l\in {\xi }_{1}\end{eqnarray} \tag{ 20.N }$$

where,

$$\begin{eqnarray}\begin{array}{ccc}{a}_{ij} & = & \displaystyle \frac{{{\rm{\Theta }}}_{j,i,l}^{(n)}}{{q}_{i,l}^{(n)}}\\ {x}_{n} & = & {R}_{l}{{\boldsymbol{f}}}_{n}\\ {b}_{n}^{{\rm{known}}} & = & {p}_{n,l}-{p}_{n,l}^{(n)}+\displaystyle \sum _{k=1}^{{\rm{N}}}\displaystyle \frac{{{\rm{\Theta }}}_{k,n,l}^{(n)}}{{q}_{n,l}^{(n)}}{R}_{l}{{{\boldsymbol{f}}}_{k}}^{(n)},n=1,2,3,\cdots ,{\rm{N}}\end{array}\end{eqnarray} \tag{ 21 }$$

where ${b}_{n}^{{\rm{known}}},n=1,2,3,\cdots ,{\rm{N}}$ are the values calculated from the measured polychromatic projection data ${p}_{n,l}$ on the current x-ray path. Here,${b}_{n}^{{\rm{known}}},n=1,2,3,\cdots ,{\rm{N}}$ and ${a}_{ij},\left(i=1\cdots {\rm{N}}\right.,$ $\left.j=1\cdots {\rm{N}}\right)$ are known items. In this paper, the projections of the basis materials are updated ray-by-ray; that is, only the polychromatic projection equation corresponding to one spectrum on the current ray path is used each time. Suppose (20.1) is used, which means ${b}_{n}^{{\rm{known}}},n=2,3,\cdots ,{\rm{N}}$ is unknown, and equations (20.1)∼(20.N) are equivalent to equations (9.1)∼(9.N). The method described in section 2.2.2 also applies to the case of consistent polychromatic projections.

We summarize the implementation steps of the iteration scheme as follows.

Step 1. Assign ${{{\boldsymbol{f}}}_{1}}^{(0)},\cdots ,{{{\boldsymbol{f}}}_{n}}^{(0)},\cdots ,{{{\boldsymbol{f}}}_{{\rm{N}}}}^{(0)}$ with some initial values.

Step 2. Suppose the estimations ${{{\boldsymbol{f}}}_{1}}^{(n)},\cdots ,{{{\boldsymbol{f}}}_{n}}^{(n)},\cdots ,{{{\boldsymbol{f}}}_{{\rm{N}}}}^{(n)}$ are calculated after $n$ iterations. For a given x-ray path along which the polychromatic projection is ${p}_{1,l},$ calculate ${b}_{1}^{{\rm{known}}}$ and ${a}_{ij},(i=1\cdots {\rm{N}},j=1\cdots {\rm{N}})$ by (8) and (10).

Step 3. First calculate $dir2$ by (11) and (12). Then calculate $dir1$ by (13), (15) and (16). Finally, calculate the direction of oblique projection modification ($dir$) by (17).

Step 4. Calculate the updated projections of basis material images ${R}_{l}{{{\boldsymbol{f}}}_{1}}^{(n+1)},\cdots ,{R}_{l}{{{\boldsymbol{f}}}_{n}}^{(n+1)},\cdots ,{R}_{l}{{{\boldsymbol{f}}}_{{\rm{N}}}}^{(n+1)}$ along the current x-ray path by (18). Then, get the updated basis material images ${{{\boldsymbol{f}}}_{1}}^{(n+1)},\cdots ,{{{\boldsymbol{f}}}_{n}}^{(n+1)},\cdots ,{{{\boldsymbol{f}}}_{{\rm{N}}}}^{(n+1)}$ by (19).

Step 5. Stop if the convergence criteria are satisfied, otherwise turn to step 2.

The parameters of the algorithm are ${\lambda }_{1}$ and ${\lambda }_{2}$ in equation (16), which are used to adjust the weights of $dir1$ and $dir2$ in the oblique projection modification direction $dir.$ By increasing the value of ${\lambda }_{1}/{\lambda }_{2}$ properly, the convergence speed of the proposed algorithm can be accelerated. The higher the condition number of the linear equation set in (9.1) to (9.N) is, the more ill-posed the MCT reconstruction problem becomes. Here, the value of ${\lambda }_{1}/{\lambda }_{2}$ should be reduced so that the noise will not be excessively amplified. In general, we choose ${\lambda }_{1}/{\lambda }_{2}\in [0.5,1.5]$ in the initial iteration, so that it can not only get high-quality basis material images, but also accelerate the convergence speed of the basis material images. With the continuous iterations, we make the parameter ${\lambda }_{1}$ decrease to 0, which can increase the accuracy of the proposed algorithm.

According to (7), we construct a system of two-variable nonlinear equations along a specific x-ray path

$$\begin{eqnarray}&&{p}_{1,l}=-\,\mathrm{ln}\,\displaystyle \sum _{m=1}^{{{\rm{M}}}_{1}}{S}_{1,m}\delta {e}^{-{\theta }_{1,m}{R}_{l}{{\boldsymbol{f}}}_{1}-{\theta }_{2,m}{R}_{l}{{\boldsymbol{f}}}_{2}},l\in {\xi }_{1}\end{eqnarray} \tag{ 22.1 }$$

$$\begin{eqnarray}&&{p}_{2,l}=-\,\mathrm{ln}\,\displaystyle \sum _{m=1}^{{{\rm{M}}}_{2}}{S}_{2,m}\delta {e}^{-{\theta }_{1,m}{R}_{l}{{\boldsymbol{f}}}_{1}-{\theta }_{2,m}{R}_{l}{{\boldsymbol{f}}}_{2}},l\in {\xi }_{1}\end{eqnarray} \tag{ 22.2 }$$

where ${S}_{1}(E)=({S}_{1,1},{S}_{1,2},\cdots ,{S}_{1,{{\rm{M}}}_{1}})$ and ${S}_{2}(E)=({S}_{2,1},{S}_{2,2},\cdots ,{S}_{2,{{\rm{M}}}_{2}})$ are different spectra. To show the iterative paths of different algorithms clearly, and by decreasing the dependency of two spectra as much as we possibly can, we divide the spectrum of 140 kV into two sections. Here, ${S}_{1}(E)$ takes 1–70 kV data from the 140 kV spectrum, and ${S}_{2}(E)$ takes 71–140 kV data from the 140 kV spectrum. Here, ${\theta }_{1}(E)=({\theta }_{1,1},{\theta }_{1,2},\cdots ,{\theta }_{1,{{\rm{M}}}_{k}})$ and ${\theta }_{2}(E)=({\theta }_{2,1},{\theta }_{2,2},\cdots ,{\theta }_{2,{{\rm{E}}}_{k}})$ are the mass attenuation coefficients of water and bone, respectively; ${p}_{1,l}$ and ${p}_{2,l}$ are polychromatic projection values obtained under ${S}_{1}(E)$ and ${S}_{2}(E),$ respectively, corresponding to the projection of water basis material ${R}_{l}{{\boldsymbol{f}}}_{1}=4$ and the projection of bone basis material ${R}_{l}{{\boldsymbol{f}}}_{2}=1.$ We assume that ${S}_{1}(E),$ ${S}_{2}(E),$ ${\theta }_{1}(E)$ and ${\theta }_{2}(E)$ are known and ${R}_{l}{{\boldsymbol{f}}}_{1},$ ${R}_{l}{{\boldsymbol{f}}}_{2}$ are to be evaluated; then equations (22.1) and (22.2) constitute a nonlinear system of equations regarding the unknown number $\left({R}_{l}{{\boldsymbol{f}}}_{1},{R}_{l}{{\boldsymbol{f}}}_{2}\right).$ In the experiment, we use the ASD-NC-POCS, E-ART and the proposed algorithm to solve the equations. Only one equation is used in each iteration, and the other polychromatic projection equation is assumed to be unknown. In the proposed algorithm, increasing the value of ${\lambda }_{1}/{\lambda }_{2}$ appropriately can accelerate the convergence speed of the reconstructed images. To show the iterative paths of the proposed algorithm clearly in the graph, the parameter is selected as ${\lambda }_{1}=1,{\lambda }_{2}=1$ in this experiment.

Since the iterative paths of the ASD-NC-POCS algorithm are similar to those of the E-ART algorithm, to show this directly in figure 3, we only draw the iterative paths and convergence curve of the E-ART algorithm and the proposed algorithm. Figure 3(a) shows the iterative paths of the E-ART algorithm and the proposed algorithm. The black solid line and black dotted line are polychromatic projection curves under ${S}_{1}(E)$ and ${S}_{2}(E)$, respectively. The blue line represents the iterative paths of the E-ART algorithm, and the red one represents the proposed algorithm.

Zoom In Zoom Out Reset image size

The error metric is the ${l}_{2}$ norm of the residual vector. The formula of the error metric is as follows:

$$\begin{eqnarray}&&{\rm{Error}}\,{\rm{metric}}={\parallel ({R}_{l}{{\boldsymbol{f}}}_{1}^{n},{R}_{l}{{\boldsymbol{f}}}_{2}^{n})-({R}_{l}{{\boldsymbol{f}}}_{1}^{* },{R}_{l}{{\boldsymbol{f}}}_{2}^{* })\parallel }_{2}\end{eqnarray} \tag{ 23 }$$

where $({R}_{l}{{\boldsymbol{f}}}_{1}^{n},{R}_{l}{{\boldsymbol{f}}}_{2}^{n})$ are the estimated values after $n(\geqslant 0)$ iterations, and $({R}_{l}{{\boldsymbol{f}}}_{1}^{* },{R}_{l}{{\boldsymbol{f}}}_{2}^{* })$ are the reference values. The convergence curve is a curve of the error metric changing with the number of iterations, which reflects the convergence trend of an algorithm. It is used for comparison of different algorithms on the convergence speed. The abscissa of the convergence curve is the number of iterations, and the ordinate of the convergence curve is the error metric. Figure 3(b) shows the convergence curves of the E-ART algorithm and the proposed algorithm. The blue curve and red curve are the convergence curves of the E-ART algorithm and the proposed algorithm, respectively.

It can be seen from figure 3(a) that the angle between the polychromatic projection curves under the two spectra is very small. The E-ART algorithm updates the iteration point in the way of orthogonal projection. The convergence speed from the initial iteration point to the true solution is very slow; moreover, it cannot converge to the true solution, even after 50 iterations. In the proposed algorithm, the OPMT is used to accelerate the convergence speed from the initial iteration point to the true solution. And it converges to the true solution after about 20 iterations. The convergence curve shows that the convergence speed of the proposed algorithm is clearly faster than that of the E-ART algorithm.

A slice of the 3D FORBILD thorax phantom is used in the simulation. Figure 4(a) shows an image of the phantom with resolution $512\times 512.$ For detailed information about the phantom please refer to http://www.imp.uni-erlangen.de/phanto ms/index.htm. The parameters of the scanning configuration are set as follows: the distance from the x-ray source to the center of the turntable (SOD) is 490 mm, and the distance from the x-ray source to the detector (SDD) is 880 mm. The linear detector is composed of 512 detector cells, and the length of each cell is 0.2mm.

Zoom In Zoom Out Reset image size

The projections corresponding to the two x-ray spectra are acquired in an alternating mode, i.e. the projections are geometrically inconsistent. A total of 720 projections are collected under each x-ray spectrum for one full turn.

In the experiments, one iteration means all projections are used once in the reconstruction. The convergences of the ASD-NC-POCS, the E-ART and the proposed OPMT algorithm are optimized by updating the reconstructed basis material images alternatively with low-energy and high-energy polychromatic projections. The sizes of all the reconstructed images were $512\times 512.$ Every pixel in the initial basis material images is set to 0.

Figure 5 shows a comparison of the reconstruction results respectively obtained with the ASD-NC-POCS, the E-ART and the proposed algorithm from noise-free polychromatic projections. The first column shows water material density images, the second column shows bone material density images and the third column shows the monochromatic images of linear attenuation coefficients at energy 70 keV. The fourth column is the absolute values of the difference between the water density image and the reference image of water density. The first and second rows are the reconstruction results obtained after 15 and 1000 iterations of the ASD-NC-POCS algorithm, respectively. The third and fourth rows are the reconstruction results obtained after 15 and 500 iterations of the E-ART algorithm, respectively. The fifth and sixth rows are the reconstruction results after 15 and 50 iterations of the proposed algorithm.

Zoom In Zoom Out Reset image size

The reconstructed results from the ASD-NC-POCS and E-ART algorithms after 15 iterations are chosen for comparison with the results from the proposed algorithm after the same iterations. The results from the ASD-NC-POCS algorithm after 1000 iterations and the E-ART algorithm after 500 iterations see no improvement with more iterations, which means they are very close to convergence. And the proposed algorithm is converged after 50 iterations. So, in the simulated data experiments, we choose the results of the ASD-NC-POCS algorithm after 1000 iterations and the E-ART algorithm after 500 iterations for comparison with the results from the proposed algorithm after 50 iterations.

In the water density and bone density images, the pixel values represent the density ratio of the current materials to that of the standard materials. The pixel values of the monochromatic images represent the linear attenuation coefficient of the materials at the current pixel at 70 keV energy. For the bone density, water density and monochromatic images, there are unified display windows for each of the three types of images. Every column shows the display window at the same level of the images reconstructed by different methods, and every row shows the different display windows of different types of images of one method. C and W represent the window level and window width, respectively.

Figure 6 shows the zoomed-in images of the red box position in the water material images of figure 5. It can be seen from figures 5 and 6 that the structures are not well reconstructed in the images reconstructed after 15 iterations of the ASD-NC-POCS and E-ART algorithms. There is a great improvement in the qualities of the reconstruction results of the ASD-NC-POCS algorithm after 1000 iterations and the E-ART algorithm after 500 iterations. Moreover, due to the constraint of total variation minimization of basis material images in the ASD-NC-POCS algorithm, the convergence speed is slightly slower than that of the E-ART algorithm. In contrast, it is obvious that the images reconstructed from the proposed algorithm after 15 iterations are much better than those from the other two algorithms after 15 iterations, even close to the ASD-NC-POCS algorithm after 1000 iterations and the E-ART algorithm after 500 iterations. The reconstructed results from the proposed algorithm after 50 iterations chosen in figure 5 are not only visually close to the results of the E-ART algorithm after 500 iterations, but also better than the results of the ASD-NC-POCS and E-ART algorithms for the quantitation, shown as peak signal to noise ratio (PSNR) and image mean-squared error (IMMSE) in table 1. Obviously, it is because, for the updating mode of basis material projections, the convergence of the OPMT adopted by the proposed algorithm is much faster than that using orthogonal projection adopted by both the ASD-NC-POCS and E-ART algorithms.

Zoom In Zoom Out Reset image size

The phantom data and scanning parameters used in this section are the same as those in section 3.1.2. Noisy polychromatic projections are simulated by adding Poisson noises to the noise-free data, i.e.

$$\begin{eqnarray}&&{p}_{noisy}=-\,\mathrm{ln}\left(\displaystyle \frac{{\rm{Poissrnd}}({{\rm{I}}}_{0}\times {e}^{-p})}{{{\rm{I}}}_{0}}\right)\end{eqnarray} \tag{ 24 }$$

where $p$ and ${p}_{noisy}$ denote noise-free and noisy projections, respectively. Here, ${{\rm{I}}}_{0}$ is the x-ray intensity for each x-ray path. We set ${{\rm{I}}}_{0}={10}^{6}$ in the noisy experiments. We use the function ${\rm{Poissrnd}}({{\rm{I}}}_{0}\times {e}^{-p})$ to generate random numbers from the Poisson distribution with the mean parameter ${{\rm{I}}}_{0}\times {e}^{-p}.$

Figure 7 shows the images of the test phantom reconstructed from noisy polychromatic projections with the ASD-NC-POCS, the E-ART and the proposed algorithm. The first row to the fourth row are the reconstruction results of the ASD-NC-POCS algorithm after 15 iterations and 1000 iterations, and the E-ART algorithm after 15 iterations and 500 iterations, respectively. The last two rows are the results of the proposed algorithm after 15 and 50 iterations. Figure 8 shows the zoomed-in images of the corresponding red box position in the water material images of figure 7.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 7 shows that when the polychromatic projection data is noisy, the results of the proposed algorithm after 15 iterations are still better than those of the ASD-NC-POCS and E-ART algorithms after 15 iterations, even close to the ASD-NC-POCS algorithm after 1000 iterations and the E-ART algorithm after 500 iterations. The reconstructed results from the proposed algorithm after 50 iterations chosen in figure 7 are visually close to the results of both the ASD-NC-POCS algorithm after 1000 iterations and the E-ART algorithm after 500 iterations. The PSNR and IMMSE in table 1 show that the results of the proposed algorithm after 50 iterations are better than those of the E-ART algorithm after 500 iterations. Due to the constraint of total variation minimization of basis material images in the ASD-NC-POCS algorithm, when it converges, it performs better in PSNR and IMMSE than the proposed algorithm, as shown in table 1. However, when the ASD-NC-POCS algorithm converges, the number of iterations is dozens of times that of the E-ART algorithm.

Figure 10 shows the images of the Al-water phantom reconstructed from polychromatic projections with the ASD-NC-POCS, the E-ART and the proposed algorithm. The first row to the fourth row are the reconstruction results of the ASD-NC-POCS algorithm after 5 and 100 iterations, and the E-ART algorithm after 5 and 100 iterations, respectively. The last row shows the results of the proposed algorithm after five iterations.

Zoom In Zoom Out Reset image size

It can be seen from figure 10 that the basis material images of the ASD-NC-POCS and E-ART algorithms after five iterations are far from convergence. And there is a great improvement in the qualities of the reconstruction results of the ASD-NC-POCS and E-ART algorithms after 100 iterations. Due to the constraint of total variation minimization of basis material images in the ASD-NC-POCS algorithm, the results of the ASD-NC-POCS algorithm are smoother than the results of the E-ART algorithm. In contrast, the proposed algorithm after five iterations can obtain similar reconstruction results compared with the results of both the ASD-NC-POCS and the E-ART algorithms after 100 iterations.

For further evaluation of the performance of the proposed algorithm, 1D profiles of the water density images and aluminum density images reconstructed by different algorithms are shown in figures 11(a) and (b), respectively. In figure 10, the dotted line indicates the 1D profiles of the images. From figure 11 we can see the proposed algorithm after five iterations has similar precision compared with the ASD-NC-POCS and E-ART algorithms after 100 iterations.

Zoom In Zoom Out Reset image size

Figure 12 illustrates a comparison of the reconstruction results obtained with the ASD-NC-POCS, the E-ART and the proposed algorithm from bone–water polychromatic projections, respectively. The first column shows water material density images, the second column shows bone material density images and the third column shows the monochromatic images of linear coefficients at energy 70 keV. The first and the second rows are the reconstruction results obtained after 5 and 100 iterations of the ASD-NC-POCS algorithm, respectively. The third and the fourth rows are the reconstruction results obtained after 5 and 100 iterations of the E-ART algorithm, respectively. The last row is the reconstruction results of the proposed algorithm after five iterations.

Zoom In Zoom Out Reset image size

As can be seen from figure 12, when the ASD-NC-POCS and E-ART algorithms are iterated five times, the water basis material and bone basis material are not separated. After 100 iterations, the quality of reconstructed water density images and bone density images is significantly improved, but there are still two details worth focusing on. Firstly, from the red arrow of the water density images, we can see that there are bright edges in the reconstruction results of the E-ART algorithm, and with the increase in iterations, the bright edges become much clearer. This bright edge structure is very similar to the profile belonging to the bone basis material. Through the analysis, we found that the reason for the bright edges is that at the beginning of the E-ART iteration, the reconstructed images of the basis material are far away from the true solution. At this time, there are some miscalculated structures. As mentioned, the MCT reconstruction problem is ill-posed, and minor corrections at each step cannot eliminate these existing pseudo structures. There are also bright edges in the results of the ASD-NC-POCS algorithm, but the structure of bright edges is less than that of the E-ART algorithm due to the total variation minimization constraint in the iterative process of the algorithm. At the same time, the water basis material images obtained by the ASD-NC-POCS algorithm are smoother than those of the E-ART algorithm, and some trabecular structure details of the ASD-NC-POCS algorithm are more blurred than those of the E-ART algorithm. Secondly, from the monochromatic images, we can see that there are obvious strip artifacts in the water material images obtained by both the ASD-NC-POCS and the E-ART algorithms. Compared with the results of the ASD-NC-POCS and E-ART algorithms after 100 iterations, the proposed algorithm obtains better reconstruction results for five iterations. There is no bright edge at the red arrow of the water density image, and no obvious strip artifacts in the monochromatic image.

It is worth noting that the results of the ASD-NC-POCS and E-ART algorithms after 100 iterations and the proposed algorithm after five iterations have all converged, but the water density images are blurred, while the bone density images are relatively clear. The reason may be that the water used in this experiment is pure water, and its mass attenuation coefficient matches the values found on the NIST website. However, the mass attenuation coefficient of the sheep bone, which is submerged in water, is different from that on the NIST website. So, pure water can be correctly divided into a water density image, while sheep bone cannot be completely separated in bone density images, and there are still residuals in the water density images. As a result, the contrast of the water density images is reduced, and the images are slightly blurred.