Filtering signal in azimuthal space

Usually, diffraction signal presents a polar symmetry, this means all pixels 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 outliers, like gaps, missing pixels, shadows, cosmic rays or reflections coming from larger crystallites. In this tutorial we will see two alternative ways to remove those unwanted signals and focus on the majority of pixels: sigma clipping and median filtering.

import os
os.environ["PYOPENCL_COMPILER_OUTPUT"]="0"
import pyFAI
print(f"pyFAI version: {pyFAI.version}")

pyFAI version: 2026.8.0-dev0

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

figsize = (10,5)

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

WARNING:pyFAI.gui.matplotlib:Matplotlib already loaded with backend `inline`, setting its backend to `QtAgg` 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, figsize=figsize)
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");

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

• The algorithm needs to be the CSR-sparse matrix multiplication: since several integrations are needed, it makes no sense to use a 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
     * integrate 

17.9 ms ± 3.39 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
     * sigma-clip

15.6 ms ± 585 μs per loop (mean ± std. dev. of 7 runs, 100 loops each)
     Python
     * integrate 

10.3 ms ± 52.1 μs per loop (mean ± std. dev. of 7 runs, 100 loops each)
     * sigma-clip

158 ms ± 1.43 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
     OpenCL
     * integrate 

685 μs ± 2.07 μs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)
     * sigma-clip

2.45 ms ± 7.75 μs per loop (mean ± std. dev. of 7 runs, 100 loops each)

Pilatus2M.poni
     Cython
     * integrate 

23.1 ms ± 2.39 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
     * sigma-clip

23.5 ms ± 415 μs per loop (mean ± std. dev. of 7 runs, 10 loops each)
     Python
     * integrate 

34.1 ms ± 474 μs per loop (mean ± std. dev. of 7 runs, 10 loops each)
     * sigma-clip

555 ms ± 11.8 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
     OpenCL
     * integrate 

1.1 ms ± 17.1 μs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)
     * sigma-clip

6.04 ms ± 31.5 μs per loop (mean ± std. dev. of 7 runs, 100 loops each)

Eiger4M.poni
     Cython
     * integrate 

33.4 ms ± 4.73 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
     * sigma-clip

32.4 ms ± 792 μs per loop (mean ± std. dev. of 7 runs, 10 loops each)
     Python
     * integrate 

60.8 ms ± 577 μs per loop (mean ± std. dev. of 7 runs, 10 loops each)
     * sigma-clip

1.15 s ± 9.65 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
     OpenCL
     * integrate 

2.22 ms ± 22.8 μs per loop (mean ± std. dev. of 7 runs, 100 loops each)
     * sigma-clip

10.9 ms ± 29 μs per loop (mean ± std. dev. of 7 runs, 100 loops each)

Pilatus6M.poni
     Cython
     * integrate 

45.5 ms ± 3.22 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
     * sigma-clip

39.6 ms ± 322 μs per loop (mean ± std. dev. of 7 runs, 10 loops each)
     Python
     * integrate 

83.7 ms ± 123 μs per loop (mean ± std. dev. of 7 runs, 10 loops each)
     * sigma-clip

1.59 s ± 6.49 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
     OpenCL
     * integrate 

3.06 ms ± 13.8 μs per loop (mean ± std. dev. of 7 runs, 100 loops each)
     * sigma-clip

14.1 ms ± 16.4 μs per loop (mean ± std. dev. of 7 runs, 100 loops each)

Eiger9M.poni
     Cython
     * integrate 

63.4 ms ± 1.45 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
     * sigma-clip

65.4 ms ± 677 μs per loop (mean ± std. dev. of 7 runs, 10 loops each)
     Python
     * integrate 

141 ms ± 702 μs per loop (mean ± std. dev. of 7 runs, 10 loops each)
     * sigma-clip

3.09 s ± 4.27 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
     OpenCL
     * integrate 

4.68 ms ± 34.2 μs per loop (mean ± std. dev. of 7 runs, 100 loops each)
     * sigma-clip

25.8 ms ± 34.3 μs per loop (mean ± std. dev. of 7 runs, 10 loops each)

Mar3450.poni
     Cython
     * integrate 

67.3 ms ± 2.43 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
     * sigma-clip

72.3 ms ± 464 μs per loop (mean ± std. dev. of 7 runs, 10 loops each)
     Python
     * integrate 

166 ms ± 205 μs per loop (mean ± std. dev. of 7 runs, 10 loops each)
     * sigma-clip

3.42 s ± 3.16 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
     OpenCL
     * integrate 

5.26 ms ± 11.7 μs per loop (mean ± std. dev. of 7 runs, 100 loops each)
     * sigma-clip

28.8 ms ± 23.5 μs per loop (mean ± std. dev. of 7 runs, 10 loops each)

Fairchild.poni
     Cython
     * integrate 

115 ms ± 18.4 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
     * sigma-clip

189 ms ± 52.3 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
     Python
     * integrate 

381 ms ± 1.31 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
     * sigma-clip

9.54 s ± 12 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
     OpenCL
     * integrate 

4.91 ms ± 43.1 μs per loop (mean ± std. dev. of 7 runs, 100 loops each)
     * sigma-clip

33.6 ms ± 82.6 μs per loop (mean ± std. dev. of 7 runs, 10 loops each)
CPU times: user 32min 22s, sys: 40.5 s, total: 33min 2s
Wall time: 5min 40s

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");

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

The biggest limitation of sigma-clipping is its incompatibility with pixel-splitting, a 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 the general case (due to the cross-correlation between bins it creates), it can be a necessary evil, especially when performing Rietveld refinement where 5 points per peak 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 contributing to an azimuthal bin and to average out all pixels 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 cumulative sum of normalization is performed in order to determine 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, figsize=figsize)
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");

Unlike the sigma-clipping, this median filter does not require any error model; but the computational 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
     * integrate   

21.4 ms ± 3.08 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
     * medianfilter

21.4 ms ± 1.62 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
     Python
     * integrate   

12.8 ms ± 45.4 μs per loop (mean ± std. dev. of 7 runs, 100 loops each)
     * medianfilter

1.16 s ± 7.68 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
     OpenCL
     * integrate   

731 μs ± 1.94 μs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)
     * medianfilter

9.45 ms ± 8.14 μs per loop (mean ± std. dev. of 7 runs, 100 loops each)
Pilatus2M.poni
     Cython
     * integrate   

27.4 ms ± 2.42 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
     * medianfilter

68.4 ms ± 842 μs per loop (mean ± std. dev. of 7 runs, 10 loops each)
     Python
     * integrate   

46.2 ms ± 264 μs per loop (mean ± std. dev. of 7 runs, 10 loops each)
     * medianfilter

3.96 s ± 13.6 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
     OpenCL
     * integrate   

1.23 ms ± 6.48 μs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)
     * medianfilter

30.9 ms ± 27 μs per loop (mean ± std. dev. of 7 runs, 10 loops each)

Eiger4M.poni
     Cython
     * integrate   

33.8 ms ± 1.52 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
     * medianfilter

114 ms ± 728 μs per loop (mean ± std. dev. of 7 runs, 10 loops each)
     Python
     * integrate   

81.4 ms ± 169 μs per loop (mean ± std. dev. of 7 runs, 10 loops each)
     * medianfilter

6.85 s ± 27.7 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
     OpenCL
     * integrate   

2.36 ms ± 10.8 μs per loop (mean ± std. dev. of 7 runs, 100 loops each)
     * medianfilter

54.5 ms ± 25.1 μs per loop (mean ± std. dev. of 7 runs, 10 loops each)

Pilatus6M.poni
     Cython
     * integrate   

46.3 ms ± 4.75 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
     * medianfilter

153 ms ± 1.5 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
     Python
     * integrate   

113 ms ± 640 μs per loop (mean ± std. dev. of 7 runs, 10 loops each)
     * medianfilter

10.1 s ± 26.7 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
     OpenCL
     * integrate   

3.29 ms ± 14.8 μs per loop (mean ± std. dev. of 7 runs, 100 loops each)
     * medianfilter

78.2 ms ± 72.4 μs per loop (mean ± std. dev. of 7 runs, 10 loops each)

Eiger9M.poni
     Cython
     * integrate   

65.6 ms ± 2.35 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
     * medianfilter

247 ms ± 4.66 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
     Python
     * integrate   

180 ms ± 3.12 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
     * medianfilter

16.8 s ± 37.1 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
     OpenCL
     * integrate   

5.11 ms ± 20.7 μs per loop (mean ± std. dev. of 7 runs, 1 loop each)
     * medianfilter

139 ms ± 214 μs per loop (mean ± std. dev. of 7 runs, 10 loops each)

Mar3450.poni
     Cython
     * integrate   

69.2 ms ± 3.06 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
     * medianfilter

259 ms ± 3.3 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
     Python
     * integrate   

206 ms ± 3.51 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
     * medianfilter

19.2 s ± 41.3 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
     OpenCL
     * integrate   

5.97 ms ± 19.7 μs per loop (mean ± std. dev. of 7 runs, 1 loop each)
     * medianfilter

163 ms ± 147 μs per loop (mean ± std. dev. of 7 runs, 10 loops each)

Fairchild.poni
     Cython
     * integrate   

106 ms ± 2 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
     * medianfilter

404 ms ± 4.45 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
     Python
     * integrate   

307 ms ± 1.07 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
     * medianfilter

23.9 s ± 103 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
     OpenCL
     * integrate   

5.35 ms ± 38.2 μs per loop (mean ± std. dev. of 7 runs, 1 loop each)
     * medianfilter

173 ms ± 265 μs per loop (mean ± std. dev. of 7 runs, 10 loops each)
CPU times: user 25min 7s, sys: 24.6 s, total: 25min 32s
Wall time: 13min 46s

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");

As one can see, the penalties 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 outliers. They are not more difficult to use but slightly slower owing to their greater complexity.