# Noise Estimation for Single-Slice Sinogram of Low-Dose X-Ray Computed Tomography Using Homogenous Patch

Authors: Zhiwu Liao; Shaoxiang Hu; Ming Li; Wufan Chen
Journal: Mathematical Problems in Engineering
 (2012)

---

## Abstract

We present a new method to estimate noise for a single-slice sinogram of low-dose CT based on the homogenous patches centered at a special pixel, called center point, which has the smallest variance among all sinogram pixels. The homogenous patch, composed by homogenous points, is formed by the points similar to the center point using similarity sorting, similarity decreasing searching, and variance analysis in a very large neighborhood (VLN) to avoid manual selection of parameter for similarity measures.Homogenous pixels in the VLN allow us find the largest number of samples, who have the highest similarities to the center point, for noise estimation, and the noise level can be estimated according to unbiased estimation.Experimental results show that for the simulated noisy sinograms, the method proposed in this paper can obtain satisfied noise estimation results, especially for  sinograms with relatively serious noises.


---

## Body

## 1. Introduction

With continued technology advancement and wider applications,  use of computed Tomography (CT) is increasing.  However  radiation exposure and associated risk of cancer for patients receiving CT examination have been an increasing concern in recent years. Thus, minimizing X-ray exposure to patients has been one of the major efforts in the CT field [1–3].A simple and cost-effective means to achieve low-dose CT applications is to lower X-ray tube current (mA) as low as achievable. However, dose reduction generally leads to an increased level of noise in the measured projection data (sinogram) and the subsequent reconstructed images.Sinograms acquired from low-dose CT are corrupted by many factors, including Poisson noise, logarithmic transformation of scaled measurements, and prereconstruction corrections for system calibration [4]. All these factors complicate noise modeling on the sinogram. Up to now, various noise models for sinogram have been developed [4–15].One common model for noise estimation is Gaussian distribution with variance depending on the sinogram [4–6]. The model is developed by repeatedly acquired projection measurements of a physical phantom at a fixed angle for 900 times by a GE spiral CT scanner [5, 6].Another model is Poisson distribution based on quantum noise for CT which is due to the limited number of photons collected by the detector [7, 8]. Although a CT detector is not based on photon-counting but energy integrated that generates a signal proportional to the total energy deposited in the detector,  a photon-counting model is still a good approximation and is widely used for characterizing noise properties of the sinogram.More accurate noise model, compound Poisson model (CPM), which takes into account both the polychromatic X-ray beam and energy integration, has been investigated in [9–11]. However, for its complex expression, the exact form of the likelihood is not amendable. Therefore, many researchers have often approximated the CPM by either the Poisson distribution or the Gaussian distribution [12, 13].Based on the above models, dose-reduction simulation using synthetic-noise generators, enables ethical studies of low-dose procedures to improve the quality of reconstruction images. Dose-reduction simulation has been reported using Poisson or Gaussian noise models [14, 15]. However, above simulation methods assume that the parameters either for Poisson or for Gaussian are known. It is unreal in clinical X-ray CT systems. Thus in order to get satisfactory reconstruction images, noise estimation becomes a key problem in low-dose CT imaging both for sinogram and for reconstructed images.In this paper, we focus on how to estimate parameters for added signal independent Gaussian noise (SIGN) on sinogram. Although it does not coincide noise models introduced in this section, it is still a valuable model in investigating the properties of the noise for sinogram and the relations between common noise models for sinogram which will be discussed in Section2.Noise estimation for a single image is an important and difficult problem in image denoising for complex structures of images and is studied in the early 90s for the last century [16–19]. The main start point for these methods is that the noise level should be estimated using the smoother versions of images. However, it leads to their high computation burden and overestimate of noise levels.Recently, noise estimation from a single image both for Gaussian and Poisson noise becomes a hot spot in computer vision [20–24]. However, most of these methods must pose complex prior such as, camera parameters and so forth [21, 22], or have many parameters chosen by hand [20, 24] or segment images previously [20, 23] which hampers their application in noise estimation for sinogram.In this paper, we insist that noise estimation should be performed in a local homogenous patch and present a method, which does not need presegment or pose complex prior. Moreover, in order to get more reliable estimated results and to improve the overestimate when noise levels are low, motivated by nonlocal means and some most recent results [25–41], a very large neighborhood is adopt to find similar points. Nonlocal means estimates each real gray level of an image based on block similarity in a nonlocal neighborhood with size 21×21 and very recently is used in low-dose CT imaging [25–32].In the very large neighborhood, the homogenous patches are determined by finding similar points to the center point through similarity sorting, decreasing similarity searching and variance analysis. The noise parameters will be estimated on this homogenous patch using unbiased estimate. Since very large neighborhood provides more reliable estimation, proposed method can get satisfied noise estimation results.The remainder of this paper is arranged as the noise models will be discussed in Section2; the motivation and method will be given in Sections 3 and 4, respectively; and then the experimental results will be given in Section 5; finally, it is the conclusion, future works, and acknowledgment.


## 2. Noise Models

For clinical X-ray system, its detected that X-ray intensity follows a compound Poisson distribution [9–11]. However, for its complex expression, it has no analytic formula for the likelihood function. Various approximations have been proposed for it. Two common models include Poisson distribution and Gaussian distribution.In this section, we will introduce signal-independent Gaussian noise (SIGN), Poisson noise, and signal-dependent Gaussian noise, as well as their relations and the reason for addressing SIGN.


### 2.1. Signal-Independent Gaussian Noise (SIGN)

SIGN is a common noise for imaging system. Poisson noise and signal-dependent Gaussian noise can be converted to SIGN using scale transforms which will be discussed in Section2.3.Let the original projection data be{xi},i=1,…,m, where i is the index of the ith bin. The signal has been corrupted by additive noise {ni},i=1,…,m, and one noisy observation
(2.1)yi=xi+ni,
where ni is an observation for the random variable Ni as normal N(0,σN2) and independent to the Gaussian random variable Xi where the uppercase letters denote the random variables and the lower-case letters denote the observations for respective variables.


### 2.2. Poisson Model and Signal-Dependent Gaussian Model

The photon noise is due to the limited number of photons collected by the detector [42]. For a given attenuating path in the imaged subject, N0(i,α) and N(i,α) denote the incident and the penetrated photon numbers, respectively. Here, i denote the index of detector channel or bin and α is the index of projection angle. In the presence of noises, the sinogram should be considered as a random process and the attenuating path is given by(2.2)ri=-ln⁡⁡[N(i,α)N0(i,α)],
where N0(i,α) is a constant and N(i,α) is Poisson distribution with mean N.Thus we have(2.3)N(i,α)=N0(i,α)exp⁡(-ri).Both its mean value and variance areN.Gaussian distributions of ployenergetic systems were assumed based on limited theorem for high-flux levels and followed many repeated experiments in [6]; we have(2.4)σi2(μi)=fiexp⁡(μiγ),
where μi is the mean and σi2 is the variance of the projection data at detector channel or bin i, γ is a scaling parameter and fi is a parameter adaptive to different detector bins.The most common conclusion for the relation between Poisson distribution and Gaussian distribution is that the photon count will obey Gaussian distribution for the case with large incident intensity and Poisson distribution with feeble intensity [6]. In addition, in [42], the authors deduce the equivalency between Poisson model and Gaussian model. Therefore, both theories indicate that these two noises have similar statistical properties and can be unified into a whole framework.


### 2.3. Scale Transformations

As a statistical method for treating nonnormally distributed or signal dependent noise, scale transformations are widely-used to stabilize the variance [6, 31]. In [4], authors indicate that the variance of their signal-dependent Gaussian model represented in (2.4) can be stabilized by(2.5)g(x)=ckln⁡(kx+b)+c1,
where c is the expected constant standard deviation for transformed data and c1 is an arbitrary constant or
(2.6)g(x)=2c4ac-b2arctan2ax+b4ac-b2+c2,
where c2 also is an arbitrary constant.Armando Manduca et al. also indicate that the Anscombe transform can convert Poisson distribution data to data with an approximately normal distribution with a constant variance. Thus we can useQ=2N+3/8 to convert data with distribution expressed in (2.3) to a normal distribution with a constant variance.Since both Poisson noise and signal-dependent Gaussian noise can be converted to SIGN, the noise estimation for low-dose CT can start from SIGN to focus on our new method itself.


## 2.1. Signal-Independent Gaussian Noise (SIGN)

SIGN is a common noise for imaging system. Poisson noise and signal-dependent Gaussian noise can be converted to SIGN using scale transforms which will be discussed in Section2.3.Let the original projection data be{xi},i=1,…,m, where i is the index of the ith bin. The signal has been corrupted by additive noise {ni},i=1,…,m, and one noisy observation
(2.1)yi=xi+ni,
where ni is an observation for the random variable Ni as normal N(0,σN2) and independent to the Gaussian random variable Xi where the uppercase letters denote the random variables and the lower-case letters denote the observations for respective variables.


## 2.2. Poisson Model and Signal-Dependent Gaussian Model

The photon noise is due to the limited number of photons collected by the detector [42]. For a given attenuating path in the imaged subject, N0(i,α) and N(i,α) denote the incident and the penetrated photon numbers, respectively. Here, i denote the index of detector channel or bin and α is the index of projection angle. In the presence of noises, the sinogram should be considered as a random process and the attenuating path is given by(2.2)ri=-ln⁡⁡[N(i,α)N0(i,α)],
where N0(i,α) is a constant and N(i,α) is Poisson distribution with mean N.Thus we have(2.3)N(i,α)=N0(i,α)exp⁡(-ri).Both its mean value and variance areN.Gaussian distributions of ployenergetic systems were assumed based on limited theorem for high-flux levels and followed many repeated experiments in [6]; we have(2.4)σi2(μi)=fiexp⁡(μiγ),
where μi is the mean and σi2 is the variance of the projection data at detector channel or bin i, γ is a scaling parameter and fi is a parameter adaptive to different detector bins.The most common conclusion for the relation between Poisson distribution and Gaussian distribution is that the photon count will obey Gaussian distribution for the case with large incident intensity and Poisson distribution with feeble intensity [6]. In addition, in [42], the authors deduce the equivalency between Poisson model and Gaussian model. Therefore, both theories indicate that these two noises have similar statistical properties and can be unified into a whole framework.


## 2.3. Scale Transformations

As a statistical method for treating nonnormally distributed or signal dependent noise, scale transformations are widely-used to stabilize the variance [6, 31]. In [4], authors indicate that the variance of their signal-dependent Gaussian model represented in (2.4) can be stabilized by(2.5)g(x)=ckln⁡(kx+b)+c1,
where c is the expected constant standard deviation for transformed data and c1 is an arbitrary constant or
(2.6)g(x)=2c4ac-b2arctan2ax+b4ac-b2+c2,
where c2 also is an arbitrary constant.Armando Manduca et al. also indicate that the Anscombe transform can convert Poisson distribution data to data with an approximately normal distribution with a constant variance. Thus we can useQ=2N+3/8 to convert data with distribution expressed in (2.3) to a normal distribution with a constant variance.Since both Poisson noise and signal-dependent Gaussian noise can be converted to SIGN, the noise estimation for low-dose CT can start from SIGN to focus on our new method itself.


## 3. Backgrounds and Motivation

Following the SIGN model discussed in Section2.1, we will introduce backgrounds and motivation for proposed method.Noise estimation for low-dose CT projection data is an estimation problem. Using the terms in Section2.1, the goal for SIGN estimation is to get unbiased estimated value σ̂N of σN from the noisy observation yi. Here σN is the parameter for population distribution of the noise.However, only one noisy observationyi is provided. Thus we have to share statistical information each other. One reasonable way for sharing statistical information is to assume independent identical distributions (iid) for nearby similar pixels among noisy projection data (sinogram). In this paper, the similar point is the point with high block similarity to the center point. This section will discuss the motivation for finding samples of unbiased estimator.


### 3.1. Unbiased Estimation

Suppose we have a statistical model parameterized byθ, and a statistic θ̂ which serves as an estimator of θ based on any samples x. That is, we assume that population data follows Gaussian distribution with a fixed and unknown constant θ, and then we construct some estimator θ̂ that maps samples to values that we hope are close to θ. Then the bias of this estimator is defined to be(3.1)Bias(θ̂)=E(θ̂)-E(θ)=E(θ̂-θ),
where E(·) denotes expected value over the distribution, that is, averaging overall possible samples x.An estimator is said to be unbiased if its bias is equal to zero for parameterθ. Suppose x={x1,…,xn} are observations for a Gaussian random variable with expectation μ and variance σ, the unbiased estimators for these samples are(3.2)X¯=1n∑i=1nxi,S2=1n-1∑i=1n(xi-X¯)2.In theory, the sample estimators shown in (3.2) are asymptotically unbiased and efficient for the sample size is moderate or large. Thus unbiased estimation requires providing enough samples to approximate the population distribution.


### 3.2. Measure Similarity

In order to measure similarity between two different pointsxi,j and xs,t of noisy sinogram where i≠s or j≠t, we should measure two 7×7 blocks centered at xi,j and xs,t, respectively. The similarity is measured as
(3.3)S[(i,j),(s,t)]=∑k=-3k=3∑m=-3m=3[xi+k,j+m-xs+k,t+m].Generally, we can predefine a thresholdT to find the similar points of xi,j in a nonlocal neighborhood with 21×21 centered at xi,j(3.4)xs,tisasimilarpointofxi,jifS{(i,j),(s,t)}≤T,xs,tisnotasimilarpointofxi,jifS{(i,j),(s,t)}>T,
where s≠i or t≠j and s,t=-10,…,10.


### 3.3. Motivation

Since only samples with moderate or large size can approximate the population distribution well, we must find enough samples to estimate the noise correctly. Moreover, iid assumption requires to share statistical information among similar points of noisy sinogram. Both requirements coincide with the start point of existing methods. That is, noise estimation should be performed at smoother versions of the images to suppress the influence of the complex structures of images.However, estimation noise using smoother versions of the images has high computation burden and overestimate noise levels [24]. In this paper, we propose a method to estimate noise in a VLN by similarity sorting and variance analysis.Just as discussed in this section, the key objective is to find enough similar points for noise estimation. Unlike existing global methods, which estimate noise levels using whole images, the proposed method try to estimate noise levels in a more local way to reduce computation burden and increase flexibility for estimation.In order to accomplish the above objectives (enough samples with a more local structure), VLN is used to find similar points (samples). Moreover, the VLN should be put in a suitable position to ensure the enough samples can be found. That is, VLN should be put in a homogenous region. Since the center of VLN can determine the position of VLN, how to put VLN in a suitable position can be converted to how to locate the center point in a homogenous region.The noise level of low-dose CT is low. In this situation, only variance is enough for describing the local homogeneity roughly. That is, large variance relates to a square near singularities while small variance relates to a homogenous square. Therefore, by comparing variances of squares with fixed size for all pixels in sinogram, the center point is defined as the pixel with the smallest variance among all pixels.After determining the center point, we propose a new method for noise estimation based on similarity sorting, similarity decreasing searching, and variance analysis. Its main start point is from how to avoid threshold setting in finding similar points to the center point since the threshold selection depends on the noise level which makes it a “chicken and egg” problem.Similarity sorting provides a similarity decreasing sequence (SDS). Thus the samples must be formed by the up at the front points in the SDS since points with smaller similarities to the center point maybe the outliers for the estimation.In order to find the largest number samples with the highest similarities in the VLN, similarity decreasing searching combining variance analysis is used. That is, the samples are formed by adding 50 points each time according to the order of the SDS. Thus it ensures that each addition adds the points with the highest similarities in residual SDS.Moreover, in order to find the largest number samples used in estimation, we must find when the outliers are added. It can be achieved by variance analysis for each addition. When the variance for an addition becomes large suddenly, it means some outliers are added. Therefore, the real samples are the last before this addition. The detailed algorithm will be introduced in Section4 and two-number examples can be found in Section 5.3.


## 3.1. Unbiased Estimation

Suppose we have a statistical model parameterized byθ, and a statistic θ̂ which serves as an estimator of θ based on any samples x. That is, we assume that population data follows Gaussian distribution with a fixed and unknown constant θ, and then we construct some estimator θ̂ that maps samples to values that we hope are close to θ. Then the bias of this estimator is defined to be(3.1)Bias(θ̂)=E(θ̂)-E(θ)=E(θ̂-θ),
where E(·) denotes expected value over the distribution, that is, averaging overall possible samples x.An estimator is said to be unbiased if its bias is equal to zero for parameterθ. Suppose x={x1,…,xn} are observations for a Gaussian random variable with expectation μ and variance σ, the unbiased estimators for these samples are(3.2)X¯=1n∑i=1nxi,S2=1n-1∑i=1n(xi-X¯)2.In theory, the sample estimators shown in (3.2) are asymptotically unbiased and efficient for the sample size is moderate or large. Thus unbiased estimation requires providing enough samples to approximate the population distribution.


## 3.2. Measure Similarity

In order to measure similarity between two different pointsxi,j and xs,t of noisy sinogram where i≠s or j≠t, we should measure two 7×7 blocks centered at xi,j and xs,t, respectively. The similarity is measured as
(3.3)S[(i,j),(s,t)]=∑k=-3k=3∑m=-3m=3[xi+k,j+m-xs+k,t+m].Generally, we can predefine a thresholdT to find the similar points of xi,j in a nonlocal neighborhood with 21×21 centered at xi,j(3.4)xs,tisasimilarpointofxi,jifS{(i,j),(s,t)}≤T,xs,tisnotasimilarpointofxi,jifS{(i,j),(s,t)}>T,
where s≠i or t≠j and s,t=-10,…,10.


## 3.3. Motivation

Since only samples with moderate or large size can approximate the population distribution well, we must find enough samples to estimate the noise correctly. Moreover, iid assumption requires to share statistical information among similar points of noisy sinogram. Both requirements coincide with the start point of existing methods. That is, noise estimation should be performed at smoother versions of the images to suppress the influence of the complex structures of images.However, estimation noise using smoother versions of the images has high computation burden and overestimate noise levels [24]. In this paper, we propose a method to estimate noise in a VLN by similarity sorting and variance analysis.Just as discussed in this section, the key objective is to find enough similar points for noise estimation. Unlike existing global methods, which estimate noise levels using whole images, the proposed method try to estimate noise levels in a more local way to reduce computation burden and increase flexibility for estimation.In order to accomplish the above objectives (enough samples with a more local structure), VLN is used to find similar points (samples). Moreover, the VLN should be put in a suitable position to ensure the enough samples can be found. That is, VLN should be put in a homogenous region. Since the center of VLN can determine the position of VLN, how to put VLN in a suitable position can be converted to how to locate the center point in a homogenous region.The noise level of low-dose CT is low. In this situation, only variance is enough for describing the local homogeneity roughly. That is, large variance relates to a square near singularities while small variance relates to a homogenous square. Therefore, by comparing variances of squares with fixed size for all pixels in sinogram, the center point is defined as the pixel with the smallest variance among all pixels.After determining the center point, we propose a new method for noise estimation based on similarity sorting, similarity decreasing searching, and variance analysis. Its main start point is from how to avoid threshold setting in finding similar points to the center point since the threshold selection depends on the noise level which makes it a “chicken and egg” problem.Similarity sorting provides a similarity decreasing sequence (SDS). Thus the samples must be formed by the up at the front points in the SDS since points with smaller similarities to the center point maybe the outliers for the estimation.In order to find the largest number samples with the highest similarities in the VLN, similarity decreasing searching combining variance analysis is used. That is, the samples are formed by adding 50 points each time according to the order of the SDS. Thus it ensures that each addition adds the points with the highest similarities in residual SDS.Moreover, in order to find the largest number samples used in estimation, we must find when the outliers are added. It can be achieved by variance analysis for each addition. When the variance for an addition becomes large suddenly, it means some outliers are added. Therefore, the real samples are the last before this addition. The detailed algorithm will be introduced in Section4 and two-number examples can be found in Section 5.3.


## 4. The Method

In this section, we will give the framework for noise estimation. Just as shown in Figure1, it includes three steps.(1)
Find the center point:  It is the key step to locate the VLN, which should be put in a homogenous region. Thus the center point is defined as the point with the smallest local variance in all sinogram pixels.(2)
Determine the homogenous patch: The homogenous patch is composed by similar points to the center point (samples), and these similar points are searched using similarity sorting, similarity decreasing searching and variance analysis.(3)
Estimate noise: The parameters are estimated using unbiased estimator by the samples on the homogenous patch.Figure 1
The flow chart for the proposed method.By these steps, the largest number of samples with the highest similarities can be obtained. These samples form a very large homogenous patch for noise estimation. Thus even using simple unbiased estimate, satisfied results can be obtained. The details for these three steps will be given in the remainder of this section.


### 4.1. Find Center Pixel

The center pixel should locate at the center of a homogenous region. However, in noise, finding a largeregular homogenous region correctly and putting a pixel of sinogram in the center of this region is not a trivial task.Motivated by [27], irregular pixel patch composed by similar pixels to the center point is a reasonable choice. Thus the choice for the center point becomes an easy task which only needs to ensure it is the center of a small homogenous square.Since the variance can describe how far the numbers lie from the mean, small variances will relate to homogenous squares while large variances will relate to squares near singularities. Thus the center point is chosen as the pixel with the smallest variance in a7×7 square centered at the point among all pixels of the sinogram. That is, the position of the center point should be
(4.1)(i,j)=argxi,j∈X(i,j)min⁡∑m=-3m=3∑n=-3n=3(xi-m,j-n-x¯i,j)2,
where X is the sinogram and
(4.2)x¯i,j=149∑s=-3s=3∑t=-3t=3xi-s,j-t
is the mean.By this way, we can find a point centered at a homogenous square and then extend it to an irregular homogenous patch.


### 4.2. Determine Similar Points (Homogenous Patch)

After finding the center pointxi,j, the irregular homogenous patch P composed by the similar points to xi,j should be determined.Each similar point is found by computing the similarity between itself andxi,j using (3.3). However, unlike existing methods which use a predefined threshold value T shown in Section 4.1, our method proposes a different scheme by similarities sorting and variance analysis in a VLN.The main advantage for this scheme is that it avoids threshold setting in finding similar points. Since the threshold should be set according to different noise level, noise estimation and threshold setting become a “chickens and eggs” question.The most straight motivation for the proposed method is that the small variance indicates only the homogenous samples, while the large variance indicates not only the homogenous samples but also outliers. Thus the variances can be considered as a sign for outliers. It reminds us that if we sort the points according to the their similarities toxi,j in a VLN and then compute the variances for increasing adding points until the big variance is reached, the most reliable noise estimation can be obtained by the as large as possible number for the samples.Thus we can add samples according to the fact that samples with bigger similarities will be added earlier,  and compute the variance after each addition. When the added variance become big suddenly, it means that there are some outliers that are mixed in estimation. Thus the real noise level is estimated using the last before this addition.In summary, the steps for finding similar points are as follows.(1)
Compute the similarity between the center pointXi,j and each point in the 41×41 square using (3.3).(2)
Sort the similarities to form a similarity descending sequence (SDS).(3)
Add 50 samples each time to form collections of samples and then estimate the standard deviations (SDs) of the collection of samples to form an SDs sequence.(4)
Compute the variance difference using thei+1th and ith, i=1,…,31, elements of the SDs sequence to form a difference variances sequence and from the 21 addition, find the variance becomes big suddenly (is bigger than 5), and denote its index k.(5)
The estimated variance of the noise isk-1th element of the SDs sequence obtained on the step 3.


## 4.1. Find Center Pixel

The center pixel should locate at the center of a homogenous region. However, in noise, finding a largeregular homogenous region correctly and putting a pixel of sinogram in the center of this region is not a trivial task.Motivated by [27], irregular pixel patch composed by similar pixels to the center point is a reasonable choice. Thus the choice for the center point becomes an easy task which only needs to ensure it is the center of a small homogenous square.Since the variance can describe how far the numbers lie from the mean, small variances will relate to homogenous squares while large variances will relate to squares near singularities. Thus the center point is chosen as the pixel with the smallest variance in a7×7 square centered at the point among all pixels of the sinogram. That is, the position of the center point should be
(4.1)(i,j)=argxi,j∈X(i,j)min⁡∑m=-3m=3∑n=-3n=3(xi-m,j-n-x¯i,j)2,
where X is the sinogram and
(4.2)x¯i,j=149∑s=-3s=3∑t=-3t=3xi-s,j-t
is the mean.By this way, we can find a point centered at a homogenous square and then extend it to an irregular homogenous patch.


## 4.2. Determine Similar Points (Homogenous Patch)

After finding the center pointxi,j, the irregular homogenous patch P composed by the similar points to xi,j should be determined.Each similar point is found by computing the similarity between itself andxi,j using (3.3). However, unlike existing methods which use a predefined threshold value T shown in Section 4.1, our method proposes a different scheme by similarities sorting and variance analysis in a VLN.The main advantage for this scheme is that it avoids threshold setting in finding similar points. Since the threshold should be set according to different noise level, noise estimation and threshold setting become a “chickens and eggs” question.The most straight motivation for the proposed method is that the small variance indicates only the homogenous samples, while the large variance indicates not only the homogenous samples but also outliers. Thus the variances can be considered as a sign for outliers. It reminds us that if we sort the points according to the their similarities toxi,j in a VLN and then compute the variances for increasing adding points until the big variance is reached, the most reliable noise estimation can be obtained by the as large as possible number for the samples.Thus we can add samples according to the fact that samples with bigger similarities will be added earlier,  and compute the variance after each addition. When the added variance become big suddenly, it means that there are some outliers that are mixed in estimation. Thus the real noise level is estimated using the last before this addition.In summary, the steps for finding similar points are as follows.(1)
Compute the similarity between the center pointXi,j and each point in the 41×41 square using (3.3).(2)
Sort the similarities to form a similarity descending sequence (SDS).(3)
Add 50 samples each time to form collections of samples and then estimate the standard deviations (SDs) of the collection of samples to form an SDs sequence.(4)
Compute the variance difference using thei+1th and ith, i=1,…,31, elements of the SDs sequence to form a difference variances sequence and from the 21 addition, find the variance becomes big suddenly (is bigger than 5), and denote its index k.(5)
The estimated variance of the noise isk-1th element of the SDs sequence obtained on the step 3.


## 5. Experimental Results

In this section, we will compare our method to some well-known noise estimation methods. Just as discussed in Section1, most of existing methods for noise estimation are not suitable for sinogram for their complex prior, segmentation for images or practical setting parameters. Only methods which estimate noise by filtering images to suppress singularities previously can be used for noise estimation of sinogram [16–19, 24]. We will introduce two of these methods firstly in Section 5.2 and then compare them with the proposed method in Section 5.3.


### 5.1. Data

Five test images are used in this paper: a thorax phantom acquired from a GE HiSpeed multislice CT scanner (see Figure2(a)), two phantom data produced on Matlab (see Figures 2(b) and 2(c)), and two images acquired with a 16 multidetector row CT unit (Somatom Sensation 16; Siemens Medical Solutions) using 120 kVp, 5 mm slice thickness (see Figures 2(d) and 2(e)).Test images. (a) Thorax phantom data acquired at 120 kVp, 150 mA. (b) A phantom data produced by Matlab. (c) Modified Shepp-Logan head phantom whose size is256×256. (d) A CT image of an abdomen of a 58-year-old man with 120 kVp, 150 mAs. (e) A CT image of a chest of a 62-year-old woman with 120 kVp, 150 mAs.
(a)(b)(c)(d)(e)For Figure2(a), the distance from the center of rotation (COR) to the curved detector is 408.075 mm. The detector array is on an arc concentric to the X-ray source with a distance between the X-ray source and the COR of 541.00 mm. The detector cell spacing is 1.0239 mm, and the slice thickness is 1.0 mm. Its projection data is shown on Figure 3(a).The relative projection data for the test images.
(a)(b)(c)(d)(e)The scanning parameters for Figures2(d) and 2(e) are: gantry rotation time, 0.5 second; table feed per gantry rotation, 24 mm; and pitch, 1 : 1. The CT doses were controlled by a fixed tube currents 150 mAs. Two respective sinogram data are shown on Figures 3(d) and 3(e).In addition, two projection data for Figures2(b) and 2(c) are shown in Figures 3(b) and 3(c), respectively. The projection data are produced in Matlab using radon transform for the equally distributed 90 angles for two phantoms.The SIG noises are added on the projection data in Figure3 directly. Since the interesting SIG noises are with small variances (SIG noises with large variances are not accepted in clinic situations), the variances are form 25 to 225 (the standard deviations are from 5 to 15).


### 5.2. Comparison Methods

Noise estimation is not a trivial task in image processing because of the complex structures for images. Some researchers suggest that the noise levels should be estimated by filtering images previously to suppress the image structures [17, 18, 24]. The methods compared with the proposed method are as follows.Wavelets
This method suppresses image structures using wavelets [18]. Since coefficients in the HH subbands of wavelets only preserving high frequency energy for images, the image structures, which mainly reside in low frequency coefficients, can be considered to be suppressed in the HH subbands. Thus based on robust median estimator, the level of the noise can be estimated as
(5.1)σN=median(|WHH(X)|)0.6745,
where X is the image, W is the wavelet operator, HH indicates the coefficients in the HH subbands, |·| is the absolute value, and median is the median estimator.Fast Estimation (FE)
This method is proposed in [17]. It only has two steps for noise estimation: the first step is filtering the image using Laplacian operator and the noise level is estimated using
(5.2)σN=π216(w-2)(h-2)∑|X⊗L|,
where w and h are the width and height of the image, respectively, ⊗ is the convolution operator, L is the Laplacian operator:
(5.3)[1-21-24-21-21].


### 5.3. Experimental Results

In this section, three methods,  proposed method (PM), wavelets, and Fast estimation (FE), are compared. In order to compare these three methods on a fair stage, the parameters used in these methods are fixed for all noise levels and all test images. According to this standard, the method proposed recently in [24] is not included because one of its parameter must be adjusted and selected with the assumption that the real noise levels are known.The wavelet used for  wavelets is Symlets with support 4. In summary, the parameters for PM are the squares used for finding the center point which are7×7 squares (see Section 4.1); the size of the square computing similarity between the considering pixel and the center point is 7×7 (see Section 3.2); the size for very large neighborhood is 41×41 (see Section 4.2 step 1); adding 50 samples each time to form the SDs sequence (see Section 4.2 step 3);  the threshold of absolute difference between adjacent standard deviations for finding the variances become ssuddenly 5; in order to avoid the influence of the outliers, the comparison for the variance is from the 21th addition for samples (see Section 4.2 step 4).Firstly, in order to show the role of variance analysis for the SDs sequence, we will use two Figures, Figures2(c) and 2(d), for analyzing the SDs sequence (see Table 1), where “Real SD” is the real standard deviation (SD) for the added noise, “No” is the times adding samples, “ESD” is the estimated standard deviation for each sample adding, “Diff” is the difference between the ith and i-1th estimated variances, the bold digits indicate the difference between two adjoint variances beyond the 5 (bigger diff), and the italic digits are the final estimated SDs.Table 1
Variance analysis to estimate the noise level.
Figure2(c)Figure2(d)Real SD5101551015NO.ESDDiffESDDiffESDDiffESDDiffESDDiffESDDiff214.930.269.821.0714.670.614.950.529.65−0.2614.423.75224.940.109.860.7814.783.004.960.159.690.9014.450.95234.91-0.299.82-0.7914.74-0.984.980.199.751.2014.501.53244.88-0.299.75-1.4214.62-3.484.980.199.751.2014.520.34254.900.229.811.1914.712.635.130.809.820.5314.684.80266.6119.6810.7619.6415.3418.715.200.699.942.4214.772.42275.250.5110.021.5214.873.17285.300.5410.040.4214.890.53295.947.2410.356.1714.982.733015.3210.38According to the steps in Subsection4.2, if the index for bigger diff is i, the final estimated SD is the i-1th SD. For example, in the second column, which is Figure 2(c) added noise with real SD 5, we can find that the bigger diff is 19.68 and its index is 26 (see the last no blank row of the first and second columns). Thus the index for final estimated SD is 25 represented by an italic digit, and the SD is 4.90 (see the first and second columns with no. 25).The final estimated SDs can be get through four steps (see Section4.2).  The following two tables (Tables 2 and 3) give the final estimated SDs for sinograms in Figure 3 added with noise whose SDs are from 5 to 15.Table 2
Estimated noise levels for the phantoms.
Figure3(a)(b)(c)Real SDPMW [18]FE [17]PMW [18]FE [17]PMW [18]FE [17]55.57*4.935.314.904.916.11*5.054.985.40*66.51*5.916.305.885.897.10*6.035.936.39*77.45*6.877.286.866.878.10*7.036.897.38*88.40*7.848.277.847.859.10*8.027.858.38*99.37*8.819.268.828.8310.08*9.028.829.37*1010.35*9.7710.259.819.8111.08*10.019.7910.37*1111.33*10.7511.2410.7910.8012.07*11.0110.7611.37*1212.2611.7212.2311.7711.7913.07*12.0011.7512.37*1313.2412.6713.2212.7512.7714.06*13.0012.7313.37*1414.2413.64*14.2113.7313.7515.06*14.0013.6914.37*1515.4114.6015.2014.7114.7416.06*14.9914.6915.38*Table 3
Estimated noise levels for the real projection data.
Figure3(c)(d)Real SDPMW [18]FE [17]PMW [18]FE [17]55.305.226.12*5.015.075.83*66.236.187.08*6.216.056.80*77.177.188.04*7.117.037.77*88.128.169.01*8.068.018.75*99.089.1410.00*9.048.999.73*1010.0410.1410.97*10.009.9910.71*1111.2511.1311.95*11.0910.9711.70*1212.1512.1212.93*12.0811.9512.69*1313.1013.1113.92*13.0412.9313.68*1414.0314.1114.91*13.9913.9114.67*1515.0915.1115.90*14.9814.9215.66*Table2 gives the final estimated SDs for phantoms (see Figures 3(a)–3(c)), “Real SD” is the real standard deviation (SD) for the added noise, and “PM,” “W,” and “FE” are the abs. of the proposed method, wavelets [18], and fast estimation [17], respectively. The bold digits represent the best estimate using three estimators while the digits with * are the worst estimate. In order to ignore small differences between two estimators, two digits will not be signed if the difference between two biases (difference of estimated SDs and the real SD) is smaller than 0.05. For example, in (b), “PM” and “W” have very similar estimated results, they are not signed.From Table2, we can see  that “PM” and “W” have very similar estimated results for (b), while “FE” has the worst estimated results in all noise levels; estimated results using “PM” are the best in most of noise levels, especially in relatively serious noises of (c) and the worst estimator is “FE” also, which overestimates noise levels for both subfigures.Table3 gives the final estimated SDs for real projection data (see Figures 3(d)-3(e)).  The all signs on Table 3 are the same as they are on Table 2.From Table3, we also can find that “PM” and “W” have very similar estimated results, while “FE” is the worst estimator in all noise levels for its overestimated SDs. For (e), “PM” has better performance in relatively serious noises, while “W” has better performance in relatively low noises.It should be indicated that the estimated results using “PM” on (a) are not satisfied especially in low noises. It reminds us, that there are also future works which can be done for improving the estimated results for the proposed method which will be discussed in Section7.


## 5.1. Data

Five test images are used in this paper: a thorax phantom acquired from a GE HiSpeed multislice CT scanner (see Figure2(a)), two phantom data produced on Matlab (see Figures 2(b) and 2(c)), and two images acquired with a 16 multidetector row CT unit (Somatom Sensation 16; Siemens Medical Solutions) using 120 kVp, 5 mm slice thickness (see Figures 2(d) and 2(e)).Test images. (a) Thorax phantom data acquired at 120 kVp, 150 mA. (b) A phantom data produced by Matlab. (c) Modified Shepp-Logan head phantom whose size is256×256. (d) A CT image of an abdomen of a 58-year-old man with 120 kVp, 150 mAs. (e) A CT image of a chest of a 62-year-old woman with 120 kVp, 150 mAs.
(a)(b)(c)(d)(e)For Figure2(a), the distance from the center of rotation (COR) to the curved detector is 408.075 mm. The detector array is on an arc concentric to the X-ray source with a distance between the X-ray source and the COR of 541.00 mm. The detector cell spacing is 1.0239 mm, and the slice thickness is 1.0 mm. Its projection data is shown on Figure 3(a).The relative projection data for the test images.
(a)(b)(c)(d)(e)The scanning parameters for Figures2(d) and 2(e) are: gantry rotation time, 0.5 second; table feed per gantry rotation, 24 mm; and pitch, 1 : 1. The CT doses were controlled by a fixed tube currents 150 mAs. Two respective sinogram data are shown on Figures 3(d) and 3(e).In addition, two projection data for Figures2(b) and 2(c) are shown in Figures 3(b) and 3(c), respectively. The projection data are produced in Matlab using radon transform for the equally distributed 90 angles for two phantoms.The SIG noises are added on the projection data in Figure3 directly. Since the interesting SIG noises are with small variances (SIG noises with large variances are not accepted in clinic situations), the variances are form 25 to 225 (the standard deviations are from 5 to 15).


## 5.2. Comparison Methods

Noise estimation is not a trivial task in image processing because of the complex structures for images. Some researchers suggest that the noise levels should be estimated by filtering images previously to suppress the image structures [17, 18, 24]. The methods compared with the proposed method are as follows.Wavelets
This method suppresses image structures using wavelets [18]. Since coefficients in the HH subbands of wavelets only preserving high frequency energy for images, the image structures, which mainly reside in low frequency coefficients, can be considered to be suppressed in the HH subbands. Thus based on robust median estimator, the level of the noise can be estimated as
(5.1)σN=median(|WHH(X)|)0.6745,
where X is the image, W is the wavelet operator, HH indicates the coefficients in the HH subbands, |·| is the absolute value, and median is the median estimator.Fast Estimation (FE)
This method is proposed in [17]. It only has two steps for noise estimation: the first step is filtering the image using Laplacian operator and the noise level is estimated using
(5.2)σN=π216(w-2)(h-2)∑|X⊗L|,
where w and h are the width and height of the image, respectively, ⊗ is the convolution operator, L is the Laplacian operator:
(5.3)[1-21-24-21-21].


## 5.3. Experimental Results

In this section, three methods,  proposed method (PM), wavelets, and Fast estimation (FE), are compared. In order to compare these three methods on a fair stage, the parameters used in these methods are fixed for all noise levels and all test images. According to this standard, the method proposed recently in [24] is not included because one of its parameter must be adjusted and selected with the assumption that the real noise levels are known.The wavelet used for  wavelets is Symlets with support 4. In summary, the parameters for PM are the squares used for finding the center point which are7×7 squares (see Section 4.1); the size of the square computing similarity between the considering pixel and the center point is 7×7 (see Section 3.2); the size for very large neighborhood is 41×41 (see Section 4.2 step 1); adding 50 samples each time to form the SDs sequence (see Section 4.2 step 3);  the threshold of absolute difference between adjacent standard deviations for finding the variances become ssuddenly 5; in order to avoid the influence of the outliers, the comparison for the variance is from the 21th addition for samples (see Section 4.2 step 4).Firstly, in order to show the role of variance analysis for the SDs sequence, we will use two Figures, Figures2(c) and 2(d), for analyzing the SDs sequence (see Table 1), where “Real SD” is the real standard deviation (SD) for the added noise, “No” is the times adding samples, “ESD” is the estimated standard deviation for each sample adding, “Diff” is the difference between the ith and i-1th estimated variances, the bold digits indicate the difference between two adjoint variances beyond the 5 (bigger diff), and the italic digits are the final estimated SDs.Table 1
Variance analysis to estimate the noise level.
Figure2(c)Figure2(d)Real SD5101551015NO.ESDDiffESDDiffESDDiffESDDiffESDDiffESDDiff214.930.269.821.0714.670.614.950.529.65−0.2614.423.75224.940.109.860.7814.783.004.960.159.690.9014.450.95234.91-0.299.82-0.7914.74-0.984.980.199.751.2014.501.53244.88-0.299.75-1.4214.62-3.484.980.199.751.2014.520.34254.900.229.811.1914.712.635.130.809.820.5314.684.80266.6119.6810.7619.6415.3418.715.200.699.942.4214.772.42275.250.5110.021.5214.873.17285.300.5410.040.4214.890.53295.947.2410.356.1714.982.733015.3210.38According to the steps in Subsection4.2, if the index for bigger diff is i, the final estimated SD is the i-1th SD. For example, in the second column, which is Figure 2(c) added noise with real SD 5, we can find that the bigger diff is 19.68 and its index is 26 (see the last no blank row of the first and second columns). Thus the index for final estimated SD is 25 represented by an italic digit, and the SD is 4.90 (see the first and second columns with no. 25).The final estimated SDs can be get through four steps (see Section4.2).  The following two tables (Tables 2 and 3) give the final estimated SDs for sinograms in Figure 3 added with noise whose SDs are from 5 to 15.Table 2
Estimated noise levels for the phantoms.
Figure3(a)(b)(c)Real SDPMW [18]FE [17]PMW [18]FE [17]PMW [18]FE [17]55.57*4.935.314.904.916.11*5.054.985.40*66.51*5.916.305.885.897.10*6.035.936.39*77.45*6.877.286.866.878.10*7.036.897.38*88.40*7.848.277.847.859.10*8.027.858.38*99.37*8.819.268.828.8310.08*9.028.829.37*1010.35*9.7710.259.819.8111.08*10.019.7910.37*1111.33*10.7511.2410.7910.8012.07*11.0110.7611.37*1212.2611.7212.2311.7711.7913.07*12.0011.7512.37*1313.2412.6713.2212.7512.7714.06*13.0012.7313.37*1414.2413.64*14.2113.7313.7515.06*14.0013.6914.37*1515.4114.6015.2014.7114.7416.06*14.9914.6915.38*Table 3
Estimated noise levels for the real projection data.
Figure3(c)(d)Real SDPMW [18]FE [17]PMW [18]FE [17]55.305.226.12*5.015.075.83*66.236.187.08*6.216.056.80*77.177.188.04*7.117.037.77*88.128.169.01*8.068.018.75*99.089.1410.00*9.048.999.73*1010.0410.1410.97*10.009.9910.71*1111.2511.1311.95*11.0910.9711.70*1212.1512.1212.93*12.0811.9512.69*1313.1013.1113.92*13.0412.9313.68*1414.0314.1114.91*13.9913.9114.67*1515.0915.1115.90*14.9814.9215.66*Table2 gives the final estimated SDs for phantoms (see Figures 3(a)–3(c)), “Real SD” is the real standard deviation (SD) for the added noise, and “PM,” “W,” and “FE” are the abs. of the proposed method, wavelets [18], and fast estimation [17], respectively. The bold digits represent the best estimate using three estimators while the digits with * are the worst estimate. In order to ignore small differences between two estimators, two digits will not be signed if the difference between two biases (difference of estimated SDs and the real SD) is smaller than 0.05. For example, in (b), “PM” and “W” have very similar estimated results, they are not signed.From Table2, we can see  that “PM” and “W” have very similar estimated results for (b), while “FE” has the worst estimated results in all noise levels; estimated results using “PM” are the best in most of noise levels, especially in relatively serious noises of (c) and the worst estimator is “FE” also, which overestimates noise levels for both subfigures.Table3 gives the final estimated SDs for real projection data (see Figures 3(d)-3(e)).  The all signs on Table 3 are the same as they are on Table 2.From Table3, we also can find that “PM” and “W” have very similar estimated results, while “FE” is the worst estimator in all noise levels for its overestimated SDs. For (e), “PM” has better performance in relatively serious noises, while “W” has better performance in relatively low noises.It should be indicated that the estimated results using “PM” on (a) are not satisfied especially in low noises. It reminds us, that there are also future works which can be done for improving the estimated results for the proposed method which will be discussed in Section7.


## 6. Conclusion

In this paper, we propose a new method to estimate noise for sinograms of low-dose CT. The proposed method can obtain estimated results both for phantoms and real projection data, especially in relatively serious noises, which demonstrate its potential for noise estimation of sinograms of low-dose CT.Based on the similarity sorting and variance analysis in a very large neighborhood whose scale is41×41, we can find enough similar samples to obtain reliable estimated results which make this proposed local method have very similar estimated results to the best existing global methods.In addition, avoiding convolution for suppressing image structures and relatively homogenous local structure makes the proposed method also be easily generalized to the more complex noises, such as, Possion noise and Gaussian compound noise. Thus the proposed method is also a promising method for real sinograms of low-dose CT.


## 7. Future Works

Although this paper proposes a new powerful method for simulated sinogram noise estimation, it may be improved as follows.(1)
How to find a center point in a large homogenous patch to ensure that there are enough points to obtain reliable estimation. We try to use multiresolution method to solve this problem.(2)
How to determine the number of the similar points according to the size of the samples and farther variance analysis.(3)
How to generalize the framework to more complex noise estimation, such as Poisson noise or Poisson and Gaussian compound noise.



---

*Source: 696212-2011-09-28.xml*