--- a/points2binaryimage.py	Wed Jun 26 08:35:56 2024 +0000
+++ b/points2binaryimage.py	Wed Sep 25 08:29:59 2024 +0000
@@ -1,36 +1,75 @@
 import argparse
 import os
 import warnings
+from typing import List
 
 import numpy as np
 import pandas as pd
+import scipy.ndimage as ndi
 import skimage.io
 
 
-def points2binaryimage(point_file, out_file, shape=[500, 500], has_header=False, invert_xy=False):
+def find_column(df: pd.DataFrame, candidates: List[str]) -> str:
+    """
+    Returns the column name present in `df` and the list of `candidates`.
 
-    img = np.zeros(shape, dtype=np.int16)
+    Raises:
+        KeyError: If there is no candidate column name present in `df`, or more than one.
+    """
+    intersection = frozenset(df.columns) & frozenset(candidates)
+    if len(intersection) == 0:
+        raise KeyError(f'No such column: {", ".join(candidates)}')
+    elif len(intersection) > 1:
+        raise KeyError(f'The column names {", ".join(intersection)} are ambiguous')
+    else:
+        return next(iter(intersection))
+
+
+def points2binaryimage(point_file, out_file, shape, has_header=False, swap_xy=False, bg_value=0, fg_value=0xffff):
+
+    img = np.full(shape, dtype=np.uint16, fill_value=bg_value)
     if os.path.exists(point_file) and os.path.getsize(point_file) > 0:
-        if has_header:
-            df = pd.read_csv(point_file, skiprows=1, header=None, delimiter="\t")
-        else:
-            df = pd.read_csv(point_file, header=None, delimiter="\t")
 
-        for i in range(0, len(df)):
-            a_row = df.iloc[i]
-            if int(a_row[0]) < 0 or int(a_row[1]) < 0:
-                raise IndexError("Point {},{} is out of image with bounds {},{}.".format(int(a_row[0]), int(a_row[1]), shape[0], shape[1]))
+        # Read the tabular file with information from the header
+        if has_header:
+            df = pd.read_csv(point_file, delimiter='\t')
+            pos_x_column = find_column(df, ['pos_x', 'POS_X'])
+            pos_y_column = find_column(df, ['pos_y', 'POS_Y'])
+            pos_x_list = df[pos_x_column].round().astype(int)
+            pos_y_list = df[pos_y_column].round().astype(int)
+            assert len(pos_x_list) == len(pos_y_list)
+            try:
+                radius_column = find_column(df, ['radius', 'RADIUS'])
+                radius_list = df[radius_column]
+            except KeyError:
+                radius_list = [0] * len(pos_x_list)
 
-            if invert_xy:
-                if img.shape[0] <= int(a_row[0]) or img.shape[1] <= int(a_row[1]):
-                    raise IndexError("Point {},{} is out of image with bounds {},{}.".format(int(a_row[0]), int(a_row[1]), shape[0], shape[1]))
-                else:
-                    img[int(a_row[1]), int(a_row[0])] = 32767
+        # Read the tabular file without header
+        else:
+            df = pd.read_csv(point_file, header=None, delimiter='\t')
+            pos_x_list = df[0].round().astype(int)
+            pos_y_list = df[1].round().astype(int)
+            assert len(pos_x_list) == len(pos_y_list)
+            radius_list = [0] * len(pos_x_list)
+
+        # Optionally swap the coordinates
+        if swap_xy:
+            pos_x_list, pos_y_list = pos_y_list, pos_x_list
+
+        # Perform the rasterization
+        for y, x, radius in zip(pos_y_list, pos_x_list, radius_list):
+
+            if y < 0 or x < 0 or y >= shape[0] or x >= shape[1]:
+                raise IndexError(f'The point x={x}, y={y} exceeds the bounds of the image (width: {shape[1]}, height: {shape[0]})')
+
+            if radius > 0:
+                mask = np.ones(shape, dtype=bool)
+                mask[y, x] = False
+                mask = (ndi.distance_transform_edt(mask) <= radius)
+                img[mask] = fg_value
             else:
-                if img.shape[0] <= int(a_row[1]) or img.shape[1] <= int(a_row[0]):
-                    raise IndexError("Point {},{} is out of image with bounds {},{}.".format(int(a_row[1]), int(a_row[0]), shape[0], shape[1]))
-                else:
-                    img[int(a_row[0]), int(a_row[1])] = 32767
+                img[y, x] = fg_value
+
     else:
         raise Exception("{} is empty or does not exist.".format(point_file))  # appropriate built-in error?
 
@@ -41,14 +80,14 @@
 
 if __name__ == "__main__":
     parser = argparse.ArgumentParser()
-    parser.add_argument('point_file', type=argparse.FileType('r'), help='label file')
+    parser.add_argument('point_file', type=argparse.FileType('r'), help='point file')
     parser.add_argument('out_file', type=str, help='out file (TIFF)')
     parser.add_argument('shapex', type=int, help='shapex')
     parser.add_argument('shapey', type=int, help='shapey')
-    parser.add_argument('--has_header', dest='has_header', default=False, help='set True if CSV has header')
-    parser.add_argument('--invert_xy', dest='invert_xy', default=False, help='invert x and y in CSV')
+    parser.add_argument('--has_header', dest='has_header', default=False, help='set True if point file has header')
+    parser.add_argument('--swap_xy', dest='swap_xy', default=False, help='Swap X and Y coordinates')
 
     args = parser.parse_args()
 
     # TOOL
-    points2binaryimage(args.point_file.name, args.out_file, [args.shapey, args.shapex], has_header=args.has_header, invert_xy=args.invert_xy)
+    points2binaryimage(args.point_file.name, args.out_file, (args.shapey, args.shapex), has_header=args.has_header, swap_xy=args.swap_xy)