Mercurial > repos > imgteam > spot_detection_2d

comparison spot_detection_2d.py @ 0:d78372040976 draft

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 481cd51a76341c0ec3759f919454e95139f0cc4e"
author imgteam
date Wed, 21 Jul 2021 19:59:00 +0000
parents
children 859dd1c11ac0
comparison
equal deleted inserted replaced
-1:000000000000 0:d78372040976
1 """
2 Copyright 2021 Biomedical Computer Vision Group, Heidelberg University.
3 Author: Qi Gao (qi.gao@bioquant.uni-heidelberg.de)
4
5 Distributed under the MIT license.
6 See file LICENSE for detail or copy at https://opensource.org/licenses/MIT
7
8 """
9
10 import argparse
11
12 import imageio
13 import numpy as np
14 import pandas as pd
15 from skimage.feature import peak_local_max
16 from skimage.filters import gaussian
17
18
19 def getbr(xy, img, nb, firstn):
20 ndata = xy.shape[0]
21 br = np.empty((ndata, 1))
22 for j in range(ndata):
23 br[j] = np.NaN
24 if not np.isnan(xy[j, 0]):
25 timg = img[xy[j, 1] - nb - 1:xy[j, 1] + nb, xy[j, 0] - nb - 1:xy[j, 0] + nb]
26 br[j] = np.mean(np.sort(timg, axis=None)[-firstn:])
27 return br
28
29
30 def spot_detection(fn_in, fn_out, frame_1st=1, frame_end=0, typ_br='smoothed', th=10, ssig=1, bd=10):
31 ims_ori = imageio.mimread(fn_in, format='TIFF')
32 ims_smd = np.zeros((len(ims_ori), ims_ori[0].shape[0], ims_ori[0].shape[1]), dtype='float64')
33 if frame_end == 0 or frame_end > len(ims_ori):
34 frame_end = len(ims_ori)
35
36 for i in range(frame_1st - 1, frame_end):
37 ims_smd[i, :, :] = gaussian(ims_ori[i].astype('float64'), sigma=ssig)
38 ims_smd_max = np.max(ims_smd)
39
40 txyb_all = np.array([]).reshape(0, 4)
41 for i in range(frame_1st - 1, frame_end):
42 tmp = np.copy(ims_smd[i, :, :])
43 tmp[tmp < th * ims_smd_max / 100] = 0
44 coords = peak_local_max(tmp, min_distance=1)
45 idx_to_del = np.where((coords[:, 0] <= bd) | (coords[:, 0] >= tmp.shape[0] - bd) |
46 (coords[:, 1] <= bd) | (coords[:, 1] >= tmp.shape[1] - bd))
47 coords = np.delete(coords, idx_to_del[0], axis=0)
48 xys = coords[:, ::-1]
49
50 if typ_br == 'smoothed':
51 intens = getbr(xys, ims_smd[i, :, :], 0, 1)
52 elif typ_br == 'robust':
53 intens = getbr(xys, ims_ori[i], 1, 4)
54 else:
55 intens = getbr(xys, ims_ori[i], 0, 1)
56
57 txyb = np.concatenate(((i + 1) * np.ones((xys.shape[0], 1)), xys, intens), axis=1)
58 txyb_all = np.concatenate((txyb_all, txyb), axis=0)
59
60 df = pd.DataFrame()
61 df['FRAME'] = txyb_all[:, 0].astype(int)
62 df['POS_X'] = txyb_all[:, 1].astype(int)
63 df['POS_Y'] = txyb_all[:, 2].astype(int)
64 df['INTENSITY'] = txyb_all[:, 3]
65 df.to_csv(fn_out, index=False, float_format='%.2f', sep="\t")
66
67
68 if __name__ == "__main__":
69 parser = argparse.ArgumentParser(description="Spot detection based on local maxima")
70 parser.add_argument("fn_in", help="Name of input image sequence (stack)")
71 parser.add_argument("fn_out", help="Name of output file to save the coordinates and intensities of detected spots")
72 parser.add_argument("frame_1st", type=int, help="Index for the starting frame to detect spots (1 for first frame of the stack)")
73 parser.add_argument("frame_end", type=int, help="Index for the last frame to detect spots (0 for the last frame of the stack)")
74 parser.add_argument("typ_intens", help="smoothed or robust (for measuring the intensities of spots)")
75 parser.add_argument("thres", type=float, help="Percentage of the global maximal intensity for thresholding candidate spots")
76 parser.add_argument("ssig", type=float, help="Sigma of the Gaussian filter for noise suppression")
77 parser.add_argument("bndy", type=int, help="Number of pixels (Spots close to image boundaries will be ignored)")
78 args = parser.parse_args()
79 spot_detection(args.fn_in,
80 args.fn_out,
81 frame_1st=args.frame_1st,
82 frame_end=args.frame_end,
83 typ_br=args.typ_intens,
84 th=args.thres,
85 ssig=args.ssig,
86 bd=args.bndy)