Lisha Wu¹, Shidi Yang^2*, Jing Huang¹, Dongjiang Ji^1*
1School of Science, Tianjin University of Technology and Education, Tianjin, China.
²Tianjin Sanying Precision Instruments Co., Ltd., Tianjin, China.
DOI: 10.4236/jsip.2025.164005   PDF    HTML   XML   60 Downloads   300 Views   Citations

Abstract

In-line X-ray phase-contrast computed tomography (IL-PCCT) can reconstruct high resolution images, and the core of IL-PCCT is the image reconstruction algorithm. As we all know, image reconstruction algorithm can suppress artifacts and produce high-quality reconstructed images using a priori constraint information. For conventional absorption X-ray CT, a common constraint is that linear attenuation coefficient has to be positive, as a negative value is physically not possible. Using the non-negative values priori information is believed to be beneficial in terms of image quality and convergence speed. Phase retrieval is a key step in achieving IL-PCCT. Due to considerations such as radiation dose, single-distance phase retrieval is typically employed. However, single-distance phase retrieval relies on limited data, and especially when the data is insufficient, it may lead to inaccurate reconstruction results. If they exceed the expected range or do not match the physical characteristics of the sample, they generally indicate errors in the retrieval process. In the image reconstruction process, when using projection data obtained from single-distance phase retrieval, whether to introduce a non-negative constraint as a regularization condition is a question worth exploring, but its validity has not yet been systematically studied. In this work, we present an investigation into this problem. The investigation is divided into various experimental studies that non-negativity constraint effect on the IL-PCCT.

Keywords

Non-Negativity Constraint, IL-PCCT, X-Ray

Share and Cite:

1. Introduction

Coherent X-ray imaging is based on wave-optical propagation of electromagnetic waves, including free space propagation and the interaction of short wavelength light with sample. As a partial coherence imaging system, in-line X-ray phase contrast imaging is assumed to be a linear shift-invariant system for objects with weak absorption. Considering an in-line phase-contrast imaging system, the sample can be described as follows:

$n\left(x,y,z\right)=1-\delta \left(x,y,z\right)+i\beta \left(x,y,z\right),$ (1)

where $\left(x,y,z\right)$ denotes the spatial coordinate, the real part $\delta \left(x,y,z\right)$ is related to the phase shift information, and the imaginary part $\beta \left(x,y,z\right)$ describes absorption information. In practical applications, in-line phase-contrast (IL-PC) is combined with CT, and IL-PCCT can reconstruct image with higher resolution and better contrast in biological soft tissues compared with traditional CT (termed as absorption-based CT) imaging techniques. IL-PCCT is a widely used imaging technique for visualizing the internal architecture of biological samples at micrometer resolution [1]. IL-PCCT reconstructs image by exploiting the fact that X-rays are not just absorbed when passing through matter but also refracted, and. it produces a more pronounced contrast compared with conventional absorption-based X-ray imaging with regard to observing biological soft tissues [2]-[4].

The image reconstruction algorithm is the core of IL-PCCT. Image reconstruction algorithms can be classified into analytic and iterative reconstruction algorithms can be generally classified into two major categories: analytical reconstruction and iterative reconstruction. The analytical reconstruction algorithm, such as the filtered back projection (FBP) algorithm, attempts to formulate the solution in a closed-form equation, while an iterative reconstruction algorithm yields a final result as the solution to either a set of solutions or an optimization problem [5]-[7]. For all we know, in order to suppress artifacts and produce high-quality reconstructed images, a solution approach therefore also needs to have regularizing properties. Among these regularizing properties, total variation (TV) has been widely used. Sidky et al. proposed the TV minimization algorithm, which was applied to CT image reconstruction with a positivity constraint [8] [9], and this algorithm is named SART-TV algorithm. The TV model falls into two categories: isotropic TV and anisotropic TV [10]-[12], the SART-TV algorithm is isotropic, the anisotropic total variation (ATV) has two advantages over the isotropic total variation [13], an anisotropic TV minimization method in in-line phase-contrast tomography without a positivity constraint is proposed, and this algorithm is named SART-ATV algorithm.

Besides enforcing regularizing properties of the solution, the appropriate solution should be restricted as far as possible by available physical a priori knowledge: support constraints (real-world samples are of finite size), homogeneous objects (the sample composed of a single material, proportionality of $\delta$ and $\beta$ ), and non-negativity (by the physics of X-rays, the $\delta$ and $\beta$ are always non-negative).

This study is to explore the influence of non-negative constraints on IL-PCCT in the following research contents: (1) The influence of non-negative constraints on analytic reconstruction algorithm; (2) The impact of the non-negativity constraint for simple and complex structures sample; (3) The influence of non-negative constraints on the absorption of background in reconstructed CT images; (4) Non-negativity constraint imposes limitations on the choice of iterative reconstruction algorithm. We view the application of SART, SART-TV and SART-ATV algorithm to research the impact of the non-negativity constraint for phase contrast tomography.

The remainder of this paper is organized as follows. In Section 2, a review of comparison methods are respectively presented, In Section 3, several experiments are carried out with real data to reveal the influence of non-negative constraints on IL-PCCT. Finally, discussion and conclusion are given in Section 4.

2. Method

2.1. Phase-Attenuation Duality Paganin Algorithm

In-line phase-contrast computed tomography (IL-PCCT), a density difference within the sample causes a change in the X-ray phase shift, which has passed through the sample [13]. During the tomography scans, the projection images are collected from various views evenly distributed over $\pi$ -view in BL13W1 of the Shanghai synchrotron radiation facility (SSRF).

A significant proportion of phase-retrieval is performed before IL-PCCT reconstruction. This study uses phase-attenuation duality Paganin (PAD-PA) method [14] to solve the problem of phase retrieval. The PAD-PA method can be described by the following formula:

${\phi }_{\theta }\left(x,y\right)=\frac{1}{2}\gamma \ln \left[{ℱ}^{-1}\left(\frac{ℱ\left(I_{\theta }^D\left(x,y\right)\right)}{\pi \gamma \lambda D\left({\tau }^2+{\upsilon }^2\right)+1}\right)\right]$ (2)

where ${\phi }_{\theta }\left(x,y\right)$ denotes the phase-shift function at the projection angle $\theta$ , $I_{\theta }^D\left(x,y\right)$ represents the intensity function measured by the detector at a distance $D$ , $\lambda$ denotes the X-ray wavelength, $\gamma$ is a constant, and $\left(\tau ,\upsilon \right)$ are the coordinates of $\left(x,y\right)$ in the Fourier domain, and $ℱ$ and ${ℱ}^{-1}$ are the 2D Fourier transform and its inverse operators.

PAD-PA is a phase retrieval technique that relies on single-distance diffraction data. This approach exploits the inherent relationship between phase and attenuation in the diffraction pattern to recover the phase information of the sample, making it especially effective for X-ray diffraction imaging.

2.2. CT Reconstructed Algorithm

A CT imaging system can be modeled as the following discrete linear system after phase retrieval:

$A\delta =\phi ,$ (3)

where $A\in R^{M\times N}$ denotes the CT system matrix, $\phi \in R^M$ denotes the projection data from detector after necessary processing, and $\delta \in R^N$ denotes the phase CT image of the object to be reconstructed. $M$ is the number of detector bins and $N$ is the number of the pixels of the CT image. The task of CT imaging aims to reconstruct $\delta$ from the projection data $\phi$ and the system matrix $A$ .

2.2.1. FBP Reconstructed Algorithm

Apply the filter to the projection data ${\phi }_{\theta }\left(x,y\right)$ to obtain ${\phi }_{\theta }^{filtered}\left(x,y\right),$ the filtered projection data will be back-projected into the image space, with the goal of mapping the intensity distribution of the projections back to the corresponding image locations. The back-projection process is typically implemented through the following integral formula [5]:

$\delta \left(x,y\right)=\underset{0}{\overset{\pi }{{\displaystyle \int }}}{\phi }_{\theta }^{filtered}\left(x,y\right)|d\theta$ (4)

2.2.2. SART Reconstructed Algorithm

SART is an iterative method for image reconstruction that updates all pixels simultaneously in each iteration. It is especially useful for situations with noisy or incomplete data and provides a more efficient approach.

The SART algorithm formula can be expressed as [5]:

${\delta }_j^{\left(K+1\right)}={\delta }_j^{\left(K\right)}+\frac{{\lambda }_K}{{{\displaystyle \sum }}_{i=1}^M{a}_{ij}}{\displaystyle \sum }_{i=1}^M\left[\frac{{\phi }_i-{{\displaystyle \sum }}_{n=1}^N{a}_{in}{\delta }_n^{\left(K\right)}}{{{\displaystyle \sum }}_{n=1}^N{a}_{in}}\right]a_{ij}$ (5)

where the relaxation parameter${\lambda }_{k }$ controls the step size of each update balancing convergence speed and stability.

2.2.3. SART Algorithm Based on Total Variation (SART-TV)

The CT reconstruction problem (3) can be solved by the CS-based reconstruction method to minimize the image isotropic TV regularized by the projections. This process is equivalent to solving the following optimization program:

$\delta =\underset{\delta }{argmin}{\left(A\delta -\phi \right)}_2^2+{\lambda }_1^{*}{\delta }_{TV}^{iso}$ (6)

The model described by Equation (4) was investigated using the alternating minimization procedure [8]:

1: a SART step to solve (3);

2: descending the TV-gradient to minimize the image TV. The gradient descent is controlled via the step size parameter ${\lambda }_1^{*}$ , If the step size of the gradient descent is too strong the image becomes uniform and inconsistent with the projection data. On the other hand, if the step size of the gradient descent is too small, the algorithm reduces to standard SART with a positivity constraint included. The SART-step loop is executed $K_{max}$ times. $N_{grad}$ denotes total number of gradient descents.

2.2.4. SART Based on Anisotropic Total Variation (SART-ATV)

The CT reconstruction problem given in Equation (1) can be solved to minimize the image anisotropic TV regularized. This process is equivalent to solving the following optimization program:

$\delta =\underset{x}{argmin}{\left(A\delta -\phi \right)}_2^2+{\lambda }_2^{*}{\delta }_{TV}^{ani}$ (7)

Equation (5) is decoupled data fidelity minimization and image regularization problems based on proximal forward–backward splitting (PFBS), the SART algorithm is used to solve the data fidelity. ADMM is used to solve the image regularization, The ADMM-step loop is executed $N_2$ times, ${\lambda }_2^{*}$ is the regularization parameter [13].

2.3. Quantitative Assessment

Three metrics, the MSE, the PSNR, and structural similarity index (SSIM), are used to evaluate the quality of reconstructed images.

The PSNR is a traditional measure of image quality, and a larger PSNR value represents better quality. It is defined as follows:

$MSE\left(x,y\right)=\frac{1}{M\times M}{\displaystyle \sum }_{i=1}^M{\displaystyle \sum }_{j=1}^N{\left(x_{i,j}-y_{i,j}\right)}^2$ (8)

$PSNR\left(x,y\right)=10{\log }_{10}\left(\frac{Pea{k}^2}{MSE}\right)$ (9)

where MSE is the mean square error function, x denotes the reference image, y denotes the reconstructed image, and $M\times M$ is the size of $x$ and $y$ ; $x_{i,j}$ and $y_{i,j}$ represent the pixel intensities of $x$ and $y$ in some pixel $\left(i,j\right)$ , respectively, and Peak represents the largest pixel intensity in the normalized image.

We use the structural similarity (SSIM) index [15] as a quantitative measure for the reconstructed image quality. The computed local similarity index is defined on windows $x$ and $y$ :

$ssim\left(x,y\right)=\frac{\left(2{\mu }_x{\mu }_y+c_1\right)\left(2{\sigma }_{xy}+c_2\right)}{\left({\mu }_x^2+{\mu }_y^2+c_1\right)\left({\sigma }_x^2+{\sigma }_y^2+c_2\right)}$ (10)

where ${\mu }_x$ and ${\mu }_y$ are the averages of $x$ and $y$ ; ${\sigma }_x^2$ , ${\sigma }_y^2$ are the variances; ${\sigma }_{xy}$ is the covariance of $x$ , $y$ ; and $c_1$ , $c_2$ are two variables to stabilize the division with weak denominator. The overall SSIM is the mean of local similarity indices:

$SSIM\left(X,Y\right)=\frac{1}{N}{\displaystyle \sum }_{i=1}^{N}ssim\left(x_i,y_i\right)$ (11)

where $X$ is a reference image; $Y$ is a reconstructed image; $x_i,y_i$ are the corresponding windows indexed by $i$ ; and $N$ is the number of windows. Here, the size of the windows is $8\times 8$ .

3. Numerical Results

The results are separated into two subsections. Each subsection provides the reconstructed images by analytic and iterative reconstruction algorithms based on the approach with and without non-negativity constraint. Projection data from coaxial contrast CT imaging was acquired at SSRF’s BL13W1 line station. The X-ray energy was set to 33 keV, and the SDD was adjusted to 0.8 m.

3.1. The Impact of the Non-Negativity Constraint for Analytic Reconstruction Algorithm

Figure 1. Phase contrast reconstruction results of the sample using FBP reconstruction of the bone model II. (a) No constraints phase information; (b) non-negative phase information; (c) Negative phase information. Enlarged phase information images of (a, b, c); (d, e, f) The enlarged images were from the same regions in the reconstructed images, as marked with the green rectangle in Figure 1(a). The grayscale images were normalized to the range [0, 255], and the display window was [0, 255].

For conventional absorption X-ray CT, a common constraint is that linear attenuation coefficient has to be positive, as a negative value is physically not possible. However, for the analytic reconstruction algorithm, algorithm, such as the filtered back projection (FBP) algorithm, because the ramp filter zeros the dc term in the frequency domain, each back projection image will have zero average value. This means that the pixels in each back projection image will have negative and positive values. When all the back projections are added to form the final image, some negative locations may become positive and the average value may not be zero, but typically, the final image will still have negative pixels [16].

In order to research the impact of the negative and positive values for result image of FBP algorithm, no constraints phase information image (negative and positive values), positive phase information image, and negative phase information image are reconstructed by FBP algorithm. The result of those images are shown in Figure 1.

Figure 1(a) shows the reconstructed image result of negative and positive values. It can be seen from the reconstruction result that all structures of the reconstructed sample have been completely reconstructed. Figure 1(b) are the reconstructed images of positive refractive index. It can be seen from these images and the local enlarged images that part of the tissue structure of the sample area has not been completely reconstructed. It can be seen from Figure 1(b) that the negative refractive index is meaningful, so we reconstructed the negative refractive index. The negative refractive index reconstruction results are shown in Figure 1(c). From these images and their local magnified view, it can be seen that the negative refractive index corresponds to some organizational structures of the object.

3.2. The Impact of the Non-Negativity Constraint for Simple and Complex Structures Sample

In this section, we will analyze the influence of non-negative constraints on reconstruction results for two groups of different samples. The first sample is liver sample, the sample structure is relatively simple, mainly includes two parts, the liver region and vascular space from Figure 2(b) and its local amplification area in Figure 2(e), you can see that for homogeneous sample, nonnegative constraints will not be lost in Figure 2(a) and its local amplification area in Figure 2(d) samples of liver and vascular structure information. As can be seen from Figure 2(c) and its enlarged view, the negative refractive index reconstruction results correspond to the vascular region. Therefore, the structure information of the reconstructed object will not be lost when the simple structure sample is subjected to non-negative constraint.

Figure 2. Phase contrast reconstruction results of the sample using FBP reconstruction. (a) No constraints phase information; (b) Non-negative phase information; (c) Negative phase information; (d, e, f) Enlarged phase information images of (a, b, c). The enlarged images were from the same regions in the reconstructed images, as marked with the green rectangle in Figure 4(a). The grayscale images were normalized to the range [0, 255], and the display window was [0, 255].

To further analyze the influence of non-negative constraints on the reconstructed images of samples with complex structures, we also selected the reconstructed image results of bone trabecular samples for analysis. Figure 3(a) did not carry out the reconstruction results of non-negative constraints, and the reconstructed image contained complete structural information of the sample. The results of Figure 3(b) and local amplification Figure 3(e) show that the important structural information of the sample will be lost when the reconstructed image with non-negative constraint is applied to the complex structure sample. It can be seen from Figure 3(c) and its local amplification Figure 3(f), that the negative refractive index also contains important structural information of the sample with complex structure.

Figure 3. Phase contrast reconstruction results of the sample using FBP reconstruction. (a) No constraints phase information; (b) Non-negative phase information; (c) Negative phase information; (d, e, f) Enlarged phase information images of (a, b, c). The enlarged images were from the same regions in the reconstructed images, as marked with the green rectangle in Figure 4(a). The grayscale images were normalized to the range [0, 255], and the display window was [0, 255].

3.3. The influence of Non-Negative Constraints on the Absorption of Background in Reconstructed CT Images

For conventional absorption CT, a simple constraint is that every linear attenuation coefficient has to be positive, as a negative value is physically not possible. This condition can be implemented directly in the back projection step of the iterative method by a simple truncation. The application of this constraint results in a reduction of noise in the reconstructed image, especially in the background. This is illustrated in Figure 4(a) and Figure 4(c) (enlarged image in Figure 4(a)), where one can clearly see how several noise in the background reconstructed image reconstructed images by SART algorithm. However, several noise can no longer be observed when applying the nonnegative constraint, this is illustrated in Figure 4(b) and Figure 4(d) (enlarged image in Figure 4(b)).

Figure 4. Reconstructed images of the iron sample using the (a) No constraints; (b) Non-negative constraints; (c, d) Enlarged images of (a, b).

3.4. The Impact of the Non-Negativity Constraint for Iterative Reconstruction Algorithm

It is common to enforce non-negative values to accommodate the physical non-negativity of x-ray attenuation. However, non-negativity constraint imposes limitations on the choice of iterative optimization algorithm for IL-PCCT. Compared with the analytic reconstruction algorithm, the main advantages of iterative methods are that regularization techniques are often used in iterative reconstruction algorithms. The choice of regularization term depends on the specific object examined and, on the desire, to preserve or emphasize particular features.

First, we view the application of SART., SART-TV, and SART-ATV algorithm to research the impact of the non-negativity constraint for IL-PCCT. Under the few-view condition, the images with a size of 1821 × 1821 pixels were reconstructed by these algorithms. The parameters were set as follows:

1) SART: $\mu =0.05$ and $K_{\text{max}}$ ; 2) SART-TV: ${\lambda }_1^{*}=0.02$ , $N_{grad}=10$ ; 3) SART-ATV: ${\lambda }_2^{*}=0.00002$ , $N_2=20$ . All parameters were obtained by optimal selection.

Figure 5(a) is assumed to represent a reference image reconstructed by the SART algorithm using 1821 projections. To study the few-view problem, 455 projections were used to achieve the reconstruction in Figures 5(b)-(d), Figures 5(e)-(h) are the enlarged images from the same regions in the reconstructed images. Figure 5(c) shows that import structures of the sample area have not been completely reconstructed by SART-TV including a positivity constraint. Figure 5(b) and Figure 5(d) show that the complete structural image reconstructed by SART and SART-ATV without a positivity constraint. Compared with the slice images in Figure 5(b) and Figure 5(d), the image is reconstructed better and streak artifacts are suppressed by the SART-ATV algorithm, while the SART reconstruction image is contaminated by streak noise.

Figure 5. Reconstruction results of the bone sample using different algorithms. The images were reconstructed by (a)SART for full-view; (b) SART for few view; (c) SART-TV for few view; (d) SART-ATV for few view; (e, f, g, h) Enlarged images of (a, b, c, d).

To better quantitively value the reconstruction results by different algorithm, the PSNR, SSIM, MSE and TIME are listed in Table 1. We see that the reconstructed image obtained by SART-ATV algorithm is of the best quality.

Table 1. Quantitative results regarding the reconstructed phantom images.

	Methods	PSNR (dB)	SSIM	MSE	TIME (s)
Few-view	SART	24.51	0.7891	15.17	120
	SART-TV	12.02	0.4465	63.89	124
	SART-ATV	25.31	0.9112	13.84	180

Note: The reconstruction time was measured based on a reconstruction slice.

4. Discussion and Conclusion

Single-distance phase retrieval relies on limited data, and when the data is insufficient, it may lead to inaccurate reconstruction results. In such cases, the choice of non-negative constraints plays a crucial role in influencing the quality of image reconstruction. This study demonstrates the impact of non-negative constraints on IL-PCCT through several experiments. First, we analyze the effect of non-negative constraint truncation on the analytic reconstruction algorithm. The experimental results show that applying non-negative constraints during sample reconstruction leads to the loss of important structural information. Second, for IL-PCCT, we investigate whether applying non-negative constraints to samples with different structures results in the loss of important structural information. We selected two samples with different complexities for experimental analysis. The results show that for simple structure samples (e.g., liver samples), whether or not non-negative constraints are applied, the structural information of the liver and vascular regions is fully restored. However, when non-negative constraints are applied to complex structure samples, important structural information may be lost. Third, we analyze the impact of non-negative constraints on conventional absorption CT. In conventional absorption CT reconstruction, the linear attenuation coefficient must be positive, as negative values are physically impossible. By directly implementing this constraint in the back projection step of the iterative method, noise in the reconstructed image, especially in the background, can be reduced. Finally, we analyze the effect of non-negative constraints on iterative reconstruction algorithms. The experimental results show that, in sparse angle reconstructions, the SART-ATV algorithm, without non-negative constraint priors, is able to fully reconstruct all structural information of the sample and effectively reduce streak artifacts. The reconstruction results are significantly better than those obtained with the traditional SART algorithm and the SART-TV algorithm, which requires the introduction of non-negative constraints.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Coan, P., Bamberg, F., Diemoz, P.C., Bravin, A., Timpert, K., Mützel, E., et al. (2010) Characterization of Osteoarthritic and Normal Human Patella Cartilage by Computed Tomography X-Ray Phase-Contrast Imaging. Investigative Radiology, 45, 437-444.[CrossRef] [PubMed]
[2]	Snigirev, A., Snigireva, I., Kohn, V., Kuznetsov, S. and Schelokov, I. (1995) On the Possibilities of X-Ray Phase Contrast Microimaging by Coherent High-Energy Synchrotron Radiation. Review of Scientific Instruments, 66, 5486-5492.[CrossRef]
[3]	Davis, T.J., Gao, D., Gureyev, T.E., Stevenson, A.W. and Wilkins, S.W. (1995) Phase-contrast Imaging of Weakly Absorbing Materials Using Hard X-Rays. Nature, 373, 595-598.[CrossRef]
[4]	Momose, A. (2020) X-Ray Phase Imaging Reaching Clinical Uses. Physica Medica, 79, 93-102.[CrossRef] [PubMed]
[5]	Kak, A.C. and Slaney, M. (2001) Principles of Computerized Tomographic Imaging. Society for Industrial and Applied Mathematics.[CrossRef]
[6]	Hounsfifield, G.N. (1972) A Method and Apparatus for Examination of a Body by Radiation Such as X-Ray or γ Radiation. Patent Specifification 1283915.
[7]	Gordon, R., Bender, R. and Herman, G.T. (1970) Algebraic Reconstruction Techniques (ART) for Three-Dimensional Electron Microscopy and X-Ray Photography. Journal of Theoretical Biology, 29, 471-481.[CrossRef] [PubMed]
[8]	Sidky, E.Y., Kao, C. and Pan, X. (2006) Accurate Image Reconstruction from Few-Views and Limited-Angle Data in Divergent-Beam Ct. Journal of X-Ray Science and Technology: Clinical Applications of Diagnosis and Therapeutics, 14, 119-139.[CrossRef]
[9]	Sidky, E.Y. and Pan, X. (2008) Image Reconstruction in Circular Cone-Beam Computed Tomography by Constrained, Total-Variation Minimization. Physics in Medicine and Biology, 53, 4777-4807.[CrossRef] [PubMed]
[10]	Rudin, L.I., Osher, S. and Fatemi, E. (1992) Nonlinear Total Variation Based Noise Removal Algorithms. Physica D: Nonlinear Phenomena, 60, 259-268.[CrossRef]
[11]	Cai, J., Dong, B., Osher, S. and Shen, Z. (2012) Image Restoration: Total Variation, Wavelet Frames, and Beyond. Journal of the American Mathematical Society, 25, 1033-1089.[CrossRef]
[12]	Bui, K., Park, F., Lou, Y. and Xin, J. (2021) A Weighted Difference of Anisotropic and Isotropic Total Variation for Relaxed Mumford—Shah Color and Multiphase Image Segmentation. SIAM Journal on Imaging Sciences, 14, 1078-1113.[CrossRef]
[13]	Ji, D., Qu, G., Hu, C., Liu, B., Jian, J. and Guo, X. (2017) Anisotropic Total Variation Minimization Approach in In-Line Phase-Contrast Tomography and Its Application to Correction of Ring Artifacts. Chinese Physics B, 26, Article ID: 060701.[CrossRef]
[14]	Barrett, R., Berry, M., Chan, T.F., Demmel, J., Donato, J., Dongarra, J., et al. (1994) Templates for the Solution of Linear Systems: Building Blocks for Iterative Methods. Society for Industrial and Applied Mathematics.[CrossRef]
[15]	Wang, Z., Bovik, A.C., Sheikh, H.R. and Simoncelli, E.P. (2004) Image Quality Assessment: From Error Visibility to Structural Similarity. IEEE Transactions on Image Processing, 13, 600-612.[CrossRef] [PubMed]
[16]	Gonzalez, R.C. and Woods, R.E. (2017) Digital Image Processing. 4th Edition, Electronic Industry Press., 368-393.