Add new centroid algorithm by Martin

src/toolbox_scs/detectors/hrixs.py

@@ -156,7 +156,7 @@ CURVE_A = 2.19042931e-02  # curvature parameters as determined elsewhere
 CURVE_B = -3.02191568e-07
 
 
-def centroid(image, threshold=THRESHOLD, curvature=(CURVE_A, CURVE_B)):
+def _esrf_centroid(image, threshold=THRESHOLD, curvature=(CURVE_A, CURVE_B)):
     gs = 2
     base = image.mean()
     cp = np.argwhere(image[gs // 2: -gs // 2, gs // 2: -gs // 2] > threshold) + np.array([gs // 2, gs // 2])
@@ -175,6 +175,43 @@ def centroid(image, threshold=THRESHOLD, curvature=(CURVE_A, CURVE_B)):
     return res
 
 
+def _new_centroid(image, threshold=THRESHOLD, curvature=(CURVE_A, CURVE_B)):
+    """find the position of photons with sub-pixel precision
+
+    A photon is supposed to have hit the detector if the intensity within a
+    2-by-2 square exceeds a threshold. In this case the position of the photon
+    is calculated as the center-of-mass in a 4-by-4 square.
+
+    Return the list of x,y coordinate pairs, corrected by the curvature.
+    """
+    base = image.mean()
+    corners = image[1:, 1:] + image[:-1, 1:] + image[1:, :-1] + image[:-1, :-1]
+    threshold = corners.mean() + 3.5 * corners.std()
+    middle = corners[1:-1, 1:-1]
+    candidates = (
+            (middle > threshold)
+            * (middle >= corners[:-2, 1:-1]) * (middle > corners[2:, 1:-1])
+            * (middle >= corners[1:-1, :-2]) * (middle > corners[1:-1, 2:])
+            * (middle >= corners[:-2, :-2]) * (middle > corners[2:, :-2])
+            * (middle >= corners[:-2, 2:]) * (middle > corners[2:, 2:]))
+    cp = np.argwhere(candidates)
+    if len(cp) > 10000:
+        raise RuntimeError(
+            "too many peaks, threshold too low or acquisition time too high")
+
+    res = []
+    for cy, cx in cp:
+        spot = image[cy: cy + 4, cx: cx + 4] - base
+        mx = np.average(np.arange(cx, cx + 4), weights=spot.sum(axis=0))
+        my = np.average(np.arange(cy, cy + 4), weights=spot.sum(axis=1))
+        my -= (curvature[0] + curvature[1] * mx) * mx
+        res.append((my, mx))
+    return res
+
+
+centroid = _new_centroid
+
+
 def decentroid(res):
     res = np.array(res)
     ret = np.zeros(shape=(res.max(axis=0) + 1).astype(int))