Signaltonoise enhancement

0Associated
Cases 
0Associated
Defendants 
0Accused
Products 
0Forward
Citations 
0
Petitions 
1
Assignment
First Claim
1. A signaltonoise 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 noisefree MNF components of the given spectrum, in which the nonlinear mathematical functions are given by S(y)=exp(−
y/s^{p})/[r exp(−
y^{2})+exp(−
y/s^{p})];
applying the inverse MNF transform to convert the estimated noisefree 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 noisefree MNF components of the multiplicity of spectra.
1 Assignment
0 Petitions
Accused Products
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 signaltonoise 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.
3 Citations
No References
Efficient near neighbor search (ENNsearch) method for high dimensional data sets with noise  
Patent #
US 20030187616A1
Filed 03/29/2002

Current Assignee
The United States of America As Represented By The Secretary of Agriculture

Sponsoring Entity
The United States of America As Represented By The Secretary of Agriculture

Method and system for increasing signaltonoise ratio  
Patent #
US 8,358,866 B2
Filed 01/29/2007

Current Assignee
Canadian Space Agency

Sponsoring Entity
Canadian Space Agency

PORTABLE HYPERSPECTRAL IMAGER  
Patent #
US 20140354868A1
Filed 04/22/2014

Current Assignee
Corning Incorporated

Sponsoring Entity
Corning Incorporated

8 Claims
 1. A signaltonoise 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 noisefree MNF components of the given spectrum, in which the nonlinear mathematical functions are given by S(y)=exp(−
y/s^{p})/[r exp(−
y^{2})+exp(−
y/s^{p})];applying the inverse MNF transform to convert the estimated noisefree 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 noisefree MNF components of the multiplicity of spectra.  View Dependent Claims (5)
 2. A signaltonoise 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 noisefree MNF components of the given spectrum, in which the nonlinear mathematical functions are given by S(y)=1/[r exp(−
y^{2})+1];applying the inverse MNF transform to convert the estimated noisefree 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 noisefree MNF components of the multiplicity of spectra.  View Dependent Claims (6)
 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 noisefree MNF components of the given spectrum, in which the nonlinear mathematical functions are given by S(y)=exp(−
y/s^{p})/[r exp(−
y^{2})+exp(−
y/s^{p})];applies the inverse MNF transform to convert the estimated noisefree 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 noisefree MNF components of the multiplicity of spectra.  View Dependent Claims (7)
 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 noisefree MNF components of the given spectrum, in which the nonlinear mathematical functions are given by S(y)=1/[r exp(−
y^{2})+1];applies the inverse MNF transform to convert the estimated noisefree 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 noisefree MNF components of the multiplicity of spectra.  View Dependent Claims (8)
1 Specification
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 signaltonoise 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.
In hyperspectral imaging (HSI) sensors, images are collected in which each spatial pixel contains a spectrum—that is, a set of measured light intensities in a large number of spectral bands, typically tens to hundreds, each defined by a wavelength range and response function. Electronic noise, which may arise from the sensor itself or from statistical fluctuations in the background radiation, adds a corrupting component to the measurements, interfering with the data analysis. In the infrared (IR) region, optical sensors are especially susceptible to background noise. Measures to reduce noise in IR sensors, such as aggressive cooling and selection of only the best detectors, can make the superior IR sensors extremely expensive. Consequently, there is a need for data postprocessing algorithms that can mitigate the effects of noise in the spectral data and provide useful data products from less expensive sensors. Certain of these algorithms might also be usefully applied to spectral datasets obtained with a nonimaging sensor that makes a large number of sequential spectral measurements, such as a laser Raman or fluorescence spectrometer, where there may be considerable noise from the scattered light background.
A number of methods for spectral imagery that reduce sensor noise, referred to as denoising algorithms, have been described in the literature, and tend to fall into two major categories. Methods in the first category are known as spectral noise reduction methods since they work in the wavelength (spectral) coordinate and provide additional applicability to nonimage datasets. Methods in the second category work in image spatial coordinates, and process one wavelength image at a time. The two methods can be combined, as in U.S. Pat. No. 8,358,866.
Qian, S., J. Lévesque and R. Neville, “Evaluation of Noise Removal of Radiance Data on Onboard Data Compression of Hyperspectral Imagery,” 2005 WSEAS Int. Conf. on Remote Sensing, Venice, Italy, pp. 3742, Nov. 24, 2005 (“Qian [2005]”) describes a spectral noise reduction method for HSI data that separates the pixel spectra into distinct classes that differ by an amount greater than the noise level. The spectra in each class are then replaced by the class averages, thereby reducing the noise. A drawback of this method is that the number of required classes increases as the diversity of spectra increases and as the noise level decreases, resulting in increasing and unpredictable computation time. Another drawback is that many spectra may wind up in classes consisting only of themselves, even though they are similar to other spectra; those spectra do not undergo any noise reduction. Another drawback is that the transformation from input to output spectra is discontinuous, potentially leading to step artifacts in the output.
An alternative, more common spectral noise reduction method is to perform Maximum Noise Fraction (MNF) truncation, as first described by Green, A. A., M. Berman, P. Switzer, and M. D. Craig, “A Transformation for Ordering Multispectral Data in Terms of Image Quality with Implications for Noise Removal,” IEEE Transactions on Geoscience and Remote Sensing, 26(1):6574, 1988 (“Green [1988]”); a recent adaptation for HSI data processing is found in Bjorgan, A., and Randeberg, L L., “RealTime Noise Removal for LineScanning Hyperspectral Devices Using a Maximum Noise FractionBased Approach,” Sensors (Basel, Switzerland), 15(2): 33623378, 2015, doi:10.3390/s150203362 (“Bjorgan [2015]”). The MNF operation transforms the data to an orthonormal coordinate system in which most of the noise is isolated from most of the signal. In this method, a variance or covariance linear operator describing the sensor noise is estimated from either the data itself or from a sensor specification. The data are then multiplied by the inverse square root of the operator, which “whitens” the noise, and the result is then processed with a Principal Component (PC) transformation. The PC transformation projects the spectra onto the eigenvectors of the spectral covariance matrix of the noisewhitened data. Each eigenvector has a corresponding eigenvalue. The combined noise whitening step and PC transformation can be considered as a single linear operation comprising the MNF transform, as described in Green [1988], Phillips, R. D., L. T. Watson, C. E. Blinn and R. H. Wynne, “An adaptive noise reduction technique for improving the utility of hyperspectral data,” Pecora 17, The Future of Land Imaging . . . Going Operational, Denver, Colo., Nov. 1820, 2008 (“Phillips [2008]”) and Bjorgan [2015]. The MNF components (i.e., the projection coefficients) that correspond to very small eigenvalues, comparable to the noise variance, are dominated by noise. To suppress the noise in the data, these MNF components are set to zero. The data are then processed with the inverse MNF transform to restore the data to their original units. A serious drawback of this method is that this processing truncates the data dimensionality. This results in distortion of “rare” or anomalous spectra in the dataset, and makes it difficult, if not impossible, to use the processed data for a rare target detection application.
An improvement on the original MNF method for hyperspectral image denoising is described by Phillips [2008]. Instead of zeroing the noisy smalleigenvalue MNF components they are processed with spatial median filtering. This reduces the noise and retains the major spatial structure in those components, but their spatial resolution is reduced. As a result, rare or unusual spectra in the dataset that occupy a small area, such a few pixels or less, are distorted.
The second general class of noise reduction methods involves performing a wavelet transform on each wavelength band image. Example methods are described in Sihag, R., R. Sharma and V. Setia, “Wavelet Thresholding for Image Denoising,” ICVCI 2011, Proc. Int. J. Computer Applic. 2011 (“Sihag [2011]”) and Simoncelli, E. P. and Adelson, “Noise Removal via Bayesian Wavelet Coring,” Proc. 3^{rd }IEEE Intl. Conf. on Image Processing, Vol I, pp. 379382, Lausanne, 1619 Sep. 1996 (“Simoncelli [1996]”). The sensor noise is preferentially found in the wavelet decomposition images that correspond to fine spatial structure. These decomposition images are processed with a threshold function. Signal absolute values below the threshold, which consist predominantly of noise, are set to zero, while signal absolute values above the threshold, which consist predominantly of true signal, are either left unchanged (with a “hard” threshold) or reduced (with a “soft” threshold). The soft threshold function is also known as a shrinkage function [Sihag, 2011]. The thresholded wavelet decomposition images are then processed with the inverse wavelet transform to yield a denoised version of the original wavelength band image. A drawback of this method is that it has difficulty distinguishing sensor noise from surface texture, and thus it may oversmooth the image and degrade its sharpness. This will in turn distort the spectra. Also, wavelet methods are applicable only to image data, not individual spectra.
U.S. Pat. No. 8,358,866 for denoising HSI data combines a wavelet transform method with a spectral filtering method. The method is also described in Othman, H. and Qian, S., “Noise reduction of hyperspectral imagery using hybrid spatialspectral derivativedomain wavelet shrinkage”, IEEE Trans. Geosci. Remote Sensing, Vol. 44 (2), pp. 397408, 2006 (“Othman [2006]”). First, wavelet denoising is performed. The output is further processed by a spectral lowpass filtering step, involving taking differences of adjacent wavelengths, smoothing the result with a running lowpass filter, and then reversing the differencing process to restore the spectra, but with noise further reduced. This method has the advantages and disadvantage of other wavelet transform methods, but is capable of additional noise reduction via the spectral filtering step. On the other hand, the spectral filtering process can distort data that contain true spectral features that are comparable in resolution to the spectral resolution of the sensor. As an example of such features, remote sensing data acquired in the long wavelength infrared (LWIR) wavelength region contains fine spectral structure from atmospheric water vapor spectral lines. Since this structure is as narrow as, or narrower than, the spacing between adjacent wavelength bands of typical LWIR hyperspectral sensors, this structure may be lost in the lowpass filtering step.
Certain types of hyperspectral imaging sensors collect the data in a staring mode, in which images of a fixed view are sequentially collected at different wavelengths by tuning a wavelengthselective optical element. In this situation, the images differ in their acquisition time, and comprise staring video imagery. Therefore, it will occur to those skilled in the art that many of the HSI denoising algorithms may also be applied to ordinary staring video imagery, in which the wavelength is not being varied. These data, in which the time dimension replaces the wavelength dimension, comprise images of waveforms rather than images of spectra.
The signaltonoise enhancement, or denoising, of the present disclosure is aimed at improving the utility of hyperspectral imagery or other large spectral datasets without introducing significant spectral distortion in or spectral smoothing of the data, even for rare or anomalous spectra. In particular, this disclosure is aimed at enhancing the imagery from affordable infrared hyperspectral sensors, such as current commercial sensors operating in the midinfrared to longwavelength infrared. Since this method provides the best noise reduction for abundant spectra in the dataset, the method should especially benefit terrain classification applications, such as in minerology, agriculture, and land use. Another object of this disclosure is to provide a denoising method for staring video imagery, which has the equivalent data structure to hyperspectral imagery. Another objective is to provide a method that is fully automated. Another objective is to provide a moderate, and predictable, computing time.
In one aspect, a signaltonoise enhancement method includes 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 noisefree MNF components of the given spectrum, applying the inverse MNF transform to convert the estimated noisefree 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 noisefree MNF components of the multiplicity of spectra.
Embodiments may include one of the above and/or below features, or any combination thereof. The nonlinear mathematical functions may be given by: S(y)=exp(−y/s^{P})/[r exp(−y^{2})+exp(−y/s^{P})]. The nonlinear mathematical functions may be given by: S(y)=1/[r exp(−y^{2})+1]. The multiplicity of spectra may comprise a hyperspectral image.
In another aspect, a system includes 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 noisefree MNF components of the given spectrum, applies the inverse MNF transform to convert the estimated noisefree 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 noisefree MNF components of the multiplicity of spectra.
Embodiments may include one of the above and/or below features, or any combination thereof. The nonlinear mathematical functions may be given by: S(y)=exp(−y/s^{P})/[r exp(−y^{2})+exp(−y/s^{P})]. The nonlinear mathematical functions may be given by: S(y)=1/[r exp(−y^{2})+1]. The spectrometer may comprise a hyperspectral imaging sensor.
In another aspect, a signaltonoise enhancement method includes receiving from an imaging device a time series of images that comprise a multiplicity of image frames, in which the number of pixels exceeds the number of image frames, estimating the covariance matrix of the background noise associated with the images, transforming the multiplicity of images to Maximum Noise Fraction (MNF) components, each component comprising the coefficients resulting from projection of the pixel waveforms 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 pixel waveform into estimated values of the corresponding noisefree MNF components of the given pixel waveform, and applying the inverse MNF transform to convert the estimated noisefree MNF components of the multiplicity of pixel waveforms back into the original units of the images. An output image dataset is generated from the converted estimated noisefree MNF components of the multiplicity of pixel waveforms.
Embodiments may include one of the above and/or below features, or any combination thereof. The nonlinear mathematical functions may be given by: S(y)=exp(−y/s^{P})/[r exp(−y^{2})+exp(−y/s^{P})]. The nonlinear mathematical functions may be given by: S(y)=1/[r exp(−y^{2})+1]. The covariance matrix of the background noise may comprise a diagonal constant matrix.
In another aspect, a system includes a camera that provides a time series of images that comprise a multiplicity of image frames, in which the number of pixels exceeds the number of image frames, and digital image processing circuitry which, in operation: estimates the covariance matrix of the background noise associated with the images, transforms the multiplicity of images to Maximum Noise Fraction (MNF) components, each component comprising the coefficients resulting from projection of the pixel waveforms 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 pixel waveform into estimated values of the corresponding noisefree MNF components of the given pixel waveform, applies the inverse MNF transform to convert the estimated noisefree MNF components of the multiplicity of pixel waveforms back into the original units of the images, and generates an output image dataset from the converted estimated noisefree MNF components of the multiplicity of pixel waveforms.
Embodiments may include one of the above and/or below features, or any combination thereof. The nonlinear mathematical functions may be given by: S(y)=exp(−y/s^{P})/[r exp(−y^{2})+exp(−y/s^{P})]. The nonlinear mathematical functions may be given by S(y)=1/[r exp(−y^{2})+1]. The covariance matrix of the background noise may comprise a diagonal constant matrix.
Other objects, features and advantages will occur to those skilled in the art from the following detailed description, and the accompanying drawings, in which:
The signaltonoise 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.
In a first example of a signaltonoise enhancement shown in
An estimate of the N×N 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.
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 noisewhitening 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 MNFtransformed data consist of the projections of the mean spectrumsubtracted data onto the eigenvectors of the PC transform. We assume without loss of generality that the eigenvectors have been unitnormalized.
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 N×N covariance matrix of the spectral data and ΣN denote the N×N covariance matrix of the noise. The matrix of N eigenvectors, V, is the solution to the equation
Σ_{S}Σ_{N}^{−1}=VΛV^{−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.
The projection of the meansubtracted data onto the N eigenvectors of the V matrix yields a set of N coefficients for each spectrum, which we denote y_{j }for a given eigenvector j. We shall refer to the y_{j }as the j^{th }MNF component. If the original data comprise a hyperspectral image (data cube), then each of the N MNF components y_{j }comprises an image. Collectively, the y_{j }comprise the MNFtransformed spectral dataset 16.
Next, we seek optimal estimates of the “true” MNF components, x_{j}, which are uncontaminated with background noise, n_{j}, based on the observed y_{j}=x_{j}+n_{j}. This is done by defining a shrinkage function, denoted S_{j}(y_{j}), which is a multiplicative factor that converts the observed noisecontaining MNF component y_{j }into an estimate of the noisefree MNF component. That is, the optimal estimate x′_{j }is given by
x′_{j}(y_{j})=S_{j}(y_{j})y_{j} (2)
The preferred embodiment of the shrinkage function estimator procedure 18 is shown in
The optimal estimate of a noisefree MNF component value, denoted x′_{j}, is its mean value in the posterior PDF, which includes the added noise. This mean value, which provides an unbiasedleast squares estimate of x_{j }given y_{j}, is computed from the Bayesian signal estimation formula 28, found in Simoncelli [1996]. We assume that the noise PDF is the unitvariance Gaussian function, G, and define the function Q_{j}(x_{j})=x_{j}P_{j}(x_{j}). Then the Bayesian signal estimation formula may be compactly written as a quotient of convolutions,
x′_{j}(y_{j})=(G*Q_{j})(y_{j})/(G*P_{j})(y_{j}) (3)
Using Eq. (2), the shrinkage function 30 is then given by
S_{j}(y_{j})=x′_{j}(y_{j})/y_{j}. (4)
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 denoising results. Some examples are described below.
For the smalleigenvalue j'"'"'s, which are dominated by noise, it is assumed that the vast majority of the spectra have MNF component values x_{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_{j}(x_{j}) as a weighted sum of a broad distribution B(x_{j}) and a Dirac delta function D(x_{j}), the latter representing the limit of an extremely narrow distribution:
P_{j}(x_{j})=w_{j}B_{j}(x_{j})+(1−w_{j})D(x_{j}) (5)
where w_{j }is the weighting factor. We further approximate the convolutions in Eq. (3) by assuming that the broad PDF component, B_{j}(x_{j}), can be represented as constant within the narrow range of x_{j }that comprises the bulk of the noise PDF (i.e., within one or two standard deviations of G(x_{j})), and accordingly we factor B_{j}(x_{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′_{j}(y)=y_{j}{w_{j}B_{j}(y_{j})/[(1−w_{j})G(y_{j})+w_{j}B_{j}(y_{j})]} (6)
The quantity in braces is the approximate shrinkage function S_{j}(y_{j}). For the remainder of this discussion we drop the j subscripts for simplicity.
For illustrative purposes, we follow Simoncelli [1996] and assume a generalized Laplacian distribution,
B(y)∝exp(−y/s^{P}) (7)
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^{P})/[r exp(−y^{2})+exp(−y/s^{P})] (8)
Simoncelli [1996] suggests a range of exponents, p, from 0.5 to 1.0. AdlerGolden, S. M., “Improved Hyperspectral Anomaly Detection in HeavyTailed Backgrounds,” IEEE, First Annual WHISPERS Conference, Grenoble, FR (2009) (“AdlerGolden [2009]”) analyzed PC transforms of HSI data, which are closely related to MNF transforms, and found that the lownoise 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).
r=qexp(−2/s^{P}) (9)
with q=1 (
S(y)=1/[r exp(−y^{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.
With typical hyperspectral datasets, the vast majority of the spectra are wellmodeled 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,
Having estimated the shrinkage functions 30 by one of the above procedures, Eq. (2) is used to determine the noisefree MNF component estimates, which comprise the denoised MNFtransformed 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 denoised spectra 22.
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 usersupervised process to estimate the parameters s, p and q in the generalized Laplacian distribution model of the prior PDF estimate 26;
 Using a usersupervised 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.
In an alternative embodiment of the method shown in
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 Σ_{S }denote the N×N covariance matrix of the image data and Σ_{N }denote the estimated N×N covariance matrix of the noise 44.
 The matrix of N eigenvectors, V, is the solution to Eq. (1).
 The projection of the meansubtracted data onto the N eigenvectors of the V matrix yields a set of N coefficients for each pixel waveform, which we denote y_{j }for a given eigenvector j. Each of the N MNF components y_{j }comprises an image; collectively these are the MNFtransformed images 46.
 Optimal estimates of the “true” MNF components, x_{j}, which are uncontaminated with background noise, n_{j}, are obtained by defining shrinkage functions, denoted S_{j}(y_{j}), according to Eqs. (2), (3) and (4).
 For the smalleigenvalue j'"'"'s, which are dominated by noise, we may assume that the vast majority of the pixel waveforms have MNF component values x_{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_{j}(x_{j}) as a weighted sum of a broad distribution B(x_{j}) and a Dirac delta function D(x_{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 noisefree MNF component estimates, which comprise the denoised MNFtransformed 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 denoised images 52.
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 imagewide motion, such as pointing jitter, is of low amplitude.
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.