Truncated total variation in fractional B-spline wavelet transform for micro-CT image denoising

Abstract

BACKGROUND:

In medical applications, computed tomography (CT) is widely used to evaluate various sample characteristics. However, image quality of CT reconstruction can be degraded due to artifacts.

OBJECTIVE:

To propose and test a truncated total variation (truncation TV) model to solve the problem of large penalties for the total variation (TV) model.

METHODS:

In this study, a truncated TV image denoising model in the fractional B-spline wavelet domain is developed to obtain the best solution. The method is validated by the analysis of CT reconstructed images of actual biological Pigeons samples. For this purpose, several indices including the peak signal-to-noise ratio (PSNR), structural similarity index (SSIM) and mean square error (MSE) are used to evaluate the quality of images.

RESULTS:

Comparing to the conventional truncated TV model that yields 22.55, 0.688 and 361.17 in PSNR, SSIM and MSE, respectively, using the proposed fractional B-spline-truncated TV model, the computed values of these evaluation indices change to 24.24, 0.898 and 244.98, respectively, indicating substantial reduction of image noise with higher PSNR and SSIM, and lower MSE.

CONCLUSIONS:

Study results demonstrate that compared with many classic image denoising methods, the new denoising algorithm proposed in this study can more effectively suppresses the reconstructed CT image artifacts while maintaining the detailed image structure.

Keywords

Digital image processing image denoising fractional B-spline wavelet truncated total variation (truncated TV) model Split-Bregman iteration

1 Introduction

Micro-computed tomography (micro-CT) is an important tool in biomedical research and preclinical applications, which can provide visual inspection and quantitative information of small animals and biological samples [1]. Currently, high-quality micro-CT images need to acquire projection data in many views, which limits the system throughput. In real application, especially the biomedical applications, to reduce radiation dose while maintaining the diagnostic performance is a major challenge. The approach to decrease radiation dose is to reduce the number of projections. However, the term of deformation or damage to the image sample that arise on account of scarcity of projections is usually called the aliasing distortions, which will create the visual noise in the reconstructed CT images [2]. However, the noise of CT images can be controlled by three main denoising: projection domain, iterative reconstruction and image denoising based on post-processing [3]. In this paper, we use post-processing methods to reduce the noise of CT images.

A traditional post-processing method to reduce noise from image is spatial filtering. Spatial filters reduce noise by applying filtering directly on original noisy CT image. In spatial domain, linear filters were adopted to denoise the images, but it was not successful for preserving edge over the images. Mean filtering was used to reduce the gaussian noise but for high noise it produces blurry images. To overcome that, Some of methods was further used, such as non-local means [4], total variation [5], and BM3D [7], while some use newer deep learning techniques [8]. The most investigated transform in denoising is the wavelet transform [9]. The basic idea of wavelet transform is to use a cluster of functions to represent or approximate a signal or function [10]. Different components in the image correspond to different frequency components. The frequency corresponding to the smooth part of the image is the low frequency part, and the frequency corresponding to the noise and detail part is the high frequency part [11]. While removing noise, image denoising preserves image details as much as possible, so it needs targeted selection when removing high frequency components [12]. The scale factor in the wavelet transform reflects the frequency information. By changing the wavelet coefficients at different scales, the high frequency is suppressed and the low frequency is retained, so as to achieve the purpose of denoising [13, 14].

Spline function plays an important role in wavelet theory. Any scaling function (or wavelet) can be written as a convolution of polynomial B-spline and singular distribution [15]. Compared with other wavelet functions, spline wavelets have clearer mathematical expressions and a higher degree of approximation. B-spline is the basis function that constitutes the spline function [17]. In 2000, Unser proposed the fractional B-spline function [18, 19], which extended the original integer-order B-spline to the fractional order, and constructed the fractional-order B-spline wavelet. Fractional B-spline wavelets have good properties such as decay, two-scale equations and fractional approximation. In the literature [20], the meaning of fractional order is explained from the perspective of filter, and it is proved that the fractional B-spline wavelet can well describe the features with fractional singularity such as fine texture in the image [21].

As an edge-preserving smoothing filter, TV is widely used due to its simplicity and effectiveness in removing noise-like structures [5]. However, it penalizes larger gradient amplitudes and may reduce the contrast during the smoothing process. Therefore, a method called truncated TV is proposed, which uses a high-sparse regularization method to establish its mathematical model [6]. Compared with TV, truncated TV attempts to partially penalize the gradient of the image to ensure that only important structures in the image are retained and uninteresting structures in the image are smoothed. In addition, the algorithm can take full advantage of the existing fast TV solver without multiple modifications, so that it has the same computational efficiency as TV [6].

Based on the characteristics of the fractional B-spline wavelet transform and the truncated TV model, this paper proposes a truncated TV image denoising model in the fractional B-spline wavelet domain (referred to as the fractional B-spline truncated TV model). The proposed truncated TV image denoising model based on fractional B-spline wavelet transform possesses three significant characteristics: 1) the fractional B-spline wavelet has the characteristics of characterizing the fractional singularity of the sample; 2) the truncated total variation (truncation TV) model solves the problem of large penalties for the total variation (TV) model; 3)Since the staircase often becomes the high-frequency components in the fractional wavelet domain, the fractional wavelet shrinkage technique can suppress the staircase caused by the truncated TV model to some extent [22, 23].

The structure of the paper is organized as follows. The second section introduces the theory of fractional B-spline wavelet transform and truncated TV model. The third part presents the truncated total variation in fractional B-spline wavelet transform for Micro-CT image denoising, and the fourth part uses this method to conduct experimental analysis on the reconstructed image of the actual biological sample pigeon. Finally, the fifth section gives some concluding remarks.

2 Methodology

2.1 Fractional B-spline wavelet transform

By analogy with the classical B-splines, we define the fractional causal B-splines by taking the (α + 1)th fractional difference of the one-sided power function [18]:

$β_{+}^{α} (x) = \frac{1}{Γ (α + 1)} Δ_{+}^{α + 1} x_{+}^{α} = \frac{1}{Γ (α + 1)} \sum_{k ⩾ 0} {(- 1)}^{k} (\begin{matrix} α + 1 \\ k \end{matrix}) {(x - k)}_{+}^{α}$ (1) where $Γ (α + 1) = \int_{0}^{+ \infty} x^{α} e^{- x} dx$ , ${(x - k)}_{+}^{α} = {(max (x - k, 0))}^{α}$ , $(\begin{matrix} α \\ k \end{matrix}) = \frac{Γ (α + 1)}{Γ (k + 1) Γ (α - k + 1)}$ .

We also introduce the reversed versions of these functions: the anticausal B-splines of degree α [18], $β_{-}^{α} (x) = Δ_{-}^{α + 1} (\frac{x_{-}^{α}}{Γ (α + 1)}) = β_{+}^{α} (- x)$ (2) where $x_{-}^{α} = {(- x)}_{+}^{α}$ and where $Δ_{-}^{α}$ denotes the fractional backward difference operator $Δ_{-}^{α} f (x) = \sum_{k ⩾ 0} {(- 1)}^{k} (\begin{matrix} α \\ k \end{matrix}) f (x + k)$ (3)

To symmetrize the construction and to be in the position to calculate fractional B-spline inner products, the centered fractional B-splines of degree α are defined as [18]: $β_{*}^{α} (x) = β_{+}^{\frac{α - 1}{2}} (x) * β_{-}^{\frac{α - 1}{2}} (x)$ (4) where * denotes convolution operator and $β_{*}^{α}$ can be calculated as:

$β_{*}^{α} (x) = \frac{1}{Γ (α + 1)} Δ_{*}^{α} {| x |}_{*}^{α} = \frac{1}{Γ (α + 1)} \sum_{k = Z} {(- 1)}^{k} | \begin{matrix} α + 1 \\ k \end{matrix} | {| x - k |}_{*}^{α}$ where ${| x |}_{*}^{α} = {\begin{matrix} \frac{{| x |}^{α}}{- 2 sin (π α / 2)}, α \neq 2 m (noteven) \\ \frac{x^{2 m} log | x |}{(- 1)^{1 + m} π}, α = 2 m (even) \end{matrix}$ $| \begin{matrix} α \\ k \end{matrix} | = (\begin{matrix} α \\ k + \frac{α}{2} \end{matrix})$

Fractional B-spline function has decay, Riesz basis, approximation, and satisfies two-scale relation. Therefore, the fractional B-spline function cluster can form a multi-resolution analysis, so that the corresponding wavelet basis function can be constructed and the image can be transformed by fractional B-spline wavelet.

In particular, the fractional B-spline function satisfies the two-scale relation $β_{+}^{α} (\frac{x}{2}) = \sum_{k \in ℤ} h_{+}^{α} (k) β_{+}^{α} (x - k)$

This equation can be established by direct manipulation of the time domain formulae. The low-pass filter corresponding to the fractional B-spline wavelet transform is obtained by performing Fourier transform on the two-scale relationship [19]: $H_{+}^{α} (e^{j ω}) = \sqrt{2} {(\frac{1 + e^{- j ω}}{2})}^{α + 1}$ (5)

Orthogonalize $H_{+}^{α}$ directly to get: $H_{⊥}^{α} (e^{j ω}) = H_{+}^{α} (e^{j ω}) \sqrt{\frac{A^{2 α + 1} (e^{j ω})}{A^{2 α + 1} (e^{2 j ω})}}$ (6) where $A^{α} (e^{j ω}) = \sum_{n \in Z} e^{- j ω} \int β_{+}^{α} (x) β_{+}^{α} (x + n) dx$ .

The corresponding high-pass filter is [19]: $G_{⊥}^{α} (e^{j ω}) = e^{- j ω} H_{⊥}^{α} (- e^{- j ω}) .$ (7)

The analysis and synthesis process of fractional B-spline wavelet transform is shown in Fig. 1.

Fig. 1

Analysis and synthesis filterbank with fractional spline filters.

2.2 Truncated Total Variation (Truncate TV)

We denote the input image by f and the smoothed result by u. Thus, the truncated TV is expressed as [6] $U (x) = \int (T (u_{x}) + T (u_{y}) dxdy$ (8) $T (u) = {\begin{matrix} | u | | u | < ɛ \\ ɛ otherwise \end{matrix}$ (9) where u_x, u_y) is the gradient of u. Truncated TV only penalizes gradients whose magnitudes are smaller than the threshold ɛ, while for those whose magnitudes are greater than ɛ, the proposed penalty doesn’t penalize them.

Truncated TV is non-convex and difficult to optimize due to its truncated” shape. However, Equation (9) can be re-expressed using L₀ regularization. Taking ɛ as a parameter, T (u) defined in Equation (9) is equivalent to $φ (u, l) = min_{l} {ɛ | l |^{0} + | u - l |}$ (10) where ||⁰ is zero power operator, i.e. |l|⁰ = 1 if l ≠ 0, else |l|⁰ = 0.

3 Proposed method

Combining the fractional B-spline wavelet transform and truncated TV, the noise reduction model in this paper is shown in Equation (11): $u = \underset{u, l_{1}, l_{2}}{arg min} {∥ W_{FBS} f - W_{FBS} u ∥}_{2}^{2} + α {∥ φ (W_{FBS} u_{x}, l_{1}) + φ (W_{FBS} u_{y}, l_{2}) ∥}_{1}$ (11)

The first term on the right side of the above formula is the data fidelity term, and the second term is the fractional B-spline-truncated TV term. Where f is the noise-containing image to be processed, u is the image after noise reduction, W_FBS is the fractional B-spline wavelet transform, and α is the parameter of the truncated TV regular term defined in the wavelet domain. ∥φ (W_FBSu_x, l₁) + φ (W_FBSu_y, l₂) ∥ ₁ defined as: $\begin{matrix} {∥ φ (W_{FBS} u_{x}, l_{1}) + φ (W_{FBS} u_{y}, l_{2}) ∥}_{1} = \sum_{i = 0}^{N - 1} \sum_{j = 0}^{N - 1} (min {ɛ {| l_{1_{i, j}} |}^{0} + | W_{FBS} u_{x_{i, j}} - l_{1_{i, j}} |} \\ + min {ɛ {| l_{2_{i, j}} |}^{0} + | W_{FBS} u_{y_{i, j}} - l_{2_{i, j}} |}) \end{matrix}$ (12) where φ (W_FBSu_x, l₁) = min{ ɛ|l₁|⁰ + |W_FBSu_x - l₁| }, φ (W_FBSu_y, l₂) = min{ ɛ|l₂|⁰ + |W_FBSu_y - l₂| }.

Directly minimizing functional Equation (11) is difficult because it involves L₁ and L₂ penalty terms. We adopt an alternating optimization strategy with split Bregman framework [24]. We introduce two dual variables b₁ and b₂, corresponding to W_FBSu_x - l₁ and W_FBSu_y - l₂ respectively, and re-express the objective function as follows: $\min {\sum_{i = 0}^{N - 1} \sum_{j = 0}^{N - 1} {(W_{F B S} f_{i, j} - W_{F B S} u_{i, j})}^{2} + α ε \sum_{i = 0}^{N - 1} \sum_{j = 0}^{N - 1} ({| l_{1_{i, j}} |}^{0} + {| l_{2_{i_{i, j}}} |}^{0}) + α \sum_{i = 0}^{N - 1} \sum_{j = 0}^{N - 1} (| b_{1_{i, j}} | + | b_{2_{i, j}} |)}$ (13)

where $b_{1} = W_{F B S} u_{x} - l_{1}, b_{2} = W_{F B S} u_{y} - l_{2} .$ By using Lagrangian multipliers, Equation (13) can be converted to an unconstraint problem:

$\begin{matrix} m i n {\sum_{i = 0}^{N - 1} \sum_{j = 0}^{N - 1} {(W_{F B S} f_{i, j} - W_{F B S} u_{i, j})}^{2} + α ε \sum_{i = 0}^{N - 1} \sum_{j = 0}^{N - 1} ({| l_{1_{i, j}} |}^{0} + {| l_{2_{i, j}} |}^{0}) + α \sum_{i = 0}^{N - 1} \sum_{j = 0}^{N - 1} (| b_{1_{i, j}} | + | b_{2_{i, j}} |) \\ + λ \sum_{i = 0}^{N - 1} \sum_{j = 0}^{N - 1} {(b_{1_{i, j}} - W_{F B S} u_{x_{i, j}} + l_{1_{i, j}})}^{2} + λ \sum_{i = 0}^{N - 1} \sum_{j = 0}^{N - 1} {(b_{2_{i, j}} - W_{F B S} u_{y_{i, j}} + l_{2_{i, j}})}^{2}} \end{matrix}$ (14)

By introducing Bregman distance, Equation (14) becomes: $\begin{matrix} min {\sum_{i = 0}^{N - 1} \sum_{j = 0}^{N - 1} {(W_{FBS} f_{i, j} - W_{FBS} u_{i, j})}^{2} + α ɛ \sum_{i = 0}^{N - 1} \sum_{j = 0}^{N - 1} ({| l_{1_{i, j}} |}^{0} + {| l_{2_{i, j}} |}^{0}) + α \sum_{i = 0}^{N - 1} \sum_{j = 0}^{N - 1} (| b_{1_{i, j}} | + | b_{2_{i, j}} |) \\ + λ \sum_{i = 0}^{N - 1} {\sum_{j = 0}^{N - 1} (b_{1_{i, j}} - W_{FBS} u_{x_{i, j}} + l_{1_{i, j}} - t_{1_{i, j}})}^{2} + λ \sum_{i = 0}^{N - 1} {\sum_{j = 0}^{N - 1} (b_{2_{i, j}} - W_{FBS} u_{y_{i, j}} + l_{2_{i, j}} - t_{2_{i, j}})}^{2}} \end{matrix}$ (15)

The above joint minimizing problem can be solved alternately by decoupling it into several subproblems, described as follows:

(1) Calculate W_FBSu subproblem with fixed l₁, l₂, b₁, b₁, t₁, t₂ $\begin{matrix} W_{FBS} u^{i + 1} = min_{W_{FBS} u} {\sum_{i = 0}^{N - 1} \sum_{j = 0}^{N - 1} {(W_{FBS} f_{i, j} - W_{FBS} u_{i, j})}^{2} \\ + λ \sum_{i = 0}^{N - 1} \sum_{j = 0}^{N - 1} ({(b_{1_{i, j}}^{i} - W_{FBS} u_{x_{i, j}} + l_{1_{i, j}}^{i} - t_{1_{i, j}}^{i})}^{2} + {(b_{2_{i, j}}^{i} - W_{FBS} u_{y_{i, j}} + l_{2_{i, j}}^{i} - t_{2_{i, j}}^{i})}^{2})} \end{matrix}$ (16)

We derive the Euler-Lagrange equation as follows: $W_{FBS} f - W_{FBS} u + λ (\frac{\partial (b_{1}^{i} + l_{1}^{i} - t_{1}^{i})}{\partial x} + \frac{\partial (b_{2}^{i} + l_{2}^{i} - t_{2}^{i})}{\partial y} - Δ W_{FBS} u) = 0$ (17) where Δ is the Laplace operator. Equation (17) can be efficient solved by using Gauss-Seidel iteration algorithm or FFT operator.∥(2) Calculate b₁ and b₂ with fixed W_FBSu, l₁, l₂, t₁, t₂

The unique minimizer of this subproblem can be obtained by applying the shrinkage operator: $b_{1}^{i + 1} = shrink (W_{FBS} u_{x}^{i + 1} - l_{1}^{i} + t_{1}^{i}, \frac{α}{λ})$ (18) $b_{2}^{i + 1} = shrink (W_{FBS} u_{y}^{i + 1} - l_{2}^{i} + t_{2}^{i}, \frac{α}{λ})$ (19) where $shrink (x, α) = \frac{x}{| x |} max (x - α, 0)$ (20)

(3) Update t₁, t₂ with fixed W_FBSu, l₁, l₂, b₁, b₂ $t_{1}^{i + 1} = W_{FBS} u_{x}^{i + 1} + t_{1}^{i} - l_{1}^{i} - b_{1}^{i + 1}$ (21) $t_{2}^{i + 1} = W_{FBS} u_{y}^{i + 1} + t_{2}^{i} - l_{2}^{i} - b_{2}^{i + 1}$ (22)

(4) Update l₁, l₂ with fixed W_FBSu, t₁, t₂, b₁, b₂ $\begin{matrix} l_{1}^{i + 1}, l_{2}^{i + 1} = \begin{matrix} arg min \\ l_{1}, l_{2} \end{matrix} {α ɛ \sum_{i = 0}^{N - 1} \sum_{j = 0}^{N - 1} ({| l_{1_{i, j}} |}^{0} + {| l_{2_{i, j}} |}^{0}) + λ \sum_{i = 0}^{N - 1} \sum_{j = 0}^{N - 1} ({(l_{1_{i, j}} - W_{FBS} u_{x_{i, j}}^{i + 1} + b_{1_{i, j}}^{i + 1} - t_{1_{i, j}}^{i + 1})}^{2} \\ + {(l_{2_{i, j}} - W_{FBS} u_{y_{i, j}}^{i + 1} + b_{2_{i, j}}^{i + 1} - t_{2_{i, j}}^{i + 1})}^{2})} \end{matrix}$ (23)

The solution of this subproblem is: $l_{1}^{i + 1} = {\begin{matrix} t_{1}^{i + 1} - b_{1}^{i + 1} + W_{FBS} u_{x}^{i + 1} if | t_{1}^{i + 1} - b_{1}^{i + 1} + W_{FBS} u_{x}^{i + 1} | < \sqrt{\frac{α ɛ}{λ}} \\ 0 otherwise \end{matrix}$ (24) $l_{2}^{i + 1} = {\begin{matrix} t_{2}^{i + 1} - b_{2}^{i + 1} + W_{FBS} u_{y}^{i + 1} if | t_{2}^{i + 1} - b_{2}^{i + 1} + W_{FBS} u_{y}^{i + 1} | < \sqrt{\frac{α ɛ}{λ}} \\ 0 otherwise \end{matrix}$ (25)

(5) Finally, perform the inverse transform of the fractional B-spline wavelet transform on W_FBSu to obtain u.∥The whole algorithm is summarized as Table 1. It can be seen from Table 1 that all sub-problems have closed form solutions, so that the algorithm can be solved very quickly.

Table 1

The steps of the fractional B-spline-truncated TV image noise reduction algorithm

Algorithm 1
Input: f, N, b₁ = b₂ = t₁ = t₂ = l₁ = l₂ = 0, α_{W _FBS}, τ_{W _FBS}, ɛ_TTV, α_TTV, λ_TTV;
Perform the fractional B-spline wavelet transform on f to get W_FBSf, W_FBSu = W_FBSf;
For i from 1 to N:
Update W_FBSu by solving (17);
Update b₁, b₂ using (18) and (19);
Update t₁, t₂ using (21) and (22);
Update l₁, l₂ using (24) and (25):
End
Perform the inverse transform of fractional B-spline wavelet transform on W_FBSu to obtain u.
Output: u

4 Experiments

4.1 Experimental data

Use MATLAB R2017a software to realize the real image noise reduction experiment on the computer side (Intel(R) Core(TM)i3-5005U CPU @2.00 GHz,4.00GB). Two groups of real biological sample pigeon CT reconstruction images are used for image denoising experiments.

The pigeon samples were scanned at a voltage of 100 kV and a current of 110μA. The source to rotation axis distance is 305 mm, the source to detector distance is 852.776 mm. As shown in Fig. 2, the left is the fully projected reconstructed image (FDK reconstruction from full projections), and the right is the sparse projection reconstructed image (FDK reconstruction from 200 projections). The size of the two groups of pictures is 1024×1024 and 2048×2048.

Fig. 2

Comparison of full projection reconstruction image and sparse projection reconstruction image. The grayscale images were normalized to the range [0, 255], and the display window was [0, 255].

Fig. 3

The PSNR value of the first denoised image varies with parameters.

Fig. 4

The PSNR value of the second denoised image varies with parameters.

As shown in the rectangular frame in the figure, after comparing the images and the partial enlargement of the image, due to the sparse projection CT reconstruction image contains obvious noise, the internal structure is not clear, the detail texture is blurred, and there are obvious strip artifacts on the periphery.

4.2 Evaluation index

Quantitative evaluation uses three indicators, peak signal-to-noise ratio (PSNR), mean square error (MSE), and structural similarity index (SSIM), to evaluate the quality of images [25, 26].

Define two values, one is MSE, and the other is PSNR, the formula is as follows: $MSE (x, y) = \frac{1}{M \times N} \sum_{i = 1}^{M} {\sum_{j = 1}^{N} (x_{i, j} - y_{i, j})}^{2}$ (26) $PSNR (x, y) = 10 {log}_{10} (\frac{{Peak}^{2}}{MSE}) = 20 {log}_{10} (\frac{Peak}{\sqrt{MSE}})$ (27)

Among them, x represents the reference image, y represents the noise reduction image, M × N is the size of x and y; x_i,j and y_i,j represent the pixel intensity of x and y in a certain pixel (i, j), and Peak represents the maximum pixel intensity in the normalized image (Usually the gray level of the image, generally the value is 255). The general value range of PSNR:20–40, the larger the value, the better the image quality. SSIM is a full-reference image quality evaluation index, which measures the similarity of two images from three aspects: brightness, contrast, and structure. The range of SSIM is 0 to 1. The larger the SSIM value, the higher the similarity. Given two images x and y, the SSIM of the two images can be calculated as follows:

$SSIM (x, y) = \frac{(2 μ_{x} μ_{y} + c_{1}) (2 σ_{xy} + c_{2})}{(μ_{x}^{2} + μ_{y}^{2} + c_{1}) (σ_{x}^{2} + σ_{y}^{2} + c_{2})}$ (28) where μ_x is the average of x, μ_y is the average of y, $σ_{x}^{2}$ is the variance of x, $σ_{x}^{2}$ is the variance of y, and σ_xy is the covariance of x and y. c₁ and c₂ are constants used to maintain stability. The quantitative results are shown in Table 2. The PSNR, MSE, and SSIM of the noisy reconstructed images to be processed for the first biological sample are 22.3229 dB, 380.8818, and 0.6228. It can be clearly seen that the fractional B-spline-truncated TV model noise reduction image has the highest PSNR, SSIM, and lower MSE, and the noise reduction effect is better, and the image structure can be better maintained.

Table 2

PSNR, MSE, and SSIM between full projection reconstruction image and sparse projection reconstruction image

	PSNR (dB)	MSE	SSIM
Sample I	22.32	380.88	0.5151
Sample II	18.61	895.19	0.5026

The quantitative results (PSNR, MSE and SSIM) between full projection reconstructed image and sparse projection reconstructed image. are shown in Table 2.

5 Results and discussion

Apply fractional B-spline-truncated TV model and existing truncated total variation (truncated TV) method [6], wavelet transform soft threshold noise reduction (WT) [9], weighted least squares (WLS) method [27], and L₀ smoothing method [28] to reduce the noise of the CT reconstructed image. In the parameter selection of the fractional B-spline-truncated TV model algorithm, there are mainly five parameters: α_{W
_FBS}, τ_{W
_FBS}, ɛ_TTV, α_TTV and λ_TTV. Where α_{W
_FBS} is the degree of a fractional B spline and τ_{W
_FBS} is dummy shift. ɛ_TTV, α_TTV and λ_TTV are the same parameters as truncated TV. These parameters have a great influence on the results of noise reduction, so it is very important to optimize their values. In this article, we use a greedy strategy to determine the best parameter values one by one. As shown in the figure, PSNR is used as the evaluation index. In order to determine the best value of a certain parameter, while keeping other parameters unchanged, select different parameter values to obtain different noise-reduction images, and then calculate the corresponding PSNR value, and determine the optimal parameter value from the highest PSNR. The parameter values of other algorithms are also determined in this way.

In actual experiments, we combined the three indicators of PSNR, MSE and SSIM to determine the best parameter sets for different algorithms. In the truncated TV algorithm, ɛ is the threshold parameter in the model, α is the parameter of the truncated TV regular term, and λ is the parameter introduced in the iterative solution of each sub-problem in the model. In the wavelet transform soft threshold denoising algorithm, the selected wavelet function coif3 is subjected to two-layer wavelet decomposition, and the best threshold vector is obtained through experiments. λ is the balance parameter between the data item and the smoothness in the WLS algorithm [27]. Increasing λ will produce a smoother image. And α gives a degree of control over the affinities by non-linear scaling the gradients. Increasing α will result in sharper preserved edges. λ is a smoothing parameter that controls the degree of smoothness in the L₀ algorithm, and κ is a parameter that controls the rate [28]. A smaller κ value results in more iterations and sharper edges.

5.1 The first biological sample

Table 3 shows the best values of all the parameters of several noise reduction algorithms corresponding to the image to achieve their best noise reduction results. From Fig. 5, when the image to be processed is degraded due to noise, several noise reduction methods can reduce the image noise, but the image after the noise reduction method in this article looks the most clear, internal details and external edges are more obvious. The frame in the figure represents the partial magnification area.

Table 3
The best parameter values of different methods in the first biological sample

Truncated TV WT WLS L₀ smoothing fractional B-spline-truncated TV model

ɛ = 0 . 001 ; α = 0 . 8 ; λ = 7 Threshold vector [26,14, 26,14] λ = 0 . 1 ; α = 1 λ = 0 . 0006 ; κ = 2 α_{W
_FBS} = 5 ; τ_{W
_FBS} = 0 . 05 ; ɛ_TTV = 0 . 001 ; α_TTV = 0 . 1 ; λ_TTV = 0.03

Truncated TV	WT	WLS	L₀ smoothing	fractional B-spline-truncated TV model
ɛ = 0 . 001 ; α = 0 . 8 ; λ = 7	Threshold vector [26,14, 26,14]	λ = 0 . 1 ; α = 1	λ = 0 . 0006 ; κ = 2	α_{W _FBS} = 5 ; τ_{W _FBS} = 0 . 05 ; ɛ_TTV = 0 . 001 ; α_TTV = 0 . 1 ; λ_TTV = 0.03

Fig. 5

Comparison of noise reduction results for the first biological sample (a) original CT image; (b) truncated TV; (c) WT; (d) WLS; (e) L₀ smoothing; (f) fractional B-spline-truncated TV model. The grayscale images were normalized to the range [0, 255], and the display window was [0, 255].

After partially magnifying the image, compare the processing effects of several methods on image details and edges. As shown by the arrow in Fig. 6, the internal structure of the image is clearer after the truncated TV denoising, and the image can retain the detailed texture after the wavelet threshold denoising. After smoothing and WLS denoising, the image structure and details are smoothed, and part of the image information is lost. The method in this paper can have the advantages of truncated TV and fractional wavelet transform in terms of noise reduction. It can effectively reduce noise while maintaining image structure and detailed texture.

Fig. 6

Comparison of internal details for the first biological sample (a) original CT image; (b) truncated TV; (c) WT; (d) WLS; (e) L₀ smoothing; (f) fractional B-spline-truncated TV model. The grayscale images were normalized to the range [0, 255], and the display window was [0, 255].

As shown by the arrow in Fig. 7, among several methods, the image edge of the method proposed in this article is the clearest, WLS is relatively blurred, and the remaining methods are relatively clear. For the obvious streak artifacts at the periphery of the image, several methods have not been able to effectively remove them. as shown by the elliptical frame in Fig. 7.

Fig. 7

Comparison of external edges for the first biological sample (a) original CT image; (b) truncated TV; (c) WT; (d) WLS; (e) L₀ smoothing; (f) fractional B-spline-truncated TV model. The grayscale images were normalized to the range [0, 255], and the display window was [0, 255].

The quantitative results are shown in Table 4. The PSNR, MSE, and SSIM of the noisy reconstructed images to be processed for the first biological sample are 22.3229 dB, 380.8818, and 0.6228. It can be clearly seen that the fractional B-spline-truncated TV model noise reduction image has the highest PSNR, SSIM, and lower MSE, and the noise reduction effect is better, and the image structure can be better maintained.

Table 5 shows the time consumption of several algorithms. When processing the first biological sample image, several methods take less time, of which WLS is relatively long.

Table 4

PSNR, MSE, and SSIM of the first biological sample processed with several methods of noise reduction

	Truncated TV	WT	WLS	L₀ smoothing	fractional B-spline-truncated TV model
PSNR (dB)	22.5537	22.5087	22.4519	22.4664	26.3019
MSE	361.1696	364.9335	369.7316	368.5037	152.3685
SSIM	0.6881	0.6837	0.6822	0.6810	0.7727

Table 5

Comparison of calculation costs of different methods in the first biological sample

Algorithm	Image size	Time (seconds)
Truncated TV	1024×1024	3.692015
WT	1024×1024	2.721092
WLS	1024×1024	7.006485
L₀ smoothing	1024×1024	5.637512
Fractional B-spline-truncated TV model	1024×1024	4.587047

5.2 The second biological sample

Table 6 shows the best values of all parameters of several noise reduction algorithms. It can be seen from Fig. 8 and the partially enlarged Figs. 9 and 10 that the proposed method enables the internal structure and detailed texture of the image to be well maintained, as well as the edges are clearer. As shown by the arrow and the elliptical frame in the figure.

The PSNR, MSE, and SSIM of the noisy reconstructed images to be processed for the second biological sample are 18.6116 dB, 895.1922 and 0.5939, adjust all the parameters of several noise reduction method algorithms to obtain high PSNR, SSIM and low MSE. It can be seen from Table 7 that PSNR, MSE, and SSIM of the proposed method are better than other methods.

Table 6
The best parameter values of different methods in the second biological sample

Truncated TV WT WLS L₀ smoothing fractional B-spline-truncated TV model

ɛ = 0 . 0045 ; α = 0 . 5 ; λ = 9 Threshold vector [20,20, 20,20] λ = 0 . 2 ; α = 1 . 2 λ = 0 . 001 ; κ = 2 τ_{W
_FBS} = -0 . 25 ; ɛ_TTV = 0.004 ; α_TTV = 0.5 ; λ_TTV = 0.02

Truncated TV	WT	WLS	L₀ smoothing	fractional B-spline-truncated TV model
ɛ = 0 . 0045 ; α = 0 . 5 ; λ = 9	Threshold vector [20,20, 20,20]	λ = 0 . 2 ; α = 1 . 2	λ = 0 . 001 ; κ = 2	τ_{W _FBS} = -0 . 25 ; ɛ_TTV = 0.004 ; α_TTV = 0.5 ; λ_TTV = 0.02

Fig. 8

Comparison of noise reduction results for the second biological sample (a) original CT image; (b) truncated TV; (c) WT; (d) WLS; (e) L₀ smoothing; (f) fractional B-spline-truncated TV model. The grayscale images were normalized to the range [0, 255], and the display window was [0, 255].

Fig. 9

Comparison of internal details for the second biological sample (a) original CT image; (b) truncated TV; (c) WT; (d) WLS; (e) L₀ smoothing; (f) fractional B-spline-truncated TV model. The grayscale images were normalized to the range [0, 255], and the display window was [0, 255].

Fig. 10

Comparison of external edges for the second biological sample (a) original CT image; (b) truncated TV; (c) WT; (d) WLS; (e) L₀ smoothing; (f) fractional B-spline-truncated TV model. The grayscale images were normalized to the range [0, 255], and the display window was [0, 255].

Table 7

PSNR, MSE, and SSIM of the second biological sample processed with several methods of noise reduction

	Truncated TV	WT	WLS	L₀ smoothing	fractional B-spline-truncated TV model
PSNR (dB)	18.7950	18.7303	18.3808	18.7390	24.2394
MSE	858.1855	871.0662	944.0603	869.3158	244.9836
SSIM	0.7041	0.6839	0.6927	0.6873	0.8980

Table 8 shows the time-consuming for several algorithms. When processing the second biological sample image, the WLS algorithm takes a long time compared to other algorithms due to the large image size, while the method in this paper takes a mid-range amount of time.

Table 8

Comparison of calculation costs of different methods in the second biological sample

Algorithm	Image size	Time (seconds)
Truncated TV	2048×2048	13.993614
WT	2048×2048	3.356346
WLS	2048×2048	457.006485
L₀ smoothing	2048×2048	21.684386
fractional B-spline-truncated TV model	2048×2048	17.549406

6 Conclusion

In this article, we propose a truncated TV image denoising model based on the fractional B-splines wavelet domain to suppress noise and maintain the structure of the image. Through the noise reduction processing of the noisy biological sample pigeon CT image, the experimental results show that, compared with other methods, this new model can effectively suppress noise, and the processed image effect is better. In addition, the model can maintain high PSNR, SSIM, and low MSE.

In our study, the truncated total variation in fractional B-spline wavelet transform domain is effective for two groups of real biological sample compared other methods (e.g., truncated TV method, WT, WLS, and L₀ smoothing method). However, the shortcoming of the proposed algorithm is that its parameters should be tuned for the best performance, and there is no proof that the new model is true for all kinds of structures. In recent years, many wavelet constructions (e.g., contourlet, curvelet, shearlet) have been developed and used for different image processing tasks. Take shearlet, for instance. The shearlet is more effective in approximating piecewise smooth images containing rich geometric information. We believe that shearlet combined with the truncated total variation can be applied for samples with complex textures in image denoising.

Footnotes

Acknowledgments

This work was supported by the Foundation of Tianjin University of Technology and Education [Grant Nos. KJ12-01, KJ17-36].

References

Han

, Bian

, Eaker

et al. Image reconstruction from a reduced number of projections in Micro-CT specimen imaging, IEEE Nuclear Science Symposuim & Medical Imaging Conference, 2010, pp. 3466–3469.

Han

, Bian

, Eaker

et al. Algorithm-enabled low-dose micro-CT imaging, IEEE Transactions on Medical Imaging 30 (2011), 606–620.

Diwakar

, Kumar

A review on CT image noise and its denoising, Biomedical Signal Processing and Control 42 (2018), 73–88.

Diwakar

, Kumar

, Singh

A.K.

CT image denoising using NLM and its method noise thresholding, Multimedia Tools and Applications 79 (2020), 14449–14464.

Cai

, Wang

, Zhong

et al. Total variation combining nonlocal means filtration for image reconstruction in X-ray computed tomography, Journal of X-Ray Science and Technology 30 (2022), 613–630.

Dou

, Song

, Gao

, Jiang

Image smoothing via truncated total variation, IEEE Access 5 (2017), 27337–27344.

Zhao

, Hoffman

, McNitt-Gray

et al. Ultra-low-dose CT image denoising using modified BM3D scheme tailored to data statistics, Medical Physics 46 (2019), 190–198.

Wolterink

J.M.

, Leiner

, Viergever

M.A.

et al. Generative adversarial networks for noise reduction in low-dose CT, IEEE Transactions on Medical Imaging 36 (2017), 2536–2545.

Mallat

S.G.

A theory for multiresolution signal decomposition: The wavelet representation, IEEE Transactions on Pattern Analysis and Machine Intelligence 11 (1989), 674–693.

10.

Livstone

Wavelets: A conceptual overview, Area I Examination Report (1994), 1–29.

11.

Chakraborty

, Shaikh

S.H.

, Chakrabarti

, Ghosh

An image denoising technique using quantum wavelet transform, International Journal of Theoretical Physics 59 (2020), 3348–3371.

12.

Luisier

, Blu

, Unser

A new SURE approach to image denoising: Interscale orthonormal wavelet thresholding, IEEE Trans Image Process 16 (2007), 593–606.

13.

Fathi

, Naghsh-Nilchi

A.R.

Efficient image denoising method based on a new adaptive wavelet packet thresholding function, IEEE Trans Image Process 21 (2012), 3981–3990.

14.

Vidakovic

Nonlinear wavelet shrinkage with Bayes rules and Bayes factors, Am Stat Assoc 93 (1998), 173–179.

15.

Unser

, Aldroubi

, Eden

B-spline image processing part I-theory, IEEE Transactions on Signal Processing 40 (1993), 821–832.

16.

Unser

, Aldroubi

, Eden

B-spline image processing part II-efficient design and applications, IEEE Transactions on Signal Processing 40 (1993), 834–848.

17.

Schoenberg

I.J.

Cardinal interpolation and spline functions, Approximation Theory 2 (1969), 167–206.

18.

Unser

, Blu

Fractional splines and wavelets, SIAM Review 42 (2000), 43–67.

19.

Blu

, Unser

The fractional spline wavelet transform: Definition and implementation, IEEE Xplore 1 (2000), 512–515.

20.

Unser

, Blu

Self-Similarity: Part1-splines and operators, IEEE Transactions on Signal Processing 55 (2007), 1352–1363.

21.

Unser

, Blu

Self-Similarity: Part II-optimal estimation of fractal processes, IEEE Transactions on Signal Processing 55 (2007), 1364–1378.

22.

Kamilov

, Bostan

, Unser

Wavelet shrinkage with consistent cycle spinning generalizes total variation denoising, IEEE Signal Process 19 (2012), 187–190.

23.

Shen

, Liu

, Lou

, Hou

Wavelet-based total variation and nonlocal similarity model for image denoising, IEEE Signal Processing Letters 24 (2017), 877–881.

24.

Goldstein

, Osher

The split Bregman method for L1-regularized problems, SIAM Journal on Imaging Sciences 2 (2009), 323–343.

25.

Wang

, Bovik

A.C.

A universal image quality index, IEEE Signal Processing Letters 9 (2002), 81–84.

26.

Wang

, Bovik

A.C.

, Sheikh

H.R.

, Simoncelli

E.P.

Image quality assessment: From error visibility to structural similarity, IEEE Trans Image Process 13 (2004), 600–612.

27.

Farbman

, Fattal

, Lischinski

, Szeliski

Edge-preserving decompositions for multi-scale tone and detail manipulation, ACM Transactions on Graphics 27 (2008), 1–10.

28.

, Lu

, Xu

, Jia

Image smoothing via L0 gradient minimization, ACM Transactions on Graphics 30 (2011), 1–11.