A Method for Quantification of Noise Non-Uniformity in Computed Tomography Images: A Computational Study

filtered noisy phantoms ( I BF ), and the magnitudes of the resulting noise for each gradient were computed. The noise-gradient dependency (NGD) curve was used to display the dependency of noise magnitude on image gradient in the non-uniform noise. It was found that for uniform noise, the magnitude of noise was constant for all gradients. However, for non-uniform noise, the measured noise was dependent on the gradient levels and on the strength of the BF for every ground truth noise ( σ G ). It was found that the noise magnitude was large for the large gradients and decreased with the magnitude of the image gradient.


Introduction
Noise is defined as stochastic fluctuations of the pixel values in an image due to uncertainties in the image production [1][2][3]. The noise can be quantum noise or electronic noise or it can result from the image generating process, i.e. from image reconstruction algorithm [4][5][6]. Fluctuations of the pixel values in an image due to anatomical variations are not considered as noise. Noise assessment generally uses the standard deviation (SD) [7] or the noise power spectrum (NPS) [8,9] of pixel values in a homogeneous region of interest (ROI). The SD of pixel values quantifies the noise magnitude [10], while the NPS characterizes the noise texture [11,12], i.e. visual impression of the noise, whether it is fine or coarse. Noise evaluation using the SD or NPS values is widely accepted by the medical physics community for uniform noise [13][14][15][16][17][18].
In earlier methods of CT reconstruction, a filtered back-projection (FBP) method [19,20] was used to produce images with uniform noise. Nowadays, an iterative reconstruction (IR) method is generally used for CT image reconstruction in order to produce a high-quality image with low dose. The IR results in images with non-uniform noise [21][22][23][24][25]. In addition, many non-linear filters, such as the bilateral filter (BF) [26,27] and the nonlocal mean (NLM) filter [28,29] are implemented in CT image, where filtering is performed aggressively in homogeneous regions with zero gradients, but less aggressively in non-homogeneous regions or edge areas with high gradients [30,31]. Image noise obtained from non-linear filters is also non-uniform noise, having characteristics of low noise in homogeneous regions but the boundaries between objects still appear sharp [32].
In non-uniform noise, the noise may be smaller in the homogeneous regions than it is in regions near the edges [33]. Thus, images may have not only one noise magnitude, but have many different noise magnitudes depending on the position and environment, whether they are in a homogeneous region or in an area near an edge. As a result, assessment of noise only in homogeneous/zero gradient regions does not comprehensively describe the actual nature of noise [33]. A new metric needs to be developed to quantify noise nonuniformity.
To date, two metrics have been introduced to quantify noise non-uniformity, and both are based on the spatial noise-map (SNM) [33]. The SNM is computed as the standard deviations of a sliding ROI within the image [34,35]. The first metric is the noise spatial non-uniformity index (NUI) [36] and is for characterizing the degree of noise fluctuation globally across the field of view (FOV). It can be measured from images of a homogeneous phantom or a special structured phantom. The NUI is computed as the SD of the SNM. The second metric is the noise inhomogeneity index (η) [37] and is for characterizing highly irregular spatial distributions of noise, found in iterative reconstruction of the structured phantoms [30,37]. This metric is computed from two distinct peaks of the histogram of the SNM of the structured phantom. Specifically, it is calculated as the separation of the two peaks divided by the height difference of the two peaks.
However, these two metrics do not completely characterize the noise fluctuations due to the magnitude of the image gradient is not considered. Therefore, we present a novel method for quantifying noise non-uniformity using a noisegradient dependency (NGD) curve, and propose 1D and 2D computational phantoms designed for this purpose. We then validate the effectiveness of the NGD in evaluating noise non-uniformity in filtered images using the non-linear BF.

Quantification of noise non-uniformity
In non-uniform noise, the magnitude of the noise changes according to the image gradient (G), i.e. the gradient in pixel values throughout the image. We refer to the curve displaying the magnitude of noise as a function of the magnitude of the image gradient as the noise-gradient spectrum (NGS). The noise for each gradient (IG) can be calculated if the magnitude of the gradient itself is known. The method is simply to subtract the gradient (G) itself from the noisy images that have a certain gradient (IN).
(1) This is illustrated in Fig 1. The first row of Fig 1 shows a homogeneous 1D phantom (i.e. zerogradient, G = 0) as is commonly used to date. Because the gradient is zero, there is no need to subtract the gradient from the noisy phantom. The second row shows noisy 1D phantom at G = 1 HU/pixel and the third row shows noisy 1D phantom at G = 2 HU/pixels. The gradients (G), shown in (b), are subtracted from the noisy 1D phantoms (IN), shown in (a), to give the noise (IG), shown in (c), from which the noise magnitude ( ) can be directly calculated using the standard deviation equation: (2) with (3) where i indicates number element of data phantom with maximum N data. The equations (1)-(3) can be extended for 2D phantom with row-element (i) and column-element (j).

Design of the phantoms
The phantoms used to realize the NGD calculations were in the form of computational phantoms developed using MATLAB software (Mathworks, MA). To easily comprehend the concept of NGD, we firstly developed 1D gradient phantom and after that we implemented the NGD in 2D gradient phantom.
The design of the 1D gradient phantom is shown in Fig 2. Each gradient consists of 512 pixels. (a) to (f) are processes to obtain the 1D gradient phantom. For an object gradient of 1 HU/pixel (Fig 2(b)), at the left end was 0 HU and at the right end was 512 HU. This was because every shift of 1 pixel to the right, the pixel values increased 1 HU. For an object gradient of 2 HU/pixel (Fig 2(c)), at the left end was 0 HU and at the right end was 1024 HU (2*512 HU). The same pattern was applied for all object gradients of 3 to 48 HU/pixel. After that all gradients are added to make 1D gradient phantoms. The design of the 2D gradient phantom is shown in Fig 3. The phantom is a circle with a diameter of 512 pixels and with a homogeneous background of 100 HU (This is assumed that phantom is made from acrylic-like material with CT number of around 100 HU). The phantom is located in air with a pixel value of -1,000 HU. There are 49 objects in the phantom, with gradient values from 0 to 48 HU/pixel. Each gradient object has a size of 30 × 30 pixels, separated from the next by a distance of 5 pixels. The initial pixel value of every gradient object (top left corner) is -500 HU. The right side of  Illustration of the concept of a gradient phantom. The first row indicates 1D phantom with a gradient of 0 HU/pixel, the second row indicates 1D phantom with a gradient of 1 HU/pixel, and the third row indicates 1D phantom with gradients of 2 HU/pixels. The gradients (b) are subtracted from the noisy 1D phantom (a) to obtain the noise (c).

Research design
The research design for calculating the NGD for the 1D phantom is shown in    IN was filtered using a non-linear filter, namely the 1D bilateral filter (1D-BF) so that a filtered noisy phantom (IBF) was generated. The 1D-BF was calculated using: (4) where k is a normalization constant and was calculated using: (5) where represents the geometric spread of the BF. By increasing , more neighboring pixels are utilized for denoising, resulting in increased filtering [26]. The σ represents the pixel intensity spread of BF [26].
x indicates x-axis and i indicates pixel kernel of BF. In this study, 1D-BF filtering was performed for each noise with 3 variations of σ, i.e. 1×, 0.5×, and 0.25× ground truth noise (σG). For example, for σG of 10 HU, the magnitudes of σ values in the 1D-BF were 10, 5, and 2.5 HU.
Filtering using the 1D-BF produced noise that varies according to the magnitude of the object gradients (IBF). For objects with gradients below the value of σ, filtering works very aggressively. As the gradient increases, filtering becomes less aggressive. For certain level of gradients, the 1D-BF was no longer able to filter because the gradient of pixel values was outside the 1D-BF range.
To calculate the noise non-uniformity, the gradient of the object (G) was subtracted from the filtered noisy image (IBF), and standard deviations in every ROI were then computed. Fig 4 showed that the unfiltered original image had uniform noise, i.e. the magnitude of noise was relatively flat for all gradients. Conversely, the filtered noisy image, had non-uniform noise, i.e. the magnitude of noise was large for large gradients, and decreased for small gradients.
Similar to the 1D phantom, Gaussian noise of various levels (i.e. noise levels (σG) of 5, 10, 25, and 50 HUs) (Fig 5(b)) was added to the 2D phantom (Fig 5(c)). In addition, to simulate the real form of noise from images obtained in a real CT scanner, a homogeneous water phantom image, the AAPM CT performance phantom for homogeneous water sections (Fig 5(d)), was used. The phantom was scanned in axial mode in a GE LightSpeed CT with two tube currents, i.e. 50 mA and 100 mA, a tube voltage of 120 kVp, and a slice thickness of 5 mm. The homogeneous phantom image was added to the 2D gradient image, in order to obtain a 2D gradient image with the real noise from the CT (Fig 5(e)).  The design for calculating the NGS in the 2D phantom is shown in Fig 6. The 2D phantom was then filtered using a 2D-BF, with 3 variations of pixel intensity spread (σ), i.e. 1×, 0.5×, and 0.25× of ground truth noises (σG). The 2D-BF was carried out using equation (6). (6) where k is the normalization constant using: (7) Similar to the 1D phantom, the original gradient phantom (G) was subtracted from the filtered noisy object (IBF) to calculate the noise non-uniformity in the 2D phantom. Fig 6 shows that uniform noise was obtained in the unfiltered image, while nonuniform noise was obtained in the filtered noise phantom. To calculate noise, a 30 × 30 pixels ROI was used within each object, and the noise magnitude was calculated using the standard deviation formula for the 49 ROIs.

Results and Discussion Uniform noise
The NGD curves for noise uniformities for various gradients are shown in Fig 7(a) for the 1D phantom, and in Fig 7(b) for the 2D phantom. Measured noises are constant for all gradients depending on the given ground truth noises (σG), i.e. 5, 10, 25, and 50 HUs. It is found that in the 1D phantom, the noise fluctuations are very small and all measured noises are within ± 5% of ground truth noises. Noise fluctuations are slightly greater in the 2D phantom, with some measured noises fluctuating more than ± 5% from the ground truth noises. All of the noise measurements are uniform for all gradients, from 0 to 48 HU/pixel for both the 1D and 2D phantoms.
The NGD curves for noise uniformities for CT images scanned with two tube currents of 50 mA and 200 mA for various gradients are shown in Fig  8. The noises are constant at about 10 HU for 50 mA, and about 5 HU for 200 mA. It is found that the noise fluctuations are relatively large and sometimes more than ± 10% of its average values.

Non-uniform noise
The NGD curves for the noisy 1D phantom filtered with the 1D-BF using various pixel intensity spreads (σ), produces non-uniform noises, which are shown in Fig 9. As expected, at high gradients, the measured noises are similar to the given ground truth noises (σG). For a gradient of equal or smaller than σ, the measured noises are much smaller than σG, while for a gradient greater than σ, the measured noises are fairly similar to the values of σG. Fig 9 shows that if the noisy 1D phantom is filtered with the 1D-BF using σ = σG, then the measured noise levels at zero-gradient are about 75% of σG. At a gradient level of around σ, it produces measured noise levels of about 85-90% of σG, and at a gradient value of around 2σ, the measured noise levels are equal to σG. If the noisy 1D phantom is filtered with the 1D-BF using the σ = 0.5 σG, the measured noise levels at zero-gradient are around 90% of σG. And if the noisy 1D-phantom is filtered with 1D-BF using σ = 0.25 σG, the measured noise levels for all gradients are equal to σG. i.e. no significant filtering.   The NGD curves for the noise non-uniformities in the 2D phantom for various gradients are shown in Fig 10. Similar noise patterns are obtained for the 2D phantom as for the 1D phantom, except that the magnitude of the noise fluctuations with the 2D phantom are greater. The same noise patterns are also found in CT images scanned with two tube currents, 50 mA and 200 mA (Fig 11).  This study introduced and validated new methods for quantifying noise non-uniformity using an NGD. To the best of our knowledge, this is the first study introducing a concept of an NGD for quantifying non-uniformity. The challenge of calculating the NGD is the availability of gradient phantoms. As the first study, we only developed 1D and 2D gradient phantoms computationally. Physical realization of the designed phantom is a challenging task. However, it may be possible to realize a physical gradient phantom using the partial volume artifact (PVA), i.e. by scanning an object in an oblique position, so that a gradient of HU can be obtained.
With these two computational gradient phantoms, various amounts of noise can be added as ground truth noise and the resulting noisy phantoms can be filtered using non-linear filters, such as bilateral filters (BF) [26], non-local mean filters [28], or others. In this study, we used the BF to produce a non-uniform noise image. The effectiveness of a non-linear filter can be comprehensively assessed with these computational phantoms. We found that with uniform noise, the measured noise is independently with the image gradient. In contrary, in nonuniform noise, the measured noise varies with the magnitude of the image gradient, increasing with the magnitude of the image gradient, as expected. The noisy phantoms filtered using the BF produce a noise magnitude that is strongly influenced by the image gradient and magnitude of the pixel intensity range (σ) of the BF.
Both the 1D and 2D filtered noisy phantoms produces similar patterns, albeit with larger noise fluctuations in the 2D phantom. If a small ground truth noise (σG) is given the 2D phantom, for example σG = 5 HU, then the filtered 2D noisy image using a pixel intensity range of σ will produce the noise magnitude equal to the σG at gradients above 10 HU/pixel. However, we found that the measured noise is about 20% smaller than σG. This is due to the limitation of the 2D gradient phantom design. For a certain pixel, not all of the neighbor pixel values differ according to the gradient. For example, in pixel (i, j), the four of the neighboring pixels, i.e. (i-1, j), (i, j-1), (i + 1, j) and (i, j + 1), differ by an amount equal to the gradient magnitude from the pixel of (i, j), while two other neighboring pixels, i.e. (i-1, j-1) and (i + 1, j + 1), differ by 2 gradient magnitudes, and the other two neighboring pixels, i.e. (i + 1, j-1) and (i -1, j + 1), have pixel values equal to the pixel at (i, j). As a result, the total noise decreases slightly. As a result, the 2D phantom is not as perfect as the 1D phantom although it can be used to obtain reliable NGD for non-uniform noise images.

Conclusions
A single metric such as SD or NPS does not comprehensively characterize noise non-uniformity, since the magnitude of the noise varies with the gradient of the image. To remedy this, we have introduced noise-gradient dependency (NGD). We have developed 1D and 2D computational phantoms to implement and verify our new methodology. We find that noisy 1D and 2D phantoms, filtered with a bilateral filter, produce non-uniform noise and the NGD curves are able to determine the dependencies of the noise magnitude on the gradient level. We found that in the filtered noisy phantoms, the measured noise varies along with the magnitude of image gradient. Using the NGD enables an evaluation of the effectiveness of the available algorithms used in noise reduction and iterative reconstruction. However, realizing a physical gradient phantom remains challenging.