Implement azimuthal medfilt with test in python

src/pyFAI/engines/CSR_engine.py

@@ -26,7 +26,7 @@
 __contact__ = "[email protected]"
 __license__ = "MIT"
 __copyright__ = "European Synchrotron Radiation Facility, Grenoble, France"
-__date__ = "13/11/2024"
+__date__ = "14/11/2024"
 __status__ = "development"
 
 from collections.abc import Iterable
@@ -395,7 +395,7 @@ def medfilt(self, data, dark=None, dummy=None, delta_dummy=None,
                 variance=None, dark_variance=None,
                 flat=None, solidangle=None, polarization=None, absorption=None,
                 safe=True, error_model=None,
-                normalization_factor=1.0, percentile=50
+                normalization_factor=1.0, quantile=0.5
                 ):
         """
         Perform a median-filter/quantile mean in azimuthal space.
@@ -429,27 +429,24 @@ def medfilt(self, data, dark=None, dummy=None, delta_dummy=None,
         :param safe: Unused in this implementation
         :param error_model: Enum or str, "azimuthal" or "poisson"
         :param normalization_factor: divide raw signal by this value
-        :param percentile: which percentile use for cutting out
-                           percentil can be a 2-tuple to specify a region to
-                           average out
+        :param quantile: which percentile/100 use for cutting out quantil.
+                         can be a 2-tuple to specify a region to average out.
+                         By default, takes the median
         :return: namedtuple with "position intensity error signal variance normalization count"
 
         """
-        if isinstance(percentile, Iterable):
-            q_start = 1e-2 * min(percentile)
-            q_stop = 1e-2 * max(percentile)
+        if isinstance(quantile, Iterable):
+            q_start = min(quantile)
+            q_stop = max(quantile)
         else:
-            q_stop = q_start = 1e-2*percentile
+            q_stop = q_start = quantile
 
-        shape = data.shape
         indptr = self._csr.indptr
         indices = self._csr.indices
-
+        csr_data = self._csr.data
+        csr_data2 = self._csr2.data
 
         error_model = ErrorModel.parse(error_model)
-        if error_model is ErrorModel.NO:
-            logger.error("No variance propagation is incompatible with sigma-clipping. Using `azimuthal` model !")
-            error_model = ErrorModel.AZIMUTHAL
 
         prep = preproc(data,
                        dark=dark,
@@ -474,9 +471,9 @@ def medfilt(self, data, dark=None, dummy=None, delta_dummy=None,
 
         work0 = numpy.zeros((indices.size,4), dtype=numpy.float32)
         work0[:, 0] = pixels[:, 0]/ pixels[:, 2]
-        work0[:, 1] = pixels[:, 0] * self._csr.data
-        work0[:, 2] = pixels[:, 0] * self._csr2.data
-        work0[:, 3] = pixels[:, 0] * self._csr.data
+        work0[:, 1] = pixels[:, 0] * csr_data
+        work0[:, 2] = pixels[:, 1] * csr_data2
+        work0[:, 3] = pixels[:, 2] * csr_data
         work1 = work0.view(mf_dtype).ravel()
 
         size = indptr.size-1
@@ -499,16 +496,18 @@ def medfilt(self, data, dark=None, dummy=None, delta_dummy=None,
             norm[i] = tmp["norm"].sum(dtype=numpy.float64)
             norm2[i] = (tmp["norm"]**2).sum(dtype=numpy.float64)
 
-        avg = signal / norm
-        std = numpy.sqrt(variance / norm2)
-        sem = numpy.sqrt(variance) / norm
+        with warnings.catch_warnings():
+            warnings.simplefilter("ignore")
+            avg = signal / norm
+            std = numpy.sqrt(variance / norm2)
+            sem = numpy.sqrt(variance) / norm
         # mask out remaining NaNs
         msk = norm <= 0
         avg[msk] = self.empty
         std[msk] = self.empty
         sem[msk] = self.empty
 
-        return Integrate1dtpl(self.bin_centers, avg, std, signal, variance, norm, cnt, std, sem, norm2)
+        return Integrate1dtpl(self.bin_centers, avg, sem, signal, variance, norm, cnt, std, sem, norm2)
 
     @property
     def check_mask(self):

src/pyFAI/test/meson.build

@@ -56,6 +56,7 @@ py.install_sources(
     'test_uncertainties.py',
     'test_watershed.py',
     'test_worker.py',
+    'test_medfilt_engine.py',
     'utilstest.py'],
   pure: false,   # Will be installed next to binaries
   subdir: 'pyFAI/test'  # Folder relative to site-packages to install to

src/pyFAI/test/test_all.py

@@ -32,7 +32,7 @@
 __contact__ = "[email protected]"
 __license__ = "MIT"
 __copyright__ = "European Synchrotron Radiation Facility, Grenoble, France"
-__date__ = "30/10/2024"
+__date__ = "14/11/2024"
 
 import sys
 import unittest
@@ -96,6 +96,7 @@
 from . import test_uncertainties
 from . import test_ring_extraction
 from . import test_fiber_integrator
+from . import test_medfilt_engine
 
 logger = logging.getLogger(__name__)
 
@@ -158,6 +159,7 @@ def suite():
     testsuite.addTest(test_uncertainties.suite())
     testsuite.addTest(test_ring_extraction.suite())
     testsuite.addTest(test_fiber_integrator.suite())
+    testsuite.addTest(test_medfilt_engine.suite())
     return testsuite
 
 

src/pyFAI/test/test_medfilt_engine.py

@@ -0,0 +1,105 @@
+#!/usr/bin/env python
+# coding: utf-8
+#
+#    Project: Azimuthal integration
+#             https://github.com/silx-kit/pyFAI
+#
+#    Copyright (C) 2015-2018 European Synchrotron Radiation Facility, Grenoble, France
+#
+#    Principal author:       Jérôme Kieffer ([email protected])
+#
+# 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.
+from pygments.lexers import solidity
+
+"""Test suites for median filtering engines"""
+
+__author__ = "Jérôme Kieffer"
+__contact__ = "[email protected]"
+__license__ = "MIT"
+__copyright__ = "European Synchrotron Radiation Facility, Grenoble, France"
+__date__ = "14/11/2024"
+
+import unittest
+import numpy
+import logging
+logger = logging.getLogger(__name__)
+from .utilstest import UtilsTest
+import fabio
+from .. import load
+
+
+class TestMedfilt(unittest.TestCase):
+    """Test Azimuthal median filtering results
+    """
+
+    @classmethod
+    def setUpClass(cls)->None:
+        super(TestMedfilt, cls).setUpClass()
+        cls.method = ["full", "csr", "python"]
+        cls.img = fabio.open(UtilsTest.getimage("mock.tif")).data
+        cls.ai = load({ "dist": 0.1,
+                        "poni1":0.03,
+                        "poni2":0.03,
+                        "detector": "Detector",
+                        "detector_config": {"pixel1": 1e-4,
+                                            "pixel2": 1e-4,
+                                            "max_shape": [500, 600],
+                                            "orientation": 3}})
+        cls.npt = 100
+
+    @classmethod
+    def tearDownClass(cls)->None:
+        super(TestMedfilt, cls).tearDownClass()
+        cls.method = cls.img =cls.ai =cls.npt =None
+
+
+    def test_python(self):
+        print(self.ai)
+        method = tuple(self.method)
+        ref = self.ai.integrate1d(self.img, self.npt, unit="2th_rad", method=method, error_model="poisson")
+        print(ref.method)
+        engine = self.ai.engines[ref.method].engine
+        obt = engine.medfilt(self.img,
+                             solidangle=self.ai.solidAngleArray(),
+                             quantile=(0,1),  # taking all Like this it works like a normal mean
+                             error_model="poisson")
+
+        self.assertTrue(numpy.allclose(ref.radial, obt.position), "radial matches")
+        self.assertTrue(numpy.allclose(ref.sum_signal, obt.signal), "signal matches")
+        self.assertTrue(numpy.allclose(ref.sum_variance, obt.variance), "variance matches")
+        self.assertTrue(numpy.allclose(ref.sum_normalization, obt.normalization), "normalization matches")
+        self.assertTrue(numpy.allclose(ref.sum_normalization2, obt.norm_sq), "norm_sq matches")
+
+        self.assertTrue(numpy.allclose(ref.intensity, obt.intensity), "intensity matches")
+        self.assertTrue(numpy.allclose(ref.sigma, obt.sigma), "sigma matches")
+        self.assertTrue(numpy.allclose(ref.std, obt.std), "std matches")
+        self.assertTrue(numpy.allclose(ref.sem, obt.sem), "sem matches")
+
+
+
+def suite():
+    loader = unittest.defaultTestLoader.loadTestsFromTestCase
+    testsuite = unittest.TestSuite()
+    testsuite.addTest(loader(TestMedfilt))
+    return testsuite
+
+
+if __name__ == '__main__':
+    runner = unittest.TextTestRunner()
+    runner.run(suite())