diff --git a/src/faim_ipa/detection/blobs.py b/src/faim_ipa/detection/blobs.py index c0e46f2a..9512c832 100644 --- a/src/faim_ipa/detection/blobs.py +++ b/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 = [ diff --git a/src/faim_ipa/detection/spots.py b/src/faim_ipa/detection/spots.py index f3267ccd..a9d0baa8 100644 --- a/src/faim_ipa/detection/spots.py +++ b/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 diff --git a/src/faim_ipa/detection/utils.py b/src/faim_ipa/detection/utils.py index 1349ae40..aa9b3c3b 100644 --- a/src/faim_ipa/detection/utils.py +++ b/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() diff --git a/tests/detection/test_detection_utils.py b/tests/detection/test_detection_utils.py index ba5f26cd..949c86b9 100644 --- a/tests/detection/test_detection_utils.py +++ b/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 + ) diff --git a/tests/detection/test_spots.py b/tests/detection/test_spots.py index 0f3b8219..4bad510e 100644 --- a/tests/detection/test_spots.py +++ b/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], + ] + ), + )