comparison spot_detection_2d.py @ 5:e8c9e104e109 draft default tip

planemo upload for repository https://github.com/BMCV/galaxy-image-analysis/tree/master/tools/spot_detection_2d/ commit f1b9207ec23c0af1681c929281bcbf1d0638368e

4:14f9986800fa

"""
Copyright 2021-2022 Biomedical Computer Vision Group, Heidelberg University.
Author: Qi Gao ([e-mail])

Distributed under the MIT license.
See file LICENSE for detail or copy at https://opensource.org/licenses/MIT

"""

import argparse

import imageio
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]
    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)")
    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)

5:e8c9e104e109

"""
Copyright 2021-2022 Biomedical Computer Vision Group, Heidelberg University.
Authors:
- Qi Gao ([e-mail])
- Leonid Kostrykin ([e-mail])

Distributed under the MIT license.
See file LICENSE for detail or copy at https://opensource.org/licenses/MIT
"""

import argparse

import giatools.io
import numpy as np
import pandas as pd
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 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,
                   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)