From ceb6a94d6b93eacf4ec8c1664efa30ae3b277399 Mon Sep 17 00:00:00 2001 From: Jerome Kieffer Date: Tue, 19 Mar 2024 10:10:39 +0100 Subject: [PATCH] Support for sparse frames without background (XPCS) --- src/fabio/ext/dense.pyx | 66 ++++++++++++++++++++------------------- src/fabio/sparseimage.py | 67 ++++++++++++++++++++++++++++------------ 2 files changed, 83 insertions(+), 50 deletions(-) diff --git a/src/fabio/ext/dense.pyx b/src/fabio/ext/dense.pyx index 8ee39ba9d..fce7b9d42 100644 --- a/src/fabio/ext/dense.pyx +++ b/src/fabio/ext/dense.pyx @@ -8,7 +8,7 @@ # Project: Fable Input/Output # https://github.com/silx-kit/fabio # -# Copyright (C) 2020-2023 European Synchrotron Radiation Facility, Grenoble, France +# Copyright (C) 2020-2024 European Synchrotron Radiation Facility, Grenoble, France # # Permission is hereby granted, free of charge, to any person obtaining a copy # of this software and associated documentation files (the "Software"), to deal @@ -31,7 +31,7 @@ """Densification of sparse frame format """ __author__ = "Jérôme Kieffer" -__date__ = "20/10/2023" +__date__ = "19/03/2024" __contact__ = "Jerome.kieffer@esrf.fr" __license__ = "MIT" @@ -249,12 +249,15 @@ def densify(cython.floating[:,::1] mask, cdef: Py_ssize_t i, j, size, pos, size_over, width, height double value, fres, fpos, idelta, start, std - bint integral, noisy, do_normalization=False + bint integral, noisy, do_normalization=False, do_background=True any_t[:, ::1] dense float[:,::1] c_normalization MT mt - size = radius.shape[0] - assert background.shape[0] == size + if radius is None: + do_background = False + else: + size = radius.shape[0] + assert background.shape[0] == size size_over = index.shape[0] assert intensity.shape[0] == size_over integral = numpy.issubdtype(dtype, numpy.integer) @@ -277,36 +280,37 @@ def densify(cython.floating[:,::1] mask, mt = MT(seed) with nogil: - start = radius[0] - idelta = (size - 1)/(radius[size-1] - start) - - #Linear interpolation - for i in range(height): - for j in range(width): - fpos = (mask[i,j] - start)*idelta - if (fpos<0) or (fpos>=size) or (not isfinite(fpos)): - dense[i,j] = dummy - else: - pos = fpos - if pos+1 == size: - value = background[pos] - fres = 0.0 + if do_background: + start = radius[0] + idelta = (size - 1)/(radius[size-1] - start) + + #Linear interpolation + for i in range(height): + for j in range(width): + fpos = (mask[i,j] - start)*idelta + if (fpos<0) or (fpos>=size) or (not isfinite(fpos)): + dense[i,j] = dummy else: - fres = fpos - pos - value = (1.0 - fres)*background[pos] + fres*background[pos+1] - if noisy: + pos = fpos if pos+1 == size: - std = background_std[pos] + value = background[pos] fres = 0.0 else: - std = (1.0 - fres)*background_std[pos] + fres*background_std[pos+1] - value = max(0.0, mt._normal_m(value, std)) - if do_normalization: - value *= c_normalization[i, j] - if integral: - dense[i,j] = (value + 0.5) #this is rounding - else: - dense[i,j] = (value) + fres = fpos - pos + value = (1.0 - fres)*background[pos] + fres*background[pos+1] + if noisy: + if pos+1 == size: + std = background_std[pos] + fres = 0.0 + else: + std = (1.0 - fres)*background_std[pos] + fres*background_std[pos+1] + value = max(0.0, mt._normal_m(value, std)) + if do_normalization: + value *= c_normalization[i, j] + if integral: + dense[i,j] = (value + 0.5) #this is rounding + else: + dense[i,j] = (value) # Assignment of outliers for i in range(size_over): j = index[i] diff --git a/src/fabio/sparseimage.py b/src/fabio/sparseimage.py index 30a7eaddb..c03d0ac82 100644 --- a/src/fabio/sparseimage.py +++ b/src/fabio/sparseimage.py @@ -37,7 +37,7 @@ __contact__ = "jerome.kieffer@esrf.fr" __license__ = "MIT" __copyright__ = "2020 ESRF" -__date__ = "15/03/2024" +__date__ = "19/03/2024" import logging logger = logging.getLogger(__name__) @@ -82,7 +82,11 @@ def densify(mask, :param seed: numerical seed for random number generator :return dense array """ - dense = numpy.interp(mask, radius, background) + dtype = intensity.dtype + if radius is None or background is None: + dense = numpy.zeros(radius.shape, dtype=dtype) + else: + dense = numpy.interp(mask, radius, background) if background_std is not None: if seed is not None: numpy.random.seed(seed) @@ -93,7 +97,7 @@ def densify(mask, flat = dense.ravel() flat[index] = intensity - dtype = intensity.dtype + if numpy.issubdtype(dtype, numpy.integer): # Foolded by banker's rounding !!!! dense +=0.5 @@ -173,9 +177,13 @@ def read(self, fname, frame=None): raise NotGoodReader("HDF5 file does not contain any default NXdata.") nx_data = entry[default_data] self.mask = nx_data["mask"][()] - self.radius = nx_data["radius"][()] - self.background_avg = nx_data["background_avg"][()] - self.background_std = nx_data["background_std"][()] + try: + self.radius = nx_data["radius"][()] + self.background_avg = nx_data["background_avg"][()] + self.background_std = nx_data["background_std"][()] + except KeyError: + logger.info("No background information found") + self.radius = self.background_avg = self.background_std = None self.frame_ptr = nx_data["frame_ptr"][()] self.index = nx_data["index"][()] self.intensity = nx_data["intensity"][()] @@ -205,19 +213,29 @@ def _generate_data(self, index=0): logger.warning("Not data have been read from disk") return start, stop = self.frame_ptr[index:index + 2] - if cython_densify is not None: - return cython_densify.densify(self.mask, - self.radius, - self.index[start:stop], - self.intensity[start:stop], - self.dummy, - self.intensity.dtype, - self.background_avg[index], - self.background_std[index] * self.noisy if self.noisy else None, - self.normalization) + if self.radius is None: + if cython_densify is None: # Numpy implementation + dense = densify(self.mask, + None, + self.index[start:stop], + self.intensity[start:stop], + self.dummy, + None, + None, + self.normalization) + else: # Cython + dense = cython_densify.densify(self.mask, + None, + self.index[start:stop], + self.intensity[start:stop], + self.dummy, + self.intensity.dtype, + None, + None, + self.normalization) else: - # Fall-back on numpy code. - return densify(self.mask, + if cython_densify is None: # Numpy + dense = densify(self.mask, self.radius, self.index[start:stop], self.intensity[start:stop], @@ -225,7 +243,18 @@ def _generate_data(self, index=0): self.background_avg[index], self.background_std[index] * self.noisy if self.noisy else None, self.normalization) - + else: + dense = cython_densify.densify(self.mask, + self.radius, + self.index[start:stop], + self.intensity[start:stop], + self.dummy, + self.intensity.dtype, + self.background_avg[index], + self.background_std[index] * self.noisy if self.noisy else None, + self.normalization) + return dense + def getframe(self, num): """ returns the frame numbered 'num' in the stack if applicable""" if self.nframes > 1: