1. Introduction

Noise observed from an imaging sensor, when evaluated with a spatially uniform irradiance, exhibits both spatial and temporal variations that can be described as a three-dimensional (3D) random process. When the random process may be considered linear and shift invariant, one method of characterization is through the 3D power spectral density (PSD) [1,2]. Typically, the 3D PSD is a complex 3D function, making it difficult to draw significant conclusions without further processing or with additional metrics.

To aid in interpretation, in the early 1990s, NVESD engineers developed a model for the 3D power spectral density known as 3D noise [3]. The 3D noise model is a white noise approximation consisting of seven independent contributions that span the 3D correlation space (one 1D uncorrelated processes, three 2D uncorrelated processes, and one 3D uncorrelated process). Approximating the 3D power spectral density with the 3D noise model gives seven variances that can be loosely correlated with some underlying imager parameters. The spatial noise terms can be indicative of poor nonuniformity correction (NUC), while a high variance of any of the 1D temporal processes can be related to flicker [6].

A measurement of 3D noise is conducted by sampling the noise process along two spatial dimensions and a temporal dimension [4–7]. The number of samples available in the system being tested for the calculation of 3D noise depends both on the parameters of the system under test as well as the test geometry. For example, some testing methodologies employ a collimator to provide a uniform flat field, occupying only a portion of the field of view. Temporal limitations also can arise from either the source or system instability; for example, if the system under test employs any automatic features on a compressed time scale (automatic NUCing). In the case of these finite sampling geometries, errors will arise when using the classical 3D noise methodology because of the underlying assumptions of infinite sampling.

Recently, efforts have been made to correct for cases with finite sampling. In [8], a sampling correction was derived working directly with the 3D PSD. This method accurately captures the sampling artifacts, resulting in minimal errors on average. However, working with the PSD requires all samples to be considered, so the PSD can be significantly affected by outliers (an outlier, for example, could be a nonresponsive pixel that violates the typical assumption of Gaussianity of the underlying process). Another method in [9,10] was proposed to address the bias correction in sample space, where the effects of defective pixels can be removed from the calculation. Although this method provides a significant improvement in errors compared to the classic method, in the calculation of some of the processes it fails to capture instances of the coupling to the 2D random processing.

In addition to capturing sampling errors for the estimate in situations with limited sampling, a measure of uncertainty is also valuable. This is especially true in the current environment where the requirement margins for system performance are small and programmatic decisions are made from single measurements. To that extent, in this correspondence we derived and demonstrated the bounding confidence intervals for the 3D noise measurement. Additional details of the methods and calculation tools can be found in the MatLab functions on the Mathworks file exchange [1].

This correspondence is organized as follows: Section 2 presents the derivation of the 3D noise directional averaging and how the variance of all seven terms couple to the variance observed. Section 3 discusses confidence intervals as applied to the 3D noise measurement problem. Section 4 then shows the results of Monte Carlo simulations for both the mean estimation as well as the associated confidence intervals. Finally, we offer additional observations and comments on model assumptions, and then discuss additional preprocessing steps that can influence results.

2. 3D Noise Variance Estimation

The 3D noise model consists of three 1D random processes in each of the orthogonal directions: two spatial (horizontal and vertical); and one temporal (time). There are three 2D random processes describing the 2D correlations (time-vertical, time-horizontal, and vertical-horizontal), and one 3D random process (time-vertical-horizontal). The 3D noise model is defined as

Each of the independent noise sources is treated as a zero mean, delta correlated Gaussian random process completely characterized by its variance: ${\sigma }_t^2$, ${\sigma }_v^2$, ${\sigma }_h^2$, ${\sigma }_{tv}^2$, ${\sigma }_{th}^2$, ${\sigma }_{vh}^2$, ${\sigma }_{\mathrm{tvh}}^2$. An example cartoon depicting the 2D random process $vh$ sampled in 3D is shown in Fig. 1.

 

Fig. 1. Example cartoon depiction of the $vh$ random process sampled in 3D across $H$, $V$, and $T$ samples. Note that the $vh$ frame is independent of time. The results of the directional averaging of this process also are shown.

Download Full Size | PDF

Experimentally, the noise is measured as a 3D cube of sampled observations, as shown in Fig. 1. To obtain the variance of the independent random processes, we calculated statistics following directional averaging. The required seven equations for each of the unknown variances can be found by calculating the variance for each of: three orthogonal directional averages applied across one dimension, three different two dimensional directional averages, and the total variance. Because the individual random processes are independent, we can independently assess the effects of directional averages on each random process.

The unbiased estimate of the sample variance for each of the noise processes is defined:

The notation of $D_{x,y}$ represents the directional average of the $y$th noise component along the $x$th direction, where $x$ can be either one or two dimensions. Figure 1 shows a pictorial representation of the results of a directional average along with the relevant notation. Following the directional averaging, the variance calculated is denoted as ${\sigma }_{D_{t,v}}^2$.

The Appendix provided a complete derivation of the resulting variance following directional averaging for each of the seven independent random processes. Because the noise sources are independent, we can calculate the moments separately. Thus, the total variance of the directional averages is the sum of the individual terms:

3. Confidence Intervals for Variance Estimation

Recall that each of the measured variances represents a single realization of a random process. If this measurement were to be repeated, a different value would likely be found, and the total variability would be determined by a governing distribution. We will approximate the distribution of ${\widehat{\sigma }}_t^2$ and the other estimates to be Gaussian distributed, as they are a mix of many random processes. A Gaussian distribution is completely determined by its mean [as found by Eq. (6)] and its variance. To obtain the variance of the variance estimators, we must calculate the cross-covariance [11]:

In the Eq. (11), corr is the correlation, which can be used to aid in calculating the covariance. Recall that the random process ${\sigma }_{Dvh}^2$ and the others are each a summation of terms from the seven original independent random processes. Therefore, we can calculate the cross covariance for each one independently, then sum to obtain the total variance. To derive the cross covariance, we must derive the distribution of calculated variance for each of the independent random processes following each of the different directional averages.

Let $X_1,X_2,\dots X_n$ be random samples from a Gaussian distributed random process with a mean $\mu$, population variance ${\sigma }^2$, and sample variance $S^2$. Then the random variable:

To calculate the total cross covariance, we must know the scaling and degrees of freedom for each variance following the directional averaging, for each of the seven random processes. Below we present the scaling and degrees of freedom for the 1D temporal random process, the full details of the remaining terms can be found in the MatLab files [1].

To calculate the cross covariance of the temporal random process, we also must take into account the correlations remaining between the different directional averaging. Below is the cross correlation of the variance following the directional averaging, for the temporal random process:

Using the above expressions, the cross covariance contribution from the temporal random is found from:

From the cross covariance in Eq. (10), the confidence intervals can be found, again assuming a Gaussian distribution. The confidence bounds for the conventional two-sided confidence of $100\%(1-\alpha )$ are found from the PDF by solving for the confidence intervals $[{\mathrm{CI}}_L,{\mathrm{CI}}_H]$:

For a Gaussian PDF, the bounds are symmetric; therefore, only the difference to the mean is required.

4. Monte Carlo Simulations

To test the derived expressions and the method of calculating confidence intervals from the above, we performed a series of Monte Carlo simulations across a variety of different sampling conditions. The underlying noise simulations used a Gaussian random number generator for the distribution creation. The notional values for the different 3D noise terms were selected as ${\sigma }_t={\sigma }_{tv}={\sigma }_{th}=10$, ${\sigma }_v={\sigma }_h=20$, and ${\sigma }_{vh}={\sigma }_{\mathrm{tvh}}=50$. The simulations were carried out by varying percentages of sampling of a notional sensor size of $640\times 480$ pixels using between 30 and 240 frames. For each sampling geometry, 10,000 independent noise cube realizations were generated and then processed with directional averaging and the 3D noise calculation. In addition to the classic method of processing and the derived method, the bias correction proposed by [10] was also calculated and reported.

Figure 2 demonstrates the variability of measured variance, the probability distribution of obtained variances using the derived mixing matrix in Eq. (6) for $V=24$, $H=32$, and $T=30$.

 

Fig. 2. Plot of distributions of observed variance compared to Gaussian assumption for $V=24$, $H=32$, and $T=30$; estimated (A) temporal variance, (B) vertical variance, (C) temporal-vertical variance, (D) vertical-horizontal variance, and (E) temporal-vertical-horizontal variance.

Download Full Size | PDF

From symmetry considerations, Fig. 2 omits the $th$ and $h$ distributions as $th$ and $tv$ have similar distributions, as do $v$ and $h$ (both differing slightly due to different number of samples along $V$ and $H$). As can be seen in the figure, the distributions of values observed vary depending on the magnitude of the noise term as well as the number of samples. The average of the distribution agrees with the theoretical prediction. Also shown is the predicted distribution using the coupled Gaussian assumption derived for the confidence interval calculation (red solid). It is this distribution that is used to provide the confidence bounds. For comparison, the assumption of no coupling (the black dashed line) using a chi-squared distribution is also shown; for these particular values, there is limited coupling and the independent assumption works well for most of the noise components.

To summarize the Monte Carlo simulations first moment results, Table 1 presents the percent error on the average across the 10,000 realization for the different mixing matrices for a variety of sampling geometries.

Table 1. Percent Error in Average Variance from the Monte Carlo Simulation. Class Refers to Mixing with Eq. (7), UB0 Refers to Eq. (9), and the New Refers to the Mixing Derived of Eq. (6)

View Table | View all tables in this article

As can be seen in Table 1, in general the percent error decreases with an increase in the sampling size. However, there still remains a bias in the classic method that is on the order of a few percent. The bias correction in [10] does improve the error as compared to the classic method; however, the new method presented captures the bias to a fraction of a percent. Although not reported here, using the method of [8] using the 3D PSD gave the same answer as the proposed mixing of Eq. (6).

To validate the assumption of a Gaussian distribution for the confidence intervals, we also can use the 10,000 realizations to test what percentage of the calculated confidence intervals enclose the true population variance. The results can be seen in Table 2 using 90% confidence intervals ($\text{alpha}=0.1$, 90%).

Table 2. Observed Percentage of Confidence Intervals Capturing the Population Statistic

View Table | View all tables in this article

As Table 2 shows, as the number of samples increases the distributions become more Gaussian and the method of calculating the confidence intervals captures the uncertainty in the measurement. Even when the sampling is limited, the confidence intervals are still within approximately 5% of the observed. It is important to note that the above table is only for a single distribution of variance for the underlying random processes.

From the expressions for the cross covariance, coupling can be observed; for example, a large amount of noise in $vh$ will leak into the $v$ and $h$ variance calculations. As a demonstration, consider the situation where all the noise terms in Eq. (9) are equal to 10 except for the $vh$ term, which has a standard deviation of 100. The 3D random process is sampled with $T=30$, $V=24$, and $H=32$. The average values (10,000 realizations) correspond to the population values, but due to the coupling the confidence intervals can be significantly larger for the terms experiencing coupling

As can be seen in Eq. (20), the $v$ and $h$ terms will have a negative value for their confidence bounds. We can further demonstrate these bounds by plotting the distributions.

As shown in Fig. 3.B, there are situations where the calculated variance is negative; this is in addition to situations where the confidence bound is negative. A negative variance should be seen as a situation where there was severe coupling, and caution should be used in interpreting the significance.

 

Fig. 3. Plot of the estimated variance of the $v$ random process in the case of limited coupling (A) ${\sigma }_v=10$, ${\sigma }_{vh}=10$, and severe coupling (B) ${\sigma }_v=10$, ${\sigma }_{vh}=100$. Blue circles correspond to Monte Carlo, black dashed is independent chi-squared distribution, and red solid is coupled variance derived Gaussian distribution.

Download Full Size | PDF

The noise is typically reported through the use of the standard deviation as opposed to the variance for 3D noise (to retain the same units as the input, typically in units of counts or temperature). To accommodate a negative variance, we propose to retain the sign and take the square root as

For example, the standard deviation and confidence intervals for the large coupling case above can be found as

The severity of the negative confidence bounds depends on the sampling as well as the strength of the coupling between the different processes. In addition to strong coupling, a negative variance value also could indicate that the underlying process is not well approximated by the 3D noise model.

5. Summary and Conclusions

We have derived and demonstrated the full sampling corrected 3D noise measurement and the corresponding confidence bounds. We accomplished this analysis by decomposing the temporal and spatial variances following the directional average into each noise component, and demonstrated it through the use of Monte Carlo simulations of noise cubes. Our results provide a significant improvement in the average (biased) error when small sampling sizes are used. In the limit of infinite sampling, the measurement matrix developed here converges to the original method.

An important property of our method is that it allows for the effects of defective pixels to be ignored. A large outlier violates the assumption of the underlying Gaussian processes and can skew the variance estimation. The method presented is conducted in the spatial domain and is based on statistical calculations of the variance. It allows for defective pixels to be ignored by simply removing them from the calculation.

The methods for both the sampling correction as well as the confidence intervals are determined by the average 3D noise variances and the corresponding sampling size; therefore, an a posteriori calculation of both the sampling correction and confidence bounds can be performed. The confidence intervals also can be used in a predictive fashion; for example, if the expectation of the noise is known, the confidence intervals can be used to provide an estimate of the necessary sampling required.

Our future work will examine the signal intensity transfer function (SiTF), which is used to convert noise in counts into noise of temperature, and its contribution of uncertainty to the 3D noise process. Additionally, the characterization of defective pixels and how to report their influence also is important, requiring adequate definitions and agreed-upon metrics. The 3D noise model we have discussed assumes a summation of seven independent spectrally white random processes; future work also could aim at characterizing nonwhite effects.

Finally, in support of this work a series of four Matlab functions have been provided to the Mathworks file exchange. These scripts allow users to simulate a noise cube, calculate the seven independent variances, and calculate the confidence intervals. An additional script provides examples of how the functions can be called under different circumstances; for example, correcting a past measurement, or iteratively calculating the noise cube by chunking the data.