Signal-to-noise enhancement

Abstract

This disclosure relates to processing a spectral dataset, such as a hyperspectral image or a large collection of individual spectra taken with the same spectrometer, to increase the signal-to-noise ratio. The methods can also be used to process a stack of images that differ by acquisition time rather than wavelength. The methods remove most of the sensor background noise with minimal corruption of image texture, anomalous or rare spectra or waveforms, and spectral or time resolution.

Claims

1. A signal-to-noise enhancement method, comprising: receiving from a spectrometer a spectral dataset that comprises a multiplicity of spectra measured with the spectrometer, in which the number of spectra exceeds the number of measured wavelength bands of the spectrometer; estimating the covariance matrix of background noise associated with the measured spectra; transforming the multiplicity of spectra to Maximum Noise Fraction (MNF) components, each component comprising the coefficients resulting from projection of the multiplicity of spectra onto an individual MNF eigenvector; calculating nonlinear mathematical functions, associated with each individual MNF eigenvector, that convert the values of the MNF components of a given spectrum into estimated values of the corresponding noise-free MNF components of the given spectrum, in which the nonlinear mathematical functions are given by S(y)=exp(|y/s|.sup.p)/[r exp(y.sup.2)+exp(|y/s|.sup.p)]; applying the inverse MNF transform to convert the estimated noise-free MNF components of the multiplicity of spectra back into the original units of the spectra; and generating an output spectral dataset from the converted estimated noise-free MNF components of the multiplicity of spectra.

2. A signal-to-noise enhancement method, comprising: receiving from a spectrometer a spectral dataset that comprises a multiplicity of spectra measured with the spectrometer, in which the number of spectra exceeds the number of measured wavelength bands of the spectrometer; estimating the covariance matrix of background noise associated with the measured spectra; transforming the multiplicity of spectra to Maximum Noise Fraction (MNF) components, each component comprising the coefficients resulting from projection of the multiplicity of spectra onto an individual MNF eigenvector; calculating nonlinear mathematical functions, associated with each individual MNF eigenvector, that convert the values of the MNF components of a given spectrum into estimated values of the corresponding noise-free MNF components of the given spectrum, in which the nonlinear mathematical functions are given by S(y)=1/[r exp(y.sup.2)+1]; applying the inverse MNF transform to convert the estimated noise-free MNF components of the multiplicity of spectra back into the original units of the spectra; and generating an output spectral dataset from the converted estimated noise-free MNF components of the multiplicity of spectra.

3. A system, comprising: a spectrometer that provides a multiplicity of spectra, in which the number of spectra exceeds the number of measured wavelength bands of the spectrometer; and digital image processing circuitry which, in operation: estimates the covariance matrix of background noise associated with the measured spectra; transforms the multiplicity of spectra to Maximum Noise Fraction (MNF) components, each component comprising the coefficients resulting from projection of the multiplicity of spectra onto an individual MNF eigenvector; calculates nonlinear mathematical functions, associated with each individual MNF eigenvector, that convert the values of the MNF components of a given spectrum into estimated values of the corresponding noise-free MNF components of the given spectrum, in which the nonlinear mathematical functions are given by S(y)=exp(|y/s|.sup.p)/[r exp(y.sup.2)+exp(|y/s|.sup.p)]; applies the inverse MNF transform to convert the estimated noise-free MNF components of the multiplicity of spectra back into the original units of the spectra; and generates an output spectral dataset from the converted estimated noise-free MNF components of the multiplicity of spectra.

4. A system, comprising: a spectrometer that provides a multiplicity of spectra, in which the number of spectra exceeds the number of measured wavelength bands of the spectrometer; and digital image processing circuitry which, in operation: estimates the covariance matrix of background noise associated with the measured spectra; transforms the multiplicity of spectra to Maximum Noise Fraction (MNF) components, each component comprising the coefficients resulting from projection of the multiplicity of spectra onto an individual MNF eigenvector; calculates nonlinear mathematical functions, associated with each individual MNF eigenvector, that convert the values of the MNF components of a given spectrum into estimated values of the corresponding noise-free MNF components of the given spectrum, in which the nonlinear mathematical functions are given by S(y)=1/[r exp(y.sup.2)+1]; applies the inverse MNF transform to convert the estimated noise-free MNF components of the multiplicity of spectra back into the original units of the spectra; and generates an output spectral dataset from the converted estimated noise-free MNF components of the multiplicity of spectra.

5. The method of claim 1, in which the multiplicity of spectra comprise a hyperspectral image.

6. The method of claim 2, in which the multiplicity of spectra comprise a hyperspectral image.

7. The system of claim 3, in which the spectrometer comprises a hyperspectral imaging sensor.

8. The system of claim 4, in which the spectrometer comprises a hyperspectral imaging sensor.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

(1) Other objects, features and advantages will occur to those skilled in the art from the following detailed description, and the accompanying drawings, in which:

(2) FIG. 1A illustrates steps in a signal-to-noise enhancement method.

(3) FIG. 1B is a functional block diagram of a system that can be used to accomplish a signal-to-noise enhancement.

(4) FIG. 2 illustrates steps in an example of the shrinkage function estimator of FIG. 1A.

(5) FIGS. 3A-3D comprises plots of the MNF component optimal-estimate noise-free absolute value, |x(y)|, as a function of the original absolute value, |y|, from Eqs. (3), (8) and (10), with s= in the generalized Laplacian prior PDF, Eq. (7).

(6) FIG. 4 comprises a plot of the value of r in Eq. (9) versus MNF eigenvector number, based on a histogram analysis of an 82-band long-wavelength infrared (LWIR) hyperspectral image taken with an airborne Telops HyperCam sensor.

(7) FIG. 5 comprises plots of LWIR spectral radiance in units of microflicks versus micrometers from the hyperspectral image of FIG. 4.

(8) FIG. 6 comprises plots of the same pixel spectra as in FIG. 5 but after processing with the present signal-to-noise enhancement method.

(9) FIG. 7. illustrates steps in another signal-to-noise enhancement method.

(10) FIG. 8 comprises a frame from a 30-frame video after Gaussian noise was added.

(11) FIG. 9 comprises the same frame as in FIG. 8 but after the video with added noise was processed with the present method.

DETAILED DESCRIPTION

(12) The signal-to-noise enhancement of this disclosure uses the MNF transform method for separating signal and noise components in multidimensional datasets in combination with the method of applying nonlinear multiplicative factors, called shrinkage functions, to the MNF components of each sample or pixel. In one example, a separate shrinkage function is defined for each MNF component using the Bayesian signal estimation formula. In essence, this method suppresses MNF transform values that are within the noise level, while retaining values that exceed the noise level, which arise from anomalous or rare samples.

First Example

(13) In a first example of a signal-to-noise enhancement shown in FIG. 1A, a plurality of spectra 12 from a spectrometer 10 are fed into the processing system, comprised of items 14, 16, 18 and 20. The number of spectra must be greater than the number of wavelength measurements, or bands, in the spectra, denoted N. The invention outputs de-noised spectra 22. The spectrometer may be a hyperspectral imager. FIG. 1B illustrates system 80 that can be used to accomplish the signal-to-noise enhancement of FIG. 1A. Spectrometer 82 provides an output (a plurality of spectra) to processing system 84. Processing system 84 may be a computer that is programmed to carry out the signal-to-noise enhancement. The output (which can be the input spectra, but de-noised) are provided to output device 86. Output device 86 may be a printer or a computer display device, for example.

(14) An estimate of the NN noise covariance matrix 14 may be obtained by a number of different methods, including by analysis of the data itself or from a specification of the RMS noise supplied by the sensor manufacturer. Specific methods are described in Green [1988], Phillips [2008] and Bjorgan [2015]. Typically, the noise in different wavelength bands is considered to be uncorrelated, in which case the noise covariance is a diagonal matrix of noise variances for each wavelength band. If only an average noise variance is available, then the elements of the diagonal matrix are identical.

(15) The MNF transform is described in numerous papers, including Phillips [2008] and Bjorgan [2015]. It can be regarded as a sequence of two steps. In the first step, which is a noise-whitening operation, the data are multiplied by the square root of the noise covariance matrix inverse. If the noise covariance matrix is diagonal, the whitening step is a spectral weighting of the data, in which the data are divided by the noise standard deviation. In the second step, a Principal Component (PC) transform is applied to the data. The MNF-transformed data consist of the projections of the mean spectrum-subtracted data onto the eigenvectors of the PC transform. We assume without loss of generality that the eigenvectors have been unit-normalized.

(16) The two mathematical steps outlined above are equivalent to a single data transformation. The first step is to subtract the mean spectrum from the data. Using the notation in Phillips [2008], let s denote the NN covariance matrix of the spectral data and N denote the NN covariance matrix of the noise. The matrix of N eigenvectors, V, is the solution to the equation
.sub.S.sub.N.sup.1=VV.sup.1(1)
where is the diagonal matrix of N eigenvalues that correspond to matrix V. As is usual, for the purposes of this discussion we assume that the eigenvectors and eigenvalues, which are positive, are ordered in decreasing eigenvalue size.

(17) The projection of the mean-subtracted data onto the N eigenvectors of the V matrix yields a set of N coefficients for each spectrum, which we denote y.sub.j for a given eigenvector j. We shall refer to the y.sub.j as the j.sup.th MNF component. If the original data comprise a hyperspectral image (data cube), then each of the N MNF components y.sub.j comprises an image. Collectively, the y.sub.j comprise the MNF-transformed spectral dataset 16.

(18) Next, we seek optimal estimates of the true MNF components, x.sub.j, which are uncontaminated with background noise, n.sub.j, based on the observed y.sub.j=x.sub.j+n.sub.j. This is done by defining a shrinkage function, denoted S.sub.j(y.sub.j), which is a multiplicative factor that converts the observed noise-containing MNF component y.sub.j into an estimate of the noise-free MNF component. That is, the optimal estimate x.sub.j is given by
x.sub.j(y.sub.j)=S.sub.j(y.sub.j)y.sub.j(2)

(19) The preferred embodiment of the shrinkage function estimator procedure 18 is shown in FIG. 2. The basic approach is described in Simoncelli [1996], who however applied it to wavelet transforms rather than MNF transforms. The first step is to estimate prior probability distribution functions (PDFs) 26, denoted P.sub.j(x.sub.j), associated with the noise-free MNF components. There a number of ways by which this may be done. In one example, Simoncelli [1996] describes a blind algorithm by which the two parameters of a generalized Laplacian prior PDF (Simoncelli [1996] Equation 2) corrupted by Gaussian white noise may be estimated from the second and fourth moments of the y.sub.j distribution across the dataset. These moments are but one example of MNF component statistics 24 that can potentially be employed to estimate the prior PDFs.

(20) The optimal estimate of a noise-free MNF component value, denoted x.sub.j, is its mean value in the posterior PDF, which includes the added noise. This mean value, which provides an unbiased-least squares estimate of x.sub.j given y.sub.j, is computed from the Bayesian signal estimation formula 28, found in Simoncelli [1996]. We assume that the noise PDF is the unit-variance Gaussian function, G, and define the function Q.sub.j(x.sub.j)=x.sub.jP.sub.j(x.sub.j). Then the Bayesian signal estimation formula may be compactly written as a quotient of convolutions,
x.sub.j(y.sub.j)=(G*Q.sub.j)(y.sub.j)/(G*P.sub.j)(y.sub.j)(3)

(21) Using Eq. (2), the shrinkage function 30 is then given by
S.sub.j(y.sub.j)=x.sub.j(y.sub.j)/y.sub.j.(4)

(22) Since analytical solutions to Eq. (3) do not exist for many useful functional forms for the prior PDF, Simoncelli [1996] proposes evaluating the convolutions numerically. However, analytical approximations are useful for developing approximate shrinkage functions 30 that provide good de-noising results. Some examples are described below.

(23) For the small-eigenvalue j's, which are dominated by noise, it is assumed that the vast majority of the spectra have MNF component values x.sub.j that are very close to zero; that is, they belong to an extremely narrow PDF. On the other hand, a small minority of the spectra, comprising anomalies and rare spectra, are associated with a relatively broad PDF. Let us therefore approximate the PDF of the true MNF coefficients P.sub.j(x.sub.j) as a weighted sum of a broad distribution B(x.sub.j) and a Dirac delta function D(x.sub.j), the latter representing the limit of an extremely narrow distribution:
P.sub.j(x.sub.j)=w.sub.jB.sub.j(x.sub.j)+(1w.sub.j)D(x.sub.j)(5)
where w.sub.j is the weighting factor. We further approximate the convolutions in Eq. (3) by assuming that the broad PDF component, B.sub.j(x.sub.j), can be represented as constant within the narrow range of x.sub.j that comprises the bulk of the noise PDF (i.e., within one or two standard deviations of G(x.sub.j)), and accordingly we factor B.sub.j(x.sub.j) out of the restricted domain convolution integral. Noting that (1) convolutions of sums are sums of convolutions, (2) the convolution of a Dirac delta function with a Gaussian is a Gaussian, and (3) the convolution of a linear function with a Gaussian is a linear function, an approximate solution to Eq. (3) is
x.sub.j(y)=y.sub.j{w.sub.jB.sub.j(y.sub.j)/[(1w.sub.j)G(y.sub.j)+w.sub.jB.sub.j(y.sub.j)]}(6)

(24) The quantity in braces is the approximate shrinkage function S.sub.j(y.sub.j). For the remainder of this discussion we drop the j subscripts for simplicity.

(25) For illustrative purposes, we follow Simoncelli [1996] and assume a generalized Laplacian distribution,
B(y)exp(|y/s|.sup.P)(7)

(26) Then we write S(y) in terms of the ratio of G to B at y=0, which we denote r. Simplifying the expression, the approximate shrinkage function becomes
S(y)=exp(|y/s|.sup.P)/[r exp(y.sup.2)+exp(|y/s|.sup.P)](8)

(27) Simoncelli [1996] suggests a range of exponents, p, from 0.5 to 1.0. Adler-Golden, S. M., Improved Hyperspectral Anomaly Detection in Heavy-Tailed Backgrounds, IEEE, First Annual WHISPERS Conference, Grenoble, FR (2009) (Adler-Golden [2009]) analyzed PC transforms of HSI data, which are closely related to MNF transforms, and found that the low-noise PC components had typical p values in a similar range. Eq. (8) is a better approximation at smaller p values, and represents the exact solution to Eq. (3) when p=0 (that is, B(y) is a flat distribution up to some maximum absolute y value).

(28) FIGS. 3A-3D compare the noise-free estimation functions x(y) using the approximate Eq. (8) for p=1 with those from exact numerical solutions to Eq. (3). For reference we include the line x=y, which corresponds to no shrinkage applied to the data (i.e., the shrinkage factor is 1). We choose a representative value of s=2, corresponding to a B(y) PDF that is twice as wide as the noise PDF at the 1/e points. We consider a wide range of r values, defined by
r=qexp(|2/s|.sup.P)(9)
with q=1 (FIG. 3A), 10 (FIG. 3B), 100 (FIG. 3C), and 1000 (FIG. 3D). As will be seen below, defining r in this way provides a good means to compare PDFs corresponding to different p values. FIGS. 3A-3D show that the approximate Eq. (8) solutions for x(y) are very similar to the exact results for these p=1 cases. The values from Eq. (8) are slightly larger at large y, meaning that its shrinkage estimates are on the conservative side (closer to 1). The agreement would be better at larger values of s (a broader Laplacian distribution).

(29) FIGS. 3A-3D include the Eq. (8) results for p=0, which are exact. In contrast to the formula for p>0, which contains three unknown parameters (s, p and r), the shrinkage function formula for p=0, which is given by
S(y)=1/[r exp(y.sup.2)+1](10)
has only a single adjustable parameter, r. The excellent agreement with the p=1 results indicates that, with r properly defined, the shrinkage function has little dependence on p, and thus can be estimated well with little or no knowledge of shape of the broad component of the prior PDF. One requires only an estimate of the relative magnitudes of the narrow (delta function) and broad components of the prior PDF, which sets the value of r.

(30) With typical hyperspectral datasets, the vast majority of the spectra are well-modeled with just the first few MNF components. Thus, for all but the lowest j values, the broad PDF component is expected to be small compared to the narrow component. In this situation, r can be estimated by histogramming the distribution of y values in the dataset and taking the ratio of the histogram populations at y=0 and y=2. As an application of this method, FIG. 4 shows some estimated r values obtained from histogram population ratios of MNF components in a thermal infrared radiance data cube, which was measured with an airborne Telops HyperCam hyperspectral sensor. At the highest eigenvector numbers (j values), the broad PDF component becomes extremely small (only a tiny fraction of the spectra have y values above the noise level), and the histogram ratio becomes poorly defined. We choose to err on the conservative side by setting a maximum r value of 1000. For the lowest j values, the broad plus narrow component PDF model may be inaccurate. However, this is of little concern, since here the shrinkage function will be close to 1 regardless of what PDF model is assumed.

(31) Having estimated the shrinkage functions 30 by one of the above procedures, Eq. (2) is used to determine the noise-free MNF component estimates, which comprise the de-noised MNF-transformed spectra 20. Applying the inverse MNF transform, which includes adding back the mean spectrum that was removed in the first step of the forward MNF transform, results in the de-noised spectra 22.

(32) FIGS. 5 and 6 show the application of this signal-to-noise enhancement to the HyperCam radiance data cube. The spectra with maximum values below 700 microflicks are from a patch of reflective material. FIG. 6 shows the same spectra as FIG. 5 obtained after de-noising using the Eq. (9) shrinkage function and the FIG. 4 r values. The noise level is dramatically reduced, especially for the non-reflective surfaces. This is seen in the greatly improved consistency among spectra of the same type of surface, and particularly in the very smooth appearance of the blackbody spectra, which are at a temperature close to that of the atmosphere. At the same time, the de-noising processing preserves the atmospheric spectral structure, which is seen in the spectra of reflective surfaces and warm surfaces. The de-noising processing also preserves the sharp 10.3 micron peak in the brightest spectrum, which appears to be due to an anomaly in the focal plane array response.

(33) Other methods of estimating shrinkage functions 30 from parameters of the MNF components of the spectra may occur to those skilled in the art. These may include: Evaluating the Bayesian signal estimation formula 28 numerically using a model for the prior PDF estimate 26, where the model is a generalized Laplacian distribution, or a model that includes a Dirac delta function, or a model that includes both; Using a user-supervised process to estimate the parameters s, p and q in the generalized Laplacian distribution model of the prior PDF estimate 26; Using a user-supervised process to estimate the parameter r in the shrinkage function equation (10); As a further approximation, using the same shrinkage function 30 for all N of the MNF components.

Second Example

(34) In an alternative embodiment of the method shown in FIG. 7, a stack of image frames 42 from a staring camera 40 are fed into the processing system, comprised of items 44, 46, 48 and 50. The number of pixels in the images must be greater than the number of image frames, denoted N. The method outputs de-noised images 52. The shrinkage function estimator 48 is identical to the shrinkage function estimator 18 shown in FIG. 1.

(35) The description of the alternative embodiment is identical to the description of the preferred embodiment, except that spectrum is replaced with pixel waveform, spectral dataset is replaced with image dataset, and wavelength band is replaced with image frame. Specifically, The first step in the data processing is to subtract the mean waveform (i.e., the set of image means) from the image frames 42. Let .sub.S denote the NN covariance matrix of the image data and .sub.N denote the estimated NN covariance matrix of the noise 44. The matrix of N eigenvectors, V, is the solution to Eq. (1). The projection of the mean-subtracted data onto the N eigenvectors of the V matrix yields a set of N coefficients for each pixel waveform, which we denote y.sub.j for a given eigenvector j. Each of the N MNF components y.sub.j comprises an image; collectively these are the MNF-transformed images 46. Optimal estimates of the true MNF components, x.sub.j, which are uncontaminated with background noise, n.sub.j, are obtained by defining shrinkage functions, denoted S.sub.j(y.sub.j), according to Eqs. (2), (3) and (4). For the small-eigenvalue j's, which are dominated by noise, we may assume that the vast majority of the pixel waveforms have MNF component values x.sub.j that are very close to zero, while a small minority of the pixel waveforms are associated with a relatively broad PDF. We may approximate the PDF of the true MNF coefficients P.sub.j(x.sub.j) as a weighted sum of a broad distribution B(x.sub.j) and a Dirac delta function D(x.sub.j), the latter representing the limit of an extremely narrow distribution. Then the same logic found in the preferred embodiment leads to Eqs. (5) through (10). Other methods of estimating shrinkage functions 48 from parameters of the MNF components may occur to those skilled in the art, and may include all of the methods described in the preferred embodiment. Having estimated the shrinkage functions 48, Eq. (2) is used to determine the noise-free MNF component estimates, which comprise the de-noised MNF-transformed data 50. Applying the inverse MNF transform, which includes adding back the mean waveform that was removed in the first step of the forward MNF transform, results in the de-noised images 52.

(36) Application of this method to a stack of image frames is most beneficial when the bulk of the pixel waveforms can be represented with a small number of independent components, i.e., a small number of dimensions. This situation may apply when there is time variation in only a small number of pixels, and/or when image-wide motion, such as pointing jitter, is of low amplitude.

(37) FIGS. 8 and 9 show the application of this method to a video from a hand-held camera. The original video is described in a paper Adler-Golden, S. M., S. C. Richtsmeier, and R. M. Shroll, Suppression of Subpixel Sensor Jitter Fluctuations using Temporal Whitening, Proc. SPIE 6969, Signal and Data Processing of Small Targets 2008, 69691D (Apr. 16, 2008); doi: 10.1117/12.780005 ([Adler-Golden, 2008]), and consists of 30 frames taken with a consumer digital camera. FIG. 8 shows a frame of the video after Gaussian white noise equivalent to 15.7% of full scale was added to the video. FIG. 9 shows the same frame after the video with noise added was processed with the present method; the noise has been dramatically reduced, by between typically one-half to one order of magnitude. The method is far superior to simply averaging the frames together, which would suffer from reduced sharpness due to the motion resulting from the camera being hand-held.

(38) A number of implementations have been described. Nevertheless, it will be understood that additional modifications may be made without departing from the scope of the inventive concepts described herein, and, accordingly, other embodiments are within the scope of the following claims.