HIA

The program Hyperspectral Image Analysis was created to analyse images with up to 3 spatial dimensions and one spectral dimension. It is designed specifically to analyse coherent Raman data but it can be used also for other types of hyperspectral imaging methods, i.e. spontaneous Raman and fluorescence. It was designed and implemented by Francesco Masia and Wolfgang Langbein (2013-2015). Please contact F. Masia (masiaf@cf.ac.uk) for further information and bug reports.

Running the software

The software is written in MATLAB, compiled for Windows 7 64bit, and requires the MATLAB run-time compiler version 9.1 (2016b) to be installed. The java class DropTargetList.class, the parallel computing profile HIA.settings and the library minilzo.dll must be in the same folder as the executable. From version 2.5, to open the Andor data *.sif, it requires the installation of the Andor Matlab SIF i/o toolbox (AndorMatlabSIFIOSetup-2.102.30000.0). The installation folder must be the Runtime compiler folder, for example YourDrive\Program Files\MATLAB\MATLAB Runtime\v91.
To run the software, launch the batch file run_HIA.bat. This will open a cmd window which can show any MATLAB error occurred during the analysis. Please use this feature to identify and report bugs.

HIA 2.5, 21/02/2017, Francesco Masia


Enter number of CPUs

Enter number of CPUs and parallel computing idle time

HIA uses parallel computing in some parts. The user has the option to select the number of threads to be used by HIA during parallel computing and the idle time (in minutes) before the parallel computing will turn off. Calling a function using parallel computing will re-start it. The default value is the total number of CPU cores present and 180 minutes for the idle time.

Enter number of CPUs

The image below indicates the functions using parallel computing.

HIA 2.0, 10/02/2015, Francesco Masia


Main panel

HIA Main panel

This is the main panel of the HIA program. It gives access to all the functionality. A display of program messages shows information from the program about the calculation stage.

HIA main panel

The individual controls are:

MultiCARS Main
  • Data # : number of data sets available;
  • number: selected data set numbered from 1 to #. It defines the active dataset where the operations are performed.
  • arrows : The arrows allow to navigate within different available datasets.
  • open: Load a data set. Opens a file selection box. The accepted file formats are tiff (8 and 16 bit), Bioformats tiffs, Multicars (ASCII and srg), ASCII, Andor sif, HIA projects (*.hia). The loaded data sets gets the number # , and # is increased by one. See Load a TIF image for more information.
  • save: Save the selected data set. Opens a file selection box. The possible file formats are tiff (8 and 16 bit), and HIA (*.hia). See save a HIA project for more information.
  • delete: Delete the selected data set. The number # is reduced by one and all data sets with numbers above are reduced in number by one.
  • view : View the selected data set. Opens a Display panel.

The menu bar choices are

  • File
    • Load
      • TIF Image
      • BioFormate image
      • ASCII and Andor sif
      • HIA project: Loads an HIA project.
      • MultiCARS (Ascii)
      • MultiCARS (srg)
Alternatively the user can drag and drop suitable files on the main panel to open them
    • Save (active when images are available)
      • TIF image
      • BioFormate image
      • HIA project
  • Image (active when images are available)
    • Type
    • Scaling
    • Subtract offset
    • Sensitivity correction
    • Limit spectral range
    • Multiple spectral selection
    • Spectral interpolation
    • Exclude Regions
    • Join Images
    • Wavelength to Wavenumber
    • Cosmic Ray Removal
    • Normalize image
    • Image binning
    • Crop image
    • Register image along energy axis
    • Register image along z axis
    • Image Calculator
    • Save region spectrum
    • Store offset spectrum
    • Fluorescence corrected image
    • OffsetMask
  • Analysis (active when images are available)
    • Sparse sampling
    • Uneven background correction
    • CARS ratio calculation
    • SVD denoising
    • Field retrieval (PCKK)
    • FSC3
    • Project spectra
    • Project spectra with Concentration constrain
    • Background subtraction (robust fit)
HIA 2.5, 09/02/2017, Francesco Masia


Data display panel

The active hyperspectral image can be displayed by clicking the view icon from the Main panel. Multiple such panels can be opened at the same time. For complex susceptibility data (i.e. after Field retrieval) the imaginary part is displayed.

The follwing data display panel is opened:

Data display panel
  • The slides at the bottom of the panel allow to navigate in the spectral dimension (left) and z dimension (right, only for 4D data). The indicators show the value of the spectral frame and z position respectively, and can also be used to enter the desired value.
  • zoom: Zoom in and out the data at the actual spectral and z position. Click on the icon and then on the image. A marquee area can be defined.
  • cursor: Selects a pixel at the actual spectral and z position. Click on the icon and then on the image. Multiple pixels can be selected. To add cursors use the right click and select 'Create new Datatip'. Alternatively Alt+left click.
  • Plot spectrum: Plots the spectra at the selected cursors. Opens a Plot spectrum panel. If more than one images is displayed, spectra from different images can be plotted in the same panel by selecting pixels at the desired panels.
  • Display real part: Only for complex susceptibility data (i.e. after Field retrieval). Opens a new data display for the real part. For 4D data the slide and control of the z dimension is locked to the selection of the imaginary part display panel.
  • Grayscale: Shows the linear grayscale with start and end values used to display the pixel values in the image.
HIA 2.0, 12/02/2015, Francesco Masia


Plot spectrum panel

The panel shows the spectra at the active cursors position of a Data display. For complex susceptibility data (created after Field retrieval) the imaginary part is displayed.

Plot spectrum panel
  • Save spectra: Saves spectra as an ascii file. Spectra from the same image will be saved in the same file.
  • Update: update the display.
  • Label: Shows the data # and positions of the displayed spectra.
HIA 2.0, 12/02/2015, Francesco Masia


Load TIF image

Loads a 8 or 16 bit stack/hyperstack TIF image. Hyperstack TIFs allow to store more than two-dimensional data as channels, slices and frames. The software assumes an hierarchy of the dimensions having the detector channels on a higher hierarchy compared to the z and spectral dimensions. If the image is a 5D data, the software assumes that the channels refer to different detectors, the slices to the z dimension and the frame to the spectral dimension. The user must indicate the channel to import. If the image is a 4D data the software interprets the higher hierarchy stack dimension as the z dimension and the other one as the spectral dimension. If the image is a 3D data the third dimension is interpreted as the spectral dimension. If the TIF was saved by HIA, the information about pixel, intensity and spectral scales will be recovered. By the default the data are assumed to be CARS data.

HIA 2.0, 10/02/2015, Francesco Masia


Load BioFormat image

Loads BioFormat stack/hyperstack TIF image. The software assumes the hierarchy XYZCT, with T being the spectral dimension. If a different hierarchy is used the function returns an error. If there is more than one series, the user will be prompted to select the series with the data. If there is more than one channel, the user will be prompted to select the channel with the data. By the default the data are assumed to be CARS data.

HIA 2.3, 21/03/2016, Francesco Masia


Load Ascii

Load ASCII/Andor data

Loads an ascii file. The data should be ordered as a column. If more than one column is present, the user must indicate the column to import.

Opens a panel to insert the parameters to load the ASCII data.

  • Load parameters: Loads the parameters for the panel from a .par file.
  • Save parameters: Saves the parameters for the panel to a .par file.
  • Reset parameters: Reset the parameters to the default values.
  • Data ordering: Indicate which is the hierarchy of the scanned dimensions. There are two options.
    • Spectrum first
    • Position first
  • Spectral # points: Insert the number of spectral points
  • Calibration: Insert the spectral calibration where pixel is the pixel number starting at 1.
  • Spectral unit: Insert the unit of the spectral dimension
  • Sensitivity correction: tick the box if data needs to be corrected by the system spectral sensitivity.
  • Spatial axis i:
    • #points: Insert number of points of the spatial axis i
    • Pixel size: Insert corresponding pixel size
  • Repetitions: Insert number of repetitions of the same rasterscan. The function assumes that the rasterscan repeats the acquisitions before changing the position of the second spatial axis (The rasterscan can be set to have zero step size in the x direction and the number of points along the x direction is the number of repetitions. The different repetitions are imported as different data.
    • For sif data acquired with the Andor CCD in linescan mode use the number of rows as first dimension and the number of frames as second dimension (or a combination of second and third dimension depending on the rasterscan settings). The number of points along the x direction in the rasterscan is the number of repetitions.
  • Spatial unit: insert spatial unit
  • Camera sensitivity: Insert camera sensitivity (in electrons/count).
  • Camera Read Noise: Insert camera read noise (in electrons).
  • Cosmic Ray Removal: Tick the box to start an automatic cosmic ray removal within the data import. (active only for camera read-noise>0). See Cosmic Ray Removal for more details
  • Laser wavelength: Insert laser wavelength in the data spectral unit.
  • Conversion wl → wvn: Tick the box to include a conversion of the spectral axis from wavelength to wavenumber. Active if the laser wavelength is greater than 0 and the Spectral unit is um or nm.
  • Offset: Indicate the type of offset to subtract. See Subtract Offset for more details. Four options are possible
    • Constant
    • Range
    • Image
    • Average Image
  • Offset value: Indicate the value of offset (in counts) to subtract (Active for Constant offset type).

The data import function workflow is shown below

Load ascii Workflow

By default the data are assumed to be Raman data.

HIA 2.5, 09/02/2017, Francesco Masia


Load MultiCARS file

Load MultiCARS Ascii file

Loads ASCII-exported regularized MultiCARS data. By the default the data are assumed to be CARS data.

In this version only linear and noise compression scalings of values are supported.

Parallel computing is used if the values are stored using the noise compression scaling.

parallelcomp

HIA 2.0, 11/02/2015, Francesco Masia


Load MultiCARS file

Load MultiCARS srg file

Loads binary regularized MultiCARS data (extension srg). By the default the data are assumed to be CARS data.

In this version only linear and noise compression scalings of values are supported.

Parallel computing is used if the values are stored using the noise compression scaling.

parallelcomp

HIA 2.4, 05/12/2016, Francesco Masia


Save TIF image

Saves a 8 or 16 bit stack/hyperstack image depending on the dimensions in the image. The spectral dimension is saved as frames, and the z dimension as slices. Some information about pixel, intensity and spectrum scales is stored in the image description (ctrl+i in imageJ/FiJi).

HIA 2.0, 10/02/2015, Francesco Masia


Load BioFormat image

Save BioFormat image

Save the data as 32bit float BioFormat stack/hyperstack TIF image. The software will usw the hierarchy XYZCT, with T being the spectral dimension.

HIA 2.4, 06/12/2016, Francesco Masia


Save HIA project

Saves an HIA project. The compression of the data is based on a SVD fitlering The software will perform an SVD on the selected data (similar to the SVD denoising function) and the user will select the number of components to be retained for the saving. The output file is a matlab structure composed by the reduced SVD spectral basis U and the matrix Σ V' saved as integer (8bit or 16bit depending on the user selection). The information about the image properties will also be saved.

HIA 2.0, 10/02/2015, Francesco Masia


Image Type

Selects the type of data of the current image. Options are CARS, SRS or Raman (use the latter for fluorescence data). The choice will determine which functions will be available and what default options will be displayed. For SRS and Raman data, the functions "CARS Ratio calculation" and "Field retrieval (PCKK)" will not be available

HIA 2.0, 10/02/2015, Francesco Masia


Scaling

The scaling of the position and spectrum dimensions and pixel values for the active image can be modified with this function. There are three submenus:

  • Pixel: Insert the pixel dimensions and unit
Position scaling
  • Spectrum: Insert the spectrum step size, origin and unit
Spectrum scaling
  • Value: Insert the maximum and minimum of the pixel values
Value scaling
HIA 2.0, 10/02/2015, Francesco Masia


Subtract offset

This function allows to subtract an offset from the selected data. There are four options:

  • Constant: Subtracts a constant offset. The user is asked to give the value of the offset.
  • Range: It calculates an offset image using the mean value in a defined spectral range. The user is asked to indicate the spectral range using a slider (see Limit spectral range for more details)
  • Image: The user is asked to load an offset image. The image must have the same spatial dimensions and no spectral dependence. If the image has different number of spatial pixels, the user has the option to subtract an sptially average offset, otherwise no action is taken.
  • Average image: The user is asked to load an image. The offset will be calculated as the mean value of the image.
  • Spectrum: The user is asked to load a file with a spectrum to be removed as an offset to all the spatial voxels of the image. The loaded spectrum must have a spectral range larger than the data. If the spectral axis is different from the data a linear interpolation of the spectrum is performed
  • Stored spectrum: Remove the spectrum stored in the memory to all the spatial voxels of the image.
The user can decide to overwrite the selected dataset or to create a new one with the modified data.
HIA 2.6, 19/10/2018, Francesco Masia


Sensitivity correction

This function allows to correct the hyperspectral data by the spectral response of the system. The user is asked to load an ascii file. The file must have the spectral positions in the first column and the corresponding value of the sensitivity in the second. The spectral range of the sensitivity must include the one of the data to be corrected. The program will interpolate the sensitivity to match the spectral positions of the data.
The user can decide to overwrite the selected dataset or to create a new one with the modified data.

HIA 2.0, 12/02/2015, Francesco Masia


Limit spectral range

This function allows to limit the spectral range of the selected data.

It opens a slider panel (shown below).

Limit spectral range

The user can type the new spectral range in the controls or use the slider to define the new range. The data outside the range will be deleted. The user can decide to overwrite the selected dataset or to create a new one with the modified data.

HIA 2.0, 12/02/2015, Francesco Masia


Limit spectral range

Multiple spectral selection

This function allows to select multiple spectral ranges of the selected data.

It opens a slider panel (shown below).

Multiple spectral ranges

  • # Sections: defines the number of sections to consider for the selections

The user can type the ranges in the controls or use the sliders to define the new ranges. The data outside the ranges (regions in gray) will be deleted. The user can decide to overwrite the selected dataset or to create a new one with the modified data. This function is active only for SRS and Raman data type and for CARS data after field retrieval.

HIA 2.3, 21/03/2016, Francesco Masia


Spectrum interpolation

This function allows to linearly interpolate the selected data into new spectral positions.

A panel allows to enter the new spectral origin and end and the step size.

Spectral interpolation

The user can decide to overwrite the selected dataset or to create a new one with the modified data.

Parallel computing

HIA 2.0, 13/02/2015, Francesco Masia


Exclude regions

This functions allow to exclude spatial regions from the analysis. This is used in the singular value decomposition for denoising and in FSC3. It can also be used to reconstruct data values which deviate from the typical spectra observed in the sample using the SVD based masking/reconstruction function.

  • Type:. The selection can be moved after drawing. To activate the selection double-click in it.
    • Mask: exclude the regions inlcuded in the selection.
    • Crop: exclude the regions outside the selection.
    • Unmask: remove exclusion the regions inlcuded in the selection.
    • Invert: exchange between excluded and non-excluded regions.
    • Reset: remove all exclusions of the image
  • Selection
    • Freehand: Select points using a freehand tool. The selection can be moved after drawing. To activate the selection double-click in it.
    • Rectangular: Select points using a rectangular tool. The selection can be moved after drawing. To activate the selection double-click in it.
    • Polygonal: Select points using a polygonal tool. The selection can be modified and moved after drawing. To activate the selection double-click in it.
    • Oval: Select points using an oval tool. The selection can be moved after drawing. To activate the selection double-click in it.
    • Regular grid: Select points using a fraction of pixels in a regular grid. The fraction of pixels to be selected is defined in the following panel
Exclude regions - equidistant

    • Real part range (active only for complex susceptibility data): Select pixels which have the value of the real part within the specified range.
    • Load mask: load a tif image which is used as mask.The image will be converted to logical.
    • Save mask: save the actual mask as a 8bit tiff.
    • SVD based masking: uses singular value decomposition to mask pixels with the option to replace the value of the masked pixel with the SVD denoised value.
  • Energy: Defines if the selected area applies to the actual spectral position or a spectral range. Select the right option before selecting the region.
  • z (only for 4D data): Defines if the selected area applies to the actual z position or to a range along z direction. Select the right option before selecting the region.
  • Ok: apply the current masking to the data.

The excluded regions will appear in red in the image. Masking of the data should be applied to the full spectral range. If only a partial range of spectral points has been masked, the user has he option to either apply the mask to the full range or to recalculate the pixel values at those spectral position using the sparse sampling algorithm.

If any point of the data has been modified using either the SVD based masking /reconstruction or the sparse sampling algorithm, the user can decide to overwrite the selected dataset or to create a new one with the modified data.

HIA 2.3, 23/03/2015, Francesco Masia


SVD masking/reconstruction

This function allows to mask (and reconstruct) datapoints using singular value decomposition.

There are two options:

Masking

SVD masking parameters

The user selectes the number of rms deviations (ξ) which will be used as a threshold for the masking.

SVD filtering is performed on data D0=F(d), where d is the active data and the function F is defined by the parameters input in the SVD panel, i.e. F(x)=sign(x)α(xβ+γ)δ, and the SVD filtered data D* is calculated from D0. The data points where abs(D0-D*)>ξσ are masked, where σ=||D0-D*||/(#D0)0.5 and ||.|| indicates the Frobenius norm of the reshaped (spatial x spectral) matrix and #D0 is the number of elements in D0. Points which were excluded outside of the function are not included in the calculation.

Reconstruction

This is an iterative function which first determines the outliers by using the approach described for the masking, and then substitutes them with the value of D*. The iteration stops when the number of points to be reconstructed chenges less than a defined relative threshold in one step.

SVD reconstruction paramters

The user defines the number of rms deviations (ξ) which will be used as a threshold for the masking, and a threshold f of the minimum relative change of number of masked points per iteration.

The workflow of the iterative algorithm is shown below

SVD reconstruction workflow

j and k indicate the indices of the spatial and spectral dimensions respectively.

HIA 2.0, 26/02/2015, Francesco Masia


Join images

This function allows to join multiple images along the spectral dimension.

The user must select the data to join. The data to join must have the same number of points in each spatial dimension.

Join images

The user must then indicate the new spectral dimension range (similar to the Spectrum interpolation). The software will limit the new range to the combined range of the selected images.

If more than one data have values at certain spectral position (see circled dot in Fig1) the software will shift spectral positions by 1e-4 of the step size (see Fig2, shift is exagerated) before joining them by linear interpolation (see Fig3).

Fig1
Fig1

Fig2
Fig2

Fig3
Fig3

A new dataset will be generated.

Parallel computing
HIA 2.0, 13/02/2015, Francesco Masia


Wavelength to wavenumber

This function allows to rescale data measured in wavelength into a linear scale in wavenumber relative to a laser wavelength. This is suited to calculate the Raman shift for Raman spectra measured with a spectrometer on a wavelengths scale. Note that this is only a axis rescaling, not a corresponding change of the values according to the change of axis differentials. The latter should be done by a sensitivity correction seperately determined if needed. The function is active if the unit of the spectrum scale is either 'um' or 'nm'.

First the user has to insert the laser wavelength in the same unit as the data (this panel does not pop-up if the function is called within the importing of ASCII files)

Laser wavelength panel

then indicate the new spectral axis settings for the interpolation.

Wavenumber interpolation

The function will limit the range to the data range, which is given as default values.
The user can decide to overwrite the selected dataset or to create a new one with the modified data.

Parallel computing

HIA 2.0, 17/02/2015, Francesco Masia


Cosmic ray removal

This function allows to remove cosmic rays from data acquired with cameras. Data values have to be in units of electrons as a poisson noise is assumed.

Cosmic ray removal panel
  • Allowed pixel to pixel dynamics: defines the maximum variation (Δ) between adjacent pixels not considered a cosmic ray
  • Check along: defines the dimension to check for the pixel to pixle varation. Options are:
    • Position: compare the actual pixel with the 8 (26) adjacent pixels for a xy (xyz) image.
    • Spectrum: compare the actual pixel with the 2 adjacent pixel in the spectral direction.
  • Replace with: defines how to replace the pixel value in case a cosmic ray has been detected. Options are
    • Minimum adjacent: minimum pixel value among the adjacent pixels.
    • Average adjacent: average pixel value among the adjacent pixels.
    • NaN: not a number. The data at the position will be set to unknown.
  • Camera Read Noise: Insert camera read noise (in electrons) (σRN).

A pixel value (d) will be interpreted as cosmic ray if d-(d+(σRN2))1/2>Δ[d+(d+(σRN2))1/2] for is any adjacent pixel value.

An example of the workings of the function is shown below. The original data (left) contains of a cosmic ray pixel, which is corrected in the image on the right after cosmic ray removal. The spectra of the untreated (Data1) and modified data (Data2) are shown. also. The pixel values not recognised as cosmic ray are not modified.

Cosmic ray removal example- before Cosmic ray removal example - after

Cosmic ray removal example - spectra
The user can decide to overwrite the selected dataset or to create a new one with the modified data.
HIA 2.0, 18/02/2015, Francesco Masia


Normalize image

This function allows to normalize the image using the image value of one of the spectral positions, or the real part value if the data are complex susceptibility.

The normalize image has unity mean. If there are excluded regions in the image, they will be excluded in the calculation of the mean value. There is the option to mask regions corresponding to a normalization factor less than a threshold.

Normalize image
  • Frame: defines the spectral position to consider for the normalization. For complex susceptibility data the real part can be also used
  • Mask region if normalization factor less than: defines the threshold value which is used to mask regions where the normalization is smaller.
  • Ok: accept the selection and execute normalization.
The user can decide to overwrite the selected dataset or to create a new one with the modified data.

HIA 2.0, 03/03/2015, Francesco Masia


Image binning

This function allows to bin the active image along the spatial dimensions

A panel allows to enter the binning size and offset.

Image binning

  • Offset: indicates the number of pixels to skip alog a certain direction
  • Bin size: indicates the size of the binning along a certain direction

The user can decide to overwrite the selected dataset or to create a new one with the modified data.

HIA 2.3, 23/03/2016, Francesco Masia


Crop Image

This function allows to crop the image to a smaller region.

The new data region can be selected using a rectangular tool. The selection can be moved after drawing. To activate the selection double-click in it.
The user can decide to overwrite the selected dataset or to create a new one with the modified data.

HIA 2.4, 06/12/2016, Francesco Masia


Registration along the energy direction

This function allows to register the image along the spectral direction to correct for motion during the acquisition.
There are two options

  • Calculate track
    • Use this option to calculate the registration track. For 4D data, individual registration for each z slide are calculated. The regions without data are filled with NaNs. The user can decide to overwrite the selected dataset or to create a new one with the modified data.The suer can decide to save the registration track as a .reg file (It is a Matlab 'mat')
  • Load track
    • Use this option to apply a pre-calculated track to the data. The registration track needs to be calculated from a dataset with identical dimension, otherwise the function will return an error. The regions without data are filled with NaNs. The user can decide to overwrite the selected dataset or to create a new one with the modified data.

HIA 2.4, 06/12/2016, Francesco Masia


Registration along the z direction

This function allows to register the image along the z direction to correct for rigid motion during the acquisition.
First the user must indicate the spectral range to use to calculate the registration track.

Individual tracks are calculate for each spectral point.The average motion is calculated and applied to the data. The regions without data are filled with NaNs. The user can decide to overwrite the selected dataset or to create a new one with the modified data.

HIA 2.4, 06/12/2016, Francesco Masia


Image calculator

This function allows to do some basics calculation between two images.

  • Image 1: Selects the first image to consider in the calculation between the available datasets.
  • Image 2: Selects the second image to consider in the calculation between the available datasets.
  • Image #: Defines the operations applied to image 1 and 2 used to create the new image
  • Continue: Performs the calculation and create the new dataset

HIA 2.4, 06/12/2016, Francesco Masia


Save region spectrum

This function allows to save a spatially averaged spectrum.

It uses the same selection tool as the Exclude regions function

HIA 2.6, 19/10/2018, Francesco Masia


Store offset spectrum

This function allows to store a spatially averaged spectrum to be used as a offset.

It uses the same selection tool as the Exclude regions function.

The spectrum can be subtracted from each pixel using the "Stored spectrum" option of the Subtract offset function.

HIA 2.6, 19/10/2018, Francesco Masia


Fluorescence corrected image

This function allows to correct CARS intensity data with a strong spectrally dependent fluorescence background (the spectral dependence can usually generate from photobleaching of the fluorescence during the IFD scan).
The following data are needed:

  1. Hyperspectral CARS intensity collected in forward direction (fC)
  2. Hyperspectral CARS intensity collected in epi direction (eC)
  3. CARS intensity collected in forward direction with no pulse overlap (fCo)
  4. CARS intensity collected in epidirection with no pulse overlap (eCo)
  5. Hyperspectral TPF intensity collected in forward direction (fT)
  6. Hyperspectral TPF intensity collected in epi direction (eT)
  7. TPF intensity collected in forward direction with no pulse overlap (fTo)
  8. TPF intensity collected in epidirection with no pulse overlap (eTo)

The user must indicate the index number of the following data

The fluorescence corrected data will be calculated as fC-eC*S*CA, where S is calculated as the ratio of the xy spatially integrated fCo and eCo over a certain ROI, and C is the square root of the ratio of (ifT2+R2*C02) and (ieT2+R2), where ifT and ieT are the spectrally integrated fT and eT, respectively, and C0 is the ratio of the xy spatially integrated fTo and eTo over a the same ROI used to calculate S.
The function calculates the spectrally integrated fT and eT, i.e. ifT and ieT. The user will select a spatial region were to calculate the noise using the same routine as the Exclude region function. The function will ask to provide the values of the regularization factor and the amplification factor.

  • Regularization factor: defines value of R (if left blank, the noise value is used)
  • Amplification factor: defines value of A
  • Noise: returns the largest value of the std of the ifT and ieT
  • Maximum signal: returns the maximum value of ieT.
The user is then prompted to define the ROI for the calculation of S and C0 using the same function as in Exclude region.
HIA 2.6, 19/10/2018, Francesco Masia


Offset and mask

Offset & mask

This function allows to define two threshold to automatically apply an offset to the data and to mask pixels.

For each pixel, the minimum value along the spectral dimension is considered. If this value is smaller than the Offset threshold the spectrum is shifted by the minimum value. If this value is smaller than the Mask threshold the pixel is exlcuded from the relevant functions.

  • Offset threshold: defines the threshold to apply the offset
  • Mask threshold: defines the threshold for the masking
HIA 2.6, 19/10/2018, Francesco Masia


Sparse sampling

This function allows to reduce the number of measured frequencies in the hyperspectral image, and consequently to increase the hyperspectral frame rate, without loosing significant spectral information. First one has to determine the important spectral components of the sample under investigation from a hyperspectral image acquired over a large number of spectral points sampled at the Nyquist limit over a small number of spatial points using singular value decomposition (SVD). For best performance the SVD should be performed on the signal which has whitened noise. The obtained SVD spectra will be used as spectral basis to reconstruct the full spectral information. Subsequently, a limited set of spectral positions (equal to the number of important spectral components) is determined to yield minimum deviation over the reconstructed spectra by using a random walk minimization approach. A hyperspectral image is then acquired over a large number of spatial points only at the few determined spectral positions. Finally a reconstruction algorithm is applied to generate at each point of the image a spectrum with the full spectral resolution.

Determination of optimal spectral positions

This part of the function determines the optimal positions for the subsequent hyperspectral image acqusition. The active data should have the full spectral information. The determination of the spectral positions is based first on the determination the important spectral components by SVD. In this case the Maximum index method is used. The maximum index is determined initially (and every time the fit range is changed) with the Threshold method using factor sqrt(2).

Determination of spectral positions

After selecting the number of components, the functions determine the optimal set of spectral positions. There are two options:

Minimization of the reconstruction error

To minimize noise in the reconstructed data, the user has the option to reduce the number of components used by applying a threshold or to weight the contribution of the different components.

  • Threshold: The function calculates the norm of the column vectors of the U-1 matrix, and the user can input a threshold (in units of norm of the first column vector of U-1) which limits the components to consider. The default threshold factor is sqrt(2). In the plot below the components in red are included in the calculation, the blue are excluded and the dashed line indicates the threshold. The user also selects the number of steps in the random walk used to optimize the spectral positions.

  • Weighting: The function calculates the norm of the column vectors of the U-1 matrix, and the weighting for the component i is defined as wi=1/(1+(|U-1i|/(α*|U-11|))β). In the plot below the symbols in blue are the norm of the column vectors of the U-1 matrix (left axis), the red symbols are the corresponding weight (right axis).

    • Weight factor: Determines value of alpha
    • Weight power: Determines value of beta

The random walk start with a set of equidistant spectral points which are moved randomly (the maximum step is adjusted to avoid crossing). During random walk, the minimum relative error reached versus optimization step is displayed.

Random walk

The set of optimal spectral positions is shown once the random walk is finished, together with the error in the reconstruction using the equidistant set and the optimized set.

Results
  • MultiCARS presets: save the list of positions as presets for the MultiCARS software.

Minimization of the norm of U-1

This option optimizes the wavenumbers to minimize the norm of U-1. The user defines the number of steps in the random walk used to optimize the spectral positions. The random walk start with a set of equidistant spectral points which are moved randomly (the maximum step is adjusted to avoid crossing). During random walk, the reached minimum relative change of the norm of U-1 is displayed versus the optimization step. Similarly to the case of the minimization of the reconstructed data error, the set of optimal spectral positions is shown once the random walk is finished and the wavenumber set can be exported as preset for the MultiCARS software.

Reconstruct data

This option allows to reconstruct the spectral information from a hyperspectral image acquired over a large number of spatial points and few number of spectrally points, i.e. over the set of spectral positions determined above. The image with the original spectral information needs to be loaded. The active image is the image to be reconstructed.

Insert the spectral positions

First the user has to insert the spectral positions of the data to reconstruct. By default, the values stored in the image are given.

  • Import from MultiCARS ppr: import the spectral positions from a preset file of the MultiCARS software.
  • Ok: record the positions and continue

Then the data with the full spectral information has to be selected.

Select data with full spectral information

The function then calculates the spectral basis to project the data using SVD on the data with the full spectral information. See the SVD description page for more details.

Similarly to the spectral position determination function, the reconstruction algorithm uses either the SVD components whose norm of the U-1 vectors is below a certain threshold or all the SVD components but with a weight reduced according to the norm of the U-1 vectors (see above).

Finally the function calculates the reconstructed data with the full spectral information.

Test

This option allows to test the sparse sampling algorithm on any dataset

  • Select the spectral positions: Either select one of the data to run the optimal spectral position determination algorithm or the 'Manual' option. The latter allows to choose either equidistant spectral positions or select them from the actual spectral axis. The positions will be used in the reconstruction part of the algorithm.

    Spectral positions

  • Select the data with the full spectral information: for testing the same dataset can be used to get the full spectral information and to be reconstructed.

    Full spectral info

  • Reduce number of spatial points: To simulate the performance of the algorithm for real application, the number of spatial points of the data consider for the full spectral resolution can be reduced.

    Reduce spatial points

  • If decided to reduce the spatial points, insert the percentage and offset of the spatial points to consider

    Percentage

The following steps are the same as for the reconstruction function

  • Spectral basis calculation: using SVD
  • Determination of the SVD components to consider: using threshold or weighting on U-1 norm
  • Data reconstruction
HIA 2.1, 10/06/2015, Francesco Masia


Uneven background correction

This function allows to remove an uneven background from data. There are the following options:

  • Background region determination: This option is based on a manual or automatic determination of regions in the image which are considered as background. After the background regions have been identified, a fit of the background intensity/amplitude is performed at every spectral position. There is a choice between two fitting functions: polynomial or exp(polynomial). In the first case the background fit will be subtracted from the data, in the second case the data will be divided by the background fit. In both cases the zero-th order term of the polynomial can be set to zero before correcting the data, allowing to correct only the inhomgeneities.
    • Manual
    • Automatic (Iterative) - (Largest error)
  • Background function: This option uses given parameters to calculate the background, which is described either as a polynomial function or as an exp(polynomial). As in the background region determination, in the first case the background fit will be subtracted from the data, in the second case the ratio between the data and the fit is considered.
    • Insert background parameters
    • Load background parameters
    • Fit parameters (Smallest error)
The user can decide to overwrite the selected dataset or to create a new one with the modified data.

HIA 2.0, 03/03/2015, Francesco Masia


Manual background region determination

The determination of spatial regions corresponding to background is done manually.

Manual background region determination
  • Function: Defines the function to use for the fit: either polynomial or exp(polynomial).
  • Polynomial order: Defines the order of the polynomial.
  • Exclude zero-th order term from fit for correction: defines if the zero-th order of the polynomial is excluded from the correction of the background.
  • Exclude zero-th order term from fit for correction (real part): (Visible only for complex susceptibility data) defines if the zero-th order of the polynomial is excluded from the correction of the background for the real part.
  • Spectrally independent correction: defines if the background correction is spectrally dependent or not.
  • Save fit: saves a tiff image with the background obtained from the fit.
  • Save fit parameters: save an ascii file with the spectrally dependent fit parameters.
  • Ok: opens the region selection panel

The region selection panel is similar to the Exclude region function. The area in the blue defines what is considered as background. The areas in red corresponds to regions previously excluded (if any) and are not considered as background even if selected. The background region determination is spectrally independent.

After the region of the background is determined, the multidimensional fit is performed. The function uses the voxel size defined in the Scale panel and uses the origin in the center of the image. If the spectrally independent correction option is active, a weighted average of the parameters along the spectral dimension is performed. The weight is defined by the inverse of the error of the fit at the corresponding spectral position. The data are finally corrected by subtracting or dividing by the fit depending if the function chosen was the polynomial or the exp(polynomial). If the "Exclude zero-th order term from fit" option is active the zero-th order term of the polynomial is not included in the calculation of the background before correction. This reflects to a background equals to 0 (1) at the center of the image for the polynomial (exp(polynomial)) function, respectively.

If the "Save fit" and/or the "Save fit parameters" option are active, the user should give the name of the file where to store the background/parameters.

This function is the preferred option if there is a clear spatial separation between background and foreground.

HIA 2.0, 03/03/2015, Francesco Masia


Automatic background region determination

Automatic background region determination (Iterative)

  • Function: Defines the function to use for the fit: either polynomial or exp(polynomial).
  • Initial polynomial order: Defines the initial order of the polynomial of the iterative algorithm.
  • Final polynomial order: Defines the final order of the polynomial of the iterative algorithm.
  • Data pre-masking: defines if a pre-masking routine is used.
  • Maximum number of iterations: defines the maximum number of iterations allowed.
  • Threshold factor (σ): defines the factor used to discriminate background and foreground.
  • Exclude zero-th order term from fit for correction: defines if the zero-th order of the polynomial is excluded from the correction of the background.
  • Exclude zero-th order term from fit for correction (real part): (Visible only for complex susceptibility data) defines if the zero-th order of the polynomial is excluded from the correction of the background for the real part.
  • Spectrally independent correction: defines if the background correction is spectrally dependent or not.
  • Save mask: saves a tiff image with the mask used to define the background.
  • Save fit: saves a tiff image with the background obtained from the fit.
  • Save fit parameters: save an ascii file with the spectrally dependent fit parameters.
  • Ok: starts the function

Iterative method

The determination of spatial regions corresponding to background is done automatically, with an iterative algorithm.
The aim of the function is to automatically mask regions of the image which are interpreted as foreground, fit the unmasked regions and correct the data by the fit. The function uses the voxel size defined in the Scale panel and assume the origin in the center of the image. This function allows to include more than one dataset in the calculation. The datasets must have been acquired in the same conditions (same spatial positions, same spectral positions). The datasets will be displaced along the x direction by an amount equal to a 10th of the step size divided by the number of images considered. Assuming a slow modulation in the background, this would not affect the fit.

The next part of the function is done at each spectral position sequentially.

Data pre-masking

If the data pre-masking option is active the function will create an histogram of the data value at the specific spectral position, and it will find the mode m. If more than one mode is found the number of bins in the histogram is halved until a single mode is found. The deviation of the data D from the mode is calculated as Δ=D-m. Data where the modulus of the deviation |Δ| exceeds σpm*rms(Δ) are masked and not considered in the next step. rms() denotes the root-mean-square deviation. The number of histogram bins and the threshold σpm are entered by the user (see panel below).

Pre-masking panel

Iterative masking

The non-masked data are fitted using the defined function with the initial polynomial order, resulting in the fit F. Data for which the modulus of the deviation |D-F| exceeds rms(D-F) by a factor larger than the threshold parameter σ are masked for the next iteration. The iteration stops when the mask is unchanged in an iteration or when the maximum number of iterations is reached. If the final polynomial order is not reached the order of the polynomial is increased by one and the iteration is continued.
The resulting mask separates the foreground (masked) from background. If "spectrally independent correction" is selected, a weighted average of the fit parameters along the spectral dimension is performed. The weight is defined by the inverse of the rms error of the fit at the corresponding spectral position. The data are corrected by subtracting the fit or dividing by the fit depending if the function chosen was the polynomial or the exp(polynomial). If "Exclude zero-th order term from fit" is selected the zero-th order term of the polynomial is not included in the calculation of the background before correction. This results in a background equal to 0 (1) at the center of the image for the polynomial (exp(polynomial)) function, respectively.

If "Save mask", "Save fit" and/or the "Save fit parameters" is selected, the user is asked for the name(s) of the file(s) to store the corresponding parameters.

This method is the preferred option if there is not a clear spatial separation between background and foreground.

Largest Error method

The determination of spatial regions corresponding to background is done automatically, with an iterative algorithm, but only at the spectral positions corresponding to the largest error in the fit. First a fit of the data will be performed using the selected function and the maximum order of the polynomial. Only the spectral position which returns the largest error will be used to determine the background position, as in the iterative method. The foreground positions are excluded from the data before performing a fit using the selected function and the maximum order of the polynomial. If "spectrally independent correction" is selected, a weighted average of the fit parameters along the spectral dimension is performed similarly to the iterative case. The data are corrected by subtracting the fit or dividing by the fit depending if the function chosen was the polynomial or the exp(polynomial). If "Exclude zero-th order term from fit" is selected the zero-th order term of the polynomial is not included in the calculation of the background before correction. This results in a background equal to 0 (1) at the center of the image for the polynomial (exp(polynomial)) function, respectively.

This method is the preferred option if there is not a clear spatial separation between background and foreground, but still a clear contrast between background and foreground.

HIA 2.0, 03/03/2015, Francesco Masia


Insert background parameters

This function allows to define the background function parameters to correct the data for an uneven background.

Background function parameters
  • Function: Defines the function to use for the fit: either polynomial or exp(polynomial).
  • Polynomial order: Defines the order of the polynomial.
  • Ok: Executes the function

The selected function and polynomial order determines the layout of the "function parameters" panel in which the parameters are entered. The function uses the voxel size defined in the Scale panel and assumes the origin in the center of the image. Once the background has been calculated, the data are corrected by subtracting or dividing by the background depending if the function chosen was the polynomial or the exp(polynomial). The correction is spectrally independent.

HIA 2.0, 04/03/2015, Francesco Masia


Insert background parameters

Load background parameters

This function allows to load the background function parameters from a file to correct the data for an uneven background. A compatible parameter file is created by saving the parameters determined by other function within the background correction. After loading, the present data is corrected.

  • Manual
  • Automatic (Iterative)
  • Automatic (largest error)
  • Fit parameters (smallest error)

Once the parameters have been loaded, the panel below pops-up:

Load parameters
  • Exclude zero-th order term from fit for correction: defines if the zero-th order of the polynomial is excluded from the correction of the background.
  • Exclude zero-th order term from fit for correction (real part): (Visible only for complex susceptibility data) defines if the zero-th order of the polynomial is excluded from the correction of the background for the real part.
  • Spectrally independent correction: defines if the background correction is spectrally dependent or not.
  • offset x,y,z: defines the offset between the origin of the system used for the determination of the parameters of the background function and the center of the image. The function uses the voxel size and unit defined in the Scale panel.
  • Ok: Execute the function

If the spectrally independent correction is selected, a weighted average of the parameters along the spectral dimension is used. Otherwise the loaded file must contain the parameters defined in a spectral range which contains the spectral range off the present data. An interpolation of the parameters is performed if the spectral positions do not match. The present data are corrected by subtracting or dividing by the fit depending if the function defined in the file is the polynomial or the exp(polynomial). If the "Exclude zero-th order term from fit" option is active the zero-th order term of the polynomial is not included in the calculation of the background before correction.

HIA 2.0, 04/03/2015, Francesco Masia


Fit background parameters (smallest error)

This function finds the background function parameters by fitting the data at the different spectral positions with no additional masking (the mask defined in the exclude regions function is retained). The spectral position resulting in the smallest fir error is used to determine the fit parameters which are then applied to all spectral positions. The function uses the voxel size defined in the Scale panel and uses the origin in the center of the image.

Fit background parameters (smallest error)
  • Function: Defines the function to use for the fit: either polynomial or exp(polynomial).
  • Polynomial order: Defines the order of the polynomial.
  • Exclude zero-th order term from fit for correction: defines if the zero-th order of the polynomial is excluded from the correction of the background.
  • Exclude zero-th order term from fit for correction (real part): (Visible only for complex susceptibility data) defines if the zero-th order of the polynomial is excluded from the correction of the background for the real part.
  • Save fit: saves a tiff image with the background obtained from the fit.
  • Save fit parameters: save an ascii file with the spectrally dependent fit parameters.
  • Ok: Execute function.

The data are corrected by subtracting or dividing by the fit depending if the function chosen was the polynomial or the exp(polynomial). If the "Exclude zero-th order term from fit" option is active the zero-th order term of the polynomial is not included in the calculation of the background before correction. This reflects to a background equals to 0 (1) at the center of the image for the polynomial (exp(polynomial)) function, respectively. The correction is spectrally independent.

If the "Save fit" and/or the "Save fit parameters" option are active, the user should give the name(s) of the file(s) where to store the background fit/parameters.

HIA 2.0, 04/03/2015, Francesco Masia


CARS ratio calculation

This function calculates the ratio between the CARS intensity of the present data and the intensity measured in a non-resonant medium, e.g. glass for wavenumbers >800/cm.

There are two options:

  1. Importing directly a spectrum measured on a non-resonant medium. An ascii (.dat or .txt) file is imported. The first column has to be the spectral position, the second column has to be the CARS intensity measured in the same conditions as the data. The range of the spectrum must contain the spectral range of the present data. If the spectral positions are not the same, the function will interpolate the imported spectrum.
  2. Importing the ratio of the spectrum of water (or another pure medium present in the image) to a non-resonant medium. The user must select the region of the data containing the pure medium (water). The selection works similarly to the Exclude regions function. The area in blue represent the spatial position on which an averaged water spectrum is calculated. Excluded regions (in red) will not be included in the calculation. The non-resonant medium will then be calculated as the product of the imported ratio and the averaged water spectrum. The range of the imported ratio must contain the data spectral range. If the spectral positions are not the same, the function will interpolate the imported spectrum.
CARS ratio import settings
  • Import:
    • Non resonant medium spectrum: import directly a spectrum measured on a non-resonant medium.
    • Non resonant medium- water ratio: import the ratio between the spectrum of a non-resonant medium and water.
  • Minimum CARS Ratio: this defines the minimum value for the CARS ratio. CARS ratio must be positive to be used in the field retrieval.
  • Ok: Imports the file and performs the ratio calculation.
The user can decide to overwrite the selected dataset or to create a new one with the modified data.
HIA 2.0, 26/02/2015, Francesco Masia


Singular value decomposition

Singular value decomposition (SVD) is a factorization operation which decompose any S x P matrix D as D=(U Σ V'), where U is a S x S unitary matrix, Σ is a S x S diagonal matrix and V' (conjugate transpose of V) is a S x P unitary matrix. The columns of the matrix U represent the relevant "spectra" which are the new basis. The elements on the diagonal of Σ are the singular values of D and give the "importance" of the new "spectra", which is their average weight in the data. Random Gaussian noise has a nearly linear dependence of Σi,i as a function of i for large value of i. The columns of V are the new coefficients to describe the spectrum of each spatial voxel in terms of the new basis. The SVD operation can be used for data denoising, compression and masking.

The function performs a SVD on a matrix D=sgn(d)α (dβ+γ)δ.

Choose α, β, γ, and δ to have data with white noise. By default, α=0, β=1, and δ=0.5 for CARS Intensity and Raman data, and α=0, β=1, and δ=1 for CARS field (after field retrieval) and SRS data. γ is set to the square of the random noise value defined when loading ascii data. This SVD filtering is designed for real data D. After decomposition, the diagonal elements of Σ are displayed as a function of the index order. The data denoising/compressing is obtained by limiting the number of elements of Σ to be used in the reconstruction D=(U Σ V'), where Σ is a square matrix as Σ setting the excluded diagonal elements with 0.

In the plot of the diagonal elements of Σ, the elements which will be excluded in the reconstruction are shown in blue, while the included elements will appear red. There are two ways to define the elements to include in the reconstruction.

  1. Threshold method: First fits Σi,i with a linear function for a range of i associated with noise (typically for i>S/2, where s is the total number of spectral points in the image). Elements which are smaller than the fit times a factor (user settable, typically sqrt(2)) are excluded from the reconstruction. The fit will appear as a blue dashed line, while factor*fit will appear as a red dashed line.
  2. Maximum index method: Manual selection of a maximum index imax for Σi,i.
SVD denoising

  • Fit range: Defines the range for the linear fit of the singular values.
  • Threshold (factor*fit): activate the threshold method.
  • factor: Defines the factor for the threshold method. The smallest index which satisfies the threshold will be shown in the "maximum index" field. If the maximum index method is selected this field is not active and it indicates the value of the ratio between the fit value and the singular value at this index.
  • Maximum index: activates the maximum index method.
  • index: Defines the index for the maximum index method. The ratio between the fit value and the singular value at this index is shown in the "factor" field. If the threshold method is selected this field is not active and it indicates the value of the smallest index which satisfy the threshold relation.
  • Manual: activates the manual method.
  • Display: Displays the SVD matrices (U and V)

    • include: (Active if the manual selection is active) Indicates if the shown component is included in the denoised reconstruction
  • Save (only for data denoising): thick the box to save the U, Σi,i and V matrices. Only the indeces considered for the reconstruction will be saved. It saves three files: two ascii files (extension .dat) for the U and Σi,i (indicated by _U and _S, respectively in the name file) and a tif image for V (8 or 16bit). The file 'name_U.dat' has the spectral position in the first column and the SVD spectra in the other columns ordered by the value of the singular values. The file 'name_S.dat' is a two columns file with the first column the index position and the second column the corresponding singular values.
  • Ok: accepts the settings and executes the data reconstruction.

The effect of SVD denoising on CARS data is shown below. The spectrum of the original data (Data1) is shown in red, while the spectrum of the SVD denoised data (Data2) is shown in cyan.

SVD spectra comparison

The user can decide to overwrite the selected dataset or to create a new one with the modified data.
HIA 2.3, 23/03/2016, Francesco Masia


Field retrieval (PCKK)

This function allows to retireve the CARS susceptibility from CARS intensity measurements. It works correctly if the actual data are the CARS ratio (see CARS ratio calculation). The phase of the CARS susceptibility is calculated as φ=Im[FT[G(t)*FT-1[ln(Ratio)]]], where FT is the Fourier transform operator, G(t) is a product between a heavyside function and a Gaussian with the width defining a time domain filtering. The width should be at least equal to the pulse duration to avoid broadening in the spectral resolution.
To correct for the error in the phase (φ0) , one can either calculate a constant phase offset using the formula above with a time domain filtering time constant much shorter than the pulse duration, or to fit the phase with a weight decreasing as the residual increases similarly to what described for the Background subtraction function. The function returns the imaginary part of the CARS susceptibility Im(χ)=(Ratio)1/2sin(φ-φ0), and a spectrally averaged real part Re(χ)=<(Ratio)1/2cos(φ-φ0)>.

Filed retrieval settings
  • Time domain filter:
    • Amplitude: tick the box to calculate the amplitude of the CARS susceptibility using a time domain filtered phase
    • Phase: thick the box to filter he phase in the time domain
    • Time domain filter width (active if one the time domain filter correction is active): defines the width of the Gaussian G(t), in ps.
  • Phase offset correction:
    • None: no correction of the phase offset will be included;
    • Constant: a constant phase offset φ0 will be calculated and removed;
    • Robust: a spectrally dependent phase error will be calculated using a robust fit approach. It will open a panel to insert the parameters for the fit. See Background subtraction
  • Phase correction time domain filter width (active if constante phase correction is active): defines the width of the Gaussian G(t) to calculate φ0, in ps.
  • Binnig for phase correction (active if phase correction is active):
    • None: no binning
    • Phase offset: φ0 is binned before remove from φ in the calculation
    • Data: Data are binned in the calculation of φ0
    • The user will be prompted to insert the binning size if selected.
The user can decide to overwrite the selected dataset or to create a new one with the modified data
The user can decide to keep the phase information. In that case, other two datasets will be created. The first one is the phase calculated before correction, the second one is the phase offset (can be wavelength dependent if the Robust option was used)
HIA 2.4, 06/12/2016, Francesco Masia


FSC3

This function factorizes the data D into susceptibilities and concentrations of chemical components. The iterative algorithm finds a set of concentrations C and susceptibilities S which minimizes the error in the factorization, i.e. the Frobenious norm of D-C*S is minimized, with nonnegativity constrain on C and S. For complex susceptibility data, there is the option to include the real part in the calculation, i.e. D={Im(χ), f*sqrt(S)*Re(χ)}, where S is the number of spectral points. Also for complex susceptibility data, it is possible to determine the matrix C such as the sum of the concentrations at any spatial point is as close a possible 1. This is obtained by minimizing the sum concentration error, i.e the norm of (ΣC)-1. The quality of the factorization is measured in terms of spectral error ES=rms(D-C*S)/||D||F*sqrt(S*P), where ||.||F indicates the Frobenius norm and P is the total number of spatial positions, and in the concentration error EC= ΣC-1.

More details can be found in the publication Masia et al, J. Raman Spectroscopy (2015)

Single FSC3 step

The FSC3 algorithm is based on a non-negative matrix factorization using either a block-pivot or active-set method. This iterative method is characterized by a minimum number of iterations, a target stop criterion, a maximum number of iterations and a maximum computation time. Random initial guesses for C and S are used. The stop criterion can be either the total error or the tolerance.

  • The total error is defined as ||D-C*S||F.
  • The tolerance measures the reduction in the norm of the projected gradients of the solution with respect to the initial step. It is defined as ||{∇Ci∇SiT}||F/||{∇C0,∇S0T}||F*n0/ni, where ∇C=(C*S-D)*ST, and ∇C=CT*(C*S-D). The subscript 0 indicates the values of the gradient calculated using the initial concentrations and spectra, and n0 is the number of elements of the matrix {∇C0,∇S0T}. The subscript i refer to the step in the iteration, and in the calculation of ||{∇Ci∇SiT}||F only the elements of the matrices ∇X (X=C,S) corresponding to either negative ∇X or positive X are included. ni is the number of elements of the resulting matrix.
The output solution is the one with the smallest stop criterion. In order to improve the reproducibility of the results there are two possibilities:

  • High/Low target method: the FSC3 is performed Nh times using random initial guesses with a target stop criterion Th. The solution which gives the minimum value of the selection criterion is used as initial guesses for a FSC3 with target stop criterion Tl.
  • Knock-out method: 2Nko calculations with target stop criterion T are performed with random initial guesses.The 2Nko-1 best solutions (see Automatic selection criterion in 'high/low' and 'knock-out' methods) are considered as initial guesses for the next round of FSC3 calculations, again with target stop criterion T. The algorithm continues until the final step. There is the option to check for dissimilarity between the last two solutions. In this case if the difference in the error exceeds a defined value a new set of 2Nko FSC3 calculations are launched and once the last two solutions are obtained, they will compared with the other two solution from the previous set of calculations. The user can define the maximum number of checks to perform before accepting the results.

Weighted FSC3 algorithm

The single FSC3 calculations return a result which minimizes the error globally. In some datasets this means that some spatial regions are factorized better than others. In order to obtain a more homogenous spectral error, the data can be weighted by a factor which is related to the spectral and concentration errors. The weighted FSC3 (wFSC3) is an iterative algorithm where the weight at the step i+1 is calculated as wi+1=(wi)(1-α)*(ES/<ES>)(1-α')*e(-β*|EC|), where <.> indicates the average value. The convergence is reached if the quantity |wi+1/<wi+1>-wi/<wi>| is not improved for three consecutive iteration steps or the maximum number of calculations is reached. If the weight exceeds the threshold γ*P^(1/4), the data at this position is excluded from the next step in the iteration.

Number of chemical components

The number of chemical components can be either entered manually or determined automatically by the algorithm. In the latter case, the FSC3 or wFSC3 is performed with a defined initial number of components and when convergence is reached the calculation is repeated for a number of components increased by 1. The change in the number of components is accepted if the improvement in the error from the previous step (defined as the difference in the norm of the spectral error divided by sqrt(P)) is better of a certain factor η*std(ES), where std(.) indicates the standard deviation. Otherwise the calculation is repeated for a number of components decreased by 1. The change in the number of components is accepted if the increase in the error from the previous step is below the threshold. If a change in the number of components is accepted the algorithm continues, otherwise stops. In case the maximum number of component is reached, the algorithm stops.

The workflow of the FSC3 algorithm is shown below for a complex susceptibility data, with included real part (factor f=1) and automatic determination of the number of components. K indicates the initial number of components

FSC<sup>3</sup> workflow

The parameters for the FSC3 calculation are defined in the panel below:

  • Load FSC3 parameters: Loads the parameters for the FSC3 panel from a .par file.
  • Save FSC3 parameters: Saves the parameters for the FSC3 panel to a .par file.
  • Reset FSC3 parameters: Reset the parameters to the default values.
  • Data selection: Defines the data to be analysed by entering their indexes, e.g. "1 4 6 8"
  • Use real part: Defines if the spectrally averaged real part is included in the calculation. Active if all selected data are complex susceptibilities.
    • Real part factor: Defines the factor f.
  • Single FSC3 settings
    • Number of components: Defines the number of components used in the factorization. This option is active if the number of components is fixed.
    • Min # iterations: minimum number of iterations in the single factorization function.
    • Max # iterations: maximum number of iterations allowed in the single factorization function.
    • Max time (s): maximum time in seconds (Stops the calculation after this time, use e.g. 1000).
    • Solver: Choose the NMF solver to be used - Choices are BP (block-pivot) or AS (active set). BP is generally faster.
    • Thermal annealing: Defines if additional noise is introduced into the calculated spectra and concentration during the factorization algorithm. This option is available only for single FSC3 method. The new spectra and concentration at the end of the i-th iteration are calculated as Si=Si+ΔS*(D-Ci*Si)/(Ci*ΔS)*t*qi and Ci=Ci+ΔC*(D-Ci*Si)/(ΔC*Si)*t*qi, where ΔS and ΔC are random generated matrices of appropriate size.
      • Thermal annealing factor: Defines parameter t.
      • Thermal annealing quenching: Defines parameter q.
    • Concentration constrain: Defines if the concentrations will be normalized to minimize their RMS deviation of their sum from unity over the spatial pixels.
  • Limit spatial points: Limits the number of spatial points used in the non-negative matrix factorization by excluding a fraction of the available points For the excluded points (e.g. in the Excluded regions function) the function will project the spectra to obtain the concentrations, similarly to the Project spectra function, using the selected non-negative constrain method. There are three options: None (all points included), Equidistant (points on a regular grid are included), Random (randomly selected points are included)
    • Fraction to consider: determines the fraction of points to include for the Equidistant and Random options. For the Equidistant option, three fractions corresponding along the three spatial directions are required, while one fraction is requested for the Random option.
  • Guided analysis: If this option is selected, the user is prompted to import a file with spectra to be used as guides. The file must have a dat extension and have the first column with the wavenumbers and the following columns with the spectra. If real part is used in the analysis, the user is asked to import a file with the real part of the spectra. In this case, the file must have a dat extension and the real parts should be arranged in a single row. The number of real part components must match the number of the imaginary part spectra. Once the spectra are uploaded, the panel below gives the option to create the guide spectra using linear combinations of the imported spectra. The new spectra will be included win the FSC3 analysis with a user-defined weight.

    • Load Spectra parameters: Loads the parameters for the guide spectra panel from a .par file.
    • Save Spectra parameters: Saves the parameters for the guide spectra panel to a .par file.
    • Reset Spectra parameters: Reset the parameters to the default values.
    • NSi: Indicate the New Spectra to be used in the guided analysis obtained as linear combinations of the imported Pure Spectra (PSi). The sum of the coefficients will be normalized to 1.
    • Weight: indicate the weight to be used in the guided analysis. It is in percentage of the total number of spatial point of the images analysed with the FSC3.
    • Concentration: (Only if the Concentration constrain option is selected) indicates the concentration corresponding to the new spectrum used for the guided analysis.
    • Include: indicates if the spectrum is included or not in the guided analysis. The color of the label corresponds to the color of the line in the plot.
  • Target stop criterion
    • Method: determines the FSC3 calculation mode between the options:
      • Single: a single FSC3 is performed at each step
      • High/Low: In a first step a set of FSC3 with high target stop criterion are calculated using random initial values. In a second step a single FSC3 calculation with low stop criterion target using the high target result with the smallest error as initial value, unless the option for the manual selection of the initial value is active.
      • Knock-out: a series of FSC3 are calculated using random initial values. The results are sorted by the error and the half with the lower error are taken as initial values for the next iteration, unless the option for the manual selection of the initial value is active. This is repeated until a single result if reached.
    • Method: defines what stop criterion to consider in the iteration between the options
      • Tolerance: the stop criterion is given by the tolerance reached in the iteration. The tolerance is the change of the concentration and spectra per iteration relative to the first iteration (FM, add the formulas). It quantifies the relative speed of convergence. In a well-behaved convergence, this “speed” is reducing with the “distance” to the minimum.
      • Error: the stop criterion is given by the absolute error of the factorization.
    • Target: target stop criterion for the Single and Knock-out method.
      • High/Low method
        • High target: target stop criterion for the first step.
        • Low target: target stop criterion for the second step.
        • Successful # high target calculations: target number of successful high target calculations. If a calculation shows singularity (which can occur depending on the initial values and the dataset), the calculation is repeated with different random initial values until the number of convergent calculations or the maximum number of calculations is reached.
        • Maximum # high target calculations: defines the maximum number of high target calculations which are performed before calculating the low target FSC3.
        • Display after high target steps: Defines if the results are displayed after the high target round. It allows to select manually the solution to consider as initial value for the low target calculation. This option is active only if the weight and the number of components are not changed during iteration. See Manual selection of intermediate solutions in "high/low" and "knock-out" methods for more details.
      • Knock-out method
        • Starting successful # calculations: defines the number of stages Nko in the knock-out method. This requires in a number of 2Nko sucessfull calculations in the first stage.
        • Maximum # calculations per stage: defines the maximum number of calculations which are performed per stage in the knock-out method. This has tio be larger than 2Nko.
        • Display after # steps: defines at which stages the results are displayed. Enter Nko+1 to display only the final results. It allows to select manually the solutions to consider as initial values for the next round. This option is active only if the weight and the number of components are not changed during iteration. See Manual selection of intermediate solutions in "high/low" and "knock-out" methods for more details.
        • Check similarity at final stage: select this option to enable the comparison of the errors of the last two solutions during the knock-out method (#stage = Nko).
          • Maximum relative error: defines the threshold in the relative error above which results are considered to be dissimilar. If they are dissimilar, an additional knock-out round starting with 2Nko calculations is performed and the last two results will be added to the results from the first run.
          • Maximum # checks: defines the maximum number of checks for similarity. If this number of checks is reached, the calculation proceeds with the last stage independent of the similarity.
    • # iterations defining divergence: Defines the number of iterations which identify a divergence. During the iteration, if the stop criterion does not improve within the defined number of steps, the iteration is stopped and the solution with the minimum stop criterion is returned.
    • Save convergence data: saves the stop criterion vs iteration step in ascii format. For High/Low and Knock-out method the stop criterion will be saved at each round. For the High/Low method, the output file will have columns labeled in the format ht-i, where i is the index of the high target calculations and the last column labeled as lt-j, where j is the index of the high target calculation used as initial solution for the low target round. For the Knock-out method, the output file will have columns labeled in the format r-t-i, where is the knock-out round number and i is the index of the knock-out calculation.
    • Component orthogonality factor: defines the factor ξ which defines the weight of the component orthogonality factor in the automatic selection criterion for the "high/low" and "knock-out" methods. See Automatic selection criterion in 'high/low' and 'knock-out' methods for more details.
  • weighted FSC3
    • Change weight: activates the weighted FSC3 routine
      • alpha: defines the α parameter for the calculation of the weight.
      • alpha': defines the α' parameter for the calculation of the weight.
      • beta: defines the β parameter for the calculation of the weight.
      • gamma: defines the γ parameter to define the threshold of the weight.
      • Minimum number of weighted calculations: defines the minimum number of calculation before the function checks for convergence.
      • Maximum number of weighted calculations: defines the maximum number of calculation before the function returns the results.
      • Use previous step results as initial guesses: Defines if the results from the previous step in the weighted algorithm are used as initial guesses. After the first step the 'single' method is used. Otherwise random initial guesses are used.
    • Weight from FSC3 components: activates the weighting routine using a FSC3 component previously calculated. After a first run of FSC3 without weighting, the user is prompted after closing the FSC3 panel to indicate a component and a power to weight the data for the next round. The data will be weighted according to the FSC3 concentration elevated to the indicated power. After each FSC3 results the user has the option to continue the weighting with a newly calculated FSC3 component.
  • Determination of number of components
      • Change number of components: activates the automatic determination of the number of components.
      • Initial number of components: defines the initial number of components used.
      • Maximum number of components: defines the maximum number of components used.
      • eta: defines the η factor for the threshold in the spectral error norm difference.
      • Display FSC3 results at change of # components: activates the display of the FSC3 results every time the number of components are changed.
  • Spatial correlation: Defines if the concentrations maps C and data D are smoothed for the calculation of the spectra during the iterations. The smoothing is performed with a Gaussian filter. The values of the kernel elements are g(i,j,k)=exp(-(i/sigma_x)2)+exp(-(j/sigma_y)2)+exp(-(z/sigma_z)2), with the central pixel of the kernel is defined as (i,j,k) = (0,0,0). The sizes of the kernel at each directions are adjusted to be three times the widths of the Gaussian.
    • sigma x: defines the width along the x axis of the Gaussian filter in units defined in the pixel scaling
    • sigma y: defines the width along the y axis of the Gaussian filter in units defined in the pixel scaling
    • sigma z: defines the width along the z axis of the Gaussian filter in units defined in the pixel scaling
  • Group spectra: activates the spectra grouping routine. In this routine, the similarity of the FSC3 spectra is checked during the iterations. Spectra are considered similar if their cross-correlation is above a certaing threshold. Similar spectra are replaced by a common spectrum centered at their original central spectral position and will be inlcuded in a common group. If the spectral positions of two spectra are too close (below a certain threshold) they will replaced by a common spectrum which will inherit the concentration of both components. When this option is active, the number of components refer to the desired number of groups. The number of groups is kept constant, i.e. if components are grouped together new random components are introduced in the iterations. Additionally to the standard concentration and spectra, this routine gives the option to save
      1. an ASCII file with the following columns
        1. group index
        2. FWHM of autocorrelation of spectrum
        3. component index
        4. Relative energy position
        5. Integrated concentration
      2. an ASCII file with the spectra for each group (if includes more than one component)
      3. an ASCII file with the probability spectra for each group
      4. an ASCII file with the number of components as a function of the NMF iteration step
      5. a TIFF file with the residuals
    • Cross correlation threshold: defines the similarity threshold
    • FHWM factor: defines the proximity threshold as the minimum distance in units of auto-correlation FWHM
    • Additional it.: defines the minimum number of NMF iterations to be considered when extra components are included.
    • Trap model: activates the trap model (work in progress, details to come)
  • Constrain broad features: activates the routine to constrain broad features. In this routine FSC3 spectra which are identified as "broad" will be replaced by their low-frequency filtered spectra. This fucntion is still under development. A spectrum is considered "broad" if the ratio of integrated FFT powers of low-frequency filtered data and of the original data is above a certain threshold
    • Target FFT power threshold: defines the threshold for the assignment of the spectrum as "broad".
    • FFT Gaussian width (ps): defines the width of the Gaussian filter to create the low-frequency filtered data
    • #iterations without constrain: defines the number of free NMF iterations before the routine starts to check the spectra
    • #iterations to reach target threshold: if larger than one the threshold starts at 1 and decrease linearly to reach the target threshold within the defined number of iterations.

Automatic selection criterion in 'high/low' and 'knock-out' methods

The selection of the solutions to be used as initial guesses after the first step of a 'high/low' or 'knock-out' calculation is based on the error committed in the factorization and in the orthogonality of the concentrations of the component. A triangular matrix will be created as M(i,j)=(εi/εj)*(si/sj)ξ, where i<j are the indeces of intermediate solutions, ε is the factorization error and M(1,1) is set to 1. ξ is defined in the settings panel as the component orthogonality factor. s is defined as Σn>m 2/(K*(K-1))*(Cn·Cm)) , with Cn·Cm indicating the scalar product between the concentrations of the components m and n, and K is the number of the factorization components. The algorithm will look for the highest values in the matrix M excluding the corresponding solutions until the correct number of best solutions are left. If only the factorization error has to be considered for the selection choose ξ=0 and the algorithm will select the solutions with the smallest ε.

Manual selection of intermediate solutions in "high/low" and "knock-out" methods

If the option "Display after high target steps" is selected in a 'high/low' calculation, or the "Display after # steps" is less than Nko+1 in a 'knock-out' run, the intermediate solutions are shown. The display is similar to the result display with a selection tick box on the bottom right. The solutions which are selected according to the automatic selection criterion to be the initial guesses for the next step have the box ticked, while for the solutions which are not selected the box is dimmed. The panel shown below is also displayed together with the results

  • Continue calculation for # steps: (only for 'knock-out' calculations) defines the number of steps which the 'knock-out' calculation will run before displaying the solutions again.
  • Continue: Continues the calculations with the selected solutions. Active when the correct number of solutions are selected.
  • Reset: Resets the selected solutions using the automatic selection criterion.
  • Stop: Stops the calculations
To manually select the solutions to be used in the next step, untick the selection boxes from the displays of the solutions which shouldn't be considered. The selection boxes of the solutions initially not selected are undimmed and can be selected. Select the correct number of solutions (1 for the 'high/low' method and half of the displayed solutions for the 'knock-out' method) and continue.

Result display

At the end of the FSC3 algorithm the resulting concentration maps, susceptibility spectra, spectral and concentration errors are displayed. If more than one image is considered, they will be displayed in different panel. The concentration maps will be displayed with grayscales ranging from the minimum to the maximum of the individual component. If more than one image is considered there is the option to display the concentration maps on a grayscale with a global range (the same for the different images) or fitting to the range of concentrations in each image. If the displayed solution was obtained from a diverging iteration, the label are shown in red.

  • Save FSC3 results: save the concentration, spectral error and concentration error maps as tif images (hyper-stacks) and the spectra as ascii files (the real part will be saved in a separate file). If some regions have been excluded from the calculations there is the option to save the tif images with a grayscale which spans over the maximum range of concentrations (Consider image for range),i.e. including the concentrations obtained from projecting the calculated spectra on the excluded regions, or to the range of concentrations which consider only the regions included in the factorization (Consider ROI for range).
  • Save FSC3 results as Powerpoint: Saves the FSC3 panel result as a Powerpoint slide. If the selected file already exists, the function appends the new slides to the end of the file.

Weighted FSC3 displays

At the end of a weighted FSC3 calculation the convergence and the final weight are displayed.

  • Weighted FSC3 convergence:
    • Save weighted FSC3 convergence: The dependence of the quantity |wi+1/<wi+1>-wi/<wi>| on the iteration step will be saved as an ASCII file
  • Weighted FSC3 results:
    • Save weight: The weight corresponding to the step where the quantity |wi+1/<wi+1>-wi/<wi>| is minimum will be saved as a tiff image.
HIA 2.3, 24/03/2016, Francesco Masia


Project spectra

This function allows to find the distribution maps of known spectra in the data. The function solves the system D=X*S, where D are the data , X are the distribution maps and S are the spectra. The user has to import an ascii file with the spectra. The file consists of a series of tab-separated columns where the first column gives the spectral positions, and the others are the spectra of the chemical components which are used for the projection.

In case the data are complex susceptibility, the real part can be included (optional) in the matrix S with a factor f*sqrt(S) where S are the number of spectral points in the data. The factor f can be entered in the input dialog box shown below

Real part weight

The real part of the spectra needs to be imported as an ascii file, with the spectrally-averaged real parts of the different components tab-separated.

The function solves the system (D*ST)=X*(S*ST) using the mldivide MATLAB function. The resulting maps X can have negative values. The user is given the option create a non-negative X by using the solution X of the system as starting values for a non-negativity constrained least squares algorithm (same as used in the FSC3) for X.

Nonnegativity constrain
  • None: no additional negativity constrain
  • BP: nonnegativity constrain using the block pivot algorithm
  • AS: nonnegativity constrain using the active set method
If the fucntion was called with the Concentration constrain option, the resulting maps X and spectra S will be rescaled to minimize the concentration error, i.e. the norm of (ΣC)-1.

The resulting maps X will be displayed similarly to the FSC3 results.

HIA 2.3, 23/03/2016, Francesco Masia


Background subtraction (robust fit)

This function allows to subtract a spectrally dependent background using a robust fit approach. It means that the spectrum is fitted from either the top or the bottom. This is achieved by using a weighted least squared fitting algorithm where the weigh of the spectral points depends on the sign and amplitude of the residuals. Weighted least-squares regression minimizes the error estimate s=Σj=1n wj(dj - fj)2, where dj are the data at the spectral point j, fj is the estimated value and wj are the weights. n is the number of spectral points. The weight is adjusted iteratively until the error is smaller than the target tolerance.

The spectrum at each voxel is fitted using either a Fourier series f(x)=a0+Σi=1k ai*sin(2*i*π*(x-xs+N*δ)/(2*N*δ+xe-xs))+ai+k*cos(2*i*π*(x-xs+N*δ)/(2*N*δ+xe-xs)) or a polynomial f(x)=Σi=0n ai*xi. xe and xs are the final and initial spectral points, respectively. k is the number of Fourier terms included in the fitting function. N defines the extension of the spectral range of the Fourier series in units of the number of the period of the highest harmonics (2*δ). δ is the minimum width of the spectral features that will be considered as background. δ is calculated as (xe-xs)/(2*k - 2*N). N must be smaller than k. n is the polynomial order. If the Fourier series is selected, the initial values for the coefficients ai are calculated as the Fourier terms. If the polynomial function is selected, an initial fit is performed on the first pixel using all coefficients equal to zero except a0 which is set to the minimum or maximum value of the spectrum depdening if σth>0 or σth<0, respectively.

In the algorithm a first fit with constant weight is performed and the sign and amplitude of the residuals are used to determine the weights (rj=(dj - fj): w(rj>0) = e-(rj/σ)2 and w(rj<=0) = 1 if σ>=0, or w(rj<=0) = e-(rj/σ)2 and w(rj>0) = 1 if σ<0). σ is calculated as the sign(σth)*rms(r). If σ>σth, a weighted fit is then performed and a new σ is calculated as sign(σth)*rms(r*w(r)). If still σ>σth, a new weighted fit is performed using the new σ to calculate the weights. The loop stops if either 3 iterations are performed or σ<=σth. A final weighted fit is performed with weights calculated using σth.

To improve the speed of the function there are the following options:

  1. the spectrum can be binned in order to reduce the number of spectral points
  2. the initial coefficients can be estimated from a coarse spectrum (for the Fourier series option)
  3. the fit results obtained from adjacent pixels can be used as starting values.

Background subtraction (robust fit)

  • Load parameters: Loads the parameters for the robust fit panel from a .par file.
  • Save parameters: Saves the parameters for the robust fit panel to a .par file.
  • Reset parameters: Reset the parameters to the default values.
  • Function: defines which function to use for the fit. There are two options
    • Fourier series
    • Polynomial
  • Number of harmonics: (Active only if the function is Fourier series) defines the number of harmonics (k) to be used in the Fourier series
  • Extension (period units): (Active only if the function is Fourier series) defines the extension of the spectral range of the fitting function in units of the period of the highest harmonic components. N must be smaller than k. If N>=k, it will chnaged automatically to k-0.5.
  • Min/Max width of spectral fetaures: displays the maximum and minimum width (δ) of the spectral features that will be considered as background according to the selected number of harmonics and extension.
  • Initial coefficients from coarse data: If this option is selected, a coarse spectrum (at 2*k+1 equidistant spectral points between xs-N*δ and xe+N*δ vith values equal to the values of the spectrum at those position + σth) is used to calculate the Fourier terms.
  • Polynomial order: (Active only if the function is Polynomial) defines the polynomial order (n).
  • Sigma: defines the target Gaussian width used in the calculation of the weigth (σth).
  • Max # Iter: defines the maximum number of iterations for the fit.
  • Termination tolerance: defines the tolerance required before accepting the fit result.
  • Use coefficient from adjacent pixel: If this option is selected the function uses as starting values the coefficients obtained from the fit of the prevoius voxel.
  • Binning: defines if the spectra have to be binned before starting the fit.
  • Bins: defines the number of the original spectral points to be binned. The suggested value is the rounded square root of a fifth of the total number of spectral points.
  • Mask spectral regions: If this option is selected the user can remove spectral regions from the data to fit. A panel similar to Multiple spectral selection function opens
  • Ok: accept the selection and start the .
The function will return two new datasets: one with the background subtracted data and one with the background itself.

HIA 2.3, 02/11/2015, Francesco Masia


Third part licences

HIA contains third party code subject to the following BSD licences

saveppt2

Ver 2.2, Copyright 2005, Mark W. Brown, mwbrown@ieee.org

Ver 3.5, Copyright 2010, Jed Frey, frey_jed@cat.com

All rights reserved.

Redistribution and use in source and binary forms, with or without

modification, are permitted provided that the following conditions are

met:

    * Redistributions of source code must retain the above copyright

      notice, this list of conditions and the following disclaimer.

    * Redistributions in binary form must reproduce the above copyright

      notice, this list of conditions and the following disclaimer in

      the documentation and/or other materials provided with the distribution

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"

AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE

IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE

ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE

LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR

CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF

SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS

INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN

CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)

ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE

POSSIBILITY OF SUCH DAMAGE.

polyfitn

Copyright (c) 2012, John D'Errico

All rights reserved.

Redistribution and use in source and binary forms, with or without

modification, are permitted provided that the following conditions are

met:

    * Redistributions of source code must retain the above copyright

      notice, this list of conditions and the following disclaimer.

    * Redistributions in binary form must reproduce the above copyright

      notice, this list of conditions and the following disclaimer in

      the documentation and/or other materials provided with the distribution

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"

AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE

IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE

ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE

LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR

CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF

SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS

INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN

CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)

ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE

POSSIBILITY OF SUCH DAMAGE.

Subaxis - Subplot

Copyright (c) 2014, Aslak Grinsted

All rights reserved.

Redistribution and use in source and binary forms, with or without

modification, are permitted provided that the following conditions are

met:

    * Redistributions of source code must retain the above copyright

      notice, this list of conditions and the following disclaimer.

    * Redistributions in binary form must reproduce the above copyright

      notice, this list of conditions and the following disclaimer in

      the documentation and/or other materials provided with the distribution

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"

AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE

IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE

ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE

LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR

CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF

SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS

INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN

CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)

ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE

POSSIBILITY OF SUCH DAMAGE.

Suptitle

Copyright (c) 2004, Ben Barrowes
All rights reserved.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:

* Redistributions of source code must retain the above copyright notice, this
  list of conditions and the following disclaimer.

* Redistributions in binary form must reproduce the above copyright notice,
  this list of conditions and the following disclaimer in the documentation
  and/or other materials provided with the distribution
* Neither the name of Barrowes Consulting / Dartmouth College / US Army CRREL nor the names of its
  contributors may be used to endorse or promote products derived from this
  software without specific prior written permission.
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE
FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

Non-negative matrix factorization

Copyright 2008-2009 by Jingu Kim and Haesun Park,
School of Computational Science and Engineering,
Georgia Institute of Technology

This code comes with no guarantee or warranty of any kind.

arb_gradient.m is a modification of gradient.m

Copyright 1984-2014 The MathWorks, Inc.

fminsearchvary.m is a modification of fminsearch.m

Copyright 1984-2012 The MathWorks, Inc.

Bioformats functions

OME Bio-Formats package for reading and converting biological file formats.
Copyright (C) 2014 - 2015 Open Microscopy Environment:
- Board of Regents of the University of Wisconsin-Madison
- Glencoe Software, Inc.
- University of Dundee
This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as
published by the Free Software Foundation, either version 2 of the
License, or (at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License along
with this program; if not, write to the Free Software Foundation, Inc.,
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.

parseArgs.m

Copyright (C) 2002-2004, Aslak Grinsted
This software may be used, copied, or redistributed as long as it is not
sold and this copyright notice is reproduced on each copy made. This
routine is provided as is without any express or implied warranties
whatsoever.

logb.m

Copyright (c) 2016, The MathWorks, Inc.
All rights reserved.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:

* Redistributions of source code must retain the above copyright notice, this
  list of conditions and the following disclaimer.

* Redistributions in binary form must reproduce the above copyright notice,
  this list of conditions and the following disclaimer in the documentation
  and/or other materials provided with the distribution.
* In all cases, the software is, and all modifications and derivatives of the
  software shall be, licensed to you solely for use in conjunction with
  MathWorks products and service offerings.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE
FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

fminsearchbnd.m

Copyright (c) 2006, John D'Errico
Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:

* Redistributions of source code must retain the above copyright notice, this
  list of conditions and the following disclaimer.

* Redistributions in binary form must reproduce the above copyright notice,
  this list of conditions and the following disclaimer in the documentation
  and/or other materials provided with the distribution.
* Neither the name of the Arctic Research Centre, Aarhus University, Denmark nor the names of its
  contributors may be used to endorse or promote products derived from this
  software without specific prior written permission.
* Neither the name of the The MathWorks, Inc. nor the names of its
  contributors may be used to endorse or promote products derived from this
  software without specific prior written permission.
* In all cases, the software is, and all modifications and derivatives of the
  software shall be, licensed to you solely for use in conjunction with
  MathWorks products and service offerings.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE
FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

Gram Schmidt normalization is a modified version of mgson.m developed by Mo Chen

Copyright (c) 2016, Mo Chen
All rights reserved.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:

* Redistributions of source code must retain the above copyright notice, this
  list of conditions and the following disclaimer.

* Redistributions in binary form must reproduce the above copyright notice,
  this list of conditions and the following disclaimer in the documentation
  and/or other materials provided with the distribution
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE
FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

kkrebook2

Copyright (c) 2005, Valerio Lucarini
All rights reserved.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:

* Redistributions of source code must retain the above copyright notice, this
  list of conditions and the following disclaimer.

* Redistributions in binary form must reproduce the above copyright notice,
  this list of conditions and the following disclaimer in the documentation
  and/or other materials provided with the distribution
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE
FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

Minilzo library

Copyright (C) 1996 - 2017
Markus F.X.J. Oberhumer
Distributed under the terms of the GNU General Public License (GPL v2+).

HIA 2.5, 14/03/2019, Francesco Masia


Licence

HIA is subject to the following BSD licence

Copyright (c) 2015, Francesco Masia, Wolfgang Langbein, Paola Borri

All rights reserved.

Redistribution and use in source and binary forms, with or without

modification, are permitted provided that the following conditions are

met:

    * Redistributions of source code must retain the above copyright

      notice, this list of conditions and the following disclaimer.

    * Redistributions in binary form must reproduce the above copyright

      notice, this list of conditions and the following disclaimer in

      the documentation and/or other materials provided with the distribution

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"

AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE

IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE

ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE

LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR

CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF

SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS

INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN

CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)

ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE

POSSIBILITY OF SUCH DAMAGE.

HIA 2.0, 12/03/2015, Francesco Masia


References

The details of the implemented algorithms are described in the following papers.

  • F. Masia, A. Glen, P. Stephens, P. Borri, and W. Langbein, Anal. Chem. 85, 10820 (2013).
  • F. Masia, P. Borri, and W. Langbein, Opt. Express 22, 4021 (2014).
  • F. Masia, A.Karuna, P. Borri, and W. Langbein, J. Raman Spect. (2015) accepted for publication

The FSC3 algorithm is based on a non negative matrix factorization implementation written by J. Kim, H. Kim and H. Park.

  • J. Kim and H. Park, SIAM Journal on Scientific Computing (SISC), 33(6), pp. 3261-3281, 2011
  • J. Kim and H. Park, In Proceedings of the 2008 Eighth IEEE International Conference on Data Mining, pp. 353-362, 2008.#
  • H. Kim and H. Park, SIAM Journal on Matrix Analysis and Applications, 30(2):713-730, 2008.
  • H. Kim and H. Park, Bioinformatics, 23-12:1495-1502, 2007.

The original code can be found here

HIA 2.0, 12/03/2015, Francesco Masia


Automatic background region determination (Largest error)

{WL: This is a dublication of the iterative method - combine into one html}

The determination of spatial regions corresponding to background is done automatically, with an iterative algorithm, but only at the spectral positions corresponding to the largest error in the fit.

Automatic background region determination (Iterative)

  • Function: Defines the function to use for the fit: either polynomial or exp(polynomial).
  • Initial polynomial order: Defines the initial order of the polynomial of the iterative algorithm.
  • Final polynomial order: Defines the final order of the polynomial of the iterative algorithm.
  • Data pre-masking: defines if a pre-masking routine is used.
  • Maximum number of iterations: defines the maximum number of iterations allowed.
  • Threshold factor (σ): defines the factor used to discriminate background and foreground.
  • Exclude zero-th order term from fit for correction: defines if the zero-th order of the polynomial is excluded from the correction of the background.
  • Exclude zero-th order term from fit for correction (real part): (Visible only for complex susceptibility data) defines if the zero-th order of the polynomial is excluded from the correction of the background for the real part.
  • Spectral independent correction: defines if the background correction is spectrally dependent or not.
  • Save mask: saves a tiff image with the mask used to define the background.
  • Save fit: saves a tiff image with the background obtained from the fit.
  • Save fit parameters: save an ascii file with the spectrally dependent fit parameters.
  • Ok: starts the function
The aim of the function is to automatically mask regions of the image which are interpreted as foreground, fit the unmasked regions and correct the data by the fit. The function uses the voxel size defined in the Scale panel and assume the origin in the center of the image.

Discrimination between background and foreground

The discrimination between background and foreground will be done at the spectral position which returns the largest error. A multidimensional fit with the chosen function will be made at the different spectral positions using the "Final polynomial order". The spectral positions with the largest error will be considered for the following steps
  • Data pre-masking: If the data pre-masking option is active the function will create an histogram of the data value at the specific spectral position, and it will find the mode m. If more than one mode is found the histogram is calculated with half of the elements until a single mode is found. The deviation of the data D from the mode is calculated as Δ=D-m. Data where the modulus of the deviation |Δ| exceeds σpm*rms(Δ) are masked and not considered in the next step. Rms indicates the root-mean-square deviation. The number of histogram elements and the threshold factor σpm are entered by the user (see panel below).
Pre-masking panel

  • Iterative masking: The non-excluded data are fitted using the defined function with the initial polynomial order, obtaining the F. Data where the modulus of the deviation |D-F| exceeds the rms(D-F) by a factor larger than the threshold parameter σ are masked and not included in the next step. The non-excluded data are again fit and a new F and masked regions are calculated. The loop stops when the masked regions don't change between to successive steps or when the maximum number of iterations is reached. The order of the polynomial is increased by one and the routine is repeated until the final polynomial order is reached.

Fitting the background and correcting the data

The final mask will be then used to fit the background at the different spectral positions. If the spectral independent correction option is active, a weighted average of the parameters along the spectral dimension is performed. The weight is defined by the inverse of the error of the fit at the corresponding spectral position. The data are finally corrected by subtracting or dividing by the fit depending if the function chosen was the polynomial or the exp(polynomial). If the "Exclude zero-th order term from fit" option is active the zero-th order term of the polynomial is not included in the calculation of the background before correction. This reflects to a background equals to 0 (1) at the center of the image for the polynomial (exp(polynomial)) function, respectively.

If the "Save mask", "Save fit" and/or the "Save fit parameters" option are active, the user must give the name(s) of the file(s) where to store the mask/background image/parameters.

This function is the preferred option if there is not a clear spatial separation between background and foreground, but still a clear contrast between background and foreground.

HIA 2.0, 04/03/2015, Francesco Masia