HSTI 0.0.36

The NEWTEC HSTI package contains fundamental functions for the data analysis of hyperspectral thermal images (HSTI).

This package contains functions used in data processing of hyperspectral images captured using a scanning Fabry-Pérot interferometer (FPI). This includes transmission simulations of the FPI itself.

---

Image handling

HSTI.import_data_cube(path, verbose = False)

This function imports the hyperspectral thermal datacube from the raw output of the camera. The path that the function uses as input must be the one containing the 'images' directory.

HSTI.import_image_acquisition_settings(path, verbose = False):

This function imports the image acquisition settings during the capturing event. The path that the function uses as input must be the one containing the 'images' directory.

HSTI.export_data_cube(data_cube, folder_name)

This function takes an HSTI numpy array and exports it as individual .ppm images to a folder given by folder_name.

HSTI.save_img(img, path)

Save an individual numpy 2D array as ppm image.

HSTI.import_mirror_sep_from_path(path, number_of_bands)

This function generates a mirror separation axis based on the .npy polynomial calibration file supplied by path. A number of bands must also be supplied matching the number of spectral bands in the datacube.

HSTI.import_camera_response(path)

This function imports the response of the thermal camera (.pkl file) as a scipy interpolate object. A 1D numpy array containing the response can then be generated by supplying an appropriate wavelength axis.

---

Preprocessing

HSTI.median_filter_cube(data_cube, kernel_size)

This function runs a median filter across the image plane. The size of the kernel must be defined.

HSTI.mean_center(data_cube, axis = 's')

This function subtracts the mean from the data, either the mean of each spectrum (axis = 's') or the mean of each band (axis = 'b').

HSTI.autoscale(data_cube, axis = 's')

This function subtracts the mean and scales with STD. Setting axis = 's' is the same as doing SNV (standard normal variate).

HSTI.norm_normalization(data_cube, order, axis = 's')

This function uses the norm of a given order for normalization. If axis = 's', then each spectrum is divided by its norm. If axis = 'b', then every band is divided by the norm of the entire band.

HSTI.normalize(data_cube, axis = 's')

This function calculates the normalised cube. By setting axis = 's' each spectrum (pixel) will span from 0 to 1 and axis = 'b' normalizes each band individually.

HSTI.remove_stuck_px(data_cube)

This function removes the dead pixels (standard deviation of 0) in the bolometer by replacing them with the average of their non-zero neighbors.

HSTI.remove_outlying_px(data_cube, cut_off)

This function removes outlying pixel measurements of values higher than the cut off value.

HSTI.msc(data_cube, ref_spec = None)

This function applies multiplicative scatter correction to entire datacube. If no reference spectrum is supplied, the mean spectrum of the entire datacube is used instead If a spectrum is supplied, it must have same lenght as axes 2 of cube.

HSTI.normalize_cube(data_cube, axis = 's')

This function normalizes the entire datacube by setting the minimum value to 0 and the maximum to 1.

HSTI.debend(data_cube, central_mirror_sep)

This function takes a single HSTI as input and returns a new spectral bending corrected cube. This does however require a vector containing the mirror separation corresponding to each band in the cube.

HSTI.debend_single_band(masks, img, interp_spec, mirror_sep)

Function used as part of the HSTI.debend() function and is not meant to be run on its own.

HSTI.subtract_band(data_cube, band)

This function subtracts the selected band (given as a integer) from the remaining bands in the datacube, effectively setting that band to 0.

HSTI.targeted_median_filter(array_2D, px_idx, kernel_size)

Median filters single numpy 2D array (array_2D). A kernel size must be supplied along with a numpy boolean array with same dimensions as array_2D. The "True" elements imply where the filtering is performed.

HSTI.array2rgb(data_cube, three_layers)

Select three layers (either as a list or numpy vector) and use these three as the channels of an rgb image. The first layer is the red channel, the second layer the green channel and the third layer is the blue channel.

HSTI.apply_NUC_cube(data_cube, sensor_temp, GSK, NUC_directory = 'default')

This function calculates and applies a NUC to the entire datacube. The NUC is dependent on the sensor temperature and the GSK settings of the camera. The NUC is calculated from camera specific calibration files from the accompanying NUC directory.

HSTI.apply_NUC_image(image, sensor_temp, GSK, NUC_directory = 'default')

This function calculates and applies a NUC to a single image. The NUC is dependent on the sensor temperature and the GSK settings of the camera. The NUC is calculated from camera specific calibration files from the accompanying NUC directory.

HSTI.remove_vignette(data_cube)

This function removes the intensity vignette of an entire data cube. It does not calculate vignette based on the image itself. It assumes that the image is taken using a specific camera and applies a predetermined vignetting correction supplied in the accompanying resources files. This function can both be used on single images as well as entire cubes.

---

Common analysis

HSTI.animation(data_cube, title_list, fps, DPI, fixed_color = True)

This function generates an animation of a datacube. It scrolls through each layer of the cube with a framerate given as (fps) and the title of each frame is supplied in the list (title_list). The (DPI) input is used to control the resolution of the figure, and the color map can be chosen as with any matplotlib figure. If (fixed_color) is set, then the color map spans the range of the entire datacube and does not change for each image. If (fixed_color) is set to False, the color scale spans the entire image of a single frame and is recalculated for each frame in the cube.

HSTI.density_scatter_plot(X, Y, n_bins = 100, colormap = 'plasma', ax = None)

This function takes in two lists (X and Y) and produces a density scatter plot, where the color is controlled by the density of points. The 'resolution' can be controlled by the number of bins given as (n_bins). If no axis is given in (ax), a new figure is created.

HSTI.flatten(data_cube)

This function flattens the datacube into a two-dimensional array.

HSTI.inflate(array_2d, n_rows, n_cols)

This function inflates the datacube into a three-dimensional array. The dimensions of each layer in the datacube must also be provided as n_rows and n_cols.

HSTI.conf95lim(x)

Calculates the upper and lower 95% confidence limits of the mean of input vector, x. The output is a list with the first element representing the lower bound and the second the upper bound.

HSTI.hottelings(X)

Calculates hottelings T^2 statistic for matrix X, where the variables are represented by each column and the samples by the rows https://learnche.org/pid/latent-variable-modelling/principal-component-analysis/hotellings-t2-statistic

HSTI.remove_zeros(array_2D)

This function is used along with resizing images from RGB and thermal camera During the resizing, some gaps might appear, where some pixels are not given a new value and therefore have a value of zero. This function finds these pixels and replaces them with the average of their non-zero neighbors using the avg_neighbors(array_2D, row, col) function.

HSTI.avg_neighbors(array_2D, row, col)

This function calculates the average value of all neighbors of the element given by row and col.

HSTI.fps(points, n_seeds)

Function which distributes n_seeds (a number of points) equally within a lists of points to obtain furthest point sampling.

The function takes in a list of points. Every entry in the list contains both the x and y coordinate of a given point. It returns the coordinates of the selected sample points.

HSTI.voronoi(array_2D, n_seeds)

This function accepts a 2-dimensional array (array_2D) and splits it up into N (n_seeds) subdomains. The partitioning is done based on furthest point sampling.

HSTI.lst_mse(lst1, lst2)

This function returns the mean square error (MSE) between two lists of same length.

HSTI.r_sq(y_fit, y_meas)

This function returns the coefficeint of determination (R²) between fit values in list, y_fit, and measured data in list, y_meas.

HSTI.fpi_sim(mirror_sep, lam, temp)

This function uses the Transfer Matrix Method (TMM) to simulate the transmittance and reflectance of the FPI at given mirror separation (mirror_sep) and wavelength (lam). temp is the substrate temperature since the refractive index of Germanium is temperature dependent. This function DOES NOT take broadening due to substrate bending into account. The functions returns the transmittance, the reflectance, as well as the numeric loss of the simulation.

HSTI.fpi_sim_matrix(mirror_sep, lam, temperature)

The same as HSTI.fpi_sim(), but instead of single values, this function accepts vectors for mirror separation and wavelength. It then returns a 2D array of transmittance values of the FPIwhere each row represents a specific mirror separation while the each column indicates individual wavelengths. This function does take broadaning into account, but only returns the transmittance matrix. Since no loss is included in the model, the reflectance can be found by subtracting the transmittance matrix from a similar-sized matrix of ones.

HSTI.fpi_sim_matrix_angular(mirror_sep, lam, temperature, angle_in_deg)

This function is similar to HSTI.fpi_sim_matrix(), but also takes the angle between the incoming ray and the FPI as an argument.

---

Classes

HSTI.PCA()

The PCA class is used to calculate the principal component analysis of a 2D data structure. An instance of the class contains the following member variables: loadings, scores, singular_values and expl_var_ratio. These are calculated using the calculate_pca(array_2D) method.

.calculate_pca(self,array_2D)

This function calculates the PCA of a two-dimensional input using the covariance matrix of the input array (array_2D). This updates the member variables of the PCA class instance.

.apply_pca(self,matrix)

This function calculates the scores of a 2D array, if the loadings of the given instance of the PCA class is applied to it.

HSTI.least_squares_methods.GLS()

Class used for calculating Generalized Least Squares which can be used to calculate the contributions of target and clutter signals. The class contains the following member variables: spectra, contributions, error, X, X_hat, and lam. spectra is a 2D numpy array, with the first column containing the target spectrum and the second column containing the clutter spectrum. contributions is 2D numpy array, which contain the contribution of the target spectrum in the first column and the contributions of the clutter spectrum in the second column. X and X_hat are the the original 2D array and its estimation respectively. error is the difference between X and X_hat. lam is the regularization parameter used to regularize the eigen values of the covariance matrix of X.

.calculate_contributions(matrix, target, regularization_parameter = 1e-5)

The data matrix is structured to have each sample represented by a row with each column representing a single wavelength. A target spectrum is supplied as a 1D numpy array, and if no clutter spectum is supplied, it is then calculated as the mean of all the data. The data covariance matrix is used to downweigh the most prominent directions within the data such that the target stands out more easily. For this downweighting, the data must be regularized. This is done by adding a constant value (regularization_parameter) to the eigen values of the clutter covariance matrix. This value is determined by the regularization_parameter The contributions contain two vectors. The first describe the contribution of the targeted while the second describe the contribution of the clutter.

HSTI.least_squares_methods.ALS()

This class is used to calculate Alternating Least squares. It contains the following member variables: .spectra, .contributions, .error, .X and .X_hat. These can be calculated using either the calculate_ALS_from_spectrum() or calculate_ALS_from_contribution() methods.The input matrix is given by .X and represents spectral measurements with each individual spectrum given in its rows. Each row represents a single spectrum. .X is factorized into a matrix of contributions and a matrix of spectra as well as the error not captured by the model. A number of constraints can be set. These include nonnegativity for both spectra and contributions and closure can be imposed on the contribution profiles meaning that they will sum to 1. Finally, unimodality can be imposed for either the contribution or the spectral profiles such that each solution/spectrum only has a single maximum. The ALS algorithm requires an initial guess which can be either spectral profiles or contribution profiles. The number of supplied spectra also determine how many components are estimated by the method.

The ALS class contains the following member variables: .spectra, .contributions, .error, .X and .X_hat. Each column of .spectra is the estimated spectral profiles. The ith element in the jth row in .contributions indicate the contribution of the ith pure spectrum (ith column in .spectra) based on the jth measured spectrum (jth row in .X).

.calculate_ALS_from_spectrum(matrix, spectra, closure = False, nonnegativity_c = False, nonnegativity_s = False, unimodality_s = False, Contrast_weight = 0.00, thresh = 1e-8, max_it = 50)

The initial guess determines the number of components, the model will search for in its solution. The spectra must be supplied in columns of input vector .spectra. If .spectra has three columns, three components are calculated. closure can be enforced meaning that the contributions of the spectra must add up to 1. nonnegativity_c and nonnegativity_s indicate whether nonnegativity constraint is enforced in contributions or spectral profiles respectively. unimodality_s enforces a single peak in the estimated pure component spectra. The contrast_weight can be used to increase the contrast in either spectral profiles or in the contribution profiles. If the difference between the original data matrix, .X and the estimate, X_hat, is smaller than thresh, the algorithm is terminated. Otherwise it runs until the maximum number of iterations (max_it) is reached.

.calculate_ALS_from_contribution(matrix, contributions, closure = False, nonnegativity_c = False, nonnegativity_s = False, unimodality_s = False, thresh = 1e-8, max_it = 50)

Same as .calculate_ALS_from_spectrum but instead of supplying spectral profiles, contribution profiles are provided.

HSTI.FPI(mirror_seps = [], lams = []. temperature = None)

Class for performing FPI calculations such as converting from FTIR to FPI spectra or vice versa. All is based on the assumption of loss less mirrors and uses transfer matrix method for calculating light/mirror interactions. Each instance of the class contain the following member variables: .trans_matrix is a 2D transmission matrix showing the transmission of the FPI at given combinations of mirror separation and wavelengths. The .mirror_seps variable contain the mirror separation axis, while the .lams contain the wavelength axis. .temperature must be provided as this affects the refractive index of the germanium layer used in the mirrors. The .camera_response variable can be set to contain the response of the camera. This can then be included in the .trans_matrix if desired.

.set_trans_matrix(trans_matrix)

This member function sets the .trans_matrix member variable to be equal to input.

.set_mirror_seps(mirror_seps)

This member function sets the .mirror_seps member variable to be equal to input.

.wavelengths(lams)

This member function sets the .lams member variable to be equal to input.

.get_trans_matrix()

This function returns the .trans_matrix member variable.

.get_mirror_seps()

This function returns the .mirror_seps member variable.

.get_wavelengths()

This function returns the .lams member variable.

.wl2ms(wavelength_spec)

This function returns a mirror separation dependent spectrum based on a wavelength dependent input. This corresponds to what the camera outputs if presented with that spectrum.

.ms2wl(mirror_sep_spec, correction_factor)

This function returns an estimated wavelength dependent spectrum based on the mirror separation dependent input. This function can be used to estimate what the original spectrum looks like based on a spectrum output by the camera. Since this estimation requires taking the inverse of a matrix, regularization can be necessary. A correction factor must therefore also be supplied.

.apply_camera_response(camera_response)

This function applies the response of the camera to the .trans_matrix variable. The camera response must be supplied either as a string pointing to the location of a pickle file or given as a numpy array with the same dimensions as the .lams member variable.

.remove_camera_response()

This function removes the camera response from the .transfer_matrix member variable.

HSTI.Stack(airgap, lam)

This class is used for calculating the electric field inside a dielectric stack. The original purpose of this class is to calculate the dimensions of a graphene radiation sensor, but it can be changed to any kind of thin film stack. By default this stack consists from the left of a ZnSe layer, Ge, ThF4, Ge, Air and Gold layer. The boundary layers are set to be infinitely thick. The thicknesses of the other layers are set to match that of the thermal mirrors used in the hyperspectral thermal camera. The member variables are .lam0 = 10.5e-6, which is the wavelength used to calculate the thicknesses of the layers in the default thin film stack. .temperature is used to calculate the refractive index of the germanium layer. .gap is the size of the airgap between the final germanium layer and the gold surface. .layers contain the refractive index of each layer in its first column and their thicknesses in the second column. .layer_material is simply a list of strings indicating what material is used in each layer. .transfer_matrix is the transfer matrix of the entire stack calculated using the transfer matrix method. .lam is the wavelength used to calculate the transfer matrix.

.update_indices(lam)

This function updates the refractive index of each layer in the stack based on the supplied wavelength (lam).

.set_airgap(gap)

This function sets the size of the airgap based on the supplied lenght (gap) - specifically of the 5th layer in the cube

.set_wavelength(lam)

This function updates the .lam member variable and subsequently calculates and updates the transfer matrix (.transfer_matrix) of the stack based on the supplied wavelength (lam).

.stack_transfer_matrix(lam)

This function updates .lam and the refractive indices of the first column of .layers. Then the transfer matrix (.trans_matrix) is updated and returned.

.transmittance(lam)

This function returns the transmittance of the entire stack - The E-field on the rightmost side of the stack

.reflectance()

This function returns the reflectance of the entire stack - The E-field on the leftmost side of the stack

.loss()

This function returns the difference between reflectance and transmittance.

.refrac_ThF4(lam)

This function returns the refractive index of ThF4 at the supplied wavelength (lam).

.refrac_Ge(lam, temp)

This function returns the refractive index of Germanium at the supplied wavelength and temperature (lam and temp).

.refrac_Au(lam)

This function returns the refractive index of gold at the supplied wavelength (lam).

.interface(n1, n2)

This function returns 2x2 matrix describing the interface between two layers of refractive index n1 and n2 respectively. n1 is the refractive index of the leftmost layer, while n2 is the refractive index of the rightmost layer.

.ac_phase(n, d, lam)

This function calculates the accumulated phase as the light of wavelength (lam) travels the distance (d) through a medium of refractive index (n).

.E_at_x(lam, x)

This function returns the complex electric field at a given wavelength (lam) and a given coordinate (x) within the stack. The returned E-field is represented as a column vector with the first element representing the E-field of the rightwards moving wave, and the second element representing the E-field of the leftward moving wave. x must lie inside the stack: 0 <= x <= total thickness of stack.

---

Contact

For bug reports or other questions please contact [e-mail] or [e-mail].