Overhaul the spot_detection_2d tool and add support for single-image processing

tools/spot_detection_2d/creators.xml

@@ -0,0 +1 @@
+../../macros/creators.xml

tools/spot_detection_2d/spot_detection_2d.py

@@ -1,88 +1,128 @@
 """
 Copyright 2021-2022 Biomedical Computer Vision Group, Heidelberg University.
-Author: Qi Gao (qi.gao@bioquant.uni-heidelberg.de)
+Authors:
+- Qi Gao (qi.gao@bioquant.uni-heidelberg.de)
+- Leonid Kostrykin (leonid.kostrykin@bioquant.uni-heidelberg.de)
 
 Distributed under the MIT license.
 See file LICENSE for detail or copy at https://opensource.org/licenses/MIT
-
 """
 
 import argparse
 
-import imageio
+import giatools.io
 import numpy as np
 import pandas as pd
-from skimage.feature import peak_local_max
-from skimage.filters import gaussian, laplace
-
-
-def getbr(xy, img, nb, firstn):
-    ndata = xy.shape[0]
-    br = np.empty((ndata, 1))
-    for j in range(ndata):
-        br[j] = np.NaN
-        if not np.isnan(xy[j, 0]):
-            timg = img[xy[j, 1] - nb - 1:xy[j, 1] + nb, xy[j, 0] - nb - 1:xy[j, 0] + nb]
-            br[j] = np.mean(np.sort(timg, axis=None)[-firstn:])
-    return br
-
-
-def spot_detection(fn_in, fn_out, frame_1st=1, frame_end=0,
-                   typ_filter='Gauss', ssig=1, th=10,
-                   typ_br='smoothed', bd=10):
-    ims_ori = imageio.mimread(fn_in, format='TIFF')
-    ims_smd = np.zeros((len(ims_ori), ims_ori[0].shape[0], ims_ori[0].shape[1]), dtype='float64')
-    if frame_end == 0 or frame_end > len(ims_ori):
-        frame_end = len(ims_ori)
-
-    for i in range(frame_1st - 1, frame_end):
-        ims_smd[i, :, :] = gaussian(ims_ori[i].astype('float64'), sigma=ssig)
-    ims_smd_max = np.max(ims_smd)
-
-    txyb_all = np.array([]).reshape(0, 4)
-    for i in range(frame_1st - 1, frame_end):
-        tmp = np.copy(ims_smd[i, :, :])
-        if typ_filter == 'LoG':
-            tmp = laplace(tmp)
-
-        tmp[tmp < th * ims_smd_max / 100] = 0
-        coords = peak_local_max(tmp, min_distance=1)
-        idx_to_del = np.where((coords[:, 0] <= bd) | (coords[:, 0] >= tmp.shape[0] - bd) |
-                              (coords[:, 1] <= bd) | (coords[:, 1] >= tmp.shape[1] - bd))
-        coords = np.delete(coords, idx_to_del[0], axis=0)
-        xys = coords[:, ::-1]
-
-        if typ_br == 'smoothed':
-            intens = getbr(xys, ims_smd[i, :, :], 0, 1)
-        elif typ_br == 'robust':
-            intens = getbr(xys, ims_ori[i], 1, 4)
-        else:
-            intens = getbr(xys, ims_ori[i], 0, 1)
-
-        txyb = np.concatenate(((i + 1) * np.ones((xys.shape[0], 1)), xys, intens), axis=1)
-        txyb_all = np.concatenate((txyb_all, txyb), axis=0)
-
-    df = pd.DataFrame()
-    df['FRAME'] = txyb_all[:, 0].astype(int)
-    df['POS_X'] = txyb_all[:, 1].astype(int)
-    df['POS_Y'] = txyb_all[:, 2].astype(int)
-    df['INTENSITY'] = txyb_all[:, 3]
+import scipy.ndimage as ndi
+from numpy.typing import NDArray
+from skimage.feature import blob_dog, blob_doh, blob_log
+
+blob_filters = {
+    'dog': blob_dog,
+    'doh': blob_doh,
+    'log': blob_log,
+}
+
+
+def mean_intensity(img: NDArray, y: int, x: int, radius: int) -> float:
+    assert img.ndim == 2
+    assert radius >= 0
+    if radius == 0:
+        return float(img[y, x])
+    else:
+        mask = np.ones(img.shape, bool)
+        mask[y, x] = False
+        mask = (ndi.distance_transform_edt(mask) <= radius)
+        return img[mask].mean()
+
+
+def spot_detection(
+    fn_in: str,
+    fn_out: str,
+    frame_1st: int,
+    frame_end: int,
+    filter_type: str,
+    min_scale: float,
+    max_scale: float,
+    abs_threshold: float,
+    rel_threshold: float,
+    boundary: int,
+) -> None:
+
+    # Load the single-channel 2-D input image (or stack thereof)
+    stack = giatools.io.imread(fn_in)
+
+    # Normalize input image so that it is a stack of images (possibly a stack of a single image)
+    assert stack.ndim in (2, 3)
+    if stack.ndim == 2:
+        stack = stack.reshape(1, *stack.shape)
+
+    # Slice the stack
+    assert frame_1st >= 1
+    assert frame_end >= 0
+    stack = stack[frame_1st - 1:]
+    if frame_end > 0:
+        stack = stack[:-frame_end]
+
+    # Select the blob detection filter
+    assert filter_type.lower() in blob_filters.keys()
+    blob_filter = blob_filters[filter_type.lower()]
+
+    # Perform blob detection on each image of the stack
+    detections = list()
+    for img_idx, img in enumerate(stack):
+        blobs = blob_filter(img, threshold=abs_threshold, threshold_rel=rel_threshold, min_sigma=min_scale, max_sigma=max_scale)
+        for blob in blobs:
+            y, x, scale = blob
+
+            # Skip the detection if it is too close to the boundary of the image
+            if y < boundary or x < boundary or y >= img.shape[0] - boundary or x >= img.shape[1] - boundary:
+                continue
+
+            # Add the detection to the list of detections
+            radius = scale * np.sqrt(2) * 2
+            intensity = mean_intensity(img, round(y), round(x), round(radius))
+            detections.append(
+                {
+                    'frame': img_idx + 1,
+                    'pos_x': round(x),
+                    'pos_y': round(y),
+                    'scale': scale,
+                    'radius': radius,
+                    'intensity': intensity,
+                }
+            )
+
+    # Build and save dataframe
+    df = pd.DataFrame.from_dict(detections)
     df.to_csv(fn_out, index=False, float_format='%.2f', sep="\t")
 
 
 if __name__ == "__main__":
+
     parser = argparse.ArgumentParser(description="Spot detection")
-    parser.add_argument("fn_in", help="Name of input image sequence (stack)")
-    parser.add_argument("fn_out", help="Name of output file to save the coordinates and intensities of detected spots")
-    parser.add_argument("frame_1st", type=int, help="Index for the starting frame to detect spots (1 for first frame of the stack)")
-    parser.add_argument("frame_end", type=int, help="Index for the last frame to detect spots (0 for the last frame of the stack)")
-    parser.add_argument("filter", help="Detection filter")
-    parser.add_argument("ssig", type=float, help="Sigma of the Gaussian for noise suppression")
-    parser.add_argument("thres", type=float, help="Percentage of the global maximal for thresholding candidate spots")
-    parser.add_argument("typ_intens", help="smoothed or robust (for measuring the intensities of spots)")
-    parser.add_argument("bndy", type=int, help="Number of pixels (Spots close to image boundaries will be ignored)")
+
+    parser.add_argument("fn_in", help="Name of input image or image stack.")
+    parser.add_argument("fn_out", help="Name of output file to write the detections into.")
+    parser.add_argument("frame_1st", type=int, help="Index for the starting frame to detect spots (1 for first frame of the stack).")
+    parser.add_argument("frame_end", type=int, help="Index for the last frame to detect spots (0 for the last frame of the stack).")
+    parser.add_argument("filter_type", help="Detection filter")
+    parser.add_argument("min_scale", type=float, help="The minimum scale to consider for multi-scale detection.")
+    parser.add_argument("max_scale", type=float, help="The maximum scale to consider for multi-scale detection.")
+    parser.add_argument("abs_threshold", type=float, help=(
+        "Filter responses below this threshold will be ignored. Only filter responses above this thresholding will be considered as blobs. "
+        "This threshold is ignored if the relative threshold (below) corresponds to a higher response.")
+    )
+    parser.add_argument("rel_threshold", type=float, help=(
+        "Same as the absolute threshold (above), but as a fraction of the overall maximal filter response of an image. "
+        "This threshold is ignored if it corresponds to a response below the absolute threshold.")
+    )
+    parser.add_argument("boundary", type=int, help="Width of image boundaries (in pixel) where spots will be ignored.")
+
     args = parser.parse_args()
     spot_detection(args.fn_in, args.fn_out,
                    frame_1st=args.frame_1st, frame_end=args.frame_end,
-                   typ_filter=args.filter, ssig=args.ssig, th=args.thres,
-                   typ_br=args.typ_intens, bd=args.bndy)
+                   filter_type=args.filter_type,
+                   min_scale=args.min_scale, max_scale=args.max_scale,
+                   abs_threshold=args.abs_threshold, rel_threshold=args.rel_threshold,
+                   boundary=args.boundary)

tools/spot_detection_2d/spot_detection_2d.xml

@@ -1,16 +1,25 @@
-<tool id="ip_spot_detection_2d" name="Perform 2D spot detection" version="0.0.3-2" profile="20.05">
+<tool id="ip_spot_detection_2d" name="Perform 2-D spot detection" version="@TOOL_VERSION@+galaxy@VERSION_SUFFIX@" profile="20.05">
     <description></description>
+    <macros>
+        <import>creators.xml</import>
+        <token name="@TOOL_VERSION@">0.1</token>
+        <token name="@VERSION_SUFFIX@">0</token>
+    </macros>
+    <creator>
+        <expand macro="creators/bmcv" />
+    </creator>
     <edam_operations>
         <edam_operation>operation_3443</edam_operation>
     </edam_operations>
     <xrefs>
         <xref type="bio.tools">galaxy_image_analysis</xref>
     </xrefs>
     <requirements>
-        <requirement type="package" version="2.9.0">imageio</requirement>
-        <requirement type="package" version="1.20.2">numpy</requirement>
+        <requirement type="package" version="0.2.0">giatools</requirement>
+        <requirement type="package" version="1.26.4">numpy</requirement>
         <requirement type="package" version="1.2.4">pandas</requirement>
-        <requirement type="package" version="0.18.1">scikit-image</requirement>
+        <requirement type="package" version="0.21">scikit-image</requirement>
+        <requirement type="package" version="2024.6.18">tifffile</requirement>
     </requirements>
     <command detect_errors="aggressive">
     <![CDATA[
@@ -19,49 +28,75 @@
          '$fn_out'
          '$frame_1st'
          '$frame_end'
-         '$filter'
-         '$ssig'
-         '$thres'
-         '$typ_intens'
-         '$bndy'
+         '$filter_type'
+         '$min_scale'
+         '$max_scale'
+         '$abs_threshold'
+         '$rel_threshold'
+         '$boundary'
     ]]>
     </command>
     <inputs>
-        <param name="fn_in" type="data" format="tiff" label="Image sequence (stack)" />
+        <param name="fn_in" type="data" format="tiff" label="Intensity image or a stack of images" />
         <param name="frame_1st" type="integer" value="1" label="Starting time point (1 for the first frame of the stack)" />
         <param name="frame_end" type="integer" value="0" label="Ending time point (0 for the last frame of the stack)" />
-        <param name="filter" type="select" label="Choose a detection filter">
-            <option value="Gauss" selected="True">Gaussian</option>
-            <option value="LoG">Laplacian-of-Gaussian, LoG (more accurate when spots have similar sizes)</option>
+        <param name="filter_type" type="select" label="Detection filter">
+            <option value="LoG" selected="True">Laplacian of Gaussian</option>
+            <option value="DoG">Difference of Gaussians</option>
+            <option value="DoH">Determinant of Hessian</option>
         </param>
-        <param name="ssig" type="float" value="1.0" min="0.5" max="10" label="Sigma of the Gaussian (for noise suppression)" />
-        <param name="thres" type="float" value="10.0" min="0" max="100" label="Percentage of the global maximal as the threshold for candidate spots" />
-        <param name="typ_intens" type="select" label="How to measure the intensities">
-            <option value="smoothed" selected="True">Smoothed</option>
-            <option value="robust">Robust</option>
-        </param>
-        <param name="bndy" type="integer" value="10" min="0" max="50" label="Width (in pixel) of image boundaries where spots will be ignored" />
+        <param name="min_scale" type="float" value="1.0" min="1.0" label="Minimum scale" />
+        <param name="max_scale" type="float" value="2.0" min="1.0" label="Maximum scale" />
+        <param name="abs_threshold" type="float" value=".25" min="0" label="Minimum filter response (absolute)" help="Filter responses below this threshold will be ignored. Only filter responses above this thresholding will be considered as blobs. This threshold is ignored if the relative threshold (below) corresponds to a higher response." />
+        <param name="rel_threshold" type="float" value="0" min="0" max="1" label="Minimum filter response (relative)" help="Same as the absolute threshold (above), but as a fraction of the overall maximum filter response of an image. This threshold is ignored if it corresponds to a response below the absolute threshold." />
+        <param name="boundary" type="integer" value="10" min="0" label="Image boundary" help="Width of image boundaries (in pixel) where spots will be ignored." />
     </inputs>
     <outputs>
         <data format="tabular" name="fn_out" />
     </outputs>
     <tests>
+        <!-- Multi-frame input -->
+        <test>
+            <param name="fn_in" value="input1.tif"/>
+            <param name="frame_1st" value="1"/>
+            <param name="frame_end" value="0"/>
+            <param name="filter_type" value="LoG"/>
+            <param name="min_scale" value="1"/>
+            <param name="max_scale" value="2"/>
+            <param name="abs_threshold" value="0"/>
+            <param name="rel_threshold" value="0.1"/>
+            <param name="boundary" value="10"/>
+            <output name="fn_out" value="output1.tsv" ftype="tabular" />
+        </test>
+        <!-- Single-frame input -->
         <test>
-            <param name="fn_in" value="test_img1.tif"/>
+            <param name="fn_in" value="input2.tif"/>
             <param name="frame_1st" value="1"/>
             <param name="frame_end" value="0"/>
-            <param name="filter" value="Gauss"/>
-            <param name="ssig" value="1"/>
-            <param name="thres" value="10"/>
-            <param name="typ_intens" value="smoothed"/>
-            <param name="bndy" value="10"/>
-            <output name="fn_out" value="spots_detected.tsv" ftype="tabular" />
+            <param name="filter_type" value="LoG"/>
+            <param name="min_scale" value="1"/>
+            <param name="max_scale" value="2"/>
+            <param name="abs_threshold" value="0.04"/>
+            <param name="rel_threshold" value="0"/>
+            <param name="boundary" value="10"/>
+            <output name="fn_out" value="output2.tsv" ftype="tabular" />
         </test>
     </tests>
     <help>
-    **What it does**
 
-    This tool detects spots and measures the intensities in a 2D image (sequence).
+        **Perform spot detection and measure the image intensities.**
+
+        This tool detects spots (blobs) and measures the image intensities in a single-channel 2-D image (or a stack of such images).
+
+        The tool produces a TSV file containing all detections, with the following columns:
+
+        - ``frame``: The frame of the image stack
+        - ``pos_x``: The horizontal coordinate of the detection
+        - ``pos_y``: The vertical coordinate of the detection
+        - ``scale``: The scale at which the detection was found
+        - ``radius``: The radius of the detected spot
+        - ``intensity``: The mean intensity of the spot
+
     </help>
     <citations>
         <citation type="doi">10.1097/j.pain.0000000000002642</citation>