Merge pull request #106 from fractal-analytics-platform/91_off_by_one

Update ROIs

Author: tcompa

fractal_tasks_core/create_zarr_structure.py

@@ -34,12 +34,10 @@
 
 
 def define_omero_channels(actual_channels, channel_parameters, bit_depth):
-    from devtools import debug
 
     omero_channels = []
     default_colormaps = ["00FFFF", "FF00FF", "FFFF00"]
     for channel in actual_channels:
-        debug(channel_parameters[channel])
 
         # Set colormap. If missing, use the default ones (for the first three
         # channels) or gray
@@ -64,7 +62,6 @@ def define_omero_channels(actual_channels, channel_parameters, bit_depth):
                 },
             }
         )
-        debug(omero_channels[-1])
 
         try:
             omero_channels[-1]["window"]["start"] = channel_parameters[

fractal_tasks_core/lib_regions_of_interest.py

@@ -28,18 +28,11 @@ def prepare_FOV_ROI_table(
     # when creating AnnData object.
     # Do this in the beginning to allow concatenation with e.g. time
     df.index = df.index.astype(str)
-    from devtools import debug
 
-    debug(df)
-    # Calculate bounding box extents in physical units
-
-    for mu in ["x", "y", "z"]:
-
-        # Reset reference values for coordinates
-        df[f"{mu}_micrometer"] -= df[f"{mu}_micrometer"].min()
-
-        # Obtain box size in physical units
-        df[f"len_{mu}_micrometer"] = df[f"{mu}_pixel"] * df[f"pixel_size_{mu}"]
+    # Obtain box size in physical units
+    df = df.assign(len_x_micrometer=df.x_pixel * df.pixel_size_x)
+    df = df.assign(len_y_micrometer=df.y_pixel * df.pixel_size_y)
+    df = df.assign(len_z_micrometer=df.z_pixel * df.pixel_size_z)
 
     # Select only the numeric positional columns needed to define ROIs
     # (to avoid) casting things like the data column to float32
@@ -51,6 +44,8 @@ def prepare_FOV_ROI_table(
         "len_x_micrometer",
         "len_y_micrometer",
         "len_z_micrometer",
+        "x_micrometer_original",
+        "y_micrometer_original",
     ]
 
     # Assign dtype explicitly, to avoid
@@ -172,6 +167,16 @@ def convert_ROI_table_to_indices(
     level: int = 0,
     coarsening_xy: int = 2,
     full_res_pxl_sizes_zyx: Iterable[float] = None,
+    cols_xyz_pos: Iterable[str] = [
+        "x_micrometer",
+        "y_micrometer",
+        "z_micrometer",
+    ],
+    cols_xyz_len: Iterable[str] = [
+        "len_x_micrometer",
+        "len_y_micrometer",
+        "len_z_micrometer",
+    ],
 ) -> List[List[int]]:
 
     # Set pyramid-level pixel sizes
@@ -180,16 +185,23 @@ def convert_ROI_table_to_indices(
     pxl_size_x *= prefactor
     pxl_size_y *= prefactor
 
+    x_pos, y_pos, z_pos = cols_xyz_pos[:]
+    x_len, y_len, z_len = cols_xyz_len[:]
+
+    origin_x = min(ROI[:, x_pos].X[:, 0])
+    origin_y = min(ROI[:, y_pos].X[:, 0])
+    origin_z = min(ROI[:, z_pos].X[:, 0])
+
     list_indices = []
     for FOV in ROI.obs_names:
 
         # Extract data from anndata table
-        x_micrometer = ROI[FOV, "x_micrometer"].X[0, 0]
-        y_micrometer = ROI[FOV, "y_micrometer"].X[0, 0]
-        z_micrometer = ROI[FOV, "z_micrometer"].X[0, 0]
-        len_x_micrometer = ROI[FOV, "len_x_micrometer"].X[0, 0]
-        len_y_micrometer = ROI[FOV, "len_y_micrometer"].X[0, 0]
-        len_z_micrometer = ROI[FOV, "len_z_micrometer"].X[0, 0]
+        x_micrometer = ROI[FOV, x_pos].X[0, 0] - origin_x
+        y_micrometer = ROI[FOV, y_pos].X[0, 0] - origin_y
+        z_micrometer = ROI[FOV, z_pos].X[0, 0] - origin_z
+        len_x_micrometer = ROI[FOV, x_len].X[0, 0]
+        len_y_micrometer = ROI[FOV, y_len].X[0, 0]
+        len_z_micrometer = ROI[FOV, z_len].X[0, 0]
 
         # Identify indices along the three dimensions
         start_x = x_micrometer / pxl_size_x

fractal_tasks_core/lib_remove_FOV_overlaps.py

@@ -0,0 +1,145 @@
+def is_overlapping_1D(line1, line2, tol=0):
+    """
+    Based on https://stackoverflow.com/a/70023212/19085332
+
+    line:
+        (xmin, xmax)
+    """
+    return line1[0] <= line2[1] - tol and line2[0] <= line1[1] - tol
+
+
+def is_overlapping_2D(box1, box2, tol=0):
+    """
+    Based on https://stackoverflow.com/a/70023212/19085332
+
+    box:
+        (xmin, ymin, xmax, ymax)
+    """
+    overlap_x = is_overlapping_1D(
+        [box1[0], box1[2]], [box2[0], box2[2]], tol=tol
+    )
+    overlap_y = is_overlapping_1D(
+        [box1[1], box1[3]], [box2[1], box2[3]], tol=tol
+    )
+    return overlap_x and overlap_y
+
+
+def get_overlapping_pair(tmp_df, tol=0):
+    # NOTE: here we use positional indices (i.e. starting from 0)
+    num_lines = len(tmp_df.index)
+    for pos_ind_1 in range(num_lines):
+        for pos_ind_2 in range(pos_ind_1):
+            if is_overlapping_2D(
+                tmp_df.iloc[pos_ind_1], tmp_df.iloc[pos_ind_2], tol=tol
+            ):
+                return (pos_ind_1, pos_ind_2)
+    return False
+
+
+def remove_FOV_overlaps(df):
+
+    # Set tolerance (this should be much smaller than pixel size or expected
+    # round-offs), and maximum number of iterations in constraint solver
+    tol = 1e-10
+    max_iterations = 200
+
+    # Create temporary columns (to streamline overlap removals), which are
+    # then removed at the end of the remove_FOV_overlaps function
+    df["xmin"] = df["x_micrometer"]
+    df["ymin"] = df["y_micrometer"]
+    df["xmax"] = df["x_micrometer"] + df["pixel_size_x"] * df["x_pixel"]
+    df["ymax"] = df["y_micrometer"] + df["pixel_size_y"] * df["y_pixel"]
+    list_columns = ["xmin", "ymin", "xmax", "ymax"]
+
+    # Create columns with the original positions (not to be removed)
+    df["x_micrometer_original"] = df["x_micrometer"]
+    df["y_micrometer_original"] = df["y_micrometer"]
+
+    # Check that tolerance is much smaller than pixel sizes
+    min_pixel_size = df[["pixel_size_x", "pixel_size_y"]].min().min()
+    if tol > min_pixel_size / 1e3:
+        raise Exception(
+            f"In remove_FOV_overlaps, {tol=} but {min_pixel_size=}"
+        )
+
+    # Loop over wells
+    wells = sorted(list(set([ind[0] for ind in df.index])))
+    for well in wells:
+        print(f"[remove_FOV_overlaps] removing FOV overlaps for {well=}")
+
+        # NOTE: these are positional indices (i.e. starting from 0)
+        pair_pos_indices = get_overlapping_pair(
+            df.loc[well][list_columns], tol=tol
+        )
+
+        # Keep going until there are no overlaps, or until iteration reaches
+        # max_iterations
+        iteration = 0
+        while pair_pos_indices:
+            iteration += 1
+
+            # Identify overlapping FOVs
+            pos_ind_1, pos_ind_2 = pair_pos_indices
+            fov_id_1 = df.loc[well].index[pos_ind_1]
+            fov_id_2 = df.loc[well].index[pos_ind_2]
+            xmin_1, ymin_1, xmax_1, ymax_1 = df.loc[well][list_columns].iloc[
+                pos_ind_1
+            ]
+            xmin_2, ymin_2, xmax_2, ymax_2 = df.loc[well][list_columns].iloc[
+                pos_ind_2
+            ]
+
+            print(
+                f"[remove_FOV_overlaps] {iteration=}, removing overlap between"
+                f" {fov_id_1=} and {fov_id_2=}"
+            )
+
+            # Two overlapping FOVs MUST share either a vertical boundary or an
+            # horizontal one
+            is_x_equal = abs(xmin_1 - xmin_2) < tol and (xmax_1 - xmax_2) < tol
+            is_y_equal = abs(ymin_1 - ymin_2) < tol and (ymax_1 - ymax_2) < tol
+            if is_x_equal + is_y_equal != 1:
+                raise Exception(
+                    "Two overlapping FOVs MUST share either a "
+                    "vertical boundary or an horizontal one, but "
+                    f"{is_x_equal=} and {is_y_equal=}"
+                )
+
+            # Compute and apply shift (either along x or y)
+            if is_y_equal:
+                min_max_x = min(xmax_1, xmax_2)
+                max_min_x = max(xmin_1, xmin_2)
+                shift_x = min_max_x - max_min_x
+                ind = df.loc[well].loc[:, "xmin"] >= max_min_x - tol
+                if not (shift_x > 0.0 and ind.to_numpy().max() > 0):
+                    raise Exception(
+                        "Something wrong in remove_FOV_overlaps\n"
+                        f"{shift_x=}\n"
+                        f"{ind.to_numpy()=}"
+                    )
+                df.loc[well].loc[ind, "xmin"] += shift_x
+                df.loc[well].loc[ind, "xmax"] += shift_x
+                df.loc[well].loc[ind, "x_micrometer"] += shift_x
+            if is_x_equal:
+                min_max_y = min(ymax_1, ymax_2)
+                max_min_y = max(ymin_1, ymin_2)
+                shift_y = min_max_y - max_min_y
+                ind = df.loc[well].loc[:, "ymin"] >= max_min_y - tol
+                if not (shift_y > 0.0 and ind.to_numpy().max() > 0):
+                    raise Exception(
+                        "Something wrong in remove_FOV_overlaps\n"
+                        f"{shift_y=}\n"
+                        f"{ind.to_numpy()=}"
+                    )
+                df.loc[well].loc[ind, "ymin"] += shift_y
+                df.loc[well].loc[ind, "ymax"] += shift_y
+                df.loc[well].loc[ind, "y_micrometer"] += shift_y
+
+            pair_pos_indices = get_overlapping_pair(
+                df.loc[well][list_columns], tol=tol
+            )
+            if iteration > max_iterations:
+                raise Exception("Something went wrong, {num_iteration=}")
+    df.drop(list_columns, axis=1, inplace=True)
+
+    return df

fractal_tasks_core/metadata_parsing.py

@@ -16,6 +16,8 @@
 import pandas as pd
 from defusedxml import ElementTree
 
+from .lib_remove_FOV_overlaps import remove_FOV_overlaps
+
 
 def parse_yokogawa_metadata(mrf_path, mlf_path):
     """
@@ -77,6 +79,8 @@ def parse_yokogawa_metadata(mrf_path, mlf_path):
     # Maybe return it here for further checks and produce a warning if it does
     # not match
 
+    site_metadata = remove_FOV_overlaps(site_metadata)
+
     return site_metadata, total_files
 
 

tests/test_unit_ROI_indices.py

@@ -27,12 +27,14 @@ def get_metadata_dataframe():
     """
     Create artificial metadata dataframe
     """
-    df = pd.DataFrame(np.zeros((4, 11)), dtype=int)
+    df = pd.DataFrame(np.zeros((4, 13)), dtype=float)
     df.index = FOV_IDS
     df.columns = [
         "x_micrometer",
         "y_micrometer",
         "z_micrometer",
+        "x_micrometer_original",
+        "y_micrometer_original",
         "x_pixel",
         "y_pixel",
         "z_pixel",
@@ -56,6 +58,8 @@ def get_metadata_dataframe():
         img_size_y_micrometer,
         img_size_y_micrometer,
     ]
+    df["x_micrometer_original"] = df["x_micrometer"]
+    df["y_micrometer_original"] = df["y_micrometer"]
     df["z_micrometer"] = [0.0, 0.0, 0.0, 0.0]
     df["x_pixel"] = [IMG_SIZE_X] * 4
     df["y_pixel"] = [IMG_SIZE_Y] * 4

tests/test_unit_parse_yokogawa_metadata.py

@@ -27,7 +27,7 @@
 mlf_path_1 = f"{path}MeasurementData_SingleWell2Sites_MultiZ.mlf"
 mrf_path_1 = f"{path}MeasurementDetail.mrf"
 expected_files_1 = 4
-expected_shape_1 = (2, 11)
+expected_shape_1 = (2, 13)
 x_mic_pos_1 = (-1448.3, -1032.3)
 y_mic_pos_1 = -1517.7
 pixel_size_z_1 = 1.0
@@ -46,7 +46,7 @@
 mlf_path_2 = f"{path}MeasurementData_2x2_well.mlf"
 mrf_path_2 = f"{path}MeasurementDetail_2x2_well.mrf"
 expected_files_2 = 120
-expected_shape_2 = (4, 11)
+expected_shape_2 = (4, 13)
 x_mic_pos_2 = (-1448.3, -1032.3)
 y_mic_pos_2 = (-1517.7, -1166.7)
 pixel_size_z_2 = 1.0
@@ -62,7 +62,7 @@
 mlf_path_3 = f"{path}MeasurementData_SingleWell2Sites_SingleZ.mlf"
 mrf_path_3 = f"{path}MeasurementDetail.mrf"
 expected_files_3 = 2
-expected_shape_3 = (2, 11)
+expected_shape_3 = (2, 13)
 x_mic_pos_3 = (-1448.3, -1032.3)
 y_mic_pos_3 = -1517.7
 pixel_size_z_3 = 1.0