A dashboard for exploring MA-XRF data

Given our precomputed .datastack files present in our Nextcloud storage we can now visualize our MA-XRF data in a Dashboard viewer after downloading them once.

Warning

The total size of these 16 files is pretty big (56.5 GB) and downloading on my computer took about half an hour!

Code
from fairdatanow import DataViewer
import matplotlib.pyplot as plt 
from matplotlib.patches import Rectangle 
import os
import numpy as np 
import moseley as mos 
import maxrf4u as mx
import re 
import imageio as iio

from falnama import get_roi_ims
Code
configuration = {
    'url': "https://laboppad.nl/falnama-project", 
    'user':    os.getenv('NC_AUTH_USER'),
    'password': os.getenv('NC_AUTH_PASS')
}
filters = {'extensions': ['.datastack']} 

dv = DataViewer(configuration, **filters)
datastack_files = dv.download_filtered()
Ready with downloading 16 selected remote files to local cache: /home/frank/.cache/fairdatanow                                                                      

Ok, let’s start our exploration by plotting a thumbnail overview of all scanned areas.

imvis_list = []
extent_list = []
for datastack_file in datastack_files: 
    ds = mx.DataStack(datastack_file)
    im = ds.read('imvis_reg')
    extent = ds.read('imvis_extent')
    imvis_list.append(im)
    extent_list.append(extent)

object_nums = [re.sub(r'.*(WM-71803-\d\d).*', r'\1', dsf) for dsf in datastack_files]

fig, axs = plt.subplots(ncols=8, nrows=2, figsize=[15, 6])
axs = axs.flatten()
for i, [im, ax, extent] in enumerate(zip(imvis_list, axs, extent_list)): 
    ax.imshow(im, extent=extent)
    ax.set_title(f'[{i}] {object_nums[i]}', fontsize=8)
    ax.axis('off')

Solving the sulfur situation

It came as a shock to me today that my two sulfur maps (i.e. a slice map and an NMF map) both are worse than the Bruker map.

filters = {'search': '.*Results/04.*71803-35_S.tif', 'use_regex': True}   
dv = DataViewer(configuration, **filters)
bruker_sulfur_tif = dv.download_filtered()[0]
Ready with downloading 1 selected remote files to local cache: /home/frank/.cache/fairdatanow                                                                      
bruker_S_map = plt.imread(bruker_sulfur_tif)
fig, ax = plt.subplots()
ax.imshow(bruker_S_map)
ax.set_title('Bruker S map');

Parse datastack

A first step now is to read relevant datasets from the datastack file and manually select regions of interest for page 35 (index 15)…

n = 15
datastack_file = datastack_files[n]
ds = mx.DataStack(datastack_file)

object_num = object_nums[n]

cube = ds.read('maxrf_cube', compute=False)
imvis_highres = ds.read('imvis_reg_highres')
imvis = ds.read('imvis_reg')
extent = ds.read('imvis_extent')
x_keVs = ds.read('maxrf_energies')
y_maxspectrum = ds.read('maxrf_maxspectrum')

nmf_elementmaps = ds.read('nmf_elementmaps') 

#ppa = mx.Peak_Pattern_Atlas(datastack_file=datastack_file, tube_keV=23)
ppa = mos.PeakPatternAtlas(x_keVs=x_keVs) 


ds.tree()
Ready building Peak Pattern Atlas!                                                                   
/
├── compton_peak_energy (1,) float64
├── hotmax_baselines (31, 4096) float64
├── hotmax_noiselines (31, 4096) float64
├── hotmax_peak_idxs_flat (34,) int64
├── hotmax_peak_idxs_list (31, 2) int64
├── hotmax_spectra (31, 4096) float32
├── hotmax_spots (31, 2) int64
├── hotmax_subpeak_idxs_list (31, 18) int64
├── imvis_extent (4,) int64
├── imvis_reg (617, 357, 3) uint8
├── imvis_reg_highres (8256, 4777, 3) uint8
├── maxrf_cube (617, 357, 4096) float32
├── maxrf_energies (4096,) float64
├── maxrf_maxspectrum (4096,) float32
├── maxrf_sumspectrum (4096,) float64
├── nmf_atomnums (18,) int64
├── nmf_elementmaps (18, 617, 357) float32
├── nmf_gausscomponents (35, 4096) float32
├── nmf_peakmaps (35, 617, 357) float32
└── nmf_peaks2elements_matrix (18, 34) float32
/home/frank/.cache/fairdatanow/falnama-project/OPENDATA/maxrf/datastacks/WM-71803-35_400_500_50.datastack:

# broken... 
# mx.get_instrument_pattern(datastack_file)

Regions of interest and their elements

fig, axs = plt.subplots(ncols=5, nrows=2, figsize=[20, 15]) axs = axs.flatten() for ax in axs: ax.imshow(imvis_highres, extent=extent)

xylims_list = [[int(ax.get_xlim()[0]), int(ax.get_xlim()[1]), int(ax.get_ylim()[0]), int(ax.get_ylim()[1])] for ax in axs]

xylims_list = [[29, 51, 106, 84], [39, 88, 136, 91], [154, 183, 242, 208], [275, 296, 249, 227], [40, 59, 289, 269],
               [72, 124, 441, 397], [146, 168, 462, 443], [139, 161, 393, 370], [245, 268, 454, 433], [197, 219, 300, 278]]

roi_ims = get_roi_ims(datastack_file, xylims_list)
fig, axs = plt.subplots(ncols=len(roi_ims)+1, figsize=[20, 4], squeeze=True)

edgecolors = ['violet', 'cyan', 'blue', 'red', 'white', 'orange', 'brown', 'green', 'maroon', 'black']

axs[0].imshow(imvis_highres, extent=extent)

for i, xylims in enumerate(xylims_list):
    x0, x1, y0, y1 = xylims
    axs[0].add_patch(Rectangle([x0, y1], x1 -x0, y0 - y1, linewidth=1, edgecolor=edgecolors[i], facecolor='none'))

for i, [ax, im] in enumerate(zip(axs[1:], roi_ims)): 
    
    ax.imshow(im)
    ax.set_title(f'[{i}] {edgecolors[i]}')

Matplotlib houtje-touwtje dashboard

Here is the latest still very rough code for an ROI tabs dashboard.

plt.close('all')

Ok, we now need some pretty advanced ROI hotmax functionality. Ehm, let’s build that here from scratch…

import falnama as fn
###### BUILD DASHBOARD FIGURES CANVAS ################### r

roi_elements_list = [['Pb', 'K', 'Ca', 'Fe'], #0 We see clear peaks for K and Ca in the spectrum and slice maps why is Ca not showing clearly in the Ca (NMF) map?
                     ['Pb_1', 'K', 'Ca', 'As_1', 'S'], #1
                     ['Pb', 'Fe', 'Ca', 'Hg'], #2
                     ['Pb', 'Pb_1', 'S'], #3
                     ['Pb', 'Fe', 'Ca'], #4
                     ['Pb', 'As', 'S', 'Fe'], #5 Here we see that the NMF maps are more clear the the alpha slices
                     ['Au', 'Pb', 'As'], #6 Clearly gold on red lead. The Pb alpha slice map is confusing, no Hg present in the Au  
                     ['Pb', 'As', 'S', 'Au'], #7 The yellow robe is clearly orpiment. Here again the the overlapping Pb alpha with As is confusing 
                     ['As', 'S', 'Au', 'Pb'], #8 According to the NMF map there should be no lead present in the green robe 
                     ['Pb']] #9

roi_cubes = [cube[l:k, i:j].compute() for i, j, k, l in xylims_list] 
roi_ims = get_roi_ims(datastack_file, xylims_list)
roi_n = 3
ncols = len(roi_elements_list[roi_n])
hspace = 0.02
wspace = 0.02 
xylims = xylims_list[roi_n]
l, r, b, t = xylims 
w = r - l 
h = b - t 
roi_extent = [0, w, h, 0]
 
fig = plt.figure(layout='constrained', figsize=(15, 15))
topfigs, midfig, botfig = fig.subfigures(nrows=3, hspace=hspace, height_ratios=[2, 1, 2], squeeze=True) 
topleftfig, topmidfig, toprightfig = topfigs.subfigures(ncols=3, wspace=wspace, width_ratios=[1, 1, 2], squeeze=True) 

ax_00 = topleftfig.subplots()
ax_01 = topmidfig.subplots()
ax_020, ax_021 = toprightfig.subplots(nrows=2, squeeze=True, sharex=True)

spectra_axs = midfig.subplots(ncols=ncols, squeeze=True, sharex=True) #, sharey=True) 

maps_axs = botfig.subplots(nrows=2, ncols=ncols, sharex=True, sharey=True)


########## TOP #################

# RGB overview and RGB ROI 
ax_00.imshow(imvis_highres, extent=extent)
for i, xylims in enumerate(xylims_list):
    x0, x1, y0, y1 = xylims
    ax_00.add_patch(Rectangle([x0, y1], x1 -x0, y0 - y1, linewidth=1, edgecolor=edgecolors[i], facecolor='none'))
    ax_00.annotate(f'ROI {i}', xy=[x0, y1], xytext=[0,5], textcoords='offset points', color='w', fontsize=8)
    
ax_00.set_title(object_num)


ax_01.imshow(roi_ims[roi_n], extent=roi_extent) 
ax_01.set_title(f'ROI {roi_n}')

# ROI spectral plots and peek pattern atlas 

# compute roi cube (core) maxspectrum 
hotmax_idxs, spectra = get_hotmax_spectra(roi_cubes[roi_n])
roi_cube = roi_cubes[roi_n]
roi_h, roi_w, roi_d = roi_cube.shape
roi_flat = roi_cube.reshape([-1, roi_d])
roi_maxspectrum = roi_flat.max(axis=0) # actually this is the roi core max spectrum  
roi_elements = roi_elements_list[roi_n]

ax_020.plot(x_keVs, roi_maxspectrum, label='ROI maxspectrum')
ax_020.plot(x_keVs, y_maxspectrum, color='r', alpha=0.2, label='full page maxspectrum')
ax_020.set_xlim([-1, 25])
ax_021.set_xlabel('energy [keV]')
#ax_020.legend() # need to fix blocking of view 

# plot peak pattern atlas 
select_elems = roi_elements_list[roi_n]
select_elems = [re.sub('_.*', '', elem) for elem in select_elems]

ppa.plot_atlas(ax=ax_021)


for i, element in enumerate(roi_elements): 

    # based on element now compute element slice map and related hotmax properties 
    slice_map, slice_idx = get_slice_map(datastack_file, element)
    hx, hy, roi_x, roi_y, hotmax_spectrum = get_roi_slice_hotmax(datastack_file, slice_idx, xylims_list[roi_n]) 

    # add hotmax locations to TOP image plots 
    ax_01.scatter(roi_x, roi_y)

    # add slice line to spectral plots 
    #spectra_axs[i].axvline(x_keVs[slice_idx], color='r')

    # add hotmax locations to slice maps 
    maps_axs[0, i].scatter(hx, hy, color='w', marker='+')
    
    # add hotmax spectra and element shadows to mid figure
    spectra_axs[i].plot(x_keVs, hotmax_spectrum)
    spectra_axs[i].scatter(x_keVs[slice_idx], hotmax_spectrum[slice_idx], marker='s', edgecolor='r', facecolor='none')
    #plot_element_shadow(element, x_keVs, flip=1, ax=spectra_axs[i], y=hotmax_spectrum)

    spectra_axs[i].plot(x_keVs, y_maxspectrum, color='r', alpha=0.15)
    spectra_axs[i].set_xlim([-1, 25])
    spectra_axs[i].set_title(f'{element} (slice {x_keVs[slice_idx]:.1f} keV)', fontsize=8)
    
    plot_element_lines(element, x_keVs, ax=spectra_axs[i], color='r', y=hotmax_spectrum)
    
    # add NMF map for element 
    nmf_map =  pick_nmf_elementmap(datastack_file, nmf_elementmaps, element) 
    
    maps_axs[0, i].imshow(slice_map)
    maps_axs[0, i].set_title(f'{element} (slice {x_keVs[slice_idx]:.1f} keV)', fontsize=8)
    maps_axs[1, i].imshow(nmf_map)
    maps_axs[1, i].set_title(f'{element} (NMF)', fontsize=8)

    # add ROI rectangle 
    x0, x1, y0, y1  = xylims_list[roi_n]
    maps_axs[0, i].add_patch(Rectangle([x0, y1], x1 -x0, y0 - y1, linewidth=1, edgecolor=edgecolors[roi_n], facecolor='none'))
    maps_axs[0, i].annotate(f'ROI {roi_n}', xy=[x0, y1], xytext=[0,5], textcoords='offset points', color='w', fontsize=10)
    maps_axs[1, i].add_patch(Rectangle([x0, y1], x1 -x0, y0 - y1, linewidth=1, edgecolor=edgecolors[roi_n], facecolor='none'))
    maps_axs[1, i].annotate(f'ROI {roi_n}', xy=[x0, y1], xytext=[0,5], textcoords='offset points', color='w', fontsize=10)

Can we fit a lead layer thickness?

Ok, Let’s extract a global hotmax spectrum for Pb_1

Pb_slice_map, Pb_slice_idx = get_slice_map(datastack_file, 'Pb_1')
yi, xi = np.argwhere(Pb_slice_map == Pb_slice_map.max())[0] 
Pb_hotmax = cube[yi, xi, :].compute()
fig, axs = plt.subplots(ncols=2, figsize=[15, 10], width_ratios=[1, 3])
ax0, ax1 = axs.flatten()

ax0.imshow(Pb_slice_map)
ax0.scatter(xi, yi, color='w', marker='+')
ax0.set_title('Pb_1')

ax1.fill_between(x_keVs, Pb_hotmax / Pb_hotmax.max(), color='r', label='Pb hotmax')

h_list = [0.0001, 0.001, 0.005, 0.01, 0.05, 0.1] 

colors = cm.viridis(np.linspace(0, 1, len(h_list)))

for i, h_mm in enumerate(h_list): 
    
    Pb_xrf = mos.ElementXRF('Pb', h_mm=h_mm, excitation_energy_keV=30, x_keVs=x_keVs, std=0.02)
    Pb_x, Pb_y = Pb_xrf.spectrum_xy
    
    ax1.plot(Pb_x, 100*Pb_y, color=colors[i], label=f'Pb ({h_mm})')
    ax1.legend()
    ax1.set_xlim([0, 18]);

I need to have access to the peak indexes…

Pb_xrf.get_pattern_dict()
{'elem': 'Pb',
 'atomic_number': 82,
 'alpha_keV': 12.609790734147028,
 'name': 'Lead',
 'peaks_xy': array([[1.26097907e+01, 3.57595202e-02],
        [1.05442700e+01, 3.31442037e-02],
        [1.47655527e+01, 3.64028035e-03],
        [1.51666247e+01, 1.29834767e-03],
        [9.18062517e+00, 1.29593281e-03],
        [2.37242804e+00, 5.87008945e-04],
        [1.13464139e+01, 4.67320502e-04]]),
 'alpha_escape_keV': 10.869790734147028}

FUNCTIONS


source

plot_roi_dashboard

 plot_roi_dashboard (datastack_file, xylims_list, roi_elements_list)

Create tabbed plot for roi spectra.


source

plot_roi_peak_patterns

 plot_roi_peak_patterns (roi_cube, prominence=2, ax=None)

source

get_roi_ims

 get_roi_ims (datastack_file, xylims_list)

Get high resolution roi images for xylims_list.


source

get_slice_map

 get_slice_map (datastack_file, element)

*Compute simple single channel maps from datastack_file cube for element alpha channel.

Add ’_1’ to element string for beta peak slice map.

Returns: slice_map, slice_idx*


source

pick_nmf_elementmap

 pick_nmf_elementmap (datastack_file, nmf_elementmaps, element)

Faster get nmf element map


source

get_nmf_map

 get_nmf_map (datastack_file, element)

Read computed NMF element map for element from datastack.


source

plot_roi_maps

 plot_roi_maps (roi_cubes, xylims_list, roi_n, elements,
                y_maxspectrum=None, exc_keV=23)

*Plot roi alpha maps for list of elements.

If y_maxspectrum is specified normalize with maximum at channel index.*


source

get_roi_maps

 get_roi_maps (roi_cube, elements, x_keVs, exc_keV=23)

*Slice roi cube at alpha channel index positions for each element in elements.

Returns: alpha_idxs, roi_maps, roi_vmax_list*


source

plot_element_lines

 plot_element_lines (element, x_keVs, y=None, ax=None, exc_keV=23,
                     color=None)

Add emission lines to spectral plot.


source

get_element_lines

 get_element_lines (element, x_keVs, exc_keV=23)

*Get sorted element peak idxs.

Returns: peak_idxs_sorted*


source

plot_hotmax_spectra

 plot_hotmax_spectra (spectra, x_keVs, hotmax_idxs, elems=None, xlim=[-1,
                      22], elem_colors={'Pb': 'k', 'Au': 'orange', 'Fe':
                      'brown', 'Ca': 'g', 'K': 'b', 'As': 'magenta', 'Na':
                      'blue', 'Ni': 'red', 'Cu': 'green', 'Ti': 'blue',
                      'Zn': 'red'}, y_maxspectrum=None)

*Plot roi hotmax spectra and add emission lines for list of elements elems.

Returns: axs*


source

get_hotmax_spectra

 get_hotmax_spectra (roi_cube)

*Get hotmax spectra for roi_cube.

Returns: hotspot_idxs,spectra`*


source

plot_roi_peak_patterns

 plot_roi_peak_patterns (roi_cube, prominence=2, ax=None)

source

get_roi_slice_hotmax

 get_roi_slice_hotmax (datastack_file, slice_idx, xylims)

*Find full map and roi coordinates (x, y, roi_x, roi_y) of slice map maximum at channel slice_idx within ROI xylims.

Returns: x, y, roi_x, roi_y, hotmax_spectrum*