import os, cv2, time, sqlite3, traceback, shutil
import numpy as np
import pandas as pd
from collections import defaultdict
from scipy.stats import pearsonr, skew, kurtosis, mode
import multiprocessing as mp
from scipy.ndimage import distance_transform_edt, generate_binary_structure, binary_dilation, gaussian_filter, center_of_mass
from skimage.measure import regionprops, regionprops_table, shannon_entropy
from skimage.exposure import rescale_intensity
from skimage.segmentation import find_boundaries
from skimage.feature import graycomatrix, graycoprops
from mahotas.features import zernike_moments
from skimage import morphology, measure, filters
from skimage.util import img_as_bool
from skimage.filters import gaussian, threshold_otsu
from skimage.measure import label as sk_label, regionprops_table
import matplotlib.pyplot as plt
from math import ceil, sqrt
#from spacr.build.lib.spacr.utils import _filter_object
#from spacr.spacr import settings
[docs]
def get_components(cell_mask, nucleus_mask, pathogen_mask):
"""
Get the components (nucleus and pathogens) for each cell in the given masks.
Args:
cell_mask (ndarray): Binary mask of cell labels.
nucleus_mask (ndarray): Binary mask of nucleus labels.
pathogen_mask (ndarray): Binary mask of pathogen labels.
Returns:
tuple: A tuple containing two dataframes - nucleus_df and pathogen_df.
nucleus_df (DataFrame): Dataframe with columns 'cell_id' and 'nucleus',
representing the mapping of each cell to its nucleus.
pathogen_df (DataFrame): Dataframe with columns 'cell_id' and 'pathogen',
representing the mapping of each cell to its pathogens.
"""
# Create mappings from each cell to its nucleus, pathogens, and cytoplasms
cell_to_nucleus = defaultdict(list)
cell_to_pathogen = defaultdict(list)
# Get unique cell labels
cell_labels = np.unique(cell_mask)
# Iterate over each cell label
for cell_id in cell_labels:
if cell_id == 0:
continue
# Find corresponding component labels
nucleus_ids = np.unique(nucleus_mask[cell_mask == cell_id])
pathogen_ids = np.unique(pathogen_mask[cell_mask == cell_id])
# Update dictionaries, ignoring 0 (background) labels
cell_to_nucleus[cell_id] = nucleus_ids[nucleus_ids != 0].tolist()
cell_to_pathogen[cell_id] = pathogen_ids[pathogen_ids != 0].tolist()
# Convert dictionaries to dataframes
nucleus_df = pd.DataFrame(list(cell_to_nucleus.items()), columns=['cell_id', 'nucleus'])
pathogen_df = pd.DataFrame(list(cell_to_pathogen.items()), columns=['cell_id', 'pathogen'])
# Explode lists
nucleus_df = nucleus_df.explode('nucleus')
pathogen_df = pathogen_df.explode('pathogen')
return nucleus_df, pathogen_df
def _calculate_zernike(mask, df, degree=8):
"""
Calculate Zernike moments for each region in the given mask image.
Args:
mask (ndarray): Binary mask image.
df (DataFrame): Input DataFrame.
degree (int, optional): Degree of Zernike moments. Defaults to 8.
Returns:
DataFrame: Updated DataFrame with Zernike features appended, if any regions are found in the mask.
Otherwise, returns the original DataFrame.
Raises:
ValueError: If the lengths of Zernike moments are not consistent.
"""
zernike_features = []
for region in regionprops(mask):
zernike_moment = zernike_moments(region.image, degree)
zernike_features.append(zernike_moment.tolist())
if zernike_features:
feature_length = len(zernike_features[0])
for feature in zernike_features:
if len(feature) != feature_length:
raise ValueError("All Zernike moments must be of the same length")
zernike_df = pd.DataFrame(zernike_features, columns=[f'zernike_{i}' for i in range(feature_length)])
return pd.concat([df.reset_index(drop=True), zernike_df], axis=1)
else:
return df
def _analyze_cytoskeleton(array, mask, channel):
"""
Analyzes and extracts skeleton properties from labeled objects in a masked image based on microtubule staining intensities.
Parameters:
image : numpy array
Intensity image where the microtubules are stained.
mask : numpy array
Mask where objects are labeled for analysis. Each label corresponds to a unique object.
Returns:
DataFrame
A pandas DataFrame containing the measured properties of each object's skeleton.
"""
image = array[:, :, channel]
properties_list = []
# Process each object in the mask based on its label
for label in np.unique(mask):
if label == 0:
continue # Skip background
# Isolate the object using the label
object_region = mask == label
region_intensity = np.where(object_region, image, 0) # Use np.where for more efficient masking
# Ensure there are non-zero values to process
if np.any(region_intensity):
# Calculate adaptive offset based on intensity percentiles within the object
valid_pixels = region_intensity[region_intensity > 0]
if len(valid_pixels) > 1: # Ensure there are enough pixels to compute percentiles
offset = np.percentile(valid_pixels, 90) - np.percentile(valid_pixels, 50)
block_size = 35 # Adjust this based on your object sizes and detail needs
local_thresh = filters.threshold_local(region_intensity, block_size=block_size, offset=offset)
cytoskeleton = region_intensity > local_thresh
# Skeletonize the thresholded cytoskeleton
skeleton = morphology.skeletonize(img_as_bool(cytoskeleton))
# Measure properties of the skeleton
skeleton_props = measure.regionprops(measure.label(skeleton), intensity_image=image)
skeleton_length = sum(prop.area for prop in skeleton_props) # Sum of lengths of all skeleton segments
branch_data = morphology.skeleton_branch_analysis(skeleton)
# Store properties
properties = {
"object_label": label,
"skeleton_length": skeleton_length,
"skeleton_branch_points": len(branch_data['branch_points'])
}
properties_list.append(properties)
else:
# Handle cases with insufficient pixels
properties_list.append({
"object_label": label,
"skeleton_length": 0,
"skeleton_branch_points": 0
})
return pd.DataFrame(properties_list)
def _intensity_measurements_v1(cell_mask, nucleus_mask, pathogen_mask, cytoplasm_mask, channel_arrays, settings, sizes=[3, 6, 12, 24], periphery=True, outside=True):
"""
Calculate various intensity measurements for different regions in the image.
Args:
cell_mask (ndarray): Binary mask indicating the cell regions.
nucleus_mask (ndarray): Binary mask indicating the nucleus regions.
pathogen_mask (ndarray): Binary mask indicating the pathogen regions.
cytoplasm_mask (ndarray): Binary mask indicating the cytoplasm regions.
channel_arrays (ndarray): Array of channel images.
settings (dict): Additional settings for the intensity measurements.
sizes (list, optional): List of sizes for the measurements. Defaults to [3, 6, 12, 24].
periphery (bool, optional): Flag indicating whether to calculate periphery intensity measurements. Defaults to True.
outside (bool, optional): Flag indicating whether to calculate outside intensity measurements. Defaults to True.
Returns:
dict: A dictionary containing the calculated intensity measurements.
"""
radial_dist = settings['radial_dist']
calculate_correlation = settings['calculate_correlation']
homogeneity = settings['homogeneity']
distances = settings['homogeneity_distances']
intensity_props = ["label", "centroid_weighted", "centroid_weighted_local", "max_intensity", "mean_intensity", "min_intensity"]
col_lables = ['region_label', 'mean', '5_percentile', '10_percentile', '25_percentile', '50_percentile', '75_percentile', '85_percentile', '95_percentile']
cell_dfs, nucleus_dfs, pathogen_dfs, cytoplasm_dfs = [], [], [], []
ls = ['cell','nucleus','pathogen','cytoplasm']
labels = [cell_mask, nucleus_mask, pathogen_mask, cytoplasm_mask]
dfs = [cell_dfs, nucleus_dfs, pathogen_dfs, cytoplasm_dfs]
for i in range(0,channel_arrays.shape[-1]):
channel = channel_arrays[:, :, i]
for j, (label, df) in enumerate(zip(labels, dfs)):
if np.max(label) == 0:
empty_df = pd.DataFrame()
df.append(empty_df)
continue
mask_intensity_df = _extended_regionprops_table(label, channel, intensity_props)
#mask_intensity_df['shannon_entropy'] = shannon_entropy(channel, base=2)
if homogeneity:
homogeneity_df = _calculate_homogeneity(label, channel, distances)
mask_intensity_df = pd.concat([mask_intensity_df.reset_index(drop=True), homogeneity_df], axis=1)
if periphery:
if ls[j] == 'nucleus' or ls[j] == 'pathogen':
periphery_intensity_stats = _periphery_intensity(label, channel)
mask_intensity_df = pd.concat([mask_intensity_df, pd.DataFrame(periphery_intensity_stats, columns=[f'periphery_{stat}' for stat in col_lables])],axis=1)
if outside:
if ls[j] == 'nucleus' or ls[j] == 'pathogen':
outside_intensity_stats = _outside_intensity(label, channel)
mask_intensity_df = pd.concat([mask_intensity_df, pd.DataFrame(outside_intensity_stats, columns=[f'outside_{stat}' for stat in col_lables])], axis=1)
blur_col = [_estimate_blur(channel[label == region_label]) for region_label in mask_intensity_df['label']]
mask_intensity_df[f'{ls[j]}_channel_{i}_blur'] = blur_col
mask_intensity_df.columns = [f'{ls[j]}_channel_{i}_{col}' if col != 'label' else col for col in mask_intensity_df.columns]
df.append(mask_intensity_df)
if isinstance(settings['distance_gaussian_sigma'], int):
if settings['distance_gaussian_sigma'] != 0:
if settings['cell_mask_dim'] != None:
if settings['nucleus_mask_dim'] != None or settings['pathogen_mask_dim'] != None:
intensity_distance_df = _measure_intensity_distance(cell_mask, nucleus_mask, pathogen_mask, channel_arrays, settings)
cell_dfs.append(intensity_distance_df)
if radial_dist:
if np.max(nucleus_mask) != 0:
nucleus_radial_distributions = _calculate_radial_distribution(cell_mask, nucleus_mask, channel_arrays, num_bins=6)
nucleus_df = _create_dataframe(nucleus_radial_distributions, 'nucleus')
dfs[1].append(nucleus_df)
if np.max(pathogen_mask) != 0:
pathogen_radial_distributions = _calculate_radial_distribution(cell_mask, pathogen_mask, channel_arrays, num_bins=6)
pathogen_df = _create_dataframe(pathogen_radial_distributions, 'pathogen')
dfs[2].append(pathogen_df)
if calculate_correlation:
if channel_arrays.shape[-1] >= 2:
for i in range(channel_arrays.shape[-1]):
for j in range(i+1, channel_arrays.shape[-1]):
chan_i = channel_arrays[:, :, i]
chan_j = channel_arrays[:, :, j]
for m, mask in enumerate(labels):
coloc_df = _calculate_correlation_object_level(chan_i, chan_j, mask, settings)
coloc_df.columns = [f'{ls[m]}_channel_{i}_channel_{j}_{col}' for col in coloc_df.columns]
dfs[m].append(coloc_df)
return pd.concat(cell_dfs, axis=1), pd.concat(nucleus_dfs, axis=1), pd.concat(pathogen_dfs, axis=1), pd.concat(cytoplasm_dfs, axis=1)
def _morphological_measurements(cell_mask, nucleus_mask, pathogen_mask, organelle_mask, cytoplasm_mask, settings, zernike=True, degree=8):
"""
Calculate morphological measurements for cells, nucleus, pathogens, organelles, and cytoplasms.
Args:
cell_mask (ndarray): Binary mask of cell labels.
nucleus_mask (ndarray): Binary mask of nucleus labels.
pathogen_mask (ndarray): Binary mask of pathogen labels.
organelle_mask (ndarray): Binary mask of organelle labels.
cytoplasm_mask (ndarray): Binary mask of cytoplasm labels.
settings (dict): Dictionary containing settings for the measurements.
zernike (bool, optional): Flag indicating whether to calculate Zernike moments. Defaults to True.
degree (int, optional): Degree of Zernike moments. Defaults to 8.
Returns:
tuple: (cell_df, nucleus_df, pathogen_df, organelle_df, cytoplasm_df)
"""
morphological_props = ['label', 'area', 'area_filled', 'area_bbox', 'convex_area', 'major_axis_length', 'minor_axis_length',
'eccentricity', 'solidity', 'extent', 'perimeter', 'euler_number', 'equivalent_diameter_area', 'feret_diameter_max']
prop_ls = []
ls = []
if settings['cell_mask_dim'] is not None:
cell_to_nucleus, cell_to_pathogen = get_components(cell_mask, nucleus_mask, pathogen_mask)
cell_props = pd.DataFrame(regionprops_table(cell_mask, properties=morphological_props))
cell_props = _calculate_zernike(cell_mask, cell_props, degree=degree)
prop_ls.append(cell_props)
ls.append('cell')
else:
prop_ls.append(pd.DataFrame())
ls.append('cell')
if settings['nucleus_mask_dim'] is not None:
nucleus_props = pd.DataFrame(regionprops_table(nucleus_mask, properties=morphological_props))
nucleus_props = _calculate_zernike(nucleus_mask, nucleus_props, degree=degree)
if settings['cell_mask_dim'] is not None:
nucleus_props = pd.merge(nucleus_props, cell_to_nucleus, left_on='label', right_on='nucleus', how='left')
prop_ls.append(nucleus_props)
ls.append('nucleus')
else:
prop_ls.append(pd.DataFrame())
ls.append('nucleus')
if settings['pathogen_mask_dim'] is not None:
pathogen_props = pd.DataFrame(regionprops_table(pathogen_mask, properties=morphological_props))
pathogen_props = _calculate_zernike(pathogen_mask, pathogen_props, degree=degree)
if settings['cell_mask_dim'] is not None:
pathogen_props = pd.merge(pathogen_props, cell_to_pathogen, left_on='label', right_on='pathogen', how='left')
prop_ls.append(pathogen_props)
ls.append('pathogen')
else:
prop_ls.append(pd.DataFrame())
ls.append('pathogen')
if settings.get('organelle_mask_dim') is not None:
organelle_props = pd.DataFrame(regionprops_table(organelle_mask, properties=morphological_props))
if len(organelle_props) > 0:
organelle_props = _calculate_zernike(organelle_mask, organelle_props, degree=degree)
# Map each organelle to its parent cell
if settings['cell_mask_dim'] is not None:
organelle_to_cell = _map_child_to_parent(organelle_mask, cell_mask, child_name='organelle', parent_name='cell')
organelle_props = pd.merge(organelle_props, organelle_to_cell, left_on='label', right_on='organelle', how='left')
prop_ls.append(organelle_props)
ls.append('organelle')
else:
prop_ls.append(pd.DataFrame())
ls.append('organelle')
if settings['cytoplasm']:
cytoplasm_props = pd.DataFrame(regionprops_table(cytoplasm_mask, properties=morphological_props))
prop_ls.append(cytoplasm_props)
ls.append('cytoplasm')
else:
prop_ls.append(pd.DataFrame())
ls.append('cytoplasm')
df_ls = []
for i, df in enumerate(prop_ls):
df.columns = [f'{ls[i]}_{col}' for col in df.columns]
df = df.rename(columns={col: 'label' for col in df.columns if 'label' in col})
df_ls.append(df)
return df_ls[0], df_ls[1], df_ls[2], df_ls[3], df_ls[4]
def _map_child_to_parent(child_mask, parent_mask, child_name='organelle', parent_name='cell'):
"""
Map each child object to its parent by maximum overlap.
Returns a DataFrame with columns [child_name, parent_name].
"""
child_labels = np.unique(child_mask)
child_labels = child_labels[child_labels != 0]
mapping = []
for child_id in child_labels:
region = child_mask == child_id
parent_ids = parent_mask[region]
parent_ids = parent_ids[parent_ids != 0]
if len(parent_ids) > 0:
parent_id = np.bincount(parent_ids).argmax()
else:
parent_id = 0
mapping.append({child_name: child_id, parent_name: parent_id})
return pd.DataFrame(mapping)
def _summarize_organelles_per_parent(organelle_mask, parent_mask, channel_arrays, parent_name='cell'):
"""
Summarize organelle measurements per parent object (e.g. per cell).
For each parent object, computes:
- organelle_count: number of organelles within the parent
- organelle_total_area: sum of organelle areas
- organelle_mean_area, organelle_std_area: mean/std of individual organelle areas
- organelle_fraction: fraction of parent area occupied by organelles
- organelle_mean_eccentricity, organelle_std_eccentricity
- organelle_mean_solidity, organelle_std_solidity
- per-channel: organelle_mean_intensity, organelle_std_intensity
Parameters
----------
organelle_mask : ndarray
Label mask for organelles.
parent_mask : ndarray
Label mask for parent objects (cells, nuclei, etc.).
channel_arrays : ndarray
Shape (H, W, C) intensity images.
parent_name : str
Name for the parent object column.
Returns
-------
pd.DataFrame
One row per parent object with summarised organelle features.
"""
parent_labels = np.unique(parent_mask)
parent_labels = parent_labels[parent_labels != 0]
morphological_props = ['label', 'area', 'eccentricity', 'solidity', 'major_axis_length', 'minor_axis_length']
# Get per-organelle morphology
organelle_props = regionprops_table(organelle_mask, properties=morphological_props)
organelle_df = pd.DataFrame(organelle_props)
# Map each organelle to its parent
organelle_to_parent = _map_child_to_parent(organelle_mask, parent_mask,
child_name='organelle_label',
parent_name=parent_name)
if len(organelle_df) > 0 and len(organelle_to_parent) > 0:
organelle_df = pd.merge(organelle_df, organelle_to_parent,
left_on='label', right_on='organelle_label', how='left')
else:
# No organelles — return empty summary for all parents
rows = []
for pid in parent_labels:
row = {'label': pid, 'organelle_count': 0, 'organelle_total_area': 0,
'organelle_fraction': 0.0}
rows.append(row)
return pd.DataFrame(rows)
# Per-channel intensity per organelle
for ch in range(channel_arrays.shape[-1]):
channel = channel_arrays[:, :, ch]
intensities = []
for org_label in organelle_df['label']:
region = organelle_mask == org_label
if np.any(region):
intensities.append(channel[region].mean())
else:
intensities.append(0.0)
organelle_df[f'organelle_ch{ch}_mean_intensity'] = intensities
# Get parent areas for fraction calculation
parent_props = pd.DataFrame(regionprops_table(parent_mask, properties=['label', 'area']))
parent_area_map = dict(zip(parent_props['label'], parent_props['area']))
# Summarise per parent
summary_rows = []
for pid in parent_labels:
org_subset = organelle_df[organelle_df[parent_name] == pid]
parent_area = parent_area_map.get(pid, 1)
row = {'label': pid}
row['organelle_count'] = len(org_subset)
row['organelle_total_area'] = org_subset['area'].sum() if len(org_subset) > 0 else 0
row['organelle_fraction'] = row['organelle_total_area'] / parent_area if parent_area > 0 else 0.0
row['organelle_mean_area'] = org_subset['area'].mean() if len(org_subset) > 0 else 0.0
row['organelle_std_area'] = org_subset['area'].std() if len(org_subset) > 1 else 0.0
row['organelle_mean_eccentricity'] = org_subset['eccentricity'].mean() if len(org_subset) > 0 else 0.0
row['organelle_std_eccentricity'] = org_subset['eccentricity'].std() if len(org_subset) > 1 else 0.0
row['organelle_mean_solidity'] = org_subset['solidity'].mean() if len(org_subset) > 0 else 0.0
row['organelle_std_solidity'] = org_subset['solidity'].std() if len(org_subset) > 1 else 0.0
row['organelle_mean_major_axis'] = org_subset['major_axis_length'].mean() if len(org_subset) > 0 else 0.0
row['organelle_mean_minor_axis'] = org_subset['minor_axis_length'].mean() if len(org_subset) > 0 else 0.0
for ch in range(channel_arrays.shape[-1]):
col = f'organelle_ch{ch}_mean_intensity'
row[f'organelle_ch{ch}_mean_intensity_per_{parent_name}'] = org_subset[col].mean() if len(org_subset) > 0 else 0.0
row[f'organelle_ch{ch}_std_intensity_per_{parent_name}'] = org_subset[col].std() if len(org_subset) > 1 else 0.0
summary_rows.append(row)
return pd.DataFrame(summary_rows)
def _intensity_measurements(cell_mask, nucleus_mask, pathogen_mask, organelle_mask, cytoplasm_mask, channel_arrays, settings, sizes=[3, 6, 12, 24], periphery=True, outside=True):
"""
Calculate various intensity measurements for different regions in the image,
including organelle masks.
Args:
cell_mask (ndarray): Binary mask indicating the cell regions.
nucleus_mask (ndarray): Binary mask indicating the nucleus regions.
pathogen_mask (ndarray): Binary mask indicating the pathogen regions.
organelle_mask (ndarray): Binary mask indicating the organelle regions.
cytoplasm_mask (ndarray): Binary mask indicating the cytoplasm regions.
channel_arrays (ndarray): Array of channel images.
settings (dict): Additional settings for the intensity measurements.
sizes (list, optional): List of sizes for the measurements.
periphery (bool, optional): Calculate periphery intensity measurements.
outside (bool, optional): Calculate outside intensity measurements.
Returns:
tuple: (cell_intensity_df, nucleus_intensity_df, pathogen_intensity_df,
organelle_intensity_df, cytoplasm_intensity_df)
"""
radial_dist = settings['radial_dist']
calculate_correlation = settings['calculate_correlation']
homogeneity = settings['homogeneity']
distances = settings['homogeneity_distances']
intensity_props = ["label", "centroid_weighted", "centroid_weighted_local", "max_intensity", "mean_intensity", "min_intensity"]
col_lables = ['region_label', 'mean', '5_percentile', '10_percentile', '25_percentile', '50_percentile', '75_percentile', '85_percentile', '95_percentile']
cell_dfs, nucleus_dfs, pathogen_dfs, organelle_dfs, cytoplasm_dfs = [], [], [], [], []
ls = ['cell', 'nucleus', 'pathogen', 'organelle', 'cytoplasm']
labels = [cell_mask, nucleus_mask, pathogen_mask, organelle_mask, cytoplasm_mask]
dfs = [cell_dfs, nucleus_dfs, pathogen_dfs, organelle_dfs, cytoplasm_dfs]
for i in range(0, channel_arrays.shape[-1]):
channel = channel_arrays[:, :, i]
for j, (label, df) in enumerate(zip(labels, dfs)):
if np.max(label) == 0:
empty_df = pd.DataFrame()
df.append(empty_df)
continue
mask_intensity_df = _extended_regionprops_table(label, channel, intensity_props)
if homogeneity:
homogeneity_df = _calculate_homogeneity(label, channel, distances)
mask_intensity_df = pd.concat([mask_intensity_df.reset_index(drop=True), homogeneity_df], axis=1)
if periphery:
if ls[j] in ('nucleus', 'pathogen', 'organelle'):
periphery_intensity_stats = _periphery_intensity(label, channel)
mask_intensity_df = pd.concat([mask_intensity_df, pd.DataFrame(periphery_intensity_stats, columns=[f'periphery_{stat}' for stat in col_lables])], axis=1)
if outside:
if ls[j] in ('nucleus', 'pathogen', 'organelle'):
outside_intensity_stats = _outside_intensity(label, channel)
mask_intensity_df = pd.concat([mask_intensity_df, pd.DataFrame(outside_intensity_stats, columns=[f'outside_{stat}' for stat in col_lables])], axis=1)
blur_col = [_estimate_blur(channel[label == region_label]) for region_label in mask_intensity_df['label']]
mask_intensity_df[f'{ls[j]}_channel_{i}_blur'] = blur_col
mask_intensity_df.columns = [f'{ls[j]}_channel_{i}_{col}' if col != 'label' else col for col in mask_intensity_df.columns]
df.append(mask_intensity_df)
if isinstance(settings['distance_gaussian_sigma'], int):
if settings['distance_gaussian_sigma'] != 0:
if settings['cell_mask_dim'] is not None:
if settings['nucleus_mask_dim'] is not None or settings['pathogen_mask_dim'] is not None:
intensity_distance_df = _measure_intensity_distance(cell_mask, nucleus_mask, pathogen_mask, channel_arrays, settings)
cell_dfs.append(intensity_distance_df)
if radial_dist:
if np.max(nucleus_mask) != 0:
nucleus_radial_distributions = _calculate_radial_distribution(cell_mask, nucleus_mask, channel_arrays, num_bins=6)
nucleus_df = _create_dataframe(nucleus_radial_distributions, 'nucleus')
dfs[1].append(nucleus_df)
if np.max(pathogen_mask) != 0:
pathogen_radial_distributions = _calculate_radial_distribution(cell_mask, pathogen_mask, channel_arrays, num_bins=6)
pathogen_df = _create_dataframe(pathogen_radial_distributions, 'pathogen')
dfs[2].append(pathogen_df)
if np.max(organelle_mask) != 0:
organelle_radial_distributions = _calculate_radial_distribution(cell_mask, organelle_mask, channel_arrays, num_bins=6)
organelle_rad_df = _create_dataframe(organelle_radial_distributions, 'organelle')
dfs[3].append(organelle_rad_df)
if calculate_correlation:
if channel_arrays.shape[-1] >= 2:
for i in range(channel_arrays.shape[-1]):
for j in range(i+1, channel_arrays.shape[-1]):
chan_i = channel_arrays[:, :, i]
chan_j = channel_arrays[:, :, j]
for m, mask in enumerate(labels):
coloc_df = _calculate_correlation_object_level(chan_i, chan_j, mask, settings)
coloc_df.columns = [f'{ls[m]}_channel_{i}_channel_{j}_{col}' for col in coloc_df.columns]
dfs[m].append(coloc_df)
return (pd.concat(cell_dfs, axis=1),
pd.concat(nucleus_dfs, axis=1),
pd.concat(pathogen_dfs, axis=1),
pd.concat(organelle_dfs, axis=1),
pd.concat(cytoplasm_dfs, axis=1))
def _morphological_measurements_v1(cell_mask, nucleus_mask, pathogen_mask, cytoplasm_mask, settings, zernike=True, degree=8):
"""
Calculate morphological measurements for cells, nucleus, pathogens, and cytoplasms based on the given masks.
Args:
cell_mask (ndarray): Binary mask of cell labels.
nucleus_mask (ndarray): Binary mask of nucleus labels.
pathogen_mask (ndarray): Binary mask of pathogen labels.
cytoplasm_mask (ndarray): Binary mask of cytoplasm labels.
settings (dict): Dictionary containing settings for the measurements.
zernike (bool, optional): Flag indicating whether to calculate Zernike moments. Defaults to True.
degree (int, optional): Degree of Zernike moments. Defaults to 8.
Returns:
tuple: A tuple containing four dataframes - cell_df, nucleus_df, pathogen_df, and cytoplasm_df.
cell_df (DataFrame): Dataframe with morphological measurements for cells.
nucleus_df (DataFrame): Dataframe with morphological measurements for nucleus.
pathogen_df (DataFrame): Dataframe with morphological measurements for pathogens.
cytoplasm_df (DataFrame): Dataframe with morphological measurements for cytoplasms.
"""
morphological_props = ['label', 'area', 'area_filled', 'area_bbox', 'convex_area', 'major_axis_length', 'minor_axis_length',
'eccentricity', 'solidity', 'extent', 'perimeter', 'euler_number', 'equivalent_diameter_area', 'feret_diameter_max']
prop_ls = []
ls = []
# Create mappings from each cell to its nucleus, pathogens, and cytoplasms
if settings['cell_mask_dim'] is not None:
cell_to_nucleus, cell_to_pathogen = get_components(cell_mask, nucleus_mask, pathogen_mask)
cell_props = pd.DataFrame(regionprops_table(cell_mask, properties=morphological_props))
cell_props = _calculate_zernike(cell_mask, cell_props, degree=degree)
prop_ls = prop_ls + [cell_props]
ls = ls + ['cell']
else:
prop_ls = prop_ls + [pd.DataFrame()]
ls = ls + ['cell']
if settings['nucleus_mask_dim'] is not None:
nucleus_props = pd.DataFrame(regionprops_table(nucleus_mask, properties=morphological_props))
nucleus_props = _calculate_zernike(nucleus_mask, nucleus_props, degree=degree)
if settings['cell_mask_dim'] is not None:
nucleus_props = pd.merge(nucleus_props, cell_to_nucleus, left_on='label', right_on='nucleus', how='left')
prop_ls = prop_ls + [nucleus_props]
ls = ls + ['nucleus']
else:
prop_ls = prop_ls + [pd.DataFrame()]
ls = ls + ['nucleus']
if settings['pathogen_mask_dim'] is not None:
pathogen_props = pd.DataFrame(regionprops_table(pathogen_mask, properties=morphological_props))
pathogen_props = _calculate_zernike(pathogen_mask, pathogen_props, degree=degree)
if settings['cell_mask_dim'] is not None:
pathogen_props = pd.merge(pathogen_props, cell_to_pathogen, left_on='label', right_on='pathogen', how='left')
prop_ls = prop_ls + [pathogen_props]
ls = ls + ['pathogen']
else:
prop_ls = prop_ls + [pd.DataFrame()]
ls = ls + ['pathogen']
if settings['cytoplasm']:
cytoplasm_props = pd.DataFrame(regionprops_table(cytoplasm_mask, properties=morphological_props))
prop_ls = prop_ls + [cytoplasm_props]
ls = ls + ['cytoplasm']
else:
prop_ls = prop_ls + [pd.DataFrame()]
ls = ls + ['cytoplasm']
df_ls = []
for i,df in enumerate(prop_ls):
df.columns = [f'{ls[i]}_{col}' for col in df.columns]
df = df.rename(columns={col: 'label' for col in df.columns if 'label' in col})
df_ls.append(df)
return df_ls[0], df_ls[1], df_ls[2], df_ls[3]
def _create_dataframe(radial_distributions, object_type):
"""
Create a pandas DataFrame from the given radial distributions.
Parameters:
- radial_distributions (dict): A dictionary containing the radial distributions.
- object_type (str): The type of object.
Returns:
- df (pandas.DataFrame): The created DataFrame.
"""
df = pd.DataFrame()
for key, value in radial_distributions.items():
cell_label, object_label, channel_index = key
for i in range(len(value)):
col_name = f'{object_type}_rad_dist_channel_{channel_index}_bin_{i}'
df.loc[object_label, col_name] = value[i]
df.loc[object_label, 'cell_id'] = cell_label
# Reset the index and rename the column that was previously the index
df = df.reset_index().rename(columns={'index': 'label'})
return df
def _extended_regionprops_table(labels, image, intensity_props):
"""
Calculate extended region properties table, adding a suite of advanced quantitative features.
"""
def _gini(array):
# Compute Gini coefficient (nan safe)
array = np.abs(array[~np.isnan(array)])
n = array.size
if n == 0:
return np.nan
array = np.sort(array)
index = np.arange(1, n + 1)
return (np.sum((2 * index - n - 1) * array)) / (n * np.sum(array)) if np.sum(array) else np.nan
props = regionprops_table(labels, image, properties=intensity_props)
df = pd.DataFrame(props)
regions = regionprops(labels, intensity_image=image)
integrated_intensity = []
std_intensity = []
median_intensity = []
skew_intensity = []
kurtosis_intensity = []
mode_intensity = []
range_intensity = []
iqr_intensity = []
cv_intensity = []
gini_intensity = []
frac_high90 = []
frac_low10 = []
entropy_intensity = []
for region in regions:
intens = region.intensity_image[region.image]
intens = intens[~np.isnan(intens)]
if intens.size == 0:
integrated_intensity.append(np.nan)
std_intensity.append(np.nan)
median_intensity.append(np.nan)
skew_intensity.append(np.nan)
kurtosis_intensity.append(np.nan)
mode_intensity.append(np.nan)
range_intensity.append(np.nan)
iqr_intensity.append(np.nan)
cv_intensity.append(np.nan)
gini_intensity.append(np.nan)
frac_high90.append(np.nan)
frac_low10.append(np.nan)
entropy_intensity.append(np.nan)
else:
integrated_intensity.append(np.sum(intens))
std_intensity.append(np.std(intens))
median_intensity.append(np.median(intens))
skew_intensity.append(skew(intens) if intens.size > 2 else np.nan)
kurtosis_intensity.append(kurtosis(intens) if intens.size > 3 else np.nan)
# Mode (use first mode value if multimodal)
try:
mode_val = mode(intens, nan_policy='omit').mode
mode_intensity.append(mode_val[0] if len(mode_val) > 0 else np.nan)
except Exception:
mode_intensity.append(np.nan)
range_intensity.append(np.ptp(intens))
iqr_intensity.append(np.percentile(intens, 75) - np.percentile(intens, 25))
cv_intensity.append(np.std(intens) / np.mean(intens) if np.mean(intens) != 0 else np.nan)
gini_intensity.append(_gini(intens))
frac_high90.append(np.mean(intens > np.percentile(intens, 90)))
frac_low10.append(np.mean(intens < np.percentile(intens, 10)))
entropy_intensity.append(shannon_entropy(intens) if intens.size > 1 else 0.0)
df['integrated_intensity'] = integrated_intensity
df['std_intensity'] = std_intensity
df['median_intensity'] = median_intensity
df['skew_intensity'] = skew_intensity
df['kurtosis_intensity'] = kurtosis_intensity
df['mode_intensity'] = mode_intensity
df['range_intensity'] = range_intensity
df['iqr_intensity'] = iqr_intensity
df['cv_intensity'] = cv_intensity
df['gini_intensity'] = gini_intensity
df['frac_high90'] = frac_high90
df['frac_low10'] = frac_low10
df['entropy_intensity'] = entropy_intensity
percentiles = [5, 10, 25, 75, 85, 95]
for p in percentiles:
df[f'percentile_{p}'] = [
np.percentile(region.intensity_image[region.image], p)
for region in regions
]
return df
def _calculate_homogeneity(label, channel, distances=[2,4,8,16,32,64]):
"""
Calculate the homogeneity values for each region in the label mask.
Parameters:
- label (ndarray): The label mask containing the regions.
- channel (ndarray): The image channel corresponding to the label mask.
- distances (list): The distances to calculate the homogeneity for.
Returns:
- homogeneity_df (DataFrame): A DataFrame containing the homogeneity values for each region and distance.
"""
homogeneity_values = []
# Iterate through the regions in label_mask
for region in regionprops(label):
region_image = (region.image * channel[region.slice]).astype(int)
homogeneity_per_distance = []
for d in distances:
rescaled_image = rescale_intensity(region_image, out_range=(0, 255)).astype('uint8')
glcm = graycomatrix(rescaled_image, [d], [0], symmetric=True, normed=True)
homogeneity_per_distance.append(graycoprops(glcm, 'homogeneity')[0, 0])
homogeneity_values.append(homogeneity_per_distance)
columns = [f'homogeneity_distance_{d}' for d in distances]
homogeneity_df = pd.DataFrame(homogeneity_values, columns=columns)
return homogeneity_df
def _periphery_intensity(label_mask, image):
"""
Calculate intensity statistics for the periphery regions in the label mask.
Args:
label_mask (ndarray): Binary mask indicating the regions of interest.
image (ndarray): Input image.
Returns:
list: List of tuples containing periphery intensity statistics for each region.
Each tuple contains the region label and the following statistics:
- Mean intensity
- 5th percentile intensity
- 10th percentile intensity
- 25th percentile intensity
- 50th percentile intensity (median)
- 75th percentile intensity
- 85th percentile intensity
- 95th percentile intensity
"""
periphery_intensity_stats = []
boundary = find_boundaries(label_mask)
for region in np.unique(label_mask)[1:]: # skip the background label
region_boundary = boundary & (label_mask == region)
intensities = image[region_boundary]
if intensities.size == 0:
periphery_intensity_stats.append((region, np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, np.nan))
else:
periphery_intensity_stats.append((region, np.mean(intensities), np.percentile(intensities,5), np.percentile(intensities,10),
np.percentile(intensities,25), np.percentile(intensities,50),
np.percentile(intensities,75), np.percentile(intensities,85),
np.percentile(intensities,95)))
return periphery_intensity_stats
def _outside_intensity(label_mask, image, distance=5):
"""
Calculate the statistics of intensities outside each labeled region in the image.
Args:
label_mask (ndarray): Binary mask indicating the labeled regions.
image (ndarray): Input image.
distance (int): Distance for dilation operation (default: 5).
Returns:
list: List of tuples containing the statistics for each labeled region.
Each tuple contains the region label and the following statistics:
- Mean intensity
- 5th percentile intensity
- 10th percentile intensity
- 25th percentile intensity
- 50th percentile intensity (median)
- 75th percentile intensity
- 85th percentile intensity
- 95th percentile intensity
"""
outside_intensity_stats = []
for region in np.unique(label_mask)[1:]: # skip the background label
region_mask = label_mask == region
dilated_mask = binary_dilation(region_mask, iterations=distance)
outside_mask = dilated_mask & ~region_mask
intensities = image[outside_mask]
if intensities.size == 0:
outside_intensity_stats.append((region, np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, np.nan))
else:
outside_intensity_stats.append((region, np.mean(intensities), np.percentile(intensities,5), np.percentile(intensities,10),
np.percentile(intensities,25), np.percentile(intensities,50),
np.percentile(intensities,75), np.percentile(intensities,85),
np.percentile(intensities,95)))
return outside_intensity_stats
def _calculate_radial_distribution(cell_mask, object_mask, channel_arrays, num_bins=6):
"""
Calculate the radial distribution of average intensities for each object in each cell.
Args:
cell_mask (numpy.ndarray): The mask representing the cells.
object_mask (numpy.ndarray): The mask representing the objects.
channel_arrays (numpy.ndarray): The array of channel images.
num_bins (int, optional): The number of bins for the radial distribution. Defaults to 6.
Returns:
dict: A dictionary containing the radial distributions of average intensities for each object in each cell.
The keys are tuples of (cell_label, object_label, channel_index), and the values are numpy arrays
representing the radial distributions.
"""
def _calculate_average_intensity(distance_map, single_channel_image, num_bins):
"""
Calculate the average intensity of a single-channel image based on the distance map.
Args:
distance_map (numpy.ndarray): The distance map.
single_channel_image (numpy.ndarray): The single-channel image.
num_bins (int): The number of bins for the radial distribution.
Returns:
numpy.ndarray: The radial distribution of average intensities.
"""
radial_distribution = np.zeros(num_bins)
for i in range(num_bins):
min_distance = i * (distance_map.max() / num_bins)
max_distance = (i + 1) * (distance_map.max() / num_bins)
bin_mask = (distance_map >= min_distance) & (distance_map < max_distance)
radial_distribution[i] = single_channel_image[bin_mask].mean()
return radial_distribution
object_radial_distributions = {}
# get unique cell labels
cell_labels = np.unique(cell_mask)
cell_labels = cell_labels[cell_labels != 0]
for cell_label in cell_labels:
cell_region = cell_mask == cell_label
object_labels = np.unique(object_mask[cell_region])
object_labels = object_labels[object_labels != 0]
for object_label in object_labels:
objecyt_region = object_mask == object_label
object_boundary = find_boundaries(objecyt_region, mode='outer')
distance_map = distance_transform_edt(~object_boundary) * cell_region
for channel_index in range(channel_arrays.shape[2]):
radial_distribution = _calculate_average_intensity(distance_map, channel_arrays[:, :, channel_index], num_bins)
object_radial_distributions[(cell_label, object_label, channel_index)] = radial_distribution
return object_radial_distributions
def _calculate_correlation_object_level(channel_image1, channel_image2, mask, settings):
"""
Calculate correlation at the object level between two channel images based on a mask.
Args:
channel_image1 (numpy.ndarray): The first channel image.
channel_image2 (numpy.ndarray): The second channel image.
mask (numpy.ndarray): The mask indicating the objects.
settings (dict): Additional settings for correlation calculation.
Returns:
pandas.DataFrame: A DataFrame containing the correlation data at the object level.
"""
thresholds = settings['manders_thresholds']
corr_data = {}
for i in np.unique(mask)[1:]:
object_mask = (mask == i)
object_channel_image1 = channel_image1[object_mask]
object_channel_image2 = channel_image2[object_mask]
total_intensity1 = np.sum(object_channel_image1)
total_intensity2 = np.sum(object_channel_image2)
if len(object_channel_image1) < 2 or len(object_channel_image2) < 2:
pearson_corr = np.nan
else:
pearson_corr, _ = pearsonr(object_channel_image1, object_channel_image2)
corr_data[i] = {f'label_correlation': i,
f'Pearson_correlation': pearson_corr}
for thresh in thresholds:
chan1_thresh = np.percentile(object_channel_image1, thresh)
chan2_thresh = np.percentile(object_channel_image2, thresh)
# boolean mask where both signals are present
overlap_mask = (object_channel_image1 > chan1_thresh) & (object_channel_image2 > chan2_thresh)
M1 = np.sum(object_channel_image1[overlap_mask]) / total_intensity1 if total_intensity1 > 0 else 0
M2 = np.sum(object_channel_image2[overlap_mask]) / total_intensity2 if total_intensity2 > 0 else 0
corr_data[i].update({f'M1_correlation_{thresh}': M1,
f'M2_correlation_{thresh}': M2})
return pd.DataFrame(corr_data.values())
def _estimate_blur(image):
"""
Estimates the blur of an image by computing the variance of its Laplacian.
Parameters:
image (numpy.ndarray): The input image.
Returns:
float: The variance of the Laplacian of the image.
"""
# Check if the image is not already in a floating-point format
if image.dtype != np.float32 and image.dtype != np.float64:
# Convert the image to float64 for processing
image_float = image.astype(np.float64)
else:
# If it's already a floating-point image, use it as is
image_float = image
# Compute the Laplacian of the image
lap = cv2.Laplacian(image_float, cv2.CV_64F)
# Compute and return the variance of the Laplacian
return lap.var()
def _measure_intensity_distance(cell_mask, nucleus_mask, pathogen_mask, channel_arrays, settings):
"""
Compute Gaussian-smoothed intensity-weighted centroid distances for each cell object.
"""
sigma = settings.get('distance_gaussian_sigma', 1.0)
cell_labels = np.unique(cell_mask)
cell_labels = cell_labels[cell_labels > 0]
dfs = []
nucleus_dt = distance_transform_edt(nucleus_mask == 0)
pathogen_dt = distance_transform_edt(pathogen_mask == 0)
for ch in range(channel_arrays.shape[-1]):
channel_img = channel_arrays[:, :, ch]
blurred_img = gaussian_filter(channel_img, sigma=sigma)
data = []
for label in cell_labels:
cell_coords = np.argwhere(cell_mask == label)
if cell_coords.size == 0:
data.append([label, np.nan, np.nan])
continue
minr, minc = np.min(cell_coords, axis=0)
maxr, maxc = np.max(cell_coords, axis=0) + 1
cell_submask = (cell_mask[minr:maxr, minc:maxc] == label)
blurred_subimg = blurred_img[minr:maxr, minc:maxc]
if np.sum(cell_submask) == 0:
data.append([label, np.nan, np.nan])
continue
masked_intensity = blurred_subimg * cell_submask
com_local = center_of_mass(masked_intensity)
if np.isnan(com_local[0]):
data.append([label, np.nan, np.nan])
continue
com_global = (com_local[0] + minr, com_local[1] + minc)
com_global_int = tuple(np.round(com_global).astype(int))
x, y = com_global_int
if not (0 <= x < cell_mask.shape[0] and 0 <= y < cell_mask.shape[1]):
data.append([label, np.nan, np.nan])
continue
nucleus_dist = nucleus_dt[x, y]
pathogen_dist = pathogen_dt[x, y]
data.append([label, nucleus_dist, pathogen_dist])
df = pd.DataFrame(data, columns=['label',
f'cell_channel_{ch}_distance_to_nucleus',
f'cell_channel_{ch}_distance_to_pathogen'])
dfs.append(df)
# Merge all channel dataframes on label
merged_df = dfs[0]
for df in dfs[1:]:
merged_df = merged_df.merge(df, on='label', how='outer')
return merged_df
[docs]
def save_and_add_image_to_grid(png_channels, img_path, grid, plot=False):
"""
Add an image to a grid and save it as PNG.
Args:
png_channels (ndarray): The array representing the image channels.
img_path (str): The path to save the image as PNG.
grid (list): The grid of images to be plotted later.
Returns:
grid (list): Updated grid with the new image added.
"""
# Save the image as a PNG
cv2.imwrite(img_path, png_channels)
if plot:
# Ensure the image is in uint8 format for cv2 functions
if png_channels.dtype == np.uint16:
png_channels = (png_channels / 256).astype(np.uint8)
# Get the filename without the extension
filename = os.path.splitext(os.path.basename(img_path))[0]
# Add the label to the image
#labeled_image = cv2.putText(png_channels.copy(), filename, (10, 30),
# cv2.FONT_HERSHEY_SIMPLEX, 1, (255, 255, 255), 2, cv2.LINE_AA)
# Add the labeled image to the grid
grid.append(png_channels)
return grid
[docs]
def img_list_to_grid(grid, titles=None):
"""
Plot a grid of images with optional titles.
Args:
grid (list): List of images to be plotted.
titles (list): List of titles for the images.
Returns:
fig (Figure): The matplotlib figure object containing the image grid.
"""
n_images = len(grid)
grid_size = ceil(sqrt(n_images))
fig, axs = plt.subplots(grid_size, grid_size, figsize=(15, 15), facecolor='black')
for i, ax in enumerate(axs.flat):
if i < n_images:
image = grid[i]
ax.imshow(cv2.cvtColor(image, cv2.COLOR_BGR2RGB))
ax.axis('off')
ax.set_facecolor('black')
if titles:
# Determine text size
img_height, img_width = image.shape[:2]
text_size = max(min(img_width / (len(titles[i]) * 1.5), img_height / 10), 4)
ax.text(5, 5, titles[i], color='white', fontsize=text_size, ha='left', va='top', fontweight='bold')
else:
fig.delaxes(ax)
# Adjust spacing
plt.subplots_adjust(wspace=0.05, hspace=0.05)
plt.tight_layout(pad=0.1)
return fig
#@log_function_call
def _measure_crop_core(index, time_ls, file, settings):
"""
Measure and crop the images based on specified settings.
Parameters:
- index: int
The index of the image.
- time_ls: list
The list of time points.
- file: str
The file path of the image.
- settings: dict
The dictionary containing the settings for measurement and cropping.
Returns:
- cropped_images: list
A list of cropped images.
"""
from .plot import _plot_cropped_arrays
from .utils import _merge_overlapping_objects, _filter_object, _relabel_parent_with_child_labels, _exclude_objects, normalize_to_dtype, filepaths_to_database
from .utils import _merge_and_save_to_database, _crop_center, _find_bounding_box, _generate_names, _get_percentiles
figs = {}
grid = []
start = time.time()
try:
source_folder = os.path.dirname(settings['src'])
file_name = os.path.splitext(file)[0]
data = np.load(os.path.join(settings['src'], file))
data_type = data.dtype
if data_type not in ['uint8','uint16']:
data_type_before = data_type
data = data.astype(np.uint16)
data_type = data.dtype
if settings['verbose']:
print(f'Converted data from {data_type_before} to {data_type}')
if settings['plot']:
if len(data.shape) == 3:
figuresize = data.shape[2]*10
else:
figuresize = 10
fig = _plot_cropped_arrays(data, file, figuresize)
figs[f'{file_name}__before_filtration'] = fig
channel_arrays = data[:, :, settings['channels']].astype(data_type)
if settings['cell_mask_dim'] is not None:
cell_mask = data[:, :, settings['cell_mask_dim']].astype(data_type)
if settings['cell_min_size'] is not None and settings['cell_min_size'] != 0:
cell_mask = _filter_object(cell_mask, settings['cell_min_size'])
else:
cell_mask = np.zeros_like(data[:, :, 0])
settings['cytoplasm'] = False
settings['uninfected'] = True
if settings['nucleus_mask_dim'] is not None:
nucleus_mask = data[:, :, settings['nucleus_mask_dim']].astype(data_type)
if settings['cell_mask_dim'] is not None:
nucleus_mask, cell_mask = _merge_overlapping_objects(mask1=nucleus_mask, mask2=cell_mask)
if settings['nucleus_min_size'] is not None and settings['nucleus_min_size'] != 0:
nucleus_mask = _filter_object(nucleus_mask, settings['nucleus_min_size'])
if settings['timelapse_objects'] == 'nucleus':
if settings['cell_mask_dim'] is not None:
cell_mask, nucleus_mask = _relabel_parent_with_child_labels(cell_mask, nucleus_mask)
data[:, :, settings['cell_mask_dim']] = cell_mask
data[:, :, settings['nucleus_mask_dim']] = nucleus_mask
save_folder = settings['src']
np.save(os.path.join(save_folder, file), data)
else:
nucleus_mask = np.zeros_like(data[:, :, 0])
if settings['pathogen_mask_dim'] is not None:
pathogen_mask = data[:, :, settings['pathogen_mask_dim']].astype(data_type)
if settings['merge_edge_pathogen_cells']:
if settings['cell_mask_dim'] is not None:
pathogen_mask, cell_mask = _merge_overlapping_objects(mask1=pathogen_mask, mask2=cell_mask)
if settings['pathogen_min_size'] is not None and settings['pathogen_min_size'] != 0:
pathogen_mask = _filter_object(pathogen_mask, settings['pathogen_min_size'])
else:
pathogen_mask = np.zeros_like(data[:, :, 0])
if settings.get('organelle_mask_dim') is not None:
organelle_mask = data[:, :, settings['organelle_mask_dim']].astype(data_type)
if settings.get('organelle_min_size') and settings['organelle_min_size'] != 0:
organelle_mask = _filter_object(organelle_mask, settings['organelle_min_size'])
else:
organelle_mask = np.zeros_like(data[:, :, 0])
# Create cytoplasm mask
if settings['cytoplasm']:
if settings['cell_mask_dim'] is not None:
# Build a combined interior mask from all subcellular objects
interior = np.zeros_like(cell_mask, dtype=bool)
if settings['nucleus_mask_dim'] is not None:
interior |= (nucleus_mask != 0)
if settings['pathogen_mask_dim'] is not None:
interior |= (pathogen_mask != 0)
if settings.get('organelle_mask_dim') is not None:
interior |= (organelle_mask != 0)
cytoplasm_mask = np.where(interior, 0, cell_mask)
else:
cytoplasm_mask = np.zeros_like(cell_mask)
else:
cytoplasm_mask = np.zeros_like(cell_mask)
if settings['cell_min_size'] is not None and settings['cell_min_size'] != 0:
cell_mask = _filter_object(cell_mask, settings['cell_min_size'])
if settings['nucleus_min_size'] is not None and settings['nucleus_min_size'] != 0:
nucleus_mask = _filter_object(nucleus_mask, settings['nucleus_min_size'])
if settings['pathogen_min_size'] is not None and settings['pathogen_min_size'] != 0:
pathogen_mask = _filter_object(pathogen_mask, settings['pathogen_min_size'])
if settings['cytoplasm_min_size'] is not None and settings['cytoplasm_min_size'] != 0:
cytoplasm_mask = _filter_object(cytoplasm_mask, settings['cytoplasm_min_size'])
if settings.get('organelle_min_size') and settings['organelle_min_size'] != 0:
organelle_mask = _filter_object(organelle_mask, settings['organelle_min_size'])
if settings['cell_mask_dim'] is not None and settings['nucleus_mask_dim'] is not None and settings['pathogen_mask_dim'] is not None:
cell_mask, nucleus_mask, pathogen_mask, cytoplasm_mask = _exclude_objects(cell_mask, nucleus_mask, pathogen_mask, cytoplasm_mask, uninfected=settings['uninfected'])
data[:, :, settings['cell_mask_dim']] = cell_mask.astype(data_type)
if settings['nucleus_mask_dim'] is not None:
data[:, :, settings['nucleus_mask_dim']] = nucleus_mask.astype(data_type)
if settings['pathogen_mask_dim'] is not None:
data[:, :, settings['pathogen_mask_dim']] = pathogen_mask.astype(data_type)
if settings['cytoplasm']:
data = np.concatenate((data, cytoplasm_mask[:, :, np.newaxis]), axis=2)
if settings['plot']:
fig = _plot_cropped_arrays(data, file, figuresize)
figs[f'{file_name}__after_filtration'] = fig
if settings['save_measurements']:
cell_df, nucleus_df, pathogen_df, organelle_df, cytoplasm_df = _morphological_measurements(cell_mask, nucleus_mask, pathogen_mask, organelle_mask, cytoplasm_mask, settings)
cell_intensity_df, nucleus_intensity_df, pathogen_intensity_df, organelle_intensity_df, cytoplasm_intensity_df = _intensity_measurements(cell_mask, nucleus_mask, pathogen_mask, organelle_mask, cytoplasm_mask, channel_arrays, settings, sizes=[1, 2, 3, 4, 5], periphery=True, outside=True)
if settings['cell_mask_dim'] is not None:
_ = _merge_and_save_to_database(cell_df, cell_intensity_df, 'cell', source_folder, file_name, settings['experiment'], settings['timelapse'])
if settings['nucleus_mask_dim'] is not None:
_ = _merge_and_save_to_database(nucleus_df, nucleus_intensity_df, 'nucleus', source_folder, file_name, settings['experiment'], settings['timelapse'])
if settings['pathogen_mask_dim'] is not None:
_ = _merge_and_save_to_database(pathogen_df, pathogen_intensity_df, 'pathogen', source_folder, file_name, settings['experiment'], settings['timelapse'])
if settings.get('summarize_organelles_by') is not None:
if "organelle" in settings['summarize_organelles_by']:
if settings.get('organelle_mask_dim') is not None:
_ = _merge_and_save_to_database(organelle_df, organelle_intensity_df, 'organelle', source_folder, file_name, settings['experiment'], settings['timelapse'])
if settings['cytoplasm']:
cytoplasm_merged_df = _merge_and_save_to_database(cytoplasm_df, cytoplasm_intensity_df, 'cytoplasm', source_folder, file_name, settings['experiment'], settings['timelapse'])
if settings.get('summarize_organelles_by') is not None:
if "cell" in settings['summarize_organelles_by']:
if settings.get('organelle_mask_dim') is not None and np.max(organelle_mask) > 0:
if settings['cell_mask_dim'] is not None:
org_per_cell = _summarize_organelles_per_parent(organelle_mask, cell_mask, channel_arrays, parent_name='cell')
org_per_cell.columns = [f'organelle_summary_{col}' if col != 'label' else col for col in org_per_cell.columns]
_merge_and_save_to_database(org_per_cell, pd.DataFrame(), 'cell_organelle_summary', source_folder, file_name, settings['experiment'], settings['timelapse'])
if "nucleus" in settings['summarize_organelles_by']:
if settings['nucleus_mask_dim'] is not None:
org_per_nucleus = _summarize_organelles_per_parent(organelle_mask, nucleus_mask, channel_arrays, parent_name='nucleus')
org_per_nucleus.columns = [f'organelle_summary_{col}' if col != 'label' else col for col in org_per_nucleus.columns]
_merge_and_save_to_database(org_per_nucleus, pd.DataFrame(), 'nucleus_organelle_summary', source_folder, file_name, settings['experiment'], settings['timelapse'])
if "pathogen" in settings['summarize_organelles_by']:
if settings['pathogen_mask_dim'] is not None:
org_per_pathogen = _summarize_organelles_per_parent(organelle_mask, pathogen_mask, channel_arrays, parent_name='pathogen')
org_per_pathogen.columns = [f'organelle_summary_{col}' if col != 'label' else col for col in org_per_pathogen.columns]
_merge_and_save_to_database(org_per_pathogen, pd.DataFrame(), 'pathogen_organelle_summary', source_folder, file_name, settings['experiment'], settings['timelapse'])
if "cytoplasm" in settings['summarize_organelles_by']:
if settings['cytoplasm_mask_dim'] is not None:
org_per_cytoplasm = _summarize_organelles_per_parent(organelle_mask, cytoplasm_mask, channel_arrays, parent_name='cytoplasm')
org_per_cytoplasm.columns = [f'organelle_summary_{col}' if col != 'label' else col for col in org_per_cytoplasm.columns]
_merge_and_save_to_database(org_per_cytoplasm, pd.DataFrame(), 'cytoplasm_organelle_summary', source_folder, file_name, settings['experiment'], settings['timelapse'])
if settings['save_png'] or settings['save_arrays'] or settings['plot']:
if isinstance(settings['dialate_pngs'], bool):
dialate_pngs = [settings['dialate_pngs'], settings['dialate_pngs'], settings['dialate_pngs']]
if isinstance(settings['dialate_pngs'], list):
dialate_pngs = settings['dialate_pngs']
if isinstance(settings['dialate_png_ratios'], float):
dialate_png_ratios = [settings['dialate_png_ratios'], settings['dialate_png_ratios'], settings['dialate_png_ratios']]
if isinstance(settings['dialate_png_ratios'], list):
dialate_png_ratios = settings['dialate_png_ratios']
if isinstance(settings['crop_mode'], str):
crop_mode = [settings['crop_mode']]
if isinstance(settings['crop_mode'], list):
crop_ls = settings['crop_mode']
size_ls = settings['png_size']
if isinstance(size_ls[0], int):
size_ls = [size_ls]
if len(crop_ls) > 1 and len(size_ls) == 1:
size_ls = size_ls * len(crop_ls)
if len(crop_ls) != len(size_ls):
print(f"Setting: size_ls: {settings['png_size']} should be a list of integers, or a list of lists of integers if crop_ls: {settings['crop_mode']} has multiple elements")
for crop_idx, crop_mode in enumerate(crop_ls):
width, height = size_ls[crop_idx]
if crop_mode == 'cell':
crop_mask = cell_mask.copy()
dialate_png = dialate_pngs[crop_idx]
dialate_png_ratio = dialate_png_ratios[crop_idx]
elif crop_mode == 'nucleus':
crop_mask = nucleus_mask.copy()
dialate_png = dialate_pngs[crop_idx]
dialate_png_ratio = dialate_png_ratios[crop_idx]
elif crop_mode == 'pathogen':
crop_mask = pathogen_mask.copy()
dialate_png = dialate_pngs[crop_idx]
dialate_png_ratio = dialate_png_ratios[crop_idx]
elif crop_mode == 'organelle':
crop_mask = organelle_mask.copy()
dialate_png = dialate_pngs[crop_idx]
dialate_png_ratio = dialate_png_ratios[crop_idx]
elif crop_mode == 'cytoplasm':
crop_mask = cytoplasm_mask.copy()
dialate_png = False
else:
print(f'Value error: Posseble values for crop_mode are: cell, nucleus, pathogen, cytoplasm')
objects_in_image = np.unique(crop_mask)
objects_in_image = objects_in_image[objects_in_image != 0]
img_paths = []
for _id in objects_in_image:
region = (crop_mask == _id)
# Use the boolean mask to filter the cell_mask and then find unique IDs
region_cell_ids = np.atleast_1d(np.unique(cell_mask[region]))
region_nucleus_ids = np.atleast_1d(np.unique(nucleus_mask[region]))
region_pathogen_ids = np.atleast_1d(np.unique(pathogen_mask[region]))
if settings['use_bounding_box']:
region = _find_bounding_box(crop_mask, _id, buffer=10)
img_name, fldr, table_name = _generate_names(file_name=file_name, cell_id = region_cell_ids, cell_nucleus_ids=region_nucleus_ids, cell_pathogen_ids=region_pathogen_ids, source_folder=source_folder, crop_mode=crop_mode, timelapse=settings['timelapse'])
if dialate_png:
region_area = np.sum(region)
approximate_diameter = np.sqrt(region_area)
dialate_png_px = int(approximate_diameter * dialate_png_ratio)
struct = generate_binary_structure(2, 2)
region = binary_dilation(region, structure=struct, iterations=dialate_png_px)
if settings['save_png']:
fldr_type = f"{crop_mode}_png/"
png_folder = os.path.join(fldr,fldr_type)
img_path = os.path.join(png_folder, img_name)
img_paths.append(img_path)
png_channels = data[:, :, settings['png_dims']].astype(data_type)
if settings['normalize_by'] == 'fov':
if not settings['normalize'] is False:
percentile_list = _get_percentiles(png_channels, settings['normalize'][0], settings['normalize'][1])
png_channels = _crop_center(png_channels, region, new_width=width, new_height=height)
if isinstance(settings['normalize'], list):
if settings['normalize_by'] == 'png':
png_channels = normalize_to_dtype(png_channels, settings['normalize'][0], settings['normalize'][1])
if settings['normalize_by'] == 'fov':
png_channels = normalize_to_dtype(png_channels, settings['normalize'][0], settings['normalize'][1], percentile_list=percentile_list)
else:
png_channels = normalize_to_dtype(png_channels, 0, 100)
os.makedirs(png_folder, exist_ok=True)
if png_channels.shape[2] == 2:
dummy_channel = np.zeros_like(png_channels[:,:,0]) # Create a 2D zero array with same shape as one channel
png_channels = np.dstack((png_channels, dummy_channel))
grid = save_and_add_image_to_grid(png_channels, img_path, grid, settings['plot'])
else:
grid = save_and_add_image_to_grid(png_channels, img_path, grid, settings['plot'])
if len(img_paths) == len(objects_in_image):
filepaths_to_database(img_paths, settings, source_folder, crop_mode)
if settings['save_arrays']:
row_idx, col_idx = np.where(region)
region_array = data[row_idx.min():row_idx.max()+1, col_idx.min():col_idx.max()+1, :]
array_folder = f"{fldr}/region_array/"
os.makedirs(array_folder, exist_ok=True)
np.save(os.path.join(array_folder, img_name), region_array)
grid = save_and_add_image_to_grid(png_channels, img_path, grid, settings['plot'])
img_paths.append(img_path)
if len(img_paths) == len(objects_in_image):
filepaths_to_database(img_paths, settings, source_folder, crop_mode)
cells = np.unique(cell_mask)
except Exception as e:
print('main',e)
cells = 0
traceback.print_exc()
end = time.time()
duration = end-start
time_ls.append(duration)
average_time = np.mean(time_ls) if len(time_ls) > 0 else 0
if settings['plot']:
fig = img_list_to_grid(grid)
figs[f'{file_name}__pngs'] = fig
return index, average_time, cells, figs
#@log_function_call
[docs]
def measure_crop(settings):
"""
Measure the crop of an image based on the provided settings.
Args:
settings (dict): The settings for measuring the crop.
Returns:
None
"""
from .io import _save_settings_to_db
from .timelapse import _timelapse_masks_to_gif
from .utils import measure_test_mode, print_progress, delete_intermedeate_files, save_settings, format_path_for_system, normalize_src_path
from .settings import get_measure_crop_settings
if settings['timelapse']:
settings['save_png'] = False
if not isinstance(settings['src'], (str, list)):
ValueError(f'src must be a string or a list of strings')
return
settings['src'] = normalize_src_path(settings['src'])
if isinstance(settings['src'], str):
settings['src'] = [settings['src']]
if isinstance(settings['src'], list):
source_folders = settings['src']
for source_folder in source_folders:
print(f'Processing folder: {source_folder}')
source_folder = format_path_for_system(source_folder)
settings['src'] = source_folder
src = source_folder
settings = get_measure_crop_settings(settings)
settings = measure_test_mode(settings)
src_fldr = settings['src']
if not os.path.basename(src_fldr).endswith('merged'):
print(f"WARNING: Source folder, settings: src: {src_fldr} should end with '/merged'")
src_fldr = os.path.join(src_fldr, 'merged')
settings['src'] = src_fldr
print(f"Changed source folder to: {src_fldr}")
if settings['cell_mask_dim'] is None:
settings['uninfected'] = True
if settings['pathogen_mask_dim'] is None:
settings['uninfected'] = True
if settings['cell_mask_dim'] is not None and settings['pathogen_min_size'] is not None:
settings['cytoplasm'] = True
elif settings['cell_mask_dim'] is not None and settings['nucleus_min_size'] is not None:
settings['cytoplasm'] = True
else:
settings['cytoplasm'] = False
spacr_cores = int(mp.cpu_count() - 6)
if spacr_cores <= 2:
spacr_cores = 1
if settings['n_jobs'] > spacr_cores:
print(f'Warning reserving 6 CPU cores for other processes, setting n_jobs to {spacr_cores}')
settings['n_jobs'] = spacr_cores
settings_save = settings.copy()
settings_save['src'] = os.path.dirname(settings['src'])
save_settings(settings_save, name='measure_crop_settings', show=True)
if settings['timelapse_objects'] == 'nucleus':
if not settings['cell_mask_dim'] is None:
tlo = settings['timelapse_objects']
print(f'timelapse object:{tlo}, cells will be relabeled to nucleus labels to track cells.')
int_setting_keys = ['cell_mask_dim', 'nucleus_mask_dim', 'pathogen_mask_dim', 'cell_min_size', 'nucleus_min_size', 'pathogen_min_size', 'cytoplasm_min_size']
if isinstance(settings['normalize'], bool) and settings['normalize']:
print(f'WARNING: to notmalize single object pngs set normalize to a list of 2 integers, e.g. [1,99] (lower and upper percentiles)')
return
if isinstance(settings['normalize'], list) or isinstance(settings['normalize'], bool) and settings['normalize']:
if settings['normalize_by'] not in ['png', 'fov']:
print("Warning: normalize_by should be either 'png' to notmalize each png to its own percentiles or 'fov' to normalize each png to the fov percentiles ")
return
if not all(isinstance(settings[key], int) or settings[key] is None for key in int_setting_keys):
print(f"WARNING: {int_setting_keys} must all be integers")
return
if not isinstance(settings['channels'], list):
print(f"WARNING: channels should be a list of integers representing channels e.g. [0,1,2,3]")
return
if not isinstance(settings['crop_mode'], list):
print(f"WARNING: crop_mode should be a list with at least one element e.g. ['cell'] or ['cell','nucleus'] or [None] got: {settings['crop_mode']}")
settings['crop_mode'] = [settings['crop_mode']]
settings['crop_mode'] = [str(crop_mode) for crop_mode in settings['crop_mode']]
print(f"Converted crop_mode to list: {settings['crop_mode']}")
_save_settings_to_db(settings)
files = [f for f in os.listdir(settings['src']) if f.endswith('.npy')]
n_jobs = settings['n_jobs']
print(f'using {n_jobs} cpu cores')
print_progress(files_processed=0, files_to_process=len(files), n_jobs=n_jobs, time_ls=[], operation_type='Measure and Crop')
def job_callback(result):
completed_jobs.add(result[0])
process_meassure_crop_results([result], settings)
files_processed = len(completed_jobs)
files_to_process = len(files)
print_progress(files_processed, files_to_process, n_jobs, time_ls=time_ls, operation_type='Measure and Crop')
if files_processed >= files_to_process:
pool.terminate()
with mp.Manager() as manager:
time_ls = manager.list()
completed_jobs = set() # Set to keep track of completed jobs
with mp.Pool(n_jobs) as pool:
for index, file in enumerate(files):
pool.apply_async(_measure_crop_core, args=(index, time_ls, file, settings), callback=job_callback)
pool.close()
pool.join()
if settings['timelapse']:
if settings['timelapse_objects'] == 'nucleus':
folder_path = settings['src']
mask_channels = [settings['nucleus_mask_dim'], settings['pathogen_mask_dim'], settings['cell_mask_dim']]
object_types = ['nucleus', 'pathogen', 'cell']
_timelapse_masks_to_gif(folder_path, mask_channels, object_types)
if settings['delete_intermediate']:
delete_intermedeate_files(settings)
print("Successfully completed run")
[docs]
def process_meassure_crop_results(partial_results, settings):
"""
Process the results, display, and optionally save the figures.
Args:
partial_results (list): List of partial results.
settings (dict): Settings dictionary.
save_figures (bool): Flag to save figures or not.
"""
for result in partial_results:
if result is None:
continue
index, avg_time, cells, figs = result
if figs is not None:
for key, fig in figs.items():
part_1, part_2 = key.split('__')
save_dir = os.path.join(os.path.dirname(settings['src']), 'results', f"{part_1}")
os.makedirs(save_dir, exist_ok=True)
fig_path = os.path.join(save_dir, f"{part_2}.pdf")
fig.savefig(fig_path)
plt.figure(fig.number)
plt.show()
plt.close(fig)
result = (index, None, None, None)
[docs]
def generate_cellpose_train_set(folders, dst, min_objects=5):
os.makedirs(dst, exist_ok=True)
os.makedirs(os.path.join(dst,'masks'), exist_ok=True)
os.makedirs(os.path.join(dst,'imgs'), exist_ok=True)
for folder in folders:
mask_folder = os.path.join(folder, 'masks')
experiment_id = os.path.basename(folder)
for filename in os.listdir(mask_folder): # List the contents of the directory
path = os.path.join(mask_folder, filename)
img_path = os.path.join(folder, filename)
newname = experiment_id + '_' + filename
new_mask = os.path.join(dst, 'masks', newname)
new_img = os.path.join(dst, 'imgs', newname)
mask = cv2.imread(path, cv2.IMREAD_UNCHANGED)
if mask is None:
print(f"Error reading {path}, skipping.")
continue
nr_of_objects = len(np.unique(mask)) - 1 # Assuming 0 is background
if nr_of_objects >= min_objects: # Use >= to include min_objects
try:
shutil.copy(path, new_mask)
shutil.copy(img_path, new_img)
except Exception as e:
print(f"Error copying {path} to {new_mask}: {e}")
[docs]
def get_object_counts(src):
database_path = os.path.join(src, 'measurements/measurements.db')
# Connect to the SQLite database
conn = sqlite3.connect(database_path)
# Read the table into a pandas DataFrame
df = pd.read_sql_query("SELECT * FROM object_counts", conn)
# Group by 'count_type' and calculate the sum of 'object_count' and the average 'object_count' per 'file_name'
grouped_df = df.groupby('count_type').agg(
total_object_count=('object_count', 'sum'),
avg_object_count_per_file_name=('object_count', 'mean')
).reset_index()
# Close the database connection
conn.close()
return grouped_df