NeXL

[Ben Buse]

Ok I'm starting to get a hang of Julia and NeXL, but can't manage to get the quant results, it says no take of angle, despite adding it
I'm working through https://pages.nist.gov/NeXLSpectrum.jl/k412refs/

So far this is what I've done

Code Select Expand

using NeXLSpectrum
using NeXLMatrixCorrection
refs = references([reference(n"Mn",joinpath(path12,"2406SPC.emsa"),mat"Mn"),reference(n"Al",joinpath(path13,"2406SPC.emsa"),mat"Al2O3"),reference(n"Mg",joinpath(path14,"2406SPC.emsa"),mat"MgO"),reference(n"Fe",joinpath(path15,"2406SPC.emsa"),mat"Fe2O3"),reference(n"Si",joinpath(path16,"2406SPC.emsa"),mat"Si59Mg107Fe11O240"),reference(n"O",joinpath(path13,"2406SPC.emsa"),mat"Al2O3"),],130)
res=fit_spectrum(s1,refs)

but couldn't figure out how to quantify single result - alternative command to

Code Select Expand

quant = quantify.(res)

So did

Code Select Expand

unks = [s1,s1]
res2= [fit_spectrum(u,refs) for u in unks]
quant=quantify.(res2)

Then no take off angle so tried

Code Select Expand

unks[1][:TaKeOffAngle]=40
unks[2][:TaKeOffAngle]=40

Still no take off angle when try quant.

Should take off angle be in degrees or radians?

Maybe I need to add take of angle to standards? - How do I do this?

But got this nice fit. What do the numbers of the plot represent?


And how do I get kratio and net intensity?

Is there a nicer application to run Julia in then it's console?

[JonF]

Thanks for sharing, NeXL isn't something I was familiar with and from what you've already shown, it seems to have potential. It would be great if we could move away from OEM-based EDS interpretation and just take the raw spectrum (regardless of manufacturer) and quant it ourselves.

I'll admit to being a bit naive on this front, but if we're looking at collecting standards and unknowns on the same EDS system (e.g. the one bolted to our probes) - as we do with WDS - would all the detector-specific bits (window compositions, fano factors etc etc) get cancelled out as we create K ratios via LLSQ fitting? 

As both Julia and NeXL are open source (MIT and Unilicense, respectively), is this the kind of thing we could build in to PfE in the future?

Is there a nicer application to run Julia in then it's console?

There's a Julia plugin for VSCode: Julia VSCode

[Ben Buse]

Yes part of consistency between instruments is being able to apply same procedures to the raw data.

My current fudge until I can work out how to add take off angle to standards - or someone can tell me, is I can add the take off angle to the input files emsa/msa by adding

Code Select Expand

#ELEVANGLE  : 40.0The file keywords being https://microscopy.org/files/galleries/MSA-MAS_FileFormatSpecification_1991.pdf but unfortunately it seems take off angle isn't mandatory or hasn't been adopted.

Giving me a quant of

Code Select Expand

Converged to [O=0.3782,Mg=0.0000,Al=0.1036,Si=0.1705,Mn=0.3207,Fe=0.0037] in 6 steps.But still question how do I get the raw kratios and the net intensities?

[Nicholas Ritchie]

Ben,  You've made great progress.  I started to answer your first iteration of this question but I had guests... Julia and NeXL aren't the easiest to get going with.  Julia is "theoretically" a wonderful language (it does some of the most amazing things) but in practice it can be a bit of a pain.
```
s[:TakeOffAngle] = deg2rad(40.0)  # The take-off angle is in radians (why?  just because...)
quant = quantify.(specs) # Is a vector operation, the dot implies apply the function quantify to all elements of the container 's'.  The function of a single argument equivalent is:
quant = quantify(spec)
```
The numbers on the fit plot are the k-ratios for each fitted ROI.

I'll provide a couple more answers after I've had a little time to work up some examples.

Take care, Nicholas

[Nicholas Ritchie]

Let me help you out by going line-by-line through Fitting K412 (simple API) I'll not replicate the commentary on the page but add context.

First load the necessary `packages` as Julia calls libraries.  Then we are going to add references.  References provide the spectral data and meta-data necessary to fit the characteristic peaks in the measurand spectra.

The `references(...)` function takes a vector of FilteredReference instances created using the reference(...)
function like this:

Code Select Expand

reference(n"Fe", joinpath(path, "Fe std.msa"), mat"Fe")
Vectors are constructed with matching `[` and `]` brackets.

The first argument to

Code Select Expand

references(...) uses the

Code Select Expand

n"Xx" syntax which can be used to construct Element, AtomicSubShell and CharXRay instances depending upon the format of the string within the quotes. 

Code Select Expand

n"Fe" produces an Element representing Fe. 

Code Select Expand

n"Fe L3" produces a AtomicSubShell representing the iron L3 sub-shell. 

Code Select Expand

Fe K-L3 produces a CharXRay. 

Code Select Expand

reference(...) requires an Element so...

Code Select Expand

refs = references(
  [
    # Specify a reference for iron (arg 1), in "Fe std.msa" (arg2) that is pure iron (arg3).
    reference(n"Fe", joinpath(path, "Fe std.msa"), mat"Fe"),
    # specify a conductive surface coating layer using the `coating` named argument
    reference(n"Si", joinpath(path, "SiO2 std.msa"), mat"SiO2", coating = Film(pure(n"C"), 10.0e-7)),
    reference(n"O", joinpath(path, "SiO2 std.msa"), mat"SiO2", coating = Film(pure(n"C"), 10.0e-7)),
    reference(n"Ca", joinpath(path, "CaF2 std.msa"), mat"CaF2", coating = Film(pure(n"C"), 10.0e-7)),
    # Read the composition from the spectrum file's ##D2STDCMP tag
    reference(n"Mg", joinpath(path, "MgO std.msa"), coating = Film(pure(n"C"), 10.0e-7)),
    # Read the conductive coating from the spectrum file's ##CONDCOATING tag
    reference(n"Al", joinpath(path, "Al2O3 std.msa"), mat"Al2O3"),
  ],
  132.0  # Detector resolution at Mn Kα (eV)
)produces references for the elements O, Mg, Al, Si, Ca and Fe.

The

Code Select Expand

mat"..." syntax produces Composition instances.  It too is flexible.

Code Select Expand

mat"0.3*Fe+0.6*Cr" produces 30% Fe by mass, 60% Cr by mass

Code Select Expand

mat"Fe2O3" produces 2 atoms of Fe and 3 atoms of O.

Code Select Expand

mat"Ca5(PO4)3OH" produces 5 atoms Ca, 3 atom of P, 13 atoms of O and 1 atom of H.

The

Code Select Expand

132.0 specifies the FWHM which is the only extra piece of information needed that can't be pulled from the spectra. (See

Code Select Expand

BasicEDS[..] below which shows the other contextual information used to describe the EDS detector.

Code Select Expand

# Note: This is output - not code
References[
    BasicEDS[4096 chs, 1.63032 + 9.99856⋅ch eV, 132.0 eV @ Mn K-L3, 1 ch LLD,
[Be,Sc,Ba,Pu]],
    k[Fe L3-M5 + 13 others, Fe],
    k[Fe K-L3 + 1 other, Fe],
    k[Fe K-M3 + 3 others, Fe],
    k[Si K-L3 + 3 others, SiO2],
    k[O K-L3 + 1 other, SiO2],
    k[Ca K-L3 + 3 others, CaF2],
    k[Mg K-L3 + 1 other, MgO],
    k[Al K-L3 + 3 others, Al2O3],
]

The we load the measurand

Code Select Expand

Spectrum object instances using

Code Select Expand

loadspectra(...).  We assign an additional piece of contextual information about the sample coating using the

Code Select Expand

:Coating Symbol.  (Note: If you provide coating meta-data for the reference also provide it for the unknown.  Otherwise, your results may vary... If the coating is the same for both, you are probably better off not defining :Coating for either.) 

Code Select Expand

Spectrum objects can be indexed using integers to access the channel data or by "Symbols" to access property data.  If `s` is a Spectrum then:

Code Select Expand

s[:ProbeCurrent] # returns the probe current in nA
s[120] # returns the 120th channel in the spectrum (remember [1] is the first channel)
s[120:130] # returns channels 120 to 130 inclusive.
sum(s[120:130]) # returns the sum over channels 120 to 130.
s[1234.0] # Returns counts in the channel at 1234.0 eV
s[n"Fe K-L3"] # returns the channel at the Fe Kalpha energy


Code Select Expand

# Now load all the unknown spectra and assign a carbon coating`
unks = map(0:4) do i
    s = loadspectrum(joinpath(path, "III-E K412[$i][4].msa"))
    # assign a carbon coating
    s[:Coating] = Film(pure(n"C"), 30.0e-7)
    s
end

Other important properties are

Code Select Expand

:LiveTime, :ProbeCurrent, :TakeOffAngle, :BeamEnergy which are in seconds, nano-amps, radians (not degrees) and electron-volts (not keV).  Note the `:` character means that these labels called Symbols not variable names.

The

Code Select Expand

plot(...) function is defined in Gadfly but extended in NeXL to understand spectra and other NeXL specific data types. (In Julia, functions can have many different implementation depending upon the argument types).  So my implementation of plot uses Gadfly as a basis but implements EPMA related items like "klms".

In NeXLSpectrum, the

Code Select Expand

fit_spectrum(...) function fits the references to the measurands and returns a FilterFitResult containing the resulting k-ratios. Thus

Code Select Expand

res= [ fit_spectrum(u,refs) for u in unks ]
produces a Vector of FilterFitResult instances one for each Spectrum.

The

Code Select Expand

quantify function takes

Code Select Expand

FilterFitResult results and performs matrix correction.

Code Select Expand

quant = quantify.(res) # The dot maps the function quantify across all FilterFitResult instances in res
q = quantify(res[1])   # Applies the function quantify to `res[1]` - the first FilterFitResult instance


The subsequent

Code Select Expand

plot(...) functions plot the results using specializations of Gadfly functions for Compositions.

[Nicholas Ritchie]

There is a lot more functionality to each of these objects that is detailed here.

A method you will often find useful is the

Code Select Expand

NeXLUncertainties.asa(...) function. To use it, load the DataFrames and NeXLUncertainties packages.

Code Select Expand

using DataFrames, NeXLUncertaintiesThen using the syntax:

Code Select Expand

asa(DataFrame, res)for various different kinds of `res` will produce a tabulation of the data items (k-ratios, compositions, etc.)  DataFrames is Julia's equivalent to Python's Pandas.  It can slice and dice data in wonderful ways.

Code Select Expand

using CSV will load the CSV package which can read and write DataFrames to CSV files.  There are alternatives for other file types like HDF5 or Excel or ...

Code Select Expand

CSV.write("data.csv", asa(DataFrame, cs))


[Nicholas Ritchie]

Code Select Expand

HyperSpectra work similar to Spectrum objects but with 2D indexes.

Code Select Expand

using Unitful  # For the units in the fov

lt = 0.72*4.0*18.0*3600.0/(1024*1024) # 18.0 hours on 4 detectors (1.0-0.72 dead time)

hs = NeXLSpectrum.compress(HyperSpectrum(
           LinearEnergyScale(0.0,10.0),
           Dict{Symbol,Any}(
             :TakeOffAngle => deg2rad(35.0),
             :ProbeCurrent => 1.0,
             :LiveTime => lt,
             :BeamEnergy => 10.0e3,
             :Name => "Mn Nodule"
           ),
           readrplraw(raw"/home/nwmr/Documents/Spectra/Fe3C_10keV_map/map[15]"),
           fov = [ 4.096u"mm", 4.096u"mm"], offset= [ 0.0u"mm", 0.0u"mm" ]
       ))

This code loads the data from a RPL/RAW file (`readrplraw(...)`) and creates a HyperSpectrum from the resulting data then compresses it into the smallest representation that won't lose any channel data (smallest data type and number of channels.)  It applies the specified :TakeOffAngle, :ProbeCurrent, :LiveTime, :BeamEnergy since this isn't typically in the RPL header.  It also specifies the dimensions of the images via the `fov` property.  The detector has a `LinearEnergyScale(...)` with 10.0 eV/ch and 0.0 eV offset.

Code Select Expand

hs[1,1] # returns the Spectrum at row 1, column 1
hs[2,3][120] # Returns the 120th channel in the Spectrum at row 2, column 3
hs[2,3][120:130] # Returns the 120 to 130 channel (inclusive) in ...
sum(float.(hs[2,3][120:130]))  # This converts channels 120 to 130 to a float-type and then sums the result.
# The conversion to float ensures that the sum doesn't overflow as frequently HyperSpectrum object's data
# is represented in short unsigned integer types u8 or u16.
sum_of_120_to_130 = [ sum(float.(hs[ci][120:130])) for ci in CartesianIndices(hs) ] # Creates a 2D matrix of sums of the channels 120 to 130


[Nicholas Ritchie]

And there is more...

Code Select Expand

plane(hs, n"Fe K-L3")
will integrate one FWHM channel range centered on Fe Kalpha producing a matrix of values.

Code Select Expand

colorize(hs, [ n"Fe K-L3", n"Fe L3-M5", n"C K-L2"])
produces a RGB image with intensities around these characteristic X-ray lines. See Images.jl for details on working with images.

To quantify HyperSpectrum objects, use the same functions to create a references (as above) which you apply in the same way:

Code Select Expand

res=quantify(hs, refs)
except it produces a Materials instance that is a more efficient way of expressing a large matrix of compositions.  Index it like:

Code Select Expand

res[2,3]
to extract Material instances or

Code Select Expand

res[2,3][n"Fe"]
for individual mass fractions.  (See also NeXLCore documentation) and NeXLMatrixCorrection documentation.

[Ben Buse]

Thanks Nicholas for all this detailed information, I'm just working through it. I might have missed it but how do I add take off angle to references if it's not contained in the msa, can you address references individually after they are created. Dataframes works a treat. I'm now going to start to look at quantifying the JEOL spectrum map as extracted via Yuji code into hyperspy as discussed here https://smf.probesoftware.com/index.php?topic=1836.0.

Looking at this

Code Select Expand

LiveTime => lt,Seems live time is fixed for every pixel in the map is it possible to use an array of live time

Using Nicholas example above, I loaded the map save from hyperspy as rpl,

Code Select Expand

hs = NeXLSpectrum.compress(HyperSpectrum(
                  LinearEnergyScale(0.0,40.0),
                  Dict{Symbol,Any}(
                    :TakeOffAngle => deg2rad(35.0),
                                :ProbeCurrent => 1.0,
                    :LiveTime => lt,
                    :BeamEnergy => 30.0e3,
                    :Name => "8 sp1 spess"
                  ),
                  readrplraw(raw"C:\\Users\\glxbb\\julia\\edstry2"),
                            fov = [ 60u"mm", 60u"mm"], offset= [ 0.0u"mm", 0.0u"mm" ]
              ))
60 × 60 HyperSpectrum{UInt16,(:Y, :X)}[8 sp1 spess, 0.0 + 40.0⋅ch eV, 512 ch]
And plotted the spectrum for a pixel

Code Select Expand

plot(hs[1,1])
plot(hs[60,60])

Ok by plotting element labels, I realise the hyperspy energy offset wasn't carried through
So here's with energy offset

Code Select Expand

hs2 = NeXLSpectrum.compress(HyperSpectrum(
                  LinearEnergyScale(-1538.3,40.0),
                  Dict{Symbol,Any}(
                    :TakeOffAngle => deg2rad(35.0),
                                :ProbeCurrent => 1.0,
                    :LiveTime => lt,
                    :BeamEnergy => 30.0e3,
                    :Name => "8 sp1 spess"
                  ),
                  readrplraw(raw"C:\\Users\\glxbb\\julia\\edstry2"),
                            fov = [ 60u"mm", 60u"mm"], offset= [ 0.0u"mm", 0.0u"mm" ]
              ))
60 × 60 HyperSpectrum{UInt16,(:Y, :X)}[8 sp1 spess, -1538.3 + 40.0⋅ch eV, 512 ch]And then max pixel plot

Code Select Expand

mp2 = maxpixel(hs2)
plot(mp2, klms=[n"C",n"O",n"Al",n"Si",n"Ca",n"Mn",n"Fe"])
To give

And looking at that using the plus key to zoom in on gadfly plot realised I needed to shift it a further 100eV
And there's also plot sum spectra

Code Select Expand

plot(sum(hs3), klms=[n"C",n"O",n"Al",n"Si",n"Ca",n"Mn",n"Fe"])I reading Nicholas paper https://pmc.ncbi.nlm.nih.gov/articles/PMC9437143/ which also explains creation of ROI maps


Now need to figure out how to quantify the map

[Ben Buse]

Can I make a map sum a standard?

Code Select Expand

MapSum = sum(hs3)
MapSum[:LiveTime]=hs3[1,1][:LiveTime]*60*60
maprefs = references([reference(n"Al",MapSum,mat"Ca3Mn57Al40Si60O240")],128)
ERROR: DomainError with -2082.385000000002:
sqrt was called with a negative real argument but will only return a complex result if called with a complex argument. Try sqrt(Complex(x)).


[Nicholas Ritchie]

Looking at this

Code Select Expand

LiveTime => lt,Seems live time is fixed for every pixel in the map is it possible to use an array of live time

Using Nicholas example above, I loaded the map save from hyperspy as rpl,

Code Select Expand

hs = NeXLSpectrum.compress(HyperSpectrum(
                  LinearEnergyScale(0.0,40.0),
                  Dict{Symbol,Any}(
                    :TakeOffAngle => deg2rad(35.0),
                                :ProbeCurrent => 1.0,
                    :LiveTime => lt,
                    :BeamEnergy => 30.0e3,
                    :Name => "8 sp1 spess"
                  ),
                  readrplraw(raw"C:\\Users\\glxbb\\julia\\edstry2"),
                            fov = [ 60u"mm", 60u"mm"], offset= [ 0.0u"mm", 0.0u"mm" ]
              ))
60 × 60 HyperSpectrum{UInt16,(:Y, :X)}[8 sp1 spess, 0.0 + 40.0⋅ch eV, 512 ch]


You can specify per pixel live-times by adding the `livetime` argument to the HyperSpectrum constructor assigning it a matrix the same dimensions as the RPL/RAW data. (Remove the `:LiveTime => lt,`)

Code Select Expand

    readrplraw(....),
    livetime=fill(get(props, :LiveTime, 1.0), (512, 512))
)


[Nicholas Ritchie]

You can also construct references in two steps.  1) Load the spectrum and set properties. 2) Replace the reference spectrum path with the loaded spectrum.

Code Select Expand

spec=pathspectrum("path\to\references\SiO2")
spec[:LiveTime] = 60.0
reference(n"Si", spec, mat"SiO2")


[Nicholas Ritchie]

Can I make a map sum a standard?

Code Select Expand

MapSum = sum(hs3)
MapSum[:LiveTime]=hs3[1,1][:LiveTime]*60*60
maprefs = references([reference(n"Al",MapSum,mat"Ca3Mn57Al40Si60O240")],128)
ERROR: DomainError with -2082.385000000002:
sqrt was called with a negative real argument but will only return a complex result if called with a complex argument. Try sqrt(Complex(x)).


I'm not exactly sure where the problem lies.  It may be a overflow in

Code Select Expand

hs3[1,1][:LiveTime]*60*60
Try:

Code Select Expand

float(hs3[1,1][:LiveTime])*60*60  # or
hs3[1,1][:LiveTime]*60.0*60.0


[Ben Buse]

Thanks Nicholas, I will give it a go.

Now I'm trying to figure out how to make my map and standards have the same detector, when the map has a large energy offset, compared to the emsa standards which have no energy offset. I've reloaded map with 2048 channels

Code Select Expand

60 × 60 HyperSpectrum{UInt16,(:Y, :X)}[8 sp1 spess, -1638.3 + 10.0⋅ch eV, 2048 ch]
refs4
References[
        BasicEDS[2048 chs, 0.0 + 10.0⋅ch eV, 130.0 eV @ Mn K-L3, 2 ch LLD, [Be,Sc,Ba,Pu]]

so

Code Select Expand

resmap4 = fit_spectrum(hs4,refs4)
ERROR: AssertionError: The detector for the hyper-spectrum must match the detector for the filtered references.
Stacktrace:

Any suggestions

[Ben Buse]

Can I make a map sum a standard?

Code Select Expand

MapSum = sum(hs3)
MapSum[:LiveTime]=hs3[1,1][:LiveTime]*60*60
maprefs = references([reference(n"Al",MapSum,mat"Ca3Mn57Al40Si60O240")],128)
ERROR: DomainError with -2082.385000000002:
sqrt was called with a negative real argument but will only return a complex result if called with a complex argument. Try sqrt(Complex(x)).


I'm not exactly sure where the problem lies.  It may be a overflow in

Code Select Expand

hs3[1,1][:LiveTime]*60*60
Try:

Code Select Expand

float(hs3[1,1][:LiveTime])*60*60  # or
hs3[1,1][:LiveTime]*60.0*60.0


Here the problem seems to be using the map spectra as a reference

Code Select Expand

maprefs = references([reference(n"Al",hs3[1,1],mat"Al")],128)
ERROR: DomainError with -2082.385000000002:
sqrt was called with a negative real argument but will only return a complex result if called with a complex argument. Try sqrt(Complex(x)).
Stacktrace:
  [1] throw_complex_domainerror(f::Symbol, x::Float64)
    @ Base.Math .\math.jl:33
  [2] sqrt
    @ .\math.jl:627 [inlined]
  [3] resolution_to_fwhm
    @ C:\Users\glxbb\.julia\packages\NeXLSpectrum\eVZfj\src\detector.jl:155 [inlined]
  [4] resolution
    @ C:\Users\glxbb\.julia\packages\NeXLSpectrum\eVZfj\src\detector.jl:165 [inlined]
  [5] resolution(eV::Float64, det::BasicEDS)
    @ NeXLSpectrum C:\Users\glxbb\.julia\packages\NeXLSpectrum\eVZfj\src\detector.jl:309
  [6] buildfilter(::Type{Float64}, ty::Type{G2Filter}, det::BasicEDS, a::Float64, b::Float64)
    @ NeXLSpectrum C:\Users\glxbb\.julia\packages\NeXLSpectrum\eVZfj\src\filter.jl:320
  [7] buildfilter
    @ C:\Users\glxbb\.julia\packages\NeXLSpectrum\eVZfj\src\filter.jl:312 [inlined]
  [8] references(refs::Vector{NeXLSpectrum.ReferencePacket}, det::BasicEDS; ftype::Type{Float64}, filter::Type)
    @ NeXLSpectrum C:\Users\glxbb\.julia\packages\NeXLSpectrum\eVZfj\src\reference.jl:156
  [9] references
    @ C:\Users\glxbb\.julia\packages\NeXLSpectrum\eVZfj\src\reference.jl:146 [inlined]
 [10] references(refs::Vector{NeXLSpectrum.ReferencePacket}, fwhm::Int64; kwargs::@Kwargs{})
    @ NeXLSpectrum C:\Users\glxbb\.julia\packages\NeXLSpectrum\eVZfj\src\reference.jl:166
 [11] top-level scope
    @ REPL[255]:1

As I post just above I'm wondering how to use the emsa standards when map has energy offset