The presence of high-density objects remains an open problem in medical CT imaging. Data of projections that passing through such objects are dominated by noise. Reconstructed images become less diagnostically conclusive because of pronounced artifacts that manifest as dark and bright streaks. A new reconstruction algorithm is proposed, which incorporates information gained from a prior image. Based on a non-local regularization, these information are used to reduce streaking artifacts. In an iterative scheme, the prior image is transformed in order to match intermediate results of the reconstruction by solving a registration problem. During iterations, temporally appearing artifacts are reduced with a bilateral filter and projection values passing through high-density objects are replaced by new calculated values, which are used further on for the reconstruction. Results show that the proposed algorithm significantly reduces streaking artifacts.
Computed tomography (CT) remains one of the key methods in medical imaging. The image quality of reconstructed CT slices is reduced by the occurrence of different artifacts, which are caused by physical phenomena such as scattering, beam hardening, noise, or total absorption. One of the main sources of artifacts in clinical practice are metal implants and surgical instruments within the area that is imaged. These objects generate dark and bright streaking artifacts, which obstruct the assessment of the anatomy of the patient and reduce the diagnostic value of the images.
Since the publishing of the PICCS algorithm by Chen et al. the integration of prior images in the reconstruction of CT images has become a fast developing topic in the field [1]. A reconstruction algorithm that reduces metal artifacts by integrating a prior image in terms of a non-local regularization term is presented. The regularization term penalizes intensity variations between the image to be reconstructed and a prior image. The prior should hold information from an image that looks similar to the image that is to be reconstructed. However, the goal of the regularization term is to keep anatomical information of the original image while reducing artifacts. In order to find a correct alignment of the prior image and the image that is to be reconstructed a iterative reconstruction scheme is chosen where in certain iteration steps a registration step is incorporated. Furthermore, during these iterations temporally appearing artifacts are reduced with a bilateral filter and projection values passing through high-density objects are replaced by new calculated values, which are used further on for the reconstruction.
Given a set of intensity measurements
where f ∈ ℝ^{N} is a vector that consists of the expected attenuation coefficients [2]. The number of photons that are detected in the absence of absorption is denoted by n_{0}, the total number of projections is denoted by M, and the number of pixels within the image is denoted by N. The forward projection step is given by
where a_{ij} is an element of the system matrix A ∈ ℝ^{M×N}. Based on these formulations, let
be the optimization problem that is to be solved in order to reconstruct an image of the tomographed object, where R(f, Γ(g, y)) is a regularization term that is added to the objective function and a regularization parameter that controls the influence of the regularization. Here a non-local regularization is proposed, which is defined as
with
where Γ(g, y) is the transformed prior image g with the transformation parameter. Within (5) _{x} denotes a patch window around pixel x, 𝒩_{x} denotes a search window around pixel x, and ||·||_{p} denotes the Minkowski distance of order p. Furthermore is λ ∈ {0, 1}^{N} a mask with
which forces the regularization to ignore all pixels where the prior image holds no information.
The optimization problem (3) is solved by the l-BFGS-b algorithm [3]. In order save computation time while preserving accuracy, the registration problem
where 𝒟 : ℝ^{2N} → ℝ denotes a distance measure, is solved only in specific iteration steps. These iteration steps are specified by a set of convergence tolerances ω^{(k + 1)} > ω^{*} with k ∈ ℕ such that ω^{(k+1)} < ω^{(k)}, where ω^{*} denotes the final convergence tolerance at which point the algorithm stops and returns the reconstructed image. While for the distance measure D mutual information Γ is chosen, the calculated transformation is based on an affine transformation.
In order to further reduce metal artifacts, an initial reconstruction based on a filtered backprojection is generated. The resulting image is thresholded in order to gain a mask μ ∈ {0, 1}^{N} that describes the position of the metal object. Based on μ the set 𝓜_{1} of projection indices is introduced, which corresponds to projection values that are not affected by metal and 𝓜_{2}, a set of projection indices which corresponds to projection values that are affected by metal, such that {1, ..., M} = 𝓜 = 𝓜_{1} ∪ 𝓜_{2} and 𝓜_{1} ∩ 𝓜_{2} = ø. Before the algorithm is started, the projection values p_{j} with j ∈ 𝓜_{2} are removed from the measurement in order to discard all projection values that are associated with the metal object. However, in the course of the algorithm, likewise the solving of the registration problem (7), the projection values p_{j} with j ∈ 𝓜_{2} are replaced by a forward projection of a filtered version of the current image f^{*}(k), which is represented by the intermediate result of (3) in iteration k. For the filtering step a bilateral filter has shown to be appropriate [4]. A bilateral filter is a non-linear local filter operator that takes into account a local intensity context and as such preserves edges. Let
and the intensity values
such that the negative log-likelihood function for iteration k > 1 is defined as
In algorithm 1 the complete algorithm is given by a pseudo-code representation.
In Figure 1 the used phantom is shown. The phantom is generated by using the XCAT software [5] and is located around the pectoral girdle. A metal implant is simulated by adding an artificial object within the left humerus.
In Figure 2 a region of interest of the ground truth image and the prior image is shown. The search window 𝓝_{x} is visualized by a dotted line. The patches f_{ρx} and Γ(g, y)_{ηy} are visualized by a red rectangle.
1: choose initial point f*(0)
2: choose convergence tolerance ω^{(1)}
3: choose final convergence tolerance ω_{*}
4: set
5: for k = 1, 2,... do
6: start with the image f*^{(k−1)} and find an approximated solution f*^{(k)} with the tolerance ω^{(k)} of the problem f*^{(k)} = arg min (l^{(k)}(f) + δR(f, Γ(g, y)))
7: if ω^{(k)} ≤ ω* then
8: stop and return f*(k)
9: else
10: solve registration problem
11: define ω^{(k+1)} such that ω^{(k+1)} < ω^{(k)}
12: calculate bilateral filtered image
13: calculate new projection values
14: end if
15: end for
In Figure 3 the weights for each patch within the transformed prior image Γ(g, y) is shown. The weights correspond to the pixel f_{x}, which is shown in Figure 2. A red pixel symbolizes a higher weight for the corresponding pixel, while a dark blue color symbolizes a small weight near zero. Note that due to the mask λ in (4), pixels within the prior image that do not get a weight that is higher than zero will have no influence on the regularization.
Figure 4 shows the progression of the reconstruction algorithm. In each iteration a new prior image is calculated. If the newly calculated prior image gets a pixel value of zero assigned due to the weights that are based on the comparison between the intermediate result of the reconstruction and the prior image (see Figure 2), the value has no influence on the regularization due to the mask in (4). Five iterations are needed until the gradient of the objective function (3) becomes smaller than 10^{−}9. The last used prior image shows only small holes where the corresponding patches are all weighted with zero.
The last raw in Figure 4 shows the linear interpolation approach compared to the proposed algorithm [6]. The final reconstruction shows significantly less artifacts compared to the linear interpolation approach. Most importantly, please note that small anatomical details are preserved and are not suppressed by the regularization term.
A metal artifact reduction algorithm is proposed, which integrates information based on a prior image. A transformation between intermediate results of the reconstruction and the prior image is found. The result of this registration step is used in order to calculate a new prior image.. Finally, a regularization term penalizes intensity variations between the prior and the intermediate results. In an additional step, new projection values that are associated with the metal object, are computed based on filtered intermediate results and combined with the original measurement. The final reconstruction result features significantly less streaking artifacts compared to the linear interpolation approach.
Conflict of interest: Authors state no conflict of interest. Material and Methods: Informed consent: Informed consent has been obtained from all individuals included in this study. Ethical approval: The research related to human use has been complied with all the relevant national regulations, institutional policies and in accordance the tenets of the Helsinki Declaration, and has been approved by the authors’ institutional review board or equivalent committee.
[1] G.-H. Chen, J. Tang, and S. Leng, “Prior image constrained compressed sensing (PICCS): a method to accurately reconstruct dynamic CT images from highly undersampled projection data sets,” Medical Physics, vol. 35, no. 2, pp. 660–663, 2008. Search in Google Scholar
[2] M. Oehler and T. M. Buzug, “Statistical image reconstruction for inconsistent CT projection data,” Methods of Information in Medicine, vol. 46, no. 3, pp. 261–269, 2007. Search in Google Scholar
[3] C. Zhu, R. H. Byrd, P. Lu, and J. Nocedal, “Algorithm 778: L-BFGSB: Fortran subroutines for large-scale bound-constrained optimization,” ACM Transactions on Mathematical Software, vol. 23, no. 4, pp. 550–560, 1997. Search in Google Scholar
[4] C. Tomasi and R. Manduchi, “Bilateral filtering for gray and color images,” in Sixth International Conference on Computer Vision. IEEE, 1998, pp. 839–846. Search in Google Scholar
[5] W. Segars, M. Mahesh, T. Beck, E. Frey, and B. Tsui, “Realistic CT simulation using the 4D XCAT phantom,” Medical Physics, vol. 35, no. 8, pp. 3800–3808, 2008. Search in Google Scholar
[6] W. Kalender and R. Hebel, “Reduction of CT artifacts caused by metallic implants,” Radiology, vol. 164, no. 2, pp. 576–577, 1987. Search in Google Scholar
© 2015 by Walter de Gruyter GmbH, Berlin/Boston
This article is distributed under the terms of the Creative Commons Attribution Non-Commercial License, which permits unrestricted non-commercial use, distribution, and reproduction in any medium, provided the original work is properly cited.