Skip to content

Commit

Permalink
Separate montage methods from vendor-specific code
Browse files Browse the repository at this point in the history 
The `montage_grid_image_YX` and `montage_stage_pos_image_YX` functions
can be used independently from any input format,
so let's keep them in a separate module.
  • Loading branch information
@imagejan
imagejan committed Oct 13, 2023
1 parent 5bd055f commit d0ec631
Show file tree
Hide file tree
Showing 3 changed files with 244 additions and 230 deletions.
244 changes: 17 additions & 227 deletions src/faim_hcs/MetaSeriesUtils.py
View file
Original file line number Diff line number Diff line change
@@ -1,11 +1,12 @@
from decimal import Decimal
from typing import Any, Callable, Optional, Union
from typing import Callable, Union

import numpy as np
import pandas as pd
from numpy._typing import ArrayLike
from numpy.typing import ArrayLike

from faim_hcs.io.MetaSeriesTiff import load_metaseries_tiff
from faim_hcs.MontageUtils import montage_grid_image_YX
from faim_hcs.UIntHistogram import UIntHistogram
from faim_hcs.utils import rgb_to_hex, wavelength_to_rgb

Expand Down Expand Up @@ -84,226 +85,6 @@ def _z_metadata(metaseries_ch_metadata: dict):
return {"z-position": metaseries_ch_metadata["z-position"]}


def _get_molecular_devices_well_bbox_2D(
data: list[tuple[ArrayLike, dict]]
) -> tuple[Optional[Any], Optional[Any], Optional[Any], Optional[Any]]:
"""Compute well-shape based on stage position metadata."""
assert "stage-position-x" in data[0][1].keys(), "Missing metaseries metadata."
assert "stage-position-y" in data[0][1].keys(), "Missing metaseries metadata."
assert "spatial-calibration-x" in data[0][1].keys(), "Missing metaseries metadata."
assert "spatial-calibration-y" in data[0][1].keys(), "Missing metaseries metadata."

min_x, max_x, min_y, max_y = None, None, None, None
for d in data:
pos_x = d[1]["stage-position-x"]
pos_y = d[1]["stage-position-y"]
res_x = d[1]["spatial-calibration-x"]
res_y = d[1]["spatial-calibration-y"]

if min_x is None:
min_x = pos_x / res_x
max_x = min_x + d[0].shape[1]
elif min_x > (pos_x / res_x):
min_x = pos_x / res_x

if max_x < (pos_x / res_x) + d[0].shape[1]:
max_x = (pos_x / res_x) + d[0].shape[1]

if min_y is None:
min_y = pos_y / res_y
max_y = min_y + d[0].shape[0]
elif min_y > pos_y / res_y:
min_y = pos_y / res_y

if max_y < (pos_y / res_y) + d[0].shape[0]:
max_y = (pos_y / res_y) + d[0].shape[0]

return min_y, min_x, max_y, max_x


def montage_stage_pos_image_YX(data):
"""Montage 2D fields based on stage position metadata.
Montages 2D fields based on stage position metadata. If the stage position
specifies overlapping images, the overlapping part is overwritten
(=> just uses the data of one image). Not well suited for regular grids,
as the stage position can show overlap, but overwriting of data at the
edge is not the intended behavior. In that case, use
`montage_grid_image_YX`.
Also calculates ROI tables for the whole well and the field of views in
the Fractal ROI table format. We only stitch the xy planes here.
Therefore, the z starting position is always 0 and the z extent is set to
1. This is overwritten downsteam if the 2D planes are assembled into a
3D stack.
:param data: list of tuples (image, metadata)
:return: img (stitched 2D np array), fov_df (dataframe with region of
interest information for the fields of view)
"""

def sort_key(d):
label = d[1]["stage-label"]

label = label.split(":")

if len(label) == 1:
return label
else:
return int(label[1].replace("Site", ""))

data.sort(key=sort_key, reverse=True)

min_y, min_x, max_y, max_x = _get_molecular_devices_well_bbox_2D(data)

shape = (int(np.round(max_y - min_y)), int(np.round(max_x - min_x)))

img = np.zeros(shape, dtype=data[0][0].dtype)

fov_rois = []

for d in data:
pos_y = int(
np.round(d[1]["stage-position-y"] / d[1]["spatial-calibration-y"] - min_y)
)
pos_x = int(
np.round(d[1]["stage-position-x"] / d[1]["spatial-calibration-x"] - min_x)
)

img[pos_y : pos_y + d[0].shape[0], pos_x : pos_x + d[0].shape[1]] = d[0]

# Create the FOV ROI table for the site in physical units
fov_rois.append(
(
_stage_label(d[1]),
pos_y * d[1]["spatial-calibration-y"],
pos_x * d[1]["spatial-calibration-x"],
0.0,
d[0].shape[0] * d[1]["spatial-calibration-y"],
d[0].shape[1] * d[1]["spatial-calibration-x"],
1.0,
)
)

roi_tables = create_ROI_tables(fov_rois, shape, calibration_dict=d[1])

return img, roi_tables


def _pixel_pos(dim: str, data: dict):
return np.round(data[f"stage-position-{dim}"] / data[f"spatial-calibration-{dim}"])


def _stage_label(data: dict):
"""Get the field of view (FOV) string for a given FOV dict"""
try:
return data["stage-label"].split(":")[-1][1:]
# Return an empty string if the metadata does not contain stage-label
except KeyError:
return ""


def montage_grid_image_YX(data):
"""Montage 2D fields into fixed grid, based on stage position metadata.
Uses the stage position coordinates to decide which grid cell to put the
image in. Always writes images into a grid, thus avoiding overwriting
partially overwriting parts of images. Not well suited for arbitarily
positioned fields. In that case, use `montage_stage_pos_image_YX`.
Also calculates ROI tables for the whole well and the field of views.
Given that Fractal ROI tables are always 3D, but we only stitch the xy
planes here, the z starting position is always 0 and the
z extent is set to 1. This is overwritten downsteam if the 2D planes are
assembled into a 3D stack.
:param data: list of tuples of (image, metadata)
:return: img (stitched 2D np array), fov_df (dataframe with region of
interest information for the fields of view)
"""
min_y = min(_pixel_pos("y", d[1]) for d in data)
min_x = min(_pixel_pos("x", d[1]) for d in data)
max_y = max(_pixel_pos("y", d[1]) for d in data)
max_x = max(_pixel_pos("x", d[1]) for d in data)

assert all([d[0].shape for d in data])
step_y = data[0][0].shape[0]
step_x = data[0][0].shape[1]

shape = (
int(np.round((max_y - min_y) / step_y + 1) * step_y),
int(np.round((max_x - min_x) / step_x + 1) * step_x),
)
img = np.zeros(shape, dtype=data[0][0].dtype)
fov_rois = []

for d in data:
pos_x = int(np.round((_pixel_pos("x", d[1]) - min_x) / step_x))
pos_y = int(np.round((_pixel_pos("y", d[1]) - min_y) / step_y))
img[
pos_y * step_y : (pos_y + 1) * step_y, pos_x * step_x : (pos_x + 1) * step_x
] = d[0]
# Create the FOV ROI table for the site in physical units
fov_rois.append(
(
_stage_label(d[1]),
pos_y * step_y * d[1]["spatial-calibration-y"],
pos_x * step_x * d[1]["spatial-calibration-x"],
0.0,
step_y * d[1]["spatial-calibration-y"],
step_x * d[1]["spatial-calibration-x"],
1.0,
)
)

roi_tables = create_ROI_tables(fov_rois, shape, calibration_dict=d[1])

return img, roi_tables


def create_ROI_tables(fov_rois, shape, calibration_dict):
columns = [
"FieldIndex",
"x_micrometer",
"y_micrometer",
"z_micrometer",
"len_x_micrometer",
"len_y_micrometer",
"len_z_micrometer",
]
roi_tables = {}
roi_tables["FOV_ROI_table"] = create_fov_ROI_table(fov_rois, columns)
roi_tables["well_ROI_table"] = create_well_ROI_table(
shape[1],
shape[0],
calibration_dict["spatial-calibration-x"],
calibration_dict["spatial-calibration-y"],
columns,
)
return roi_tables


def create_well_ROI_table(shape_x, shape_y, pixel_size_x, pixel_size_y, columns):
well_roi = [
"well_1",
0.0,
0.0,
0.0,
shape_x * pixel_size_x,
shape_y * pixel_size_y,
1.0,
]
well_roi_table = pd.DataFrame(well_roi).T
well_roi_table.columns = columns
well_roi_table.set_index("FieldIndex", inplace=True)
return well_roi_table


def create_fov_ROI_table(fov_rois, columns):
roi_table = pd.DataFrame(fov_rois, columns=columns).set_index("FieldIndex")
return roi_table


def verify_integrity(field_metadata: list[dict]):
metadata = field_metadata[0]
for fm in field_metadata:
Expand Down Expand Up @@ -361,7 +142,7 @@ def get_well_image_CZYX(
well_files: pd.DataFrame,
channels: list[str],
assemble_fn: Callable = montage_grid_image_YX,
) -> tuple[ArrayLike, list[UIntHistogram], list[dict], dict]:
) -> tuple[ArrayLike, list[UIntHistogram], list[dict], dict, dict]:
"""Assemble image data for the given well-files."""
planes = well_files["z"].unique()

Expand Down Expand Up @@ -437,11 +218,11 @@ def get_well_image_CYX(
well_files: pd.DataFrame,
channels: list[str],
assemble_fn: Callable = montage_grid_image_YX,
) -> tuple[ArrayLike, list[UIntHistogram], list[dict], dict]:
"""Assemble image data for the given well-files.
) -> tuple[ArrayLike, list[UIntHistogram], list[dict], dict, dict]:
"""Assemble image data for the given well files.
For each channel a single 2D image is computed. If the well has multiple
fields per channel the `assemble_fn` has to montage or stitch the fields
For each channel, a single 2D image is computed. If the well has multiple
fields per channel, the `assemble_fn` has to montage or stitch the fields
accordingly.
:param well_files: all files corresponding to the well
Expand Down Expand Up @@ -490,6 +271,15 @@ def get_well_image_CYX(


def get_img_YX(assemble_fn, files):
"""Assemble single 2D image for all fields.
general_metadata: spatial-calibration-x, spatial-calibration-y, spatial-calibration-units, pixel-type
img: 2D pixel array (yx)
metadata: list[dict] (wavelength, power, exposure-time, exposure-time-unit, shading-correction,
channel-name, objective-numerical-aperture, objective, display-color)
z_position
roi_tables: dict[DataFrame]
"""
imgs = []
field_metadata = []
z_positions = []
Expand Down
Loading

0 comments on commit d0ec631

Please sign in to comment.