Add method to calibrate dataset/ refine spectral line poisition

A common issue when running espm is that the peak poisitions of the G matrix do not perfectly align with where the peaks are in the dataset, as well as their shape.

A good solution would be to perform a "precalibration" step, where the average spectrum is fit by the G matrix but alowing shifting the positions of the gaussians.

The idea would be to fit the average data with a set of free gaussians with initial positions and intensity ratios set by the kV_xrays.json tables and then fit with least squares, allowing small energy shifts and small changes in width. The result of this fit should then be used to set the energy scale, energy offset, gaussian width and intercept (as is called in espm) and small shifts in the line positions).

Ideally, this can set a proper G matrix that will improve the results of the subsequent NMF decompositon.

The addition of this functionality will be part of Perry's summer intership with us.

I think this is a crucial feature that needs to be implemented as soon as possible. I was planning to do this but I couldn't find the time.

I have a few suggestions :

I don't recommend setting the ratios of the gaussian but to use fully free gaussians (except that they should be in a certain energy range). The current energy calibration method is not good but is not too far away either so it should be good starting point for the fit.
Originally I was thinking of implementing the fit such that it produces a table similar (same format) to the kV_xrays.json table but with modified energies. It would then be stored internally in the .hspy file. The espm modelling would then use the generated tables for performing decomposition, etc ...
It would then be easier for the user to come back to the result of the fit and even tune the table if needed. It would also prepare a next step : we could provide the user with functions to visualize or browse espm tables and even apply modifications (such as applying experimentally determined k factors).

By the way :

It could be interesting not to use gaussian to fit but 2nd order polynomials. It only works well over a limited number of points (less than 10) around the maximum of the peak. It is basically similar to performing a taylor development, that's why it is only valid very locally. The advantage are :

It doesn't care about the actual correct theoretical model
It is very fast to calculate
It is very stable, if you have peaks around, a background, etc ... It is not an issue.

It is a well known tricks in Xray diffraction.

in fact i have tried fitting a degree N-1 polynomial using N points, where N is the number after filtering out bad lines. the problem is that low energy lines can act quite funny.

And did you try 2nd order on N points ?

yes actually i plotted it out vs different degrees -- all degrees of polynomial agree that the energy misalignment is linear (they start diverging after 15 keV since there are no sample lines after 15)

fitting peaks individually independent of each other (dashed lines are theoretical lines, and dotted lines are calibrated lines)

the problem with this approach is that, for example grey lines are Re, which has a very low concentration, thus its lines basically disappeared into other bell curves or the background, making its naive fitted lines unrealistically far away from where it should be.

calibrated peaks using polynomial of degree 3

using polynomial, the problem for low concentration or low cs lines is stablised, but the polynomial calibration is actually even worse then the theoretical table for large peaks.

i will use some other strategies. more coming next monday.

Could you, if possible, detail a bit more the last two graphs ?

What model did you choose to model the peak position (with equation) ? On how many points ? Did you put boundaries to the fit ? Could you show a list of results (peak name, peak energy (tabulated), peak energy (fitted), peak intensity, etc ...) ?

Shall we organise a remote meeting ? (if necessary)

I am not sure what you mean by model, but i read the 200keV table and the dashed (- - - -) lines are the lines present in the table. The dotted (.....) lines are calibrated lines.

For the first plot, it was done by the following algorithm:

def gauss(A, x, x0, sigma, C, m):
    return A * np.exp(-((x - x0) ** 2) / (2 * sigma**2)) + C + m * (x - x0)

$A\cdot\mathcal N(x_0, \sigma) + (m (x-x_0)+C)$ where $m(x-x_0)+C$ is the background noise.

def fit_single_peak(E_theoretical, sigma_expected, window):
    mask = (energy_axis >= E_theoretical - window) & (
        energy_axis <= E_theoretical + window
    )
    xdata = energy_axis[mask]
    ydata = average_spectrum[mask]

    idx_max = np.argmax(ydata)
    x0_guess = xdata[idx_max]
    A_guess = np.max(ydata) - np.min(ydata)
    sigma_guess = sigma_expected
    C_guess = np.min(ydata)

    p0 = [A_guess, x0_guess, sigma_guess, C_guess, 0.0]
    bounds = (
        [0.0, E_theoretical - window, 0.9 * sigma_expected, 0.0, -1.0],
        [np.inf, E_theoretical + window, 1.1 * sigma_expected, np.inf, 1.0],
    )
    popt, _ = curve_fit(
        gauss,
        xdata,
        ydata,
        p0,
        bounds=bounds,
    )
    return popt

The above function is ran either for all lines in the table, or a selected subset of lines, with sigma_expected = model.width_slope * E_theoretical + model.width_intercept / 2.3548 and window = 3 * sigma_expected. It basically tries to fit gauss in the ball of radius window around E_theoretical, where E_theoretical is the line energy given by the table.

For the second plot, we collect the E_theoretical and the returned calibrated energies, poly fit a degree $d$ polynomial $p$, and set calibrated energy = p(E_theoretical).

Let's look at the polynomial $p$ fitted using different degrees:

I imagine we want the resulting polynomial to look linear, so any $p$ of degree $>8$ is not good. Each dot in this plot is the theoretical energy vs calibrated energy. Notice that there are a lot of pink and grey dots below the black line:

They are Pt and Re, which are of very low concentration in the sample, meaning they get fitted into another peak on the left, which in turns drag the fitted $p$ away from what it should be.

To combat this, I have added a weight for each line in the poly fit by the peak magnitude, the result is way better compared to plot 2 from the previous comment:

I cannot get over the fact that it can happen that the lines of an element with low concentration aligns closely to the lines of another abundant element, making it "very important" during the fitting. Also, I just don't see how this is better than the existing affine calibration using 2 lines.

In fact, if we know which element has the highest concentration, this problem would have been much easier.

Ok, I see one main issue here. I think that instead of fitting the peaks with gaussians we should fit each peak with a 2nd order polynomial on a severly reduced amount of points for each peak. I believe, that second order polynomial results can be sorted into correct ones :

the maximum of the polynomial sits in the selected range
and discard incorrect ones :
The maximum of the polynomial sits outside (or on the edge) of the selected range.
You don't need all peaks to participate to the calibration.
You can also constraint the polynomial to negative 2nd degree coeff, so that it has a bell shape.

For the energy scale, so far we should keep the 1st order polynomial. For quite a few datasets there are peaks beyond 15keV.

You can also constraint the polynomial to negative 2nd degree coeff, so that it has a bell shape.

Is there a physical reason for this? I honestly also expect a 2nd degree polynomical to be largely enough, but why just negative coefficient?

@perrymok

In fact, if we know which element has the highest concentration, this problem would have been much easier.

Why can know by runing the quan spectrum function on the average spectrum.
Aside from that I'd like to know how (est.G@[email protected]).mean(1) compares to s.data.mean((0,1)) once you run decompositions on calibrated/uncalibrated data to see how well the calibration works.
@perrymok do you think you can produce this plots in on a couple of the datasets you have?

refer to the calibration branch for the polyfit calibration and an example notebook in the notebook folder

as a side note, i think we can include the plot_table function in the new calibration.ipynb in the library as a utility function.

In fact, if we know which element has the highest concentration, this problem would have been much easier.

Why can know by runing the quan spectrum function on the average spectrum. Aside from that I'd like to know how (est.G@[email protected]).mean(1) compares to s.data.mean((0,1)) once you run decompositions on calibrated/uncalibrated data to see how well the calibration works. @perrymok do you think you can produce this plots in on a couple of the datasets you have?

pushing soon

@PauTorru

I am sorry this comment is coming a bit late. I believe we side-tracked a little bit.

I searched a bit the litterature on the topic and technically the detectors are made so that there is a linear relationship between the bin number (1) and the energy. I couldn't find a more refined analytical formula (I hope I didn't miss an obvious paper). The calibration can already be done using hyperspy's calibration method (and that's why eds_spim inherits from EDSTEMSpectrum object). I don't think we need more than that in terms of calibration.

The issue with espm NMF is that we need the peak position very precisely for each peak (if possible). The theoretical tables are not accurate enough and their errors seem random. That's why I think we should design a method that finds the experimental position of peaks in the spectrum (based on the already avaible detector calibration). Then we can produce a corrected table of x ray peaks. In your graphs, the idea would be to project the peak energies on the straight line.

I think the peak position correction can be done in two steps

Sort out the peak that are easy to find from the rest. Often, some peaks overlaps and thus it is probably hard to find their position accurately.
1. @PauTorru If we use the 2nd order polynomial to find peaks, we can levrage the fact that the 2nd order coefficient has to be negative for the polynomial to have a bell shap. If not, it means that the peak is ill-defined.
2. From the theoretical table you can probably already determine troublesome overlaps (Ok) and low intensity peaks (Trickier, it still depends on concentration).
From a good enough quantity of calibrated peaks we can probably interpolate the positions of the other peaks. It needs to be tested though.
We may want to raise a warning if the method didn't find enough peaks.

@perrymok I hope I did not misunderstood your previous efforts. But it seemed to me we were heading in the wrong direction.

(1) Just to clarify : Each X ray detected is assigned to a bin from 0 up to 2048 (or 4096) bins. Then using the usually already known calibration given by the manufacturer we can assign each bin to an energy in keV.

So the manufacturers provide an affine map from the bin number to the energy that it thinks it is observing (call this empirical energy), and we need a calibration map from the empirical energy to the actual energy. i was working on the latter (specifically we focus on mappings that are polynomial of degree 2 or 3), so I don't think we were heading in the wrong direction.

ii. From the theoretical table you can probably already determine troublesome overlaps (Ok) and low intensity peaks (Trickier, it still depends on concentration).

I agree, I have done some of it.

I don't entirely understand why we need to fit the peaks using deg 2 polynomials (with negative 2nd coefficient) instead of bell curves, though.

@adriente

I believe Perry already has what you expect Adrien, using the flag use_poly=False. That is, just redifining the peak position to be the one that has been fitted.

@adriente i think the second order polynomial to find peak positios is not that important, as the gaussian works well to find the position of strong lines. It runs sufficently fast and accurately.

@perrymok I do agree with adrien and believe there is some misconception and we are straying to far away.

Let's say we identify a series of peaks P and find through fitting that their center is on the energy channel x_P. There is also rough energy calibration E_empirical(x) = offset+scale*x, that gives a good rough estimate of the energy a which the peak lies. In an ideal world ther may exist a function E(x_P)=E_P_theoretical, that puts all peaks in their theoretical energy (and all energy channels in between), and ideally this could be a low order polynomial. However, we have discussed a few times that it is enough to just say tell the library that E_P_theortical = E_empirical(x_P) and leave all points in between with their rough estimate. We can also improve E_empirical(x) with the fit of the main peak positions, and we can even say that E_empirical is an order two polynomial to better adress energy channels in between peaks. However, we still need to to correct the table that generates the G matrix by setting E_P_theortical = E_empirical(x_P), because that is the most important part for espm to work properly. I hope the nomenclature I introduced clarifies things a little.

There is extensive practical and theoretical justification on why this needs to be done (for example the same line for the same element can be at slightly different positions for differnt samples depending on the chemistry of the compound). So the "hacky" step of setting the peak position always needs to be done.

We have discussed this but I think leaving it written here will make things more clear.

In parallel, I have starte testing the code and I am finding some bugs.
fit_single_peak gives sometimes gives an error. I believe it is when the maks for the peak falls out of bounds. @perrymok pleas address.

Navigation Menu

Repository navigation

Add method to calibrate dataset/ refine spectral line poisition #6

Metadata

Assignees

Labels

Projects

Milestone

Relationships

Development

Add method to calibrate dataset/ refine spectral line poisition #6

Description

Activity

adriente commented on Jul 2, 2026

adriente commented on Jul 9, 2026

perrymok commented on Jul 9, 2026

adriente commented on Jul 9, 2026

perrymok commented on Jul 10, 2026

adriente commented on Jul 13, 2026

perrymok commented on Jul 13, 2026

adriente commented on Jul 13, 2026

PauTorru commented on Jul 13, 2026

PauTorru commented on Jul 13, 2026

perrymok commented on Jul 13, 2026

perrymok commented on Jul 13, 2026

adriente commented on Jul 16, 2026

perrymok commented on Jul 16, 2026

PauTorru commented on Jul 16, 2026

Metadata

Metadata

Assignees

Labels

Projects

Milestone

Relationships

Development

Issue actions