From 4eb78216789fb125f983ca57ddaa11a044b44da4 Mon Sep 17 00:00:00 2001
From: David Hammer <dhammer@mailbox.org>
Date: Tue, 27 Jun 2023 15:59:21 +0200
Subject: [PATCH] Use faster common mode (row / col) for ePix100, too

---
 src/calng/corrections/Epix100Correction.py | 19 ++++++++-----------
 src/calng/corrections/PnccdCorrection.py   |  5 +++--
 src/calng/kernels/common_cpu.pyx           | 12 ++++++------
 src/tests/test_pnccd_kernels.py            |  1 -
 4 files changed, 17 insertions(+), 20 deletions(-)

diff --git a/src/calng/corrections/Epix100Correction.py b/src/calng/corrections/Epix100Correction.py
index 98411ecd..9c7df5c2 100644
--- a/src/calng/corrections/Epix100Correction.py
+++ b/src/calng/corrections/Epix100Correction.py
@@ -21,6 +21,7 @@ from .. import (
     utils,
 )
 from .._version import version as deviceVersion
+from ..kernels import common_cython
 
 
 class Constants(enum.Enum):
@@ -251,11 +252,11 @@ class Epix100CpuRunner(base_kernel_runner.BaseKernelRunner):
 
         if flags & CorrectionFlags.COMMONMODE:
             # per rectangular block that looks like something is going on
-            masked = np.ma.masked_array(
-                data=output,
-                mask=(self._q_bad_pixel_map[q] != 0)
-                | (output > self._q_noise_map[q] * cm_noise_sigma),
-            )
+            cm_mask=(
+                (self._q_bad_pixel_map[q] != 0)
+                | (output > self._q_noise_map[q] * cm_noise_sigma)
+            ).astype(np.uint8, copy=False)
+            masked = np.ma.masked_array(data=output, mask=cm_mask)
             if cm_block:
                 for block in np.hsplit(masked, 4):
                     if block.count() < block.size * cm_min_frac:
@@ -263,14 +264,10 @@ class Epix100CpuRunner(base_kernel_runner.BaseKernelRunner):
                     block.data[:] -= np.ma.median(block)
 
             if cm_row:
-                subset_rows = masked.count(axis=1) >= masked.shape[1] * cm_min_frac
-                output[subset_rows] -= np.ma.median(
-                    masked[subset_rows], axis=1, keepdims=True
-                )
+                common_cython.cm_fs(output, cm_mask, cm_noise_sigma, cm_min_frac)
 
             if cm_col:
-                subset_cols = masked.count(axis=0) >= masked.shape[0] * cm_min_frac
-                output[:, subset_cols] -= np.ma.median(masked[:, subset_cols], axis=0)
+                common_cython.cm_ss(output, cm_mask, cm_noise_sigma, cm_min_frac)
 
         if flags & CorrectionFlags.RELGAIN:
             output *= self._q_rel_gain_map[q]
diff --git a/src/calng/corrections/PnccdCorrection.py b/src/calng/corrections/PnccdCorrection.py
index 36b8ee15..4a8c1491 100644
--- a/src/calng/corrections/PnccdCorrection.py
+++ b/src/calng/corrections/PnccdCorrection.py
@@ -259,8 +259,9 @@ class PnccdCpuRunner(base_kernel_runner.BaseKernelRunner):
             output -= self._q_offset_map[q]
 
         if flags & CorrectionFlags.COMMONMODE:
-            cm_mask = (self._q_bad_pixel_map[q] != 0) | (
-                output > self._q_noise_map[q] * cm_noise_sigma
+            cm_mask = (
+                (self._q_bad_pixel_map[q] != 0)
+                | (output > self._q_noise_map[q] * cm_noise_sigma)
             ).astype(np.uint8, copy=False)
             if cm_row:
                 common_cython.cm_fs(
diff --git a/src/calng/kernels/common_cpu.pyx b/src/calng/kernels/common_cpu.pyx
index fc61d6d0..893015f2 100644
--- a/src/calng/kernels/common_cpu.pyx
+++ b/src/calng/kernels/common_cpu.pyx
@@ -48,7 +48,7 @@ def cm_ss(
     float ratio,
 ):
     # transposing to give std::nth_element contiguous input
-    cdef float[:, ::1] a = np.empty_like(np.transpose(data), order="C")
+    cdef float[:, ::1] aux = np.empty_like(np.transpose(data), order="C")
     cdef int num_dark, min_dark, i, j, k
     cdef float median
     min_dark = <int>(<float>data.shape[1] * ratio)
@@ -56,15 +56,15 @@ def cm_ss(
         num_dark = 0
         for i in range(data.shape[0]):
             if mask[i, j] == 0:
-                a[j, num_dark] = data[i, j]
+                aux[j, num_dark] = data[i, j]
                 num_dark = num_dark + 1
         if num_dark < min_dark:
             continue
         k = num_dark // 2
-        nth_element(&a[j, 0], (&a[j, 0]) + k, (&a[j, 0]) + <int>num_dark)
-        median = a[j, k]
+        nth_element(&aux[j, 0], (&aux[j, 0]) + k, (&aux[j, 0]) + <int>num_dark)
+        median = aux[j, k]
         if num_dark % 2 == 0:
-            nth_element(&a[j, 0], (&a[j, 0]) + k - 1, (&a[j, 0]) + <int>num_dark)
-            median = median / 2 + a[j, k - 1] / 2
+            nth_element(&aux[j, 0], (&aux[j, 0]) + k - 1, (&aux[j, 0]) + <int>num_dark)
+            median = median / 2 + aux[j, k - 1] / 2
         for i in range(data.shape[0]):
             data[i, j] = data[i, j] - median
diff --git a/src/tests/test_pnccd_kernels.py b/src/tests/test_pnccd_kernels.py
index 470eb025..18f438fa 100644
--- a/src/tests/test_pnccd_kernels.py
+++ b/src/tests/test_pnccd_kernels.py
@@ -20,7 +20,6 @@ def test_common_mode():
             subset_ss = masked.count(axis=0) >= masked.shape[0] * cm_min_frac
             data[:, subset_ss] -= np.ma.median(masked[:, subset_ss], axis=0)
 
-
     data_shape = (1024, 1024)
     data = (np.random.random(data_shape) * 1000).astype(np.uint16)
     runner = PnccdCpuRunner(1024, 1024, 1, 1)
-- 
GitLab