Skip to main content
  • Research article
  • Open access
  • Published:

An improved nonlinear diffusion in Laplacian pyramid domain for cone beam CT denoising during image-guided vascular intervention

Abstract

Background

Cone-beam computed tomography (CBCT) acquisition during endovascular aneurysm repair is an emergent technology with more and more applications. It may provide 3-D information to achieve guidance of intervention. However, there is growing concern on the overall radiation doses delivered to patients, thus a low dose protocol is called when scanning. But CBCT images with a low dose protocol are degraded, resulting in streak artifacts and decreased contrast-to-noise ratio (CNR). In this paper, a Laplacian pyramid-based nonlinear diffusion is proposed to improve the quality of CBCT images.

Method

We first transform the CBCT image into its pyramid domain, then a modified nonlinear diffusion is performed in each level to remove noise across edges while keeping edges as far as possible. The improved diffusion coefficient is a function of the gradient magnitude image; the threshold in the modified diffusion function is estimated using the median absolute deviation (MAD) estimator; the time step is automatically determined by iterative image changes and the iteration is stopped according to mean absolute error between two adjacent diffusions. Finally, we reconstruct the Laplacian pyramid using the processed pyramid images in each level.

Result

Results from simulation show that the filtered image from the proposed method has the highest peak signal-noise ratio (81.92), the highest correlation coefficient (99.77%) and the lowest mean square error (27.61), compared with the other four methods. In addition, it has highest contrast-to-noise ratio and sharpness in ROIs. Results from real CBCT images show that the proposed method shows better smoothness in homogeneous regions meanwhile keeps bony structures clear.

Conclusion

Simulation and patient studies show that the proposed method has a good tradeoff between noise/artifacts suppression and edge preservation.

Peer Review reports

Background

With the progressive development of hybrid operating rooms in vascular surgery, cone-beam CT (CBCT) mounted on a C-arm becomes an increasingly commonly used imaging technology in vascular interventions, such as endovascular aneurysm repair (EVAR), thanks to its capability of 3D imaging of arterial structures and less radiation in comparison with multi-slice CT [1,2,3,4,5]. There are at least two advantages of 3D CBCT acquisition during EVAR. First, detecting any potential complication (endoleaks, stentgraft kinking) that can be treated immediately after the procedure [2,3,4,5], i.e., after the deployment of the stentgraft. Secondly, 3D CBCT acquisition (without contrast media injection) was described to fuse the pre-operative CT-scan by means of a 3D/3D rigid bone registration [6,7,8,9,10,11]. With this fusion imaging technique catheterization and stentgraft deployment can be achieved with a 3D visualization of the vascular tree and leads to a significant decrease of contrast media volume injection.

Although CBCT produces high spatial resolution and offers benefits for patients, there is growing concern on the overall radiation doses delivered to patients due to pre-, intra- and post-operative X-ray imaging during endovascular procedures [12]. As a consequence, it is very important to reduce the radiation dose related to CBCT, not only for patient care, but also for the medical staff to avoid or reduce potential determinist and stochastic risks from radiological procedure [13,14,15]. However, lowering radiation dose inevitably produces more noise, thus leading to degraded CBCT images with streak artifacts, and decreased contrast-to-noise ratio (CNR) [16]. Noise and streak artifacts suppression is therefore called as a preprocessing step to access a cleaner CBCT image, thus significantly improves the accuracy of subsequent image segmentation and registration.

Generally speaking, there are three categories for reducing noise and artifacts in CBCT: processed before reconstruction [17,18,19], during reconstruction procedure [20,21,22] and after reconstruction [23,24,25,26,27,28,29]. Since post-processing methods don’t have necessary to access to projection data, various sophisticated filters were applied on the reconstructed images, such as bilateral filtering [23], nonlocal means filtering [24], and nonlinear diffusion filtering [25, 26], most of them consider the strong structural and statistical properties of objects in image space. Recently, dictionary learning and sparse representation were used for reconstruction and enhancement for low-dose X-ray imaging [27,28,29]. However, they have the limitation of computation time.

The nonlinear diffusion filtering is a useful technique that uses an edge seeking function to encourage diffusion within regions and prohibit it across strong edges. Hence edges can be preserved while removing noise from the image. However, the gradient dependent diffusion, such as the anisotropic diffusion proposed by Perona and Malik [30], cannot effectively distinguish between edges and ramps [31], causing the staircase effect when smoothing ramp regions. Thus, variations of the traditional diffusion filtering (PM) have been proposed to overcome its shortcoming of staircase effect, such as quaternion diffusion [32] and nonlinear complex diffusion [33]. However, most of them worked in a simple scale thus may have the limitation of retaining subtle features. This leads us to consider diffusion in multi-scale.

In this paper, we propose a multi-scale modified nonlinear diffusion (MND) filter using Laplacian pyramid decomposition for CBCT with low radiation dose, with the aim of improving the quality of intra-operative CBCT in EVAR procedures, so that the denoised result is good enough for image interpretation or further processing (e.g. registration with pre-operative CT or post-operative CT). In the modified diffusion filter, the diffusion function is constructed using the gradient magnitude image, rather than gradient values in neighborhood, so that edges can be well preserved and smoothed without introducing obvious staircase effect. And then the proposed Laplacian pyramid-based nonlinear diffusion (LPMND) can reduce noise and streak artifacts while preserving edges and detailed features and completely eliminate the staircase effect.

Methods

Laplacian pyramid-based modified nonlinear diffusion

The proposed Laplacian pyramid-based modified nonlinear diffusion (LPMND) method is illustrated in Fig. 1. The Laplacian Pyramid is a decomposition of the original image into a hierarchy of images so that each level corresponds to a different band of image frequencies [34]. More details are reported in Additional file 1: Appendix A1. In our method, the degraded image is decomposed into 3 levels with the approximation image as the highest level.Actually noise and useful signal components of an image will be reflected in different levels after decomposition using the Laplacian pyramid. Regard to a degraded image, as noise has high frequency, it mainly exists in the lower pyramid level. However, we need notice that although noise mainly exists in lower level, some also exists in higher level and should be discarded. Similarly, some useful information of structure may exist in the lower level and should be retained. The proposed modified nonlinear diffusion (MND) is therefore applied on the image in each level. Obviously, there are three steps as shown in Fig. 1: (1) transforming the image to be processed into its pyramid domain, (2) restoring all pyramid images by performing MND filter, and (3) reconstructing the Laplacian pyramid using the processed pyramid images in each level.

Fig. 1
figure 1

Outline of the proposed LPMND method

The proposed modified nonlinear diffusion

Although the anisotropic diffusion proposed by Perona and Malik (PM diffusion), more details about it is in Additional file 1: Appendix A2, has been proved to be useful in removing noise (except salt-pepper noise) in homogeneous region, it limits the smoothing at the edge pixels due to permitting more diffusion along the edge than across it and easily generates staircase effect due to the incorrect diffusion at ramp. The reason is that the gradient operator is not a proper measure to detect the ramp features (endpoints) and less effective in noise reduction within the ramp-edge. Several authors suggested using the second derivative to replace the gradient, so that edges and ramps can be effectively distinguished [35, 36]. Gilboa et al. [33] also pointed that the second derivative was a more suitable choice than the gradient because it has a high magnitude near the endpoints and low magnitude elsewhere. That’s, in the diffusion process, the diffusion coefficient should be small near the endpoints and large within the ramp so that noise over the ramp can be reduced and ramp edge can be preserved at the endpoints. However, using the second derivative as the edge indicator may introduce a numerical problem when third order derivatives are computed.

According to the above-mentioned principle, we modified the diffusion function to remove noise while preserving ramp edges. The modified diffusion coefficient c is a function of the gradient magnitude image which does not involve the second derivative. It is defined as

$$ c\left(g(I)\right)=\frac{1}{1+{\left(g(I)/k\right)}^2} $$
(1)

where I is the original image and g(I) is its gradient magnitude image and given by

$$ g(I)=\sqrt{{G_x}^2+{G_y}^2} $$
(2)

Gx = (I(i, j + 1) − I(i, j − 1))/2, Gy = (I(i + 1, j) − I(i − 1, j))/2.

In this paper, we adopt the explicit numerical scheme, and then the modified nonlinear diffusion (MND) process is

$$ {I}_{i,j}^{\left(n+1\right)}={I}_{i,j}^{(n)}+\Delta {t}^{(n)}{\sum}_{Z\in \left\{N,S,W,E\right\}}c\left({g}_Z^{(n)}\right)\kern0.5em {\nabla}_Z{I}_{i,j}^{(n)} $$
(3)

Where Δt is the time step and \( {\nabla}_Z{I}_{i,j}^{(n)} \) means directional derivative. This modified diffusion function not only involves the gradient component of four directions, but also involves far pixels when the corresponding neighboring pixel gzin the gradient magnitude image is computed.

Figure 2 illustrates the difference between PM model and MND model. Let the green point be the central pixel to be processed. In PM model, it involves four neighboring pixels (blue points) when calculating the directional diffusion coefficients. In comparison, MND model involves another 8 points (orange points) to calculate the corresponding four neighboring pixels in the gradient magnitude image. The MND model therefore takes more information into consideration when estimating the central pixel. Besides, MND diffusion function only refers to the first derivative, avoiding numerical problems when using the second derivative as the edge indicator. In addition, it is simple to be performed even if it involves more points.

Fig. 2
figure 2

The nonlinear diffusion models. a PM model, b MND model

Let’s consider the example shown in Fig. 3. Figure 3a shows a sample of ramp edge, Fig. 3b shows the calculated gradients in neighborhood, and Fig. 3c shows the corresponding value in the gradient magnitude version by convolved with Gx. With respect to the ramp point marked by A, with value 5, the diffusion coefficient in each direction (just West and East in this case) by PM model is cPM = 1/(1 + (3/k)2), whereas the corresponding diffusion coefficient in each direction is cMND = 1/(1 + (1.5/k)2) by (1). cMND > cPM, therefore, MND model has a larger flux at the ramp point than PM model as they have the same difference in gray. That’s to say, it has a stronger diffusion to reduce the noise across the ramp. Endpoints of the ramp were marked by B and C, with value 2 and 8, respectively. For endpoint B, the diffusion coefficients in the two models are both cW = 1 and cE = 1/(1 + (3/k)2), thus they have the same flux. The same result can be observed at the endpoint C. It is obvious that the two models have the same diffusion strength at endpoints, indicating that MND model has the same capability of preserving endpoints as that in PM model in this example, but has a stronger capability of reducing noise over the ramp. In Fig. 4, we plotted the filtered results of a one-dimensional ramp signal by PM model and MND model. Figure 4a shows the original ramp signal and noisy version contaminated by white Gaussian with SNR = 15.7 dB. Figure 4b shows the denoised results from the PM process and the proposed MND process with 10, 50, and 100 iterations. We observe that in the result from MND process less noise appears with less iterations (10 iterations) and the denoised ramp still keeps its shape with more iterations (100 iterations). This illustrates that ramp edges can be effectively preserved while noise is effectively removed by the MND process in comparison with PM process.

Fig. 3
figure 3

The response of operators to a vertical ramp edge. a Ramp edge, b response of the gradient operator, c response of the gradient magnitude operator

Fig. 4
figure 4

Performs of PM and the MND processes on a ramp edge. a Original (top) and noisy ramp signal (bottom), b denoised signal from PM process (blue line) and the proposed MND process (red line) with 10, 50, 100 iterations from top to bottom, respectively

Additionally, Bernardes et al. [37] pointed that using a Gaussian filter for c was beneficial to remove speckle noise. We also used a Gaussian filter to smooth the diffusion coefficient because noise in CBCT is complex and may contain speckle noise [38]. Thus, the diffusion coefficient in (1) is redefined as

$$ c\left(g(I)\right)={\frac{1}{1+{\left(g(I)/k\right)}^2}}^{\ast }{G}_{\sigma } $$
(4)

where is the convolution operator and Gσ is a Gaussian kernel of 3 × 3 size and standard deviation σ. When σ is small, the filtering of diffusion function in (4) does not change c dramatically at edges, whereas it turns diffusion less conservative at speckle points.

Parameters determination

Since the selection of threshold k is very important in the diffusion process, we prefer an automatic k estimated using the median absolute deviation (MAD) estimator proposed by Donoho [38].

$$ k=\frac{MAD}{0.6745} $$
(5)

where MAD is the median absolute deviation of the wavelet coefficients at the nest level in which most noise exist.

The time step Δt is another important parameter to control the speed of diffusion. It is usually set to a constant that closes to the time step limit of the convergence of the iterative update process. For explicit 2D schemes the maximum time step to achieve stability of the iterative update is 0.25 s. In this work, we adopted an adaptive time step in (6) so that it is small at the initial iterations in which higher values of c can be found due to noise and increases until to a steady condition, in which changes of c over time are small.

$$ \varDelta {t}^{(n)}=\frac{1}{4}\left[a+b\exp \left\{-\left|\mathrm{sum}\left(\raisebox{1ex}{$\partial {I}^{(n)}$}\!\left/ \!\raisebox{-1ex}{$\partial t$}\right.\right)/\mathrm{sum}\left({I}^{(n)}\right)\right|\right\}\right] $$
(6)

where \( \mathrm{sum}\left(\raisebox{1ex}{$\partial {I}^{(n)}$}\!\left/ \!\raisebox{-1ex}{$\partial t$}\right.\right)/\mathrm{sum}\left({I}^{(n)}\right) \)is the fraction of change of the image at iteration n, parameters a and b control the time step with a + b ≤ 1. We set a = 0.25 and b = 0.75 in this study. An example evolution of Δt(n) over the iterative process is shown in Fig. 5.

Fig. 5
figure 5

An evolution of the time step using (15)

The stopping criterion is one challenge of diffusion filtering and simply can be stopped manually by setting a fixed number of iterations. However, in the pyramid domain, pyramid images contain different noise levels. Thus, it is difficult to assign a unified value for each level. In this study, we don’t further discuss the stopping criterion, but stop the diffusion using the mean absolute error (MAE) between two adjacent diffusions

$$ MAE\left({I}^{(n)}\right)=\sum \limits_{i=1,j=1}^{i=M,j=N}\left|{I}_{i,j}^{(n)}-{I}_{i,j}^{\left(n-1\right)}\right|/\left(M\times N\right) $$
(7)

where \( {I}_{i,j}^{\left(n-1\right)} \) and \( {I}_{i,j}^{(n)} \) are the filtered value at pixel at (n-1)th and (n)th iterations, respectively, and M and N are height and width of the processed image in pyramid domain, respectively. For each level image, the diffusion process stops automatically when the value of MAE is smaller than a preset threshold.

Figure 6 shows one slice of CBCT with poor quality in the left column and its 3 levels of Laplacian pyramid in the right top row. We can observe from the right top row that noise and artifacts play a dominant role in the lowest level and much decrease with the level increases, while progressively coarser features such as edges and structures are prominent in images of higher level. It is obvious that the distributions of noise and artifacts are very different in each level. We therefore should consider different parameters including parameterk, time step Δt, standard deviation σ in (4) and the number of iterations N when applying the MND filter on each pyramid level.

Fig. 6
figure 6

Noise and artifacts suppression in 3 levels of Laplacian pyramid for the degraded CBCT image. Left column shows the original image; right top row shows its 3 levels of Laplacian pyramid; right bottom row shows the restored Laplacian pyramid levels by MND filter

The parameter k and time step Δt could be adaptively determined for each pyramid image still using (5) and (6). With regard to σ in (4), it changes in different levels. Therefore, we should do two things: (i) set a high σ for the lowest level, since it contains most noise and artifacts, and (ii) reduce σ when the level increases and set a small value in the highest level so that it will not change c significantly.

The right bottom row in Fig. 6 shows the processed images by the proposed MND filter in each pyramid level. σ from level L0 to L2was set to 0.8, 0.5, and 0.1, respectively. We can see that the MND filter leads to artifacts suppression effectively in the lower pyramid levels, particularly in the lowest level, yet modifies little in the coarsest level.

CBCT simulation and acquisition

In order to evaluate the performance of the proposed LPMND, we performed experiments on simulated and real abdominal CBCT images. The simulation experiment was conducted by using a numerical phantom (shown in Fig. 7a) which imitates the structures of a real abdominal CBCT image, and thus we can evaluate the proposed method quantitatively. A large ellipse object in the center of the phantom was used as reference “background tissue”. Object A was used to mimic bony structures and elliptical objects B, C, D were used to mimic soft tissues. Object E in the center simulates the guidewire inserted into the artery and was also chosen for a spatial resolution test. This phantom was simulated for the central slice in the CBCT imaging.

Fig. 7
figure 7

CBCT images reconstructed by FDK for digital phantom. a Clean CBCT image (as the reference image), b noisy CBCT image

This phantom was then back-projected according to the following scanning: The distance of X-ray source to the center of rotation was 541 and 984 angular samples evenly spanned on a circular orbit of 360°. The central slice was then reconstructed from this noise-free sinogram using the well-known FDK algorithm [39], as shown in Fig. 7a, and can be as the reference image in the performance analysis. To simulate a FDK reconstructed CBCT image acquired by a low-dose protocol, we generated a noisy sinogram by adding signal-dependent Gaussian noise to the noise-free sinogram, according to the noise model in (8). This noise model was presented by Li et al... [40] and Wang et al [41] with the content that the projection data after system calibration and logarithm transformation was approximately Gaussian distributed with a non-linear dependency between sample mean and variance as follows.

$$ {\sigma^2}_{p_s}=f\exp \left({\overline{p}}_s/\eta \right) $$
(8)

where \( {\overline{p}}_s \) and \( {\sigma}_{p_s}^2 \) is the mean and variance at detector bin s. η and f are parameters determined by system settings. In our study, η and f were set to 22,000 and 200, respectively. The CBCT image was then reconstructed by FDK, shown in Fig. 7(b). Note that the simulated image contains obvious noise and streak artifacts that is the characteristic of low-dose CBCT.

For experiments using patient data, the study was approved by the ethics committee of Rennes University Hospital (France) on April 10, 2016. Patient informed consent was obtained for being registered anonymously in the database. We used a Siemens Artis zeego C-arm system in the operation room to obtain a CBCT for a patient during an EVAR. Four protocols are preset in this system: they are 5sDR protocol, 5sDSA protocol, 8sDR protocol and 20sDR protocol. Low-dose intra-operative CBCT was realized using a 5sDR protocol. The 5sDR protocol has 5-s acquisition time capturing 133 frames at 30 frames/second (f/s). The 5sDSA protocol is used for acquiring a digital subtraction angiography (DSA). The 5sDSA protocol has the same characteristics of the 5sDR but it realized two rotations, the first one is used to acquire the mask, and the second one is synchronized with the injection with the injection of contrast medium. Other acquisition parameters were: the distance from the X-ray source to detector is 1199 mm, the field of view (FOV) is 100 mm × 100 mm, the slice thickness is 0.4804 mm, the size of one slice of the CBCT is 512 × 512 and pixel spacing is 0.4804 mm × 0.4804 mm. The related scanning parameters were 91KV, 243mAs.

Performance evaluation

To evaluate the performance of noise reduction, we computed several well-known image performance metrics: mean square error (MSE), peak signal-to-noise (PSNR), and correlation coefficient (CC). MSE and the PSNR are error metrics used to test the quality of filtered image. MSE represents the cumulative squared error between the filtered and the original image, while PSNR represents a measure of the peak error. CC indicates the degree of spatial similarity between the original and the filtered image. Definitions of these metrics are in (9)–(11).

$$ MSE=\sum \limits_{i=1,j=1}^{i=M,j=N}{\left({\widehat{I}}_{i,j}-{I}_{i,j}\right)}^2/\left(M\times N\right) $$
(9)
$$ PSNR=20{\log}_{10}\left(\frac{2^B-1}{\sqrt{MSE}}\right) $$
(10)
$$ CC=\sum \limits_{i=1,j=1}^{i=M,j=N}\left({\widehat{I}}_{i,j}-{\mu}_{\widehat{I}}\right)\left({I}_{i,j}-{\mu}_I\right)/\sqrt{\sum \limits_{i=1,j=1}^{i=M,j=N}{\left({\widehat{I}}_{i,j}-{\mu}_{\widehat{I}}\right)}^2{\left({I}_{i,j}-{\mu}_I\right)}^2} $$
(11)

where I and μI are the original image and its mean value, respectively. \( \widehat{I} \) and \( {\mu}_{\widehat{I}} \) are the filtered image and its mean value, respectively. B represents the bits per sample. A good result should have small MSE, high PSNR and CC.

Besides, to quantify localized differences, we evaluated the contrast-to-noise ratio (CNR) over specific regions of interest (ROIs). The CNR measures the contrast between a ROI and the background region. It is generally defined as

$$ CNR=\left|{\mu}_R-{\mu}_B\right|/\sqrt{\sigma_R^2+{\sigma}_B^2} $$
(12)

where μR and σR are the mean and standard deviation of the ROI, including the homogeneous region and μB and σB are the mean and standard deviation of the background region, respectively. To evaluate the performance of resolution preservation, we also analyzed the spatial resolution by using the modulation transfer function (MTF), which was approximately obtained by using the object E in Fig. 7(a). A ROI placed at the center of the object was extracted and the line spread function (LSF) then was obtained by integrating the ROI in y-direction, finally MTF can be obtained by calculating the Fourier transform of the normalized and zero-padded LSF. In addition, the sharpness measurement tool proposed by Taubmann et al. [42] was used to support comparison of edge sharpness in images. This tool not only gives a plot of sharpness estimates along an edge for images to be tested, but also computed sharpness values.

Results

In this section, we evaluate the proposed LPMND filter by using both simulated image and real abdominal CBCT data. The simulation study gives quantitative performance analysis and the patient study demonstrates the practical applicability of the proposed filter in the procedure of EVARs. In each study, we compared the LPMND filter with other filters including PM filter with diffusion function in (1), nonlinear complex diffusion filter (NCDF), bilateral filter (BF) and nonlocal means filter (NLMF). For LPMND, αin Gaussian kernel for REDUCE and EXPAND operators was set to 0.375 and the images were decomposed into three pyramid levels as more levels didn’t guarantee improved performances. “db2” wavelet was used to calculate the parameter k in (5). The standard deviation σ of Gaussian filter in (4) was set to 0.8, 0.5, and 0.1 from the lowest level to the highest level, respectively. These parameters were hand-tuned set to obtain the possible optimization based on trial errors.

Results from simulated image

The five different filters were applied on the simulated image in Fig. 7b, in which plenty of noise and artifacts exist. The time step was 0.15 and 0.055 for PM and NCDF, respectively. The parameter k of PM and NCDF was also determined by (5). LPMND, PM and NCDF followed the same stop criterion, i.e., MAE criterion in (7), and the thresholds for the MAE in PM, NCDF and LPMND were all set to 0.11. The search window and similarity window in NLMF ware set to 21 × 21 and 7 × 7, and the filtering parameter was 18, whereas in the BF, the standard deviations of closeness function and similarity function were 3 and 100, respectively. All of the parameters in each filter were appropriately set to achieve a matched noise level. PSNR values in the filtered images from PM, NCDF, BF, NLMF, and LPMND were 79.98, 78.89, 80.37, 80.97, and 81.92, respectively. Figure 8 shows the comparison of different filters. PM doesn’t give a satisfactory result due to the fact that obvious staircase effect exists in the filtered image, although it has a good capability of preserving edge. NCDF behaves well on avoiding staircase effect but is less satisfactory on suppressing enough noise and streak artifacts. Like NCDF, BF don’t perform well on suppressing noise and streak artifacts, while NLMF generates wavy artifacts in the filtered image, as shown in Fig. 8e. In comparison, the proposed LPMND gives a better tradeoff between noise/artifacts reduction and edge preservation. From Fig. 8f, we can observe that weak edges are well enhanced.

Fig. 8
figure 8

Comparison of the simulated image. a Simulated low-dose CBCT image, b-f images filtered by PM, NCDF, BF, NLMF, and LPMND, respectively

To quantitatively compare the performances of different filters, we then calculated the quality parameters MSE and CC at the matched PSNR, as shown in Table 1. Furthermore, we calculated the CNR and ENR to quantify the local performance in regions of interest. ROI1 indicated by a rectangle in Fig. 8a was seen as the background region, and three homogeneous regions (ROI2, ROI3, and ROI4) were used to calculate CNR and ENR values. Table 2 lists the CNR value for each region, respectively. It is obvious from Table 1 that LPMND has the lowest value of MSE and highest values of CC in comparison with the other filters, meaning it is most closed to the original image. Furthermore, it has highest values of CNR for each homogeneous region, as shown in Table 2.

Table 1 MSE, PSNR, and CORR values of different filtered images
Table 2 CNR values of three ROIs

To measure the resolution of the filtered images, ROI5 indicated in Fig. 8a was used to calculate the MTF. Figure 9 depicts the MTF curves of the original image and the filtered images in Fig. 8b-f. We observe that the MTF curve of LPMND filtered image is most closed to that of the original image, compared with results from PM, NCDF, BF, and NLMF. It means that the spatial resolution in the LPMND filtered result is higher than those filtered by the other four filters. In order to evaluate the performances of edge preservation, we showed an enlarged low-contrast region (ROI4) and analyzed its edge by using the sharpness measurement tool. The corresponding enlarged images related to Fig. 8b-f are shown in Fig. 10b-f with the enlarged one of original image shown in Fig. 10a. The tested edge is indicated by the red curve in Fig. 10a. Figure 10g plots the sharpness estimates along the edge for different filtered image. The computed edge sharpness estimates are summarized in Table 3. In addition, Fig. 10h shows the profiles passing ROI4 indicated by a blue line in Fig. 10a. Through those profiles, it can be observed that NLMF failed on preserving edges of the low-contrast region; PM, BF and NCDF don’t have good behaviors of smoothing the homogeneous region, seeing the part marked by a red circle; By contrast, LPMND likely reaches a compromise between smoothing and edge preserving.

Fig. 9
figure 9

Comparison of MTF curves from the original image and the filtered images in Fig. 8b-f

Fig. 10
figure 10

Comparison of sharpness for a low-contrast edge. a-f Enlarged ROI4 of the original image and the filtered image by PM, NCDF, BF, NLMF, and LPMND, respectively. g The plot shows the corresponding sharpness estimates along the edge (indicated by red curve) for different filtered image. h Comparison of horizontal profiles, indicated by the blue line in Fig. 10a

Table 3 Edge sharpness of low-contrast edge in different filtered images

Results from real CBCT images

In this section, we used the intra-operative CBCT data of an AAA (abdominal aortic aneurysm) patient to evaluate the effectiveness of the proposed filter. The original image of slice 152 is shown in Fig. 11(a1) and the filtered results obtained by PM, NCDF, BF, NLMF and LPMND are shown in Fig. 11(b1-f1), respectively. The thresholds for the MAE in PM, NCDF and LPMND were all set to 2. The time step was 0.24 and 0.1 for PM and NCDF, respectively. The search window and similarity window in NLMF were set to 7 × 7 and 3 × 3, and the filtering parameter was 150, whereas in the BF, the standard deviations of closeness function and similarity function were 3 and 500, respectively. All of parameters in each filter were hand-tuned set to obtain possible best results based on trial errors. For a better illustration, the enlarged regions of bony structure and renal structure, marked by a white block in Fig. 11 (a1), are illustrated in Fig. 11(a2-f2). We can see that the filtered images obtained by PM, BF and NLMF have sharp edges, but PM yields staircase effect while BF and NLMF (in particular) are less helpful for artifacts suppression. NCDF behaves well on smoothing homogeneous region, but makes the bony structures blurred. In comparison, LPMND outperforms other filters on suppressing artifacts and preserving edges. Furthermore, a vertical profile (indicated by a red line in Fig. 11(a1) is plotted in Fig. 12 for the original image and each filtered result. The filtered profiles are overlapped by the original profile for comparison. Through these profiles, we can observe that LPMND filtered image shows better smoothness in homogeneous regions.

Fig. 11
figure 11

One slice of the CBCT data from one patient and its enlarged region which is indicatd by a white square in (a1). a1 the original image, b1-f1 filtered image by PM, NCDF, BF, NLMF and LPMND, respectively. a2-f2 Corresponding enlargements in (a1-f1). Display window: [− 1024 HU, 1976 HU]

Fig. 12
figure 12

Comparison of vertical profiles, indicated by the red line in Fig. 11 (a1), for different filtered images as shown in Fig. 11 (b1-f1). (a) Original image, (b)-(f) images filtered by PM, NCDF, BF, NLMF, and LPMND, respectively

In practice, the noise level changes across slices even in the same CBCT data, we therefore selected another slice (slice 29) that contains more artifacts to test the proposed filter. Figure 13 shows the comparison of processed images by different filters. In order to evaluate the consistent performance of different filters, the parameters in each filter were kept the same as those used for slice 152. It can be observed that PM, NLMF and BF cannot suppress artifacts effectively compared with NCDF and LPMND. On the other hand, the bony structures in the NCDF filtered image are a little blurred, whereas the LPMND filtered image has clear bony structures. Thus, the LPMND filter may have superiority for practical applications.

Fig. 13
figure 13

Comparison of slice 29. Display window: [− 1024 HU, 1976 HU]. (a) Original image, (b)-(f) images filtered by PM, NCDF, BF, NLMF, and LPMND, respectively

We also calculated the average time for one CBCT data of 365 slices under the environment that OS: 64-bit; Windows 7; CPU: Intel Core(TM) i7, 8G RAM. The computation time for PM, NCDF, BF, NLMF and the proposed LPMND is about 6 min, 43 min, 19 min, 59 min, and 4 min, respectively.

Discussion

The number of endovascular aneurysm repair is yearly growing and concerns about radiation are more topical than ever. Works dealing with low-dose X-ray imaging protocols are appropriate for optimization of the clinical use of this acquisition. Therefore, noise and artifacts caused by low dose protocols in intra-operative CBCT acquisition need to be addressed. Our study aimed to improve the quality of intra-operative CBCT acquired with low dose protocols so that a more accurate segmentation or registration can be carried out, thus to guide and monitor the insertion procedure. Furthermore a post-procedural CBCT with good noise reduction could better help assess early complications compared to classical 2D angiography. In practice, we can restore the projection image or the reconstructed images, or obtain the final reconstructed image using a statistical image reconstruction (SIR) algorithm. Projection denoising takes noise properties in projections into account, yet has the potential disadvantage that the definition of edge in projection data is not definite, resulting in sharpness loss in image domain. SIR methods more focus on reconstruction rather than denoising techniques, utilizing the noise properties in projections and image prior information. While more sophisticated physical modeling and image reconstruction can be viewed as a superior approach, they always have high computation load and are highly dependent on special scanner model, i.e., requiring more detailed information such as scanning geometry, correction physics, etc. This limitation appeals a more broadly used denoising method that can perform on different systems, and leads us to think more about denoising after reconstruction. Image post-processing techniques, working in the image space alone, are retrospectively applied and relatively simple to implement. For the reasons mentioned above, we choose to implement noise and artifacts suppression in image space.

Laplacian pyramid is a very useful tool for multi-scale analysis. Unlike sub-band decomposition in wavelet transform, REDUCE and EXPAND operators generate a series of low-pass and band-pass images that contain approximation and detailed images without requiring quadrature mirror filters. The pyramid offers a useful image representation for denoising tasks, since it can represent information that is localized in both space and spatial frequency in each level [43]. Since CBCT noise in image space is difficult to model accurately and has strong correlations, the Laplacian pyramid is considered to remove much of the pixel-to-pixel correlation.

In this study, we have proposed a modified nonlinear diffusion filter in Laplacian pyramid to reduce noise and artifacts without damaging image features. It consists of the following three steps: (1) Laplacian pyramid decomposition, (2) performance of the modified nonlinear diffusion filter in each pyramid level, and (3) reconstruction from the processed pyramid images. We modified the diffusion function of PM model so that it can remove noise across ramp edges, to a large extent avoiding transforming ramp edges in to stairs. The staircase effect is then remedied by multi-scale decomposition through Laplacian pyramid. There are several advantages of this proposed method beyond that. First, the threshold parameter k in diffusion function is adaptive and automatically determined by (5). Secondly, we adopted an adaptive time step which changes according to iterative image changes. And then the speed of diffusion process is accelerated as well. Thirdly, the diffusion is automatically stopped using the MAE between two adjacent diffusions, rather than manually stopped by setting a fixed number of iterations. Furthermore, the value of MAE for different level image in pyramid domain could be fixed the same according to our experiments, avoiding multi-choice in multi-scale decomposition. In addition, in view of the computation cost, the proposed method is fast and easy to implement, thus is possible to be used for clinical applications.

In this work, our effort was focused on the noise and streak artifacts suppression for the intra-operative CBCT during an EVAR. We have performed simulated and real CBCT images to evaluate the proposed method. Noise and artifacts lead to a degraded image, especially in low-contrast regions. The proposed method can suppress noise and streak artifacts meanwhile protect low-contrast regions. We also compared the proposed method with other edge-preserving filters through quantitative analysis in simulation study and qualitative analysis in patient study. The MTF study shows that LPMND can produce good image resolution (Fig. 9), and the sharpness estimation shows that it works well on preserving edges of low-contrast regions (Fig. 10 and Table 3). Actually, we have verified that the proposed method performed well on clinical CBCT data with full slices. This method could be used to improve results of further processing (e.g. CT – CBCT registration).

There are also some drawbacks of this proposed method. One is the selection of parameter k. Actually, k plays an important role in diffusion procedures and should correspond to the noise level. However, the noise level in CBCT is very complex, resulting in the difficulty to establish a noise model. The performance of the proposed filter may be further improved with an accurate noise model to determine k. Furthermore, although the proposed method performs well for the streak artifacts caused by noise, it is not very effective for the metal artifacts caused by high-attenuation objects, e.g., the inserted guidewire. In fact, metal artifacts remove is another challenge for CBCT imaging. Finding an effective solution to suppress metal artifacts is necessary and is considered in our next step. In addition, the proposed method is presented on 2D case, making the information from the successive slices not to be concerned. Although some other methods were proposed in three dimension, such as the KL-PWLS method in [20], the computational load is an issue to be resolved. We therefore processed the CBCT data slice by slice, and in fact the proposed method was satisfactory when we used it in practical applications.

Conclusion

We have presented a Laplacian pyramid-based nonlinear diffusion filter for intra-operative CBCT in endovascular aneurysm repair. The multi-scale diffusion can effectively remove noise and artifacts, improve the CNR for low-contrast regions, and meanwhile preserve edges and detailed features. Compared with some other edge-preserving filters, the proposed method has shown better performance for the CBCT in EVAR procedures.

Abbreviations

AAA:

Abdominal aortic aneurysm

BF:

Bilateral filter

CBCT:

Cone-beam computed tomography

CC:

Correlation coefficient

CNR:

Contrast-to-noise ratio

DSA:

Digital subtraction angiography

EVAR:

Endovascular aneurysm repair

FOV:

Field of view

LPMND:

Laplacian pyramid-based modified nonlinear diffusion

LSF:

Line spread function

MAD:

Median absolute deviation

MND:

Modified nonlinear diffusion

MSE:

Mean square error

MTF:

Modulation transfer function

NCDF:

Nonlinear complex diffusion filter

NLMF:

Nonlocal means filter

PSNR:

Peak signal-to-noise ratio

ROI:

Region of interest

SIR:

Statistical image reconstruction

References

  1. Herana NS, Songa JK, Nambaa K, Smitha W, Niimia Y, Berensteina A. The utility of DynaCT in neuroendovascular procedures. AJNR Am J Neuroradiol. 2006;27:330–2.

    Google Scholar 

  2. Biasi L, Ali T, Arumagum LA, Loftus I, Morgan R, Thompson M. Intra-operative DynaCT improves technical success of endovascular repair of abdominal aortic aneurysms. J Vasc Surg. 2009;49:288–95.

    Article  PubMed  Google Scholar 

  3. Biasi L, Ali T, Hinchliffe LR, Morgan R, Loftus I, Thompson M. Intraoperative DynaCT detection and immediate correction of a type la endoleak following endovascular repair of abdominal aortic aneurysm. Cardiovasc Intervent Radiol. 2009;32:535–8.

    Article  PubMed  Google Scholar 

  4. Dijkstra ML, Eagleton MJ, Greenberg RK, Mastracci T, Hernandez A. Intraoperative C-arm cone-beam computed tomography in fenestrated /branched aortic endografting. J Vasc Surg. 2011;53:583–90.

    Article  PubMed  Google Scholar 

  5. Törnqvist P, Dias N, Sonesson B, Kristmundsson T, Resch T. Intra-operative cone beam computed tomography can help avoid Reinterventions and reduce CT follow up after Infrarenal EVAR. Eur J Vasc Endovasc Surg. 2015;49:390–5.

    Article  PubMed  Google Scholar 

  6. Tacher V, Lin M, Desgranges P, Deux JF, Grünhagen T, Becquemin JP, Luciani A, Rahmouni A, Kobeiter H. Image guidance for endovascular repair of complex aortic aneurysms: comparison of two-dimensional and three-dimensional angiography and image fusion. JVIR. 2013;24:1698–706.

    Article  PubMed  Google Scholar 

  7. Kobeiter H, Nahum J, Becquemin JP. Zero-contrast thoracic endovascular aortic repair using image fusion. Circulation. 2011;124:e280–2.

    Article  PubMed  Google Scholar 

  8. Hertault A, Maurel B, Sobocinski J, Gonzalez TM, Roux ML, Azzaoui R, Midull M, Haulon S. Impact of hybrid rooms with image fusion on radiation exposure during endovascular aortic repair. Eur J Vasc Endovasc Surg. 2014;48:382–90.

    Article  PubMed  CAS  Google Scholar 

  9. Sailer AM, de Haan MW, Peppelenbosch AG, Jacobs MJ, Wildberger JE, Schurink GWH. CTA with fluoroscopy image fusion guidance in endovascular complex aortic aneurysm repair. Eur J Vasc Endovasc Surg. 2014;47:349–56.

    Article  PubMed  CAS  Google Scholar 

  10. McNally MM, Scali ST, Feezor RJ, Neal D, Huber TS, Beck AW. Three-dimensional fusion computed tomography decreases radiation exposure, procedure time, and contrast use during fenestrated endovascular aortic repair. J Vasc Surg. 2015;61:309–16.

    Article  PubMed  Google Scholar 

  11. Kaladji A, Dumenil A, Castro M, Haigron P, Heautot JF, Haulon S. Endovascular aortic repair of a post-dissecting Thoracoabdominal aneurysm using intraoperative fusion imaging. J Vasc Surg. 2013;57:1109–12.

    Article  PubMed  Google Scholar 

  12. Hertault A, Maurel B, Midulla M, Bordier C, Desponds L, Kilani MS, Sobocinski J, Haulon S. Minimizing radiation exposure during endovascular procedures: basic knowledge, literature review, and Reporting Standards. Eur J Vasc Endovasc Surg. 2015;50:21–36.

    Article  PubMed  CAS  Google Scholar 

  13. Islam MK, Purdie TG, Norrlinger BD, Alasti H, Moseley DJ, Sharpe MB, Siewerdsen JH, Jaffray DA. Patient dose from kilovoltage cone beam computed tomography imaging in radiation therapy. Med Phys. 2006;33:1573–82.

    Article  PubMed  Google Scholar 

  14. Wen N, Guan H, Hammoud R, Pradhan D, Nurushev T, Li S, Movsas B. Dose delivered from Varian's CBCT to patients receiving IMRT for prostate cancer. Phys Med Biol. 2007;52:2267–76.

    Article  PubMed  Google Scholar 

  15. Wang J, Li T, Xing L. Iterative image reconstruction for CBCT using edge-preserving prior. Med Phys. 2009;36:252–60.

    Article  PubMed  Google Scholar 

  16. Lee RD. Common image artifacts in cone beam CT. In: Summer 2008 AADMRT Newsletter. American Association of Dental Maxillofacial Radiographic Technicians. 2008.  http://www.aadmrt.com/article-1---2008.html. Accessed 19 Aug 2018.

  17. Kachelrieß M, Watzke O, Kalender WA. Generalized multi-dimensional adaptive filtering for conventional and spiral single-slice, multi-slice, and cone-beam CT. Med Phys. 2001;28:475–90.

    Article  PubMed  Google Scholar 

  18. Manduca A, Yu L, Trzasko JD, Khaylova N, Kofler JM, McCollough CM, Fletcher JG. Projection space denoising with bilateral filtering and CT noise modeling for dose reduction in CT. Med Phys. 2009;36:4911–9.

    Article  PubMed  PubMed Central  Google Scholar 

  19. Yu L, Manduca A, Trzasko J, Khaylova N, Kofler JM, McCollough CH, Fletcher JG. Sinogram smoothing with bilateral filtering for low-dose CT. Proc SPIE. 2008;6913:691329.

    Article  CAS  Google Scholar 

  20. Ouyang L, Solberg T, Wang J. Noise reduction in low-dose cone beam CT by incorporating prior volumetric image information. Med Phys. 2012;39:2569–77.

    Article  PubMed  Google Scholar 

  21. Sidky EY, Pan X. Image reconstruction in circular cone-beam computed tomography by constrained, total-variation minimization. Phys Med Biol. 2008;53:4777–807.

    Article  PubMed  PubMed Central  Google Scholar 

  22. Li Z, Yu L, Trzasko J, Fletcher J, McCollough C, Manduca A. Adaptive nonlocal means filtering based on local noise level for CT denoising. Med Phys. 2014;41:011908.

    Article  PubMed  Google Scholar 

  23. Barash D, Comaniciu D. A common framework for nonlinear diffusion, adaptive smoothing, bilateral filtering and mean shift. Image Vison Comput. 2004;22:73–81.

    Article  Google Scholar 

  24. Chen Y, Chen W, Yin X, Ye X, Bao X, Luo L, Feng Q, Li Y, Yu X. Improving low-dose abdominal CT images by weighted intensity averaging over large-scale neighborhoods. Eur J Radiol. 2011;80:e42–9.

    Article  PubMed  Google Scholar 

  25. Buades A, Coll B, Morel JM. A review of image denoising algorithms, with a new one. Multiscale Model Simul. 2005;4:490–530.

    Article  Google Scholar 

  26. Kroon DJ, Slump CH, Maal TJJ. Optimized anisotropic rotational invariant diffusion scheme on cone-beam CT. In: 13th International Conference on Medical Image Computing and Computer-Assisted Intervention (MICCAI), Springer LNCS 6363, Beijing, China; 2010. p. 221-228.

  27. Bai T, Yan H, Shi F, Jia X, Lou Y, Xu Q, Jiang S, Mou X. 3D dictionary learning based iterative cone beam CT reconstruction. Int J Cancer Ther Oncol. 2014;2:020240.

    Article  Google Scholar 

  28. Ghadrdan S, Alirezaie J, Dillenseger J, Babyn P. Low-dose computed tomography image denoising based on joint wavelet and sparse representation. Conf Proc IEEE Eng Med Biol Soc. 2014;2014:3325–8.

    PubMed  Google Scholar 

  29. Shi L, Chen Y, Shu H, Luo L, Toumoulin C, Coatrieux J. Low-dose CT image processing using artifact suppressed dictionary learning. In: 11th International Symposium on Biomedical Imaging (ISBI), IEEE, Beijing, China; 2014. p.1127–1130.

  30. Perona P, Malik J. Scale-space and edge detection using anisotropic diffusion. IEEE Trans Pattern Anal Mach Intell. 1990;12:629–39.

    Article  Google Scholar 

  31. Chen Q, Montesinos P, Sun QS. Ramp preserving PM model. Signal Process. 2010;90:1963–75.

    Article  Google Scholar 

  32. Liu ZX, Lian SG, Ren Z. Quaternion diffusion for color image filtering. J Comput Sci Technol. 2006;21:126–36.

    Article  Google Scholar 

  33. Gilboa G, Sochen N, Zeevi YY. Image enhancement and denoising by complex diffusion processes. IEEE Trans Pattern Anal Mach Intell. 2004;26:1020–36.

    Article  PubMed  Google Scholar 

  34. Burt PJ, Adelson EH. The Laplacian pyramid as a compact image code. IEEE Trans Commun. 1983;31:532–40.

    Article  Google Scholar 

  35. Chan T, Marquina A, Mulet P. High order total variation-based image restoration. SIAM J Sci Comput. 2000;22:503–16.

    Article  Google Scholar 

  36. You YL, Kaveh M. Fourth-order partial differential equations for noise removal. IEEE Trans Med Process. 2000;10:1723–30.

    Article  Google Scholar 

  37. Bernardes R, Maduro C, Serranho P, Araújo A, Barbeiro S, Cunha-Vaz J. Improved adaptive complex diffusion despeckling filter. Opt Express. 2010;18:24048–59.

    Article  PubMed  Google Scholar 

  38. Donoho DL. De-noising by soft-thresholding. IEEE Trans Inf Theory. 1995;4:613–27.

    Article  Google Scholar 

  39. Miao H, Zhao H, Gao F, Gong S. Implementation of FDK reconstruction algorithm in cone-beam CT based on the 3D Shepp-Logan model. In: 2nd International Conference on Biomedical Engineering and Informatics, IEEE, Tianjin, China; 2009. p. 1-5.

  40. Li T, Li X, Wang J, Wen J, Lu H, Hsieh J, Liang Z. Nonlinear sinogram smoothing for low-dose x-ray CT. IEEE Trans Nucl Sci. 2004;51:2505–13.

    Article  Google Scholar 

  41. Wang J, Lu H, Liang Z, Eremina D, Zhang G, Wang S, Chen J, Manzione J. An experimental study on the noise properties of x-ray CT sinogram data in radon space. Phys Med Biol. 2008;53:3327–41.

    Article  PubMed  PubMed Central  Google Scholar 

  42. Taubmann O, Wetzl J, Lauritsch G, Maier A, Hornegger J. Sharp as a tack: Measuring and Comparing Edge Sharpness in Motion-Compensated Medical Image Reconstruction. In: Bildverarbeitung für die Medizin; 2015. p. 425–30.

    Google Scholar 

  43. Adelson EH, Anderson CH, Bergen JR, Burt PJ, Ogden JM. Pyramid methods in image processing. RCA engineer. 1984;29:33–41.

    Google Scholar 

Download references

Funding

This work was partially conducted in the experimental platform TherA-Image (Rennes, France) supported by Europe FEDER. This work was partially supported by the French National Research Agency (ANR) in the context of Endosim project (grant n° ANR-13-TECS-0012), and in the framework of the Investissement d’Avenir Program through Labex CAMI (grant n° ANR-11-LABX-0004-01). The funders had no role in the design of the study, the collection, analysis and interpretation of the data, and in writing the manuscript.

Availability of data and materials

The datasets used and/or analysed during the current study are available from the corresponding author on reasonable request.

Author information

Authors and Affiliations

Authors

Contributions

YL designed the study and completed drafting of the manuscript. MC and AK acquired the data. MC completed the data collection. ML and AK provided medical guidance and interpreted and discussed the results. PH contributed to design the study and revised the manuscript critically for important intellectual content. All authors read and approved the final version of the manuscript.

Corresponding author

Correspondence to Pascal Haigron.

Ethics declarations

Ethics approval and consent to participate

The study was approved by the Ethics Committee of the University Hospital of Rennes on April 10, 2016 (decision n°16.38) and declared to the French data protection agency (CNIL) on February 19, 2016 (declaration n°1,932,135). All patients received an information leaflet and provided written consent.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Additional file

Additional file 1:

Appendix. The Laplacian Pyramid (A1), Anisotropic Nonlinear Diffusion (A2). (PDF 62 kb)

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Liu, Y., Castro, M., Lederlin, M. et al. An improved nonlinear diffusion in Laplacian pyramid domain for cone beam CT denoising during image-guided vascular intervention. BMC Med Imaging 18, 25 (2018). https://doi.org/10.1186/s12880-018-0269-1

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12880-018-0269-1

Keywords