Skip to content

Example 2: A powder pattern from the same pass: pyFAI's sparse matrix in the decoder

chunk2sparseCSC(mask, csc, dtype) is as in 0.0.15 and 0.0.16: chunk2sparse plus one more output. While a chunk is decoded, every pixel is also multiplied into a sparse matrix (bins x pixels) and added into the bins it falls in. pyFAI's CSC integration engine is exactly such a matrix, so the powder pattern comes out of the decompression with no dense frame. A call returns

npx            the number of pixels above the cut
values, indices  those pixels (values and their flat positions; the first npx entries)
powder_sum     for every bin, the sum over its pixels of weight x pixel value: pyFAI's
               sum_signal, float64

pyFAI divides powder_sum by sum_normalization, the same sum for an image of ones with the corrections (solid angle, polarisation, flat field) folded in, to get the intensity. sum_normalization does not depend on the data: it is computed once. The matrix is used as it is given, as in the released versions. The last part of this example shows what is new in 0.0.22: fixed_point_matrix prepares the matrix for exact integer sums.

Source: examples/ex2_pyfai_powder.py, run as python examples/ex2_pyfai_powder.py, or download the notebook (figures and output as run when this page was built).

Output, from running it when this page was built, on a synthetic stand-in for your files (powder rings with noise on an Eiger-like detector, tools/standin.py; with your own data the numbers are yours):

import h5py
import hdf5plugin  # registers the bitshuffle filter, so that h5py can decompress too
import numpy as np
import pyFAI
from pyFAI.method_registry import IntegrationMethod

import bslz4_to_sparse

h5filename = "eiger_0000.h5"  # your file: bitshuffle-compressed, one frame per chunk
h5dataset = "/entry_0000/ESRF-ID11/eiger/data"  # the (nframes, ni, nj) dataset inside it
mask_filename = (
    "mask.h5"  # your detector mask: a dataset "mask", pyFAI's convention (non-zero = masked)
)
poni_filename = "calib.poni"  # your pyFAI geometry

with h5py.File(h5filename, "r") as h5:
    ds = h5[h5dataset]
    dtype = ds.dtype
    frame0 = ds[0]
    chunk0 = ds.id.read_direct_chunk((0, 0, 0))[1]
with h5py.File(mask_filename, "r") as h5:
    bad = h5["mask"][()] != 0  # the masked pixels
mask = (~bad).astype(np.uint8)

pyFAI makes its sparse matrix once, when it first integrates with a sparse method (bbox pixel splitting, CSC layout). Integrating an image of ones also gives sum_normalization.

ai = pyFAI.load(poni_filename)
method = IntegrationMethod.select_one_available(("bbox", "csc", "cython"), dim=1)
ones = ai.integrate1d(
    np.ones(frame0.shape, np.float32), 400, unit="q_nm^-1", method=method, mask=bad
)
engine = ai.engines[method].engine  # pyFAI's CSC matrix
sum_normalization = np.asarray(ones.sum_normalization).ravel()

WARNING:silx.opencl.common:Unable to import pyOpenCl. Please install it from: https://pypi.org/project/pyopencl

Our object: decode, collect the sparse pixels and multiply into the matrix, all in one pass.

integ = bslz4_to_sparse.chunk2sparseCSC(mask, engine, dtype)
cut = 3
npx, (values, indices), powder_sum = integ(chunk0, cut)
print("pipeline dot:", bslz4_to_sparse.describe(integ.pipeline)["dot"])
print("sparse pixels above the cut of %d: %d" % (cut, npx))
print("powder_sum: %d bins of %s" % (powder_sum.size, powder_sum.dtype))
sum_signal = powder_sum

pipeline dot: padded-avx2
sparse pixels above the cut of 3: 7023
powder_sum: 400 bins of float64

the same frame, integrated by pyFAI from the decompressed pixels

ref = ai.integrate1d(frame0, 400, unit="q_nm^-1", method=method, mask=bad)
assert npx == np.count_nonzero((frame0 > cut) & ~bad)
signal = np.asarray(ref.sum_signal).ravel()
print(
    "largest difference from pyFAI's sum_signal: %.3g (largest bin %.3g)"
    % (abs(sum_signal - signal).max(), signal.max())
)
assert np.allclose(sum_signal, signal, rtol=1e-9, atol=0)

largest difference from pyFAI's sum_signal: 0 (largest bin 5.61e+03)

intensity = np.divide(
    sum_signal, sum_normalization, out=np.zeros_like(sum_signal), where=sum_normalization > 0
)
assert np.allclose(intensity, ref.intensity, rtol=1.2e-7, atol=0)
print("intensity = powder_sum / sum_normalization equals pyFAI's integrate1d")

intensity = powder_sum / sum_normalization equals pyFAI's integrate1d

A picture: why powder_sum is divided by sum_normalization. The raw sum in a bin grows with the number of pixels that fall in it (more pixels at large radius, none in the module gaps); the normalisation counts them (weighted), and the ratio is the intensity.

import matplotlib.pyplot as plt

q = ones.radial
fig, axes = plt.subplots(1, 3, figsize=(12, 3.5))
axes[0].plot(q, sum_signal)
axes[0].set_title("powder_sum: sum of pixel values in each bin")
axes[1].plot(q, sum_normalization)
axes[1].set_title("sum_normalization: pixels in each bin")
axes[2].plot(q, intensity, label="powder_sum / sum_normalization")
axes[2].plot(q, ref.intensity, "--", label="pyFAI integrate1d")
axes[2].set_title("intensity")
axes[2].legend()
for ax in axes:
    ax.set_xlabel("q (nm$^{-1}$)")
plt.tight_layout()
plt.show()

figure 1

Only the powder is wanted? Then no pixel needs collecting: a cut at the dtype's maximum finds none (and howmany=0 would keep none of any that were found).

nothing = np.iinfo(dtype).max
only_powder = bslz4_to_sparse.chunk2sparseCSC(mask, engine, dtype)
npx, _, powder_only = only_powder(chunk0, nothing)
assert npx == 0 and np.array_equal(powder_only, powder_sum)
print("with the cut at the dtype maximum: no pixels, the same powder_sum")

with the cut at the dtype maximum: no pixels, the same powder_sum

2D (radial x azimuth): the matrix has nrad * nazim rows, in radial-major order

method2 = IntegrationMethod.select_one_available(("bbox", "csc", "cython"), dim=2)
ref2 = ai.integrate2d(frame0, 100, 18, unit="q_nm^-1", method=method2, mask=bad)
integ2 = bslz4_to_sparse.chunk2sparseCSC(mask, ai.engines[method2].engine, dtype)
_, _, powder_sum2 = integ2(chunk0, nothing)
cake = powder_sum2.reshape(100, 18).T  # pyFAI's (azimuth, radial) layout
assert np.allclose(cake, ref2.sum_signal, rtol=1e-9, atol=0)
print("2D cake", cake.shape, "matches pyFAI's integrate2d sum_signal")

2D cake (18, 100) matches pyFAI's integrate2d sum_signal

New in 0.0.22: exact sums. fixed_point_matrix folds the mask into the matrix and stores its weights as fixed-point integers; chunk2sparseCSC given it returns the int64 sums of integer weight x pixel, and dividing by powder_scale (an integer power of two) turns them into float. Integer pixels only; for 32-bit pixels the weights get fewer bits, unless maxval says how large the data can be.

fixed = bslz4_to_sparse.fixed_point_matrix(engine, mask, dtype)
exact = bslz4_to_sparse.chunk2sparseCSC(mask, fixed, dtype)
_, _, powder_int = exact(chunk0, nothing)
print("fixed point: powder %s, scale %d" % (powder_int.dtype, fixed.scale))
assert np.allclose(powder_int * (1.0 / exact.powder_scale), powder_sum, rtol=1e-6, atol=0)

fixed point: powder int64, scale 2147483648

Generated by tools/generate_docs.py from bslz4_to_sparse 0.0.22a2 (git v0.0.22-13-g610ed27, sources sha256 6929bd75ab7a).