Restoration

Welcome to the main page of this documentation. The inpystem library is nothing else than a pluggin to HyperSpy to allow reconstruction.

Some words about pre and post-processing steps

For 2D data

2D data are centered and normalized before reconstruction. The reason is to avoid highly variable reconstruction methods parameters.

# data numpy array is ready to be reconstructed.
#

# Let's get data mean and standard devition.
data_mean, data_std = data.mean(), data.std()

# Let's center and normalize the data
data_ready = (data - data_mean) / data_std

#
# Reconstruction is performed
#

# Let's perform the inverse transformation
data_out = reconstructed_data * data_std + data_mean

# data_out is returned

For 3D data

3D data have the save pre and post-processing as for 2D data. But in addition to centering and normalization, a thresholded principal component analysis (PCA) is performed to reduce the data dimension along the bands axis and to ensure the low rank assumption. This one states that most multi-band data result in the mixing of R basic data with R small compared to the data dimension).

The default behavior is to perform PCA with an automatically estimated threshold (which can be really over-estimated in case of data starvation situations, i.e. if you have almost as many samples as the data size). Though, the user can set both parameters and choose if PCA should be performed and which value the threshold should have.

To sum up, the steps are:

  • perform thresholded PCA if required,

  • center and normalize the data,

  • perform reconstruction,

  • re-set the data mean and standard deviation values,

  • perform inverse PCA.

How to reconstruct my data

To reconstruct the data, the user should use the restore() method. Both Stem2D and Stem3D classes need the following arguments to restore their data:

  • method which is the method name (default is 'interpolation'),

  • parameters which should be a dictionary with input parameters,

  • verbose which allows the method to display information on the console (default is True).

In addition to these common parameters, the Stem3D class has the foolowing inputs:

  • PCA_transform which controlls the PCA execution (defualt is True for PCA execution),

  • PCA_th which states the PCA threshold.

A common reconstruction task will then look like this.

>>> data = inpystem.load_example('HR-sample', ndim=2, scan_ratio=0.2)
Reading configuration file ...
Generating data ...
Creating STEM acquisition...
Correcting STEM acquisition...

>>> rec, info = data.restore()
Restoring the 2D STEM acquisition...
-- Interpolation reconstruction algorithm --
Done in 0.01s.
---
>>> rec
<Signal2D, title: HR-sample, dimensions: (|113, 63)>
>>> info
{'time': 0.011758089065551758}

The reconstruction methods available

All you need to know for each method is:

  • what the method do (of course you need to know a little about it),

  • his nickname to give to restore(),

  • his parameters,

  • what informations are returned.

Restoration cheet sheet

Additional info in case PCA_transform is True is PCA_info which stores the following keys:

  • H: the truncated PCA basis,

  • PCA_th: the PCA threshold,

  • Ym: the data mean.

Interpolation

The interpolation method calls linear, cubic or nearest neighbor interpolation.

The method to give to the restore() method is interpolation. The associated function is resp. interpolate().

The input parameters are:

  • method: (optional, str) The interpolation method (among nearest, linear and cubic). Default is nearest neighbor.

The output dictionary stores the following informations:

  • time: the execution time (in sec.),

  • PCA_info: in case of 3D data with PCA pre-processing, it stores info about PCA.

L1

This regularized least-square method solves the following optimization problem:

\[\gdef \x {\mathbf{x}} \gdef \y {\mathbf{y}} \hat{\x} = \mathrm{arg}\min_{ \x\in\mathbb{R}^{m \times n} } \frac{1}{2} ||(\x-\y)\cdot \Phi||_F^2 + \lambda ||\x\Psi||_1\]

where \(\mathbf{y}\) are the corrupted data, \(\Phi\) is a subsampling operator and \(\Psi\) is a 2D DCT operator.

The method to give to the restore() method is L1. The associated function is resp. L1_LS().

The input parameters are:

  • Lambda: (float) The regularization parameter,

  • init: (optional, numpy array) An initial point for the gradient descent algorithm which should have the same shape as the input data,

  • Nit: (optional, int) The number of iterations.

The output dictionary stores the following informations:

  • E: The evolution of the functional value,

  • Gamma: The set of all pixel positions which coefficient in the DCT basis is non-zero,

  • nnz-ratio: The ratio of non-zero coefficients over the number of DCT coefficients,

  • time: the execution time (in sec.).

Smoothed SubSpace

The 3S algorithm denoise or reconstructs a multi-band image possibly spatially sub-sampled in the case of spatially smooth images. It is well adapted to intermediate scale images.

This algorithm performs a PCA pre-processing operation to estimate:

  • the data subspace basis \(\mathbf{H}\),

  • the subspace dimension \(R\),

  • the associated eigenvalues in decreasing order \(\mathbf{d}\),

  • the noise level \(\hat{\sigma}\).

After this estimation step, the algorithm solves the folowing regularization problem in the PCA space:

\[\gdef \S {\mathbf{S}} \gdef \Y {\mathbf{Y}} \gdef \H {\mathbf{H}} \gdef \I {\mathcal{I}} \begin{aligned} \hat{\S} &= \underset{\S\in\mathbb{R}^{m \times n \times R}}{\arg\min} \frac{1}{2R}\left\|\S \mathbf{D}\right\|_\mathrm{F}^2 + \frac{\lambda}{2}\sum_{m=1}^{R} w_{m} |\S_{m,:}|_2^2\\ &\textrm{s.t.}\quad \frac{1}{R}|\H_{1:R}^T\Y_{\I(n)}-\S_{\mathcal{I}(n)}|^2_2 \leq\alpha\hat{\sigma}^2,\ \forall n \in \{1, \dots,\ m*n\} \end{aligned}\]

where \(\mathbf{Y}\) are the corrupted data, \(\mathbf{D}\) is a spatial finite difference operator and \(\mathcal{I}\) is the set of all sampled pixels. The coefficient \(\alpha\) is a coefficient which scales the power of the data fidelity term.

For more details, see [BMonierOberlinBrun+18].

The method to give to the restore() method is 3S. The associated function is resp. SSS().

The input parameters are:

  • Lambda: (float) The regularization parameter,

  • scale: (optional, float) The spectr

  • init: (optional, numpy array) An initial point for the gradient descent algorithm which should have the same shape as the input data,

  • Nit: (optional, int) The number of iterations.

The output dictionary stores the following informations:

  • E: The evolution of the functional value,

  • time: the execution time (in sec.),

  • PCA_info: in case of 3D data with PCA pre-processing, it stores info about PCA.

Smoothed Nuclear Norm

The SNN algorithm denoise or reconstructs a multi-band image possibly spatially sub-sampled in the case of spatially smooth images. It is well adapted to intermediate scale images.

This algorithm solves the folowing optimization problem:

\[\gdef \X {\mathbf{X}} \gdef \Y {\mathbf{Y}} \gdef \H {\mathbf{H}} \gdef \I {\mathcal{I}} \hat{\X} = \underset{\X\in\mathbb{R}^{m \times n \times B}}{\arg\min} \frac{1}{2}||\Y_\I - \X_\I||_\mathrm{F}^2 + \frac{\lambda}{2}\left\|\X \mathbf{D}\right\|_\mathrm{F}^2 + \mu ||\X||_*\]

where \(\mathbf{Y}\) are the corrupted data, \(\mathbf{D}\) is a spatial finite difference operator and \(\mathcal{I}\) is the set of all sampled pixels.

For more details, see [BMonierOberlinBrun+18].

The method to give to the restore() method is SNN. The associated function is resp. SNN().

The input parameters are:

  • Lambda: (float) The \(\lambda\) regularization parameter,

  • Mu: (float) The \(\mu\) regularization parameter,

  • init: (optional, numpy array) An initial point for the gradient descent algorithm which should have the same shape as the input data,

  • Nit: (optional, int) The number of iterations.

The output dictionary stores the following informations:

  • E: The evolution of the functional value,

  • time: the execution time (in sec.),

  • PCA_info: in case of 3D data with PCA pre-processing, it stores info about PCA.

Cosine Least Square

The CLS algorithm denoises or reconstructs a multi-band image possibly spatially sub-sampled in the case of spatially sparse content in the DCT basis. It is well adapted to periodic data.

This algorithm solves the folowing optimization problem:

\[\gdef \X {\mathbf{X}} \gdef \Y {\mathbf{Y}} \gdef \H {\mathbf{H}} \gdef \I {\mathcal{I}} \hat{\X} = \underset{\X\in\mathbb{R}^{m \times n \times B}}{\arg\min} \frac{1}{2}||\Y_\I - \X_\I||_\mathrm{F}^2 + \lambda ||\X \Psi||_{2, 1}\]

where \(\mathbf{Y}\) are the corrupted data, \(\mathbf{D}\) is a spatial finite difference operator and \(\mathcal{I}\) is the set of all sampled pixels.

The method to give to the restore() method is CLS. The associated function is resp. CLS().

The input parameters are:

  • Lambda: (float) The \(\lambda\) regularization parameter,

  • init: (optional, numpy array) An initial point for the gradient descent algorithm which should have the same shape as the input data,

  • Nit: (optional, int) The number of iterations.

The output dictionary stores the following informations:

  • E: The evolution of the functional value,

  • Gamma: The set of all pixel positions which coefficient in the DCT basis is non-zero,

  • nnz-ratio: The ratio of non-zero coefficients over the number of DCT coefficients,

  • time: the execution time (in sec.),

  • PCA_info: in case of 3D data with PCA pre-processing, it stores info about PCA.

Post-Lasso CLS algorithm

This algorithms consists in applying CLS to restore the data and determine the data support in DCT basis. A post-least square optimization is performed to reduce the coefficients bias.

The method to give to the restore() method is Post_LS_CLS. The associated function is resp. Post_LS_CLS().

The input parameters are:

  • Lambda: (float) The \(\lambda\) regularization parameter,

  • init: (optional, numpy array) An initial point for the gradient descent algorithm which should have the same shape as the input data,

  • Nit: (optional, int) The number of iterations.

The output dictionary stores the following informations:

  • E_CLS: The evolution of the functional value for the CLS optimization step,

  • E_post_ls: The evolution of the functional value for the post-LS optimization step,

  • Gamma: The set of all pixel positions which coefficient in the DCT basis is non-zero,

  • nnz-ratio: The ratio of non-zero coefficients over the number of DCT coefficients,

  • time: the execution time (in sec.),

  • PCA_info: in case of 3D data with PCA pre-processing, it stores info about PCA.

ITKrMM and wKSVD

Weighted K-SVD (see [BMairalEladSapiro08]) and Iterative Thresholding and K residual Means for Masked data (see [BNS18]) methods.

The wKSVD and ITKrMM algorithms share a lots of their code so that their input and output are the same. Though, two implementations exist to run these algorithms: one with python (ITKrMM and wKSVD methods) and one with maltab (ITKrMM_matlab and wKSVD_matlab methods). The original Matlab codes are broadcasted by Karin Schnass. They were translated afterwards into python. Nothing distinguish them but for wKSVD where matlab is faster. The only problem is that you should have the matlab command in your system path.

The methods to give to the restore() method are ITKrMM, wKSVD, ITKrMM_matlab or wKSVD_matlab. The associated functions are resp. ITKrMM(), wKSVD(), ITKrMM_matlab() and wKSVD_matlab().

The input parameters are:

  • Patchsize: (optional, int) The patch width,

  • K: (optional, int) The dictionary size (incl. low-rank component),

  • L: (optional, int) The number of low-rank components to estimate,

  • S: (optional, int) The sparsity level,

  • Nit: (optional, int) The number of iterations for the dictionary estimation.

  • Nit_lr: (optional, int) The number of iterations for the low-rank estimation.

The output dictionary stores the following informations:

  • dico: The dictionary,

  • E: The evolution of the error,

  • time: the execution time (in sec.),

  • PCA_info: in case of 3D data with PCA pre-processing, it stores info about PCA.

BPFA

Beta Process Factor Analysis algorithm (see [BXZC+12]).

As for wKSVD and ITKrMM, BPFA is based on a Matlab code from Zhengming Xing (these codes were broadcasted without any license). The python code just calls it, so matlab should be in the path system so that the matlab command could be called from the command line.

The method to give to the restore() method is BPFA_matlab. The associated function is resp. BPFA_matlab().

The input parameters are:

  • Patchsize: (optional, int) The patch width,

  • K: (optional, int) The dictionary size,

  • step: (optional, int) That’s the pixel space between two consecutive patches (if 1, full overlap),

  • Nit: (optional, int) The number of iterations for the dictionary estimation.

The output dictionary stores the following informations:

  • dico: The dictionary,

  • time: the execution time (in sec.),

  • PCA_info: in case of 3D data with PCA pre-processing, it stores info about PCA.

That’s all folks !

This was the main content of the documentation. Congrats, you understood 90% of this library :)

References

BNS18

Valeriya Naumova and Karin Schnass. Fast dictionary learning from incomplete data. EURASIP journal on advances in signal processing, 2018(1):12, 2018.

BXZC+12

Z. Xing, M. Zhou, A. Castrodad, G. Sapiro, and L. Carin. Dictionary learning for noisy and incomplete hyperspectral images. SIAM J. Imag. Sci., 5(1):33–56, 2012.

BMairalEladSapiro08

J. Mairal, M. Elad, and G. Sapiro. Sparse representation for color image restoration. IEEE Trans. Image Process., 17(1):53–69, Jan 2008.

BMonierOberlinBrun+18(1,2)

É. Monier, T. Oberlin, N. Brun, M. Tencé, M. de Frutos, and N. Dobigeon. Reconstruction of partially sampled multiband images—application to stem-eels imaging. IEEE Trans. Comput. Imag., 4(4):585–598, dec. 2018.