From 52b946a6ec7236f9c7c82273422580582fa46381 Mon Sep 17 00:00:00 2001 From: Jerome Kieffer Date: Thu, 14 Nov 2024 11:55:24 +0100 Subject: [PATCH] Implement azimuthal medfilt with test in python --- src/pyFAI/engines/CSR_engine.py | 41 +++++----- src/pyFAI/test/meson.build | 1 + src/pyFAI/test/test_all.py | 4 +- src/pyFAI/test/test_medfilt_engine.py | 105 ++++++++++++++++++++++++++ 4 files changed, 129 insertions(+), 22 deletions(-) create mode 100644 src/pyFAI/test/test_medfilt_engine.py diff --git a/src/pyFAI/engines/CSR_engine.py b/src/pyFAI/engines/CSR_engine.py index 7a5facc12..62b436428 100644 --- a/src/pyFAI/engines/CSR_engine.py +++ b/src/pyFAI/engines/CSR_engine.py @@ -26,7 +26,7 @@ __contact__ = "Jerome.Kieffer@ESRF.eu" __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): diff --git a/src/pyFAI/test/meson.build b/src/pyFAI/test/meson.build index 61c66b51b..fd0090226 100644 --- a/src/pyFAI/test/meson.build +++ b/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 diff --git a/src/pyFAI/test/test_all.py b/src/pyFAI/test/test_all.py index 41b6de431..d738acab3 100644 --- a/src/pyFAI/test/test_all.py +++ b/src/pyFAI/test/test_all.py @@ -32,7 +32,7 @@ __contact__ = "jerome.kieffer@esrf.eu" __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 diff --git a/src/pyFAI/test/test_medfilt_engine.py b/src/pyFAI/test/test_medfilt_engine.py new file mode 100644 index 000000000..12accfb18 --- /dev/null +++ b/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 (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. +from pygments.lexers import solidity + +"""Test suites for median filtering engines""" + +__author__ = "Jérôme Kieffer" +__contact__ = "Jerome.Kieffer@ESRF.eu" +__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())