Mercurial > repos > imgteam > points_association_nn

changeset 0:04e692ee53a8 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/points_association_nn/ commit db4c2a87a21f32e5d12d11e68f32773bfc06fcfd"
author imgteam
date Thu, 22 Jul 2021 22:29:47 +0000
parents
children fd4293bee0dc
files points_association_nn.py points_association_nn.xml test-data/spots_detected.tsv test-data/spots_linked.xlsx
diffstat 4 files changed, 672 insertions(+), 0 deletions(-) [+]
line wrap: on
line diff
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/points_association_nn.py	Thu Jul 22 22:29:47 2021 +0000
@@ -0,0 +1,168 @@
+"""
+Copyright 2021 Biomedical Computer Vision Group, Heidelberg University.
+Author: Qi Gao (qi.gao@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 numpy as np
+import openpyxl  # noqa: F401
+import pandas as pd
+import skimage.util
+
+
+def disk_mask(imsz, ir, ic, nbpx):
+    ys, xs = np.ogrid[-nbpx:nbpx + 1, -nbpx:nbpx + 1]
+    se = xs ** 2 + ys ** 2 <= nbpx ** 2
+    mask = np.zeros(imsz, dtype=int)
+    if ir - nbpx < 0 or ic - nbpx < 0 or ir + nbpx + 1 > imsz[0] or ic + nbpx + 1 > imsz[1]:
+        mask = skimage.util.pad(mask, nbpx)
+        mask[ir:ir + 2 * nbpx + 1, ic:ic + 2 * nbpx + 1] = se
+        mask = skimage.util.crop(mask, nbpx)
+    else:
+        mask[ir - nbpx:ir + nbpx + 1, ic - nbpx:ic + nbpx + 1] = se
+    return mask
+
+
+def find_nn(cim, icy, icx, nim, nbpx):
+    mask = disk_mask(cim.shape, icy, icx, nbpx)
+    iys_nim, ixs_nim = np.where(nim * mask)
+    if iys_nim.size == 0:
+        return np.NaN, np.NaN
+
+    d2 = (icy - iys_nim) ** 2 + (icx - ixs_nim) ** 2
+    I1 = np.argsort(d2)
+    iy_nim = iys_nim[I1[0]]
+    ix_nim = ixs_nim[I1[0]]
+
+    mask = disk_mask(cim.shape, iy_nim, ix_nim, nbpx)
+    iys_cim, ixs_cim = np.where(cim * mask)
+    d2 = (iy_nim - iys_cim) ** 2 + (ix_nim - ixs_cim) ** 2
+    I2 = np.argsort(d2)
+    if not iys_cim[I2[0]] == icy or not ixs_cim[I2[0]] == icx:
+        return np.NaN, np.NaN
+
+    return iy_nim, ix_nim
+
+
+def points_linking(fn_in, fn_out, nbpx=6, th=25, minlen=50):
+    data = pd.read_csv(fn_in, delimiter="\t")
+    all_data = np.array(data)
+    assert all_data.shape[1] in [3, 4], 'unknow collum(s) in input data!'
+
+    coords = all_data[:, :3].astype('int64')
+
+    frame_1st = np.min(coords[:, 0])
+    frame_end = np.max(coords[:, 0])
+    assert set([i for i in range(frame_1st, frame_end + 1)]).issubset(set(coords[:, 0].tolist())), "spots missing at some time point!"
+
+    nSlices = frame_end
+    stack_h = np.max(coords[:, 2]) + nbpx
+    stack_w = np.max(coords[:, 1]) + nbpx
+    stack = np.zeros((stack_h, stack_w, nSlices), dtype='int8')
+    stack_r = np.zeros((stack_h, stack_w, nSlices), dtype='float64')
+
+    for i in range(all_data.shape[0]):
+        iyxz = tuple(coords[i, ::-1] - 1)
+        stack[iyxz] = 1
+        stack_r[iyxz] = all_data[i, -1]
+
+    tracks_all = np.array([], dtype=float).reshape(0, nSlices, 4)
+    maxv = np.max(stack_r)
+    br_max = maxv
+    idx_max = np.argmax(stack_r)
+    while 1:
+        iyxz = np.unravel_index(idx_max, stack.shape)
+
+        spot_br = np.empty((nSlices, 1))
+        track = np.empty((nSlices, 3))
+        for i in range(nSlices):
+            spot_br[i] = np.NaN
+            track[i, :] = np.array((np.NaN, np.NaN, np.NaN))
+
+        spot_br[iyxz[2]] = maxv
+        track[iyxz[2], :] = np.array(iyxz[::-1]) + 1
+
+        # forward
+        icy = iyxz[0]
+        icx = iyxz[1]
+        for inz in range(iyxz[2] + 1, nSlices):
+            iny, inx = find_nn(stack[:, :, inz - 1], icy, icx, stack[:, :, inz], nbpx)
+            if np.isnan(iny) and not inz == nSlices - 1:
+                iny, inx = find_nn(stack[:, :, inz - 1], icy, icx, stack[:, :, inz + 1], nbpx)
+                if np.isnan(iny):
+                    break
+                else:
+                    iny = icy
+                    inx = icx
+                    stack[iny, inx, inz] = 1
+                    stack_r[iny, inx, inz] = stack_r[iny, inx, inz - 1]
+            elif np.isnan(iny) and inz == nSlices - 1:
+                break
+
+            track[inz, :] = np.array((inz, inx, iny)) + 1
+            spot_br[inz] = stack_r[iny, inx, inz]
+            icy = iny
+            icx = inx
+
+        # backward
+        icy = iyxz[0]
+        icx = iyxz[1]
+        for inz in range(iyxz[2] - 1, -1, -1):
+            iny, inx = find_nn(stack[:, :, inz + 1], icy, icx, stack[:, :, inz], nbpx)
+            if np.isnan(iny) and not inz == 0:
+                iny, inx = find_nn(stack[:, :, inz + 1], icy, icx, stack[:, :, inz - 1], nbpx)
+                if np.isnan(iny):
+                    break
+                else:
+                    iny = icy
+                    inx = icx
+                    stack[iny, inx, inz] = 1
+                    stack_r[iny, inx, inz] = stack_r[iny, inx, inz + 1]
+            elif np.isnan(iny) and inz == 0:
+                break
+
+            track[inz, :] = np.array((inz, inx, iny)) + 1
+            spot_br[inz] = stack_r[iny, inx, inz]
+            icy = iny
+            icx = inx
+
+        for iz in range(nSlices):
+            if not np.isnan(track[iz, 0]):
+                stack[track[iz, 2].astype(int) - 1, track[iz, 1].astype(int) - 1, iz] = 0
+                stack_r[track[iz, 2].astype(int) - 1, track[iz, 1].astype(int) - 1, iz] = 0
+
+        # discard short trajectories
+        if np.count_nonzero(~np.isnan(spot_br)) > minlen * (frame_end - frame_1st) / 100:
+            tmp = np.concatenate((track, spot_br), axis=1)
+            tracks_all = np.concatenate((tracks_all, tmp.reshape(1, -1, 4)), axis=0)
+
+        maxv = np.max(stack_r)
+        idx_max = np.argmax(stack_r)
+        if maxv < th * br_max / 100:
+            break
+
+    with pd.ExcelWriter(fn_out, engine="openpyxl") as writer:
+        for i in range(tracks_all.shape[0]):
+            df = pd.DataFrame()
+            df['FRAME'] = tracks_all[i, :, 0]
+            df['POS_X'] = tracks_all[i, :, 1]
+            df['POS_Y'] = tracks_all[i, :, 2]
+            df['INTENSITY'] = tracks_all[i, :, 3]
+            df.to_excel(writer, sheet_name='spot%s' % (i + 1), index=False, float_format='%.2f')
+        writer.save()
+
+
+if __name__ == "__main__":
+    parser = argparse.ArgumentParser(description="Association of points in consecutive frames using the nearest neighbor algorithm")
+    parser.add_argument("fn_in", help="Name of input file (tsv tabular)")
+    parser.add_argument("fn_out", help="Name of output file (xlsx)")
+    parser.add_argument("nbpx", type=int, help="Neighborhood size in pixel")
+    parser.add_argument("thres", type=float, help="Percentage of the global maximal intensity for thresholding some event")
+    parser.add_argument("minlen", type=float, help="Minimum length of tracks (percentage of senquence length)")
+    args = parser.parse_args()
+    points_linking(args.fn_in, args.fn_out, args.nbpx, args.thres, args.minlen)
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/points_association_nn.xml	Thu Jul 22 22:29:47 2021 +0000
@@ -0,0 +1,42 @@
+<tool id="ip_points_association_nn" name="Association of points" version="0.0.1" profile="20.05"> 
+    <description>in consecutive frames (slices) using the nearest neighbor algorithm</description>
+    <requirements>
+        <requirement type="package" version="1.20.2">numpy</requirement>
+        <requirement type="package" version="3.0.7">openpyxl</requirement>
+        <requirement type="package" version="1.2.4">pandas</requirement>
+        <requirement type="package" version="0.18.1">scikit-image</requirement>
+    </requirements>
+    <command>
+    <![CDATA[
+    python '$__tool_directory__/points_association_nn.py'
+         '$fn_in'
+         ./output.xlsx
+         '$nbpx'
+         '$thres'
+         '$minlen'
+    ]]>
+    </command>
+    <inputs>
+        <param name="fn_in" type="data" format="tabular" label="Name of input file (tsv tabular)" />
+        <param name="nbpx" type="integer" value="6" label="Neighborhood size in pixel" />
+        <param name="thres" type="float" value="25" label="Percentage (%) of the global maximal intensity for thresholding some event" />
+        <param name="minlen" type="float" value="50" label="Minimum length of tracks (% of sequence length)" />
+    </inputs>
+    <outputs>
+        <data format="xlsx" name="fn_out" from_work_dir="output.xlsx" />
+    </outputs>
+    <tests>
+        <test>
+            <param name="fn_in" value="spots_detected.tsv"/>
+            <param name="nbpx" value="6"/>
+            <param name="thres" value="25"/>
+            <param name="minlen" value="50"/>
+            <output name="fn_out" value="spots_linked.xlsx" ftype="xlsx" compare="sim_size"/>
+        </test>
+    </tests>
+    <help>
+    **What it does**
+
+    This tool associates points in consecutive frames (slices) using the nearest neighbor algorithm.
+    </help>
+</tool>
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/test-data/spots_detected.tsv	Thu Jul 22 22:29:47 2021 +0000
@@ -0,0 +1,462 @@
+FRAME	POS_X	POS_Y	INTENSITY
+1	126	202	55939.36
+1	218	227	41590.82
+1	132	201	41849.38
+1	120	199	45491.74
+1	113	177	27103.64
+1	95	135	33206.37
+1	130	209	36872.91
+1	128	170	30994.56
+1	315	240	38767.19
+1	123	195	32988.47
+1	117	181	40031.76
+1	118	184	36812.77
+1	127	181	29651.66
+1	134	192	28531.76
+1	137	189	35093.95
+1	195	223	30206.57
+1	78	87	30367.18
+1	146	195	26779.82
+1	175	225	28515.74
+1	213	227	30812.06
+1	114	158	31229.97
+1	171	223	26982.37
+1	87	117	27143.89
+1	129	194	29475.93
+1	115	192	28510.64
+1	248	230	22615.54
+1	154	208	26394.30
+1	109	153	26268.78
+1	328	234	26387.65
+1	138	179	20215.40
+1	159	212	24400.22
+1	335	233	21698.80
+1	387	198	24933.29
+1	121	154	18710.45
+1	308	239	22328.27
+1	258	230	18680.93
+1	340	229	15937.72
+1	78	78	18151.49
+1	81	80	24704.83
+1	117	161	24293.65
+1	120	164	24242.78
+1	160	215	23759.44
+1	86	141	10949.28
+1	222	230	22277.61
+1	370	215	18379.39
+1	109	199	18850.01
+1	83	109	15871.12
+1	117	155	15703.14
+1	100	182	17281.91
+1	168	227	17431.67
+1	392	197	15897.67
+1	184	227	11346.34
+1	294	237	10236.84
+1	377	210	15642.67
+1	152	226	9851.10
+1	93	145	12584.04
+1	147	212	14735.83
+1	80	101	11021.47
+1	394	193	14063.18
+1	105	139	12039.57
+1	301	237	11466.69
+1	273	233	8421.28
+1	111	170	14406.46
+1	237	228	9466.36
+1	106	206	14522.55
+1	83	75	11115.23
+1	290	238	12636.07
+1	205	229	12808.78
+1	381	207	11531.86
+1	92	196	9024.83
+1	185	222	11196.17
+1	180	221	8339.93
+1	146	185	7207.70
+1	142	209	11730.52
+1	88	101	10415.08
+1	235	231	11264.22
+1	85	103	5911.48
+1	104	183	8101.94
+1	85	148	7290.63
+1	145	189	9776.95
+1	361	219	7009.70
+1	282	236	8853.80
+1	121	169	7733.35
+1	107	180	8135.17
+1	266	235	7470.55
+1	79	96	6920.43
+1	152	219	5623.45
+1	149	220	5889.31
+2	126	202	59098.11
+2	126	205	53383.42
+2	218	227	46745.07
+2	120	199	43055.88
+2	116	180	45577.81
+2	127	181	42874.19
+2	130	209	41387.22
+2	123	195	33569.00
+2	113	177	31339.11
+2	128	170	35552.32
+2	135	191	30546.65
+2	95	136	34169.11
+2	315	238	31430.32
+2	314	240	30015.35
+2	116	201	26928.89
+2	196	223	29118.84
+2	175	225	29053.82
+2	213	227	34704.93
+2	115	159	29750.71
+2	146	195	27597.86
+2	328	234	28969.08
+2	113	191	25325.73
+2	79	87	27998.19
+2	340	229	19788.22
+2	154	207	25348.79
+2	81	81	25643.11
+2	172	224	25762.85
+2	87	117	24830.46
+2	169	222	19415.52
+2	325	232	23879.48
+2	248	230	23270.45
+2	110	153	24396.90
+2	387	198	23947.67
+2	87	143	16954.62
+2	118	162	26061.81
+2	78	79	23543.41
+2	258	230	16996.33
+2	121	165	22822.04
+2	387	201	25260.27
+2	335	233	19553.33
+2	223	230	22795.66
+2	159	212	20956.68
+2	370	216	19989.36
+2	320	233	13806.24
+2	308	239	17964.97
+2	105	139	18000.56
+2	109	200	17896.86
+2	83	109	16421.63
+2	105	151	17984.13
+2	137	179	12811.40
+2	110	169	14111.65
+2	228	231	18325.92
+2	392	196	14572.82
+2	153	228	14935.83
+2	122	154	11708.73
+2	83	75	13067.82
+2	153	213	16283.67
+2	295	237	9677.50
+2	377	210	13457.71
+2	273	233	9148.80
+2	205	229	13083.06
+2	381	206	10603.41
+2	80	101	9485.05
+2	152	202	15278.04
+2	185	226	11093.24
+2	161	216	15543.22
+2	146	212	14364.79
+2	235	231	12509.71
+2	290	238	13411.96
+2	93	145	9119.49
+2	142	210	13054.23
+2	239	232	13043.98
+2	92	196	10749.96
+2	185	221	9922.40
+2	299	238	13577.34
+2	278	235	11194.38
+2	100	182	11485.41
+2	145	189	11315.24
+2	237	228	11272.61
+2	107	189	8051.49
+2	87	101	8814.45
+2	146	186	7585.02
+2	106	207	8720.64
+2	109	207	7676.81
+2	402	188	7813.91
+2	88	108	7921.85
+2	360	220	7726.74
+2	86	148	7184.91
+2	106	180	8401.49
+2	95	173	4093.76
+2	266	235	5678.52
+2	122	172	5731.76
+2	127	218	4787.41
+3	126	206	56658.14
+3	126	202	54234.74
+3	218	226	35022.37
+3	119	200	40576.02
+3	133	201	44149.96
+3	115	178	43820.79
+3	124	196	41934.96
+3	117	180	43996.85
+3	130	209	42016.60
+3	136	190	39500.63
+3	316	237	31963.21
+3	97	137	36182.78
+3	127	181	31990.16
+3	119	184	34251.83
+3	114	158	29844.61
+3	146	196	37064.84
+3	128	186	34532.76
+3	176	225	30746.00
+3	196	223	29807.10
+3	213	226	25523.49
+3	128	171	28664.10
+3	152	204	24687.37
+3	329	234	28277.06
+3	78	87	26725.30
+3	81	87	30473.18
+3	172	225	26643.89
+3	120	163	26474.56
+3	387	197	21111.91
+3	86	117	20367.28
+3	335	233	20895.73
+3	324	232	23086.43
+3	339	229	15275.48
+3	129	194	24033.17
+3	123	185	24073.76
+3	77	79	19152.47
+3	169	222	17195.61
+3	132	175	25123.86
+3	81	81	24895.65
+3	109	153	21172.27
+3	117	161	24299.24
+3	307	240	19257.85
+3	370	215	18554.16
+3	222	230	22017.20
+3	104	139	20194.29
+3	248	230	23396.01
+3	257	230	17273.50
+3	116	154	15659.69
+3	114	190	21449.57
+3	158	211	18454.95
+3	83	109	17776.43
+3	84	75	15659.83
+3	158	215	19800.41
+3	109	199	15805.24
+3	101	140	20715.26
+3	143	210	20769.18
+3	88	144	20132.54
+3	93	145	12670.17
+3	228	231	19053.86
+3	393	195	13463.33
+3	80	101	13305.58
+3	105	151	15913.95
+3	137	179	15719.37
+3	167	226	17465.42
+3	293	238	15459.81
+3	377	209	10830.81
+3	185	222	13711.26
+3	108	189	16955.73
+3	185	226	12648.79
+3	153	227	14010.45
+3	101	186	14180.61
+3	298	238	16257.78
+3	121	189	14906.88
+3	151	210	12848.29
+3	274	233	11026.17
+3	93	197	14149.59
+3	234	231	13284.28
+3	236	228	8838.32
+3	106	207	10942.68
+3	381	206	8704.24
+3	238	231	13112.03
+3	109	169	9573.93
+3	88	101	9934.59
+3	104	183	7882.07
+3	85	103	8202.08
+3	361	219	7171.66
+3	278	235	10260.85
+3	268	235	8025.01
+3	146	188	9096.57
+3	283	236	9090.54
+3	399	190	6873.99
+3	145	185	6134.20
+3	105	179	6141.86
+3	97	176	6082.70
+3	110	207	4505.63
+3	127	219	5028.42
+3	94	171	6146.80
+4	126	202	55939.36
+4	218	227	41590.82
+4	132	201	41849.38
+4	120	199	45491.74
+4	113	177	27103.64
+4	95	135	33206.37
+4	130	209	36872.91
+4	128	170	30994.56
+4	315	240	38767.19
+4	123	195	32988.47
+4	117	181	40031.76
+4	118	184	36812.77
+4	127	181	29651.66
+4	134	192	28531.76
+4	137	189	35093.95
+4	195	223	30206.57
+4	78	87	30367.18
+4	146	195	26779.82
+4	175	225	28515.74
+4	213	227	30812.06
+4	114	158	31229.97
+4	171	223	26982.37
+4	87	117	27143.89
+4	129	194	29475.93
+4	115	192	28510.64
+4	248	230	22615.54
+4	154	208	26394.30
+4	109	153	26268.78
+4	328	234	26387.65
+4	138	179	20215.40
+4	159	212	24400.22
+4	335	233	21698.80
+4	387	198	24933.29
+4	121	154	18710.45
+4	308	239	22328.27
+4	258	230	18680.93
+4	340	229	15937.72
+4	78	78	18151.49
+4	81	80	24704.83
+4	117	161	24293.65
+4	120	164	24242.78
+4	160	215	23759.44
+4	86	141	10949.28
+4	222	230	22277.61
+4	370	215	18379.39
+4	109	199	18850.01
+4	83	109	15871.12
+4	117	155	15703.14
+4	100	182	17281.91
+4	168	227	17431.67
+4	392	197	15897.67
+4	184	227	11346.34
+4	294	237	10236.84
+4	377	210	15642.67
+4	152	226	9851.10
+4	93	145	12584.04
+4	147	212	14735.83
+4	80	101	11021.47
+4	394	193	14063.18
+4	105	139	12039.57
+4	301	237	11466.69
+4	273	233	8421.28
+4	111	170	14406.46
+4	237	228	9466.36
+4	106	206	14522.55
+4	83	75	11115.23
+4	290	238	12636.07
+4	205	229	12808.78
+4	381	207	11531.86
+4	92	196	9024.83
+4	185	222	11196.17
+4	180	221	8339.93
+4	146	185	7207.70
+4	142	209	11730.52
+4	88	101	10415.08
+4	235	231	11264.22
+4	85	103	5911.48
+4	104	183	8101.94
+4	85	148	7290.63
+4	145	189	9776.95
+4	361	219	7009.70
+4	282	236	8853.80
+4	121	169	7733.35
+4	107	180	8135.17
+4	266	235	7470.55
+4	79	96	6920.43
+4	152	219	5623.45
+4	149	220	5889.31
+5	126	206	56658.14
+5	126	202	54234.74
+5	218	226	35022.37
+5	119	200	40576.02
+5	133	201	44149.96
+5	115	178	43820.79
+5	124	196	41934.96
+5	117	180	43996.85
+5	130	209	42016.60
+5	136	190	39500.63
+5	316	237	31963.21
+5	97	137	36182.78
+5	127	181	31990.16
+5	119	184	34251.83
+5	114	158	29844.61
+5	146	196	37064.84
+5	128	186	34532.76
+5	176	225	30746.00
+5	196	223	29807.10
+5	213	226	25523.49
+5	128	171	28664.10
+5	152	204	24687.37
+5	329	234	28277.06
+5	78	87	26725.30
+5	81	87	30473.18
+5	172	225	26643.89
+5	120	163	26474.56
+5	387	197	21111.91
+5	86	117	20367.28
+5	335	233	20895.73
+5	324	232	23086.43
+5	339	229	15275.48
+5	129	194	24033.17
+5	123	185	24073.76
+5	77	79	19152.47
+5	169	222	17195.61
+5	132	175	25123.86
+5	81	81	24895.65
+5	109	153	21172.27
+5	117	161	24299.24
+5	307	240	19257.85
+5	370	215	18554.16
+5	222	230	22017.20
+5	104	139	20194.29
+5	248	230	23396.01
+5	257	230	17273.50
+5	116	154	15659.69
+5	114	190	21449.57
+5	158	211	18454.95
+5	83	109	17776.43
+5	84	75	15659.83
+5	158	215	19800.41
+5	109	199	15805.24
+5	101	140	20715.26
+5	143	210	20769.18
+5	88	144	20132.54
+5	93	145	12670.17
+5	228	231	19053.86
+5	393	195	13463.33
+5	80	101	13305.58
+5	105	151	15913.95
+5	137	179	15719.37
+5	167	226	17465.42
+5	293	238	15459.81
+5	377	209	10830.81
+5	185	222	13711.26
+5	108	189	16955.73
+5	185	226	12648.79
+5	153	227	14010.45
+5	101	186	14180.61
+5	298	238	16257.78
+5	121	189	14906.88
+5	151	210	12848.29
+5	274	233	11026.17
+5	93	197	14149.59
+5	234	231	13284.28
+5	236	228	8838.32
+5	106	207	10942.68
+5	381	206	8704.24
+5	238	231	13112.03
+5	109	169	9573.93
+5	88	101	9934.59
+5	104	183	7882.07
+5	85	103	8202.08
+5	361	219	7171.66
+5	278	235	10260.85
+5	268	235	8025.01
+5	146	188	9096.57
+5	283	236	9090.54
+5	399	190	6873.99
+5	145	185	6134.20
+5	105	179	6141.86
+5	97	176	6082.70
+5	110	207	4505.63
+5	127	219	5028.42
+5	94	171	6146.80
Binary file test-data/spots_linked.xlsx has changed