# A projection-domain iterative algorithm for metal artifact reduction by minimizing the total-variation norm and the negative-pixel energy – Visual Computing for Industry, Biomedicine, and Art

#### ByGengsheng L. Zeng

Jan 2, 2022

This paper is inspired by the two observations that the metal artifacts often result in dark/bright streaks radiating from the metals [25] and negative image pixels around the metals [26]. The streaking artifacts increase the TV norm of the non-metal regions of the image. The attenuation coefficients cannot be negative. These two effects were dealt with in two papers, separately [16, 26]. The goal of this current paper is to combine these two methods into one.

Let the measured sinogram (i.e., line integrals) be P, which can be divided into two parts: the metal affected part PM and the metal unaffected part PnotM. This separation can be readily achieved by performing a raw filtered backprojection (FBP) reconstruction, thresholding, and forward projection. The objective function is hence

$${displaystyle begin{array}{c}mathrm{Objectivefunction}=T(P)={beta}_1{T}_1+{beta}_2{T}_2\ {}={beta}_1{leftVert Gleft{F(P)right}rightVert}_{TV}+{beta}_2{leftVert min Big{0,F(P)Big}rightVert}^2end{array}}$$

(1)

where the first term is the TV norm, and the second term is the squared L2 norm. In Formula (1), F denotes the FBP reconstruction operator, G is the operator that removes the metals from the image, and β1 and β2 are the parameters balancing the two norms. In the gradient descent algorithm, β1 and β2 also act as the step sizes (also, known as the relaxation factors).

We minimize this objective function with respect to unknowns PM with a gradient descent algorithm. It is noticed that neither of the two terms in Formula (1) is differentiable. Sub-gradients are used to combat this difficulty [27].

A sub-gradient is a ‘surrogate’ of the gradient: it coincides with the gradient whenever a gradient exists, and it generalizes the notion of gradient at points where the function is non-differentiable. Let us consider a point of interest and the function is non-differentiable at this point. A sub-gradient is evaluated as follows. We find the left-gradient by taking the limit from the left and the right-gradient by taking the limit from the right, which are also known as the directional gradients. Any values between these two left and right values can be used as the sub-gradient.

We now use the absolute function |x| to illustrate how to find the sub-gradient. When x > 0, the sub-gradient is the same as the gradient: 1. When x < 0, the sub-gradient is the same as the gradient: −1. When x = 0, the function |x| is non-differentiable. The left-gradient is −1 and the right-gradient is 1. The sub-gradient at x = 0 can be any value in the interval [− 1, 1].

In our implementation, three values of β1 were tested: 0, 0.002 and 0.004; two values of β2 were tested: 0 and 5.

### The FBP algorithm

In Formula (1), F(P) is a linear algorithm that can be deposed into two steps. The first step is to convolve the sinogram P along the detector dimension with a convolution kernel h(n):

$$h(n)=left{begin{array}{ccc}frac{1}{4}& if& n=0,\ {}-frac{1}{{left( npi right)}^2}& if& n is odd,\ {}0& otherwise.& end{array}right.$$

(2)

The second step is the backprojection. Let the FBP reconstruction from P be X. If both P and X are represented in the vector forms, F(P) can be written as matrix multiplication

where F is a matrix.

For a given threshold value t, a metal image is obtained by modifying the image X. If a pixel in X is greater than t, its value is set to 1, otherwise its value is set to 0. The forward projection of X generates a ‘shadow’ in the projection domain. If the shadow values are positive, the corresponding sinogram values are called metal-affected projections, PM. If the shadow values are zero, the corresponding sinogram values are called metal-not-affected projections, PnotM. The sinogram P is thus divided into two parts: PM and PnotM. The projection values in PM are treated as variables in the proposed algorithm.

In Formula (1), the metal removed FBP image, G{F(P)} = G{X}, is almost the same as X, except that if a pixel value is greater than t, its value is set to 0.

In Formula (1), min{0, F(P)} = min {0, X} is almost the same as X, except that if a pixel value is positive, its value is set to 0. Here, the ‘min’ function is an element-wise function.

### Optimization of the objective function (Formula 1) by the gradient descent algorithm

A gradient descent algorithm to minimize the objective function T(P) in Formula (1) can be expressed as

$${P}_M^{left(k+1right)}={P}_M^{(k)}-Dleft[{beta}_1nabla {T}_1(P)+{beta}_2nabla {T}_2(P)right]$$

(4)

where the super script (k) is the iteration index. The projection vector P consists of two parts: the metal affected part PM and the metal not-affected part PnotM. The metal not-affected part PnotM does not get updated from iteration to iteration. The matrix D in Formula (4) is a diagonal matrix with 0’s and 1’s, which discards the entries in PnotM. The parameter β1 and β2 in Formula (4) controls the step sizes of the gradient descent algorithm for the two different norms, respectively. In Formula (4), T1 is the gradient of the TV norm of the image with the metals removed, and T2 is the gradient of the energy of the negative pixels in the image. These two gradients are calculated with respect to the projections P. These two gradients are discussed further in the following two sub sections.

### Minimization of the first term in the objective function (Formula 1)

Let Y = G{F(P)} = G{X} be the metal-removed FBP reconstruction. The image Y is represented in a two-dimensional array and yi,j is its pixel value at the ith row and jth column. The TV norm of Y is defined as

$${T}_1={sum}_{i,j}sqrt{{left({y}_{i,j}-{y}_{i,j+1}right)}^2+{left({y}_{i,j}-{y}_{i+1,j}right)}^2}$$

(5)

The partial derivative of T1 with respect to pixel (i, j) (if exists) is readily calculated as

$${displaystyle begin{array}{c}{u}_{i,j}=frac{partial {T}_1}{partial {y}_{i,j}}\ {}begin{array}{c}=frac{left({y}_{i,j}-{y}_{i,j+1}right)+left({y}_{i,j}-{y}_{i+1,j}right)}{sqrt{{left({y}_{i,j}-{y}_{i,j+1}right)}^2+{left({y}_{i,j}-{y}_{i+1,j}right)}^2}}\ {}begin{array}{c}+frac{y_{i,j}-{y}_{i,j-1}}{sqrt{{left({y}_{i,j-1}-{y}_{i,j}right)}^2+{left({y}_{i,j-1}-{y}_{i+1,j-1}right)}^2}}\ {}+frac{y_{i,j}-{y}_{i-1,j}}{sqrt{{left({y}_{i-1,j}-{y}_{i-1,j+1}right)}^2+{left({y}_{i-1,j}-{y}_{i,j}right)}^2}}end{array}end{array}end{array}}$$

(6)

When the quantity under a square root is zero, the gradient ui,j is not well-defined. In this case, we use the sub-gradient and assume the sub-gradient ui,j to be zero

The gradient of T1 given in Formula (6) is with respect to the image pixels yi,j. However, the gradient decent algorithm (Formula 4) requires the gradient of T1 be with respect to the projections P.

Realizing that Y is in the image domain, the algorithm variables are in the projection domain related by the mapping

where, A represents the forward projection matrix. The matrix A maps the image into its projections. This same matrix A can map the gradient U into the projection domain as AU, where the entries of U are ui,j. We have

$$nabla {T}_1(P)= AU$$

(9)

### Minimization of the second term in the objective function (Formula 1)

Let Z = min {0, F(P)} = min {0, X} be the negative pixels in the FBP reconstruction X. The energy of the negative image pixels is the square of the L2 norm of Z given as

$${T}_2={leftVert ZrightVert}_2^2={leftVert min left{0, FPright}rightVert}_2^2$$

(10)

To find the gradient T2(P) , that is, ∂T2/∂pj, is not straightforward, because the ‘min’ function in Formula (10) makes T2 non-differentiable. We can use the subdifferential concept to find the gradient of ∂T2/∂pj as [27]

$$nabla {T}_2=2{F}^Tmin left{0, FPright}=2{F}^TZ$$

(11)

where FT represents forward projection followed by ramp filtering. The operator FT maps an image into the projection domain.

### Proposed gradient descent algorithm to estimate PM

Combining the gradients above, the proposed gradient descent iterative algorithm (Formula 4) to estimate PM can given as

$${P}_M^{left(k+1right)}={P}_M^{(k)}-Dleft[{beta}_1tanh (AU)+{beta}_2 {F}^TZright]$$

(12)

where tanh(T1) is used in place of T1 in Formula (4). The purpose of the hyperbolic tangent function ‘tanh’ is to hard limit the TV gradient values, making the iterative algorithm more stable.

The proposed algorithm is different from the commonly used iterative algorithms in the sense that the sinogram values in a common iterative reconstruction algorithm are never altered. Also, the proposed algorithm does not have a data fidelity term in the objective function (Formula 1). Our final image is obtained by the FBP algorithm, after PM is modified. The initial value of the metal affected part PM is the measured value.

We tested the proposed algorithm with some unknown airport bags, which contained some unknown metals. The projection sinograms were provided by the US Department of Homeland Security. The data sets were acquired with unknown kVp and unknown current settings. The original sinograms were converted into parallel-beam geometry with 180 views over 180° and with 597 detection channels (i.e., detector bins).