Merge pull request #558 from kif/557_sparse

Support for sparse frames without background (XPCS among other)

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__ = "[email protected]"
 __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 = <double>(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 = <uint32_t> fpos
-                    if pos+1 == size:
-                        value = background[pos]
-                        fres = 0.0
+        if do_background:
+            start = radius[0]
+            idelta = <double>(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 = <uint32_t> 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] =  <any_t>(value + 0.5) #this is rounding
-                    else:
-                        dense[i,j] =  <any_t>(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] =  <any_t>(value + 0.5) #this is rounding
+                        else:
+                            dense[i,j] =  <any_t>(value) 
         # Assignment of outliers
         for i in range(size_over):
             j = index[i]

src/fabio/sparseimage.py

@@ -37,7 +37,7 @@
 __contact__ = "[email protected]"
 __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,27 +213,48 @@ 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],
                            self.dummy,
                            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: