7105148444

Committed 05 Dec 2023 06:42PM UTC coverage: 93.707% (+0.06%) from 93.651%

Build # 7105148444

Build Type

Pull #1092

github

Committed by

web-flow

Commit Message

Merge 2575a214a into 99a9abbda

Pull Request Pull Request #1092: Generalized masking

Run Details

2052 of 2287 branches covered (0.0%)

Branch coverage included in aggregate %.

47 of 47 new or added lines in 3 files covered. (100.0%)

7 existing lines in 1 file now uncovered.

4008 of 4180 relevant lines covered (95.89%)

5.72 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

80.77

/src/ark/segmentation/ez_seg/ez_object_segmentation.py

import pathlib
from typing import Optional, Tuple, Union, Literal
import numpy as np
from skimage import measure, filters, morphology
from skimage.util import map_array
import pandas as pd
from alpineer import misc_utils, load_utils, image_utils, io_utils
from ark.segmentation.ez_seg.ez_seg_utils import log_creator
import xarray as xr
import warnings


def create_object_masks(
        image_data_dir: Union[str, pathlib.Path],
        img_sub_folder: Optional[str],
        fov_list: list[str],
        mask_name: str,
        channel_to_segment: str,
        masks_dir: Union[str, pathlib.Path],
        log_dir: Union[str, pathlib.Path],
        object_shape_type: str = "blob",
        sigma: int = 1,
        thresh: Optional[np.float32] = None,
        hole_size: Optional[int] = None,
        fov_dim: int = 400,
        min_object_area: int = 100,
        max_object_area: int = 100000,
) -> None:
    """
    Calculates a mask for each channel in the FOV for circular or 'blob'-like objects such as: single large cells or amyloid
    plaques. It will blur the input image, then threshold the blurred image on either a given
    fixed value, or an adaptive thresholding method. In addition, it removes small holes using
    that same thresholding input and filters out objects which are either too small or too large.

    Args:
        image_data_dir (Union[str, pathlib.Path]): The directory to pull images from to perform segmentation on.
        img_sub_folder (str): A name for sub-folders within each fov in the image_data location.
        fov_list: A list of fov names to segment on.
        mask_name (str): The name of the masks you are creating.
        channel_to_segment: The channel on which to perform segmentation.
        masks_dir (Union[str, pathlib.Path]): The directory to save segmented images to.
        object_shape_type (str, optional): Specify whether the object is either "blob" or
        "projection" shaped. Defaults to "blob".
        sigma (int): The standard deviation for Gaussian kernel, used for blurring the
        image. Defaults to 1.
        thresh (np.float32, optional): The global threshold value for image thresholding if
        desired. Defaults to None.
        hole_size (int, optional): A specific area to close small holes over in object masks.
        Defaults to None.
        fov_dim (int): The dimension in μm of the FOV.
        min_object_area (int): The minimum size (area) of an object to capture in
        pixels. Defaults to 100.
        max_object_area (int): The maximum size (area) of an object to capture in
        pixels. Defaults to 100000.
        log_dir (Union[str, pathlib.Path]): The directory to save log information to.
    """

    # Input validation
    io_utils.validate_paths([image_data_dir, masks_dir, log_dir])

    misc_utils.verify_in_list(
        object_shape=[object_shape_type], object_shape_options=["blob", "projection"]
    )

    for fov in fov_list:
        fov_xr: xr.DataArray = load_utils.load_imgs_from_tree(
            data_dir=image_data_dir, img_sub_folder=img_sub_folder, fovs=fov
        ).squeeze()

        # handles folders where only 1 channel is loaded, often for composites
        try:
            len(fov_xr.channels)
            channel: xr.DataArray = fov_xr.sel({"channels": channel_to_segment}).astype(
                int
            )
        except TypeError:
            channel: xr.DataArray = fov_xr.astype(int)

        object_masks: np.ndarray = _create_object_mask(
            input_image=channel,
            object_shape_type=object_shape_type,
            sigma=sigma,
            thresh=thresh,
            hole_size=hole_size,
            fov_dim=fov_dim,
            min_object_area=min_object_area,
            max_object_area=max_object_area,
        )

        # save the channel overlay
        save_path = pathlib.Path(masks_dir) / f"{fov}_{mask_name}.tiff"
        image_utils.save_image(fname=save_path, data=object_masks)

    # Write a log saving ez segment info
    variables_to_log = {
        "image_data_dir": image_data_dir,
        "fov_list": fov_list,
        "mask_name": mask_name,
        "channel_to_segment": channel_to_segment,
        "masks_dir": masks_dir,
        "object_shape_type": object_shape_type,
        "sigma": sigma,
        "thresh": thresh,
        "hole_size": hole_size,
        "fov_dim": fov_dim,
        "min_object_area": min_object_area,
        "max_object_area": max_object_area,
    }
    log_creator(variables_to_log, log_dir, f"{mask_name}_segmentation_log.txt")
    print("ez masks built and saved")


def _create_object_mask(
        input_image: xr.DataArray,
        object_shape_type: Union[Literal["blob"], Literal["projection"]] = "blob",
        sigma: int = 1,
        thresh: Union[int, Literal["auto"]] = None,
        hole_size: Union[int, Literal["auto"]] = "auto",
        fov_dim: int = 400,
        min_object_area: int = 10,
        max_object_area: int = 100000,
) -> np.ndarray:
    """
    Calculates a mask for circular or 'blob'-like objects such as: single large cells or amyloid
    plaques. It will blur the input image, then threshold the blurred image on either a given
    fixed value, or an adaptive thresholding method. In addition, it removes small holes using
    that same thresholding input and filters out objects which are either too small or too large.

    Args:
        input_image (xr.DataArray): The numpy array (image) to perform segmentation on.
        object_shape_type (str, optional): Specify whether the object is either "blob" or
        "projection" shaped. Defaults to "blob".
        sigma (int): The standard deviation for Gaussian kernel, used for blurring the
        image. Defaults to 1.
        thresh (int, str,  optional): The global threshold value for image thresholding if
        desired. Defaults to "auto".
        hole_size (int, str, optional): A specific area to close small holes over in object masks.
        Defaults to None.
        fov_dim (int): The dimension in μm of the FOV.
        min_object_area (int): The minimum size (area) of an object to capture in
        pixels. Defaults to 100.
        max_object_area (int): The maximum size (area) of an object to capture in
        pixels. Defaults to 100000.

    Returns:
        np.ndarray: The object mask.
    """

    # Do not display any UserWarning msg's about boolean arrays here.
    warnings.filterwarnings(
        "ignore", message="Any labeled images will be returned as a boolean array. Did you mean to use a boolean array?")

    # Input validation
    misc_utils.verify_in_list(object_shape_type=[object_shape_type], object_shape_options=["blob", "projection"])

    # Copy the input image, and get its shape
    img2mask: np.ndarray = input_image.copy().to_numpy()
    img2mask_shape: Tuple[int, int] = img2mask.shape

    # Blur the input mask using given sigma value
    if sigma is None:
        img2mask_blur = img2mask
    else:
        img2mask_blur: np.ndarray = filters.gaussian(
            img2mask, sigma=sigma, preserve_range=True
        )

    # Apply binary thresholding to the blurred image
    if isinstance(thresh, int):
        # Find the threshold value based on the given percentile number
        img_nonzero = img2mask_blur[img2mask_blur != 0]
        thresh_percentile = np.percentile(img_nonzero, thresh)

        # Threshold the array where values below the threshold are set to 0
        img2mask_thresh = np.where(img2mask_blur < thresh_percentile, 0, img2mask_blur)

    elif thresh == "auto":
        local_thresh_block_size: int = get_block_size(
            block_type="local_thresh", fov_dim=fov_dim, img_shape=img2mask_shape[0]
        )
        img2mask_thresh: np.ndarray = img2mask_blur > filters.threshold_local(
            img2mask_blur, block_size=local_thresh_block_size
        )
    elif thresh is None:
        img2mask_thresh = img2mask_blur
    else:
        raise ValueError(
            f"Invalid `threshold` value: {thresh}. Must be either `auto`, `None` or an integer."
        )

    # Binarize the image in prep for removing holes.
    img2mask_thresh_binary = img2mask_thresh > 0
    img2mask_thresh[img2mask_thresh_binary] = 1
    img2mask_thresh = img2mask_thresh.astype(int)

    # Remove small holes within the objects
    if isinstance(hole_size, int):
        img2mask_rm_holes: np.ndarray = morphology.remove_small_holes(
            img2mask_thresh, area_threshold=hole_size
        )
    elif hole_size == "auto":
        small_holes_block_size: int = get_block_size(
            block_type="small_holes", fov_dim=fov_dim, img_shape=img2mask_shape[0]
        )
        img2mask_rm_holes: np.ndarray = morphology.remove_small_holes(
            img2mask_thresh, area_threshold=small_holes_block_size
        )
    elif hole_size is None:
        img2mask_rm_holes = img2mask_thresh
    else:
        raise ValueError(
            f"Invalid `hole_size` value: {hole_size}. Must be either `auto`, `None` or an integer."
        )

    # Filter projections
    if object_shape_type == "projection":
        img2mask_filtered: np.ndarray = filters.meijering(
            img2mask_rm_holes, sigmas=range(1, 5, 1), black_ridges=False
        )
    else:
        img2mask_filtered: np.ndarray = img2mask_rm_holes

    # Binarize the image in prep for labeling.
    img2mask_filtered_binary = img2mask_filtered > 0
    img2mask_filtered[img2mask_filtered_binary] = 1

    # Extract `label` and `area` from regionprops
    labeled_object_masks = measure.label(img2mask_filtered, connectivity=2)

    # Convert dictionary of region properties to DataFrame
    object_masks_df: pd.DataFrame = pd.DataFrame(
        measure.regionprops_table(
            label_image=labeled_object_masks,
            cache=True,
            properties=[
                "label",
                "area",
            ],
        )
    )

    # zero out objects not meeting size requirements
    keep_labels_bool = (object_masks_df["area"] >= min_object_area) & (
            object_masks_df["area"] <= max_object_area
    )
    all_labels = object_masks_df["label"]
    labels_to_keep = all_labels * keep_labels_bool

    # Map filtered objects to the object mask
    objects_filtered_by_size = map_array(
        labeled_object_masks, all_labels.to_numpy(), labels_to_keep.to_numpy()
    )

    return objects_filtered_by_size.astype(np.int32)


def get_block_size(block_type: str, fov_dim: int, img_shape: int) -> int:
    """
    Computes the approximate local otsu threshold based on fov size (in μm) and pixel resolution.

    Args:
        block_type (str): Either "small_holes" or "local_thresh"
        fov_dim (int, optional): The size in μm for the FOV.
        img_shape (int, optional): The shape of the image.

    Returns:
        int: Returns the approximate block area
    """

    # Input validation
    misc_utils.verify_in_list(
        block_type=[block_type], block_types=["small_holes", "local_thresh"]
    )
    # Get the size of the pixel

    pixel_size = fov_dim / img_shape
    # grab block size for removing small holes
    if block_type == "small_holes":
        size = (np.pi * 5) ** 2 / pixel_size
        # round up above value to make it into an integer
        area: int = round(size)
    # grab local threshold block size
    else:
        # use this to calculate out how many pixels it takes to get to roughly 10 μm
        # (roughly a cell soma diameter)
        size: float = 10 / pixel_size

        # round the area up to the nearest odd number
        area: int = round(size)
        if area % 2 == 0:
            area += 1
    return area

1	import pathlib	6✔
2	from typing import Optional, Tuple, Union, Literal	6✔
3	import numpy as np	6✔
4	from skimage import measure, filters, morphology	6✔
5	from skimage.util import map_array	6✔
6	import pandas as pd	6✔
7	from alpineer import misc_utils, load_utils, image_utils, io_utils	6✔
8	from ark.segmentation.ez_seg.ez_seg_utils import log_creator	6✔
9	import xarray as xr	6✔
10	import warnings	6✔
11
12
13	def create_object_masks(	6✔
14	image_data_dir: Union[str, pathlib.Path],
15	img_sub_folder: Optional[str],
16	fov_list: list[str],
17	mask_name: str,
18	channel_to_segment: str,
19	masks_dir: Union[str, pathlib.Path],
20	log_dir: Union[str, pathlib.Path],
21	object_shape_type: str = "blob",
22	sigma: int = 1,
23	thresh: Optional[np.float32] = None,
24	hole_size: Optional[int] = None,
25	fov_dim: int = 400,
26	min_object_area: int = 100,
27	max_object_area: int = 100000,
28	) -> None:
29	"""
30	Calculates a mask for each channel in the FOV for circular or 'blob'-like objects such as: single large cells or amyloid
31	plaques. It will blur the input image, then threshold the blurred image on either a given
32	fixed value, or an adaptive thresholding method. In addition, it removes small holes using
33	that same thresholding input and filters out objects which are either too small or too large.
34
35	Args:
36	image_data_dir (Union[str, pathlib.Path]): The directory to pull images from to perform segmentation on.
37	img_sub_folder (str): A name for sub-folders within each fov in the image_data location.
38	fov_list: A list of fov names to segment on.
39	mask_name (str): The name of the masks you are creating.
40	channel_to_segment: The channel on which to perform segmentation.
41	masks_dir (Union[str, pathlib.Path]): The directory to save segmented images to.
42	object_shape_type (str, optional): Specify whether the object is either "blob" or
43	"projection" shaped. Defaults to "blob".
44	sigma (int): The standard deviation for Gaussian kernel, used for blurring the
45	image. Defaults to 1.
46	thresh (np.float32, optional): The global threshold value for image thresholding if
47	desired. Defaults to None.
48	hole_size (int, optional): A specific area to close small holes over in object masks.
49	Defaults to None.
50	fov_dim (int): The dimension in μm of the FOV.
51	min_object_area (int): The minimum size (area) of an object to capture in
52	pixels. Defaults to 100.
53	max_object_area (int): The maximum size (area) of an object to capture in
54	pixels. Defaults to 100000.
55	log_dir (Union[str, pathlib.Path]): The directory to save log information to.
56	"""
57
58	# Input validation
59	io_utils.validate_paths([image_data_dir, masks_dir, log_dir])	6✔
60
61	misc_utils.verify_in_list(	6✔
62	object_shape=[object_shape_type], object_shape_options=["blob", "projection"]
63	)
64
65	for fov in fov_list:	6✔
66	fov_xr: xr.DataArray = load_utils.load_imgs_from_tree(	6✔
67	data_dir=image_data_dir, img_sub_folder=img_sub_folder, fovs=fov
68	).squeeze()
69
70	# handles folders where only 1 channel is loaded, often for composites
71	try:	6✔
72	len(fov_xr.channels)	6✔
73	channel: xr.DataArray = fov_xr.sel({"channels": channel_to_segment}).astype(	6✔
74	int
75	)
76	except TypeError:	×
77	channel: xr.DataArray = fov_xr.astype(int)	×
78
79	object_masks: np.ndarray = _create_object_mask(	6✔
80	input_image=channel,
81	object_shape_type=object_shape_type,
82	sigma=sigma,
83	thresh=thresh,
84	hole_size=hole_size,
85	fov_dim=fov_dim,
86	min_object_area=min_object_area,
87	max_object_area=max_object_area,
88	)
89
90	# save the channel overlay
91	save_path = pathlib.Path(masks_dir) / f"{fov}_{mask_name}.tiff"	6✔
92	image_utils.save_image(fname=save_path, data=object_masks)	6✔
93
94	# Write a log saving ez segment info
95	variables_to_log = {	6✔
96	"image_data_dir": image_data_dir,
97	"fov_list": fov_list,
98	"mask_name": mask_name,
99	"channel_to_segment": channel_to_segment,
100	"masks_dir": masks_dir,
101	"object_shape_type": object_shape_type,
102	"sigma": sigma,
103	"thresh": thresh,
104	"hole_size": hole_size,
105	"fov_dim": fov_dim,
106	"min_object_area": min_object_area,
107	"max_object_area": max_object_area,
108	}
109	log_creator(variables_to_log, log_dir, f"{mask_name}_segmentation_log.txt")	6✔
110	print("ez masks built and saved")	6✔
111
112
113	def _create_object_mask(	6✔
114	input_image: xr.DataArray,
115	object_shape_type: Union[Literal["blob"], Literal["projection"]] = "blob",
116	sigma: int = 1,
117	thresh: Union[int, Literal["auto"]] = None,
118	hole_size: Union[int, Literal["auto"]] = "auto",
119	fov_dim: int = 400,
120	min_object_area: int = 10,
121	max_object_area: int = 100000,
122	) -> np.ndarray:
123	"""
124	Calculates a mask for circular or 'blob'-like objects such as: single large cells or amyloid
125	plaques. It will blur the input image, then threshold the blurred image on either a given
126	fixed value, or an adaptive thresholding method. In addition, it removes small holes using
127	that same thresholding input and filters out objects which are either too small or too large.
128
129	Args:
130	input_image (xr.DataArray): The numpy array (image) to perform segmentation on.
131	object_shape_type (str, optional): Specify whether the object is either "blob" or
132	"projection" shaped. Defaults to "blob".
133	sigma (int): The standard deviation for Gaussian kernel, used for blurring the
134	image. Defaults to 1.
135	thresh (int, str, optional): The global threshold value for image thresholding if
136	desired. Defaults to "auto".
137	hole_size (int, str, optional): A specific area to close small holes over in object masks.
138	Defaults to None.
139	fov_dim (int): The dimension in μm of the FOV.
140	min_object_area (int): The minimum size (area) of an object to capture in
141	pixels. Defaults to 100.
142	max_object_area (int): The maximum size (area) of an object to capture in
143	pixels. Defaults to 100000.
144
145	Returns:
146	np.ndarray: The object mask.
147	"""
148
149	# Do not display any UserWarning msg's about boolean arrays here.
150	warnings.filterwarnings(	6✔
151	"ignore", message="Any labeled images will be returned as a boolean array. Did you mean to use a boolean array?")
152
153	# Input validation
154	misc_utils.verify_in_list(object_shape_type=[object_shape_type], object_shape_options=["blob", "projection"])	6✔
155
156	# Copy the input image, and get its shape
157	img2mask: np.ndarray = input_image.copy().to_numpy()	6✔
158	img2mask_shape: Tuple[int, int] = img2mask.shape	6✔
159
160	# Blur the input mask using given sigma value
161	if sigma is None:	6✔
UNCOV 162	img2mask_blur = img2mask	×
163	else:
164	img2mask_blur: np.ndarray = filters.gaussian(	6✔
165	img2mask, sigma=sigma, preserve_range=True
166	)
167
168	# Apply binary thresholding to the blurred image
169	if isinstance(thresh, int):	6✔
170	# Find the threshold value based on the given percentile number
171	img_nonzero = img2mask_blur[img2mask_blur != 0]	6✔
172	thresh_percentile = np.percentile(img_nonzero, thresh)	6✔
173
174	# Threshold the array where values below the threshold are set to 0
175	img2mask_thresh = np.where(img2mask_blur < thresh_percentile, 0, img2mask_blur)	6✔
176
177	elif thresh == "auto":	6✔
UNCOV 178	local_thresh_block_size: int = get_block_size(	×
179	block_type="local_thresh", fov_dim=fov_dim, img_shape=img2mask_shape[0]
180	)
UNCOV 181	img2mask_thresh: np.ndarray = img2mask_blur > filters.threshold_local(	×
182	img2mask_blur, block_size=local_thresh_block_size
183	)
184	elif thresh is None:	6✔
185	img2mask_thresh = img2mask_blur	6✔
186	else:
UNCOV 187	raise ValueError(	×
188	f"Invalid `threshold` value: {thresh}. Must be either `auto`, `None` or an integer."
189	)
190
191	# Binarize the image in prep for removing holes.
192	img2mask_thresh_binary = img2mask_thresh > 0	6✔
193	img2mask_thresh[img2mask_thresh_binary] = 1	6✔
194	img2mask_thresh = img2mask_thresh.astype(int)	6✔
195
196	# Remove small holes within the objects
197	if isinstance(hole_size, int):	6✔
198	img2mask_rm_holes: np.ndarray = morphology.remove_small_holes(	6✔
199	img2mask_thresh, area_threshold=hole_size
200	)
201	elif hole_size == "auto":	3!
202	small_holes_block_size: int = get_block_size(	3!
203	block_type="small_holes", fov_dim=fov_dim, img_shape=img2mask_shape[0]
204	)
205	img2mask_rm_holes: np.ndarray = morphology.remove_small_holes(	3✔
206	img2mask_thresh, area_threshold=small_holes_block_size
207	)
UNCOV 208	elif hole_size is None:	×
UNCOV 209	img2mask_rm_holes = img2mask_thresh	×
210	else:
211	raise ValueError(	×
212	f"Invalid `hole_size` value: {hole_size}. Must be either `auto`, `None` or an integer."
213	)
214
215	# Filter projections
216	if object_shape_type == "projection":	6✔
217	img2mask_filtered: np.ndarray = filters.meijering(	6✔
218	img2mask_rm_holes, sigmas=range(1, 5, 1), black_ridges=False
219	)
220	else:
221	img2mask_filtered: np.ndarray = img2mask_rm_holes	6✔
222
223	# Binarize the image in prep for labeling.
224	img2mask_filtered_binary = img2mask_filtered > 0	6✔
225	img2mask_filtered[img2mask_filtered_binary] = 1	6✔
226
227	# Extract `label` and `area` from regionprops
228	labeled_object_masks = measure.label(img2mask_filtered, connectivity=2)	6✔
229
230	# Convert dictionary of region properties to DataFrame
231	object_masks_df: pd.DataFrame = pd.DataFrame(	6✔
232	measure.regionprops_table(
233	label_image=labeled_object_masks,
234	cache=True,
235	properties=[
236	"label",
237	"area",
238	],
239	)
240	)
241
242	# zero out objects not meeting size requirements
243	keep_labels_bool = (object_masks_df["area"] >= min_object_area) & (	6✔
244	object_masks_df["area"] <= max_object_area
245	)
246	all_labels = object_masks_df["label"]	6✔
247	labels_to_keep = all_labels * keep_labels_bool	6✔
248
249	# Map filtered objects to the object mask
250	objects_filtered_by_size = map_array(	6✔
251	labeled_object_masks, all_labels.to_numpy(), labels_to_keep.to_numpy()
252	)
253
254	return objects_filtered_by_size.astype(np.int32)	6✔
255
256
257	def get_block_size(block_type: str, fov_dim: int, img_shape: int) -> int:	6✔
258	"""
259	Computes the approximate local otsu threshold based on fov size (in μm) and pixel resolution.
260
261	Args:
262	block_type (str): Either "small_holes" or "local_thresh"
263	fov_dim (int, optional): The size in μm for the FOV.
264	img_shape (int, optional): The shape of the image.
265
266	Returns:
267	int: Returns the approximate block area
268	"""
269
270	# Input validation
271	misc_utils.verify_in_list(	6✔
272	block_type=[block_type], block_types=["small_holes", "local_thresh"]
273	)
274	# Get the size of the pixel
275
276	pixel_size = fov_dim / img_shape	6✔
277	# grab block size for removing small holes
278	if block_type == "small_holes":	6✔
279	size = (np.pi * 5) ** 2 / pixel_size	6✔
280	# round up above value to make it into an integer
281	area: int = round(size)	6✔
282	# grab local threshold block size
283	else:
284	# use this to calculate out how many pixels it takes to get to roughly 10 μm
285	# (roughly a cell soma diameter)
286	size: float = 10 / pixel_size	6✔
287
288	# round the area up to the nearest odd number
289	area: int = round(size)	6✔
290	if area % 2 == 0:	6✔
UNCOV 291	area += 1	×
292	return area	6✔

angelolab / ark-analysis / 7105148444

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous