# HG changeset patch # User imgteam # Date 1626992987 0 # Node ID 04e692ee53a83e29f7bd2580908ceef499dfb37f "planemo upload for repository https://github.com/BMCV/galaxy-image-analysis/tree/master/tools/points_association_nn/ commit db4c2a87a21f32e5d12d11e68f32773bfc06fcfd" diff -r 000000000000 -r 04e692ee53a8 points_association_nn.py --- /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) diff -r 000000000000 -r 04e692ee53a8 points_association_nn.xml --- /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 @@ + + in consecutive frames (slices) using the nearest neighbor algorithm + + numpy + openpyxl + pandas + scikit-image + + + + + + + + + + + + + + + + + + + + + + + + **What it does** + + This tool associates points in consecutive frames (slices) using the nearest neighbor algorithm. + + diff -r 000000000000 -r 04e692ee53a8 test-data/spots_detected.tsv --- /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 diff -r 000000000000 -r 04e692ee53a8 test-data/spots_linked.xlsx Binary file test-data/spots_linked.xlsx has changed