Remains at least a bug since one fragment of a pixel is missing

src/pyFAI/opencl/azim_csr.py

@@ -179,6 +179,7 @@ def __init__(self, lut, image_size, checksum=None,
         self.buffer_dtype = {i.name:numpy.dtype(i.dtype) for i in self.buffers}
         if numpy.allclose(self._data, numpy.ones(self._data.shape)):
             self.cl_mem["data"] = None
+            self.cl_kernel_args["csr_medfilt"]["data"] = None
             self.cl_kernel_args["csr_sigma_clip4"]["data"] = None
             self.cl_kernel_args["csr_integrate"]["data"] = None
             self.cl_kernel_args["csr_integrate4"]["data"] = None
@@ -1263,6 +1264,7 @@ def medfilt(self, data, dark=None, dummy=None, delta_dummy=None,
             ev = pyopencl.enqueue_copy(self.queue, merged, self.cl_mem["merged8"])
             events.append(EventDescription("copy D->H merged8", ev))
         self.profile_multi(events)
+
         res = Integrate1dtpl(self.bin_centers, avgint, sem, merged[:, 0], merged[:, 2], merged[:, 4], merged[:, 6],
                              std, sem, merged[:, 7])
         return res

src/pyFAI/opencl/test/test_ocl_azim_csr.py

@@ -33,7 +33,7 @@
 __contact__ = "[email protected]"
 __license__ = "MIT"
 __copyright__ = "2019-2021 European Synchrotron Radiation Facility, Grenoble, France"
-__date__ = "28/06/2022"
+__date__ = "05/12/2024"
 
 import logging
 import numpy
@@ -43,10 +43,9 @@
 if ocl:
     import pyopencl.array
 from ...test.utilstest import UtilsTest
-from silx.opencl.common import _measure_workgroup_size
 from ...integrator.azimuthal import AzimuthalIntegrator
 from ...method_registry import IntegrationMethod
-from scipy.ndimage import gaussian_filter1d
+from ...containers import ErrorModel
 logger = logging.getLogger(__name__)
 
 
@@ -81,10 +80,12 @@ def tearDownClass(cls):
         cls.ai = None
 
     @unittest.skipUnless(ocl, "pyopencl is missing")
-    def integrate_ng(self, block_size=None):
+    def integrate_ng(self, block_size=None, method_called="integrate_ng", extra=None):
         """
         tests the 1d histogram kernel, with variable workgroup size
         """
+        if extra is None:
+            extra={}
         from ..azim_csr import OCL_CSR_Integrator
         data = numpy.ones(self.ai.detector.shape)
         npt = 500
@@ -95,14 +96,14 @@ def integrate_ng(self, block_size=None):
                                                             dim=1, default=None, degradable=False)
 
         # Retrieve the CSR array
-        cpu_integrate = self.ai._integrate1d_legacy(data, npt, unit=unit, method=csr_method)
+        cpu_integrate = self.ai._integrate1d_ng(data, npt, unit=unit, method=csr_method)
         r_m = cpu_integrate[0]
         csr_engine = list(self.ai.engines.values())[0]
         csr = csr_engine.engine.lut
-        ref = self.ai._integrate1d_ng(data, npt, unit=unit, method=method)
-        integrator = OCL_CSR_Integrator(csr, data.size, block_size=block_size)
+        ref = self.ai._integrate1d_ng(data, npt, unit=unit, method=method, error_model="poisson")
+        integrator = OCL_CSR_Integrator(csr, data.size, block_size=block_size, empty=-1)
         solidangle = self.ai.solidAngleArray()
-        res = integrator.integrate_ng(data, solidangle=solidangle)
+        res = integrator.__getattribute__(method_called)(data, solidangle=solidangle, error_model=ErrorModel.POISSON, **extra)
         # for info, res contains: position intensity error signal variance normalization count
 
         # Start with smth easy: the position
@@ -112,23 +113,25 @@ def integrate_ng(self, block_size=None):
         if "AMD" in integrator.ctx.devices[0].platform.name:
             logger.warning("This test is known to be complicated for AMD-GPU, relax the constrains for them")
         else:
-            self.assertLessEqual(delta.max(), 1, "counts are almost the same")
-            self.assertEqual(delta.sum(), 0, "as much + and -")
+            if method_called=="integrate_ng":
+                self.assertLessEqual(delta.max(), 1, "counts are almost the same")
+                self.assertEqual(delta.sum(), 0, "as much + and -")
+            elif method_called=="medfilt":
+                pix = csr[2][1:]-csr[2][:-1]
+                for i,j in zip(pix,res.count):
+                    print(j-i,i,j)
 
             # Intensities are not that different:
             delta = ref.intensity - res.intensity
-            self.assertLessEqual(abs(delta.max()), 1e-5, "intensity is almost the same")
+            self.assertLessEqual(abs(delta).max(), 1e-5, "intensity is almost the same")
 
             # histogram of normalization
-            ref = self.ai._integrate1d_ng(solidangle, npt, unit=unit, method=method).sum_signal
-            sig = res.normalization
-            err = abs((sig - ref).max())
+            err = abs((res.normalization - ref.sum_normalization).max())
             self.assertLess(err, 5e-4, "normalization content is the same: %s<5e-5" % (err))
 
             # histogram of signal
-            ref = self.ai._integrate1d_ng(data, npt, unit=unit, method=method).sum_signal
-            sig = res.signal
-            self.assertLess(abs((sig - ref).sum()), 5e-5, "signal content is the same")
+            print(res.signal - ref.sum_signal)
+            self.assertLess(abs((res.signal - ref.sum_signal).sum()), 5e-5, "signal content is the same")
 
     @unittest.skipUnless(ocl, "pyopencl is missing")
     def test_integrate_ng(self):
@@ -144,6 +147,20 @@ def test_integrate_ng_single(self):
         """
         self.integrate_ng(block_size=1)
 
+    @unittest.skipUnless(ocl, "pyopencl is missing")
+    def test_sigma_clip(self):
+        """
+        tests the sigma-clipping kernel, default block size
+        """
+        self.integrate_ng(None, "sigma_clip",{"cutoff":100.0, "cycle":0,})
+
+    @unittest.skipUnless(ocl, "pyopencl is missing")
+    def test_medfilt(self):
+        """
+        tests the median filtering kernel, default block size
+        """
+        self.integrate_ng(None, "medfilt", {"quant_min":0, "quant_max":1})
+
 
 def suite():
     loader = unittest.defaultTestLoader.loadTestsFromTestCase

src/pyFAI/resources/openCL/medfilt.cl

@@ -129,13 +129,15 @@ csr_medfilt    (  const   global  float4  *data4,
     float sum=0.0f, ratio=1.3f;
     float2 acc_sig, acc_nrm, acc_var, acc_nrm2;
 
+//    if (tid==0)printf("bin %d start %d stop %d\n",bin_num,start, stop);
     // first populate the work4 array from data4
     for (int i=start+tid; i<stop; i+=wg)
     {
         float4 r4, w4;
         int idx = indices[i];
-        float coef = coefs[i];
+        float coef = (coefs == ZERO)?1.0f:coefs[i];
         r4 = data4[idx];
+//        if (tid==0)printf("bin %d start %d stop %d i=%d idx=%d coef=%6.4f\n",bin_num,start, stop,i,idx,coef);
 
         w4.s0 = r4.s0 / r4.s2;
         w4.s1 = r4.s0 * coef;
@@ -203,8 +205,9 @@ csr_medfilt    (  const   global  float4  *data4,
         }
     }
 //  Now parallel reductions, one after the other :-/
+//    if (tid==0)printf("bin %d cnt %d acc_sig %6.1f acc_var %6.1f acc_nrm %6.1f acc_nrm2 %6.1f\n",bin_num, cnt, acc_sig.s0, acc_var.s0, acc_nrm.s0, acc_nrm2.s0);
     shared_int[tid] = cnt;
-    cnt += sum_int_reduction(shared_int);
+    cnt = sum_int_reduction(shared_int);
 
     shared_float[2*tid] = acc_sig.s0;
     shared_float[2*tid+1] = acc_sig.s1;
@@ -225,7 +228,7 @@ csr_medfilt    (  const   global  float4  *data4,
 
     // Finally store the accumulated value
 
-    if (get_local_id(0) == 0) {
+    if (tid == 0) {
         summed[bin_num] = (float8)(acc_sig.s0, acc_sig.s1,
                                    acc_var.s0, acc_var.s1,
                                    acc_nrm.s0, acc_nrm.s1,