Merge pull request #2467 from kif/2644_automatic_dynamic_mask

automatic dynamic mask for Dectris detetors

Author: kif

src/pyFAI/detectors/_common.py

@@ -949,6 +949,25 @@ def calc_mask(self):
 #        logger.debug("Detector.calc_mask is not implemented for generic detectors")
         return None
 
+    def get_dummies(self, img):
+        """Calculate the actual dummy value from dtype of the img
+
+        :param img: numpy array (or actually its dtype)
+        :return: actual (dummy, delta_dummy) values as data_d (i.e. float32)
+        """
+        if self.dummy is None:
+            return None, None
+        if isinstance(img, numpy.ndarray):
+            dtype = numpy.dtype(img.dtype)
+        else:
+            dtype = numpy.dtype(img)
+        actual_dummy = numpy.float32(numpy.dtype(img.dtype).type(numpy.int64(self.dummy)))
+        if self.delta_dummy is None:
+            actual_delta_dummy = numpy.finfo("float32").eps
+        else:
+            actual_delta_dummy = numpy.float32(self.delta_dummy)
+        return actual_dummy, actual_delta_dummy
+
     def dynamic_mask(self, img):
         """Calculate the dynamic mask for the given image.
 

src/pyFAI/ext/CSR_common.pxi

@@ -3,7 +3,7 @@
 #    Project: Azimuthal integration
 #             https://github.com/silx-kit/pyFAI
 #
-#    Copyright (C) 2015-2024 European Synchrotron Radiation Facility, Grenoble, France
+#    Copyright (C) 2015-2025 European Synchrotron Radiation Facility, Grenoble, France
 #
 #    Principal author:       Jérôme Kieffer ([email protected])
 #
@@ -29,7 +29,7 @@
 
 __author__ = "Jérôme Kieffer"
 __contact__ = "[email protected]"
-__date__ = "05/12/2024"
+__date__ = "12/03/2025"
 __status__ = "stable"
 __license__ = "MIT"
 
@@ -347,7 +347,7 @@ cdef class CsrIntegrator(object):
             bint do_azimuthal_variance = error_model is ErrorModel.AZIMUTHAL
             bint do_variance = error_model is not ErrorModel.NO
         assert weights.size == self.input_size, "weights size"
-        empty = dummy if dummy is not None else self.empty
+        empty = self.empty if dummy is None else dummy
         #Call the preprocessor ...
         preproc4 = preproc(weights.ravel(),
                            dark=dark,

src/pyFAI/integrator/azimuthal.py

@@ -30,7 +30,7 @@
 __contact__ = "[email protected]"
 __license__ = "MIT"
 __copyright__ = "European Synchrotron Radiation Facility, Grenoble, France"
-__date__ = "10/03/2025"
+__date__ = "12/03/2025"
 __status__ = "stable"
 __docformat__ = 'restructuredtext'
 
@@ -128,7 +128,12 @@ def integrate1d(self, data, npt, filename=None,
         method = self._normalize_method(method, dim=1, default=self.DEFAULT_METHOD_1D)
         assert method.dimension == 1
         unit = units.to_unit(unit)
-        empty = numpy.float32(dummy) if dummy is not None else self._empty
+        if dummy is None:
+            dummy, delta_dummy = self.detector.get_dummies(data)
+            empty = self._empty if dummy is None else dummy
+        else:
+            dummy = empty = numpy.float32(dummy)
+            delta_dummy = None if delta_dummy is None else numpy.float32(delta_dummy)
         shape = data.shape
         pos0_scale = unit.scale
 
@@ -605,7 +610,12 @@ def integrate_radial(self, data, npt, npt_rad=100,
         sum_normalization = res._sum_normalization.sum(axis=-1)
 
         mask = numpy.where(count == 0)
-        empty = numpy.float32(dummy) if dummy is not None else self._empty
+        if dummy is None:
+            dummy, delta_dummy = self.detector.get_dummies(data)
+            empty = self._empty if dummy is None else dummy
+        else:
+            dummy = empty = numpy.float32(dummy)
+            delta_dummy = None if delta_dummy is None else numpy.float32(delta_dummy)
         intensity = sum_signal / sum_normalization
         intensity[mask] = empty
 
@@ -703,7 +713,13 @@ def integrate2d_ng(self, data, npt_rad, npt_azim=360,
         space = (radial_unit.space, azimuth_unit.space)
         pos0_scale = radial_unit.scale
         pos1_scale = azimuth_unit.scale
-        empty = numpy.float32(dummy) if dummy is not None else self._empty
+        if dummy is None:
+            dummy, delta_dummy = self.detector.get_dummies(data)
+            empty = self._empty if dummy is None else dummy
+        else:
+            dummy = empty = numpy.float32(dummy)
+            delta_dummy = None if delta_dummy is None else numpy.float32(delta_dummy)
+
         if mask is None:
             has_mask = "from detector"
             mask = self.mask
@@ -1103,7 +1119,7 @@ def integrate2d_ng(self, data, npt_rad, npt_azim=360,
                                          dummy=dummy,
                                          delta_dummy=delta_dummy,
                                          normalization_factor=normalization_factor,
-                                         empty=self._empty,
+                                         empty=empty,
                                          variance=variance,
                                          dark_variance=None,
                                          error_model=error_model,

src/pyFAI/test/test_azimuthal_integrator.py

@@ -390,6 +390,22 @@ def test_guess_polarization(self):
         ai = AzimuthalIntegrator.sload(UtilsTest.getimage("Eiger4M.poni"))
         self.assertLess(abs(ai.guess_polarization(img) - 0.5), 0.1)
 
+    def test_detector_dynamic_mask(self):
+        """Check if the dummy value from detector is actually used #2466"""
+        det = detector_factory("Pilatus CdTe 300k")
+        corners = det.get_pixel_corners()
+        config = {"detector": det.name,
+          "dist": 1,
+          "poni1": corners[...,1].max()/4,
+          "poni2": corners[...,2].max()/4}
+        ai1 = AzimuthalIntegrator.sload(config)
+        img1 = det.mask.astype("uint16")*(65535-1)+1
+        img1[300:400, 300:400] = 65535
+        result = ai1.integrate1d(img1, 100, unit="r_mm",
+                                 method=("no", "histogram", "cython"),
+                                 correctSolidAngle=False)
+        self.assertTrue(numpy.allclose(result.intensity, numpy.ones_like(result.radial)))
+
 
 class TestSaxs(unittest.TestCase):
     saxsPilatus = "bsa_013_01.edf"

src/pyFAI/test/test_csr.py

@@ -195,16 +195,16 @@ def test_2d_nosplit(self):
         result_nosplit = self.ai.integrate2d(self.data, self.N, unit="2th_deg", method=("no", "csr", "cython"))
         self.assertTrue(numpy.allclose(result_histo.radial, result_nosplit.radial), " 2Th are the same")
         self.assertTrue(numpy.allclose(result_histo.azimuthal, result_nosplit.azimuthal, atol=1e-5), " Chi are the same")
-        if False:
-            print(result_histo.method, result_histo.method_called)
-            print(result_histo.sum_signal.min(), result_histo.sum_signal.max(), result_histo.sum_signal.mean(), result_histo.sum_signal.std())
-            print(result_histo.sum_normalization.min(), result_histo.sum_normalization.max(), result_histo.sum_normalization.mean(), result_histo.sum_normalization.std())
-            print(result_histo.intensity.min(), result_histo.intensity.max(), result_histo.intensity.mean(), result_histo.intensity.std())
-            print(result_nosplit.intensity.min(), result_nosplit.intensity.max(), result_nosplit.intensity.mean(), result_nosplit.intensity.std())
-            print(result_nosplit.intensity)
+        # print("result_histo", result_histo.method, result_histo.method_called)
+        # print("result_histo", result_histo.sum_signal.min(), result_histo.sum_signal.max(), result_histo.sum_signal.mean(), result_histo.sum_signal.std())
+        # print("result_histo",result_histo.sum_normalization.min(), result_histo.sum_normalization.max(), result_histo.sum_normalization.mean(), result_histo.sum_normalization.std())
+        # print("result_histo",result_histo.intensity.min(), result_histo.intensity.max(), result_histo.intensity.mean(), result_histo.intensity.std())
+        # print("result_nosplit",result_nosplit.intensity.min(), result_nosplit.intensity.max(), result_nosplit.intensity.mean(), result_nosplit.intensity.std())
+        # print("result_nosplit",result_nosplit.method)
         error = (result_histo.intensity - result_nosplit.intensity)
         logger.debug("ref: %s; obt: %s", result_histo.intensity.shape, result_nosplit.intensity.shape)
         logger.debug("error mean: %s, std: %s", error.mean(), error.std())
+
         self.assertLess(error.mean(), 1e-3, "img are almost the same")
         self.assertLess(error.std(), 3, "img are almost the same")
 

src/pyFAI/test/test_geometry.py

@@ -34,7 +34,7 @@
 __contact__ = "[email protected]"
 __license__ = "MIT"
 __copyright__ = "European Synchrotron Radiation Facility, Grenoble, France"
-__date__ = "20/02/2025"
+__date__ = "12/03/2025"
 
 import unittest
 import random
@@ -559,18 +559,21 @@ class TestCalcFrom(unittest.TestCase):
     def test_calcfrom12d(self):
         det = detector_factory("pilatus300k")
         ai = AzimuthalIntegrator(0.1, 0.05, 0.04, detector=det)
-        prof_1d = ai.integrate1d_ng(UtilsTest.get_rng().random(det.shape), 200, unit="2th_deg")
-        sig = numpy.sinc(prof_1d.radial * 10) ** 2
-        img1 = ai.calcfrom1d(prof_1d.radial, sig, dim1_unit="2th_deg", mask=det.mask, dummy=-1)
+        img0 = UtilsTest.get_rng().random(det.shape)
+        prof_1d = ai.integrate1d_ng(img0, 200, unit="2th_deg")
+        sig = 1e6*numpy.sinc(prof_1d.radial / 10) ** 2
+        img1 = ai.calcfrom1d(prof_1d.radial, sig, dim1_unit="2th_deg",
+                            mask=det.mask, dummy=-1)
         new_prof_1d = ai.integrate1d_ng(img1, 200, unit="2th_deg")
         delta = abs((new_prof_1d.intensity - sig)).max()
-        self.assertLess(delta, 2e-3, "calcfrom1d works delta=%s" % delta)
+        self.assertLess(delta, 600, "calcfrom1d works delta=%s" % delta)
         prof_2d = ai.integrate2d(img1, 400, 360, unit="2th_deg")
         img2 = ai.calcfrom2d(prof_2d.intensity, prof_2d.radial, prof_2d.azimuthal,
                              mask=det.mask,
                              dim1_unit="2th_deg", correctSolidAngle=True, dummy=-1)
-        delta2 = abs(img2 - img1).max()
-        self.assertLess(delta2, 1e-3, "calcfrom2d works delta=%s" % delta2)
+        delta2 = img2 - img1
+        self.assertLess(abs(delta2.mean()), 2, "calcfrom2d works delta.mean=%s" % abs(delta2.mean()))
+        self.assertLess(delta2.std(), 100, "calcfrom2d works delta.std=%s" % delta2.std())
 
 
 class TestBugRegression(unittest.TestCase):

src/pyFAI/test/test_integrate.py

@@ -4,7 +4,7 @@
 #    Project: Azimuthal integration
 #             https://github.com/silx-kit/pyFAI
 #
-#    Copyright (C) 2015-2022 European Synchrotron Radiation Facility, Grenoble, France
+#    Copyright (C) 2015-2025 European Synchrotron Radiation Facility, Grenoble, France
 #
 #    Principal author:       Jérôme Kieffer ([email protected])
 #
@@ -31,7 +31,7 @@
 __contact__ = "[email protected]"
 __license__ = "MIT"
 __copyright__ = "European Synchrotron Radiation Facility, Grenoble, France"
-__date__ = "21/05/2024"
+__date__ = "12/03/2025"
 
 import contextlib
 import os
@@ -190,7 +190,12 @@ def setUpClass(cls):
         cls.img = UtilsTest.getimage("Pilatus1M.edf")
         with fabio.open(cls.img) as fimg:
             cls.data = fimg.data
-        cls.ai = AzimuthalIntegrator(1.58323111834, 0.0334170169115, 0.0412277798782, 0.00648735642526, 0.00755810191106, 0.0, detector=Pilatus1M())
+        class DummyLessPilatus(Pilatus1M):
+            DUMMY = None
+            DELTA_DUMMY = None
+
+        cls.ai = AzimuthalIntegrator(1.58323111834, 0.0334170169115, 0.0412277798782, 0.00648735642526, 0.00755810191106, 0.0,
+                                     detector=DummyLessPilatus())
         cls.ai.wavelength = 1e-10
         cls.Rmax = 30
         cls.delta_pos_azim_max = 0.28

src/pyFAI/test/test_mask.py

@@ -4,7 +4,7 @@
 #    Project: Azimuthal integration
 #             https://github.com/silx-kit/pyFAI
 #
-#    Copyright (C) 2015-2018 European Synchrotron Radiation Facility, Grenoble, France
+#    Copyright (C) 2015-2025 European Synchrotron Radiation Facility, Grenoble, France
 #
 #    Principal author:       Jérôme Kieffer ([email protected])
 #
@@ -32,7 +32,7 @@
 __contact__ = "[email protected]"
 __license__ = "MIT"
 __copyright__ = "European Synchrotron Radiation Facility, Grenoble, France"
-__date__ = "21/05/2024"
+__date__ = "12/03/2024"
 
 import unittest
 import numpy
@@ -48,21 +48,25 @@
 
 
 class TestMask(unittest.TestCase):
-    dataFile = "testMask.edf"
-    poniFile = "Pilatus1M.poni"
 
-    def setUp(self):
+    @classmethod
+    def setUpClass(cls):
         """Download files"""
-        self.dataFile = UtilsTest.getimage(self.__class__.dataFile)
-        self.poniFile = UtilsTest.getimage(self.__class__.poniFile)
-        self.ai = load(self.poniFile)
-        with fabio.open(self.dataFile) as fimg:
-            self.data = fimg.data
-        self.mask = self.data < 0
-
-    def tearDown(self):
-        unittest.TestCase.tearDown(self)
-        self.dataFile = self.data = self.ai = self.mask = self.poniFile = None
+        cls.dataFile = UtilsTest.getimage("testMask.edf")
+        cls.poniFile = UtilsTest.getimage("Pilatus1M.poni")
+        cls.ai = load(cls.poniFile)
+        # hack to prevent dynamic masking: this is a Pilatus without dummy values
+        class MyPilatus(cls.ai.detector.__class__):
+            DUMMY=None
+            DELTA_DUMMY=None
+        cls.ai.detector = MyPilatus()
+        with fabio.open(cls.dataFile) as fimg:
+            cls.data = fimg.data
+        cls.mask = cls.data < 0
+
+    @classmethod
+    def tearDownClass(cls):
+        cls.dataFile = cls.data = cls.ai = cls.mask = cls.poniFile = None
 
     def test_mask_hist(self):
         """

src/pyFAI/test/test_utils_mathutil.py

@@ -170,7 +170,7 @@ def test_quality_of_fit(self):
         cal = calibrant.get_calibrant("AgBh")
         cal.wavelength = ai.wavelength
         res = mathutil.quality_of_fit(img, ai, cal, rings=[0,1], npt_azim=36, npt_rad=100)
-        self.assertLess(res, 0.3, "Fit of good quality")
+        self.assertLess(res, 0.31, "Fit of good quality")
 
     def test_nan_equal(self):
         nan_equal = mathutil.nan_equal