Merge pull request #181 from fmi-faim/2d-spot-detection

Add detect_spots_2d.

src/faim_ipa/detection/blobs.py

@@ -56,7 +56,7 @@ def detect_blobs(
     )
 
     rescale_factor = estimate_log_rescale_factor(
-        axial_sigma=axial_sigma, lateral_sigma=lateral_sigma
+        (axial_sigma, lateral_sigma, lateral_sigma),
     )
 
     sigmas = [

src/faim_ipa/detection/spots.py

@@ -2,10 +2,12 @@
 
 import numpy as np
 from scipy.ndimage import gaussian_laplace
-from skimage.morphology import ball, h_maxima
+from skimage.morphology import ball, h_maxima, disk
 from skimage.util import img_as_float32
 
-from faim_ipa.detection.utils import estimate_log_rescale_factor
+from faim_ipa.detection.utils import (
+    estimate_log_rescale_factor,
+)
 
 
 def detect_spots(
@@ -46,7 +48,7 @@ def detect_spots(
     )
 
     rescale_factor = estimate_log_rescale_factor(
-        axial_sigma=axial_sigma, lateral_sigma=lateral_sigma
+        sigmas=(axial_sigma, lateral_sigma, lateral_sigma)
     )
     log_img = (
         -gaussian_laplace(image, sigma=(axial_sigma, lateral_sigma, lateral_sigma))
@@ -56,3 +58,47 @@ def detect_spots(
     h_ = img_as_float32(np.array(h, dtype=img.dtype))
     h_detections = h_maxima(log_img, h=h_, footprint=ball(1))
     return np.array(np.where(h_detections)).T
+
+
+def detect_spots_2d(
+    img: np.ndarray,
+    lateral_sigma: float,
+    h: int,
+    background_img: np.ndarray | None = None,
+) -> np.ndarray:
+    """Detect diffraction limited spots.
+
+    The spot detection uses a Laplacian of Gaussian filter to detect
+    spots of a given size. These detections are intensity filtered with
+    a h-maxima filter.
+
+    Parameters
+    ----------
+    img :
+        Image containing spot signal.
+    lateral_sigma :
+        YX extension of the spots.
+    h :
+        h-maxima threshold.
+    background_img :
+        Estimated background image. This is subtracted before the
+        spot detection.
+
+    Returns
+    -------
+    Detected spots.
+    """
+    image = (
+        img_as_float32(img) - img_as_float32(background_img)
+        if background_img is not None
+        else img_as_float32(img)
+    )
+
+    rescale_factor = estimate_log_rescale_factor(sigmas=(lateral_sigma, lateral_sigma))
+    log_img = (
+        -gaussian_laplace(image, sigma=(lateral_sigma, lateral_sigma)) * rescale_factor
+    )
+
+    h_ = img_as_float32(np.array(h, dtype=img.dtype))
+    h_detections = h_maxima(log_img, h=h_, footprint=disk(1))
+    return np.array(np.where(h_detections)).T

src/faim_ipa/detection/utils.py

@@ -52,30 +52,28 @@ def compute_lateral_sigma(
     return wavelength / (2 * NA) / (2 * np.sqrt(2 * np.log(2))) / lateral_spacing
 
 
-def estimate_log_rescale_factor(axial_sigma: float, lateral_sigma: float) -> float:
+def estimate_log_rescale_factor(sigmas: tuple[float, ...]) -> float:
     """
     Estimate the rescale factor for a LoG filter response, such that
     the LoG filter response intensities are equal to the input image
-    intensities for spots of size equal to a Gaussian with sigma
+    intensities for spots of size equal to a Gaussian with sigmas
     (axial_sigma, lateral_sigma, lateral_sigma).
 
+    Note: Arbitrary number of sigmas is possible.
+
     Parameters
     ----------
-    axial_sigma :
-        Sigma in axial direction. Usually along Z.
-    lateral_sigma :
-        Sigma in lateral direction. Usually along Y and X.
+    sigmas :
+        Tuple of sigmas used in LoG detection.
 
     Returns
     -------
     rescale_factor
     """
-    extend = int(max(axial_sigma, lateral_sigma) * 7)
-    img = np.zeros((extend,) * 3, dtype=np.float32)
-    img[extend // 2, extend // 2, extend // 2] = 1
-    img = gaussian_filter(img, (axial_sigma, lateral_sigma, lateral_sigma))
+    extend = int(max(sigmas) * 7)
+    img = np.zeros((extend,) * len(sigmas), dtype=np.float32)
+    img[(extend // 2,) * len(sigmas)] = 1
+    img = gaussian_filter(img, sigmas)
     img = img / img.max()
-    img_log = -gaussian_laplace(
-        input=img, sigma=(axial_sigma, lateral_sigma, lateral_sigma)
-    )
+    img_log = -gaussian_laplace(input=img, sigma=sigmas)
     return 1 / img_log.max()

tests/detection/test_detection_utils.py

@@ -25,7 +25,7 @@ def test_compute_lateral_sigma():
 
 
 def test_estimate_log_rescale_factor():
-    rescale_factor = estimate_log_rescale_factor(2.07, 0.75)
+    rescale_factor = estimate_log_rescale_factor((2.07, 0.75, 0.75))
 
     # Create image with diffraction limited spot with intensity = 10
     img = np.zeros((101, 101, 101), dtype=np.float32)
@@ -37,3 +37,18 @@ def test_estimate_log_rescale_factor():
     assert_almost_equal(
         np.max(-gaussian_laplace(img, (2.07, 0.75, 0.75)) * rescale_factor), 10, 2
     )
+
+
+def test_estimate_log_rescale_factor_2d():
+    rescale_factor = estimate_log_rescale_factor((0.75, 0.75))
+
+    # Create image with diffraction limited spot with intensity = 10
+    img = np.zeros((101, 101), dtype=np.float32)
+    img[50, 50] = 1
+    img = gaussian_filter(img, (0.75, 0.75))
+    img = 10 * img / img.max()
+
+    # Max value of the rescaled LoG filter response should be 10
+    assert_almost_equal(
+        np.max(-gaussian_laplace(img, (0.75, 0.75)) * rescale_factor), 10, 2
+    )

tests/detection/test_spots.py

@@ -2,7 +2,7 @@
 from numpy.testing import assert_array_equal
 from scipy.ndimage import gaussian_filter
 
-from faim_ipa.detection.spots import detect_spots
+from faim_ipa.detection.spots import detect_spots, detect_spots_2d
 
 
 def test_detect_spots():
@@ -115,3 +115,111 @@ def test_detect_spots():
             ]
         ),
     )
+
+
+def test_detect_spots_2d():
+    np.random.seed(0)
+    img = np.zeros((101, 101), dtype=np.float32)
+
+    # Create 3 spots with intensities of 100, 200, 300
+    img[25, 25] = 1
+    img[50, 50] = 2
+    img[75, 75] = 3
+    img = gaussian_filter(img, (0.75, 0.75))
+    img = 300 * img / img.max()
+
+    # Create 1 larger spot
+    large_spot = np.zeros((101, 101), dtype=np.float32)
+    large_spot[75, 75] = 1
+    large_spot = gaussian_filter(large_spot, (3 * 0.75, 3 * 0.75))
+    large_spot = 100 * large_spot / large_spot.max()
+
+    # Create random background noise
+    background_img = np.random.normal(10, 2, (101, 101)).astype(np.float32)
+
+    # Create hot-pixel
+    hot_pixels = np.zeros((101, 101), dtype=np.float32)
+    hot_pixels[25, 50] = 500
+
+    # Combine to final image
+    img_final = (img + large_spot + hot_pixels + background_img).astype(np.uint16)
+
+    # Fake estimated background with hot-pixels
+    estimated_bg = (np.ones_like(background_img) * background_img.mean()).astype(
+        np.uint16
+    )
+    estimated_bg[hot_pixels > 0] = 500
+
+    # Detect spots without estimated background
+    spots = detect_spots_2d(
+        img=img_final,
+        lateral_sigma=0.75,
+        h=90,
+        background_img=None,
+    )
+    assert spots.shape[0] == 4
+    assert_array_equal(
+        spots,
+        np.array(
+            [
+                [25, 25],
+                [25, 50],
+                [50, 50],
+                [75, 75],
+            ]
+        ),
+    )
+
+    # Detect spot with estimated background
+    spots = detect_spots_2d(
+        img=img_final,
+        background_img=estimated_bg,
+        lateral_sigma=0.75,
+        h=90,
+    )
+    assert spots.shape[0] == 3
+    assert_array_equal(
+        spots,
+        np.array(
+            [
+                [25, 25],
+                [50, 50],
+                [75, 75],
+            ]
+        ),
+    )
+
+    # Detect spots brighter than 190
+    spots = detect_spots_2d(
+        img=img_final,
+        background_img=estimated_bg,
+        lateral_sigma=0.75,
+        h=190,
+    )
+    assert spots.shape[0] == 2
+    assert_array_equal(
+        spots,
+        np.array(
+            [
+                [50, 50],
+                [75, 75],
+            ]
+        ),
+    )
+
+    # Detect spots brighter than 290
+    spots = detect_spots_2d(
+        img=img_final,
+        background_img=estimated_bg,
+        lateral_sigma=0.75,
+        h=290,
+    )
+    assert spots.shape[0] == 1
+    assert_array_equal(
+        spots,
+        np.array(
+            [
+                [75, 75],
+            ]
+        ),
+    )