Mercurial > repos > imgteam > points2labelimage

comparison points2label.py @ 3:de611b3b5ae8 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/points2labelimage/ commit 6fc9ab8db9ef72ac7ded30d7373768feeae9390d
author imgteam
date Fri, 27 Sep 2024 17:41:21 +0000
parents 30ca5d5d03ec
children
comparison
equal deleted inserted replaced
2:30ca5d5d03ec 3:de611b3b5ae8
1 import argparse 1 import argparse
2 import sys 2 import os
3 import warnings 3 import warnings
4 4
5 import giatools.pandas
5 import numpy as np 6 import numpy as np
6 import pandas as pd 7 import pandas as pd
8 import scipy.ndimage as ndi
7 import skimage.io 9 import skimage.io
10 import skimage.segmentation
8 11
9 12
10 def points2label(labels, shape, output_file=None, has_header=False, is_TSV=True): 13 def rasterize(point_file, out_file, shape, has_header=False, swap_xy=False, bg_value=0, fg_value=None):
11 labelimg = np.zeros([shape[0], shape[1]], dtype=np.int32)
12 14
13 if is_TSV: 15 img = np.full(shape, dtype=np.uint16, fill_value=bg_value)
16 if os.path.exists(point_file) and os.path.getsize(point_file) > 0:
17
18 # Read the tabular file with information from the header
14 if has_header: 19 if has_header:
15 df = pd.read_csv(labels, sep='\t', skiprows=1, header=None) 20 df = pd.read_csv(point_file, delimiter='\t')
21
22 pos_x_column = giatools.pandas.find_column(df, ['pos_x', 'POS_X'])
23 pos_y_column = giatools.pandas.find_column(df, ['pos_y', 'POS_Y'])
24 pos_x_list = df[pos_x_column].round().astype(int)
25 pos_y_list = df[pos_y_column].round().astype(int)
26 assert len(pos_x_list) == len(pos_y_list)
27
28 try:
29 radius_column = giatools.pandas.find_column(df, ['radius', 'RADIUS'])
30 radius_list = df[radius_column]
31 assert len(pos_x_list) == len(radius_list)
32 except KeyError:
33 radius_list = [0] * len(pos_x_list)
34
35 try:
36 label_column = giatools.pandas.find_column(df, ['label', 'LABEL'])
37 label_list = df[label_column]
38 assert len(pos_x_list) == len(label_list)
39 except KeyError:
40 label_list = list(range(1, len(pos_x_list) + 1))
41
42 # Read the tabular file without header
16 else: 43 else:
17 df = pd.read_csv(labels, sep='\t', header=None) 44 df = pd.read_csv(point_file, header=None, delimiter='\t')
45 pos_x_list = df[0].round().astype(int)
46 pos_y_list = df[1].round().astype(int)
47 assert len(pos_x_list) == len(pos_y_list)
48 radius_list = [0] * len(pos_x_list)
49 label_list = list(range(1, len(pos_x_list) + 1))
50
51 # Optionally swap the coordinates
52 if swap_xy:
53 pos_x_list, pos_y_list = pos_y_list, pos_x_list
54
55 # Perform the rasterization
56 for y, x, radius, label in zip(pos_y_list, pos_x_list, radius_list, label_list):
57 if fg_value is not None:
58 label = fg_value
59
60 if y < 0 or x < 0 or y >= shape[0] or x >= shape[1]:
61 raise IndexError(f'The point x={x}, y={y} exceeds the bounds of the image (width: {shape[1]}, height: {shape[0]})')
62
63 # Rasterize circle and distribute overlapping image area
64 if radius > 0:
65 mask = np.ones(shape, dtype=bool)
66 mask[y, x] = False
67 mask = (ndi.distance_transform_edt(mask) <= radius)
68
69 # Compute the overlap (pretend there is none if the rasterization is binary)
70 if fg_value is None:
71 overlap = np.logical_and(img > 0, mask)
72 else:
73 overlap = np.zeros(shape, dtype=bool)
74
75 # Rasterize the part of the circle which is disjoint from other foreground.
76 #
77 # In the current implementation, the result depends on the order of the rasterized circles if somewhere
78 # more than two circles overlap. This is probably negligable for most applications. To achieve results
79 # that are invariant to the order, first all circles would need to be rasterized independently, and
80 # then blended together. This, however, would either strongly increase the memory consumption, or
81 # require a more complex implementation which exploits the sparsity of the rasterized masks.
82 #
83 disjoint_mask = np.logical_xor(mask, overlap)
84 if disjoint_mask.any():
85 img[disjoint_mask] = label
86
87 # Distribute the remaining part of the circle
88 if overlap.any():
89 dist = ndi.distance_transform_edt(overlap)
90 foreground = (img > 0)
91 img[overlap] = 0
92 img = skimage.segmentation.watershed(dist, img, mask=foreground)
93
94 # Rasterize point (there is no overlapping area to be distributed)
95 else:
96 img[y, x] = label
97
18 else: 98 else:
19 if has_header: 99 raise Exception("{} is empty or does not exist.".format(point_file)) # appropriate built-in error?
20 df = pd.read_csv(labels, skiprows=1, header=None)
21 else:
22 df = pd.read_csv(labels, header=None)
23 100
24 for i in range(0, len(df)): 101 with warnings.catch_warnings():
25 a_row = df.iloc[i] 102 warnings.simplefilter("ignore")
26 labelimg[a_row[0], a_row[1]] = i + 1 103 skimage.io.imsave(out_file, img, plugin='tifffile') # otherwise we get problems with the .dat extension
27
28 if output_file is not None:
29 with warnings.catch_warnings():
30 warnings.simplefilter("ignore")
31 skimage.io.imsave(output_file, labelimg, plugin='tifffile')
32 else:
33 return labelimg
34 104
35 105
36 if __name__ == "__main__": 106 if __name__ == "__main__":
37 parser = argparse.ArgumentParser() 107 parser = argparse.ArgumentParser()
38 parser.add_argument('label_file', type=argparse.FileType('r'), default=sys.stdin, help='label file') 108 parser.add_argument('point_file', type=argparse.FileType('r'), help='point file')
39 parser.add_argument('out_file', type=argparse.FileType('w'), default=sys.stdin, help='out file') 109 parser.add_argument('out_file', type=str, help='out file (TIFF)')
40 parser.add_argument('org_file', type=argparse.FileType('r'), default=sys.stdin, help='input original file') 110 parser.add_argument('shapex', type=int, help='shapex')
41 parser.add_argument('--has_header', dest='has_header', type=bool, default=False, help='label file has header') 111 parser.add_argument('shapey', type=int, help='shapey')
42 parser.add_argument('--is_tsv', dest='is_tsv', type=bool, default=True, help='label file is TSV') 112 parser.add_argument('--has_header', dest='has_header', default=False, help='set True if point file has header')
113 parser.add_argument('--swap_xy', dest='swap_xy', default=False, help='Swap X and Y coordinates')
114 parser.add_argument('--binary', dest='binary', default=False, help='Produce binary image')
115
43 args = parser.parse_args() 116 args = parser.parse_args()
44 117
45 original_shape = skimage.io.imread(args.org_file.name, plugin='tifffile').shape 118 rasterize(
46 119 args.point_file.name,
47 points2label(args.label_file.name, original_shape, args.out_file.name, args.has_header, args.is_tsv) 120 args.out_file,
121 (args.shapey, args.shapex),
122 has_header=args.has_header,
123 swap_xy=args.swap_xy,
124 fg_value=0xffff if args.binary else None,
125 )