diff --git a/src/pyFAI/test/test_utils_mathutil.py b/src/pyFAI/test/test_utils_mathutil.py index 6529aa7a7..2e6465ee3 100644 --- a/src/pyFAI/test/test_utils_mathutil.py +++ b/src/pyFAI/test/test_utils_mathutil.py @@ -4,7 +4,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 (Jerome.Kieffer@ESRF.eu) # @@ -32,7 +32,7 @@ __contact__ = "Jerome.Kieffer@ESRF.eu" __license__ = "MIT" __copyright__ = "European Synchrotron Radiation Facility, Grenoble, France" -__date__ = "16/05/2024" +__date__ = "21/01/2025" import unittest import numpy @@ -40,7 +40,12 @@ import logging from . import utilstest logger = logging.getLogger(__name__) +import fabio from .. import utils +from .. import load +from .. import calibrant +from ..utils import mathutil + import scipy.ndimage @@ -148,19 +153,24 @@ def test_interp_filter(self): 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) + ref = mathutil.is_far_from_group_python(pt, pts, dst2) + res = mathutil.is_far_from_group_cython(pt, pts, dst2) self.assertEqual(ref, res, "cython implementation matches *is_far_from_group*") def test_allclose_mod(self): - from ..utils.mathutil import allclose_mod - self.assertTrue(allclose_mod(numpy.arctan2(+1e-10, -1), numpy.arctan2(-1e-10, -1)),"angles matches modulo 2pi") - + self.assertTrue(mathutil.allclose_mod(numpy.arctan2(+1e-10, -1), numpy.arctan2(-1e-10, -1)),"angles matches modulo 2pi") + + def test_quality_of_fit(self): + img = fabio.open(utilstest.UtilsTest.getimage("Pilatus1M.edf")).data + ai = load(utilstest.UtilsTest.getimage("Pilatus1M.poni")) + 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") def suite(): loader = unittest.defaultTestLoader.loadTestsFromTestCase diff --git a/src/pyFAI/utils/mathutil.py b/src/pyFAI/utils/mathutil.py index e1c00e7e8..e852ece36 100644 --- a/src/pyFAI/utils/mathutil.py +++ b/src/pyFAI/utils/mathutil.py @@ -34,7 +34,7 @@ __contact__ = "Jerome.Kieffer@ESRF.eu" __license__ = "MIT" __copyright__ = "European Synchrotron Radiation Facility, Grenoble, France" -__date__ = "13/01/2025" +__date__ = "21/01/2025" __status__ = "production" import logging @@ -44,6 +44,7 @@ import numpy import time import scipy.ndimage +from scipy.signal import peak_widths from .decorators import deprecated try: @@ -939,3 +940,65 @@ def allclose_mod(a, b, modulo=2*numpy.pi, **kwargs): """ di = numpy.minimum((a-b)%modulo, (b-a)%modulo) return numpy.allclose(modulo*0.5, (di+modulo*0.5), **kwargs) + + +def quality_of_fit(img, ai, calibrant, + npt_rad=1000, npt_azim=360, + unit="q_nm^-1", + method=("full", "csr", "cython"), + empty = numpy.nan, rings=None): + """Provide an indicator for the quality of fit of a given geometry for an image + + :param img: 2D image with a calibration image (containing rings) + :param ai: azimuthal integrator object (instance of pyFAI.integrator.azimuthal.AzimuthalIntegrator) + :param calibrant: calibration object, instance of pyFAI.calibrant.Calibrant + :param npt_rad: int with the number of radial bins + :param npt_azim: int with the number of azimuthal bins + :param unit: typically "2th_deg" or "q_nm^-1", the quality of fit should be largely independant from the space. + :param method: integration method + :param empty: value of the empy bins, discarded values + :param rings: list of rings to evaluate (0-based) + :return: QoF indicator, similar to reduced χ², the smaller, the better + """ + + ai.empty = empty + q_theo = calibrant.get_peaks(unit=unit) + res = ai.integrate2d(img, npt_rad, npt_azim, method=method, unit=unit) + if rings is None: + rings = list(range(len(calibrant.get_2th()))) + q_theo = q_theo[rings] + idx_theo = abs(numpy.add.outer(res.radial,-q_theo)).argmin(axis=0) + idx_maxi = numpy.empty((res.azimuthal.size, q_theo.size))+numpy.nan + idx_fwhm = numpy.empty((res.azimuthal.size, q_theo.size))+numpy.nan + signal = res.intensity + gradient = numpy.gradient(signal, axis=-1) + minima = numpy.where(numpy.logical_and(gradient[:,:-1]<0, gradient[:,1:]>=0)) + maxima = numpy.where(numpy.logical_and(gradient[:,:-1]>0, gradient[:,1:]<0)) + for idx in range(res.azimuthal.size): + for ring in rings: + q_th = q_theo[ring] + idx_th = idx_theo[ring] + if (q_th<=res.radial[0]) or (q_th>=res.radial[-1]): + continue + maxi = maxima[1][maxima[0]==idx] + mini = minima[1][minima[0]==idx] + idx_max = maxi[abs(maxi-idx_th).argmin()] + idx_inf = mini[miniidx_max] + if idx_sup.size: + idx_sup = idx_sup[0] + if idx_inf< idx_th< idx_sup: + sub = signal[idx, idx_inf:idx_sup+1] - numpy.linspace(signal[idx, idx_inf],signal[idx, idx_sup], 1+idx_sup-idx_inf) + com = (sub*numpy.linspace(idx_inf, idx_sup, 1+idx_sup-idx_inf)).sum()/sub.sum() + if numpy.isfinite(com): + width = peak_widths(sub, [numpy.argmax(sub)])[0][0] + if width==0: + print(f" #{idx},{ring}: {idx_inf} < th:{idx_th} max:{idx_max} com:{com:.3f} < {idx_sup}; fwhm={width}") + print(signal[idx, idx_inf:idx_sup+1]) + print(sub) + else: + idx_fwhm[idx, ring] = width + idx_maxi[idx, ring] = idx_max + return numpy.nanmean((2.355*(idx_maxi-idx_theo)/idx_fwhm)**2)