PDS-MAR: a fine-grained Projection-Domain Segmentation-based Metal Artifact Reduction method for intraoperative CBCT images with guidewires

This study focuses primarily on the segmentation of metallic guidewires. Guidewires are tubular objects in the projection domain, so techniques for enhancing vascular objects can be applied. Here we refer to the neurite enhancement filtering method proposed in [32] to enhance guidewires. This approach relies on a modified Hessian matrix $H^{\prime}(x)$ defined with the following equation

$$H^{\prime}_{P*G_{\sigma}}(x)=H_{P*G{\sigma}}(x)+\alpha R^{T}_{\pi/2}H_{P*G{\sigma}}(x)R_{\pi/2}$$

(1)

where $P*G_{\sigma}$ denotes the result of the input projection view $P$ Gaussian-smoothed with $\sigma$ variance, $H_{f*G{\sigma}}$ represents the Hessian matrix of the smoothed image, $R_{\pi/2}$ stands for the rotation matrix with angle $\pi/2$. Let $\lambda^{\sigma}_{1}(x)$ and $\lambda^{\sigma}_{2}(x)$ be the two eigenvalues of the Hessian matrix $H^{\prime}_{P*G_{\sigma}}(x)$ and we have $|\lambda^{\sigma}_{1}(x)|>|\lambda^{\sigma}_{2}(x)|$, the output vesselness is then defined by

$$V_{\sigma}(x)=\begin{cases}\lambda^{\sigma}_{1}(x)/\lambda^{\sigma}_{max},&\lambda^{\sigma}_{max}>0\\ 0,&\lambda^{\sigma}_{max}\leq 0\end{cases}$$

(2)

where $\lambda^{\sigma}_{max}$ refers to the largest $\lambda^{\sigma}_{1}(x)$ among all pixels in the image. Eq. 2 normalizes vesselness image to the range of $[0,1]$. To better model guidewires with different cross-sectional radii, vesselnesses under a set of different Gaussian variances $\Sigma={\sigma_{i}},i=1,2,...,N$ are calculated, and the enhancement result is set as the maximal vesselness value on each pixel.

$$T(x)=max_{\sigma\in\Sigma}V_{\sigma}(x)$$

(3)

Even though the tubular enhancement filtering step substantially enhances guidewire-affected regions, it also enhances other image boundaries, such as bone and body boundaries. A simple threshold $Th_{enh}$ will include a lot of undesired regions. Fortunately, the metal traces from the image domain segmentation are still available for guidance. The image domain metals are segmented with a threshold $Th_{metal}$ and then forward-projected into the image domain. Points with values greater than 0 are considered seed points for region-growing, and the resultant binary masks are considered metal-affected region masks.

Since guidewires are typically inserted horizontally, conventional 1-D linear interpolation in rows is likely to lose a great deal of tissue information. To enhance inpainting precision, we use a Delaunay triangulation-based 2-D interpolation method.

To reduce computation costs, the input binary mask is first subdivided into connected components. The outer boundary of each connected component is extracted by subtracting itself from the morphologically dilated mask. The area within the outer boundary is then triangulated using the Delauney method [33]. For each point $x$ inside a triangle facet $\Delta_{x_{1}x_{2}x_{3}}$, the interpolated value $P(x)$ is defined as a linear weighting of value on the vertices.

$$P(x)=w_{x_{1}}P(x_{1})+w_{x_{2}}P(x_{2})+w_{x_{3}}P(x_{3})$$

(4)

where the weights are calculated with the following equations.

$$\begin{cases}w_{x_{1}}=\frac{\vec{x_{3}x}\times\vec{x_{3}x_{2}}}{\vec{x_{3}x_{1}}\times\vec{x_{3}x_{2}}}\\ w_{x_{2}}=\frac{\vec{x_{3}x}\times\vec{x_{3}x_{1}}}{\vec{x_{3}x_{2}}\times\vec{x_{3}x_{1}}}\\ w_{x_{3}}=1-w_{x_{1}}-w_{x_{2}}\end{cases}$$

(5)

With metal traces interpolated in the projection domain, metal-free images can be reconstructed with the widely-accepted FDK algorithm [34].

Despite that we have already reconstructed a set of metal-free images, all information regarding the locations of metallic guidewires has been lost. Here, we present a method to restore image-domain metal masks $I_{M}(x)$ from projection-domain metal traces $P_{M}(u,v,\theta)$ based on multiplicative-form backprojection. Assuming that we have the perfect metal traces in the projection domain, forward-projecting an image-domain metal voxel will always find a projection-domain point inside the metal traces, whereas projecting a voxel outside the metallic convex hull will find points outside the metal traces in at least some views. The metal mask can therefore be reconstructed using the following equation:

$$I_{M}(x)=\prod_{\theta=0}^{2\pi}P_{M}(u(x,\theta),v(x,\theta),\theta)$$

(6)

where $u(x,\theta)$ and $v(x,\theta)$ are the detector coordinate at view angle $\theta$ corresponding to image point $x$.

Considering false negative points in metal trace segmentation which may introduce a large number of false negatives in $I_{M}$, we replace $P_{M}$ in Eq. 6 with a soft mask $P^{soft}_{M}$ defined as follows

$$P^{soft}_{M}=(1-\gamma)P_{M}+\gamma(1-P_{M})$$

(7)

where $\gamma$ is a weighting parameter in $(0,1)$. We usually choose $\gamma$ values close to 1 (like 0.9). After soft-masking, metal points get large $I_{M}$ values even if affected by false negatives on limited views while non-metal points get values close to 0. Metal masks are therefore segmented with a threshold $Th_{mask}$.

All algorithms are implemented and tested in MATLAB R2021 with some operators implemented in C++/CUDA and compiled with mex/mexcuda in MATLAB. As for parameter settings, we set $\alpha=1/3$, $\Sigma=\{1,3,5,7,9\}$, $Th_{enh}=0$, $\gamma=0.9$ and $Th_{mask}=0.5$ for all experiments, as the results are not sensitive to these parameters. Parameter $Th_{metal}$ was set to 3000HU by default, and its settings will be discussed in 3.5.

To quantitatively evaluate the performance of PDS-MAR, we simulated CBCT scans under a clinical intraoperative C-arm CBCT geometry with metallic guidewires inserted. XCIST [35], an open-source X-ray CT simulation toolkit, is used to simulate the projection data using a cone beam geometry. The detailed geometrical parameters are presented in Table 1, where $DSD$ represents the distance between source and detector, $DSO$ stands for distance between source and rotation axis. We used the default 110kV spectrum provided in XCIST without any pre-filtration for simulations. For the digital phantom, the voxelized female chest phantom provided along with XCIST (https://github.com/xcist/phantoms-voxelized) is adopted for simulation, which is also a part of the XCAT phantom dataset by Segars et al. [36].

To simulate CBCT metal artifact, an additional material mask is manually outlined from a set of intraoperative clinical CBCT images using guidewires and applied to the digital phantom as the distribution of material ’iron’. Additionally, we have modified the phantom position offset to simulate CBCT data truncation during spine procedures.

As the performance of MAR algorithms mainly depends on the projection-domain interpolation, here we first evaluate the projection-domain accuracy. Classical MAR methods (MAR-LI [17] and NMAR [18]) are also compared here to show the performance of the proposed algorithm. To evaluate the effectiveness of each component, ablation analysis with two different configurations is also included here: 1) MAR-tri, MAR with image-domain metal trace segmentation and triangulation-based interpolation; 2) PDS-MAR-LI, MAR with projection-domain metal trace segmentation and linear interpolation.

Fig. 3 provides a visual comparison of different algorithms, and results on 3 different views are depicted. The reference projection views are simulated using the digital phantom without metal implants under the same geometry. Compared to other algorithms, the proposed algorithm generate corrected projection views closest to the reference views. To better view the details, Fig. 3 also shows the enlarged projection-domain blocks corresponding to the red rectangles. The errors of MAR algorithms mainly come from two aspects, false segmentation and non-ideal interpolation. MAR-LI, MAR-tri and NMAR adopt image-domain segmentation, which introduces obvious false segmentations (white arrows in Fig. 3). Despite the fact that metals can be easily segmented with a threshold on reconstructed phantom images, there is no way for image-domain methods to get metal traces outside the reconstruction FoV. Besides, 1-D linear interpolation produces over-smoothed rows (green arrows in Fig. 3). This characteristic is particularly harmful when metal implants are long in the horizontal direction (View 0 in Fig. 3), which is, however, common for guidewires. Triangulation-based interpolation, on the other hand, inpaints the void regions according to 2-D information and results in interpolated views much closer to references. There are still some defects shown in the results of the proposed algorithm (orange arrows in Fig. 3). Still, generally, the proposed algorithm produces the best-interpolated projection data among all algorithms.

To quantitatively evaluate the performance of the proposed algorithm, several metrics are calculated between the outputs and the references regarding both segmentation accuracy and final interpolation quality. For segmentation, we compute Precision, Recall, and Dice index between the generated metal trace and the simulated projection data with the metal object only. Commonly used metrics including Root Mean Square Error (RMSE), Peak-Signal-to-Noise Ratio (PSNR), and Structural Similarity Index Metric (SSIM) are also calculated between interpolated projection data and reference data to compare interpolation similarity. The results are displayed in Table 2. The proposed projection-domain segmentation method achieves a dice index of 0.9188, whereas conventional image-domain segmentation only gets 0.8141. The projection-domain method achieves a $\sim 4\%$ improvement in segmentation precision and a $\sim 14\%$ improvement in recall over the image-domain method. The significant increase in recall is primarily attributable to the successful segmentation of metal traces outside the reconstruction field of view. The proposed method also obtains the highest metric values for final interpolation similarity compared to all other comparing methods. From the results, we get the following observations. (1) PDS-MAR-LI/PDS-MAR achieves better interpolation than MAR-LI/MAR-tri. This proves that projection-domain segmentation is more appropriate than image-domain segmentation. (2) Triangulation-based interpolation methods (MAR-tri/PDS-MAR) achieve better results than linear interpolation methods (MAR-LI/PDS-MAR-LI). (3) NMAR gets the worst results on all metrics. This may result from poor tissue/bone segmentation in the image domain.

Fig. 4 gives a visual comparison between the reconstruction results of different algorithms. Generally, image domain segmentation-based methods (MAR-LI, MAR-tri and NMAR) suffer from severe streak artifacts resulting from metal objects out of scanning FoV. Linear interpolation-based algorithms (MAR-LI, PDS-MAR-LI) suffer from over-smoothing near metal objects, and tissue contrast details get lost. The proposed algorithm shows much-reduced streak artifacts and sharper tissue boundaries around metal implants than other algorithms. However, the proposed algorithm does not eliminate streak artifacts totally, since there are still interpolation failures in the projection domain (orange arrows in Fig. 3).

Quantitative evaluation results on the image domain are provided in Table 3. We first compared the metal mask accuracy between image domain thresholding and the proposed metal mask reconstruction method. The image domain thresholding method achieves a pretty good dice score in the phantom study. However, the success in thresholding mainly results from low image noise, low bone value and much-suppressed image artifacts in our simulated data, and cannot be reproduced on authentic CBCT images (see Fig. 1). The proposed algorithm achieves a dice score $\sim 1\%$ higher than image domain thresholding does. As for artifact reduction, we compute RMSE, PSNR and SSIM metrics between MAR results and reference images from metal-free data. The proposed algorithm shows significant improvements in all metrics compared to other algorithms.

According to the results of MAR-tri at view 0 in Fig. 3, metal traces are nearly eliminated at this view. The near-perfect interpolation results from a near-perfect segmentation at this view, indicating that metal objects within the scanning FoV in the image domain are nearly accurately segmented. However, the results of the reconstruction indicate that metal artifacts cannot be eliminated under these conditions. There are metal objects outside the scanning FoV that can be observed from other perspectives, and the artifacts introduced by these metal objects cannot be corrected using image-domain segmentation.

To prove this conjecture, we conducted an experiment on an ideal image-domain segmentation-based MAR method (ID-ideal). In ID-ideal, the image-domain metal mask is generated by masking the ground truth with scanning FoV. The metal mask is then forward-projected into the projection to produce metal traces. For interpolation, we directly replace the values within the metal traces with corresponding values in the reference projection data (simulated without metal), which should be considered perfect interpolation. The resulting projection data at a view and reconstruction outputs on 3 slices are presented in 5. From the results, we see that metal traces are not completely corrected in the projection domain, which further introduces severe streak artifacts to the image domain. There are also streak artifacts introduced by boundary precision issues (orange arrows in Fig. 5 (b,d)). Quantitative evaluations of the image domain are also provided in Table 3. The metrics attained by PDS-MAR are superior to those attained by ID-ideal.

The proposed algorithm is evaluated on a set of animal body data. The back of a deceased sheep was scanned with Kirschner wires (K-wires) inserted by an expert. These projection data were acquired using a commercial intraoperative CBCT system (B51s) from First-Imaging Medical Equipment Co., Ltd. All original data were initially pre-processed with commercial data correction software from First-Imaging (including field correction, denoising, and water beam hardening correction).

Fig. 6 provides the reconstructed results from sheep body data. Unlike in digital phantom data, results on animal body data suffer from severe segmentation failures in the image domain with a threshold. Some K-wires are partially segmented (the left one in Fig. 6 (a1-a4)) some are not found at all (the left one in Fig. 6 (b1-b4), the one in Fig. 6 (c1-c4)). The missegmentation in the image domain causes a significant metal trace missing in the projection domain which further leads to uncorrected metal artifacts in the reconstructed images. Generally, the proposed algorithm reconstructs image slices with substantially reduced metal artifacts and tissue details that are preserved. In addition, the proposed metal mask reconstruction method is able to preserve the correct image domain metal mask even when the thresholding strategy fails.

Some K-wires are nearly correctly segmented in the image domain (the right one in Fig. 6 (b1-b4)), but still get slight deviations from the real metal traces in the projection domain after the forward projection. These deviations sometimes come from the geometrical randomness of the entire CBCT system. Triangulation-based interpolation is ineffective in this situation and results in saw-toothed inpainting results (red arrows in Fig. 7 (a)). The saw-toothed projection data further introduces severe streak artifacts to the reconstructed images as a result of ramp filtering (Fig. 6 (b3) as an example). As segmentation in the projection domain is more precise than in the image domain, the proposed method generates interpolation output that is substantially smoother.

Two human-body phantoms (an abdominal phantom and a head phantom) are scanned with guidewires inserted for educational purposes. The projection data used in this section were also collected with a commercial intraoperative CBCT system from First-Imaging Medical Equipment Co., Ltd. All original data were pre-processed as described in 3.4.

The reconstruction results for the abdominal phantom are depicted in Fig. 8. Similar to the animal body results, conventional image domain segmentation suffers from segmentation faults and results in incomplete metal artifact suppression (see Fig. 8 (c2-c3)) and secondary artifacts (see Fig. 8 (a3)). The proposed algorithm shows superior ability in metal artifact reduction, and the metal masks are well-reconstructed with our multiplicative-form backprojection method. Furthermore, streak artifacts introduced by metal objects outside FoV are still observed even on slices without metal in Fov for image domain segmentation methods (bottom of Fig. 8 (b2-b3)) while being totally eliminated by the proposed method (Fig. 8 (b5)).

The reconstruction results for the head phantom are depicted in Fig. 9. These results reveal another shortcoming of image domain segmentation, namely the difficulty in threshold choice. Referring to Fig. 9 as an example, we compared results under $Th_{metal}=3000HU$ (b-e) and $Th_{metal}=6000HU$ (f-i) for different methods. Under a lower threshold, image domain methods can segment metal better. However, tooth regions are identified as metal objects due to their high attenuation properties, which introduce secondary (see Fig. 9 (b,c)). Under the higher threshold, tooth regions are not included; however, guidewires are not recognized and metal artifacts are not corrected. The proposed also suffers from segmentation faults on teeth under the threshold of 3000HU but performs well under the threshold of 6000HU. A slice in sagittal view and a slice in coronal view are displayed in Fig. 10. Besides guidewires, there are also metal screws inserted in the head phantom and the proposed algorithm can effectively deal with metal artifacts introduced by these screws.