Source code for pysptools.noise.dnoise_int

#
#------------------------------------------------------------------------------
# Copyright (c) 2013-2014, Christian Therien
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#    http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
#------------------------------------------------------------------------------
#
# dnoise_int.py - This file is part of the PySptools package.
#

from __future__ import division

import os.path as osp
import numpy as np
from scipy.signal import fftconvolve
from math import factorial
from . import dnoise
from .inval import *


[docs]class SavitzkyGolay(object): """Apply a Savitzky Golay low pass filter.""" def __init__(self): self.denoised = None self.dbands = None @DenoiseSpectraInputValidation('SavitzkyGolay') def denoise_spectra(self, M, window_size, order, deriv=0, rate=1): """ Apply the Savitzky Golay filter on each spectrum. Smooth (and optionally differentiate) data with a Savitzky-Golay filter. The Savitzky-Golay filter removes high frequency noise from data. It has the advantage of preserving the original shape and features of the signal better than other types of filtering approaches, such as moving averages techniques. Parameters: M: `numpy array` A HSI cube (m x n x p). window_size: `int` the length of the window. Must be an odd integer number. order: `int` the order of the polynomial used in the filtering. Must be less then `window_size` - 1. deriv: `int [default 0]` the order of the derivative to compute (default = 0 means only smoothing). Returns: `numpy array` the smoothed signal (or it's n-th derivative) (m x n x p). Code source: The scipy Cookbook, SavitzkyGolay section. This class is not under the copyright of this file. """ h, w, numBands = M.shape M = np.reshape(M, (w*h, numBands)) self.denoised = self._denoise1d(M, window_size, order, deriv, rate) self.denoised = np.reshape(self.denoised, (h, w, numBands)) return self.denoised @DenoiseBandsInputValidation('SavitzkyGolay') def denoise_bands(self, M, window_size, order, derivative=None): """ Apply the Savitzky Golay filter on each band. Parameters: M: `numpy array` A HSI cube (m x n x p). window_size: `int` the length of the window. Must be an odd integer number. order: `int` the order of the polynomial used in the filtering. Must be less then `window_size` - 1. derivative: `string [default None]` direction of the derivative to compute, can be None, 'col', 'row', 'both'. Returns: `numpy array` the smoothed signal (or it's n-th derivative) (m x n x p). Code source: The scipy Cookbook, SavitzkyGolay section. This class is not under the copyright of this file. """ h, w, numBands = M.shape self.dbands = np.ones((h, w, numBands), dtype=np.float) for i in range(numBands): self.dbands[:,:,i] = self._denoise2d(M[:,:,i], window_size, order, derivative) return self.dbands def _denoise1d(self, M, window_size, order, deriv, rate): try: window_size = np.abs(np.int(window_size)) order = np.abs(np.int(order)) except ValueError as msg: raise ValueError("in SavitzkyGolay.denoise_spectra(), window_size and order have to be of type int") if window_size % 2 != 1 or window_size < 1: raise TypeError("in SavitzkyGolay.denoise_spectra(), window_size size must be a positive odd number") if window_size < order + 2: raise TypeError("in SavitzkyGolay.denoise_spectra(), window_size is too small for the polynomials order") order_range = range(order+1) half_window = (window_size -1) // 2 # precompute coefficients b = np.mat([[k**i for i in order_range] for k in range(-half_window, half_window+1)]) m = np.linalg.pinv(b).A[deriv] * rate**deriv * factorial(deriv) # pad the signal at the extremes with # values taken from the signal itself N, p = M.shape dn = np.ones((N,p), dtype=np.float) long_signal = np.ndarray(p+2, dtype=np.float) for i in range(N): y = M[i] firstvals = y[0] - np.abs( y[1:half_window+1][::-1] - y[0] ) lastvals = y[-1] + np.abs(y[-half_window-1:-1][::-1] - y[-1]) long_signal = np.concatenate((firstvals, y, lastvals)) dn[i] = fftconvolve(long_signal, m, mode='valid') return dn def _denoise2d(self, z, window_size, order, derivative=None): # number of terms in the polynomial expression n_terms = ( order + 1 ) * ( order + 2) / 2.0 if window_size % 2 == 0: raise ValueError('in SavitzkyGolay.denoise_bands(), window_size must be odd') if window_size**2 < n_terms: raise ValueError('in SavitzkyGolay.denoise_bands(), order is too high for the window size') half_size = window_size // 2 # exponents of the polynomial. # p(x,y) = a0 + a1*x + a2*y + a3*x^2 + a4*y^2 + a5*x*y + ... # this line gives a list of two item tuple. Each tuple contains # the exponents of the k-th term. First element of tuple is for x # second element for y. # Ex. exps = [(0,0), (1,0), (0,1), (2,0), (1,1), (0,2), ...] exps = [ (k-n, n) for k in range(order+1) for n in range(k+1) ] # coordinates of points ind = np.arange(-half_size, half_size+1, dtype=np.float64) dx = np.repeat( ind, window_size ) dy = np.tile( ind, [window_size, 1]).reshape(window_size**2, ) # build matrix of system of equation A = np.empty( (window_size**2, len(exps)) ) for i, exp in enumerate( exps ): A[:,i] = (dx**exp[0]) * (dy**exp[1]) # pad input array with appropriate values at the four borders new_shape = z.shape[0] + 2*half_size, z.shape[1] + 2*half_size Z = np.zeros( (new_shape) ) # top band band = z[0, :] Z[:half_size, half_size:-half_size] = band - np.abs( np.flipud( z[1:half_size+1, :] ) - band ) # bottom band band = z[-1, :] Z[-half_size:, half_size:-half_size] = band + np.abs( np.flipud( z[-half_size-1:-1, :] ) -band ) # left band band = np.tile( z[:,0].reshape(-1,1), [1,half_size]) Z[half_size:-half_size, :half_size] = band - np.abs( np.fliplr( z[:, 1:half_size+1] ) - band ) # right band band = np.tile( z[:,-1].reshape(-1,1), [1,half_size] ) Z[half_size:-half_size, -half_size:] = band + np.abs( np.fliplr( z[:, -half_size-1:-1] ) - band ) # central band Z[half_size:-half_size, half_size:-half_size] = z # top left corner band = z[0,0] Z[:half_size,:half_size] = band - np.abs( np.flipud(np.fliplr(z[1:half_size+1,1:half_size+1]) ) - band ) # bottom right corner band = z[-1,-1] Z[-half_size:,-half_size:] = band + np.abs( np.flipud(np.fliplr(z[-half_size-1:-1,-half_size-1:-1]) ) - band ) # top right corner band = Z[half_size,-half_size:] Z[:half_size,-half_size:] = band - np.abs( np.flipud(Z[half_size+1:2*half_size+1,-half_size:]) - band ) # bottom left corner band = Z[-half_size:,half_size].reshape(-1,1) Z[-half_size:,:half_size] = band - np.abs( np.fliplr(Z[-half_size:, half_size+1:2*half_size+1]) - band ) # solve system and convolve if derivative == None: m = np.linalg.pinv(A)[0].reshape((window_size, -1)) return fftconvolve(Z, m, mode='valid') elif derivative == 'col': c = np.linalg.pinv(A)[1].reshape((window_size, -1)) return fftconvolve(Z, -c, mode='valid') elif derivative == 'row': r = np.linalg.pinv(A)[2].reshape((window_size, -1)) return fftconvolve(Z, -r, mode='valid') elif derivative == 'both': c = np.linalg.pinv(A)[1].reshape((window_size, -1)) r = np.linalg.pinv(A)[2].reshape((window_size, -1)) return fftconvolve(Z, -r, mode='valid'), scipy.signal.fftconvolve(Z, -c, mode='valid')
[docs] def plot_spectrum_sample(self, M, path, x, y, suffix=None): """ Plot a spectrum sample with the original and the filtered signal. Parameters: M: 'numpy array' The original cube (m x n x p). path: `string` The path where to put the plot. x: `int` The x coordinate. y: `int` The y coordinate. suffix: `string [default None]` Add a suffix to the file name. """ import matplotlib.pyplot as plt plt.ioff() plt.plot(M[x,y]) plt.plot(self.denoised[x,y], color='r') plt.xlabel('Wavelength') plt.ylabel('Brightness') if suffix == None: fout = osp.join(path, 'SavitzkyGolay_x{0}_y{1}.png'.format(x,y)) else: fout = osp.join(path, 'SavitzkyGolay_x{0}_y{1}_{2}.png'.format(x,y,suffix)) plt.savefig(fout) plt.close()
[docs] def plot_bands_sample(self, path, band_no, suffix=None): """ Plot a filtered band. Parameters: path: `string` The path where to put the plot. band_no: `int or string` The band index. If band_no == 'all', plot all the bands. suffix: `string [default None]` Add a suffix to the file name. """ import matplotlib.pyplot as plt plt.ioff() if band_no == 'all': for i in range(self.dbands.shape[2]): plt.imshow(self.dbands[:,:,i], interpolation='none') if suffix == None: fout = osp.join(path, 'SavitzkyGolay_band_{0}.png'.format(i)) else: fout = osp.join(path, 'SavitzkyGolay_band_{0}_{1}.png'.format(i, suffix)) plt.savefig(fout) plt.close() else: plt.imshow(self.dbands[:,:,band_no], interpolation='none') if suffix == None: fout = osp.join(path, 'SavitzkyGolay_band_{0}.png'.format(band_no)) else: fout = osp.join(path, 'SavitzkyGolay_band_{0}_{1}.png'.format(band_no, suffix)) plt.savefig(fout) plt.close()
[docs] def display_spectrum_sample(self, M, x, y, suffix=None): """ Display a spectrum sample with the original and the filtered signal. Parameters: M: 'numpy array' The original cube (m x n x p). x: `int` The x coordinate. y: `int` The y coordinate. suffix: `string [default None]` Add a suffix to the title. """ import matplotlib.pyplot as plt plt.plot(M[x,y]) plt.plot(self.denoised[x,y], color='r') plt.xlabel('Wavelength') plt.ylabel('Brightness') if suffix == None: plt.title('SG spectrum sample, x={0}, y={1}'.format(x,y)) else: plt.title('SG spectrum sample, x={0}, y={1} - {2}'.format(x,y,suffix)) plt.show() plt.close()
[docs] def display_bands_sample(self, band_no, suffix=None): """ Display a filtered band. Parameters: band_no: `int or string` The band index. If band_no == 'all', plot all the bands. suffix: `string [default None]` Add a suffix to the title. """ import matplotlib.pyplot as plt if band_no == 'all': for i in range(self.dbands.shape[2]): plt.imshow(self.dbands[:,:,i], interpolation='none') if suffix == None: plt.title('SG band {0}'.format(i)) else: plt.title('SG band {0} - {1}'.format(i, suffix)) plt.show() plt.close() else: plt.imshow(self.dbands[:,:,band_no], interpolation='none') if suffix == None: plt.title('SG band {0}'.format(band_no)) else: plt.title('SG band {0} - {1}'.format(band_no, suffix)) plt.show() plt.close()
[docs]class Whiten(object): """Whiten the cube.""" def __init__(self): self.dM = None @ApplyInputValidation('Whiten') def apply(self, M): """ Whitens a HSI cube. Use the noise covariance matrix to decorrelate and rescale the noise in the data (noise whitening). Results in transformed data in which the noise has unit variance and no band-to-band correlations. Parameters: M: `numpy array` A HSI cube (m x n x p). Returns: `numpy array` A whitened HSI cube (m x n x p). """ h, w, numBands = M.shape M = np.reshape(M, (w*h, numBands)) dM = dnoise.whiten(M) self.dM = np.reshape(dM, (h, w, numBands)) return self.dM
[docs] def get(self): """ Returns: `numpy array` The whitened HSI cube (m x n x p). """ return self.dM
[docs]class MNF(object): """Transform a HSI cube.""" def __init__(self): self.mnf = None self.transform = None self.wdata = None # temp @ApplyInputValidation('MNF') def apply(self, M): """ A linear transformation that consists of a noise whitening step and one PCA rotation. This process is designed to * determine the inherent dimensionality of image data, * segregate noise in the data, * allow efficient elimination and/or reduction of noise, and * reduce the computational requirements for subsequent processing. Parameters: M: `numpy array` A HSI cube (m x n x p). Returns: `numpy array` A MNF transformed cube (m x n x p). References: C-I Change and Q Du, "Interference and Noise-Adjusted Principal Components Analysis," IEEE TGRS, Vol 36, No 5, September 1999. """ from sklearn.decomposition import PCA w = Whiten() wdata = w.apply(M) self.wdata = wdata #temp h, w, numBands = wdata.shape X = np.reshape(wdata, (w*h, numBands)) self.transform = PCA() mnf = self.transform.fit_transform(X) self.mnf = np.reshape(mnf, (h, w, numBands)) return self.mnf @XInputValidation('MNF') def inverse_transform(self, X): """ Inverse the PCA rotation step. The cube stay whitened. Usefull if you want to denoise noisy bands before the rotation. X: `numpy array` A transformed (MNF) cube (m x n x p). Return: `numpy array` A inverted cube (m x n x p). """ h, w, numBands = X.shape X = np.reshape(X, (w*h, numBands)) M = self.transform.inverse_transform(X) M = np.reshape(M, (h, w, numBands)) return M
[docs] def get_components(self, n): """ Return: `numpy array` Return the first n bands (maximum variance bands). """ return self.mnf[:,:,:n]
[docs] def plot_components(self, path, n_first=None, colorMap='jet', suffix=None): """ Plot some bands. Parameters: path: `string` The path where to put the plot. n_first: `int [default None]` Print the first n components. colorMap: `string [default jet]` A matplotlib color map. suffix: `string [default None]` Suffix to add to the title. """ import matplotlib.pyplot as plt if n_first != None: n = min(n_first, self.mnf.shape[2]) else: n = self.mnf.shape[2] plt.ioff() for i in range(n): plt.imshow(self.mnf[:,:,i], interpolation='none', cmap=colorMap) if suffix == None: fout = osp.join(path, 'MNF_bandno_{0}.png'.format(i+1)) else: fout = osp.join(path, 'MNF_bandno_{0}_{1}.png'.format(i+1, suffix)) plt.savefig(fout) plt.clf() plt.close()
[docs] def display_components(self, n_first=None, colorMap='jet', suffix=None): """ Display some bands. Parameters: n_first: `int [default None]` Display the first n components. colorMap: `string [default jet]` A matplotlib color map. suffix: `string [default None]` Suffix to add to the title. """ import matplotlib.pyplot as plt if n_first != None: n = min(n_first, self.mnf.shape[2]) else: n = self.mnf.shape[2] for i in range(n): plt.imshow(self.mnf[:,:,i], interpolation='none', cmap=colorMap) if suffix == None: plt.title('MNF bandno {0}'.format(i+1)) else: plt.title('MNF bandno {0} - {1}.png'.format(i+1, suffix)) plt.show() plt.clf() plt.close()