Poisson-Gaussian Noise Reduction for X-Ray Images Based on Local Linear Minimum Mean Square Error Shrinkage in Nonsubsampled Contourlet Transform Domain

Noise reduction is important for X-ray images because it can reduce radiation exposure to patients. X-ray image noise has a Poisson-Gaussian distribution, and recently, noise analysis and removal in multiscale transformations have been widely implemented. The nonsubsampled contourlet transform (NSCT) is a multiscale transformation suitable for medical images that separates the scale and direction. This study proposes a Poisson-Gaussian noise-removal method using NSCT shrinkage that is based on the characteristics of Poisson-Gaussian noise in NSCT domain. It has the structure of a block-matching 3D filtering algorithm in the form of basic estimation and noise removal process; however, the main processes are modified to consider Poisson-Gaussian noise characteristics. In the basic estimation process, an NSCT shrinkage method that is suitable for Poisson-Gaussian noise characteristics is developed by optimizing the local linear minimum mean square error estimator in the NSCT domain. In the denoising step, the noise term of the Wiener filter is determined using the result of the NSCT shrinkage, and finally, the denoised image is obtained. The proposed method is applied to simulated and real X-ray images and is compared with other state-of-the-art Poisson-Gaussian noise removal methods; it exhibits excellent results in both quantitative and qualitative aspects.


I. INTRODUCTION
Image acquisition devices acquire digital images by converting light into electrical signals via complementary metal oxide semiconductors or charged coupled devices image sensors [1], [2]. The intensity of the image is determined by the number of photons incident on the sensor, and the number of photons follows a Poisson distribution. This is the primary cause of noise in the image, which also has Poisson characteristics. However, the noise of an image acquired in a general situation (medium, high-illumination situation) primarily follows a Gaussian distribution because the number of photons incident on the sensor is sufficiently large; therefore, The associate editor coordinating the review of this manuscript and approving it for publication was Shiqi Wang. the Poisson distribution is approximated to a Gaussian distribution. However, in low-light, ultra-low-light conditions, or when the number of photons incident on the sensor is small, the image noise cannot be approximated by a Gaussian distribution, and it follows the Poisson distribution [3], [4]. In addition to the Poisson noise generated during the image acquisition process, noise from the electric equipment also appears in the image. Noise generated by sensors or other electronic devices follows a Gaussian distribution, and it is modeled using a combination of Poisson and Gaussian distribution, called a Poisson-Gaussian distribution [5].
X-ray image noise follows a Poisson-Gaussian distribution [6], which causes the image to be degraded; thus, a method to overcome this degradation is necessary. There are two ways to increase the signal-to-noise ratio (SNR) of an VOLUME 9, 2021 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ image. The first method is to increase the signal power while maintaining the intensity of the image noise. The second method is to reduce the image noise while maintaining the signal power. For the first method, the power of the signal can be increased by increasing the photon energy using stronger X-ray generator or by increasing the number of photons via longer exposure time. However, X-rays use high-energy; they can pass through the human body to obtain an image inside the body, and these high-energy photons expose the human body to radiation. Thus, these methods increase the amount of radiation exposure, which may adversely affects the human body. Therefore, alternative methods are required to obtain a clean X-ray image while reducing the amount of exposure [7]. The second method to increase the SNR is to reduce the power of the noise while maintaining the power of the signal, i.e., noise removal. Noise removal has been extensively studied, and traditional denoising methods mainly focus on removing additive Gaussian noise. Robbins proposed an empirical Bayesian framework to estimate Gaussian noise [8], and Lee proposed a two-step empirical Bayesian estimation [9], which estimates the signal variance from local statistics and applies a minimum mean square error estimator to obtain the noise filtering algorithms. Buades et al. proposed nonlocal means (NLM) filtering [10], which is a method to remove noise by calculating the weighting coefficients of neighboring pixels using the similarity between patches. Block-matching and 3D filtering (BM3D) [11], which is a state-of-the-art technology, effectively removes noise by grouping similar 2D patches into 3D data and then implementing thresholding and collaborative Wiener filtering.
Poisson noise is signal dependent: therefore, it is difficult to apply an additive Gaussian noise removal method because it does not have a constant noise variance. To solve this problem, variance stabilization transformation (VST), such as Anscombe transformation [12] or Fisz transformation [13] has been introduced. Signal-dependent Poisson noise that has undergone the Anscombe transform has a constant variance; thus, noise can be removed using a Gaussian noise removal algorithm. The image from which Poisson noise has been removed can be obtained through an inverse transform [14]- [16]. The Fisz transform can be combined with the Haar wavelet [17], and they are widely used as a Haar-Fisz transform [18]. Moreover, because these VSTs take into account when noise has a Poisson distribution, a generalized Anscombe transform (GAT) was introduced [19] that developed Anscombe transform into Poisson-Gaussian noise. However, because these transforms are nonlinear transforms, bias errors occur during inverse transforms. Therefore, noise removal may not be performed properly or images may be damaged [20]- [22]. To overcome this problem, Makitalo and Foi proposed the exact unbiased inverse of the GAT [23]- [26].
There are also techniques to remove noise using the signal-dependent noise characteristic without VSTs. Kuan et al. proposed a local linear minimum mean square error (LLMMSE) filter [27] that estimates the variance of signal and noise using the local mean and variance; they also design a noise removal filter using the estimated value. Le et al. proposed a noise removal method based on total variation normalization considering Poisson noise [28], and Bindilatti and Mascarenhas proposed an NLM filter that uses stochastic distances for Poisson noise instead of Euclidean distances for Gaussian noise when calculating the similarity between patches [29]. Another way to effectively remove Poisson noise is to use multiscale transformation. The wavelet transform [30], which is a representative multiscale transform, divides images into vertical and horizontal directions for each band. An expectation-maximization image restoration technique based on a maximum penalized likelihood estimator [31] in the wavelet domain was proposed. A noise removal method based on Poisson-Gaussian unbiased risk estimator (PG-URE) [32], [33] is also performed in the wavelet domain, in which Stein's unbiased risk estimator [34] is extended to Poisson-Gaussian noise. The wavelet transform is excellent for separating scales; however, it does not consider directions other than vertical and horizontal. Thus, it is ineffective for separating curves included in natural images. To overcome this limitation, ridgelet [35], curvelet [36], and contourlet transforms (CT) [37], [38] that considered scale and directionality were introduced.
As mentioned above, many studies have been conducted on methods that remove Poisson noise; however, there have been few studies on the removal Poisson-Gaussian noise using its own characteristics without using VSTs. Because X-ray images, particularly on low-dose condition, have Poisson-Gaussian noise [5], it is necessary to effectively remove them. Medical images, including X-rays, are primarily composed of curves; therefore, the nonsubsampled contourlet transform (NSCT), which is a type of CT, effectively separates the images into their scale and direction among multiscale transformations. In our previous study, we performed Poisson-Gaussian noise analysis in the NSCT domain [39]. Based on the analyzed results, an LLMMSE-based shrinkage method is applied through an NSCT, and a BM3D-based noise removal method is introduced. Experiments are conducted using simulated images with artificially generated noise as well as images acquired with actual X-rays, and the denoising results are compared with those of other recent studies. The main contributions of this paper are as follows: 1) the NSCT is used to separate bands suitable for curves; 2) LLMMSE-based NSCT shrinkage is developed based on the noise relationship between the low-band layer and the detail layer of the same scale level; 3) finer level noise is removed by inheriting the denoised result of coarser level; 4) the noise removal performance is maximized using BM3D-based method. The remainder of this paper is organized as follows. In Section II, we briefly review the Poisson-Gaussian noise analysis in NSCT that we performed in our previous study and the noise reduction method, which is the basis of the proposed method. In Section III, we describe the proposed LLMMSE-based inherited NSCT shrinkage noise removal method. Finally, we demonstrate the performance of the proposed method on both simulated and real X-ray images in Section IV, and we present our conclusion in Section V.

II. RELATED WORKS
This section summarizes the key ideas referenced in this study from previously studies. We briefly describe Poisson-Gaussian noise analysis in the NSCT domain, LLMMSE filtering, and BM3D filtering. The reader is referred to [11], [27], [39] for further details. In our previous study, we established a Poisson-Gaussian noise model and analyzed its suitability using real X-ray images acquired under low-dose conditions [39]. The observation model of X-ray and Poisson-Gaussian noise model are expressed as follows [5]: where (i, j) are spatial coordinates, y (i, j) is the acquired image, x (i, j) is the original image, η (x (i, j)) is the standard deviation of the Poisson-Gaussian noise, δ is the independent Gaussian noise with zero mean and a standard deviation equal to one, α is the Poisson noise parameter, and β is the standard deviation of the Gaussian noise. A large number of stationary X-ray images were acquired to confirm that model (1) is suitable. The average and variance of the acquired images were assumed to be a noiseless image and the noise variance of images, respectively. These two images were then compared and analyzed. As a result, the noise variance exhibited a linear relationship with the signal power, which can be expressed as a linear equation (1). Thus, it was demonstrated that the observation and noise model were properly modeled.

2) NONSUBSAMPLED CONTOURLET TRANSFORM
The CT is a transformation that separates bands by scale and direction using a Laplacian pyramid (LP) [40], [41] and directional filter bank (DFB) [42]. As shown in Fig. 1(a), the process of the CT is to first use an LP to separate the scale and then use a DFB to separate the band-pass and the high-pass subbands by direction, excluding the low-pass subband. However, the LP and DFB used in CT both contain a subsampling process; thus, the size of the subbands after transformation is different from the original image size. That is, the CT has spatial variant property, which is unsuitable for noise removal. NSCT [43] was proposed to overcome these shortcomings. The NSCT replaces the LP and DFB of CT with nonsubsampled pyramid (NSP) and nonsubsampled DFB (NSDFB), respectively, to perform multiscale, multidirectional, and shift-invariant image decomposition. The decomposition process of the NSCT is shown in the Fig. 1(b). We describe the NSCT separation process in detail. First, the image is separated into a low-frequency band and a high-frequency band by using the low-pass (H 0 (z)) and high-pass (H 1 (z)) filters of the NSP. The separated high-frequency band is used as a detail layer, and the low-frequency band passes through a low-pass (H 0 z 2I ) and high-pass filter (H 1 z 2I ) of next level, and it is divided into lower-and band-frequency regions. The separated band-frequency region becomes the detail layer of the next level, and the low-frequency band is repeatedly used for the next level of scale transformation. The filter that separates the m-level detail layer is referred to as m , and it is expressed as follows: Assuming that the image is divided into a total of M scale levels, each scale level is denoted by m (m = 1, 2, · · · , M ), where m = 1 and m = M are the finest and the coarsest levels, respectively. The decomposed low-band layer is denoted as G (·), the detail layer is denoted as L (·), the original image VOLUME 9, 2021 is denoted as G 0 (·), and the low-band and detail layers at the m-th scale level is denoted as G m (·) and L m (·), respectively.
Next, the detail layer, the scale of which is separated by NSP, is separated for each direction by the NSDFB. The number of directions separated by NSDFB is expressed as N m , which must be a power of two. Each direction level is denoted by n (n = 1, 2, · · · , N m ), and the n-th direction subband of the m-th scale is denoted by L m,n (·).

3) POISSON-GAUSSIAN NOISE ANALYSIS IN NSCT DOMAIN
In our previous study, to analyze Poisson-Gaussian noise in the NSCT domain, we conducted a step-wise analysis of the noise distribution through the NSP and NSDFB. First, as in the Section II-A1, 100 still images were decomposed by the NSP, and the Poisson-Gaussian noise in the NSP domain was analyzed using the detail layer variance. This shows that Poisson-Gaussian noise still maintains its characteristics, even after NSP decomposition. The noise of the NSP detail layer has a Poisson-Gaussian distribution, which is dependent on the low-band layer at the same level and is expressed as following equations: where η m (G m (x (i, j))), α m , and β m are the standard deviation of the noise, Poisson noise parameter, and standard deviation of the Gaussian noise of NSP detail layer of m-th scale, respectively. At the same m-th scale, α and α m as well as β 2 and β 2 m have the same ratio, E m is equal to the energy of the filters passing through the NSP decomposition. The energy of these filters can be calculated as E m = m (z) 2 2 . Next, the noise distribution after applying the NSDFB was analyzed, and finally, the Poisson-Gaussian noise in the NSCT domain was analyzed. The analysis was performed in the same manner as the noise analysis process in NSP, and NSDFB decomposition was also demonstrated to retain the Poisson-Gaussian distribution when decomposing Poisson-Gaussian noise. Therefore, it was confirmed that Poisson-Gaussian noise has a Poisson-Gaussian distribution even after NSCT decomposition. The noise of the NSCT detail layer is also dependent on the NSP low-band signal, as shown in the following equations: where η m,n (G m (x (i, j))) is the standard deviation of the noise of NSCT subband of m-th scale and n-th direction, and α m,n and β m,n are the Poisson noise parameter and the standard deviation of the Gaussian noise of the NSCT detail layer of m-th scale and n-th direction, respectively. NSDFB also changes the noise parameter at the same rate for the same direction level and distributes 1/N m equally for the same scale level. This can be expressed by the formula: Because the filter of NSDFB distributes the entire band equally, the energy of the filter in each direction is 1/N m , and the noise parameters are divided at an equal ratio.
The minimum mean square error (MMSE) estimate of x, given observation y, is the conditional mean estimate of In boldface, the signal is expressed in the form of a vector, and E (x) means the ensemble mean of x. In general, the MMSE estimate is nonlinear, and it depends on the probability density function of x and η; therefore, in most cases, it is difficult to obtain an explicit form of the MMSE estimator. Thus, the following linear minimum mean square error (LMMSE) [44] estimator is used by applying a linear constraint to the structure of the estimator: where C xy is the cross-covariance matrix of x and y, and C y is the covariance matrix of y. Kuan et al. proposed the LLMMSE noise removal filter [27] for images without degradation, except for nonstationary mean and nonstationary variance. In the case where there is no degradation owing to blurring and the noise is uncorrelated, (8) is expressed by the following scalar processing filter: (9) where σ 2 x (i, j) and σ 2 η (i, j) are the ensemble variances of the original image and nonstationary noise, respectively. Because the noise is assumed to be nonstationary, the ensemble statistics can be replaced with local spatial statistics; thus, (9) can be written aŝ where x (i, j) and y (i, j) are the local means of x (i, j) and y (i, j), respectively, and v x (i, j) is the local spatial variance of x (i, j). The local mean and variance can be calculated using the uniform moving average window of size (2r + 1)× (2r + 1). Then, and The method used to obtain σ 2 η (i, j) differs according to the type of noise, and the LLMMSE filter can be designed using the relationship between v x (i, j), v y (i, j), and σ 2 η (i, j).
2) BLOCK-MATCHING AND 3D FILTERING BM3D filter [11] is one of the most advanced noise removal methods, and it is excellent for Gaussian noise removal. It has the features of LMMSE, nonlocal means, and transformation domain-based filtering, and it works by synthesizing them efficiently. It consists of two steps: a basic estimation step and a denoising step. In the basic estimation step, noise is coarsely removed to obtain a relatively clean image, and in the denoising step, reliable statistics are obtained using basic estimation, and the actual noise removal is performed. All processes are performed block-wise, similar to nonlocal approach, and not pixel-wise using neighboring pixels. Filtering is performed in the transform domain, and the final result is derived through the Wiener filter, which is one of the optimum LMMSE estimators. Each step of BM3D consists of three stages: grouping, filtering, and aggregation. In the grouping stage, blocks that are similar to the reference block are identified in the image and grouped in three dimensions. Because the noise of the image is assumed to be Gaussian noise, blocks around the reference block are searched based on the Euclidean distance. In the basic estimation step, a noisy image is used for searching and grouping. In the denoising step, the search uses the basic estimation result of the first step, and in the grouping stage, both the noisy image and basic estimated image in the same position are grouped. In the filtering stage, a 3D transform is applied to the generated group, noise is removed in the transform domain, and an inverse transform is then performed. For 3D transformation, a biorthogonal spline wavelet or discrete cosine transform (DCT) is typically used. Because similar blocks are grouped together, the 3D transform domain has a sparse property. Using this property, noise is coarsely removed via hard thresholding in the basic estimation step, and it is effectively removed by applying a Wiener filter in the denoising step. Finally, in the aggregation stage, each block from which noise has been removed is returned to its original position, and the pixel value is determined using a weighted average. Because this is a block-wise process,, if there are many similar blocks, there are many overlapping blocks at the same pixel position. The pixel value is determined by the weighted average of the blocks that overlap at one location. The weight is inversely proportional to the noise variance; therefore, so when the noise is strong, a small weight is given, and when the noise is weak, a large weight is given.
Because BM3D is a Gaussian noise removal filter, it is not suitable for Poisson-Gaussian noise removal. However, BM3D can be applied by converting Poisson-Gaussian noise into Gaussian noise using VST. Makitalo and Foi proposed GAT [24]- [26], the VST of Poisson-Gaussian noise, and they confirmed that the algorithm that combines GAT and BM3D performs better than other methods, as described in [26].

III. PROPOSED METHOD
This section introduces the removal of Poisson-Gaussian noise by applying LLMMSE-based shrinkage in the NSCT domain and modifying the BM3D algorithm. Because there are more curves than straight lines or patterns in medical images, it is necessary to use a transformation that is suitable for this condition. To remove noise, a transformation with a shift-invariant feature must be used; therefore, NSCT was used among many multiscale transformations. LLMMSE shrinkage in the NSCT domain was proposed based on the Poisson-Gaussian noise analysis in the NSCT domain. We briefly explain the proposed algorithm. The input image y (i, j) is degraded by Poisson-Gaussian noise. It is first separated into its scale and direction components through NSCT decomposition. From the subbands of the coarsest level, noise is removed by LLMMSE shrinkage. After removing the noise of all subbands at the same scale level, the low-band of the finer level is reconstructed using NSDFB and NSP reconstruction. The reconstructed low-band from which noise has been removed is used again for the next finer level LLMMSE shrinkage, and this process is repeated up to the finest level to obtain an image from which the noise has been primarily removed. Using this image as a basic estimation, the modified BM3D algorithm is performed with the noisy input image to obtain a clean image with the noise removed. An overall block diagram of the proposed algorithm is shown in Fig. 2.

A. LLMMSE BASED NSCT SHRINKAGE
Low-dose X-ray images have Poisson-Gaussian noise as shown in (1). The input image is first separated into subbands for each scale and directional band using NSCT decomposition. Because there is no low-frequency component in the subband, E L m,n (x (i, j)) = E L m,n (y (i, j)) = 0; thus, the LLMMSE filter in (10) becomes where r m is the size of the moving window. Because NSCT does not have a subsampling process; therefore, to obtain accurate local statistics, the size of the window needs to VOLUME 9, 2021 increase as the level increases. Thus, the window size was applied differently according to the scale level as r m = 2 m . We set up the noise model in (1) in the form of an additive stationary and this property does not change after NSCT decomposition, as shown in (5); therefore, the relationship between v L m,n (y) (i, j) and v L m, The variance of the noise σ 2 L m,n (η) (i, j) can be estimated based on the Poisson-Gaussian noise analysis in the NSCT domain using (5), and the noise parameters α m,n , β 2 m,n can be calculated using (4) and (6), as follows: Substituting (14)- (17) into (13) yields the LLMMSE shrinkage filter as per the following equation: .
Both the numerator and denominator must have a value of zero or higher in (13). However, there is a possibility that the local variance of the subband may have a smaller value than the estimated power of noise; taking this into consideration, the final shrinkage type equation is obtained as follows: The problem is now reduced to the estimation of G m (x (i, j)) in the aforementioned formulas. In (17), the noise parameter of each subband decreases as the level coarsens. The noise that remains after being distributed to each subbands exists in the low-band of the coarsest level, and the ratio of the low-band noise to the total noise can be calculated as 1 − M m=1 E m . Using this to calculate the noise ratio of the low-band of the coarsest level when M = 3, for E 1 = 0.769155, E 2 = 0.178771, and E 3 = 0.044181, a value of 0.007893 can be obtained. That is, when M ≥ 3, the noise existing in the low-band of the coarsest level is less than 1% of the total noise, and it can be assumed that this is a noiseless image, The low-bands of noiseless and noisy images for M = 1, 2, and 3 are shown in Fig. 3. As M increases, the noise appears less in the low-band of the noisy image and appears similar to the low-band of the noiseless image. It is recommended that M ≥ 3 is set for the noise reduction performance.
Substituting (20) into (19), the noise of the subbands of the coarsest level can be removed, and the denoised low-band of the finer level is then obtained by NSDFB and NSP reconstruction. The denoised low-band is used to estimate the noise variance of the finer level, and this process is repeated up to the finest level to obtainx shr (i, j) with the noise removed.

B. MODIFIED BM3D
The basic estimation, which is the result of the first step of BM3D, was obtained using NSCT shrinkage; thus, the second step of the BM3D algorithm is now performed. In the grouping stage, similar blocks are identified by calculating the difference between the blocks inx shr (i, j). Because noise has been removed once inx shr (i, j), the distance between the blocks is calculated using the Euclidean distance. Let X i,j be a W B × W B size block extracted from x (i, j), where the subscripts i, j are the coordinates of the top-left corner. When the top-left coordinate of the reference block is (i R , j R ), the distance between the reference block and other blocks in After calculating this distance, blocks within a particular range are grouped into similar blocks as follows: After obtaining S i R ,j R , blocks of the noisy image y (i, j) and basic estimationx shr (i, j) that are located at these coordinates are stacked and grouped into Y S i R ,j R andX shr S i R ,j R , respectively.
In the second step, Y S i R ,j R andX shr S i R ,j R are transformed via 3D transformation T 3D . Collaborative filtering W S i R ,j R is then performed, and an inverse 3D transformation T −1 3D is performed to obtainX Wie , where noise is removed. This can be expressed aŝ The collaborative filter operates as an element-wise multiplication of the filter and a 3D transform. Collaborative filtering uses a Wiener filter, which is also based on LLMMSE. The LLMMSE filter (9) is an adaptive Wiener shrinkage filter in the transform domain [45]. The empirical Wiener shrinkage defined in BM3D is where σ 2 η shr i R ,j R represents the noise variance of the block with the top-left coordinates (i R , j R ). The original BM3D experiences Gaussian noise; therefore, σ 2 However, in the case of Poisson-Gaussian noise, the value of the signal changes according to the pixel position; thus, the power of the noise also changes. Therefore, the noise variance values for each group must be set. Noise is estimated using the difference between the noisy image and the basic estimation as follows: Because E η shr = 0, the noise variance of the reference block is E η shr 2 , which is expressed as After removing the noise of the group, the final denoised result is obtained by decomposing the group and distributing NSDFB decomposition: L m (y) → L m,1 , L m,2 , · · · , L m,N m (y).

5:
for n = 1 to N m do 6: Obtain noise removed subband L m,n x shr using (19). 7: end for 8: NSDFB reconstruction: L m x shr ← L m,1 , L m,2 , · · · , L m,N m x shr . 9: NSP reconstruction:  (23), (24), and (26). 14: Obtainx using (27) and (28). the blocks. A pixel can belong to more than one group; therefore, it can be estimated multiple times and have different values each time. These values should be averaged using the appropriate weights, Here,x Wie U (i, j) is the estimated value obtained through decomposition in group U , w U is the corresponding weight, U is all groups including a pixel at positions (i, j), and V = U ∈U (i,j) w U is a normalization factor. As shown in [11], w U is set in inverse proportion to the noise and Wiener shrinkage coefficient as follows: The procedure of the proposed Poisson-Gaussian noise removal method based on the LLMMSE NSCT shrinkage is shown in Algorithm 1, and the result of each algorithm step is shown in Fig. 4. After decomposing the noisy image using NSCT, it is assumed that there is no noise in the low-band of the coarsest level, and noise is removed from the subbands of the coarsest level. The low-band of the finer level is reconstructed from the subbands from which the noise has been removed, as shown in Fig. 4(b) and 4(c). Fig. 4(d) shows the basic estimation obtained by an iterative process up to the finest level. The final result is obtained by performing BM3D which is modified to fit the Poisson-Gaussian noise using basic estimation, as shown in Fig. 4(e).

IV. EXPERIMENTAL RESULTS AND DISCUSSION
In this section, various experiments are conducted to prove the efficiency and robustness of the proposed Poisson-Gaussian noise algorithm. It is difficult to obtain a ground truth image for an X-ray image; therefore, experiments were conducted by artificially adding noise to noiseless images to quantitatively evaluate the noise removal effect. In addition, the proposed algorithm was applied to an actual X-ray image to confirm the noise-removal performance. The parameters of the proposed method used in the experiment are listed in Table. 1.
As a measure of quality, we used the peak SNR (PSNR), which is defined as PSNR = 10 log 10 where I max is the maximum value admitted by the data format and the mean-square error (MSE) is defined as where Row and Col are the sizes of the rows and columns of the image, respectively. In addition, the structural similarity index (SSIM) was used, which is defined as where µ x and µ y are the averages of x and y, respectively, σ 2 x and σ 2 y are the variances of x and y, respectively, σ xy is the covariance of x and y, c 1 = (0.01L) 2 , and c 2 = (0.03L) 2 , where L is the dynamic range of the pixel values.

A. DATA SET
The original images used for the synthetic images are shown in Fig. 5. Experiments using synthetic images were conducted using Barbara image, shown in Fig. 5(a), with many straight lines and patterns, Lena image, shown in Fig. 5(b), with flat areas and curves. As shown in Fig. 5(c)-(f), relatively clean real X-ray images of patients [46] were also synthesized. Images have a range of 0 to 255, and noisy images were synthesized by changing the α and β values from low to high.
α was set to values of 0.1, 0.5, and 1, and β was set to 1, 5, and 10. The noise conditions were classified as Poissondominant, Gaussian-dominant, moderately dense, and highly dense noise, and the combinations of the noise parameter values (α, β) are as follows: (1, 1), (0.1, 10), (0.5, 5), and (1, 10). The actual X-ray images to be used in the experiments are shown in Fig. 6. The images with improved contrast were shown together with the original images because it is difficult to check the details in the images that had not been post-processed. The actual X-ray images were acquired using a clinical angiography prototype system supported by Samsung Electronics and a chest phantom (Multipurpose Chest Phantom N1 ''LUNGMAN'', Kyoto Kagaku, Kyoto, Japan), that yields life-like radiographs that are very close to actual clinical images. The image intensity had a 12-bit range. To obtain an image that included Poisson-Gaussian noise, a low-dose X-ray was used. The detailed image acquisition environment was as follows: the source-to-image-receptor distance was 120 cm, and the source-to-object distance was 70 cm. The radiation exposure level was 1.94 µGy/pulse and the scan parameters were 62 kVp and 40 mA for Real1. And the radiation exposure level was 1.49 µGy/pulse and the scan parameters were 63 kVp and 20 mA for Real2.

B. COMPARING METHODS
The performance of the proposed algorithm was compared using GAT+BM3D [26], the Poisson-Gaussian unbiased risk estimator linear expansion of thresholds (PURE-LET) [32], and PG-URE [33], which are state-of-the-art algorithms that remove Poisson-Gaussian noise. The GAT is the most reliable VST for Poisson-Gaussian noise, and BM3D is a representative Gaussian noise removal algorithm with the best performance; therefore, it was selected as a comparison method. The performance of BM3D was maximized through two steps of basic estimation and noise removal, and the loss of high-frequency components was minimized using patch-based 3D block filtering. In addition, PURE-LET and PG-URE were used for the comparison of results because they performed in the multiscale transform domain. Moreover, they did not use VST, but the characteristics of noise itself, similar to the proposed method. In PURE-LET, the Poisson-Gaussian model was built by changing the maximum value of the pixel intensity without setting the Poisson noise parameter. That is, the noise parameter α was fixed to 1, and the intensity of Poisson noise was determined by adjusting the maximum value of the pixel value. It was performed on a pixel-wise; therefore, similarity with neighboring pixels is not considered, unlike BM3D. PG-URE is an extended and developed algorithm based on PURE-LET. Unlike PURE-LET, which was performed only in the multiscale transform domain, it had been extended based on the wavelet transform, total variation, and nonlocal means, and VST can be applied to each method. A model that included a gain value that normalized after adjusting the maximum value of the signal was used. In this experiment, a method that performed in the wavelet transform domain without VST was selected because this method is similar to the proposed method.

C. QUANTITATIVE EVALUATION ON SYNTHETIC IMAGES
After synthesizing various noise environments to the images in Fig. 5, each method was applied and the results were compared. Table. 2 Table. 2 shows the PSNR results for each experimental image by noise condition and noise removal method. For moderately dense noise condition, GAT+BM3D exhibits best performance, and for Gaussian-dominant noise condition, the proposed method exhibits good performance in the all test images. In Chest 2-4 images, PURE-LET shows the best performance in Poisson-dominant and highly dense noise condition. Because those two conditions had the noise parameter α of 1; therefore, they were suitable for the noise model of PURE-LET. PG-URE did not record the highest value in any image, but exhibits a sense of stability showing consistent PSNR values in any noise situation. The proposed method exhibits the best performance on average for all images except Barbara image, for which the GAT+BM3D performance is the best on average. Because Barbara image has more straight lines and patterns than curves, it seems that patch-based method of BM3D is more advantageous than the shrinkage in the NSCT domain of the proposed method. Additionally, even if the proposed method did not record the best value when compared with other methods, it recorded the second highest value. In terms of quantitative evaluation, GAT+BM3D and the proposed method had a stable performance in all noise situations. PURE-LET exhibited a difference in performance for each noise situation. Particularly, when the Poisson noise is small, the performance of this method was significantly degraded. In PURE-LET, when the Poisson noise intensity in the image was small, the maximum value of the signal had to be normalized to a very small value. As the maximum value of the signal decreased, the difference between the signals became subtle, and the intensity of the estimated noise was too small, or the Gaussian noise was misrecognized as Poisson noise. For this reason, PURE-LET exhibits poor results for images with low Poisson noise intensity, but excellent performance for images with high Poisson noise intensity. Owing to the extended method and model, PG-URE exhibits a similar performance in various noisy environments. The variation in performance for each situation was small, and the method showed satisfactory performance, even in Gaussian-dominant noise, which was a limitation of PURE-LET. However, owing to the complicated model, it was difficult to obtain both noise reduction and high-frequency component preservation because the method is sensitive to the setting of parameter values. Table. 3 shows the SSIM results for synthetic images. In most cases, methods with high values in PSNR also exhibit high values in SSIM. For Poisson-dominant noise condition, GAT+BM3D or PURE-LET exhibit the best SSIM values, and for Gaussian-dominant noise condition, the proposed method exhibits the best SSIM scores. For the results of the moderately dense noise, GAT+BM3D exhibits the best values and for the results of the highly dense noise, the proposed method shows the highest SSIM values except for Lena image. However, despite the highest PSNR value of PURE-LET, there are cases where the SSIM value of GAT+BM3D or the proposed method is the highest. In the Poisson-dominant noise condition of the Chest2 image, where the PSNR was best in PURE-LET, the SSIM was the best in GAT+BM3D, and in the highly dense noise condition of the chest3 and chest4 images, the SSIM of the proposed method was the highest. GAT+BM3D, PG-URE, and the proposed method show uniform improvement in SSIM values for all noise situations, and PURE-LET has poor SSIM results in Gaussian-dominated noise condition. As for the average SSIM value, the proposed method was the best for all of the images except for the Barbara image where GAT+BM3D is the best. Fig. 7-10 show the results of noise removal from an image degraded by Poisson-dominant noise, Gaussian-dominant noise, moderately dense noise, and highly dense noise, respectively. The red boxes are enlarged and displayed in the insets of the images. In Fig. 7, the results of GAT+BM3D and the proposed method show that the noise is well suppressed, and the edge and detail, which are high-frequency components, are also well preserved. In the case of PURE-LET, the noise was removed; however, the noise removal performance was inferior to that of the other methods. The high-frequency components were not damaged and were very well preserved, thereby resulting in sharp results. PG-URE removed most of the noise,  but it also damaged the high-frequency components and exhibited a blurry result. Fig. 8 shows the noise removal results of Chest1 images dominated by Gaussian noise. The GAT+BM3D method and the proposed method show stable results as the case of Poisson-dominant noise. However, PURE-LET hardly suppressed Gaussian-dominant noise. As mentioned in Section IV-C, PURE-LET did not have the gain of Poisson noise in the noise model; thus, when the Poisson noise parameter was small, the noise removal performance of PURE-LET decreased. PG-URE suppressed the noise excessively and damaged many details in the process, resulting in a blurry image, as in the case of Poisson-dominant noise removal. Fig. 9 shows the noise removal results of Lena images in which the Poisson and Gaussian noise were mixed by an appropriate amount. All methods effectively removed the noise, and GAT+BM3D and the proposed method also yielded stable results that preserved high-frequency components with the noise removal. The results of PURE-LET showed that high-frequency components were maintained, and they looked sharp. However, thin edges were damaged and disappeared, but the only strong edges were preserved. In the PG-URE results, all the weak edges disappeared, and the strong edges were also damaged, thereby indicating that the image was blurred overall. Fig. 10 shows the results of Chest2-4 image sets for highly dense noise. GAT+BM3D, PG-URE, and the proposed method show a similar tendency to the results of other noise environments. PURE-LET performed well for suppressing the noise, but did not effectively recover the high-frequency components from the Chest image. As shown in Table. 2 and 3, PURE-LET removed the noise well and shows high PSNR, which indicates the fidelity of data, but SSIM is lower than the proposed method because the structure of the image is not well reconstructed. Because only the differential in the vertical and horizontal direction was used in PURE-LET, the accuracy in the curve was relatively low, and the structural part of the result is restored in poorer performance than the patch-based method because PURE-LET was based on the pixel-wise method. The proposed method considered the curve structure of image and applied the patch-base method to restore high-frequency components. Although the PSNR of the proposed method shows lower than that of PURE-LET, but the SSIM of the proposed method exhibits higher than that of PURE-LET.

E. EVALUATION ON REAL X-RAY IMAGE
The algorithms were applied to actual X-ray images. Noise parameters were estimated from the acquired low-dose X-ray images, and their values were estimated to be α = 0.4663 and β = 7.2111 for Real1, and α = 0.0337 and β = 10.9457 for Real2. The noise environments of Real1 and  Real2 images were similar to that of moderately dense noise and Gaussian-dominant noise in the synthetic images, respectively. The noise removal results of real X-ray raw images are shown in Fig. 11. To confirm the damage and preservation of the high-frequency components, the contrasts were enhanced for the denoised images. The performance of GAT+BM3D was still good, and PURE-LET exhibited little noise, which was properly removed, and the high-frequency components were less damaged, thereby making it easier to identify the edge compared with other methods. PG-URE did not suppress much noise for Real1 images, however, effectively removes noise while preserving high-frequency components for Real2 image. Finally, the proposed method removed the noise cleanly, and it effectively restored the original signal that was corrupted by noise.

V. CONCLUSION
We proposed the LLMMSE shrinkage method in the NSCT based on noise characteristics using an experimental analysis and a theoretical approach. The proposed method has a similar structure to BM3D, and noise is removed through basic estimation using shrinkage and a Poisson-Gaussian-oriented Wiener filter. The NSCT has been used for medical images that consist of many curves, and the LLMMSE filter was designed and derived in the form of shrinkage, based on the noise analysis in the NSCT. The reliability of the algorithm was improved through a two-step process, and high-frequency components that may be damaged during noise removal were preserved as much as possible via patch-based block filtering.
The results for the artificially synthesized images were very satisfactory, and the proposed method exhibited the best or second highest PSNR and SSIM values. It showed consistent performance in all noise environments and recorded the most evaluated value on average when compared with other noise removal methods. Experiments on real X-ray images are also encouraging, because proposed method can better preserve relevant details while suppressing noise and smoothing out homogeneous areas.
The NSCT has the advantage of separating the image into the scale and direction components, but has the disadvantage that the size of the decomposition and reconstruction filter increases as the scale level increases, owing to the elimination of the subsampling process. For a stable noise removal performance, the maximum scale level should be set to M ≥ 3. If the size of the image is small, there is a possibility that the NSP may not work well, owing to the large size of the decomposition filter. Future work will include the use of transforms other that the NSCT or those with subsampling to remove noise in small images.