diff --git a/src/faim_ipa/detection/__init__.py b/src/faim_ipa/detection/__init__.py
new file mode 100644
index 00000000..e69de29b
diff --git a/src/faim_ipa/detection/blobs.py b/src/faim_ipa/detection/blobs.py
new file mode 100644
index 00000000..ae6a36f5
--- /dev/null
+++ b/src/faim_ipa/detection/blobs.py
@@ -0,0 +1,96 @@
+from typing import Optional
+
+import numpy as np
+from scipy.ndimage import gaussian_laplace
+from skimage.feature import peak_local_max
+from skimage.morphology import h_maxima, ball
+from skimage.util import img_as_float32
+from skimage.feature.blob import _prune_blobs
+
+from faim_ipa.detection.utils import estimate_log_rescale_factor
+
+
+def detect_blobs(
+    img: np.ndarray,
+    axial_sigma: float,
+    lateral_sigma: float,
+    h: int,
+    n_scale_levels: int,
+    overlap: float,
+    background_img: Optional[np.ndarray] = None,
+) -> np.ndarray:
+    """Detect blobs of different sizes.
+
+    The blob detection finds blobs of different sizes by applying a
+    Laplacian of Gaussian with increasing sigmas followed by h-maxima
+    filtering. The blob detection starts with the provided sigmas and
+    then doubles them for `n_scale_levels`.
+
+    Parameters
+    ----------
+    img :
+        Image containing spot signal.
+    axial_sigma :
+        Z extension of the spots.
+    lateral_sigma :
+        YX extension of the spots.
+    h :
+        h-maxima threshold.
+    n_scale_levels :
+        Number of upscaling rounds.
+    overlap :
+        A value between 0 and 1. If the fraction of area overlapping for 2
+        blobs is greater than `overlap` the smaller blob is eliminated.
+    background_img :
+        Estimated background image. This is subtracted before the
+        blob detection.
+
+    Returns
+    -------
+    Detected spots.
+    """
+    if background_img is not None:
+        image = img_as_float32(img) - img_as_float32(background_img)
+    else:
+        image = img_as_float32(img)
+
+    rescale_factor = estimate_log_rescale_factor(
+        axial_sigma=axial_sigma, lateral_sigma=lateral_sigma
+    )
+
+    sigmas = [
+        (axial_sigma * 2**i, lateral_sigma * 2**i, lateral_sigma * 2**i)
+        for i in range(n_scale_levels)
+    ]
+
+    scale_cube = np.empty(image.shape + (len(sigmas),), dtype=np.uint8)
+
+    h_ = img_as_float32(np.array(h, dtype=img.dtype))
+    for i, sigma in enumerate(sigmas):
+        log_img = (
+            -gaussian_laplace(image, sigma=sigma)
+            * rescale_factor
+            * (np.mean(sigma) / np.mean(sigmas[0])) ** 2
+        )
+        scale_cube[..., i] = h_maxima(log_img, h=h_, footprint=ball(1))
+
+    maxima = peak_local_max(
+        scale_cube,
+        threshold_abs=0.1,
+        exclude_border=False,
+        footprint=np.ones((3,) * scale_cube.ndim),
+    )
+
+    # Convert local_maxima to float64
+    lm = maxima.astype(np.float64)
+
+    # translate final column of lm, which contains the index of the
+    # sigma that produced the maximum intensity value, into the sigma
+    sigmas_of_peaks = np.array(sigmas)[maxima[:, -1]]
+
+    # Remove sigma index and replace with sigmas
+    lm = np.hstack([lm[:, :-1], sigmas_of_peaks])
+
+    sigma_dim = sigmas_of_peaks.shape[1]
+
+    return _prune_blobs(np.array(lm), overlap=overlap, sigma_dim=sigma_dim)
diff --git a/src/faim_ipa/detection/spots.py b/src/faim_ipa/detection/spots.py
new file mode 100644
index 00000000..9ff97c8c
--- /dev/null
+++ b/src/faim_ipa/detection/spots.py
@@ -0,0 +1,57 @@
+from typing import Optional
+
+import numpy as np
+from scipy.ndimage import gaussian_laplace
+from skimage.morphology import h_maxima, ball
+from skimage.util import img_as_float32
+
+from faim_ipa.detection.utils import estimate_log_rescale_factor
+
+
+def detect_spots(
+    img: np.ndarray,
+    axial_sigma: float,
+    lateral_sigma: float,
+    h: int,
+    background_img: Optional[np.ndarray] = 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.
+    axial_sigma :
+        Z extension of the spots.
+    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.
+    """
+    if background_img is not None:
+        image = img_as_float32(img) - img_as_float32(background_img)
+    else:
+        image = img_as_float32(img)
+
+    rescale_factor = estimate_log_rescale_factor(
+        axial_sigma=axial_sigma, lateral_sigma=lateral_sigma
+    )
+    log_img = (
+        -gaussian_laplace(image, sigma=(axial_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=ball(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
new file mode 100644
index 00000000..a12140b9
--- /dev/null
+++ b/src/faim_ipa/detection/utils.py
@@ -0,0 +1,77 @@
+import numpy as np
+from scipy.ndimage import gaussian_laplace, gaussian_filter
+
+
+def compute_axial_sigma(wavelength: float, NA: float, axial_spacing: float):
+    """
+    Sigma which produces a Gaussian with the same full width
+    half maximum as described by Abbe's diffraction formula for axial resolution.
+
+    R = 2 * lambda / (NA**2)
+
+    Parameters
+    ----------
+    wavelength :
+        Emission wavelength
+    NA :
+        Numerical Aperture
+    axial_spacing :
+        Spacing of z planes in nanometers
+
+    Returns
+    -------
+    theoretical sigma in pixels
+    """
+    return 2 * wavelength / (NA**2) / (2 * np.sqrt(2 * np.log(2))) / axial_spacing
+
+
+def compute_lateral_sigma(wavelength: float, NA: float, lateral_spacing: float):
+    """
+    Sigma which produces a Gaussian with the same full width
+    half maximum as the theoretical resolution limit in Y/X described by E. Abbe.
+
+    d = lambda / (2*NA)
+
+    Parameters
+    ----------
+    wavelength :
+        Emission wavelength
+    NA :
+        Numerical Aperture
+    lateral_spacing :
+        Pixel size in YX.
+
+    Returns
+    -------
+    theoretical sigma in pixels
+    """
+    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:
+    """
+    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
+    (axial_sigma, lateral_sigma, lateral_sigma).
+
+    Parameters
+    ----------
+    axial_sigma :
+        Sigma in axial direction. Usually along Z.
+    lateral_sigma :
+        Sigma in lateral direction. Usually along Y and X.
+
+    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))
+    img = img / img.max()
+    img_log = -gaussian_laplace(
+        input=img, sigma=(axial_sigma, lateral_sigma, lateral_sigma)
+    )
+    return 1 / img_log.max()
diff --git a/tests/detection/__init__.py b/tests/detection/__init__.py
new file mode 100644
index 00000000..e69de29b
diff --git a/tests/detection/test_blobs.py b/tests/detection/test_blobs.py
new file mode 100644
index 00000000..2870c62f
--- /dev/null
+++ b/tests/detection/test_blobs.py
@@ -0,0 +1,86 @@
+import numpy as np
+from numpy.testing import assert_array_equal
+from scipy.ndimage import gaussian_filter
+
+from faim_ipa.detection.blobs import detect_blobs
+
+
+def test_detect_blobs():
+    np.random.seed(0)
+    img = np.zeros((101, 101, 101), dtype=np.float32)
+
+    # Create 3 spots with intensities of 100, 200, 300
+    img[25, 25, 25] = 1
+    img[50, 50, 50] = 1
+    img[75, 75, 75] = 1
+    img = gaussian_filter(img, (2.07, 0.75, 0.75))
+    img = 300 * img / img.max()
+
+    # Create 1 larger spot
+    large_spot = np.zeros((101, 101, 101), dtype=np.float32)
+    large_spot[50, 75, 75] = 1
+    large_spot = gaussian_filter(large_spot, (2 * 2.07, 2 * 0.75, 2 * 0.75))
+    large_spot = 300 * large_spot / large_spot.max()
+
+    # Create random background noise
+    background_img = np.random.normal(10, 2, (101, 101, 101)).astype(np.float32)
+
+    # Create hot-pixel
+    hot_pixels = np.zeros((101, 101, 101), dtype=np.float32)
+    hot_pixels[25, 25, 50] = 1000
+
+    # 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] = 1000
+
+    # Detect spots without estimated background
+    blobs = detect_blobs(
+        img=img_final,
+        axial_sigma=2.07,
+        lateral_sigma=0.75,
+        h=200,
+        n_scale_levels=2,
+        overlap=0.875,
+        background_img=None,
+    )
+    assert blobs.shape[0] == 5
+    assert_array_equal(
+        blobs,
+        np.array(
+            [
+                [25, 25, 25, 2.07, 0.75, 0.75],
+                [25, 25, 50, 2.07, 0.75, 0.75],
+                [50, 50, 50, 2.07, 0.75, 0.75],
+                [50, 75, 75, 2 * 2.07, 2 * 0.75, 2 * 0.75],
+                [75, 75, 75, 2.07, 0.75, 0.75],
+            ]
+        ),
+    )
+
+    # Detect spots with estimated background
+    blobs = detect_blobs(
+        img=img_final,
+        axial_sigma=2.07,
+        lateral_sigma=0.75,
+        h=200,
+        n_scale_levels=2,
+        overlap=0.875,
+        background_img=estimated_bg,
+    )
+    assert blobs.shape[0] == 4
+    assert_array_equal(
+        blobs,
+        np.array(
+            [
+                [25, 25, 25, 2.07, 0.75, 0.75],
+                [50, 50, 50, 2.07, 0.75, 0.75],
+                [50, 75, 75, 2 * 2.07, 2 * 0.75, 2 * 0.75],
+                [75, 75, 75, 2.07, 0.75, 0.75],
+            ]
+        ),
+    )
diff --git a/tests/detection/test_spots.py b/tests/detection/test_spots.py
new file mode 100644
index 00000000..0f3b8219
--- /dev/null
+++ b/tests/detection/test_spots.py
@@ -0,0 +1,117 @@
+import numpy as np
+from numpy.testing import assert_array_equal
+from scipy.ndimage import gaussian_filter
+
+from faim_ipa.detection.spots import detect_spots
+
+
+def test_detect_spots():
+    np.random.seed(0)
+    img = np.zeros((101, 101, 101), dtype=np.float32)
+
+    # Create 3 spots with intensities of 100, 200, 300
+    img[25, 25, 25] = 1
+    img[50, 50, 50] = 2
+    img[75, 75, 75] = 3
+    img = gaussian_filter(img, (2.07, 0.75, 0.75))
+    img = 300 * img / img.max()
+
+    # Create 1 larger spot
+    large_spot = np.zeros((101, 101, 101), dtype=np.float32)
+    large_spot[50, 75, 75] = 1
+    large_spot = gaussian_filter(large_spot, (2 * 2.07, 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, 101)).astype(np.float32)
+
+    # Create hot-pixel
+    hot_pixels = np.zeros((101, 101, 101), dtype=np.float32)
+    hot_pixels[25, 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(
+        img=img_final,
+        axial_sigma=2.07,
+        lateral_sigma=0.75,
+        h=90,
+        background_img=None,
+    )
+    assert spots.shape[0] == 4
+    assert_array_equal(
+        spots,
+        np.array(
+            [
+                [25, 25, 25],
+                [25, 25, 50],
+                [50, 50, 50],
+                [75, 75, 75],
+            ]
+        ),
+    )
+
+    # Detect spot with estimated background
+    spots = detect_spots(
+        img=img_final,
+        background_img=estimated_bg,
+        axial_sigma=2.07,
+        lateral_sigma=0.75,
+        h=90,
+    )
+    assert spots.shape[0] == 3
+    assert_array_equal(
+        spots,
+        np.array(
+            [
+                [25, 25, 25],
+                [50, 50, 50],
+                [75, 75, 75],
+            ]
+        ),
+    )
+
+    # Detect spots brighter than 190
+    spots = detect_spots(
+        img=img_final,
+        background_img=estimated_bg,
+        axial_sigma=2.07,
+        lateral_sigma=0.75,
+        h=190,
+    )
+    assert spots.shape[0] == 2
+    assert_array_equal(
+        spots,
+        np.array(
+            [
+                [50, 50, 50],
+                [75, 75, 75],
+            ]
+        ),
+    )
+
+    # Detect spots brighter than 290
+    spots = detect_spots(
+        img=img_final,
+        background_img=estimated_bg,
+        axial_sigma=2.07,
+        lateral_sigma=0.75,
+        h=290,
+    )
+    assert spots.shape[0] == 1
+    assert_array_equal(
+        spots,
+        np.array(
+            [
+                [75, 75, 75],
+            ]
+        ),
+    )
diff --git a/tests/detection/test_utils.py b/tests/detection/test_utils.py
new file mode 100644
index 00000000..bd7adab3
--- /dev/null
+++ b/tests/detection/test_utils.py
@@ -0,0 +1,40 @@
+import numpy as np
+from numpy.testing import assert_almost_equal
+from scipy.ndimage import gaussian_filter, gaussian_laplace
+
+from faim_ipa.detection.utils import (
+    compute_axial_sigma,
+    compute_lateral_sigma,
+    estimate_log_rescale_factor,
+)
+
+
+def test_compute_axial_sigma():
+    wl = 514
+    NA = 1.45
+    axial_spacing = 100
+
+    assert_almost_equal(compute_axial_sigma(wl, NA, axial_spacing), 2.08, 2)
+
+
+def test_compute_lateral_sigma():
+    wl = 514
+    NA = 1.45
+    lateral_spacing = 100
+    print(compute_lateral_sigma(wl, NA, lateral_spacing))
+    assert_almost_equal(compute_lateral_sigma(wl, NA, lateral_spacing), 0.75, 2)
+
+
+def test_estimate_log_rescale_factor():
+    rescale_factor = estimate_log_rescale_factor(2.07, 0.75)
+
+    # Create image with diffraction limited spot with intensity = 10
+    img = np.zeros((101, 101, 101), dtype=np.float32)
+    img[50, 50, 50] = 1
+    img = gaussian_filter(img, (2.07, 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, (2.07, 0.75, 0.75)) * rescale_factor), 10, 2
+    )