Filtering signal in azimuthal space

Usually, diffraction signal presents a polar symmetry, this means all pixel with the same azimuthal angle (χ) have similar intensities. The best way to exploit this is to take the mean, what is called azimuthal average. But the average is very sensitive to outlier, like gaps, missing pixels, shadows, cosmic rays or reflection coming from larger crystallite. In this tutorial we will see two alternative ways to remove those unwanted signal and focus on the majority of pixels: sigma clipping and median filtering.

import os
os.environ["PYOPENCL_COMPILER_OUTPUT"]="0"

%matplotlib inline
from matplotlib.pyplot import subplots
from pyFAI.gui import jupyter
import numpy, fabio, pyFAI
from pyFAI import benchmark
from pyFAI.test.utilstest import UtilsTest

ai = pyFAI.load(UtilsTest.getimage("Pilatus6M.poni"))
img = fabio.open(UtilsTest.getimage("Pilatus6M.cbf")).data
fig, ax = subplots(1, 2)
jupyter.display(img, ax=ax[1])
jupyter.plot1d(ai.integrate1d(img, 1000), ax=ax[0])
ax[1].set_title("With a few Bragg peaks")
pass

WARNING:pyFAI.gui.matplotlib:matplotlib already loaded, setting its backend may not work

Azimuthal sigma-clipping

The idea is to discard pixels which look like outliers in the distribution of all pixels contributing to a single azimuthal bin. It requires an error model like poisson but it has been proven to be better to use the variance in the given azimuthal ring. All details are available in this publication: https://doi.org/10.1107/S1600576724011038 also available at https://doi.org/10.48550/arXiv.2411.09515

fig, ax = subplots(1, 2)
jupyter.display(img, ax=ax[1])
jupyter.plot1d(ai.sigma_clip(img, 1000, error_model="hybrid", method=("no", "csr", "cython")), ax=ax[0])
ax[1].set_title("With a few Bragg peaks")
ax[0].set_title("Sigma_clipping")
pass

Of course, sigma-clip takes several extra parameters like the number of iterations to perform, the cut-off, the error model, … There are alo a few limitation:

• The algorithm needs to be the CSR-sparse matrix multiplication: since several integration are needed, it makes no sense to use an histogram based algorithm.

• The algorithm is available with any implementation: Python (using scipy.saprse), Cython and OpenCL, and it runs just fine on GPU.

• Sigma-clipping is incompatible with any kind of pixel splitting: With pixel splitting, a single pixel can contribute to several azimuthal bins and discarding a pixel in one ring could disable it in the neighboring ring (or not, since bins are processed in parallel).

Sigma-clipping performances:

method = ["no", "csr", "cython"]

%%time
perfs_integrate_python = {}
perfs_integrate_cython = {}
perfs_integrate_opencl = {}
perfs_sigma_clip_python = {}
perfs_sigma_clip_cython = {}
perfs_sigma_clip_opencl = {}

for ds in pyFAI.benchmark.PONIS:
    ai = pyFAI.load(UtilsTest.getimage(ds))
    if ai.wavelength is None: ai.wavelength=1.54e-10
    img = fabio.open(UtilsTest.getimage(pyFAI.benchmark.datasets[ds])).data
    size = numpy.prod(ai.detector.shape)
    print(ds)
    print("     Cython")
    meth = tuple(method)
    nbin = max(ai.detector.shape)
    print("     * integrate ", end="")
    perfs_integrate_cython[size] = %timeit -o ai.integrate1d(img, nbin, method=meth)
    print("     * sigma-clip", end="")
    perfs_sigma_clip_cython[size] = %timeit -o ai.sigma_clip(img, nbin, method=meth, error_model="azimuthal")
    print("     Python")
    meth = tuple(method[:2]+["python"])
    print("     * integrate ", end="")
    perfs_integrate_python[size] = %timeit -o ai.integrate1d(img, nbin, method=meth)
    print("     * sigma-clip", end="")
    perfs_sigma_clip_python[size] = %timeit -o ai.sigma_clip(img, nbin, method=meth, error_model="azimuthal")

    print("     OpenCL")
    meth = tuple(method[:2]+["opencl"])
    print("     * integrate ", end="")
    perfs_integrate_opencl[size] = %timeit -o ai.integrate1d(img, nbin, method=meth)
    print("     * sigma-clip", end="")
    perfs_sigma_clip_opencl[size] = %timeit -o ai.sigma_clip(img, nbin, method=meth, error_model="azimuthal")

Pilatus1M.poni
     Cython
4.3 ms ± 6.97 μs per loop (mean ± std. dev. of 7 runs, 100 loops each)
4.93 ms ± 186 μs per loop (mean ± std. dev. of 7 runs, 100 loops each)
     Python
8.13 ms ± 1.67 μs per loop (mean ± std. dev. of 7 runs, 100 loops each)
     * sigma-clip

WARNING:pyFAI.opencl.azim_csr:Please upgrade to silx v2.2+

110 ms ± 221 μs per loop (mean ± std. dev. of 7 runs, 10 loops each)
     OpenCL
466 μs ± 459 ns per loop (mean ± std. dev. of 7 runs, 1,000 loops each)
2.26 ms ± 654 ns per loop (mean ± std. dev. of 7 runs, 100 loops each)
Pilatus2M.poni
     Cython
11.8 ms ± 1.22 ms per loop (mean ± std. dev. of 7 runs, 100 loops each)
26.1 ms ± 398 μs per loop (mean ± std. dev. of 7 runs, 100 loops each)
     Python
26 ms ± 35.4 μs per loop (mean ± std. dev. of 7 runs, 10 loops each)
     * sigma-clip

WARNING:pyFAI.opencl.azim_csr:Please upgrade to silx v2.2+

391 ms ± 2.07 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
     OpenCL
722 μs ± 320 ns per loop (mean ± std. dev. of 7 runs, 1,000 loops each)
5.01 ms ± 1.47 μs per loop (mean ± std. dev. of 7 runs, 100 loops each)
Eiger4M.poni
     Cython
20.5 ms ± 5.49 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
25.9 ms ± 6.18 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
     Python
53 ms ± 67.4 μs per loop (mean ± std. dev. of 7 runs, 10 loops each)
     * sigma-clip

WARNING:pyFAI.opencl.azim_csr:Please upgrade to silx v2.2+

860 ms ± 2.2 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
     OpenCL
1.18 ms ± 296 ns per loop (mean ± std. dev. of 7 runs, 1,000 loops each)
8.9 ms ± 1.05 μs per loop (mean ± std. dev. of 7 runs, 100 loops each)
Pilatus6M.poni
     Cython
26.4 ms ± 174 μs per loop (mean ± std. dev. of 7 runs, 1 loop each)
32 ms ± 1.18 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
     Python
75 ms ± 19.3 μs per loop (mean ± std. dev. of 7 runs, 10 loops each)
     * sigma-clip

WARNING:pyFAI.opencl.azim_csr:Please upgrade to silx v2.2+

1.24 s ± 5.02 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
     OpenCL
1.75 ms ± 960 ns per loop (mean ± std. dev. of 7 runs, 1,000 loops each)
12.1 ms ± 2.2 μs per loop (mean ± std. dev. of 7 runs, 100 loops each)
Eiger9M.poni
     Cython
48.2 ms ± 1.34 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
61.7 ms ± 1.3 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
     Python
140 ms ± 453 μs per loop (mean ± std. dev. of 7 runs, 10 loops each)
     * sigma-clip

WARNING:pyFAI.opencl.azim_csr:Please upgrade to silx v2.2+

2.5 s ± 38.9 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
     OpenCL
2.92 ms ± 931 ns per loop (mean ± std. dev. of 7 runs, 100 loops each)
21.4 ms ± 6.34 μs per loop (mean ± std. dev. of 7 runs, 10 loops each)
Mar3450.poni
     Cython
62.6 ms ± 1.07 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
69.6 ms ± 1.46 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
     Python
143 ms ± 800 μs per loop (mean ± std. dev. of 7 runs, 10 loops each)
     * sigma-clip

WARNING:pyFAI.opencl.azim_csr:Please upgrade to silx v2.2+

2.67 s ± 19.4 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
     OpenCL
3.33 ms ± 4.78 μs per loop (mean ± std. dev. of 7 runs, 100 loops each)
23.9 ms ± 5.11 μs per loop (mean ± std. dev. of 7 runs, 10 loops each)

---------------------------------------------------------------------------
KeyError                                  Traceback (most recent call last)
File <timed exec>:11

File /nobackup/scratch/kieffer/py310/lib/python3.10/site-packages/silx/utils/ExternalResources.py:216, in ExternalResources.getfile(self, filename)
    214 with open(fullfilename, mode="rb") as fd:
    215     h.update(fd.read())
--> 216 if h.hexdigest() != self.all_data[filename]:
    217     logger.warning(f"Detected corruped file {fullfilename}")
    218     self.all_data.pop(filename)

KeyError: 'Fairchild.edf'

fig, ax = subplots()
ax.set_xlabel("Image size (Mpix)")
ax.set_ylabel("Frames per seconds")
sizes = numpy.array(list(perfs_integrate_python.keys()))/1e6
ax.plot(sizes, [1/i.best for i in perfs_integrate_python.values()], label="Integrate/Python", color='green', linestyle='dashed', marker='1')
ax.plot(sizes, [1/i.best for i in perfs_integrate_cython.values()], label="Integrate/Cython", color='orange', linestyle='dashed', marker='1')
ax.plot(sizes, [1/i.best for i in perfs_integrate_opencl.values()], label="Integrate/OpenCL", color='blue', linestyle='dashed', marker='1')
ax.plot(sizes, [1/i.best for i in perfs_sigma_clip_python.values()], label="Sigma-clip/Python", color='green', linestyle='dotted', marker='2')
ax.plot(sizes, [1/i.best for i in perfs_sigma_clip_cython.values()], label="Sigma-clip/Cython", color='orange', linestyle='dotted', marker='2')
ax.plot(sizes, [1/i.best for i in perfs_sigma_clip_opencl.values()], label="Sigma-clip/OpenCL", color='blue', linestyle='dotted', marker='2')
ax.set_yscale("log")
ax.legend()

ax.set_title("Performance of Sigma-clipping vs integrate")

Text(0.5, 1.0, 'Performance of Sigma-clipping vs integrate')

The penalties is very limited in Cython, much more in Python.

The biggest limitation of sigma-clipping is its incompatibility with pixel-splitting, feature needed when oversampling, i.e. taking many more points than the size of the diagonal of the image. While oversampling is not recommended in general case (due to the cross-corelation between bins it creates), it can be a nessessary evil, especially when performing Rietveld refinement where 5 points per peaks are needed, resolution that cannot be obtained with the pixel-size/distance couple accessible by the experimental setup.

Median filter in Azimuthal space

The idea is to sort all pixels contibuting to an azimuthal bin and to average out all pixel between the lower and upper quantile. When those two thresholds are at one half, this filter provides actually the median. In order to be compatible with pixel splitting, each pixel is duplicated as many times as it contributes to different bins. After sorting fragments of pixels according to their normalization corrected signal, the cummulative sum of normalization is performed in order to detemine which fragments to average out.

ai = pyFAI.load(UtilsTest.getimage("Pilatus6M.poni"))
img = fabio.open(UtilsTest.getimage("Pilatus6M.cbf")).data

method = ["full", "csr", "cython"]
percentile=(40,60)
pol=0.99

fig, ax = subplots(1, 2)
jupyter.display(img, ax=ax[1])
jupyter.plot1d(ai.medfilt1d_ng(img, 1000, method=method, percentile=percentile, polarization_factor=pol), ax=ax[0])
ax[1].set_title("With a few Bragg peaks")
ax[0].set_title("Median filtering")
pass

Unlike the sigma-clipping, this median filter does not equires any error model; but the computationnal cost induced by the sort is huge. In addition, the median is very sensitive and requires a good geometry and modelisation of the polarization.

%%time
perf2_integrate_python = {}
perf2_integrate_cython = {}
perf2_integrate_opencl = {}
perf2_medfilt_python = {}
perf2_medfilt_cython = {}
perf2_medfilt_opencl = {}

for ds in pyFAI.benchmark.PONIS:
    ai = pyFAI.load(UtilsTest.getimage(ds))
    if ai.wavelength is None: ai.wavelength=1.54e-10
    img = fabio.open(UtilsTest.getimage(pyFAI.benchmark.datasets[ds])).data
    size = numpy.prod(ai.detector.shape)
    print(ds)
    print("     Cython")
    meth = tuple(method)
    nbin = max(ai.detector.shape)
    print("     * integrate   ", end="")
    perf2_integrate_cython[size] = %timeit -o ai.integrate1d(img, nbin, method=meth)
    print("     * medianfilter", end="")
    perf2_medfilt_cython[size] = %timeit -o ai.medfilt1d_ng(img, nbin, method=meth)
    print("     Python")
    meth = tuple(method[:2]+["python"])
    print("     * integrate   ", end="")
    perf2_integrate_python[size] = %timeit -o ai.integrate1d(img, nbin, method=meth)
    print("     * medianfilter", end="")
    perf2_medfilt_python[size] = %timeit -o ai.medfilt1d_ng(img, nbin, method=meth)

    print("     OpenCL")
    meth = tuple(method[:2]+["opencl"])
    print("     * integrate   ", end="")
    perf2_integrate_opencl[size] = %timeit -o ai.integrate1d(img, nbin, method=meth)
    print("     * medianfilter", end="")
    perf2_medfilt_opencl[size] = %timeit -o ai.medfilt1d_ng(img, nbin, method=meth)

Pilatus1M.poni
     Cython
4.51 ms ± 198 μs per loop (mean ± std. dev. of 7 runs, 100 loops each)
15.5 ms ± 129 μs per loop (mean ± std. dev. of 7 runs, 100 loops each)
     Python
11.3 ms ± 8.24 μs per loop (mean ± std. dev. of 7 runs, 100 loops each)
     * medianfilter

WARNING:pyFAI.opencl.azim_csr:Please upgrade to silx v2.2+

913 ms ± 5.34 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
     OpenCL
495 μs ± 285 ns per loop (mean ± std. dev. of 7 runs, 1,000 loops each)
4.43 ms ± 903 ns per loop (mean ± std. dev. of 7 runs, 100 loops each)
Pilatus2M.poni
     Cython
11.7 ms ± 669 μs per loop (mean ± std. dev. of 7 runs, 1 loop each)
82.2 ms ± 2.94 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
     Python
39.3 ms ± 157 μs per loop (mean ± std. dev. of 7 runs, 10 loops each)
     * medianfilter

WARNING:pyFAI.opencl.azim_csr:Please upgrade to silx v2.2+

3.08 s ± 35.9 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
     OpenCL
833 μs ± 473 ns per loop (mean ± std. dev. of 7 runs, 1,000 loops each)
12.6 ms ± 7.77 μs per loop (mean ± std. dev. of 7 runs, 100 loops each)
Eiger4M.poni
     Cython
19.8 ms ± 197 μs per loop (mean ± std. dev. of 7 runs, 1 loop each)
138 ms ± 1.45 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
     Python
75.8 ms ± 364 μs per loop (mean ± std. dev. of 7 runs, 10 loops each)
     * medianfilter

WARNING:pyFAI.opencl.azim_csr:Please upgrade to silx v2.2+

5.32 s ± 15.7 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
     OpenCL
1.49 ms ± 490 ns per loop (mean ± std. dev. of 7 runs, 1,000 loops each)
22.2 ms ± 7.54 μs per loop (mean ± std. dev. of 7 runs, 10 loops each)
Pilatus6M.poni
     Cython
30.1 ms ± 3.73 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
172 ms ± 1.39 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
     Python
108 ms ± 38.6 μs per loop (mean ± std. dev. of 7 runs, 10 loops each)
     * medianfilter

WARNING:pyFAI.opencl.azim_csr:Please upgrade to silx v2.2+

7.9 s ± 36.2 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
     OpenCL
2 ms ± 1.2 μs per loop (mean ± std. dev. of 7 runs, 100 loops each)
31.1 ms ± 11.8 μs per loop (mean ± std. dev. of 7 runs, 10 loops each)
Eiger9M.poni
     Cython
46.4 ms ± 858 μs per loop (mean ± std. dev. of 7 runs, 1 loop each)
266 ms ± 4.98 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
     Python
190 ms ± 418 μs per loop (mean ± std. dev. of 7 runs, 1 loop each)
     * medianfilter

WARNING:pyFAI.opencl.azim_csr:Please upgrade to silx v2.2+

13.1 s ± 47.7 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
     OpenCL
3.2 ms ± 799 ns per loop (mean ± std. dev. of 7 runs, 100 loops each)
57.1 ms ± 18.7 μs per loop (mean ± std. dev. of 7 runs, 10 loops each)
Mar3450.poni
     Cython
65.6 ms ± 185 μs per loop (mean ± std. dev. of 7 runs, 1 loop each)
291 ms ± 3.11 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
     Python
223 ms ± 1.54 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
     * medianfilter

WARNING:pyFAI.opencl.azim_csr:Please upgrade to silx v2.2+

14.9 s ± 44.5 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
     OpenCL
3.63 ms ± 2.76 μs per loop (mean ± std. dev. of 7 runs, 100 loops each)
65.8 ms ± 18.4 μs per loop (mean ± std. dev. of 7 runs, 10 loops each)

---------------------------------------------------------------------------
KeyError                                  Traceback (most recent call last)
File <timed exec>:11

File /nobackup/scratch/kieffer/py310/lib/python3.10/site-packages/silx/utils/ExternalResources.py:216, in ExternalResources.getfile(self, filename)
    214 with open(fullfilename, mode="rb") as fd:
    215     h.update(fd.read())
--> 216 if h.hexdigest() != self.all_data[filename]:
    217     logger.warning(f"Detected corruped file {fullfilename}")
    218     self.all_data.pop(filename)

KeyError: 'Fairchild.edf'

fig, ax = subplots()
ax.set_xlabel("Image size (Mpix)")
ax.set_ylabel("Frames per seconds")
sizes = numpy.array(list(perf2_integrate_python.keys()))/1e6
ax.plot(sizes, [1/i.best for i in perf2_integrate_python.values()], label="Integrate/Python", color='green', linestyle='dashed', marker='1')
ax.plot(sizes, [1/i.best for i in perf2_integrate_cython.values()], label="Integrate/Cython", color='orange', linestyle='dashed', marker='1')
ax.plot(sizes, [1/i.best for i in perf2_integrate_opencl.values()], label="Integrate/OpenCL", color='blue', linestyle='dashed', marker='1')
ax.plot(sizes, [1/i.best for i in perf2_medfilt_python.values()], label="Medfilt/Python", color='green', linestyle='dotted', marker='2')
ax.plot(sizes, [1/i.best for i in perf2_medfilt_cython.values()], label="Medfilt/Cython", color='orange', linestyle='dotted', marker='2')
ax.plot(sizes, [1/i.best for i in perf2_medfilt_opencl.values()], label="Medfilt/OpenCL", color='blue', linestyle='dotted', marker='2')
ax.set_yscale("log")
ax.legend()
ax.set_title("Performance of Median filtering vs integrate")
pass

As one can see, the penalities are much larger for OpenCL and Python than for Cython.

Conclusion

Sigma-clipping and median-filtering are alternatives to azimuthal integration and offer the ability to reject outlier. They are mot more difficult to use but slightly slower owing to their greater complexity.