Abstract
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.
1 Introduction
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.
2 Methods
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 n0, 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 aij 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 pj 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 pj 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.
3 Results
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.
Algorithm 1: Pseudo-code representation of the proposed algorithm
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
The used phantom. Ground truth image of the image that is to be reconstructed, which shows a region around the pectoral girdle. Segmentation of the metal object. The used prior image, which is located approx. 1 cm underneath the ground truth image. The difference of the ground truth image and the prior image.
In Figure 3 the weights for each patch within the transformed prior image Γ(g, y) is shown. The weights correspond to the pixel fx, 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.
Comparison of different patches between the image to be reconstructed, on the left side, and the prior image, on the right side. The search window is visualized by a dotted line. The red rectangle shows the patches fx and Γ(g, y)ηy, respectively.
Weights for the different patches within the search window which is visualized in Fig. 2.
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.
Progression of the reconstruction. In each iteration a new non-local prior image is calculated. In the last raw a comparison of the proposed algorithm and the linear interpolation approach is shown.
4 Conclusion
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.
Author’s Statement
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.
References
[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.10.1118/1.2836423Search in Google Scholar PubMed PubMed Central
[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.10.1160/ME9041Search in Google Scholar PubMed
[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.10.1145/279232.279236Search 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.10.1109/ICCV.1998.710815Search 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.10.1118/1.2955743Search in Google Scholar PubMed PubMed Central
[6] W. Kalender and R. Hebel, “Reduction of CT artifacts caused by metallic implants,” Radiology, vol. 164, no. 2, pp. 576–577, 1987.10.1148/radiology.164.2.3602406Search in Google Scholar PubMed
© 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.
