A Generalized Construction Model for CT Projection-Wise Filters on the SDBP Technique

Abstract

Nowadays, with the rapid development of computed tomography (CT), the theoretical expansion of CT reconstruction has become the crucial ingredient to break through technical bottlenecks. This article supplements relevant knowledge of the method for constructing CT projection-wise filters. From the perspective of the mathematical principle of CT reconstruction, the second-order divided-difference back projection (SDBP) technique is firstly proposed, which is an implementation of accurate and efficient inverse Radon transform. On the basis of the SDBP technique, a brand new computational model for filter expressions is derived. The model reveals the correlation between the convolution kernel for data restoration and the filter expression and also specifies the principle of filter construction. On the proposed filter construction model, the decomposition form of filters is discovered for the first time. A series of basic filters are acquired, which can be used to compose filters in actual need. According to the superposition principle, the properties of a filter depend on the decomposed basic filters. The selection of the kernel function on the properties of basic filters clarifies the general rule of filter construction. These theories proposed in this article are of reference value for the filter optimization research in the CT reconstruction field.

1. Introduction

Computed tomography (CT) is a commonly used nondestructive diagnostic imaging technology [1,2]. Reconstruction-algorithm-related technology has always been the research focus. From the perspective of applied mathematics, CT reconstruction is essentially an inverse problem of the Radon transform [3]. Therefore, the properties of the Radon transform have always been the research topic in the CT reconstruction field [4,5,6,7]. Among all algorithms based on the inverse Radon transform, the filtered back projection (FBP) technique is a typical case [8], which is convenient to implement. However, due to the divergence of the ideal ramp function, the FBP technique cannot be equivalent to the rigorous Radon inverse transform. An accurate and efficient inverse Radon transform method is essential to the theoretical research of CT reconstruction.

For discrete projection data, the performance of the CT projection-wise filter is an important factor that affects the CT image quality. Presently, the most commonly used filters in FBP are the R-L filter [9] and S-L filter [10]. The CT image reconstructed by the R-L filter has a clear contour and high spatial resolution, but the image is affected by an obvious Gibbs phenomenon [11]. Compared with the R-L filter, the CT image reconstructed by the S-L filter has a lighter Gibbs phenomenon and less pixel value error, but the sharpness of the contour edge is lower. Additionally, the reconstructed image quality usually needs to be quantified by evaluating indexes [12,13]. On the other hand, not all filters can be obtained by the traditional design method, which is realized by the inverse Fourier transform after windowing the ramp function in the frequency domain [14]. For that reason, scholars have also studied some filter designing methods, which have achieved remarkable results [15,16,17,18]. Recently, a class of neural network and deep learning-based filter optimizing methods have been proposed and shown some advantages of machine learning [19,20,21]. However, most of the proposed methods are aimed at certain targets without widespread applicability. Filter optimizing methods based on machine learning can be used to obtain some local optimal solutions, but in most cases, they cannot be proved to be a global optimal solution. The supplement of a theoretical basis is necessary.

This article has three main modules related to filter construction: the SDBP technique, the computational model for filters, and the basic composition of filters. Among them, the SDBP technique contributes a new approach for inverse Radon transform; the computational model for filters reveals the relationship between filter function and convolution kernel of data preprocessing; the basic composition of filters shows the decomposition and essential elements of filters. Their corresponding theorems are proposed and proved in the order of logical recurrence.

As an extension of the inverse Radon transform in the form of Hilbert transform, this article introduces an inverse Radon transform formula in the form of second-order divided difference back projection (SDBP). Different from the FBP technique, the SDBP technique calculates the second-order divided difference instead of filtering during data preprocessing. Theoretically, the SDBP technique is an approach for accurate CT reconstruction for dense sampling.

On the SDBP technique, a generalized construction model for CT projection-wise filters is proposed. It begins with the discrete projection data being restored into a continuous function by convolution with kernel function. By applying the SDBP formula to the restored second-order difference projection function, the relationship between filter function and convolution kernel is revealed. As an application, a variety of filters can be acquired by applying different kernel functions to the proposed computational model. In particular, the R-L filter and S-L filter can be acquired, respectively, with sinc kernel and rectangular kernel. As an image evaluation index, the root mean square error (RMSE) is used to evaluate the reconstruction accuracy of different filters. The result shows that the delta filter has the minimum reconstruction error among the tested filters.

By further applying the filter construction model, it is proved that any kind of filter can be decomposed into basic filters. Through simulations of CT reconstruction for the Shepp–Logan phantom and the digital resolution chart, some characteristics of the basic filters are demonstrated, which provide guidance for the filter design using other techniques. This study provides sufficient theoretical basis for the analysis of filters and makes a foundation for filter optimizing methods.

2. Inverse Radon Transform with SDBP Technique

The SDBP technique is an approach for the inverse Radon transform. For the projection data of dense sets, the SDBP technique is a practicable way to achieve accurate analytical reconstruction of the target function. The SDBP technique also provides an approach for CT reconstruction with a discrete sinogram, even if the sampling is not equidistant.

2.1. Radon Transform

Each $•$ refers to a variable, and each $\otimes$ refers to a function. Given a target real function $f\left(•\right)$, for $\overrightarrow{x}={\left[\mathit{x},y\right]}^T\in {ℝ}^2$, assume that $f\left(\overrightarrow{\mathit{x}}\right)$ satisfies the following conditions:

Hypothesis 1.

• $f\left(\overrightarrow{\mathit{x}}\right)$is continuous;
• The double integral${\displaystyle \iint \frac{\left|f\left(\overrightarrow{\mathit{x}}\right)\right|}{\sqrt{{\overrightarrow{\mathit{x}}}^T\overrightarrow{\mathit{x}}}}{d}^2\overrightarrow{\mathit{x}}}$extending over the whole plane, converges;
• For any arbitrary point$\overrightarrow{\mathit{x}}$on the plane, $\underset{\epsilon \to 0_{+}}{\lim }{\displaystyle {\int }_0^{2\pi }f\left(\overrightarrow{\mathit{x}}+\frac{\overrightarrow{\mathit{\phi }}}{\epsilon }\right)}d\phi =0$, where$\overrightarrow{\mathit{\phi }}={\left[\cos \phi ,\sin \phi \right]}^T$.

Definition 1.

The Radon transform of$f\left(\overrightarrow{\mathit{x}}\right)$is a linear integral transform in the plane, which is defined as

$$\mathcal{R}\left[f\right]\left(r,\phi \right):={\displaystyle {\int }_{ℝ}f\left(r\overrightarrow{\mathit{\phi }}+\tau {\overrightarrow{\mathit{\phi }}}_{\perp }\right)d\tau },$$

(1)

where$r$is the distance of straight line$\left\{L:{\overrightarrow{\mathit{\phi }}}^T\overrightarrow{\mathit{x}}-r=0\right\}$from the origin;${\overrightarrow{\mathit{\phi }}}_{\perp }={\left[\sin \phi ,-\cos \phi \right]}^T$.

The Radon-transformed functions defined in Equation (1) and the original functions constitute a function space of bijections.

2.2. SDBP Technique

For $r\in ℝ$ and $\phi \in \left(0,2\pi \right]$, given the Radon transform $\mathcal{R}\left[f\right]\left(r,\phi \right)$ of the target function $f\left(\overrightarrow{\mathit{x}}\right)$ to be solved, assume that $\mathcal{R}\left[f\right]\left(r,\phi \right)$ satisfies the following conditions:

Hypothesis 2.

• The integral${\displaystyle {\int }_{ℝ}\left|\mathcal{R}\left[f\right]\left(r,\phi \right)\right|}dr$converges;
• $\underset{s\to 0}{\lim }{\Delta }_r^2\left[\mathit{ℛ}\left[f\right]\left(r,\phi \right)\right]\left(s\right)=0$, where${\Delta }_r^2\left[\otimes \right]\left(s\right)$of$s\in ℝ$is the second-order difference operator with respect to variable$r$, and more specifically

  $${\Delta }_r^2\left[\mathcal{R}\left[f\right]\left(r,\phi \right)\right]\left(s\right):=\mathcal{R}\left[f\right]\left(r-s,\phi \right)-2\mathcal{R}\left[f\right]\left(r,\phi \right)+\mathcal{R}\left[f\right]\left(r+s,\phi \right);$$
• The first-order partial derivative$\frac{\partial }{\partial r}\mathcal{R}\left[f\right]\left(r,\phi \right)$is continuous, or only has removable discontinuities;
• The second-order partial derivative$\frac{{\partial }^2}{\partial r^2}\mathcal{R}\left[f\right]\left(r,\phi \right)$is continuous, or only has removable or jump discontinuities.

On the premise of the Radon transform $\mathcal{R}\left[f\right]\left(r,\phi \right)$ meeting Hypothesis 2, the formula of the SDBP technique for the inverse Radon transform is shown as follows:

Theorem 1.

For$\overrightarrow{\mathit{x}}\in {ℝ}^2$,$r\in ℝ$,$s\in ℝ$, and$\phi \in \left[0,2\pi \right)$,

$$f\left(\overrightarrow{\mathit{x}}\right)=-\frac{1}{4{\pi }^2}\underset{\epsilon \to 0_{+}}{\lim }{\displaystyle {\int }_0^{2\pi }{\displaystyle {\int }_{\epsilon }^{+\infty }\frac{{\Delta }_r^2\left[\mathcal{R}\left[f\right]\left(r,\phi \right)\right]\left(s\right)}{s^2}\left|{}_{r={\overrightarrow{\mathit{\phi }}}^T\overrightarrow{\mathit{x}}}\right.ds}d\phi }.$$

(2)

If $s\to 0$, the integrand is equal to the mean value of second-order partial derivatives of $\mathcal{R}\left[f\right]\left(r,\phi \right)$ on both sides of the index position, as

$$\underset{s\to 0}{\lim }\frac{{\Delta }_r^2\left[\mathcal{R}\left[f\right]\left(r,\phi \right)\right]\left(s\right)}{s^2}=\frac{1}{2}\left(\frac{{\partial }^2}{\partial r^2}\mathcal{R}\left[f\right]\left(r_{+},\phi \right)+\frac{{\partial }^2}{\partial r^2}\mathcal{R}\left[f\right]\left(r_{-},\phi \right)\right).$$

Compared with the inverse Radon transform in the form of the Hilbert transform, the SDBP technique is operable, for it avoids the derivation process and the divergence of the Hilbert-transformed function at singularity. The SDBP technique is not same as the second-order differential operation of lambda tomography [22]. Compared with FBP, SDBP is also quite different in the data processing of projection values. SDBP does not need to introduce a filter function. Its realization process can be summarized as the following three steps:

• Calculation of the second-order divided difference of the projection function under a certain projection angle;
• Linear integration of the second-order divided difference at the projection position of the target point, which is calculated in step 1;
• Rotational integration of the linear integral value calculated in step 2, also known as the back projection progress.

Proof of Theorem 1.

In order to facilitate the subsequent proofs, some properties of Fourier transform are introduced primarily. The rules of functional relations of Fourier transform are shown in Appendix A, in which the rules cited in this article are represented by ‘Rule Number’.

For the target $f\left(\overrightarrow{\mathit{x}}\right)$, on the promise that $f\left(\overrightarrow{\mathit{x}}\right)$ satisfies Hypothesis 1, then according to the properties of Fourier transform, there is

$$f\left(\overrightarrow{\mathit{x}}\right)\hspace{0.17em}=\hspace{0.17em}{\mathcal{F}}_{\overrightarrow{\mathit{\xi }}}^{-1}\left[{\mathcal{F}}_{\overrightarrow{\mathit{x}}}\left[f\right]\right]\left(\overrightarrow{\mathit{x}}\right)\hspace{0.17em}=\hspace{0.17em}{\displaystyle {\iint }_{{ℝ}^2}{\mathcal{F}}_{\overrightarrow{\mathit{x}}}\left[f\right]\left(\overrightarrow{\mathit{\xi }}\right)\cdot \exp \left(2\pi \mathit{i}{\overrightarrow{\mathit{\xi }}}^T\overrightarrow{\mathit{x}}\right)d^2\overrightarrow{\mathit{\xi }}},$$

where $\overrightarrow{\mathit{\xi }}={\left(\xi ,\psi \right)}^T$ are the independent variables in frequency domain corresponding to $\overrightarrow{\mathit{x}}$; $\mathit{i}$ represents the imaginary unit; ${\mathcal{F}}_{\overrightarrow{\mathit{x}}}\left[\otimes \right]\left(\overrightarrow{\mathit{\xi }}\right)$ is the 2D Fourier transform for independent variables $\overrightarrow{\mathit{x}}$; and ${\mathcal{F}}_{\overrightarrow{\mathit{\xi }}}^{-1}\left[\otimes \right]\left(\overrightarrow{\mathit{x}}\right)$ is the inverse 2D Fourier transform for independent variables $\overrightarrow{\mathit{\xi }}$. Replace variable $\overrightarrow{\mathit{\phi }}\rho \leftarrow \overrightarrow{\mathit{\xi }}$, and combined with the Jacobian determinant, there is

$$d^2\overrightarrow{\mathit{\xi }}=\left|\begin{array}{cc}\cos \phi & -\rho \sin \phi \\ \sin \phi & \rho \cos \phi \end{array}\right|d\rho d\phi =\rho d\rho d\phi .$$

The orthogonal direction integral of the rectangular coordinate system is transformed into the rotation direction integral of the polar coordinate system. According to the central slice theorem [23], there is

$$f\left(\overrightarrow{\mathit{x}}\right)\hspace{0.17em}=\hspace{0.17em}\frac{1}{2}{\displaystyle {\int }_0^{2\pi }{\displaystyle {\int }_{ℝ}{\mathcal{F}}_r\left[\mathcal{R}\left[f\right]\left(r,\phi \right)\right]\left(\rho \right)\cdot \exp \left(2\pi \mathit{i}\rho {\overrightarrow{\mathit{\phi }}}^T\overrightarrow{\mathit{x}}\right)\cdot \left|\rho \right|d\rho d\phi }},$$

(3)

where ${\mathcal{F}}_r\left[\otimes \right]\left(\rho \right)$ is the Fourier transform with respect to variable $r$. For $\left|\rho \right|$ in Equation (3), there is

$$\left|\rho \right|\hspace{0.17em}=\hspace{0.17em}\rho \mathrm{sgn}\left(\rho \right)\hspace{0.17em}=\hspace{0.17em}-\mathit{i}\mathrm{sgn}\left(\rho \right)\cdot \mathit{i}\rho ,$$

where $\mathrm{sgn}\left(•\right)$ represents the sign function. According to Rule 6, there is

$$\delta \left(r\right)\hspace{0.17em}=\hspace{0.17em}\mathcal{F}\left[1\right]\left(\rho \right)\hspace{0.17em}=\hspace{0.17em}{\mathcal{F}}_r\left[r\cdot \frac{1}{r}\right]\left(\rho \right)\hspace{0.17em}=\hspace{0.17em}\frac{\mathit{i}}{2\pi }\cdot \frac{d}{d\rho }\left({\mathcal{F}}_r\left[\frac{1}{r}\right]\left(\rho \right)\right),$$

where $\delta \left(•\right)$ represents the Dirac delta function. Combined with the principle of parity invariance,

$${\mathcal{F}}_r\left[\frac{1}{r}\right]\left(\rho \right)\hspace{0.17em}=\hspace{0.17em}-\pi \mathit{i}\mathrm{sgn}\left(\rho \right).$$

By the inverse Fourier transform, there is

$${\mathcal{F}}_{\rho }^{-1}\left[-\mathit{i}\mathrm{sgn}\left(\rho \right)\right]\left(r\right)\hspace{0.17em}=\hspace{0.17em}\frac{1}{\pi r},$$

(4)

where ${\mathcal{F}}_{\rho }^{-1}\left[\otimes \right]\left(r\right)$ is the inverse Fourier transform with respect to variable $\rho$ in the frequency domain. According to Rule 5, there is

$${\mathcal{F}}_{\rho }^{-1}\left[\mathit{i}\rho \cdot {\mathcal{F}}_r\left[\mathcal{R}\left[f\right]\left(r,\phi \right)\right]\left(\rho \right)\right]\left(r\right)\hspace{0.17em}=\hspace{0.17em}\frac{1}{2\pi }\cdot \frac{\partial }{\partial r}\mathcal{R}\left[f\right]\left(r,\phi \right).$$

(5)

By substituting Equations (4) and (5) into Equation (3), and according to Rule 1 and Rule 7, the inverse Radon transform formula in the form of the Hilbert transform is obtained, as

$$f\left(\overrightarrow{\mathit{x}}\right)=\frac{1}{4\pi }{\displaystyle {\int }_0^{2\pi }{\mathbf{H}}_{\tau }\left[\frac{\partial }{\partial \tau }\mathcal{R}\left[f\right]\left(\tau ,\phi \right)\right]\left(r\right)\left|{}_{r={\overrightarrow{\mathit{\phi }}}^T\overrightarrow{\mathit{x}}}\right.d\phi },$$

(6)

where ${\mathbf{H}}_{\tau }\left[\otimes \right]\left(r\right)$ is the Hilbert transform with respect to variable $\tau$.

In Equation (6), the derivative term in the Hilbert transform is sensitive to the minor errors in the projections, which is not conducive to accurate reconstruction. The Hilbert transform is defined using the Cauchy principal value (marked with $\mathbf{p}\mathbf{.}\mathbf{v}\mathbf{.}$) as

$${\mathbf{H}}_{\tau }\left[\frac{\partial }{\partial \tau }\mathcal{R}\left[f\right]\left(\tau ,\phi \right)\right]\left(r\right):=\mathbf{p}\mathbf{.}\mathbf{v}\mathbf{.}{\displaystyle {\int }_{ℝ}\frac{1}{\pi \left(r-\tau \right)}\cdot \frac{\partial }{\partial \tau }\mathcal{R}\left[f\right]\left(\tau ,\phi \right)d\tau }.$$

(7)

On the basis that the anomalous integral in Equation (7) converges, the right side of Equation (7) can be transformed by the partial integration as

$${\displaystyle \int \frac{1}{r-\tau }\cdot \frac{\partial }{\partial \tau }\mathcal{R}\left[f\right]\left(\tau ,\phi \right)d\tau }=\frac{\mathcal{R}\left[f\right]\left(\tau ,\phi \right)+C}{r-\tau }-{\displaystyle \int \frac{\mathcal{R}\left[f\right]\left(\tau ,\phi \right)+C}{{\left(r-\tau \right)}^2}d\tau },$$

(8)

where $C$ is a constant. To ensure the validity of the partial integration, the convergence of the integral on the right side of Equation (8) must be guaranteed, so there is $C=-\mathcal{R}\left[f\right]\left(r,\phi \right)$. Equation (6) is further equivalent to

$$f\left(\overrightarrow{\mathit{x}}\right)=\frac{1}{4{\pi }^2}\cdot \mathbf{p}\mathbf{.}\mathbf{v}\mathbf{.}{\displaystyle {\int }_0^{2\pi }{\displaystyle {\int }_{ℝ}\frac{\mathcal{R}\left[f\right]\left(r,\phi \right)-\mathcal{R}\left[f\right]\left(\tau ,\phi \right)}{{\left(r-\tau \right)}^2}\left|{}_{r={\overrightarrow{\mathit{\phi }}}^T\overrightarrow{\mathit{x}}}\right.d\tau d\phi }}.$$

(9)

Eliminate the singularity of the integrand in Equation (9), then replace variable $s\leftarrow r-\tau$. Equation (2) is obtained, which is the formula of SDBP technique. Theorem 1 is proved. □

3. Computational Model for Filter Expressions

On the basis of the SDBP technique, a filter construction method is proposed. Implemented in the spatial domain, the method reveals the relationship between the filter function and convolution kernel, which can be used to restore discrete sampling data into continuous functions by convolution.

3.1. Projection Data Preprocess

For the target $f\left(\overrightarrow{\mathit{x}}\right)$, on the premise that the projection of $f\left(\overrightarrow{\mathit{x}}\right)$ at any angle is systematic sampling data, it is assumed that the known conditions are as follows:

Hypothesis 3.

• The sampling interval of the detector is$\mathit{d}$,  $\mathit{d}\in {ℝ}_{+}$;
• The projection angles are discrete at equal intervals, and$M$is the number of total projections;
• For$m\in \mathbb{N}$and$n\in \mathbb{Z}$, the projection value is defined as

  $$P_m\left(n\right):={\displaystyle {\int }_{ℝ}\mathcal{R}\left[f\right]\left(r,{\phi }_m\right)\cdot \delta \left(r-n\mathit{d}\right)dr},$$

where${\phi }_m=2\pi m/M$.

In order to acquire the continuous second-order difference projection function at any angle whose domain is $ℝ$, the discrete projection signal should be converted into a continuous function. The continuous second-order difference projection function ${\Delta }_n^2\left[{\tilde{P}}_m\right]\left(\widehat{s}\right)$ of $\widehat{s}\in ℝ$ is constructed by projection values $P_m\left(n\right)$, where $\widehat{s}$ is the normalized variable of $s$ as $\widehat{s}=s/\mathit{d}$. The constructed ${\Delta }_n^2\left[{\tilde{P}}_m\right]\left(\widehat{s}\right)$ is expressed as follows:

Definition 2.

For$m\in \mathbb{N}$,$n\in \mathbb{Z}$,$\nu \in \mathbb{Z}$, and$\widehat{s}\in ℝ$,

$${\Delta }_n^2\left[{\tilde{P}}_m\right]\left(\widehat{s}\right):={\displaystyle \sum _{\nu \in \mathbb{Z}}\left({\Delta }_n^2\left[P_m\right]\left(\nu \right)\cdot k\left(\widehat{s}-\nu \right)\right)},$$

(10)

where$k\left(•\right)$represents the convolution kernel function, and$k\left(\widehat{s}\right)$of$\widehat{s}\in ℝ$satisfies the following conditions:

Hypothesis 4.

For$\nu \in \mathbb{Z}$,$\widehat{s}\in ℝ$,

• $k\left(\widehat{s}\right)$is piecewise continuous, and the integral${\displaystyle {\int }_{ℝ}k\left(\widehat{s}\right)d\widehat{s}}=1$;
• $k\left(\widehat{s}\right)\hspace{0.17em}=\hspace{0.17em}k\left(-\widehat{s}\right)$,$k\left(\nu \right)\left|{}_{\nu \ne 0}\right.=0$, and$\underset{t\to \infty }{\lim }k\left(\widehat{s}\right)=0$;
• $k\left(\widehat{s}\right)$is differentiable at any$\widehat{s}=\nu \ne 0$;
• The second-order derivative$\frac{d^{2}k\left(\widehat{s}\right)}{d{\widehat{s}}^2}$is continuous in the deleted neighborhood of any$\widehat{s}=\nu \ne 0$, and$\frac{d^{2}k\left(\widehat{s}\right)}{d{\widehat{s}}^2}$converges on both sides of any$\widehat{s}=\nu \ne 0$.

3.2. Conversion Relation between Convolution Kernels and Filter Functions

The reconstructed value $f\left(\overrightarrow{\mathit{x}}\right)$ of $\overrightarrow{\mathit{x}}\in {ℝ}^2$ can be approximated by the following formula [14]:

$$\tilde{f}\left(\overrightarrow{\mathit{x}}\right)\hspace{0.17em}=\hspace{0.17em}\frac{\pi \mathit{d}}{M}{\displaystyle \sum _{m=0}^{M-1}{\displaystyle \sum _{n\in \mathbb{Z}}\left(P_m\left(⌊{\overrightarrow{\mathit{\phi }}}_m^T\overrightarrow{\mathit{x}}/\mathit{d}+1/2⌋-n\right)\cdot h\left(n\right)\right)}},$$

(11)

where ${\overrightarrow{\mathit{\phi }}}_m={\left(\cos {\phi }_m,\sin {\phi }_m\right)}^T$; $⌊•⌋$ represents the floor function, which gives the largest integer no greater than the input real number; $h\left(n\right)$ is the expression of filter, which can be obtained by $k\left(\widehat{s}\right)$ as follows:

Theorem 2.

For$n\in \mathbb{Z}\backslash \left\{0\right\}$,$\tau \in ℝ$, there is

$$h\left(n\right)\left|{}_{n\ne 0}\right.=\frac{1}{2\pi {\mathit{d}}^2}{\mathbf{H}}_{\tau }\left[\frac{dk\left(\tau \right)}{d\tau }\right]\left(n\right);$$

(12)

if$n=0$, there is

$$h\left(0\right)\hspace{0.17em}=\hspace{0.17em}-{\displaystyle \sum _{n\in \mathbb{Z}\backslash \left\{0\right\}}h\left(n\right)}.$$

(13)

For any given kernel function, its corresponding filter is uniquely determined. On the other hand, the kernel function corresponding to any given $h\left(n\right)$ is not unique.

Given the analytic function (not unique) of $h\left(n\right)\left|{}_{n\ne 0}\right.$, which is $h\left(\widehat{r}\right)$ of $\widehat{r}\in ℝ$, where $\widehat{r}$ is the normalized variable of $r$ as $\widehat{r}=r/\mathit{d}$, the corresponding convolution kernel of $h\left(\widehat{r}\right)$ is as follows:

Theorem 3.

For$\widehat{s}\in ℝ$,$\tau \in ℝ$, there is

$$k\left(\widehat{s}\right)\hspace{0.17em}=\hspace{0.17em}2\pi {\mathit{d}}^2{\mathbf{H}}_{\tau }\left[{\displaystyle \int h\left(\tau \right)d\tau }\right]\left(\widehat{s}\right).$$

(14)

By applying the proposed computational model, a variety of new filters can be obtained by designing the interpolation functions. Conversely, according to Equation (14), the kernel functions corresponding to some commonly used filters are revealed. Therefore, the filter expression is closely related to the restoration process of the projection signal.

Proof of Theorem 2.

For the discrete projections, according to the SDBP technique, the inverse transform can be approximately expressed as

$$f\left(\overrightarrow{\mathit{x}}\right)\approx -\frac{1}{2\pi M}\underset{\epsilon \to 0+}{\lim }{\displaystyle \sum _{m=0}^{M-1}{\displaystyle {\int }_{\epsilon }^{+\infty }\frac{{\Delta }_r^2\left[\mathcal{R}\left[f\right]\left(r,{\phi }_m\right)\right]\left(s\right)}{s^2}\left|{}_{r={\overrightarrow{\mathit{\phi }}}_m^T\overrightarrow{\mathit{x}}}\right.ds}}.$$

Replace variable $\widehat{s}\leftarrow s/\mathit{d}$, then replace the second-order difference with ${\Delta }_n^2\left[{\tilde{P}}_m\right]\left(\widehat{s}\right)$ in Equation (10). The estimation formula of reconstructed $f\left(\overrightarrow{\mathit{x}}\right)$ can be acquired as

$$\tilde{f}\left(\overrightarrow{\mathit{x}}\right)\hspace{0.17em}=\hspace{0.17em}-\frac{1}{2\pi \mathit{d}M}\underset{\epsilon \to 0+}{\lim }{\displaystyle \sum _{m=0}^{M-1}{\displaystyle {\int }_{\epsilon }^{+\infty }\frac{{\Delta }_n^2\left[{\tilde{P}}_m\right]\left(\widehat{s}\right)}{{\widehat{s}}^2}\left|{}_{n=⌊\frac{{\overrightarrow{\mathit{\phi }}}_m^T\overrightarrow{\mathit{x}}}{\mathit{d}}+\frac{1}{2}⌋}\right.d\widehat{s}}},$$

(15)

According to Equation (10), Equation (15) can be further transformed into

$$\tilde{f}\left(\overrightarrow{\mathit{x}}\right)\hspace{0.17em}=\hspace{0.17em}-\frac{1}{2\pi \mathit{d}M}\underset{\epsilon \to 0+}{\lim }{\displaystyle \sum _{m=0}^{M-1}{\displaystyle \sum _{\nu \in {\mathbb{N}}_{+}}\left({\Delta }_n^2\left[P_m\right]\left(\nu \right)\cdot \mathbf{p}\mathbf{.}\mathbf{v}\mathbf{.}{\displaystyle {\int }_{ℝ}\frac{k\left(\nu -\tau \right)}{{\tau }^2}d\tau }\right)\left|{}_{n=⌊\frac{{\overrightarrow{\mathit{\phi }}}_m^T\overrightarrow{\mathit{x}}}{\mathit{d}}+\frac{1}{2}⌋}\right.}}.$$

(16)

In light of the operational characteristics of discrete convolution, Equation (16) can be transformed into Equation (11) in which

$$h\left(n\right)\left|{}_{n\ne 0}\right.=-\frac{1}{2{\pi }^2{\mathit{d}}^2}\cdot \mathbf{p}\mathbf{.}\mathbf{v}\mathbf{.}{\displaystyle {\int }_{ℝ}\frac{k\left(n-\tau \right)}{{\tau }^2}d\tau }.$$

(17)

By performing partial integral in Equation (17), there is

$$h\left(n\right)\left|{}_{n\ne 0}\right.=\frac{1}{2\pi {\mathit{d}}^2}\cdot \mathbf{p}\mathbf{.}\mathbf{v}\mathbf{.}{\displaystyle {\int }_{ℝ}\frac{1}{\pi \tau }\cdot \frac{dk\left(n-\tau \right)}{d\tau }d\tau }.$$

Additionally, if $k\left(\widehat{s}\right)$ satisfies Hypothesis 4 and is not differentiable at countable discrete points, $\frac{dk\left(\widehat{s}\right)}{d\widehat{s}}$ refers to the generalized derivative. Combined with the definition of the Hilbert transform, Equation (12) can be acquired. Since the sum of all the $P_m\left(n\right)$ coefficients in Equation (16) is zero, Equation (13) holds. Theorem 2 is proved. □

Proof of Theorem 3.

Inversely, if the analytic function $h\left(\widehat{s}\right)$ of $\widehat{s}\in ℝ$ is given, then $h\left(\widehat{s}\right)$ satisfies Equation (17) as

$$h\left(\widehat{s}\right)\hspace{0.17em}=\hspace{0.17em}-\frac{1}{2{\pi }^2{\mathit{d}}^2}\cdot \mathbf{p}\mathbf{.}\mathbf{v}\mathbf{.}{\displaystyle {\int }_{ℝ}\frac{k\left(\widehat{s}-\tau \right)}{{\tau }^2}d\tau }.$$

By the Fourier transform and according to Rule 7, there is

$${\mathcal{F}}_{\tau }\left[h\right]\left(\widehat{\rho }\right)\hspace{0.17em}=\hspace{0.17em}-\frac{1}{2{\pi }^2{\mathit{d}}^2}\cdot {\mathcal{F}}_{\tau }\left[k\right]\left(\widehat{\rho }\right)\cdot {\mathcal{F}}_{\tau }\left[\frac{1}{{\tau }^2}\right]\left(\widehat{\rho }\right).$$

(18)

For ${\mathcal{F}}_{\tau }\left[\frac{1}{{\tau }^2}\right]\left(\widehat{\rho }\right)$, according to Rule 5, there is

$${\mathcal{F}}_{\tau }\left[\frac{1}{{\tau }^2}\right]\left(\widehat{\rho }\right)\hspace{0.17em}=\hspace{0.17em}{\mathcal{F}}_{\tau }\left[-\frac{d}{d\tau }\left(\frac{1}{\tau }\right)\right]\left(\widehat{\rho }\right)\hspace{0.17em}=\hspace{0.17em}-2\pi \mathit{i}\widehat{\rho }\cdot {\mathcal{F}}_{\tau }\left[\frac{1}{\tau }\right]\left(\widehat{\rho }\right)\hspace{0.17em}=\hspace{0.17em}-2{\pi }^2\left|\widehat{\rho }\right|.$$

(19)

By substituting Equation (19) into Equation (18), the following equation is acquired after performing the inverse Fourier transform as

$$k\left(\widehat{s}\right)\hspace{0.17em}=\hspace{0.17em}{\mathit{d}}^2{\mathcal{F}}_{\widehat{\rho }}^{-1}\left(\frac{{\mathcal{F}}_{\tau }\left[h\right]\left(\widehat{\rho }\right)}{\left|\widehat{\rho }\right|}\right)\left(\widehat{s}\right).$$

(20)

In Equation (20), there is

$$\frac{{\mathcal{F}}_{\tau }\left[h\right]\left(\widehat{\rho }\right)}{\left|\widehat{\rho }\right|}=-2\pi \mathit{i}\mathrm{sgn}\left(\widehat{\rho }\right)\cdot \frac{{\mathcal{F}}_{\tau }\left[h\right]\left(\widehat{\rho }\right)}{2\pi \mathit{i}\widehat{\rho }}.$$

(21)

According to Rule 5,

$${\mathcal{F}}_{\tau }\left[h\right]\left(\widehat{\rho }\right)\hspace{0.17em}=\hspace{0.17em}2\pi \mathit{i}\widehat{\rho }\cdot {\mathcal{F}}_{\tau }\left[{\displaystyle \int h\left(\tau \right)d\tau }\right]\left(\widehat{\rho }\right).$$

(22)

By substituting Equation (22) into Equation (21), it can be deducted that

$$\frac{{\mathcal{F}}_{\tau }\left[h\right]\left(\widehat{\rho }\right)}{\left|\widehat{\rho }\right|}=-2\pi \mathit{i}\mathrm{sgn}\left(\widehat{\rho }\right)\cdot {\mathcal{F}}_{\tau }\left[{\displaystyle \int h\left(\tau \right)d\tau }\right]\left(\widehat{\rho }\right).$$

(23)

According to Rule 7 and combined with Equation (4), the following equation is acquired after performing the Fourier transform on both sides of Equation (23) as

$${\mathcal{F}}_{\widehat{\rho }}^{-1}\left(\frac{{\mathcal{F}}_{\tau }\left[h\right]\left(\widehat{\rho }\right)}{\left|\widehat{\rho }\right|}\right)\left(\widehat{s}\right)\hspace{0.17em}=\hspace{0.17em}2\pi \cdot \mathbf{p}\mathbf{.}\mathbf{v}\mathbf{.}{\displaystyle {\int }_{ℝ}\frac{{\displaystyle \int h\left(\tau \right)d\tau }}{\pi \left(\widehat{s}-\tau \right)}d\tau }.$$

(24)

By substituting Equation (24) into Equation (20), and combined with the definition of the Hilbert transform, Equation (14) can be acquired. Theorem 3 is proved. □

It can be clearly observed that there is a reciprocal relationship between filter function $h\left(•\right)$ and kernel function $k\left(•\right)$, which can be realized by the Hilbert transform, shown as follows:

$$h\left(•\right)\underset{\frac{1}{2\pi d^2}{\mathbf{H}}_{\tau }\left[\frac{dk\left(\tau \right)}{d\tau }\right]\left(•\right)}{\overset{2\pi {\mathit{d}}^2{\mathbf{H}}_{\tau }\left[{\displaystyle \int h\left(\tau \right)d\tau }\right]\left(•\right)}{\rightleftarrows }}k\left(•\right).$$

Based on the principles above, the mutual conversion between the filter function and kernel function can be realized. It shows the correlation between the filter expression and convolution kernel, which is not available in the frequency domain windowing method and the Riesz transform method [24].

4. Illustration with Examples

To verify the validity of the proposed method for constructing filters, we apply some common kernels to the proposed model as examples, by which some known filters are derived. We also compare the reconstruction errors of different filters.

4.1. R-L Filter: $\sin \mathrm{c}\left(•\right)$ Kernel

If ${\Delta }_n^2\left[{\tilde{P}}_m\right]\left(\widehat{s}\right)$ of $\widehat{s}\in ℝ$ is the Whittaker–Shannon interpolation of ${\Delta }_n^2\left[P_m\right]\left(\nu \right)$, then the convolution kernel is the normalized sinc function as

$$k\left(\widehat{s}\right)\hspace{0.17em}=\hspace{0.17em}\sin \mathrm{c}\left(\pi \widehat{s}\right).$$

The derivative of the normalized sinc function is

$$\frac{d}{d\tau }\sin \mathrm{c}\left(\pi \tau \right)\hspace{0.17em}=\hspace{0.17em}\frac{\tau \cos \left(\pi \tau \right)-\sin \left(\pi \tau \right)}{\pi {\tau }^2}.$$

According to Equation (12), there is

$$h_{\sin \mathrm{c}}\left(n\right)\left|{}_{n\ne 0}\right.\hspace{0.17em}=\hspace{0.17em}\frac{1}{2\pi {\mathit{d}}^2}\cdot \mathbf{p}\mathbf{.}\mathbf{v}\mathbf{.}{\displaystyle {\int }_{ℝ}\frac{\tau \cos \left(\pi \tau \right)-\sin \left(\pi \tau \right)}{{\pi }^2{\tau }^2\left(n-\tau \right)}d\tau }=-\frac{1-{\left(-1\right)}^n}{2{\pi }^2{n}^2{\mathit{d}}^2},$$

By Equation (13),

$$h_{\sin \mathrm{c}}\left(0\right)\hspace{0.17em}=\hspace{0.17em}-\frac{1}{{\pi }^2{\mathit{d}}^2}{\displaystyle \sum _{n\in {\mathbb{N}}_{+}}\frac{{\left(-1\right)}^n-1}{n^2}}=\frac{3\varsigma \left(2\right)}{2{\pi }^2{\mathit{d}}^2}=\frac{1}{4{\mathit{d}}^2},$$

where $\varsigma \left(•\right)$ represents the Riemann zeta function, and $\varsigma \left(2\right)\hspace{0.17em}=\hspace{0.17em}\frac{{\pi }^2}{6}$. Therefore, the filter corresponding to the normalized sinc kernel is

$$h_{\sin \mathrm{c}}\left(n\right)\hspace{0.17em}=\hspace{0.17em}\left\{\begin{array}{l}\frac{1}{4{\mathit{d}}^2},\hspace{1em}\hspace{1em}\hspace{1em}\hspace{0.17em}\hspace{0.17em}n=0\\ -\frac{1-{\left(-1\right)}^n}{2{\pi }^2{n}^2{\mathit{d}}^2},\hspace{1em}\hspace{0.17em}\hspace{0.17em}n\ne 0\end{array}\right..$$

(25)

Equation (25) is the expression of the R-L filter [9].

4.2. S-L Filter: $\Pi \left(•\right)$ Kernel

If ${\Delta }_n^2\left[{\tilde{P}}_m\right]\left(\widehat{s}\right)$ of $\widehat{s}\in ℝ$ is the proximal interpolation of ${\Delta }_n^2\left[P_m\right]\left(\nu \right)$, then the kernel function is the unit rectangular function, as

$$k\left(\widehat{s}\right)\hspace{0.17em}=\hspace{0.17em}\Pi \left(\widehat{s}\right),$$

where $\Pi \left(•\right)$ represents the unit rectangular function. The generalized derivative of unit rectangular function is

$$\frac{d}{d\tau }\Pi \left(\tau \right):=\delta \left(\tau +1/2\right)-\delta \left(\tau -1/2\right).$$

According to Equation (12), there is

$$h_{\Pi }\left(n\right)\left|{}_{n\ne 0}\right.\hspace{0.17em}=\hspace{0.17em}\frac{1}{2\pi {\mathit{d}}^2}{\displaystyle {\int }_{ℝ}\frac{\delta \left(\tau +1/2\right)-\delta \left(\tau -1/2\right)}{\pi \left(n-\tau \right)}d\tau }=-\frac{2}{{\pi }^2{\mathit{d}}^2\left(4n^2-1\right)}.$$

By Equation (13),

$$h_{\Pi }\left(0\right)=\frac{4}{{\pi }^2{\mathit{d}}^2}{\displaystyle \sum _{n\in {\mathbb{N}}_{+}}\frac{1}{4n^2-1}}=\frac{2}{{\pi }^2{\mathit{d}}^2}.$$

Therefore, the filter corresponding to the unit rectangular kernel is

$$h_{\Pi }\left(n\right)=\left\{\begin{array}{l}\frac{2}{{\pi }^2{\mathit{d}}^2},\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{0.17em}n=0\\ -\frac{2}{{\pi }^2{\mathit{d}}^2\left(4n^2-1\right)},\hspace{1em}\hspace{1em}\hspace{0.17em}n\ne 0\end{array}\right..$$

(26)

Equation (26) is the expression of the S-L filter [10].

4.3. Delta Filter: $\delta \left(•\right)$ Kernel

If ${\Delta }_n^2\left[{\tilde{P}}_m\right]\left(\widehat{s}\right)$ of $\widehat{s}\in ℝ$ is the sampling function with each amplitude ${\Delta }_n^2\left[P_m\right]\left(\nu \right)$, then the convolution kernel is the unit delta function, as

$$k\left(\widehat{s}\right)\hspace{0.17em}=\hspace{0.17em}\delta \left(\widehat{s}\right).$$

According to Equation (12), there is

$$h_{\delta }\left(n\right)\left|{}_{n\ne 0}\right.\hspace{0.17em}=\hspace{0.17em}\frac{1}{2\pi {\mathit{d}}^2}{\displaystyle {\int }_{ℝ}\frac{d\delta \left(\tau \right)}{d\tau }\cdot \frac{1}{\pi \left(n-\tau \right)}d\tau }=-\frac{1}{2\pi {\mathit{d}}^2}{\displaystyle {\int }_{ℝ}\frac{\delta \left(\tau \right)}{\pi {\left(n-\tau \right)}^2}d\tau }=-\frac{1}{2{\pi }^2{n}^2{\mathit{d}}^2}.$$

By Equation (13),

$$h_{\delta }\left(0\right)\hspace{0.17em}=\hspace{0.17em}\frac{1}{{\pi }^2{\mathit{d}}^2}{\displaystyle \sum _{n\in {\mathbb{N}}_{+}}\frac{1}{n^2}}=\frac{\varsigma \left(2\right)}{{\pi }^2{\mathit{d}}^2}=\frac{1}{6{\mathit{d}}^2}.$$

Therefore, the filter corresponding to the unit delta kernel is

$$h_{\delta }\left(n\right)\hspace{0.17em}=\hspace{0.17em}\left\{\begin{array}{l}\frac{1}{6{\mathit{d}}^2},\hspace{1em}\hspace{1em}\hspace{1em}\hspace{0.17em}\hspace{0.17em}n=0\\ -\frac{1}{2{\pi }^2{n}^2{\mathit{d}}^2},\hspace{1em}\hspace{0.17em}\hspace{0.17em}n\ne 0\end{array}\right..$$

(27)

The expression of the delta filter shown in Equation (27) was firstly proposed by Yu-chuan Wei et al. [25], and was introduced as a practical filter to approximate the ideal filter.

4.4. Imaging Tests

The Shepp–Logan phantom of 1024 × 1024 pixels was used to test the CT imaging performance of different example filters. The number of projections was set to 720, and the sampling interval of the detector was equal to the pixel size of the phantom. For further comparison, Gaussian white noise of standard deviation $\sigma$ was added to the sinogram to test the performance of the filters in noise. The reconstructed images are shown in Figure 1.

To reflect the reconstruction accuracy more intuitively, we used the root mean square error (RMSE) to evaluate the error between the reconstructed image and the original image. The RMSE of CT image is defined as follows:

$$\mathrm{RMSE}=\sqrt{\frac{{\displaystyle \sum _{u=1}^N{\displaystyle \sum _{v=1}^N{\left({\tilde{f}}_{u,v}-f_{u,v}\right)}^2}}}{{\displaystyle \sum _{u=1}^N{\displaystyle \sum _{v=1}^N{f}_{u,v}^2}}}},$$

where $f_{u,v}$ is the gray value of the original image pixel in row $u$ and column $v$; ${\tilde{f}}_{u,v}$ is the gray value of the reconstructed image pixel in row $u$ and column $v$; $N$ is the maximum value of $u$ and column $v$. By calculation, the RMSE table of reconstructed images is shown in

Table 1. RMSE of reconstructed images.

Noise Density	R-L Filter	S-L Filter	Delta Filter
$\sigma =0$	0.2672	0.2508	0.2431
$\sigma =1$	0.3231	0.2919	0.2784
$\sigma =5$	0.6544	0.5886	0.5332

, which shows that the reconstructed image of delta filter has the smallest RMSE under each noise intensity.

5. Analysis of Basic Composition of Filters

In order to make further explorations in the rules of the filter properties changing with the kernels, the basic components of filters are analyzed such that an arbitrary filter can be decomposed into a combination of basic filters, which is similar to the Fourier transform. The properties of basic filters provide guidance for the design of the kernel functions.

5.1. Decomposition of Kernel Functions

A generalized function can be described by Schwartz distribution. As for the kernel function $k\left(•\right)$, its Schwartz distribution is defined by

$$〈k,\mathfrak{T}〉:={\displaystyle {\int }_{ℝ}k\left(\tau \right)\cdot \mathfrak{T}\left(\tau \right)d\tau },$$

(28)

where $\mathfrak{T}\left(•\right)$ represents the test function.

By performing discretization, the local function in a definite interval is regarded as a delta function whose amplitude is equal to the interval integral. In this way, the kernel function can be decomposed to a combination of delta functions. Assume that the interval between any adjacent integers is equally divided by piecewise intervals, then the discretized kernel function can be defined as follows:

Definition 3.

For$\nu \in \mathbb{Z}$,$\Lambda \in {ℝ}_{+}$, and$\tau \in ℝℝ$,

$${\tilde{k}}_{\Lambda }\left(\tau \right):={\displaystyle \sum _{\nu \in \mathbb{Z}}\left(\delta \left(\tau -\frac{2\nu +1}{2\Lambda }\right)\cdot {\displaystyle {\int }_{ℝ}k\left(\tau \right)\cdot \Pi \left(\Lambda \tau -\frac{2\nu +1}{2}\right)d\tau }\right)}.$$

(29)

Considering when $\Lambda \to \infty$, by applying Equation (28) to Equation (29), there is

$$〈{\tilde{k}}_{\Lambda },\mathfrak{T}〉\left|{}_{\Lambda \to \infty }\right.=\underset{\Lambda \to \infty }{\lim }{\displaystyle \sum _{\nu \in \mathbb{Z}}\left(\mathfrak{T}\left(\frac{2\nu +1}{2\Lambda }\right)\cdot {\displaystyle {\int }_{ℝ}k\left(\tau \right)\cdot \Pi \left(\Lambda \tau -\frac{2\nu +1}{2}\right)d\tau }\right)}.$$

If $\mathfrak{T}\left(•\right)$ is continuous within its domain, according to the first mean value theorem for definite integrals, the discretized kernel ${\tilde{k}}_{\Lambda }\left(•\right)$ when $\Lambda \to \infty$ and the original kernel $k\left(•\right)$ have the identical Schwartz distribution,

$$〈{\tilde{k}}_{\Lambda },\mathfrak{T}〉\left|{}_{\Lambda \to \infty }\right.=\underset{\Lambda \to \infty }{\lim }{\displaystyle \sum _{\nu \in \mathbb{Z}}\left({\displaystyle {\int }_{ℝ}k\left(\tau \right)\cdot \mathfrak{T}\left(\tau \right)\cdot \Pi \left(\Lambda \tau -\frac{2\nu +1}{2}\right)d\tau }\right)}=〈k,\mathfrak{T}〉.$$

(30)

To describe the basic elements of kernel functions, the delta_2 function ${\delta }_2\left(•,\lambda \right)$ is defined, which is described as follows:

Definition 4.

For$\lambda \in ℝ$,$\tau \in ℝ$,

$${\delta }_2\left(\tau ,\lambda \right):=\frac{1}{2}\left(\delta \left(\tau -\lambda \right)+\delta \left(\tau +\lambda \right)\right),$$

(31)

where$\lambda$represents the displacement of the delta functions from the origin.

According to the conditions in Hypothesis 4, ${\tilde{k}}_{\Lambda }\left(\tau \right)$ can be expressed in the following form as

$${\tilde{k}}_{\Lambda }\left(\tau \right)\hspace{0.17em}=\hspace{0.17em}{\displaystyle \sum _{\nu \in \mathbb{Z}}\left({\mathfrak{a}}_{\Lambda }\left(\nu \right)\cdot {\delta }_2\left(\tau ,\frac{2\nu +1}{2\Lambda }\right)\right)},$$

(32)

where ${\mathfrak{a}}_{\Lambda }\left(\nu \right)$ is the amplitude of delta_2 function,

$${\mathfrak{a}}_{\Lambda }\left(\nu \right)\hspace{0.17em}=\hspace{0.17em}{\displaystyle {\int }_{ℝ}k\left(\tau \right)\cdot \Pi \left(\Lambda \tau -\frac{2\nu +1}{2}\right)d\tau }.$$

According to Hypothesis 4, $\mathfrak{a}\left(\nu \right)$ satisfies the following:

• ${\displaystyle \sum _{\nu \in \mathbb{Z}}{\mathfrak{a}}_{\Lambda }\left(\nu \right)}\hspace{0.17em}=\hspace{0.17em}1$;
• ${\mathfrak{a}}_{\Lambda }\left(\nu \right)\hspace{0.17em}=\hspace{0.17em}{\mathfrak{a}}_{\Lambda }\left(-\nu \right)$;
• For a finite fixed value of $\Lambda$, $\underset{\nu \to \infty }{\lim }{\mathfrak{a}}_{\Lambda }\left(\nu \right)\hspace{0.17em}=\hspace{0.17em}0$.

5.2. Filters of Basic Type

In the definition of Schwartz distribution, the delta_2 function is the basis function of any kernel. According to Equations (12) and (31), if ${\delta }_2\left(\tau ,\lambda \right)$ is the kernel function of $\tau \in ℝ$ and $\lambda \in ℝ\backslash \{\mathbb{Z}\backslash \{0\left\}\right\}$, there is

$$h_{{\delta }_2}\left(n,\lambda \right)\left|{}_{n\ne 0}\right.=-\frac{1}{4\pi {\mathit{d}}^2}{\displaystyle {\int }_{ℝ}\frac{\delta \left(\tau -\lambda \right)+\delta \left(\tau +\lambda \right)}{\pi {\left(n-\tau \right)}^2}d\tau }=-\frac{1}{4{\pi }^2{\mathit{d}}^2}\left(\frac{1}{{\left(n-\lambda \right)}^2}+\frac{1}{{\left(n+\lambda \right)}^2}\right).$$

(33)

By Equation (13), the following equation can be acquired, as

$$h_{{\delta }_2}\left(0,\lambda \right)\hspace{0.17em}=\hspace{0.17em}-\frac{1}{2{\pi }^2{\lambda }^2{\mathit{d}}^2}-\frac{1}{2{\pi }^2{\mathit{d}}^2}\cdot \left(\frac{d}{d\tau }{\displaystyle \sum _{n\in \mathbb{Z}}\left(\frac{1}{n+\tau }\right)}\right)\left|{}_{\tau =\lambda }\right..$$

(34)

The partial fraction decomposition of cotangent function can be used to calculate the exact value of $h\left(0\right)$. For $\tau \in ℝ\backslash \mathbb{Z}$, the sum of the sequence in Equation (34) is

$${\displaystyle \sum _{n\in \mathbb{Z}}\frac{1}{n+\tau }}=\pi \cot \left(\pi \tau \right).$$

Therefore,

$$h_{{\delta }_2}\left(0,\lambda \right)=\frac{1}{2{\mathit{d}}^2}\left(\frac{1}{{\sin }^2\left(\pi \lambda \right)}-\frac{1}{{\pi }^2{\lambda }^2}\right).$$

(35)

Combined with Equation (33) and Equation (35), the filter expression corresponding to delta_2 function kernel is

$$h_{{\delta }_2}\left(n,\lambda \right)=\left\{\begin{array}{l}\frac{1}{2{\mathit{d}}^2}\left(\frac{1}{{\sin }^2\left(\pi \lambda \right)}-\frac{1}{{\pi }^2{\lambda }^2}\right),\hspace{1em}\hspace{1em}\hspace{1em}\hspace{0.17em}\hspace{0.17em}n=0\\ -\frac{1}{4{\pi }^2{\mathit{d}}^2}\left(\frac{1}{{\left(n-\lambda \right)}^2}+\frac{1}{{\left(n+\lambda \right)}^2}\right),\hspace{1em}\hspace{0.17em}\hspace{0.17em}n\ne 0\end{array}\right..$$

(36)

Equation (36) is the expression of basic filter series. In particular, $h_{{\delta }_2}\left(n,0\right)$ is the delta filter as $\lambda =0$.

A filter can be decomposed into a linear combination of basic filters, which is expressed as follows:

Theorem 4.

For$n\in \mathbb{Z}$,$\lambda \in ℝ\backslash \{\mathbb{Z}\backslash \{0\left\}\right\}$,

$$h\left(n\right)={\displaystyle {\int }_{ℝ}k\left(\lambda \right)\cdot h_{{\delta }_2}\left(n,\lambda \right)d\lambda }.$$

(37)

Proof of Theorem 4.

The operation in Equation (12) is a special case of Schwartz distribution expressed in Equation (28). Combined with Equation (30), it can be concluded that ${\tilde{k}}_{\Lambda }\left(•\right)$ and $k\left(•\right)$ correspond to the same filter when $\Lambda \to \infty$, which means

$$\underset{\Lambda \to \infty }{\lim }{\tilde{h}}_{\Lambda }\left(n\right)\hspace{0.17em}=\hspace{0.17em}h\left(n\right),$$

where ${\tilde{h}}_{\Lambda }\left(n\right)$ corresponds to ${\tilde{k}}_{\Lambda }\left(\tau \right)$. On the basis of the linearity in Equation (12), it can be derived from Equation (32) that

$${\tilde{h}}_{\Lambda }\left(n\right)\hspace{0.17em}=\hspace{0.17em}{\displaystyle \sum _{\nu \in \mathbb{Z}}\left({\mathfrak{a}}_{\Lambda }\left(\nu \right)\cdot h_{{\delta }_2}\left(n,\frac{2\nu +1}{2\Lambda }\right)\right)}.$$

(38)

If $k\left(\tau \right)$ is continuous for $\tau \in \left(\frac{\nu }{\Lambda },\frac{\nu +1}{\Lambda }\right)$, then

$$\underset{\Lambda \to \infty }{\lim }{\mathfrak{a}}_{\Lambda }\left(\nu \right)\hspace{0.17em}=\hspace{0.17em}\underset{\Lambda \to \infty }{\lim }\left(\frac{1}{\Lambda }\cdot k\left(\frac{2\nu +1}{2\Lambda }\right)\right).$$

(39)

Suppose that $k\left(\tau \right)$ is continuous for $\tau \in ℝ$. According to Equations (38) and (39), there is

$$\underset{\Lambda \to \infty }{\lim }{\tilde{h}}_{\Lambda }\left(n\right)\hspace{0.17em}=\hspace{0.17em}\underset{\Lambda \to \infty }{\lim }{\displaystyle \sum _{\nu \in \mathbb{Z}}\left(\frac{1}{\Lambda }\cdot k\left(\frac{2\nu +1}{2\Lambda }\right)\cdot h_{{\delta }_2}\left(n,\frac{2\nu +1}{2\Lambda }\right)\right)}.$$

(40)

By means of the definition of integral, Equation (37) can be obtained by the conversion of Equation (40). Under the constraint of Hypothesis 4, it can be proved that Equation (37) is still valid in the case that the kernel function contains finite discontinuities, which include delta function components. Theorem 4 is proved. □

5.3. Imaging Tests

The Shepp–Logan phantom of 1024×1024 pixels was used to test the CT imaging performance of different basic filters. The number of projections was set to 720, and the sampling interval of the detector was equal to the pixel size of the phantom. The images of the Shepp–Logan phantom were reconstructed by a number of basic filters with different $\lambda$, and some results are shown in Figure 2.

In addition to RMSE, we used the average gradient module (AGM) to measure the sharpness of the CT images. The image sharpness has a positive correlation with AGM. The computing formula of the CT image AGM is as follows:

$$\mathrm{AGM}=\frac{{\displaystyle \sum _{u=1}^N{\displaystyle \sum _{v=1}^{N-1}\left|{\tilde{f}}_{u,v}-{\tilde{f}}_{u,v+1}\right|}}+{\displaystyle \sum _{u=1}^{N-1}{\displaystyle \sum _{v=1}^N\left|{\tilde{f}}_{u,v}-{\tilde{f}}_{u+1,v}\right|}}}{2N\left(N-1\right)}.$$

The calculated RMSE and AGM of CT images reconstructed by basic filters are shown in Figure 3 and Figure 4, and the corresponding data set is provided in Appendix B.

If the parameter $\lambda$ is divided into two parts, integer and decimal, as $\lambda ={\lambda }_0+{\lambda }_1$ (${\lambda }_0\in \mathbb{Z}$ and ${\lambda }_1\in \left(-1/2,1/2\right]$), it can be summarized from Figure 3 and Figure 4 that the reconstructed image gradually blurs with the increase in ${\lambda }_0$ while ${\lambda }_1$ is constant, and the sharpness of the image increases as ${\lambda }_1$ is closer to 0 while ${\lambda }_0$ is constant. The RMSE is minimal only when $\lambda =0$. The phenomenon verifies that the delta filter has the characteristic of minimal reconstruction error.

The simulated digital resolution chart of 1024 × 1024 pixels was reconstructed by basic filters with different $\lambda$ to highlight the phenomenon mentioned above. The number of projections was set to 1440, and the sampling interval of the detector was equal to the pixel size of the phantom. Some results are shown in Figure 5.

By observing the local centers of the reconstructed images shown in Figure 6, it can be noticed that the stripes near center become clearer as ${\lambda }_1$ is closer to 0 while ${\lambda }_0$ is constant, indicating improved resolution.

6. Discussion

The convolution kernel mentioned in this article is basically used for projection data restoration. Meanwhile, the kernel is the inverse Fourier transform of the filtering window function in the frequency domain. In that sense, the proposed method indicates that the applicable window function has more possibilities, and even that the window function does not necessarily have to converge. It shows from a side view that the traditional filter design method does have some limitations.

As a conclusion, any CT image reconstructed through FBP is a linear superposition of the images reconstructed by basic filters. However, there are more properties of basic filters that need to be revealed. Additionally, the optimal composition of the basic filters remains a subject to be studied. This article mainly describes the relevant theories of filter construction, about the rationale extension and basic filter decomposition. The possible applications of the theories need to be further explored, which are not discussed here.

7. Conclusions

The SDBP technique proposed in this article is an approach for an accurate inverse Radon transform, which can be realized through a series of well-defined calculation processes. The transformation between the filter function and convolution kernel is derived from the SDBP formula. The conversion relationship indicates a systematic filter design method, showing that a desired filter can be acquired by an appropriate convolution kernel. By extension, a decomposition method of filters is proposed. For all filters used for CT reconstruction, any of them can be decomposed into basic filters. Some characteristics of the basic filters are demonstrated through simulations of CT reconstruction for the Shepp–Logan phantom and the digital resolution chart. By comparing RMSE and AGM, it is concluded that the error and resolution of the basic filter reconstructed image change periodically with parameter $\lambda$. As an inference, the properties of basic filters can be exploited to construct filters that meet some certain index requirements.