Handle the case "median filter"

Tests remain noisy
silx-kit · Dec 5, 2024 · 7cfa9ff · 7cfa9ff
1 parent 9a3ff3d
commit 7cfa9ff
Show file tree

Hide file tree

Showing 4 changed files with 15 additions and 6 deletions.
diff --git a/src/pyFAI/ext/CSR_common.pxi b/src/pyFAI/ext/CSR_common.pxi
@@ -29,7 +29,7 @@
 
 __author__ = "Jérôme Kieffer"
 __contact__ = "[email protected]"
-__date__ = "19/11/2024"
+__date__ = "05/12/2024"
 __status__ = "stable"
 __license__ = "MIT"
 
@@ -859,7 +859,8 @@ cdef class CsrIntegrator(object):
                 for i in range(start, stop):
                     former_element = element
                     element = work[i]
-                    if (qmin<=former_element.s0) and (element.s0 <= qmax):
+                    if ((qmin<=former_element.s0) and (element.s0 <= qmax)) or \
+                       ((qmin>=former_element.s0)  and (element.s0 >= qmax)):   #specific case where qmin==qmax
                         acc_sig = acc_sig + element.s1
                         acc_var = acc_var + element.s2
                         acc_norm = acc_norm + element.s3

diff --git a/src/pyFAI/integrator/azimuthal.py b/src/pyFAI/integrator/azimuthal.py
@@ -1148,7 +1148,7 @@ def integrate2d_ng(self, data, npt_rad, npt_azim=360,
 
     integrate2d = _integrate2d_ng = integrate2d_ng
 
-    @deprecated(since_version="2024.12", only_once=True, replacement="medfilt1d_ng", deprecated_since="2024.12")
+    @deprecated(since_version="2024.12.0", only_once=True, replacement="medfilt1d_ng", deprecated_since="2024.12.0")
     def medfilt1d_legacy(self, data, npt_rad=1024, npt_azim=512,
                   correctSolidAngle=True,
                   radial_range=None, azimuth_range=None,

diff --git a/src/pyFAI/opencl/azim_csr.py b/src/pyFAI/opencl/azim_csr.py
@@ -1221,7 +1221,7 @@ def medfilt(self, data, dark=None, dummy=None, delta_dummy=None,
             events.append(EventDescription(kernel_correction_name, ev))
 
             kw_int["quant_min"] = numpy.float32(quant_min)
-            kw_int["quant_max"] = numpy.float32(quant_max)*EPS32
+            kw_int["quant_max"] = numpy.float32(quant_max)#*EPS32
 
             wg_min = max(self.workgroup_size["csr_medfilt"])
             kw_int["shared_int"] = pyopencl.LocalMemory(4 * wg_min)

diff --git a/src/pyFAI/resources/openCL/medfilt.cl b/src/pyFAI/resources/openCL/medfilt.cl
@@ -105,7 +105,7 @@ csr_medfilt    (  const   global  float4  *data4,
                   const   global  int     *indices,
                   const   global  int     *indptr,
                   const           float    quant_min,
-                  const           float    quant_max,
+                                  float    quant_max,
                   const           char     error_model,
                   const           float    empty,
                           global  float8  *summed,
@@ -129,6 +129,10 @@ csr_medfilt    (  const   global  float4  *data4,
     float sum=0.0f, ratio=1.3f;
     float2 acc_sig, acc_nrm, acc_var, acc_nrm2;
 
+    // ensure the last element is always taken
+    if (quant_max == 1.0f)
+        quant_max = 1.1f;
+
     // first populate the work4 array from data4
     for (int i=start+tid; i<stop; i+=wg)
     {
@@ -160,6 +164,7 @@ csr_medfilt    (  const   global  float4  *data4,
     // Then perform the cumsort of the weights to s0
     // In blelloch scan, one workgroup can process 2wg in size.
 
+    barrier(CLK_GLOBAL_MEM_FENCE);
 
     sum = 0.0f;
     for (int i=0; i<(size + 2*wg-1)/(2*wg); i++)
@@ -179,6 +184,8 @@ csr_medfilt    (  const   global  float4  *data4,
 
     }
 
+    barrier(CLK_GLOBAL_MEM_FENCE);
+
     // Perform the sum for accumulator of signal, variance, normalization and count
 
     cnt = 0;
@@ -194,7 +201,8 @@ csr_medfilt    (  const   global  float4  *data4,
     {
         float q_last = (i>start)?work4[i-1].s0:0.0f;
         float4 w = work4[i];
-        if ((q_last>=qmin) && (w.s0<=qmax)){
+        if (((q_last>=qmin) && (w.s0<=qmax))||
+            ((q_last<=qmin) && (w.s0>=qmax)))   { // case qmin==qmax
             cnt ++;
             acc_sig = dw_plus_fp(acc_sig, w.s1);
             acc_var = dw_plus_fp(acc_var, w.s2);