diff --git a/src/pyFAI/integrator/azimuthal.py b/src/pyFAI/integrator/azimuthal.py index 0c209a2ba..c53bf839e 100644 --- a/src/pyFAI/integrator/azimuthal.py +++ b/src/pyFAI/integrator/azimuthal.py @@ -30,7 +30,7 @@ __contact__ = "Jerome.Kieffer@ESRF.eu" __license__ = "MIT" __copyright__ = "European Synchrotron Radiation Facility, Grenoble, France" -__date__ = "10/10/2024" +__date__ = "05/12/2024" __status__ = "stable" __docformat__ = 'restructuredtext' @@ -49,6 +49,7 @@ from ..io.ponifile import PoniFile error = None from ..method_registry import IntegrationMethod +from ..utils.decorators import deprecated # from .load_engines import ocl_sort #ocl_azim_csr, ocl_azim_lut, , histogram, splitBBox, \ @@ -1147,8 +1148,8 @@ def integrate2d_ng(self, data, npt_rad, npt_azim=360, integrate2d = _integrate2d_ng = integrate2d_ng - - def medfilt1d(self, data, npt_rad=1024, npt_azim=512, + @deprecated(since_version="2024.12", only_once=True, replacement="medfilt1d_ng", deprecated_since="2024.12") + def medfilt1d_legacy(self, data, npt_rad=1024, npt_azim=512, correctSolidAngle=True, radial_range=None, azimuth_range=None, polarization_factor=None, dark=None, flat=None, @@ -1285,6 +1286,305 @@ def medfilt1d(self, data, npt_rad=1024, npt_azim=512, result._set_normalization_factor(normalization_factor) return result + medfilt1d = medfilt1d_legacy + + def medfilt1d_ng(self, data, + npt=1024, + correctSolidAngle=True, + polarization_factor=None, + variance=None, + error_model=ErrorModel.NO, + radial_range=None, + azimuth_range=None, + dark=None, + flat=None, + absorption=None, + method=("full", "csr", "cython"), + unit=units.Q, + percentile=50, + dummy=None, + delta_dummy=None, + mask=None, + normalization_factor=1.0, + metadata=None, + safe=True, + **kwargs): + """Performs iteratively the 1D integration with variance propagation + and performs a sigm-clipping at each iteration, i.e. + all pixel which intensity differs more than thres*std is + discarded for next iteration. + + Keep only pixels with intensty: + + ``|I - | < thres * σ(I)`` + + This enforces a symmetric, bell-shaped distibution (i.e. gaussian-like) + and is very good at extracting background or amorphous isotropic scattering + out of Bragg peaks. + + :param data: input image as numpy array + :param npt_rad: number of radial points + :param bool correctSolidAngle: correct for solid angle of each pixel if True + :param float polarization_factor: polarization factor between: + -1 (vertical) + +1 (horizontal). + - 0 for circular polarization or random, + - None for no correction, + - True for using the former correction + :param radial_range: The lower and upper range of the radial unit. If not provided, range is simply (data.min(), data.max()). Values outside the range are ignored. + :type radial_range: (float, float), optional + :param azimuth_range: The lower and upper range of the azimuthal angle in degree. If not provided, range is simply (data.min(), data.max()). Values outside the range are ignored. + :type azimuth_range: (float, float), optional + + :param ndarray dark: dark noise image + :param ndarray flat: flat field image + :param ndarray absorption: Detector absorption (image) + :param ndarray variance: the variance of the signal + :param str error_model: can be "poisson" to assume a poissonian detector (variance=I) or "azimuthal" to take the std² in each ring (better, more expenive) + :param unit: unit to be used for integration + :param IntegrationMethod method: IntegrationMethod instance or 3-tuple with (splitting, algorithm, implementation) + :param percentile: which percentile use for cutting out. + percentil can be a 2-tuple to specify a region to average out, + like: (25,75) to average the second and third quartile. + :param mask: masked out pixels array + :param float normalization_factor: Value of a normalization monitor + :param metadata: any other metadata, + :type metadata: JSON serializable dict + :param safe: set to False to skip some tests + :return: Integrate1D like result like + + The difference with the previous `medfilt_legacy` implementation is that there is no 2D regrouping. + """ + for k in kwargs: + if k == "npt_azim": + logger.warning("'npt_azim' argument is not used in sigma_clip_ng as not 2D intergration is performed anymore") + else: + logger.warning("Got unknown argument %s %s", k, kwargs[k]) + + error_model = ErrorModel.parse(error_model) + if variance is not None: + assert variance.size == data.size + error_model = ErrorModel.VARIANCE + + unit = units.to_unit(unit) + if radial_range: + radial_range = tuple(radial_range[i] / unit.scale for i in (0, -1)) + if azimuth_range is not None: + azimuth_range = self.normalize_azimuth_range(azimuth_range) + try: + quant_min = min(percentile)/100 + quant_max = max(percentile)/100 + except: + quant_min = quant_max = percentile/100.0 + + method = self._normalize_method(method, dim=1, default=self.DEFAULT_METHOD_1D) + + if mask is None: + has_mask = "from detector" + mask = self.mask + mask_crc = self.detector.get_mask_crc() + if mask is None: + has_mask = False + mask_crc = None + else: + has_mask = "user provided" + mask = numpy.ascontiguousarray(mask) + mask_crc = crc32(mask) + + if dark is None: + dark = self.detector.darkcurrent + if dark is None: + has_dark = False + else: + has_dark = "from detector" + else: + has_dark = "provided" + + if flat is None: + flat = self.detector.flatfield + if dark is None: + has_flat = False + else: + has_flat = "from detector" + else: + has_flat = "provided" + + if correctSolidAngle: + solidangle = self.solidAngleArray(data.shape, correctSolidAngle) + else: + solidangle = None + + if polarization_factor is None: + polarization = polarization_crc = None + else: + polarization, polarization_crc = self.polarization(data.shape, polarization_factor, with_checksum=True) + + if (method.algo_lower == "csr"): + "This is the only method implemented for now ..." + # Prepare LUT if needed! + # initialize the CSR integrator in Cython as it may be needed later on. + cython_method = IntegrationMethod.select_method(method.dimension, method.split_lower, method.algo_lower, "cython")[0] + if cython_method not in self.engines: + cython_engine = self.engines[cython_method] = Engine() + else: + cython_engine = self.engines[cython_method] + with cython_engine.lock: + # Validate that the engine used is the proper one + cython_integr = cython_engine.engine + cython_reset = None + if cython_integr is None: + cython_reset = "of first initialization" + if (not cython_reset) and safe: + if cython_integr.unit != unit: + cython_reset = "unit was changed" + if cython_integr.bins != npt: + cython_reset = "number of points changed" + if cython_integr.size != data.size: + cython_reset = "input image size changed" + if cython_integr.empty != self._empty: + cython_reset = "empty value changed " + if (mask is not None) and (not cython_integr.check_mask): + cython_reset = "mask but CSR was without mask" + elif (mask is None) and (cython_integr.check_mask): + cython_reset = "no mask but CSR has mask" + elif (mask is not None) and (cython_integr.mask_checksum != mask_crc): + cython_reset = "mask changed" + if (radial_range is None) and (cython_integr.pos0_range is not None): + cython_reset = "radial_range was defined in CSR" + elif (radial_range is not None) and cython_integr.pos0_range != (min(radial_range), max(radial_range)): + cython_reset = "radial_range is defined but not the same as in CSR" + if (azimuth_range is None) and (cython_integr.pos1_range is not None): + cython_reset = "azimuth_range not defined and CSR had azimuth_range defined" + elif (azimuth_range is not None) and cython_integr.pos1_range != (min(azimuth_range), max(azimuth_range)): + cython_reset = "azimuth_range requested and CSR's azimuth_range don't match" + if cython_reset: + logger.info("AI.sigma_clip_ng: Resetting Cython integrator because %s", cython_reset) + split = method.split_lower + if split == "pseudo": + split = "full" + try: + cython_integr = self.setup_sparse_integrator(data.shape, npt, mask=mask, + mask_checksum=mask_crc, + unit=unit, split=split, algo="CSR", + pos0_range=radial_range, + pos1_range=azimuth_range, + empty=self._empty, + scale=False) + except MemoryError: # CSR method is hungry... + logger.warning("MemoryError: falling back on forward implementation") + cython_integr = None + self.reset_engines() + method = self.DEFAULT_METHOD_1D + else: + cython_engine.set_engine(cython_integr) + if method not in self.engines: + # instanciated the engine + engine = self.engines[method] = Engine() + else: + engine = self.engines[method] + with engine.lock: + # Validate that the engine used is the proper one + integr = engine.engine + reset = None + # This whole block uses CSR, Now we should treat all the various implementation: Cython, OpenCL and finally Python. + + # Validate that the engine used is the proper one + if integr is None: + reset = "of first initialization" + if (not reset) and safe: + if integr.unit != unit: + reset = "unit was changed" + if integr.bins != npt: + reset = "number of points changed" + if integr.size != data.size: + reset = "input image size changed" + if integr.empty != self._empty: + reset = "empty value changed " + if (mask is not None) and (not integr.check_mask): + reset = "mask but CSR was without mask" + elif (mask is None) and (integr.check_mask): + reset = "no mask but CSR has mask" + elif (mask is not None) and (integr.mask_checksum != mask_crc): + reset = "mask changed" + # TODO + if (radial_range is None) and (integr.pos0_range is not None): + reset = "radial_range was defined in CSR" + elif (radial_range is not None) and integr.pos0_range != (min(radial_range), max(radial_range)): + reset = "radial_range is defined but not the same as in CSR" + if (azimuth_range is None) and (integr.pos1_range is not None): + reset = "azimuth_range not defined and CSR had azimuth_range defined" + elif (azimuth_range is not None) and integr.pos1_range != (min(azimuth_range), max(azimuth_range)): + reset = "azimuth_range requested and CSR's azimuth_range don't match" + + if reset: + logger.info("ai.sigma_clip_ng: Resetting ocl_csr integrator because %s", reset) + csr_integr = self.engines[cython_method].engine + if method.impl_lower == "opencl": + try: + integr = method.class_funct_ng.klass(csr_integr.lut, + image_size=data.size, + checksum=csr_integr.lut_checksum, + empty=self._empty, + unit=unit, + mask_checksum=csr_integr.mask_checksum, + bin_centers=csr_integr.bin_centers, + platformid=method.target[0], + deviceid=method.target[1]) + except MemoryError: + logger.warning("MemoryError: falling back on default forward implementation") + self.reset_engines() + method = self.DEFAULT_METHOD_1D + else: + # Copy some properties from the cython integrator + integr.pos0_range = csr_integr.pos0_range + integr.pos1_range = csr_integr.pos1_range + engine.set_engine(integr) + elif method.impl_lower in ("python", "cython"): + integr = method.class_funct_ng.klass(lut=csr_integr.lut, + image_size=data.size, + empty=self._empty, + unit=unit, + mask_checksum=csr_integr.mask_checksum, + bin_centers=csr_integr.bin_centers) + # Copy some properties from the cython integrator + integr.pos0_range = csr_integr.pos0_range + integr.pos1_range = csr_integr.pos1_range + engine.set_engine(integr) + else: + logger.error(f"Implementation {method.impl_lower} not supported") + else: + integr = self.engines[method].engine + kwargs = {"dark":dark, "dummy":dummy, "delta_dummy":delta_dummy, + "variance":variance, "dark_variance":None, + "flat":flat, "solidangle":solidangle, "polarization":polarization, "absorption":absorption, + "error_model":error_model, "normalization_factor":normalization_factor, + "quant_min":quant_min, "quant_max":quant_max} + + intpl = integr.medfilt(data, **kwargs) + else: + raise RuntimeError("Not yet implemented. Sorry") + result = Integrate1dResult(intpl.position * unit.scale, intpl.intensity, intpl.sem) + result._set_method_called("sigma_clip_ng") + result._set_method(method) + result._set_compute_engine(str(method)) + result._set_percentile(percentile) + result._set_unit(unit) + result._set_has_mask_applied(has_mask) + result._set_has_dark_correction(has_dark) + result._set_has_flat_correction(has_flat) + result._set_metadata(metadata) + result._set_sum_signal(intpl.signal) + result._set_sum_normalization(intpl.normalization) + result._set_sum_normalization2(intpl.norm_sq) + result._set_std(intpl.std) + result._set_sem(intpl.sem) + result._set_sum_variance(intpl.variance) + result._set_count(intpl.count) + result._set_polarization_factor(polarization_factor) + result._set_normalization_factor(normalization_factor) + result._set_error_model(error_model) + return result + def sigma_clip_legacy(self, data, npt_rad=1024, npt_azim=512, correctSolidAngle=True, polarization_factor=None, radial_range=None, azimuth_range=None, diff --git a/version.py b/version.py index 324cf9028..dbeb16fc1 100755 --- a/version.py +++ b/version.py @@ -47,7 +47,7 @@ __authors__ = ["Jérôme Kieffer", "V. Valls"] __license__ = "MIT" __copyright__ = "European Synchrotron Radiation Facility, Grenoble, France" -__date__ = "10/10/2024" +__date__ = "05/12/2024" __status__ = "production" __docformat__ = 'restructuredtext' __all__ = ["date", "version_info", "strictversion", "hexversion", "debianversion", @@ -61,7 +61,7 @@ "final": 15} MAJOR = 2024 -MINOR = 11 +MINOR = 12 MICRO = 0 RELEV = "dev" # <16 SERIAL = 0 # <16