Mercurial > repos > imgteam > spot_detection_2d

comparison spot_detection_2d.py @ 5:e8c9e104e109 draft default tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
planemo upload for repository https://github.com/BMCV/galaxy-image-analysis/tree/master/tools/spot_detection_2d/ commit f1b9207ec23c0af1681c929281bcbf1d0638368e
author imgteam
date Wed, 25 Sep 2024 08:19:30 +0000
parents 4645b356bf3b
children
comparison
equal deleted inserted replaced
4:14f9986800fa 5:e8c9e104e109
1 """ 1 """
2 Copyright 2021-2022 Biomedical Computer Vision Group, Heidelberg University. 2 Copyright 2021-2022 Biomedical Computer Vision Group, Heidelberg University.
3 Author: Qi Gao (qi.gao@bioquant.uni-heidelberg.de) 3 Authors:
4 - Qi Gao (qi.gao@bioquant.uni-heidelberg.de)
5 - Leonid Kostrykin (leonid.kostrykin@bioquant.uni-heidelberg.de)
4 6
5 Distributed under the MIT license. 7 Distributed under the MIT license.
6 See file LICENSE for detail or copy at https://opensource.org/licenses/MIT 8 See file LICENSE for detail or copy at https://opensource.org/licenses/MIT
7
8 """ 9 """
9 10
10 import argparse 11 import argparse
11 12
12 import imageio 13 import giatools.io
13 import numpy as np 14 import numpy as np
14 import pandas as pd 15 import pandas as pd
15 from skimage.feature import peak_local_max 16 import scipy.ndimage as ndi
16 from skimage.filters import gaussian, laplace 17 from numpy.typing import NDArray
18 from skimage.feature import blob_dog, blob_doh, blob_log
19
20 blob_filters = {
21 'dog': blob_dog,
22 'doh': blob_doh,
23 'log': blob_log,
24 }
17 25
18 26
19 def getbr(xy, img, nb, firstn): 27 def mean_intensity(img: NDArray, y: int, x: int, radius: int) -> float:
20 ndata = xy.shape[0] 28 assert img.ndim == 2
21 br = np.empty((ndata, 1)) 29 assert radius >= 0
22 for j in range(ndata): 30 if radius == 0:
23 br[j] = np.NaN 31 return float(img[y, x])
24 if not np.isnan(xy[j, 0]): 32 else:
25 timg = img[xy[j, 1] - nb - 1:xy[j, 1] + nb, xy[j, 0] - nb - 1:xy[j, 0] + nb] 33 mask = np.ones(img.shape, bool)
26 br[j] = np.mean(np.sort(timg, axis=None)[-firstn:]) 34 mask[y, x] = False
27 return br 35 mask = (ndi.distance_transform_edt(mask) <= radius)
36 return img[mask].mean()
28 37
29 38
30 def spot_detection(fn_in, fn_out, frame_1st=1, frame_end=0, 39 def spot_detection(
31 typ_filter='Gauss', ssig=1, th=10, 40 fn_in: str,
32 typ_br='smoothed', bd=10): 41 fn_out: str,
33 ims_ori = imageio.mimread(fn_in, format='TIFF') 42 frame_1st: int,
34 ims_smd = np.zeros((len(ims_ori), ims_ori[0].shape[0], ims_ori[0].shape[1]), dtype='float64') 43 frame_end: int,
35 if frame_end == 0 or frame_end > len(ims_ori): 44 filter_type: str,
36 frame_end = len(ims_ori) 45 min_scale: float,
46 max_scale: float,
47 abs_threshold: float,
48 rel_threshold: float,
49 boundary: int,
50 ) -> None:
37 51
38 for i in range(frame_1st - 1, frame_end): 52 # Load the single-channel 2-D input image (or stack thereof)
39 ims_smd[i, :, :] = gaussian(ims_ori[i].astype('float64'), sigma=ssig) 53 stack = giatools.io.imread(fn_in)
40 ims_smd_max = np.max(ims_smd)
41 54
42 txyb_all = np.array([]).reshape(0, 4) 55 # Normalize input image so that it is a stack of images (possibly a stack of a single image)
43 for i in range(frame_1st - 1, frame_end): 56 assert stack.ndim in (2, 3)
44 tmp = np.copy(ims_smd[i, :, :]) 57 if stack.ndim == 2:
45 if typ_filter == 'LoG': 58 stack = stack.reshape(1, *stack.shape)
46 tmp = laplace(tmp)
47 59
48 tmp[tmp < th * ims_smd_max / 100] = 0 60 # Slice the stack
49 coords = peak_local_max(tmp, min_distance=1) 61 assert frame_1st >= 1
50 idx_to_del = np.where((coords[:, 0] <= bd) | (coords[:, 0] >= tmp.shape[0] - bd) | 62 assert frame_end >= 0
51 (coords[:, 1] <= bd) | (coords[:, 1] >= tmp.shape[1] - bd)) 63 stack = stack[frame_1st - 1:]
52 coords = np.delete(coords, idx_to_del[0], axis=0) 64 if frame_end > 0:
53 xys = coords[:, ::-1] 65 stack = stack[:-frame_end]
54 66
55 if typ_br == 'smoothed': 67 # Select the blob detection filter
56 intens = getbr(xys, ims_smd[i, :, :], 0, 1) 68 assert filter_type.lower() in blob_filters.keys()
57 elif typ_br == 'robust': 69 blob_filter = blob_filters[filter_type.lower()]
58 intens = getbr(xys, ims_ori[i], 1, 4)
59 else:
60 intens = getbr(xys, ims_ori[i], 0, 1)
61 70
62 txyb = np.concatenate(((i + 1) * np.ones((xys.shape[0], 1)), xys, intens), axis=1) 71 # Perform blob detection on each image of the stack
63 txyb_all = np.concatenate((txyb_all, txyb), axis=0) 72 detections = list()
73 for img_idx, img in enumerate(stack):
74 blobs = blob_filter(img, threshold=abs_threshold, threshold_rel=rel_threshold, min_sigma=min_scale, max_sigma=max_scale)
75 for blob in blobs:
76 y, x, scale = blob
64 77
65 df = pd.DataFrame() 78 # Skip the detection if it is too close to the boundary of the image
66 df['FRAME'] = txyb_all[:, 0].astype(int) 79 if y < boundary or x < boundary or y >= img.shape[0] - boundary or x >= img.shape[1] - boundary:
67 df['POS_X'] = txyb_all[:, 1].astype(int) 80 continue
68 df['POS_Y'] = txyb_all[:, 2].astype(int) 81
69 df['INTENSITY'] = txyb_all[:, 3] 82 # Add the detection to the list of detections
83 radius = scale * np.sqrt(2) * 2
84 intensity = mean_intensity(img, round(y), round(x), round(radius))
85 detections.append(
86 {
87 'frame': img_idx + 1,
88 'pos_x': round(x),
89 'pos_y': round(y),
90 'scale': scale,
91 'radius': radius,
92 'intensity': intensity,
93 }
94 )
95
96 # Build and save dataframe
97 df = pd.DataFrame.from_dict(detections)
70 df.to_csv(fn_out, index=False, float_format='%.2f', sep="\t") 98 df.to_csv(fn_out, index=False, float_format='%.2f', sep="\t")
71 99
72 100
73 if __name__ == "__main__": 101 if __name__ == "__main__":
102
74 parser = argparse.ArgumentParser(description="Spot detection") 103 parser = argparse.ArgumentParser(description="Spot detection")
75 parser.add_argument("fn_in", help="Name of input image sequence (stack)") 104
76 parser.add_argument("fn_out", help="Name of output file to save the coordinates and intensities of detected spots") 105 parser.add_argument("fn_in", help="Name of input image or image stack.")
77 parser.add_argument("frame_1st", type=int, help="Index for the starting frame to detect spots (1 for first frame of the stack)") 106 parser.add_argument("fn_out", help="Name of output file to write the detections into.")
78 parser.add_argument("frame_end", type=int, help="Index for the last frame to detect spots (0 for the last frame of the stack)") 107 parser.add_argument("frame_1st", type=int, help="Index for the starting frame to detect spots (1 for first frame of the stack).")
79 parser.add_argument("filter", help="Detection filter") 108 parser.add_argument("frame_end", type=int, help="Index for the last frame to detect spots (0 for the last frame of the stack).")
80 parser.add_argument("ssig", type=float, help="Sigma of the Gaussian for noise suppression") 109 parser.add_argument("filter_type", help="Detection filter")
81 parser.add_argument("thres", type=float, help="Percentage of the global maximal for thresholding candidate spots") 110 parser.add_argument("min_scale", type=float, help="The minimum scale to consider for multi-scale detection.")
82 parser.add_argument("typ_intens", help="smoothed or robust (for measuring the intensities of spots)") 111 parser.add_argument("max_scale", type=float, help="The maximum scale to consider for multi-scale detection.")
83 parser.add_argument("bndy", type=int, help="Number of pixels (Spots close to image boundaries will be ignored)") 112 parser.add_argument("abs_threshold", type=float, help=(
113 "Filter responses below this threshold will be ignored. Only filter responses above this thresholding will be considered as blobs. "
114 "This threshold is ignored if the relative threshold (below) corresponds to a higher response.")
115 )
116 parser.add_argument("rel_threshold", type=float, help=(
117 "Same as the absolute threshold (above), but as a fraction of the overall maximal filter response of an image. "
118 "This threshold is ignored if it corresponds to a response below the absolute threshold.")
119 )
120 parser.add_argument("boundary", type=int, help="Width of image boundaries (in pixel) where spots will be ignored.")
121
84 args = parser.parse_args() 122 args = parser.parse_args()
85 spot_detection(args.fn_in, args.fn_out, 123 spot_detection(args.fn_in, args.fn_out,
86 frame_1st=args.frame_1st, frame_end=args.frame_end, 124 frame_1st=args.frame_1st, frame_end=args.frame_end,
87 typ_filter=args.filter, ssig=args.ssig, th=args.thres, 125 filter_type=args.filter_type,
88 typ_br=args.typ_intens, bd=args.bndy) 126 min_scale=args.min_scale, max_scale=args.max_scale,
127 abs_threshold=args.abs_threshold, rel_threshold=args.rel_threshold,
128 boundary=args.boundary)