Sparse to sparse integration?

[jadball]

Hello all,

At ID11 we are collecting many millions of extremely sparse Eiger frames currently. In the past with Frelon data, we have just integrated dense -> dense images and saved them, but the data volumes were much lower. In ImageD11, we have a SparseFrame format which is basically just a COO format (I think!). Would it be possible in some way to feed this into pyFAI to convert sparse detector pixels into sparse 2d (tth, chi) values?

[kif]

Hi James,

If you are just interested in the 2th and chi values, those array are always initialized and readily available. Have a look at ai._cached_arrays ... there is plenty of things hidden there inside.

But if you want to perform the integration, such integrator does not yet exists but it should not be too complicated to adapt the sparse-CSC engine (probably the Cython implementation is the most effective one to start with) to be able to work on such COO data. The code should be single threaded (per frame) but can be parallelized between frames.

Basically, the COO format is a list of (x, y, signal) tuples. The first thing is to convert the (x, y) into an index position idx = x + y*width

The CSC-sparse array is composed of 3 arrays called data, index and indptr

Initialize 2 output containers called sum_signal and sum_normalization.

*Then loop over all pixels in your sparse frame, calculate its idx value.

• Then loop over i taking all values from indptr[idx] to indptr[idx+1]
  • Append signal * data[i] to sum_signal[index[i]]
  • Append data[i] * norm[i] to sum_normalization[index[i]] where norm is is product of the polarization, solid-angle and flatfield of the detector.

The average value is then the ratio between sum_signal/sum_normalization

[jonwright]

There is an example in https://github.com/jonwright/bslz4_to_sparse/blob/main/test/test_dot.py

It calls on a bslz4_to_sparse.chunk2sparseCSC object that is defined in bslz4_to_sparse/init.py

This is installed in the ESRF conda. It does need a more user-friendly wrapper....!

[jadball]

Thanks both for taking a look.
I experimented with this a bit last week and made something that sort of works.
It doesn't produce results that exactly match the regular integration in terms of pixel-by-pixel values, but it's certainly close.
@kif what do you think about my normalisation attempt?

import pyFAI
from pyFAI.integrator import load_engines
import fabio
import numpy as np
import scipy.sparse
from ImageD11.sinograms import dataset
import ImageD11.sparseframe


# prepare files needed
ponifile = 'path/to/poni'
maskfile = 'path/to/mask'
dsfile = 'path/to/ImageD11/ds'

# prepare pyFAI integrator
ai = pyFAI.load(ponifile)
msk = fabio.open(ds.maskfile).data

npt_rad = 3000
npt_azim = 5000

engine = ai.setup_sparse_integrator(algo='CSR', shape=ai.detector.shape, npt=(npt_rad,npt_azim), mask=msk, unit='2th_deg', split='full')

# get the lut of the engine
lut = (engine._data, engine._indices, engine._indptr)

# convert to CSR matrix format
csr_transformer = scipy.sparse.csr_array(lut)

# now we make the sparse CSR from the ImageD11 data
ds = ImageD11.sinograms.dataset.load(dsfile)
sparse_scan = ImageD11.sparseframe.SparseScan(ds.sparsefile, ds.scans[0])

# sparse_scan is a list of COO over (n_frames, rows, cols)
# convert (rows, cols) to flat index of (total_pixels,)
width = msk.shape[1]
flat_indices = (sparse_scan.row.astype(np.int64) * width + sparse_scan.col).astype(np.int64)

# nnz is "active pixels per frame" in ImageD11 COO list
# convert to indptr
indptr = np.zeros(len(sparse_scan.nnz) + 1, dtype=np.int32)
np.cumsum(sparse_scan.nnz, out=indptr[1:])

# shape of data CSR - (n_frames, total_pixels)
total_pixels = msk.size
shape = (len(sparse_scan.nnz), total_pixels)

# make scipy CSR
csr_data = scipy.sparse.csr_array((sparse_scan.intensity, flat_indices, indptr), shape=shape)

# perform the transform - very very fast!!!
# this happens over all frames at once.
result_sparse = csr_transformer.dot(csr_data.T)

# normalise the transform
# avoiding dividing by zero weights
# dot-product array of ones
weights_vec = csr_transformer.dot(np.ones((msk.size, 1)))
inv_w = np.zeros_like(weights_vec)
np.divide(1.0, weights_vec, out=inv_w, where=weights_vec != 0)
norm_sparse = result_sparse.multiply(inv_w)

[jonwright]

It should be just to change weights_vec to use the polarisation and solid angle arrays (masked) that you have in the ai.

[kif]

I agree with Jon about the normalization, it is the denominator of eq1 in
https://journals.iucr.org/j/issues/2025/01/00/jo5109/index.html#FD1

Everything is available in the cached_array's dictonary.

Beside this, if you don't want to densify your image (stored as COO), you have to work with a linear transformation which is stored as a CSC sparse matrix, not as a CSR sparse matrix.

• The former (CSC) offers you an access pixel per pixel, useful with your COO image, but incompatible with parallelization.
• The later (CSR) offers you an access bin per bin, useful when parallelizing the processing, since writing does not require locks.

[jadball]

When I initialise an integrator and setup a sparse integrator, I actually don't see much in the cached arrays. I am now interested in the 2th, chi values given my specified units and npts, but I don't seem to be able to easily access them:

ai = pyFAI.load(ponifile)

engine = ai.setup_sparse_integrator(algo='CSR', shape=ai.detector.shape, npt=(npt_rad,npt_azim), mask=msk, unit='2th_deg', split='full')

# get the lut of the engine
lut = (engine._data, engine._indices, engine._indptr)

print(ai._cached_array.keys())
# gives:
dict_keys(['2th_corner'])

[jonwright]

Just integrate a frame and everything gets populated?

[kif]

Calculate the arrays like this: * `ai.array_from_unit(unit="2th_deg")` * `ai.array_from_unit(unit="chi_deg")` Cheers,

[jadball]

Thanks - a related comment - once I apply the transformation, I'm left with a sparse result with some sort of data layout. I know which axis is radial and which is azimuthal (different shapes) but I'm not aware about how the axis indices align with the units. In my testing, it seems that the result images are tth-aligned (along the i axis) but chi is decreasing as the j axis increases. Is that expected?

In the end I did this to get flat arrays:

tth_axis = pyFAI.units.to_unit(engine.unit).scale * engine.bin_centers0
chi_axis = np.degrees(engine.bin_centers1)

[kif]

To answer your question: pyFAI does most of the operation in SI-units (exception: q is in inverse nm), so engine.bin_centers[01] are in internal units and need to be converted to the user's requested unit by multiplication of the scale factor of the unit.

To come back on different sparse matrix representation (for images or linear transformations):
If one COO frame is a triplet of 1D vector (x,y,I):

pos = x+y*witdth
tth = ai.array_from_unit(unit="2th_deg")[pos]
chi = ai.array_from_unit(unit="chi_deg")[pos]

Then your pixels are at tth and chi values, the intensity I remains the same.

But this is not azimuthal integration because neither the rebinning, nor the averaging is performed.

I could write for you a rebinning engine taking as input a COO sparse matrix, but are you actually using this format ? The main issue is that this is not a bug, so it has no priority.

Performing the double loop in python / cython / numba is the way to go (probably not vectorisable in a simple manner, maybe investigate einsum, it has magic features).