Source-linked AI summary
PyFAI: a Python library for high performance azimuthal integration on GPU
Jérôme Kieffer, Giannis Ashiotis
TL;DR
Fast pixel detectors created a need for Python azimuthal integration that existing programs could not meet. pyFAI implements and accelerates the algorithm across CPU and many-core approaches; its conclusion reports 15x speed-up over binary code and 400x over NumPy for OpenCL on GPUs.
Problem
High-speed pixel detectors impose demands that available Fortran or C programs could not meet, motivating a convenient Python azimuthal-integration routine.
Source: §1 · p. 1Method
Results
15x speed-up was obtained by switching from binary code to OpenCL on GPUs, and the conclusion reports 400x versus NumPy.
Source: §8 · p. 7Takeaways & Limitations
The authors report that Python-interfaced parallel programming can match fast-detector readout speeds, with strong GPU performance including a mid-range GeForce 750Ti.
Source: §1 · p. 1, §8 · pp. 7–8Cython OpenMP parallelization produced wrong results, with 2% errors on eight cores, because atomic operations were unavailable.
Source: §5.2 · p. 4
Abstract
from arXiv · showhide
The pyFAI package has been designed to reduce X-ray diffraction images into powder diffraction curves to be further processed by scientists. This contribution describes how to convert an image into a radial profile using the Numpy package, how the process was accelerated using Cython. The algorithm was parallelised, needing a complete re-design to benefit from massively parallel devices like graphical processing units or accelerators like the Intel Xeon Phi using the PyOpenCL library.
1 INTRODUCTION
The introduction motivates a convenient Python azimuthal-integration routine for a community facing detector-speed demands. It presents pyFAI as a Python interface whose implementation is accelerated for fast detectors.
1 INTRODUCTION
- High-speed pixel detectors, reaching 3000 frames per second, imposed demands that available Fortran or C programs could not meet.The synchrotron community requested a convenient azimuthal-integration routine for Python, widely adopted in crystallography.Source: §1 · p. 1
- pyFAI aims to provide an adaptable Pythonic interface for azimuthal integration and accelerate the algorithm to match fast-detector readout speeds.The project started in 2011.Source: §1 · p. 1
- The contribution covers NumPy vectorisation, Cython acceleration, accuracy enhancement, PyOpenCL parallelisation, and OpenMP/OpenCL benchmarks.Source: §1 · p. 1
2 DESCRIPTION OF THE EXPERIMENT
The section explains how X-ray diffraction from crystalline powders produces rings whose intensity is measured as a function of aperture. Azimuthal integration averages the signal along each ring.
2 DESCRIPTION OF THE EXPERIMENT
- X-rays have wavelengths near inter-atomic distances, making them a probe of atomic and molecular structures.Source: §2 · p. 1
- Constructive interference occurs when crystal-scattered X-ray paths differ by a multiple of the wavelength, as expressed by Bragg’s law.The passage gives the relation 2dsin(θ) = nλ and defines d, θ, λ, and n.Source: §2 · p. 1
- Randomly oriented crystallites scatter X-rays along concentric cones, and powder diffraction measures intensity by cone aperture through azimuthal averaging.Source: §2 · p. 1
3 AZIMUTHAL INTEGRATION
The section develops azimuthal integration from a Python-facing library interface and simulated test case through NumPy histograms and Cython implementations. The naive repeated-slicing method is slow, while histogram-based methods reduce execution time; OpenMP gains are limited and require a redesign for greater speed.
3 AZIMUTHAL INTEGRATION
- pyFAI provides a Pythonic interface for image integration, with selectable output space and units, and this contribution limits its experiment to 1D full-azimuthal integration on a planar detector.The example uses FabIO for image reading and a pyFAI geometry file.Source: §3 · pp. 1–2
3.1 Test case
- The test case uses a simulated 2048×2048 gold-powder diffraction image with equal-intensity rings and a planar detector orthogonal to the beam.Its radial coordinate is radius, although crystallographers would prefer 2θ or scattering-vector length q.Source: §3.1 · p. 2
3.2 Naive implementation
- 40 seconds per image made the naive method unacceptably slow because slicing took tens of milliseconds and was repeated thousands of times.Source: §3.2 · p. 2
3.3 Numpy histograms
- 800 ms is the reported execution time for NumPy histogram integration, faster than the Python loop and open to further optimization through a single radius-array read.The implementation divides the weighted-value histogram by the unweighted radius histogram.Source: §3.3 · pp. 2–3
3.4 Cython implementation
- 150 ms on one core is the reported Cython histogram time after generating weighted and unweighted histograms with one radius-array read.The authors attribute the improvement to better CPU-cache use.Source: §3.4 · p. 3
3.4.1 OpenMP support in Cython
- OpenMP histogramming gains were minor, especially beyond two threads, and write conflicts without atomic_add operations produced errors of a few percent.The authors conclude that faster execution requires redesigning the algorithm for many threads from the beginning.Source: §3.4.1 · p. 3
4 PIXEL SPLITTING
Pixel splitting distributes a detector pixel’s contribution across the histogram bins it spans, aiming to reduce noise when bin statistics are low. The paper describes both a bounding-box approximation and a more precise full-pixel scheme, while noting effects on signal blur and bin statistics.
4 PIXEL SPLITTING
- With pixel splitting, high-frequency patterns near the beam center in the 2D integration disappear, addressing defects associated with low statistics.The paper cautions that the transformation resembles interpolation but does not guarantee conservation of the signal or pixel count.Source: §4 · p. 3, Figure 3 · p. 3
4.1 Bounding Box
- The bounding-box method approximates each pixel by a rectangle and distributes contributions among bins using the bounding-box segment areas.A pixel wholly within one bin is processed without splitting; otherwise, the outer-bin contributions are calculated first and the remainder is distributed across internal bins.Source: §4.1 · p. 4
- The bounding-box approximation overestimates pixel size, causing slight signal blurring.Source: §4.1 · p. 4
4.2 Full Pixel Splitting
- Full pixel splitting integrates the functions describing the pixel’s sides over each bin’s bounds to calculate segment areas, adding computational complexity and slowing processing.The method assumes straight pixel edges and requires careful integral limits to avoid excessive conditionals.Source: §4.2 · p. 4
4.3 Discussion on the statistics
- Either splitting scheme creates cross-correlation between neighboring histogram bins, potentially affecting the statistics of the results.Source: §4.3 · p. 4
5 MORE PARALLELISATION
The paper redesigns azimuthal integration for parallel execution, addressing atomic-operation constraints by changing the computation from scatter to gather and using sparse-matrix methods. It reports faster execution with CSR and up to 750 MPixel/s on recent computers with a mid-range graphics card.
5.1 Parallel algorithms
- Parallel algorithms are assembled from building blocks such as map, scatter, gather, reduction, scan, and sort.OpenCL execution on processors, graphics cards, or accelerators depends on memory access patterns, while many-core devices run numerous threads simultaneously.Source: §5.1 · p. 4, §5 · p. 4
5.2 Parallel azimuthal integration
- Azimuthal integration is a scatter operation requiring double-precision atomic operations; the Cython OpenMP implementation produced wrong results, with 2% errors on 8 cores.Source: §5.2 · p. 4
- The scatter-to-gather transformation stores pixel-to-bin correspondences and weights in a lookup table, enabling a backward algorithm with random reads and linear writes.The lookup-table array may be transposed to suit CPU or GPU memory access patterns.Source: §5.2 · pp. 4–5
5.2.1 Optimisation of the sparse matrix multiplication
- CSR reduces LUT storage needs by a factor of two to three and enables larger images on the same hardware; its OpenCL implementation uses parallel reductions, with workgroup threads collaborating to calculate each bin.Source: §5.2.1 · p. 5
5.2.2 About numerical precision
- Kahan summation enables numerically accurate results with single-precision floating-point arithmetic on consumer-grade GPUs, at a cost of four times as many arithmetic operations.The paper says that cost is hidden by smaller data types, more single-precision units, and typical GPU memory-bandwidth limits.Source: §5.2.2 · p. 5
- 750 MPixel/s is reported for parallel azimuthal integration on recent computers with a mid-range graphics card.On multi-socket servers with high-end GPUs such as Tesla cards, performance is equivalent, with the added benefit of processing multiple detectors simultaneously.Source: §5.2.2 · p. 5
6 BENCHMARKS
Benchmarks compare integration algorithms and OpenMP/OpenCL implementations across CPUs and many-core accelerators, then profile the OpenCL pipeline. The reported comparisons cover speed, memory use, and execution-time components under the stated benchmark conditions.
6 BENCHMARKS
- Benchmarks use equivalent-quality profiles, bounding-box pixel splitting, and images from actual diffraction experiments; the 12 Mpixel cast ran on one CPU thread.Times were averaged over 10 iterations, taking the best of 3; tested image sizes ranged from 1 to 16 Mpixel.Source: §6 · p. 5
6.1 Choice of the algorithm
- On two Xeon E2667s, both LUT and CSR outperform serial code and behave similarly because CSR's smaller transfer volume offsets its array-of-struct penalty.LUT's array-of-struct layout favors CPU caches; CSR transfers a smaller data chunk.Source: §6.1 · p. 5, Figure 4 · p. 6
6.2 OpenMP vs OpenCL
- OpenCL outperforms OpenMP on every tested CPU, while the dual Xeon E5520's OpenMP performance compares well with newer systems.The E5520 was the only tested computer with hyper-threading enabled.Source: Figure 5 · p. 6, §6.2 · p. 6
- On Ivy Bridge, the Intel OpenCL driver clearly outperforms AMD's, whereas the driver-speed difference is much smaller on Nehalem.The authors attribute the Ivy Bridge gap to AVX instructions supported by Intel's driver but not AMD's.Source: Figure 6 · p. 6, §6.2 · p. 6
6.3 GPUs and Xeon Phi
- CSR runs much faster than LUT on the tested GPUs and processes 12 Mpixel images where LUT is limited to 4 Mpixel on the 1 GB ATI GPU.The authors attribute CSR's speed to collaborative parallel reductions and its larger processable image size to lower memory use.Source: Figure 7 · p. 6, §6.3 · p. 6
- Xeon Phi matches the 2010 AMD GPU, while the consumer GeForce 750Ti performs well against higher-priced Kepler cards.The comparison includes the Titan and Tesla K20.Source: §6.3 · p. 6, Figure 8 · p. 7
6.4 Kernel timings
- On the GeForce Titan, Python-level execution takes 2 ms versus 1.4 ms summed OpenCL-profiler time; image transfer takes 0.8 ms and integration 0.4 ms.The authors report that Python overhead is under 40% and integration is 20% of total runtime.Source: Table 3 · p. 7, §6.4 · p. 7
6.5 Configuration and Drivers used
- The benchmark configurations span three Debian systems with Xeon CPUs and a range of Nvidia and AMD GPUs, plus an Intel Xeon Phi accelerator.The systems use Debian 7 or Debian 8/Jessie and vendor-specific OpenCL drivers.Source: §6.5 · p. 7
7 PROJECT DESCRIPTION
PyFAI is open-source software for azimuthal integration, with detector-reading, numerical, acceleration, and graphical components and support across several operating systems.
7 PROJECT DESCRIPTION
- PyFAI is open-source software released under the GPL and available on GitHub.Source: §7 · p. 7
- The package depends on Python 2.6 or 2.7 and NumPy, and uses FabIO to read images from different X-ray detectors.Source: §7 · p. 7
- Optional OpenCL acceleration uses PyOpenCL; graphical applications rely on matplotlib, PyQt, and SciPy, with a C compiler needed to build Cython code.Source: §7 · p. 7
8 CONCLUSIONS
The contribution implements azimuthal integration in Python, optimizes it with Cython, and ports it to many-core architectures using PyOpenCL. OpenCL on GPUs achieves substantial speedups, while a consumer GeForce 750Ti performs strongly against high-end and accelerator hardware.
8 CONCLUSIONS
- 15x speed-up over binary code and 400x over NumPy were obtained with OpenCL code running on GPUs.Source: §8 · p. 7
- The consumer GeForce 750Ti runs as fast as high-end Kepler graphics such as Titan and outperforms AMD GPUs and Xeon Phi.Source: §8 · p. 8
- Programming CPUs in parallel is described as being as easy as programming GPUs using PyOpenCL interfaced with Python.Source: §8 · p. 8