Merge pull request #2691 from kif/2634_2pi_gi

Ensure 2pi periodicity is not used in 2D integration with fiber integration

Author: kif

src/pyFAI/engines/histogram_engine.py

@@ -26,10 +26,11 @@
 __contact__ = "Jerome.Kieffer@ESRF.eu"
 __license__ = "MIT"
 __copyright__ = "European Synchrotron Radiation Facility, Grenoble, France"
-__date__ = "06/10/2025"
+__date__ = "18/11/2025"
 __status__ = "development"
 
 import logging
+from math import pi
 import numpy
 from ..utils.mathutil import EPS32
 from .preproc import preproc as preproc_np
@@ -42,7 +43,7 @@
     preproc = preproc_np
 else:
     preproc = preproc_cy
-
+twopi = 2.0*pi
 
 
 def histogram1d_engine(radial, npt,
@@ -175,7 +176,7 @@ def histogram2d_engine(radial, azimuthal, bins,
                        azimuth_range=None,
                        allow_radial_neg=False,
                        chiDiscAtPi=True,
-                       clip_pos1=True
+                       pos1_period=twopi
                        ):
     """Implementation of 2D rebinning engine using pure numpy histograms
 

src/pyFAI/ext/histogram.pyx

@@ -38,7 +38,7 @@ Can be replaced by ``silx.math.histogramnd``.
 """
 
 __author__ = "Jérôme Kieffer"
-__date__ = "25/04/2024"
+__date__ = "18/11/2025"
 __license__ = "MIT"
 __copyright__ = "2011-2022, ESRF"
 __contact__ = "jerome.kieffer@esrf.fr"
@@ -583,7 +583,7 @@ def histogram2d_engine(radial, azimuthal,
                        azimuth_range=None,
                        bint allow_radial_neg=False,
                        bint chiDiscAtPi=1,
-                       bint clip_pos1=True
+                       position_t pos1_period=twopi
                        ):
     """Implementation of 2D rebinning engine using pure numpy histograms
 
@@ -609,7 +609,8 @@ def histogram2d_engine(radial, azimuthal,
     :param azimuth_range: enforce boundaries in azimuthal dimension, 2tuple with lower and upper bound
     :param allow_radial_neg: clip negative radial position (can a dimension be negative ?)
     :param chiDiscAtPi: set the azimuthal discontinuity at π (True) or at 0/2π (False)
-    :param clip_pos1: clip the azimuthal range to [-π π] (or [0 2π]), set to False to deactivate behavior.
+    :param pos1_period: periodicity of dim1, 2π, or 0 to non-periodic dimension
+    If pos1_period: clip_pos1 is enforced, the azimuthal range is set to [-π π] (or [0 2π] depending on chiDiscAtPi)
 
     NaN are always considered as invalid values
 
@@ -659,6 +660,7 @@ def histogram2d_engine(radial, azimuthal,
         bint do_dark = False, do_flat = False, do_polarization = False, do_solidangle = False, do_absorption=False
         bint do_variance=error_model, do_dark_variance=False
         preproc_t value
+        bint clip_pos1 = True if pos1_period>0 else False
 
     if variance is not None:
         assert variance.size == size, "variance size matches"

src/pyFAI/ext/regrid_common.pxi

@@ -32,7 +32,7 @@ Some are defined in the associated header file .pxd
 
 __author__ = "Jérôme Kieffer"
 __contact__ = "Jerome.kieffer@esrf.fr"
-__date__ = "12/11/2024"
+__date__ = "18/11/2025"
 __status__ = "stable"
 __license__ = "MIT"
 
@@ -373,16 +373,20 @@ cdef inline floating area4(floating a0,
 def _sp_area4(floating a0, floating a1, floating b0, floating b1, floating c0, floating c1, floating d0, floating d1):
     return area4(a0, a1, b0, b1, c0, c1, d0, d1)
 
-cdef inline position_t _recenter_helper(position_t azim, bint chiDiscAtPi) noexcept nogil:
+cdef inline position_t _recenter_helper(position_t azim,
+                                        position_t period,
+                                        bint chiDiscAtPi=1) noexcept nogil:
     """Helper function
     """
-    if (chiDiscAtPi and azim<0) or (not chiDiscAtPi and azim<pi):
-        return azim + twopi
+    if (chiDiscAtPi and azim<0) or (not chiDiscAtPi and azim<0.5*period):
+        return azim + period
     else:
         return azim
 
 
-cdef inline position_t _recenter(position_t[:, ::1] pixel, bint chiDiscAtPi) noexcept nogil:
+cdef inline position_t _recenter(position_t[:, ::1] pixel,
+                                 position_t pos1_period=twopi,
+                                 bint chiDiscAtPi=1) noexcept nogil:
     cdef position_t a0, a1, b0, b1, c0, c1, d0, d1, center1, area, hi
     a0 = pixel[0, 0]
     a1 = pixel[0, 1]
@@ -393,27 +397,29 @@ cdef inline position_t _recenter(position_t[:, ::1] pixel, bint chiDiscAtPi) noe
     d0 = pixel[3, 0]
     d1 = pixel[3, 1]
     area = area4p(a0, a1, b0, b1, c0, c1, d0, d1) # check if the quad is crossed
-    if area>0:
+    if pos1_period>0.0 and area>0:
         # area are expected to be negative except for pixel on the boundary
-        a1 = _recenter_helper(a1, chiDiscAtPi)
-        b1 = _recenter_helper(b1, chiDiscAtPi)
-        c1 = _recenter_helper(c1, chiDiscAtPi)
-        d1 = _recenter_helper(d1, chiDiscAtPi)
+        a1 = _recenter_helper(a1, pos1_period, chiDiscAtPi)
+        b1 = _recenter_helper(b1, pos1_period, chiDiscAtPi)
+        c1 = _recenter_helper(c1, pos1_period, chiDiscAtPi)
+        d1 = _recenter_helper(d1, pos1_period, chiDiscAtPi)
         center1 = 0.25 * (a1 + b1 + c1 + d1)
-        hi = pi if chiDiscAtPi else twopi
+        hi = 0.5*pos1_period if chiDiscAtPi else pos1_period
         if center1>hi:
-            a1 -= twopi
-            b1 -= twopi
-            c1 -= twopi
-            d1 -= twopi
+            a1 -= pos1_period
+            b1 -= pos1_period
+            c1 -= pos1_period
+            d1 -= pos1_period
         pixel[0, 1] = a1
         pixel[1, 1] = b1
         pixel[2, 1] = c1
         pixel[3, 1] = d1
         area = area4p(a0, a1, b0, b1, c0, c1, d0, d1)
     return area
 
-def recenter(position_t[:, ::1] pixel, bint chiDiscAtPi=1):
+def recenter(position_t[:, ::1] pixel,
+             position_t pos1_period=twopi,
+             bint chiDiscAtPi=True):
     """This function checks the pixel to be on the azimuthal discontinuity
     via the sign of its algebraic area and recenters the corner coordinates in a
     consistent manner to have all azimuthal coordinate in
@@ -424,7 +430,7 @@ def recenter(position_t[:, ::1] pixel, bint chiDiscAtPi=1):
     :param chiDiscAtPi: set to 0 to indicate the range goes from 0-2π instead of the default -π:π
     :return: signed area (approximate & negative)
     """
-    return _recenter(pixel, chiDiscAtPi)
+    return _recenter(pixel, pos1_period, chiDiscAtPi)
 
 
 cdef inline any_t _clip(any_t value, any_t min_val, any_t max_val) noexcept nogil:

src/pyFAI/ext/splitBBox.pyx

@@ -36,7 +36,7 @@ Splitting is done on the pixel's bounding box similar to fit2D
 
 __author__ = "Jérôme Kieffer"
 __contact__ = "Jerome.kieffer@esrf.fr"
-__date__ = "15/06/2024"
+__date__ = "18/11/2025"
 __status__ = "stable"
 __license__ = "MIT"
 
@@ -124,6 +124,7 @@ def histoBBox1d(weights,
         acc_t inv_area, delta_right, delta_left
         acc_t[::1] sum_data, sum_count
 
+
     cdata = numpy.ascontiguousarray(weights.ravel(), dtype=data_d)
     cpos0 = numpy.ascontiguousarray(pos0.ravel(), dtype=position_d)
     dpos0 = numpy.ascontiguousarray(delta_pos0.ravel(), dtype=position_d)
@@ -526,11 +527,11 @@ def histoBBox2d(weights,
                 solidangle=None,
                 polarization=None,
                 bint allow_pos0_neg=0,
+                position_t pos1_period=twopi,
                 bint chiDiscAtPi=1,
                 empty=0.0,
                 double normalization_factor=1.0,
                 int coef_power=1,
-                bint clip_pos1=1,
                 **back_compat_kwargs):
     """
     Calculate 2D histogram of pos0(tth),pos1(chi) weighted by weights
@@ -557,8 +558,8 @@ def histoBBox2d(weights,
     :param empty: value of output bins without any contribution when dummy is None
     :param normalization_factor: divide the result by this value
     :param coef_power: set to 2 for variance propagation, leave to 1 for mean calculation
-    :param clip_pos1: clip the azimuthal range to -pi/pi (or 0-2pi), set to False to deactivate behavior
-
+    :param pos1_period: periodicity of dim1, 2π, or 0 to non-periodic dimension
+        If pos1_period: clip_pos1 is enforced, the azimuthal range is set to [-π π] (or [0 2π] depending on chiDiscAtPi)
 
     :return: I, bin_centers0, bin_centers1, weighted histogram(2D), unweighted histogram (2D)
     """
@@ -611,6 +612,7 @@ def histoBBox2d(weights,
         Py_ssize_t  bin0_max, bin0_min, bin1_max, bin1_min
         bint check_mask = False, check_dummy = False
         bint do_dark = False, do_flat = False, do_polarization = False, do_solidangle = False
+        bint clip_pos1 = True if pos1_period>0 else False
 
     if mask is not None:
         assert mask.size == size, "mask size"
@@ -810,12 +812,11 @@ def histoBBox2d_engine(weights,
                        polarization=None,
                        absorption=None,
                        bint allow_pos0_neg=False,
+                       position_t pos1_period=twopi,
                        bint chiDiscAtPi=1,
                        data_t empty=0.0,
                        double normalization_factor=1.0,
-                       bint weighted_average=True,
-                       bint clip_pos1=True
-                       ):
+                       bint weighted_average=True):
     """
     Calculate 2D histogram of pos0(tth),pos1(chi) weighted by weights
 
@@ -845,7 +846,8 @@ def histoBBox2d_engine(weights,
     :param empty: value of output bins without any contribution when dummy is None
     :param normalization_factor: divide the result by this value
     :param bool weighted_average: set to False to use an unweighted mean (similar to legacy) instead of the weighted average.
-    :param clip_pos1: clip the azimuthal range to [-pi pi] (or [0 2pi]), set to False to deactivate behavior
+    :param pos1_period: periodicity of dim1, 2π, or 0 to non-periodic dimension
+        If pos1_period: clip_pos1 is enforced, the azimuthal range is set to [-π π] (or [0 2π] depending on chiDiscAtPi)
     :return: Integrate2dtpl namedtuple: "radial azimuthal intensity error signal variance normalization count"
     """
 
@@ -887,6 +889,7 @@ def histoBBox2d_engine(weights,
         bint do_dark = False, do_flat = False, do_polarization = False, do_solidangle = False
         bint do_absorption = False
         preproc_t value
+        bint clip_pos1 = True if pos1_period>0 else False
 
     if variance is not None:
         assert variance.size == size, "variance size"

src/pyFAI/ext/splitBBoxCSC.pyx

@@ -37,7 +37,7 @@ Serial implementation based on a sparse CSC matrix multiplication
 
 __author__ = "Jérôme Kieffer"
 __contact__ = "Jerome.kieffer@esrf.fr"
-__date__ = "26/01/2024"
+__date__ = "18/11/2025"
 __status__ = "stable"
 __license__ = "MIT"
 
@@ -79,7 +79,7 @@ class HistoBBox1d(CscIntegrator, SplitBBoxIntegrator):
                  unit="undefined",
                  empty=None,
                  bint chiDiscAtPi=True,
-                 bint clip_pos1=False):
+                 position_t pos1_period=twopi):
         """
         :param pos0: 1D array with pos0: tth or q_vect or r ...
         :param delta_pos0: 1D array with delta pos0: max center-corner distance
@@ -93,14 +93,16 @@ class HistoBBox1d(CscIntegrator, SplitBBoxIntegrator):
         :param unit: can be 2th_deg or r_nm^-1 ...
         :param empty: value for bins without contributing pixels
         :param chiDiscAtPi: tell if azimuthal discontinuity is at 0° or 180°
-        :param clip_pos1: clip the azimuthal range to [-π π] (or [0 2π] depending on chiDiscAtPi), set to False to deactivate behavior
+        :param pos1_period: periodicity of dim1, 2π, or 0 to non-periodic dimension
+            If pos1_period: clip_pos1 is enforced, the azimuthal range is set to [-π π] (or [0 2π] depending on chiDiscAtPi)
         """
         self.unit = unit
         self.space = tuple(str(u).split("_")[0] for u in unit) if isinstance(unit, (tuple, list)) else  str(unit).split("_")[0]
         SplitBBoxIntegrator.__init__(self, pos0, delta_pos0, pos1, delta_pos1,
                                      bins, pos0_range, pos1_range,
                                      mask, mask_checksum,
-                                     allow_pos0_neg, chiDiscAtPi, clip_pos1=clip_pos1)
+                                     allow_pos0_neg, chiDiscAtPi,
+                                     pos1_period=pos1_period)
 
 
         self.delta = (self.pos0_max - self.pos0_min) / (<position_t> (self.bins))
@@ -155,7 +157,7 @@ class HistoBBox2d(CscIntegrator, SplitBBoxIntegrator):
                  unit="undefined",
                  empty=None,
                  bint chiDiscAtPi=True,
-                 bint clip_pos1=True):
+                 position_t pos1_period=twopi):
         """
         :param pos0: 1D array with pos0: tth or q_vect
         :param delta_pos0: 1D array with delta pos0: max center-corner distance
@@ -169,11 +171,12 @@ class HistoBBox2d(CscIntegrator, SplitBBoxIntegrator):
         :param unit: can be 2th_deg or r_nm^-1 ...
         :param empty: value for bins where no pixels are contributing
         :param chiDiscAtPi: tell if azimuthal discontinuity is at 0 (0° when False) or π (180° when True)
-        :param clip_pos1: clip the azimuthal range to [-π π] (or [0 2π] depending on chiDiscAtPi), set to False to deactivate behavior
+        :param pos1_period: periodicity of dim1, 2π, or 0 to non-periodic dimension
+            If pos1_period: clip_pos1 is enforced, the azimuthal range is set to [-π π] (or [0 2π] depending on chiDiscAtPi)
         """
         SplitBBoxIntegrator.__init__(self, pos0, delta_pos0, pos1, delta_pos1,
                                      bins, pos0_range, pos1_range, mask, mask_checksum, allow_pos0_neg, chiDiscAtPi,
-                                     clip_pos1)
+                                     pos1_period=pos1_period)
         self.unit = unit
         self.space = tuple(str(u).split("_")[0] for u in unit) if isinstance(unit, (tuple, list)) else  str(unit).split("_")[0]
         self.bin_centers = None

src/pyFAI/ext/splitBBoxCSR.pyx

@@ -37,7 +37,7 @@ reverse implementation based on a sparse matrix multiplication
 
 __author__ = "Jérôme Kieffer"
 __contact__ = "Jerome.kieffer@esrf.fr"
-__date__ = "04/10/2023"
+__date__ = "18/11/2025"
 __status__ = "stable"
 __license__ = "MIT"
 
@@ -79,7 +79,7 @@ class HistoBBox1d(CsrIntegrator, SplitBBoxIntegrator):
                  unit="undefined",
                  empty=None,
                  bint chiDiscAtPi=True,
-                 bint clip_pos1=False):
+                 position_t pos1_period=twopi):
         """
         :param pos0: 1D array with pos0: tth or q_vect or r ...
         :param delta_pos0: 1D array with delta pos0: max center-corner distance
@@ -93,14 +93,15 @@ class HistoBBox1d(CsrIntegrator, SplitBBoxIntegrator):
         :param unit: can be 2th_deg or r_nm^-1 ...
         :param empty: value for bins without contributing pixels
         :param chiDiscAtPi: tell if azimuthal discontinuity is at 0° or 180°
-        :param clip_pos1: clip the azimuthal range to [-π π] (or [0 2π] depending on chiDiscAtPi), set to False to deactivate behavior
+        :param pos1_period: periodicity of dim1, 2π, or 0 to non-periodic dimension
+            If pos1_period: clip_pos1 is enforced, the azimuthal range is set to [-π π] (or [0 2π] depending on chiDiscAtPi)
         """
         self.unit = unit
         self.space = tuple(str(u).split("_")[0] for u in unit) if isinstance(unit, (list, tuple)) else  str(unit).split("_")[0]
         SplitBBoxIntegrator.__init__(self, pos0, delta_pos0, pos1, delta_pos1,
                                      bins, pos0_range, pos1_range,
                                      mask, mask_checksum,
-                                     allow_pos0_neg, chiDiscAtPi, clip_pos1=clip_pos1)
+                                     allow_pos0_neg, chiDiscAtPi, pos1_period=pos1_period)
 
 
         self.delta = (self.pos0_max - self.pos0_min) / (<position_t> (self.bins))
@@ -153,7 +154,7 @@ class HistoBBox2d(CsrIntegrator, SplitBBoxIntegrator):
                  unit="undefined",
                  empty=None,
                  bint chiDiscAtPi=True,
-                 bint clip_pos1=True):
+                 position_t pos1_period=twopi):
         """
         :param pos0: 1D array with pos0: tth or q_vect
         :param delta_pos0: 1D array with delta pos0: max center-corner distance
@@ -167,11 +168,12 @@ class HistoBBox2d(CsrIntegrator, SplitBBoxIntegrator):
         :param unit: can be 2th_deg or r_nm^-1 ...
         :param empty: value for bins where no pixels are contributing
         :param chiDiscAtPi: tell if azimuthal discontinuity is at 0 (0° when False) or π (180° when True)
-        :param clip_pos1: clip the azimuthal range to [-π π] (or [0 2π] depending on chiDiscAtPi), set to False to deactivate behavior
+        :param pos1_period: periodicity of dim1, 2π, or 0 to non-periodic dimension
+            If pos1_period: clip_pos1 is enforced, the azimuthal range is set to [-π π] (or [0 2π] depending on chiDiscAtPi)
         """
         SplitBBoxIntegrator.__init__(self, pos0, delta_pos0, pos1, delta_pos1,
                                      bins, pos0_range, pos1_range, mask, mask_checksum, allow_pos0_neg, chiDiscAtPi,
-                                     clip_pos1)
+                                     pos1_period=pos1_period)
         self.unit = unit
         self.space = tuple(str(u).split("_")[0] for u in unit) if isinstance(unit, tuple) else  str(unit).split("_")[0]
         self.bin_centers = None