Merge pull request #2081 from kif/2078_optimize_extract_cp

Improve performances for ring extraction

doc/source/usage/tutorial/Goniometer/Rotation-Pilatus100k/ROBL_v1.json

@@ -1,15 +1,17 @@
 {
   "content": "Goniometer calibration v2",
   "detector": "Pilatus 100k",
-  "detector_config": {},
+  "detector_config": {
+    "orientation": 3
+  },
   "wavelength": 7.7490120575e-11,
   "param": [
-    0.7994908637172136,
-    -0.0248850431840963,
-    0.12722802891692944,
-    0.10448630616575466,
-    0.048840045382446656,
-    0.01745706471069648
+    0.7999015765795351,
+    -0.042536815083314417,
+    0.096018037587983,
+    0.06541469598422345,
+    0.0709179790175431,
+    0.0174568035572029
   ],
   "param_names": [
     "dist",

doc/source/usage/tutorial/Goniometer/Rotation-Pilatus100k/ROBL_v2.json

@@ -1,16 +1,18 @@
 {
   "content": "Goniometer calibration v2",
   "detector": "Pilatus 100k",
-  "detector_config": {},
+  "detector_config": {
+    "orientation": 3
+  },
   "wavelength": 7.7490120575e-11,
   "param": [
-    0.8028183874123671,
-    -0.030497146146421752,
-    0.06660620717040736,
-    0.029413083810715633,
-    -9.57985973026547e-05,
-    0.0558376290842621,
-    0.017457086319994403
+    0.8004534770437763,
+    -0.044802777930768284,
+    0.07173709328138,
+    0.03557894996410206,
+    -6.461701150618871e-05,
+    0.07378133590152441,
+    0.017456966023298034
   ],
   "param_names": [
     "dist",

doc/source/usage/tutorial/Goniometer/Translation-Pilatus6M/ID29.json

@@ -1,7 +1,9 @@
 {
   "content": "Goniometer calibration v2",
   "detector": "Pilatus 6M",
-  "detector_config": {},
+  "detector_config": {
+    "orientation": 3
+  },
   "wavelength": 9.72386e-11,
   "param": [
     -0.001187934615261645,

src/pyFAI/ext/invert_geometry.pyx

@@ -37,8 +37,8 @@ coordinate.
 
 __author__ = "Jerome Kieffer"
 __license__ = "MIT"
-__date__ = "03/03/2023"
-__copyright__ = "2018-2018, ESRF"
+__date__ = "20/02/2024"
+__copyright__ = "2018-2024, ESRF"
 __contact__ = "jerome.kieffer@esrf.fr"
 
 include "regrid_common.pxi"

src/pyFAI/ext/mathutil.pyx

@@ -0,0 +1,120 @@
+# -*- coding: utf-8 -*-
+#cython: embedsignature=True, language_level=3, binding=True
+#cython: boundscheck=False, wraparound=False, cdivision=True, initializedcheck=False,
+## This is for developping:
+##cython: profile=True, warn.undeclared=True, warn.unused=True, warn.unused_result=False, warn.unused_arg=True
+#
+#    Project: Fast Azimuthal integration
+#             https://github.com/silx-kit/pyFAI
+#
+#    Copyright (C) 2012-2020 European Synchrotron Radiation Facility, Grenoble, France
+#
+#    Principal author:       Jérôme Kieffer (Jerome.Kieffer@ESRF.eu)
+#
+#  Permission is hereby granted, free of charge, to any person obtaining a copy
+#  of this software and associated documentation files (the "Software"), to deal
+#  in the Software without restriction, including without limitation the rights
+#  to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+#  copies of the Software, and to permit persons to whom the Software is
+#  furnished to do so, subject to the following conditions:
+#  .
+#  The above copyright notice and this permission notice shall be included in
+#  all copies or substantial portions of the Software.
+#  .
+#  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+#  IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+#  FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+#  AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+#  LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+#  OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+#  THE SOFTWARE.
+
+
+"""
+Extenstion with Cython implementation of pyFAI.utils.mathutil
+"""
+
+__author__ = "Jerome Kieffer"
+__license__ = "MIT"
+__date__ = "20/02/2024"
+__copyright__ = "2024-2024, ESRF"
+__contact__ = "jerome.kieffer@esrf.fr"
+
+from libc.stdint cimport int8_t, uint8_t, int16_t, uint16_t, \
+                         int32_t, uint32_t, int64_t, uint64_t
+import numpy
+
+def is_far_from_group_cython(pt, list lst_pts, double d2):
+    """
+    Tells if a point is far from a group of points, distance greater than d2 (distance squared)
+
+    :param pt: point of interest
+    :param lst_pts: list of points
+    :param d2: minimum distance squarred
+    :return: True If the point is far from all others.
+
+    """
+    cdef:
+        double p0, p1, q0, q1, dsq, d0, d1
+    p0, p1 = pt
+    for (q0, q1) in lst_pts:
+        d0 = p0 - q0
+        d1 = p1 - q1
+        dsq = d0*d0 + d1*d1
+        if dsq <= d2:
+            return False
+    return True
+
+
+def build_qmask(tth_array,
+                tth_min,
+                tth_max,
+                mask=None):
+    """Build a qmask, array with the index of the reflection for each pixel.
+    qmask==-1 is for uninteresting pixels
+    qmask==-2 is for masked ones
+    The count array (same size as tth_min or tth_max) holds the number of pixel per reflection.
+
+    :param tth_array: 2D array with 2th positions
+    :param tth_min: 1D array with lower bound for each reflection, same size as tth_max
+    :param tth_max: 1D array with upper bound for each reflection, same size as tth_min
+    :param mask: marked pixel are invalid
+    :return: qmask, count
+    """
+    cdef Py_ssize_t nref = tth_min.size
+    assert tth_max.size == nref, "tth_max.size == tth_min.size"
+    cdef:
+        int i, r, n = tth_array.size
+        double tthv, lb, ub
+        double[::1] ctth = numpy.ascontiguousarray(tth_array.ravel(), dtype=numpy.float64)
+        double[::1] ctth_min = numpy.ascontiguousarray(tth_min.ravel(), dtype=numpy.float64)
+        double[::1] ctth_max = numpy.ascontiguousarray(tth_max.ravel(), dtype=numpy.float64)
+        int32_t[::1] qmask = numpy.zeros(n, dtype=numpy.int32)
+        int32_t[::1] count = numpy.zeros(nref, dtype=numpy.int32)
+        int8_t[::1] cmask
+        bint do_mask=False
+    if mask is not None:
+        assert mask.size == n, "mask.size == tth_array.size"
+        cmask = numpy.ascontiguousarray(mask.ravel(), dtype=numpy.int8)
+        do_mask = True
+
+    with nogil:
+        for i in range(n):
+            if do_mask and cmask[i]:
+                qmask[i] = -2
+                continue
+            tthv = ctth[i]
+            for r in range(nref):
+                lb = ctth_min[r]
+                ub = ctth_max[r]
+                if tthv < lb:
+                    qmask[i] = -1
+                    break
+                if lb <= tthv < ub:
+                    qmask[i] = r
+                    count[r] += 1
+                    break
+            else:
+                qmask[i] = -1
+
+    return numpy.asarray(qmask).reshape(tth_array.shape), numpy.asarray(count)

src/pyFAI/ext/meson.build

@@ -101,3 +101,6 @@ py.extension_module('_distortion', '_distortion.pyx',
 
 py.extension_module('fastcrc', ['fastcrc.pyx','src/crc32.c'],
         dependencies : [py_dep], install: true, subdir: 'pyFAI/ext')
+
+py.extension_module('mathutil', ['mathutil.pyx'],
+        dependencies : [py_dep], install: true, subdir: 'pyFAI/ext')

src/pyFAI/goniometer.py

@@ -34,7 +34,7 @@
 __contact__ = "Jerome.Kieffer@ESRF.eu"
 __license__ = "MIT"
 __copyright__ = "European Synchrotron Radiation Facility, Grenoble, France"
-__date__ = "25/04/2023"
+__date__ = "20/02/2024"
 __status__ = "development"
 __docformat__ = 'restructuredtext'
 
@@ -54,7 +54,7 @@
 from .utils import StringTypes
 from .multi_geometry import MultiGeometry
 from .units import CONST_hc, CONST_q
-
+from .ext.mathutil import build_qmask
 logger = logging.getLogger(__name__)
 
 try:
@@ -680,27 +680,26 @@ def extract_cp(self, max_rings=None, pts_per_deg=1.0, Imin=0):
         if max_rings is None:
             max_rings = tth.size
 
+        qmask, count = build_qmask(ttha, tth_min, tth_max, self.geometry_refinement.detector.mask)
+        mask2 = numpy.empty(qmask.shape, dtype=bool)
         ms = marchingsquares.MarchingSquaresMergeImpl(ttha,
                                                       mask=self.geometry_refinement.detector.mask,
                                                       use_minmax_cache=True)
         for i in range(tth.size):
             if rings >= max_rings:
                 break
-            mask = numpy.logical_and(ttha >= tth_min[i], ttha < tth_max[i])
-            if self.detector.mask is not None:
-                mask = numpy.logical_and(mask, numpy.logical_not(self.geometry_refinement.detector.mask))
-            size = mask.sum(dtype=int)
-            if (size > 0):
+            if count[i]:
                 rings += 1
-                sub_data = self.image.ravel()[numpy.where(mask.ravel())]
+                mask = qmask==i
+                sub_data = self.image[mask]
                 mean = sub_data.mean(dtype=numpy.float64)
                 std = sub_data.std(dtype=numpy.float64)
                 upper_limit = mean + std
-                mask2 = numpy.logical_and(self.image > upper_limit, mask)
+                numpy.logical_and(self.image > upper_limit, mask, out=mask2)
                 size2 = mask2.sum(dtype=int)
                 if size2 < 1000:
                     upper_limit = mean
-                    mask2 = numpy.logical_and(self.image > upper_limit, mask)
+                    numpy.logical_and(self.image > upper_limit, mask, out=mask2)
                     size2 = mask2.sum()
                 # length of the arc:
                 points = ms.find_pixels(tth[i])

src/pyFAI/massif.py

@@ -30,7 +30,7 @@
 __contact__ = "Jerome.Kieffer@ESRF.eu"
 __license__ = "MIT"
 __copyright__ = "European Synchrotron Radiation Facility, Grenoble, France"
-__date__ = "25/04/2023"
+__date__ = "20/02/2024"
 __status__ = "production"
 
 import sys
@@ -221,21 +221,26 @@ def peaks_from_area(self, mask, Imin=numpy.finfo(numpy.float64).min,
         :param seed: list of good guesses to start with
         :return: list of peaks [y,x], [y,x], ...]
         """
-        all_points = numpy.vstack(numpy.where(mask)).T
+
         res = []
         cnt = 0
         dmin2 = dmin * dmin
+
+        all_points = numpy.vstack(numpy.where(mask)).T
         if len(all_points) > 0:
             numpy.random.shuffle(all_points)
+
         if seed:
             seeds = numpy.array(list(seed))
             if len(seeds) > 0:
                 numpy.random.shuffle(seeds)
             all_points = numpy.concatenate((seeds, all_points))
+
+        msg = "[ %3i, %3i ] -> [ %.1f, %.1f ]"
         for idx in all_points:
             out = self.nearest_peak(idx)
             if out is not None:
-                msg = "[ %3i, %3i ] -> [ %.1f, %.1f ]"
+
                 logger.debug(msg, idx[1], idx[0], out[1], out[0])
                 p0, p1 = int(round(out[0])), int(round(out[1]))
                 if mask[p0, p1]:

src/pyFAI/test/test_utils_mathutil.py

@@ -32,7 +32,7 @@
 __contact__ = "Jerome.Kieffer@ESRF.eu"
 __license__ = "MIT"
 __copyright__ = "European Synchrotron Radiation Facility, Grenoble, France"
-__date__ = "05/09/2023"
+__date__ = "20/02/2024"
 
 import unittest
 import numpy
@@ -147,6 +147,15 @@ def test_interp_filter(self):
         z[w] = numpy.NaN
         self.assertLess(abs(y - utils.mathutil.interp_filter(z)).max(), 0.01, "error is small")
 
+    def test_is_far_from_group_cython(self):
+        from ..utils.mathutil import is_far_from_group_python, is_far_from_group_cython
+        rng = utilstest.UtilsTest.get_rng()
+        pt = rng.uniform(size=2)
+        pts = list(rng.uniform(size=(10,2)))
+        dst2 = rng.uniform()**2
+        ref = is_far_from_group_python(pt, pts, dst2)
+        res = is_far_from_group_cython(pt, pts, dst2)
+        self.assertEqual(ref, res, "cython implementation matches *is_far_from_group*")
 
 def suite():
     loader = unittest.defaultTestLoader.loadTestsFromTestCase

src/pyFAI/utils/mathutil.py

@@ -34,7 +34,7 @@
 __contact__ = "Jerome.Kieffer@ESRF.eu"
 __license__ = "MIT"
 __copyright__ = "European Synchrotron Radiation Facility, Grenoble, France"
-__date__ = "29/08/2023"
+__date__ = "20/02/2024"
 __status__ = "production"
 
 import logging
@@ -716,7 +716,7 @@ def roundfft(*args, **kwargs):
     return round_fft(*args, **kwargs)
 
 
-def is_far_from_group(pt, lst_pts, d2):
+def is_far_from_group_python(pt, lst_pts, d2):
     """
     Tells if a point is far from a group of points, distance greater than d2 (distance squared)
 
@@ -732,6 +732,12 @@ def is_far_from_group(pt, lst_pts, d2):
             return False
     return True
 
+try:
+    from ..ext.mathutil import is_far_from_group_cython
+except ImportError:
+    is_far_from_group = is_far_from_group_python
+else:
+    is_far_from_group = is_far_from_group_cython
 
 def rwp(obt, ref, scale=1.0):
     """Compute :math:`\\sqrt{\\sum \\frac{4\\cdot(obt-ref)^2}{(obt + ref)^2}}`.