   1# -*- coding: utf-8 -*-
      2"""
      3Created on Wed Feb 22 13:46:15 2023
      4
      5@author: Bart Steemans. Govers Lab. 
      6"""
      7import os
      8import shutil
      9import tifffile
     10import logging
     11
     12import numpy as np
     13import cv2 as cv
     14import networkx as nx
     15import math
     16from shapely import Polygon
     17from .logging_utils import tqdm_if_verbose
     18import pandas as pd
     19from itertools import combinations
     20import matplotlib.pyplot as plt
     21
     22from scipy import interpolate, ndimage, spatial
     23from scipy.ndimage import gaussian_filter1d, fourier_shift, gaussian_laplace
     24from scipy.signal import find_peaks
     25from scipy.stats import kurtosis, skew
     26
     27from skimage.feature import graycomatrix, graycoprops
     28from skimage.filters import threshold_otsu
     29from skimage.draw import polygon2mask
     30from skimage.morphology import skeletonize, closing
     31from skimage.measure import regionprops_table, regionprops, label, shannon_entropy
     32from skimage.registration import phase_cross_correlation
     33from scipy.optimize import curve_fit
     34
     35bactoscoop_logger = logging.getLogger("logger")
     36
     37
     38# -------------------------------------------------------------------------------------------------------
     39# FUNCTIONS
     40# %% Cellular mesh creation functions
     41def smallest_path(x_coords, y_coords, start_x, start_y):
     42    G = nx.Graph()
     43    for i in range(len(x_coords)):
     44        G.add_node(i, x=x_coords[i], y=y_coords[i])
     45    for i in range(len(x_coords)):
     46        for j in range(i + 1, len(x_coords)):
     47            G.add_edge(
     48                i,
     49                j,
     50                weight=np.sqrt(
     51                    (x_coords[i] - x_coords[j]) ** 2 + (y_coords[i] - y_coords[j]) ** 2
     52                ),
     53            )
     54    start = None
     55    for i in range(len(x_coords)):
     56        if x_coords[i] == start_x and y_coords[i] == start_y:
     57            start = i
     58            break
     59    T = nx.minimum_spanning_tree(G)
     60    path = nx.dfs_postorder_nodes(T, source=start)
     61    sorted_points = list(path)
     62    return [x_coords[i] for i in sorted_points], [y_coords[i] for i in sorted_points]
     63
     64
     65def skeleton_endpoints(skel):
     66
     67    skel = skel.copy()
     68    skel[skel != 0] = 1
     69    skel = np.uint8(skel)
     70
     71    # Apply the convolution.
     72    kernel = np.uint8([[1, 1, 1], [1, 10, 1], [1, 1, 1]])
     73
     74    filtered = cv.filter2D(skel, -1, kernel)
     75
     76    out = np.zeros_like(skel)
     77    out[np.where(filtered == 11)] = 1
     78    return out
     79
     80
     81def get_ordered_list(points, x, y):
     82    distance = ((x - points[1]) ** 2 + (y - points[0]) ** 2) ** 0.5
     83    args = np.argsort(distance)
     84    return args
     85
     86
     87# @jit(nopython=True, cache=True)
     88# MOMIA
     89def line_intersect(x1, y1, x2, y2, x3, y3, x4, y4):  # Ax+By = C
     90    A1 = y2 - y1
     91    B1 = x1 - x2
     92    C1 = A1 * x1 + B1 * y1
     93    A2 = y4 - y3
     94    B2 = x3 - x4
     95    C2 = A2 * x3 + B2 * y3
     96    intersect_x = (B1 * C2 - B2 * C1) / (A2 * B1 - A1 * B2)
     97    intersect_y = (A1 * C2 - A2 * C1) / (B2 * A1 - B1 * A2)
     98    return intersect_x, intersect_y
     99
    100
    101# MOMIA
    102def line_contour_intersection(p1, p2, contour):
    103    v1, v2 = contour[:-1], contour[1:]
    104    x1, y1 = v1.T
    105    x2, y2 = v2.T
    106    x3, y3 = p1
    107    x4, y4 = p2
    108    xy = np.array(line_intersect(x1, y1, x2, y2, x3, y3, x4, y4)).T
    109    dxy_v1 = xy - v1
    110    dxy_v2 = xy - v2
    111    dxy = dxy_v1 * dxy_v2
    112    intersection_points = xy[np.where(np.logical_and(dxy[:, 0] < 0, dxy[:, 1] < 0))]
    113    if len(intersection_points) > 2:
    114        dist = np.sum(
    115            np.square(np.tile(p1, (len(intersection_points), 1)) - intersection_points),
    116            axis=1,
    117        )
    118        intersection_points = intersection_points[np.argsort(dist)[0:2]]
    119    return intersection_points
    120
    121
    122# MOMIA
    123def unit_perpendicular_vector(data, closed=True):
    124
    125    p1 = data[1:]
    126    p2 = data[:-1]
    127    dxy = p1 - p2
    128    ang = np.arctan2(dxy.T[1], dxy.T[0]) + 0.5 * np.pi
    129    dx, dy = np.cos(ang), np.sin(ang)
    130    unit_dxy = np.array([dx, dy]).T
    131    if not closed:
    132        unit_dxy = np.concatenate([[unit_dxy[0]], unit_dxy])
    133    else:
    134        unit_dxy = np.concatenate([unit_dxy, [unit_dxy[-1]]])
    135    return unit_dxy
    136
    137
    138def get_major_axis_line_coordinates(masks):
    139    # Find contours in the mask
    140    contours, _ = cv.findContours(
    141        masks.astype(np.uint8), cv.RETR_EXTERNAL, cv.CHAIN_APPROX_NONE
    142    )
    143
    144    ellipse = cv.fitEllipse(contours[0])
    145
    146    # Extract major axis parameters
    147    # Extract major axis parameters
    148    center, (major_axis, minor_axis), angle = ellipse
    149
    150    # Calculate endpoints of the major axis with reversed orientation
    151    angle_rad = np.deg2rad(angle)
    152    cos_angle = np.cos(angle_rad)
    153    sin_angle = np.sin(angle_rad)
    154
    155    # Center coordinates of the major axis
    156    mid_x = center[0]
    157    mid_y = center[1]
    158
    159    # Endpoint coordinates of the major axis with reversed orientation
    160    endpoint1_x = int(mid_x + 0.5 * major_axis * sin_angle)
    161    endpoint1_y = int(mid_y - 0.5 * major_axis * cos_angle)
    162    endpoint2_x = int(mid_x - 0.5 * major_axis * sin_angle)
    163    endpoint2_y = int(mid_y + 0.5 * major_axis * cos_angle)
    164
    165    # Generate coordinates along the major axis line
    166    line_x = np.linspace(endpoint1_x, endpoint2_x, num=200)
    167    line_y = np.linspace(endpoint1_y, endpoint2_y, num=200)
    168
    169    # Combine x and y coordinates into an array of shape (n, 2) with (x, y) convention
    170    line_coordinates = np.column_stack((line_y, line_x))
    171    return line_coordinates
    172
    173
    174def extract_skeleton(masks):
    175    skeleton = skeletonize(masks, method="lee")
    176
    177    # Retrieve the skeleton coordinates in the original image space
    178    skeleton_coords = np.column_stack(np.where(skeleton)).astype(int)
    179
    180    x_skel = skeleton_coords[:, 0]
    181    y_skel = skeleton_coords[:, 1]
    182
    183    skel_mask = np.zeros_like(masks)
    184    skel_mask[x_skel, y_skel] = 1
    185
    186    # Sort coordinates based on distance
    187    skel_end = skeleton_endpoints(skel_mask)
    188    try:
    189
    190        end_coords = np.argwhere(skel_end == 1)[0]
    191
    192    except IndexError:
    193
    194        return get_major_axis_line_coordinates(masks)
    195
    196    sorted_skelx, sorted_skely = smallest_path(
    197        x_skel, y_skel, end_coords[0], end_coords[1]
    198    )
    199
    200    if len(y_skel) <= 3:
    201        return get_major_axis_line_coordinates(masks)
    202
    203    elif 4 <= len(y_skel) <= 13:
    204
    205        fkc, u = interpolate.splprep([sorted_skelx, sorted_skely], k=1, s=1, per=0)
    206    else:
    207        fkc, u = interpolate.splprep([sorted_skelx, sorted_skely], k=2, s=10, per=0)
    208    smoothed_skel = np.asarray(
    209        interpolate.splev(np.linspace(u.min(), u.max(), 200), fkc, der=0)
    210    ).T
    211    return smoothed_skel
    212
    213
    214# MOMIA
    215def find_poles(smoothed_skeleton, smoothed_contour, find_pole1=True, find_pole2=True):
    216    # find endpoints and their nearest neighbors on a midline
    217    length = len(smoothed_skeleton)
    218    extended_pole1 = [smoothed_skeleton[0]]
    219    extended_pole2 = [smoothed_skeleton[-1]]
    220    i = 0
    221    j = 0
    222    if find_pole1:
    223        for i in range(10):
    224            p1 = smoothed_skeleton[i]
    225            p2 = smoothed_skeleton[i + 1]
    226            # find the two intersection points between
    227            # the vectorized contour and line through pole1
    228            intersection_points_pole1 = line_contour_intersection(
    229                p1, p2, smoothed_contour
    230            )
    231            # find the interesection point with the same direction as the outward pole vector
    232            dxy_1 = p1 - p2
    233            p1_tile = np.tile(p1, (len(intersection_points_pole1), 1))
    234            p1dot = (intersection_points_pole1 - p1_tile) * dxy_1
    235            index_1 = np.where((p1dot[:, 0] > 0) & (p1dot[:, 1] > 0))[0]
    236            if len(index_1) > 0:
    237                extended_pole1 = intersection_points_pole1[index_1]
    238                break
    239    else:
    240        i = 1
    241
    242    if find_pole2:
    243        for j in range(10):
    244            p3 = smoothed_skeleton[-1 - j]
    245            p4 = smoothed_skeleton[-2 - j]
    246            # find the two intersection points between
    247            # the vectorized contour and line through pole1
    248            intersection_points_pole2 = line_contour_intersection(
    249                p3, p4, smoothed_contour
    250            )
    251            # find the interesection point with the same direction as the outward pole vector
    252            dxy_2 = p3 - p4
    253            p3_tile = np.tile(p3, (len(intersection_points_pole2), 1))
    254            p3dot = (intersection_points_pole2 - p3_tile) * dxy_2
    255            index_2 = np.where((p3dot[:, 0] > 0) & (p3dot[:, 1] > 0))[0]
    256            if len(index_2) > 0:
    257                extended_pole2 = intersection_points_pole2[index_2]
    258                break
    259    else:
    260        j = 1
    261    trimmed_midline = smoothed_skeleton[i : length - j]
    262    return extended_pole1, extended_pole2, trimmed_midline
    263
    264
    265# MOMIA
    266def extend_skeleton(
    267    smoothed_skeleton, smoothed_contour, find_pole1=True, find_pole2=True
    268):
    269    # initiate approximated tip points
    270    new_pole1, new_pole2, trimmed_skeleton = find_poles(
    271        smoothed_skeleton,
    272        smoothed_contour,
    273        find_pole1=find_pole1,
    274        find_pole2=find_pole2,
    275    )
    276    extended_skeleton = np.concatenate([new_pole1, trimmed_skeleton, new_pole2])
    277
    278    extended_skeleton = spline_approximation(
    279        extended_skeleton, n=200, smooth_factor=1, closed=False
    280    )
    281    return extended_skeleton, new_pole1, new_pole2
    282
    283
    284# MOMIA
    285def straighten_by_orthogonal_lines(contour, midline, length, width, unit_micron=0.5):
    286    # estimate profile mesh size
    287    median_width = np.median(width)
    288    N_length = int(round(length / unit_micron))
    289    N_width = int(round(median_width / unit_micron))
    290    midline = spline_approximation(midline, n=N_length, smooth_factor=1, closed=False)
    291
    292    half_contour_1, half_contour_2 = divide_contour_by_midline(
    293        midline, contour
    294    )  # [:-1]
    295    # infer orthogonal vectors
    296    ortho_unit_vectors = unit_perpendicular_vector(midline)
    297    # generate orthogonal profile lines for all midline points except for the polar ones
    298    l1 = orthogonal_intersection_point(
    299        midline, half_contour_1, precomputed_orthogonal_vector=ortho_unit_vectors
    300    )
    301
    302    l2 = orthogonal_intersection_point(
    303        midline, half_contour_2, precomputed_orthogonal_vector=ortho_unit_vectors
    304    )
    305
    306    dl = (l2 - l1) / N_width
    307    mult_mat = np.tile(np.arange(N_width + 1), (len(l1), 1))
    308    mat_x = l1[:, 0][:, np.newaxis] + mult_mat * dl[:, 0][:, np.newaxis]
    309    mat_y = l1[:, 1][:, np.newaxis] + mult_mat * dl[:, 1][:, np.newaxis]
    310    profile_mesh = np.array([mat_x, mat_y])
    311    return l1, l2, profile_mesh, midline
    312
    313
    314def radial_intensities(
    315    cropped_signal_interp2d,
    316    contour,
    317    width,
    318    min_val,
    319    max_val,
    320    num_erosions,
    321    erosion_scale=None,
    322):
    323
    324    erosion = width / (2 * num_erosions)
    325    if erosion_scale is not None:
    326        erosion = erosion_scale
    327
    328    eroded_contours = [contour]
    329    for i in range(num_erosions):
    330        eroded_contour = erode_contour(eroded_contours[-1], scale=erosion)
    331        eroded_contours.append(eroded_contour)
    332
    333    eroded_contours.reverse()
    334
    335    intensities_per_mesh = []
    336    for contour in eroded_contours:
    337        intensities = measure_smoothened_intensity(
    338            contour, cropped_signal_interp2d, width=erosion, subpixel=0.1
    339        )
    340        normalized_intensities = (intensities - min_val) / (max_val - min_val)
    341        average_intensity = np.average(normalized_intensities)
    342        intensities_per_mesh.append(average_intensity)
    343
    344    return intensities_per_mesh
    345
    346
    347# MOMIA
    348def spline_approximation(init_contour, n=200, smooth_factor=1, closed=True):
    349    if closed:
    350        tck, u = interpolate.splprep(init_contour.T, u=None, s=smooth_factor, per=1)
    351    else:
    352        tck, u = interpolate.splprep(init_contour.T, u=None, s=smooth_factor)
    353    u_new = np.linspace(u.min(), u.max(), n)
    354    x_new, y_new = interpolate.splev(u_new, tck, der=0)
    355    return np.array([x_new, y_new]).T
    356
    357
    358# MOMIA
    359def divide_contour_by_midline(midline, contour):
    360    dist1 = distance_matrix(contour, midline[0]).flatten()
    361    dist2 = distance_matrix(contour, midline[-1]).flatten()
    362
    363    id1, id2 = np.argsort(dist1)[:2]
    364    id3, id4 = np.argsort(dist2)[:2]
    365
    366    contour_cp = contour.copy()
    367    if max(id1, id2) < max(id3, id4):
    368        term_p1 = max(id1, id2)
    369        if abs(id3 - id4) == 1:
    370            term_p2 = max(id3, id4) + 1
    371        elif abs(id3 - id4) > 1:
    372            term_p2 = max(id3, id4) + 2
    373        contour_cp = np.insert(contour_cp, term_p1, midline[0], axis=0)
    374        contour_cp = np.insert(contour_cp, term_p2, midline[-1], axis=0)
    375
    376    else:
    377        term_p1 = max(id3, id4)
    378        if abs(id1 - id2) == 1:
    379            term_p2 = max(id1, id2) + 1
    380        elif abs(id1 - id2) > 1:
    381            term_p2 = max(id1, id2) + 2
    382        contour_cp = np.insert(contour_cp, term_p1, midline[-1], axis=0)
    383        contour_cp = np.insert(contour_cp, term_p2, midline[0], axis=0)
    384
    385    if term_p1 == term_p2:
    386        raise ValueError("Two endpoints are identical!")
    387    else:
    388        pos1, pos2 = sorted([term_p1, term_p2])
    389        # print(id1, id2, id3, id4, pos1, pos2, len(contour_cp))
    390        half_contour_1 = contour_cp[pos1 : min(pos2 + 1, len(contour_cp) - 1)]
    391        half_contour_2 = np.concatenate([contour_cp[pos2:], contour_cp[: pos1 + 1]])
    392    return half_contour_1, half_contour_2
    393
    394
    395# MOMIA
    396def distance_matrix(data1, data2):
    397    x1, y1 = data1.T
    398    x2, y2 = data2.T
    399    dx = x1[:, np.newaxis] - x2
    400    dy = y1[:, np.newaxis] - y2
    401    dxy = np.sqrt(dx**2 + dy**2)
    402    return dxy
    403
    404
    405# MOMIA
    406def orthogonal_intersection_point(
    407    midline, outerline, precomputed_orthogonal_vector=None, min_dist=1e-100
    408):
    409    v1, v2 = outerline[:-1], outerline[1:]
    410    skel_x, skel_y = midline.T
    411    if precomputed_orthogonal_vector is None:
    412        intersect_x, intersect_y = intersect_matrix(midline, outerline)
    413    else:
    414        intersect_x, intersect_y = intersect_matrix(
    415            midline, outerline, orthogonal_vectors=precomputed_orthogonal_vector
    416        )
    417    dx_v1 = intersect_x - v1.T[0][:, np.newaxis]
    418    dx_v2 = intersect_x - v2.T[0][:, np.newaxis]
    419    dy_v1 = intersect_y - v1.T[1][:, np.newaxis]
    420    dy_v2 = intersect_y - v2.T[1][:, np.newaxis]
    421    dx = dx_v1 * dx_v2
    422    dy = dy_v1 * dy_v2
    423
    424    dist_x = skel_x[np.newaxis, :] - intersect_x
    425    dist_y = skel_y[np.newaxis, :] - intersect_y
    426    # influence on extreme points
    427    non_bounadry_points = np.where(np.logical_and(dy >= 0, dx >= 0))
    428    dist_matrix = np.sqrt(dist_x**2 + dist_y**2)
    429    dist_matrix[non_bounadry_points] = np.inf
    430    dist_matrix[dist_matrix <= min_dist] = np.inf
    431
    432    nearest_id_x = np.argsort(dist_matrix, axis=0)[:1]
    433    nearest_id_y = np.linspace(
    434        0, dist_matrix.shape[1] - 1, dist_matrix.shape[1]
    435    ).astype(int)
    436    pos_list = np.array(
    437        [
    438            intersect_x[nearest_id_x[0], nearest_id_y],
    439            intersect_y[nearest_id_x[0], nearest_id_y],
    440        ]
    441    ).T
    442    return pos_list
    443
    444
    445# MOMIA
    446def intersect_matrix(line, contour, orthogonal_vectors=None):
    447    if orthogonal_vectors is None:
    448        dxy = unit_perpendicular_vector(line, closed=False)
    449    else:
    450        dxy = orthogonal_vectors
    451    v1, v2 = contour[:-1], contour[1:]
    452    x1, y1 = v1.T
    453    x2, y2 = v2.T
    454    x3, y3 = line.T
    455    perp_xy = line + dxy
    456    x4, y4 = perp_xy.T
    457    A1 = y2 - y1
    458    B1 = x1 - x2
    459    C1 = A1 * x1 + B1 * y1
    460    A2 = y4 - y3
    461    B2 = x3 - x4
    462    C2 = A2 * x3 + B2 * y3
    463
    464    A1B2 = A1[:, np.newaxis] * B2
    465    A1C2 = A1[:, np.newaxis] * C2
    466    B1A2 = B1[:, np.newaxis] * A2
    467    B1C2 = B1[:, np.newaxis] * C2
    468    C1A2 = C1[:, np.newaxis] * A2
    469    C1B2 = C1[:, np.newaxis] * B2
    470    # can give outliers when B1A2 = A1B2, division by zero or close to zero
    471    intersect_x = (B1C2 - C1B2) / (B1A2 - A1B2)
    472    intersect_y = (A1C2 - C1A2) / (A1B2 - B1A2)
    473    return intersect_x, intersect_y
    474
    475
    476def add_poles(coords1, coords2, pole1, pole2):
    477    coords1 = np.vstack([pole1, coords1[1:-1], pole2])
    478    coords2 = np.vstack([pole1, coords2[1:-1], pole2])
    479    result = np.stack(
    480        [coords1[:, 0], coords1[:, 1], coords2[:, 0], coords2[:, 1]], axis=1
    481    )
    482    return result, coords1, coords2
    483
    484
    485# MOMIA
    486def line_length(line):
    487    v1 = line[:-1]
    488    v2 = line[1:]
    489    d = v2 - v1
    490    return np.sum(np.sqrt(np.sum(d**2, axis=1)))
    491
    492
    493# MOMIA
    494def direct_intersect_distance(skeleton, contour):
    495
    496    v1, v2 = contour[:-1], contour[1:]
    497    skel_x, skel_y = skeleton[1:-1].T
    498    intersect_x, intersect_y = intersect_matrix(skeleton[1:-1], contour)
    499    dx_v1 = intersect_x - v1.T[0][:, np.newaxis]
    500    dx_v2 = intersect_x - v2.T[0][:, np.newaxis]
    501    dy_v1 = intersect_y - v1.T[1][:, np.newaxis]
    502    dy_v2 = intersect_y - v2.T[1][:, np.newaxis]
    503    dx = dx_v1 * dx_v2
    504    dy = dy_v1 * dy_v2
    505    dist_x = skel_x[np.newaxis, :] - intersect_x
    506    dist_y = skel_y[np.newaxis, :] - intersect_y
    507
    508    non_boundry_points = np.where(np.logical_and(dy > 0, dx > 0))
    509    dist_matrix = np.sqrt(dist_x**2 + dist_y**2)
    510    dist_matrix[non_boundry_points] = np.inf
    511    nearest_id_x = np.argsort(dist_matrix, axis=0)[:2]
    512    nearest_id_y = np.linspace(
    513        0, dist_matrix.shape[1] - 1, dist_matrix.shape[1]
    514    ).astype(int)
    515    dists = (
    516        dist_matrix[nearest_id_x[0], nearest_id_y]
    517        + dist_matrix[nearest_id_x[1], nearest_id_y]
    518    )
    519    return np.concatenate([[0], dists, [0]])
    520
    521
    522def contour2mesh(contour, cropped_mask):
    523    try:
    524        skeleton = extract_skeleton(cropped_mask)
    525    except Exception:
    526        bactoscoop_logger.debug("Skeletonization failed")
    527        skeleton = []
    528    if np.any(skeleton):
    529        try:
    530            extended_skeleton, pole1, pole2 = extend_skeleton(
    531                skeleton, contour, find_pole1=True, find_pole2=True
    532            )
    533
    534            length = line_length(extended_skeleton)
    535            width = direct_intersect_distance(extended_skeleton, contour)
    536            if np.any(width):
    537                l1, l2, profile_mesh, midline = straighten_by_orthogonal_lines(
    538                    contour, extended_skeleton, length, width, unit_micron=0.5
    539                )
    540                result, l1, l2 = add_poles(l1, l2, pole1, pole2)
    541                return contour, result, midline
    542        except Exception:
    543            bactoscoop_logger.debug("width calculations for mesh creation failed")
    544            width = []
    545    return [], [], []
    546
    547
    548def draw_mask(contour, image):
    549    dilated_contour = expand_contour(contour, 0.5)
    550    mask_dilated = polygon2mask(image.shape, dilated_contour)
    551    return mask_dilated
    552
    553
    554def crop_image(image=None, contour=None, mask_to_crop=None, phase_img=None, pad=10):
    555    if contour is not None:
    556        cropped_contour = contour.copy()
    557        contour = contour.astype(int)
    558        x, y, w, h = crop(contour, pad, image)
    559        cropped_img = image[x : x + w, y : y + h]
    560
    561        if phase_img is not None:
    562            cropped_phase_img = phase_img[x : x + w, y : y + h]
    563            # Perform registration between cropped_img and cropped_phase_img
    564            shift, error, _diff = phase_cross_correlation(
    565                cropped_phase_img, cropped_img, upsample_factor=20
    566            )
    567            if max(np.abs(shift)) <= 15:
    568
    569                shift = np.array(shift)  # Ensure shift is in correct format
    570                cropped_img = shift_image(cropped_img, shift)
    571            else:
    572                cropped_img = image[x : x + w, y : y + h]
    573        cropped_mask = None
    574        cropped_contour -= [x, y]
    575
    576    elif mask_to_crop is not None:
    577        x, y, w, h = crop(mask_to_crop, pad, mask_to_crop, xy_coordinates=True)
    578        cropped_mask = mask_to_crop[y : y + h, x : x + w]
    579        cropped_img = None
    580        cropped_contour = contour
    581
    582    return cropped_img, cropped_mask, cropped_contour, x, y
    583
    584
    585def crop(to_be_cropped, pad, image, xy_coordinates=False):
    586
    587    padding = pad  # Default padding value
    588    x, y, w, h = cv.boundingRect(
    589        to_be_cropped
    590    )  # Get the smallest bounding rectangle of the contour
    591
    592    # Contours use internal (row, column) coordinates; mask bounding boxes
    593    # returned by OpenCV use (column, row) coordinates.
    594    x_limit, y_limit = (
    595        (image.shape[1], image.shape[0])
    596        if xy_coordinates
    597        else (image.shape[0], image.shape[1])
    598    )
    599    x = np.clip(x - padding, 0, x_limit)
    600    y = np.clip(y - padding, 0, y_limit)
    601    w = np.clip(w + 2 * padding, 0, x_limit - x)
    602    h = np.clip(h + 2 * padding, 0, y_limit - y)
    603    return x, y, w, h
    604
    605
    606def get_object_contours(mask, smoothing=0.1):  # x_offset, y_offset,
    607    contours, _ = cv.findContours(mask, cv.RETR_EXTERNAL, cv.CHAIN_APPROX_NONE)
    608    interpolated_contours = []
    609
    610    for contour in contours:
    611        contour = np.array(contour)
    612        contourx = contour[:, 0, 1]
    613        contoury = contour[:, 0, 0]
    614        s = int(len(contourx) * smoothing)
    615        try:
    616            tck, u = interpolate.splprep([contourx, contoury], u=None, s=s, per=1)
    617            u_new = np.linspace(u.min(), u.max(), len(contourx))
    618            outx, outy = interpolate.splev(u_new, tck, der=0)
    619            interpolated_contour = np.array([outx, outy]).T
    620        except Exception as e:
    621            bactoscoop_logger.debug(f"{e}")
    622            continue
    623        interpolated_contours.append(interpolated_contour)
    624
    625    return interpolated_contours
    626
    627def maybe_reverse_orientation(contour, result, midline):
    628    if np.random.rand() < 0.5:
    629
    630        contour = contour[::-1]
    631        result = result[::-1]
    632        midline = midline[::-1]
    633    return contour, result, midline
    634
    635
    636def _collect_touching_label_pairs(label_slice_a, label_slice_b):
    637    touching = (
    638        (label_slice_a != label_slice_b)
    639        & (label_slice_a != 0)
    640        & (label_slice_b != 0)
    641    )
    642    if not np.any(touching):
    643        return np.empty((0, 2), dtype=np.int64)
    644
    645    pairs = np.stack((label_slice_a[touching], label_slice_b[touching]), axis=1)
    646    pairs = np.sort(pairs, axis=1)
    647    return np.unique(pairs, axis=0)
    648
    649
    650def get_label_neighbor_counts(labels, connectivity=1):
    651    labels = np.asarray(labels)
    652    if labels.ndim != 2:
    653        raise ValueError("labels must be a 2D labeled mask.")
    654
    655    if connectivity not in (1, 4, 2, 8):
    656        raise ValueError("connectivity must be one of 1, 4, 2, or 8.")
    657
    658    unique_labels = np.unique(labels)
    659    unique_labels = unique_labels[unique_labels != 0]
    660    neighbor_map = {int(label): set() for label in unique_labels}
    661
    662    pair_arrays = [
    663        _collect_touching_label_pairs(labels[:, :-1], labels[:, 1:]),
    664        _collect_touching_label_pairs(labels[:-1, :], labels[1:, :]),
    665    ]
    666
    667    if connectivity in (2, 8):
    668        pair_arrays.extend(
    669            [
    670                _collect_touching_label_pairs(labels[:-1, :-1], labels[1:, 1:]),
    671                _collect_touching_label_pairs(labels[:-1, 1:], labels[1:, :-1]),
    672            ]
    673        )
    674
    675    non_empty_pairs = [pairs for pairs in pair_arrays if len(pairs) > 0]
    676    if non_empty_pairs:
    677        all_pairs = np.unique(np.vstack(non_empty_pairs), axis=0)
    678        for left_label, right_label in all_pairs:
    679            left_label = int(left_label)
    680            right_label = int(right_label)
    681            neighbor_map[left_label].add(right_label)
    682            neighbor_map[right_label].add(left_label)
    683
    684    return {label: len(neighbors) for label, neighbors in neighbor_map.items()}
    685
    686
    687def filter_labels_by_neighbor_count(labels, max_neighbors, connectivity=1):
    688    if max_neighbors is None:
    689        raise ValueError("max_neighbors must be provided.")
    690
    691    labels = np.asarray(labels)
    692    if labels.ndim != 2:
    693        raise ValueError("labels must be a 2D labeled mask.")
    694
    695    neighbor_counts = get_label_neighbor_counts(labels, connectivity=connectivity)
    696    removed_labels = sorted(
    697        label for label, count in neighbor_counts.items() if count > max_neighbors
    698    )
    699
    700    filtered_labels = labels.copy()
    701    if removed_labels:
    702        filtered_labels[np.isin(filtered_labels, removed_labels)] = 0
    703
    704    remaining_labels = np.unique(filtered_labels)
    705    remaining_labels = remaining_labels[remaining_labels != 0]
    706
    707    if remaining_labels.size == 0:
    708        relabeled = np.zeros_like(filtered_labels, dtype=np.int32)
    709    else:
    710        lookup = np.zeros(int(remaining_labels.max()) + 1, dtype=np.int32)
    711        lookup[remaining_labels] = np.arange(
    712            1, remaining_labels.size + 1, dtype=np.int32
    713        )
    714        relabeled = lookup[filtered_labels]
    715
    716    stats = {
    717        "connectivity": 4 if connectivity in (1, 4) else 8,
    718        "max_neighbors": int(max_neighbors),
    719        "total_label_count": int(len(neighbor_counts)),
    720        "removed_label_count": int(len(removed_labels)),
    721        "kept_label_count": int(len(remaining_labels)),
    722        "max_neighbor_count": int(max(neighbor_counts.values(), default=0)),
    723        "mean_neighbor_count": float(
    724            np.mean(list(neighbor_counts.values())) if neighbor_counts else 0.0
    725        ),
    726        "removed_labels": removed_labels,
    727    }
    728    return relabeled, stats
    729
    730def get_cellular_mesh(masks, smoothing):
    731    meshdata = {}
    732    cell_labels = np.unique(masks)
    733    cell_labels = cell_labels[cell_labels != 0]
    734
    735    for cell in tqdm_if_verbose(cell_labels):
    736        temp_mask = np.where(masks == cell, 1, 0)
    737        temp_mask = cv.convertScaleAbs(temp_mask)
    738
    739        _, cropped_mask, _, x, y = crop_image(mask_to_crop=temp_mask)
    740
    741        cell_contours = get_object_contours(cropped_mask, smoothing=smoothing)
    742        for cn in cell_contours:
    743            contour, result, midline = contour2mesh(cn, cropped_mask)
    744            contour, result, midline = maybe_reverse_orientation(contour, result, midline)
    745            if np.any(contour):
    746
    747                contour += [y, x]
    748                result += [y, x, y, x]
    749                midline += [y, x]
    750
    751                meshdata[cell] = {
    752                    "contour": contour,
    753                    "mesh": result,
    754                    "midline": midline,
    755                }
    756    meshdataframe = pd.DataFrame.from_dict(
    757        meshdata, orient="index", columns=["contour", "mesh", "midline"]
    758    )
    759    meshdataframe.reset_index(drop=True, inplace=True)
    760    return meshdataframe
    761
    762
    763# MOMIA
    764def distance(v1, v2):
    765    # Euclidean distance of two points
    766    return np.sqrt(np.sum((np.array(v1) - np.array(v2)) ** 2))
    767
    768
    769# MOMIA
    770def measure_width(extended_skeleton, smoothed_contour):
    771    length = line_length(extended_skeleton)
    772    d_perp = unit_perpendicular_vector(extended_skeleton, closed=False)
    773    width_list = []
    774    for i in range(1, len(extended_skeleton) - 1):
    775        xy = line_contour_intersection(
    776            extended_skeleton[i], d_perp[i] + extended_skeleton[i], smoothed_contour
    777        )
    778        coords = np.average(xy, axis=0)
    779        if (len(xy) == 2) and (np.isnan(coords).sum() == 0):
    780            width_list.append(distance(xy[0], xy[1]))
    781        else:
    782            raise ValueError(
    783                "Error encountered while computing line intersection points!"
    784            )
    785    return np.array([0] + width_list + [0]), length
    786
    787
    788def get_cell_pairs(pole_1, pole_2, cell_ID, maxdist=None):
    789    dm11 = distance_matrix(pole_1, pole_1)
    790    dm22 = distance_matrix(pole_2, pole_2)
    791    dm12 = distance_matrix(pole_1, pole_2)
    792    for i in cell_ID:
    793        dm11 = np.insert(
    794            dm11, i - 1, values=np.zeros((1, dm11.shape[1]), dtype=int), axis=0
    795        )
    796        dm11 = np.insert(dm11, i - 1, values=np.zeros(dm11.shape[0], dtype=int), axis=1)
    797        dm22 = np.insert(
    798            dm22, i - 1, values=np.zeros((1, dm22.shape[1]), dtype=int), axis=0
    799        )
    800        dm22 = np.insert(dm22, i - 1, values=np.zeros(dm22.shape[0], dtype=int), axis=1)
    801        dm12 = np.insert(
    802            dm12, i - 1, values=np.zeros((1, dm12.shape[1]), dtype=int), axis=0
    803        )
    804        dm12 = np.insert(dm12, i - 1, values=np.zeros(dm12.shape[0], dtype=int), axis=1)
    805    if maxdist != None:
    806        maxdistance = maxdist
    807        neighboring_cells_dm11 = np.argwhere(
    808            np.logical_and(dm11 > 0, dm11 <= maxdistance)
    809        )
    810        neighboring_cells_dm22 = np.argwhere(
    811            np.logical_and(dm22 > 0, dm22 <= maxdistance)
    812        )
    813        neighboring_cells_dm12 = np.argwhere(
    814            np.logical_and(dm12 > 0, dm12 <= maxdistance)
    815        )
    816        nc_concat = np.concatenate(
    817            (neighboring_cells_dm11, neighboring_cells_dm22, neighboring_cells_dm12),
    818            axis=0,
    819        )
    820        unique_nc1 = np.unique(nc_concat, axis=0)
    821        unique_nc2 = unique_nc1[unique_nc1[:, 0] != unique_nc1[:, 1]]
    822    else:
    823        # get index list of cell couples. Loop decides the distance for when more than 2 cells from couples with other cells
    824        maxdistance = 1
    825        while True:  # or maxdist = 8
    826            neighboring_cells_dm11 = np.argwhere(
    827                np.logical_and(dm11 > 0, dm11 <= maxdistance)
    828            )
    829            neighboring_cells_dm22 = np.argwhere(
    830                np.logical_and(dm22 > 0, dm22 <= maxdistance)
    831            )
    832            neighboring_cells_dm12 = np.argwhere(
    833                np.logical_and(dm12 > 0, dm12 <= maxdistance)
    834            )
    835            nc_concat = np.concatenate(
    836                (
    837                    neighboring_cells_dm11,
    838                    neighboring_cells_dm22,
    839                    neighboring_cells_dm12,
    840                ),
    841                axis=0,
    842            )
    843            unique_nc1 = np.unique(nc_concat, axis=0)
    844            unique_nc2 = unique_nc1[unique_nc1[:, 0] != unique_nc1[:, 1]]
    845            unique_values, inverse_indices = np.unique(unique_nc2, return_inverse=True)
    846            counts = np.bincount(inverse_indices)
    847            # Check if the current neighboring cells are already in the set of unique cells
    848            if counts == []:
    849                print("No neighboring cells present in the image")
    850                break
    851            try:
    852                if np.max(counts) > 2:
    853                    print(
    854                        f"Found triple neighboring cells after {maxdistance} distance units"
    855                    )
    856                    maxdistance = maxdistance
    857                    break
    858            except ValueError:
    859                pass
    860
    861            # Increment the maximum distance to search for
    862            maxdistance += 1
    863            if maxdistance > 10:
    864                print("pole distance to reach 4 stitched cells exceeds 20")
    865                break
    866            if counts == []:
    867                unique_nc2 = []
    868            else:
    869                neighboring_cells_dm11 = np.argwhere(
    870                    np.logical_and(dm11 > 0, dm11 <= maxdistance)
    871                )
    872                neighboring_cells_dm22 = np.argwhere(
    873                    np.logical_and(dm22 > 0, dm22 <= maxdistance)
    874                )
    875                neighboring_cells_dm12 = np.argwhere(
    876                    np.logical_and(dm12 > 0, dm12 <= maxdistance)
    877                )
    878                nc_concat = np.concatenate(
    879                    (
    880                        neighboring_cells_dm11,
    881                        neighboring_cells_dm22,
    882                        neighboring_cells_dm12,
    883                    ),
    884                    axis=0,
    885                )
    886                unique_nc1 = np.unique(nc_concat, axis=0)
    887                unique_nc2 = unique_nc1[unique_nc1[:, 0] != unique_nc1[:, 1]]
    888    return unique_nc2
    889
    890
    891def mesh2contour(x1, y1, x2, y2):
    892    x2f = np.flip(x2)
    893    y2f = np.flip(y2)
    894    # Concatenate the x and y coordinates
    895    xspline = np.concatenate((x2f[1:], x1[1:]))
    896    yspline = np.concatenate((y2f[1:], y1[1:]))
    897
    898    tck, u = interpolate.splprep(np.array([xspline, yspline]), k=3, s=2, per=1)
    899    u_new = np.linspace(u.min(), u.max(), 200)
    900    outx, outy = interpolate.splev(u_new, tck)
    901
    902    return np.array([outx, outy]).T
    903
    904
    905def split_point(x1, y1, x2, y2, ctpos):
    906    dx = x2 - x1
    907    dy = y2 - y1
    908    dist = np.sqrt(dx**2 + dy**2)
    909    # ctpos is the index at which the cell is most constricted and where the algorith will split the cell
    910    xc1 = x1[ctpos]
    911    yc1 = y1[ctpos]
    912    xc2 = x2[ctpos]
    913    yc2 = y2[ctpos]
    914    newpole_x = (xc1 + xc2) / 2
    915    newpole_y = (yc1 + yc2) / 2
    916
    917    xt1, yt1 = x1[ctpos + 1], y1[ctpos + 1]
    918    distt = dist[ctpos + 1]
    919    # Use indexing to get the x and y coordinates separately
    920    xt13, yt13 = (
    921        np.array([xt1, yt1])
    922        + (distt / 3) * np.array([dx[ctpos + 1], dy[ctpos + 1]]) / distt
    923    )
    924    xt23, yt23 = (
    925        np.array([xt1, yt1])
    926        + (distt * (2 / 3)) * np.array([dx[ctpos + 1], dy[ctpos + 1]]) / distt
    927    )
    928
    929    xs1, ys1 = x1[ctpos - 1], y1[ctpos - 1]
    930    dists = dist[ctpos - 1]
    931
    932    # Use indexing to get the x and y coordinates separately
    933    xs13, ys13 = (
    934        np.array([xs1, ys1])
    935        + (dists / 3) * np.array([dx[ctpos - 1], dy[ctpos - 1]]) / dists
    936    )
    937    xs23, ys23 = (
    938        np.array([xs1, ys1])
    939        + (dists * (2 / 3)) * np.array([dx[ctpos - 1], dy[ctpos - 1]]) / dists
    940    )
    941
    942    return (
    943        newpole_x,
    944        newpole_y,
    945        np.array([[xt13, yt13, xt23, yt23]]),
    946        np.array([[xs13, ys13, xs23, ys23]]),
    947    )
    948
    949
    950# calculates the constriction degree in relative and absolute value,
    951# and the relative position along the length of the cell
    952# absolute constr degree needs to be multiplied by the width
    953def constr_degree_single_cell_min(
    954    intensity,
    955    new_step_length,
    956    upper_limit=0.75,
    957    lower_limit=0.25,
    958    interp_factor=10,
    959):
    960    """
    961    Compute constriction degree and return both coarse and fine minima positions.
    962
    963    Parameters
    964    ----------
    965    intensity : array-like
    966        1D intensity profile.
    967    new_step_length : array-like
    968        Currently unused here, kept for API compatibility.
    969    upper_limit : float, optional
    970        Upper bound of region to search, as fraction of profile length.
    971    lower_limit : float, optional
    972        Lower bound of region to search, as fraction of profile length.
    973    interp_factor : int, optional
    974        Upsampling factor for fine-scale interpolation.
    975
    976    Returns
    977    -------
    978    minsize : float
    979        Relative constriction degree.
    980    ctposr : float
    981        Relative position of the constriction in [0, 1].
    982    minsizeabs : float
    983        Absolute constriction depth.
    984    ctpos_coarse : int or list
    985        Constriction position in original intensity index space.
    986    ctpos_fine : int or list
    987        Constriction position in interpolated fine index space.
    988    """
    989    minsize = 0
    990    minsizeabs = 0
    991    ctpos_coarse = []
    992    ctpos_fine = []
    993    ctposr = np.nan
    994
    995    intensity = np.asarray(intensity).squeeze()
    996
    997    if intensity.ndim != 1 or len(intensity) < 3:
    998        return minsize, ctposr, minsizeabs, ctpos_coarse, ctpos_fine
    999
   1000    # Original and fine coordinate grids
   1001    x = np.linspace(0, 1, len(intensity))
   1002    x_fine = np.linspace(0, 1, len(intensity) * interp_factor)
   1003
   1004    try:
   1005        f_interp = interpolate.interp1d(
   1006            x, intensity, kind="cubic", bounds_error=False, fill_value="extrapolate"
   1007        )
   1008        intensity_fine = f_interp(x_fine)
   1009    except Exception:
   1010        return minsize, ctposr, minsizeabs, ctpos_coarse, ctpos_fine
   1011
   1012    # Search window in fine space
   1013    start_idx = int(lower_limit * len(x_fine))
   1014    end_idx = int(upper_limit * len(x_fine))
   1015
   1016    # Keep valid bounds
   1017    start_idx = max(0, min(start_idx, len(x_fine) - 1))
   1018    end_idx = max(start_idx + 1, min(end_idx, len(x_fine)))
   1019
   1020    intensity_section = intensity_fine[start_idx:end_idx]
   1021
   1022    try:
   1023        minima = np.concatenate(
   1024            (
   1025                [False],
   1026                (intensity_section[1:-1] < intensity_section[:-2])
   1027                & (intensity_section[1:-1] < intensity_section[2:]),
   1028                [False],
   1029            )
   1030        )
   1031    except Exception:
   1032        return minsize, ctposr, minsizeabs, ctpos_coarse, ctpos_fine
   1033
   1034    if (not np.any(minima)) or np.sum(intensity_section) == 0:
   1035        return minsize, ctposr, minsizeabs, ctpos_coarse, ctpos_fine
   1036
   1037    index, dh, dhi, hgt = identify_minima(intensity_section, minima)
   1038    fix = np.argmax(dh)
   1039
   1040    minsizeabs = dhi[fix]
   1041    hgt_fix = hgt[fix]
   1042
   1043    if hgt_fix != 0:
   1044        minsize = minsizeabs / hgt_fix
   1045
   1046    # Fine-space global index
   1047    ctpos_fine = int(start_idx + index[fix])
   1048
   1049    # Relative position [0, 1]
   1050    ctposr = float(x_fine[ctpos_fine])
   1051
   1052    # Map fine position back to original/coarse intensity index
   1053    ctpos_coarse = int(round(ctposr * (len(intensity) - 1)))
   1054    ctpos_coarse = max(0, min(ctpos_coarse, len(intensity) - 1))
   1055
   1056    return minsize, ctposr, minsizeabs, ctpos_coarse, ctpos_fine
   1057
   1058
   1059def interp2d(image):
   1060    ff = interpolate.RectBivariateSpline(
   1061        range(image.shape[0]), range(image.shape[1]), image, kx=1, ky=1
   1062    )
   1063    return ff
   1064
   1065
   1066def get_image_for_frame(df_nr, images):
   1067    return images[df_nr]
   1068
   1069
   1070def get_cell_mask(contour, midline, image):
   1071    cropped_img, cropped_mask, cropped_contour, x, y = crop_image(
   1072        image=image, contour=contour, mask_to_crop=None
   1073    )
   1074    cropped_midline = midline.copy()
   1075    cropped_midline -= [x, y]
   1076    mask = draw_mask(cropped_contour, cropped_img)
   1077    return mask, cropped_img, cropped_contour, cropped_midline, x, y
   1078
   1079
   1080# %% Object detection and mesh creation functions
   1081def keep_significant_masks(
   1082    nucleoid_mask, cropped_mask, min_overlap_ratio=0.01, max_external_ratio=0.1
   1083):
   1084    # Find connected components (blobs) in the nucleoid_mask
   1085    num_labels, labels, stats, _ = cv.connectedComponentsWithStats(
   1086        nucleoid_mask.astype(np.uint8)
   1087    )
   1088
   1089    # Create a mask to keep blobs with significant overlap
   1090    significant_mask = np.zeros_like(nucleoid_mask, dtype=np.uint8)
   1091
   1092    # Get the total area of the cropped_mask
   1093    total_cropped_mask_area = np.sum(cropped_mask)
   1094
   1095    for label in range(1, num_labels):  # Skip label 0 as it represents the background
   1096        # Extract the blob region for the current label
   1097        blob_mask = (labels == label).astype(np.uint8)
   1098
   1099        # Calculate the overlap with the cropped_mask
   1100        overlap_area = np.sum(np.logical_and(blob_mask, cropped_mask))
   1101
   1102        # Calculate the overlap ratio
   1103        overlap_ratio = overlap_area / total_cropped_mask_area
   1104
   1105        # Calculate the external ratio (portion outside of cropped_mask)
   1106        external_area = np.sum(np.logical_and(blob_mask, ~cropped_mask))
   1107        external_ratio = external_area / np.sum(blob_mask)
   1108        # Keep blobs that have a significant overlap with the cropped_mask
   1109        # and don't exceed the maximum external ratio
   1110        if overlap_ratio >= min_overlap_ratio and external_ratio <= max_external_ratio:
   1111            significant_mask |= blob_mask
   1112
   1113    return significant_mask
   1114
   1115
   1116def get_object_masks(contour, object_contours, object_midlines, object_meshes, signal):
   1117
   1118    masks = []
   1119    cropped_img, cropped_mask, cropped_contour, x, y = crop_image(
   1120        image=signal, contour=contour, mask_to_crop=None
   1121    )
   1122    cropped_object_contours, cropped_object_midlines, cropped_object_meshes = (
   1123        get_cropped_object_data(object_contours, object_midlines, object_meshes, x, y)
   1124    )
   1125    for cropped_object_contour in cropped_object_contours:
   1126
   1127        temp_masks = draw_mask(cropped_object_contour, cropped_img)
   1128        masks.append(temp_masks)
   1129
   1130    return (
   1131        masks,
   1132        cropped_img,
   1133        cropped_object_contours,
   1134        cropped_object_midlines,
   1135        cropped_object_meshes,
   1136        x,
   1137        y,
   1138    )
   1139
   1140
   1141def get_subcellular_objects(
   1142    contour,
   1143    signal,
   1144    image,
   1145    cell_id,
   1146    smoothing,
   1147    log_sigma,
   1148    kernel_width,
   1149    min_overlap_ratio=0.01,
   1150    max_external_ratio=0.1,
   1151):
   1152
   1153    cropped_img, _, cropped_contour, x_offset, y_offset = crop_image(
   1154        image=signal, contour=contour, mask_to_crop=None, phase_img=image
   1155    )
   1156
   1157
   1158    cropped_mask = draw_mask(cropped_contour, cropped_img)
   1159    bgr = median_background(cropped_img, cropped_mask)
   1160    cropped_img_bgr = np.maximum(cropped_img - bgr, 0)
   1161
   1162    # Apply Laplacian of Gaussian (LoG) filter to the cropped_img
   1163    log_filtered_img = gaussian_laplace(cropped_img_bgr, sigma=log_sigma)
   1164
   1165    threshold_value = threshold_otsu(log_filtered_img)
   1166
   1167    object_mask = log_filtered_img < threshold_value
   1168
   1169    kernel = np.ones(
   1170        (kernel_width, kernel_width), np.uint8
   1171    )  # You can adjust the kernel size as needed
   1172    object_mask = cv.dilate(object_mask.astype(np.uint8), kernel, iterations=1)
   1173
   1174    object_mask = keep_significant_masks(
   1175        object_mask, cropped_mask, min_overlap_ratio, max_external_ratio
   1176    )
   1177
   1178    if np.any(object_mask == 1):
   1179        object_contours = get_object_contours(
   1180            object_mask, smoothing
   1181        )  # , x_offset, y_offset
   1182    else:
   1183        object_contours = []
   1184
   1185    return object_contours, cropped_img, object_mask, x_offset, y_offset
   1186
   1187
   1188def get_object_mesh(
   1189    contour,
   1190    signal,
   1191    image,
   1192    cell_id,
   1193    smoothing,
   1194    log_sigma,
   1195    kernel_width,
   1196    min_overlap_ratio,
   1197    max_external_ratio,
   1198):
   1199    meshdata = {
   1200        "cell_id": cell_id,
   1201        "object_contour": [],
   1202        "cropped_object_contour": [],
   1203        "object_mesh": [],
   1204        "object_midline": [],
   1205        "mesh_construction_errors": [],
   1206    }
   1207    object_contours, cropped_img, object_mask, x, y = get_subcellular_objects(
   1208        contour,
   1209        signal,
   1210        image,
   1211        cell_id,
   1212        smoothing,
   1213        log_sigma,
   1214        kernel_width,
   1215        min_overlap_ratio,
   1216        max_external_ratio,
   1217    )
   1218
   1219    if object_contours:
   1220
   1221        for object_index, cn in enumerate(object_contours):
   1222            meshdata["object_contour"].append(cn + [x, y])
   1223            meshdata["cropped_object_contour"].append(cn)
   1224
   1225            try:
   1226                temp_mask = draw_mask(cn, object_mask)
   1227                _, result, midline = contour2mesh(cn, temp_mask)
   1228                if np.any(result) and result.shape[0] < 800:
   1229                    result += [x, y, x, y]
   1230                    midline += [x, y]
   1231                    meshdata["object_mesh"].append(result)
   1232                    meshdata["object_midline"].append(midline)
   1233                    continue
   1234            except Exception as exc:
   1235                meshdata["mesh_construction_errors"].append(
   1236                    {
   1237                        "object_index": object_index,
   1238                        "error_type": type(exc).__name__,
   1239                        "error": str(exc),
   1240                    }
   1241                )
   1242
   1243            # Keep the contour and a matching placeholder for an invalid mesh.
   1244            meshdata["object_mesh"].append(np.array([]))
   1245            meshdata["object_midline"].append(np.array([]))
   1246
   1247    return meshdata
   1248
   1249
   1250# %% File readers/handlers
   1251from natsort import natsorted
   1252import re
   1253
   1254
   1255def image_field_key(filename, suffix="", is_mask=False):
   1256    """Return a normalized acquisition/field key from a TIFF filename.
   1257
   1258    Channel suffixes such as ``C1``/``C2`` and mask suffixes such as
   1259    ``cp_masks`` are removed so corresponding files can be matched by field.
   1260    """
   1261    stem = os.path.splitext(os.path.basename(os.fspath(filename)))[0]
   1262    if is_mask:
   1263        stem = re.sub(r"(?:[_-]?cp)?[_-]?masks?$", "", stem, flags=re.IGNORECASE)
   1264        # A mask tagged for C2 must not be paired with a C1 phase image.
   1265        if suffix:
   1266            suffix = str(suffix)
   1267            if stem.casefold().endswith(suffix.casefold()):
   1268                stem = stem[:-len(suffix)].rstrip(" _-")
   1269        else:
   1270            stem = re.sub(r"[_-]?C\d+$", "", stem, flags=re.IGNORECASE)
   1271    elif suffix:
   1272        suffix = str(suffix)
   1273        if stem.casefold().endswith(suffix.casefold()):
   1274            stem = stem[: -len(suffix)].rstrip(" _-")
   1275    else:
   1276        stem = re.sub(r"[_-]?C\d+$", "", stem, flags=re.IGNORECASE)
   1277    return stem.strip().casefold()
   1278
   1279
   1280def read_channels(folder_path, images_dict, channel_list=None):
   1281    if channel_list is None:
   1282        # If no suffixes are provided, assume an empty list
   1283        channel_list = []
   1284
   1285    missing_channels = []
   1286    for channel in channel_list:
   1287        matching_files = natsorted([
   1288            f
   1289            for f in os.listdir(folder_path)
   1290            if f.casefold().endswith(f"{channel}.tif".casefold())
   1291            or f.casefold().endswith(f"{channel}.tiff".casefold())
   1292        ])
   1293
   1294        if matching_files:
   1295            images_dict[channel] = [
   1296                tifffile.imread(os.path.join(folder_path, tiff_file))
   1297                for tiff_file in matching_files
   1298            ]
   1299            bactoscoop_logger.info(
   1300                f"Images with specified suffixes {channel} found in folder: {folder_path}"
   1301            )
   1302        else:
   1303            missing_channels.append(channel)
   1304    if missing_channels:
   1305        raise ValueError(
   1306            f"No images found for requested channel(s) {missing_channels} in folder: {folder_path}"
   1307        )
   1308
   1309    return images_dict
   1310
   1311
   1312def read_tiff_folder(folder_path, suffix="", include_paths=False):
   1313    # Get a list of all TIFF files in the folder
   1314    tiff_files = natsorted([
   1315        f
   1316        for f in os.listdir(folder_path)
   1317        if f.casefold().endswith(f"{suffix}.tif".casefold())
   1318        or f.casefold().endswith(f"{suffix}.tiff".casefold())
   1319    ])
   1320    if not tiff_files:
   1321        # Raise an exception if no TIFF files are found in the folder
   1322        raise ValueError("\nNo TIFF files found in folder: " + folder_path)
   1323
   1324    images = []
   1325    file_names = []
   1326    paths = [] if include_paths else None
   1327
   1328    for tiff_file in tiff_files:
   1329        tiff_path = os.path.join(folder_path, tiff_file)
   1330        image = tifffile.imread(tiff_path)
   1331        images.append(image)
   1332        file_names.append(tiff_file)
   1333        if include_paths:
   1334            paths.append(tiff_path)
   1335
   1336    if len(images) == 1:
   1337        bactoscoop_logger.info(f"Only one TIFF file found in folder: {folder_path}")
   1338        images = np.expand_dims(images[0], axis=0)
   1339        file_names = [file_names[0]]
   1340        if include_paths:
   1341            paths = [paths[0]]
   1342
   1343    else:
   1344        bactoscoop_logger.info(f"{len(images)} TIFF files found in folder: {folder_path}")
   1345
   1346    if include_paths:
   1347        return images, file_names, paths
   1348    else:
   1349        return images, file_names
   1350
   1351
   1352# %% Contour and mesh operations
   1353def separate_meshdata(df_nr, cell_nr, data):
   1354    x1 = data[df_nr]["mesh"][cell_nr][:, 0]
   1355    y1 = data[df_nr]["mesh"][cell_nr][:, 1]
   1356    x2 = data[df_nr]["mesh"][cell_nr][:, 2]
   1357    y2 = data[df_nr]["mesh"][cell_nr][:, 3]
   1358    contour = data[df_nr]["contour"][cell_nr]
   1359    return x1, y1, x2, y2, contour
   1360
   1361
   1362def separate_singleframe_meshdata(cell_id, data):
   1363    row = data[data.index == cell_id].iloc[0]
   1364    x1 = row["mesh"][:, 0]
   1365    y1 = row["mesh"][:, 1]
   1366    x2 = row["mesh"][:, 2]
   1367    y2 = row["mesh"][:, 3]
   1368    contour = row["contour"]
   1369    midline = row["midline"]
   1370    return x1, y1, x2, y2, contour, midline
   1371
   1372
   1373def get_cropped_object_data(cnts, mdls, mshs, x, y):
   1374    cropped_cnts = []
   1375    cropped_mdls = []
   1376    cropped_mshs = []
   1377
   1378    for cn, ml, msh in zip(cnts, mdls, mshs):
   1379        crop_cn = cn.copy()
   1380        crop_ml = ml.copy()
   1381        crop_msh = msh.copy()
   1382
   1383        crop_cn -= [x, y]
   1384        if crop_ml.size:
   1385            crop_ml -= [x, y]
   1386        if crop_msh.size:
   1387            crop_msh -= [x, y, x, y]
   1388
   1389        cropped_cnts.append(crop_cn)
   1390        cropped_mdls.append(crop_ml)
   1391        cropped_mshs.append(crop_msh)
   1392
   1393    return cropped_cnts, cropped_mdls, cropped_mshs
   1394
   1395
   1396def get_cropped_cell_data(cnt, mdl, msh, x, y):
   1397    cropped_cn = cnt.copy()
   1398    cropped_ml = mdl.copy()
   1399    cropped_msh = msh.copy()
   1400
   1401    cropped_cn -= [x, y]
   1402    cropped_ml -= [x, y]
   1403    cropped_msh -= [x, y, x, y]
   1404
   1405    return cropped_cn, cropped_ml, cropped_msh
   1406
   1407
   1408def get_uncropped_cell_data(cnt, mdl, msh, x, y):
   1409    uncropped_cn = cnt.copy()
   1410    uncropped_ml = mdl.copy()
   1411    uncropped_msh = msh.copy()
   1412
   1413    uncropped_cn += [x, y]
   1414    uncropped_ml += [x, y]
   1415    uncropped_msh += [x, y, x, y]
   1416
   1417    return uncropped_cn, uncropped_ml, uncropped_msh
   1418
   1419
   1420def mesh2midline(x1, y1, x2, y2):
   1421    x = (x1 + x2) / 2
   1422    y = (y1 + y2) / 2
   1423    line = np.array([x, y]).T
   1424    line = spline_approximation(line, n=len(x1), smooth_factor=3, closed=False)
   1425    return line
   1426
   1427
   1428def get_weighted_intprofile(intensity, width):
   1429
   1430    if intensity.shape[0] == width.shape[0]:
   1431        return intensity * width
   1432
   1433    else:
   1434        return []
   1435
   1436
   1437def get_width(x1, y1, x2, y2):
   1438    width_not_ordered = np.sqrt((x1 - x2) ** 2 + (y1 - y2) ** 2)
   1439    return width_not_ordered
   1440
   1441
   1442def _validate_max_daughter_cell_mesh_rows(value):
   1443    if (
   1444        isinstance(value, (bool, np.bool_))
   1445        or not isinstance(value, (int, np.integer))
   1446        or value < 4
   1447    ):
   1448        raise ValueError("max_daughter_cell_mesh_rows must be an integer >= 4")
   1449    return int(value)
   1450
   1451
   1452def split_mesh2mesh(x1, y1, x2, y2, max_daughter_cell_mesh_rows=800):
   1453    """Rebuild a daughter, rejecting oversized resampling before allocation.
   1454
   1455    The limit applies to straightened rows here and to final daughter rows in
   1456    Image.split_cells. Pole extension can add rows, so the minimum stays there.
   1457    """
   1458    max_daughter_cell_mesh_rows = _validate_max_daughter_cell_mesh_rows(
   1459        max_daughter_cell_mesh_rows
   1460    )
   1461
   1462    length = np.sum(get_step_length_no_px(x1, y1, x2, y2))
   1463    mesh_spacing_px = 0.5
   1464    # Use the same spacing/rounding as straighten_by_orthogonal_lines,
   1465    # before contour reconstruction or its dense resampling allocations.
   1466    predicted_rows = int(round(length / mesh_spacing_px)) if np.isfinite(length) else None
   1467    if predicted_rows is None or predicted_rows > max_daughter_cell_mesh_rows:
   1468        bactoscoop_logger.debug(
   1469            "Skipping daughter before reconstruction: predicted_rows=%s, max_rows=%s",
   1470            predicted_rows, max_daughter_cell_mesh_rows,
   1471        )
   1472        return np.empty((0, 4)), np.empty((0, 2)), np.empty((0, 2))
   1473
   1474    contour = mesh2contour(x1, y1, x2, y2)
   1475    try:
   1476        width, widthno = get_avg_width_no_px(x1, y1, x2, y2)
   1477
   1478        x = (x1 + x2) / 2
   1479        y = (y1 + y2) / 2
   1480        line = np.array([x, y]).T
   1481        midline, pole1, pole2 = extend_skeleton(line[4:-4], contour)
   1482        l1, l2, profile_mesh, midline = straighten_by_orthogonal_lines(
   1483            contour, midline, length, width, unit_micron=mesh_spacing_px
   1484        )
   1485        result, l1, l2 = add_poles(l1, l2, pole1, pole2)
   1486    except Exception as e:
   1487        bactoscoop_logger.debug(f"Exception occurred: {e}")
   1488        pass  # needs logging
   1489    return result, contour, midline
   1490
   1491
   1492# MOMIA
   1493def get_profile_mesh(mesh, width, micron_unit=1):
   1494    l1 = mesh[:, :2]
   1495    l2 = mesh[:, 2:]
   1496    sorted_width = sorted(width, reverse=True)
   1497    width = sum(sorted_width[: math.floor(len(sorted_width) / 3)]) / math.floor(
   1498        len(sorted_width) / 3
   1499    )
   1500    N_width = int(round(np.median(width) / 0.05))
   1501    dl = (l2 - l1) / N_width
   1502    mult_mat = np.tile(np.arange(N_width + 1), (len(l1), 1))
   1503    mat_x = l1[:, 0][:, np.newaxis] + mult_mat * dl[:, 0][:, np.newaxis]
   1504    mat_y = l1[:, 1][:, np.newaxis] + mult_mat * dl[:, 1][:, np.newaxis]
   1505    profile_mesh = np.array([mat_x, mat_y])
   1506    return profile_mesh
   1507
   1508
   1509# MOMIA
   1510def expand_contour(contour, scale=1):
   1511    """
   1512    enlarges contour
   1513    :param contour:
   1514    :param scale:
   1515    :return:
   1516    """
   1517    area = 0.5 * np.sum(np.diff(contour[:, 0]) * (contour[:-1, 1] + contour[1:, 1]))
   1518    if area < 0:
   1519        # If the area is negative, flip the sign of the unit perpendicular vector
   1520        dxy = unit_perpendicular_vector(contour, closed=True)
   1521    else:
   1522        # Otherwise, use the unit perpendicular vector as is
   1523        dxy = unit_perpendicular_vector(contour, closed=True) * (-1)
   1524    # dxy = unit_perpendicular_vector(contour, closed=True)
   1525
   1526    return contour - (scale * dxy)
   1527
   1528
   1529def erode_contour(contour, scale=1):
   1530    """
   1531    shrinks contour
   1532    :param contour:
   1533    :param scale:
   1534    :return:
   1535    """
   1536    area = 0.5 * np.sum(np.diff(contour[:, 0]) * (contour[:-1, 1] + contour[1:, 1]))
   1537    if area < 0:
   1538        # If the area is negative, flip the sign of the unit perpendicular vector
   1539        dxy = unit_perpendicular_vector(contour, closed=True)
   1540    else:
   1541        # Otherwise, use the unit perpendicular vector as is
   1542        dxy = unit_perpendicular_vector(contour, closed=True) * (-1)
   1543    # dxy = unit_perpendicular_vector(contour, closed=True)
   1544
   1545    return contour + (scale * dxy)
   1546
   1547
   1548### ------------------------------------------------------------------------------------------------
   1549### Background subtraction
   1550# calculates the average pixel value of the background region in an image,
   1551# given a binary cell mask.
   1552def mean_background(image, cell_mask):
   1553    # Create a mask for the background region
   1554    background_mask = np.logical_not(cell_mask)
   1555
   1556    # Calculate the average pixel value of the background
   1557    # =>indexing the image array with the background_mask, which selects
   1558    #  only the pixel values
   1559    bgr = np.mean(image[background_mask])
   1560    return bgr
   1561
   1562
   1563def median_background(image, cell_mask):
   1564    # Create a mask for the background region
   1565    background_mask = np.logical_not(cell_mask)
   1566
   1567    # Calculate the average pixel value of the background
   1568    # =>indexing the image array with the background_mask, which selects
   1569    #  only the pixel values
   1570    bgr = np.median(image[background_mask])
   1571    return bgr
   1572
   1573
   1574def poly_func(coords, a, b, c, d, e, f):
   1575    x, y = coords
   1576    return a + b * x + c * y + d * x**2 + e * x * y + f * y**2
   1577
   1578
   1579def construct_bg_function(image, seg_mat, sampling_interval=None, plotF=False):
   1580    """
   1581    Reconstructs the background of an image using 2D polynomial fitting.
   1582
   1583    Parameters:
   1584        image (ndarray): Input grayscale image.
   1585        seg_mat (ndarray): Binary segmentation mask (non-zero = foreground).
   1586        sampling_interval (int): Interval between samples. Default: height // 64.
   1587        plotF (bool): Whether to plot the background fit.
   1588
   1589    Returns:
   1590        bg_function (function): Background function F(x, y).
   1591        bg_matrix (ndarray): Estimated background image.
   1592    """
   1593    H, W = image.shape
   1594    if sampling_interval is None:
   1595        sampling_interval = max(1, H // 64)
   1596
   1597    # Get background pixel coordinates using sparse sampling
   1598    ys, xs = np.mgrid[0:H:sampling_interval, 0:W:sampling_interval]
   1599    xs_flat = xs.ravel()
   1600    ys_flat = ys.ravel()
   1601    mask_flat = seg_mat[ys_flat, xs_flat] == 0
   1602
   1603    x_bg = xs_flat[mask_flat]
   1604    y_bg = ys_flat[mask_flat]
   1605    z_bg = image[y_bg, x_bg]
   1606
   1607    # Fit polynomial model to background samples
   1608    popt, _ = curve_fit(poly_func, (x_bg, y_bg), z_bg, p0=(0, 0, 0, 0, 0, 0))
   1609
   1610    # Vectorized background matrix generation
   1611    xx, yy = np.meshgrid(np.arange(W), np.arange(H))
   1612    bg_matrix = poly_func((xx, yy), *popt)
   1613
   1614    # Optional plotting
   1615    if plotF:
   1616        plot_interval = max(1, sampling_interval)
   1617        xx_sample, yy_sample = np.meshgrid(
   1618            np.arange(0, W, plot_interval), np.arange(0, H, plot_interval)
   1619        )
   1620        zz_sample = poly_func((xx_sample, yy_sample), *popt)
   1621
   1622        fig = plt.figure()
   1623        ax = fig.add_subplot(111, projection="3d")
   1624        ax.plot_surface(xx_sample, yy_sample, zz_sample, alpha=0.5, cmap="viridis")
   1625        ax.scatter(x_bg, y_bg, z_bg, color="r", label="Sampled BG")
   1626        ax.set_xlabel("X")
   1627        ax.set_ylabel("Y")
   1628        plt.legend()
   1629        plt.show()
   1630
   1631    # Return as function and matrix
   1632    bg_function = lambda x, y: poly_func((x, y), *popt)
   1633    return bg_function, bg_matrix
   1634
   1635
   1636###  -------------------------------------------------------------------------------------------------
   1637# %% Morphological features
   1638def get_step_length(x1, y1, x2, y2, px):
   1639    dx = x1[1:] + x2[1:] - x1[:-1] - x2[:-1]
   1640    dy = y1[1:] + y2[1:] - y1[:-1] - y2[:-1]
   1641    return (np.sqrt(dx**2 + dy**2) / 2) * px
   1642
   1643
   1644def get_step_length_no_px(x1, y1, x2, y2):
   1645    dx = x1[1:] + x2[1:] - x1[:-1] - x2[:-1]
   1646    dy = y1[1:] + y2[1:] - y1[:-1] - y2[:-1]
   1647    return np.sqrt(dx**2 + dy**2) / 2
   1648
   1649
   1650def get_length(step_length):
   1651    return np.sum(step_length)
   1652
   1653
   1654def get_avg_width(x1, y1, x2, y2, px):
   1655    width_not_ordered = np.sqrt((x1 - x2) ** 2 + (y1 - y2) ** 2)
   1656    sorted_width = sorted(width_not_ordered, reverse=True)
   1657    width = sum(sorted_width[: math.floor(len(sorted_width) / 3)]) / math.floor(
   1658        len(sorted_width) / 3
   1659    )
   1660    return width * px, width_not_ordered * px
   1661
   1662
   1663def get_avg_width_no_px(x1, y1, x2, y2):
   1664    width_not_ordered = np.sqrt((x1 - x2) ** 2 + (y1 - y2) ** 2)
   1665    sorted_width = sorted(width_not_ordered, reverse=True)
   1666    width = sum(sorted_width[: math.floor(len(sorted_width) / 3)]) / math.floor(
   1667        len(sorted_width) / 3
   1668    )
   1669    return width, width_not_ordered
   1670
   1671
   1672def get_area(contour, px):
   1673    poly = Polygon(contour)
   1674    area = poly.area
   1675    return area * px * px
   1676
   1677
   1678def get_volume(x1, y1, x2, y2, step_length, px):
   1679    d = np.sqrt((x1 - x2) ** 2 + (y1 - y2) ** 2)
   1680    volume = np.trapz((np.pi * (d / 2) ** 2), dx=step_length)
   1681    return volume * px * px
   1682
   1683
   1684def get_surface_area(width_no, step_length):
   1685    widths = width_no[1:]
   1686    surface_areas = 2 * np.pi * (widths / 2) * step_length
   1687    total_surface_area = np.sum(surface_areas)
   1688    return total_surface_area
   1689
   1690
   1691def get_surface_area_over_volume(sa, vol):
   1692    return sa / vol
   1693
   1694
   1695def get_cell_perimeter_measurements(contour, area, px):
   1696    v1 = contour[:-1]
   1697    v2 = contour[1:]
   1698    d = v2 - v1
   1699    perimeter = np.sum(np.sqrt(np.sum(d**2, axis=1))) * px
   1700    circularity = (4 * np.pi * area) / (perimeter) ** 2
   1701    compactness = (perimeter**2) / area
   1702    sphericity = (np.pi * 1.5 * (perimeter / (2 * np.pi)) ** 1.5) / area
   1703    return perimeter, circularity, compactness, sphericity
   1704
   1705
   1706# cell width variability calculated based on the 50% highest cell widths
   1707def get_cell_width_variability(width_no):
   1708    sorted_width = sorted(width_no, reverse=True)
   1709    half_idx = len(sorted_width) // 2
   1710    half_width = sorted_width[:half_idx]
   1711    width_var = np.std(half_width) / np.mean(half_width)
   1712    return width_var
   1713
   1714
   1715def get_curvature_characteristics(contour, px):
   1716    dx = np.gradient(contour[:, 0] * px)
   1717    dy = np.gradient(contour[:, 1] * px)
   1718    d2x = np.gradient(dx)
   1719    d2y = np.gradient(dy)
   1720    curvature = np.abs(dx * d2y - dy * d2x) / (dx**2 + dy**2) ** 1.5
   1721    max_c = np.round(np.max(curvature), 8)
   1722    min_c = np.round(np.min(curvature), 8)
   1723    mean_c = np.round(np.nanmean(curvature), 8)
   1724    std_c = np.round(np.std(curvature), 8)
   1725    return curvature, max_c, min_c, mean_c, std_c
   1726
   1727
   1728def get_total_phaco_intensity(contour, shape, interp2d):
   1729
   1730    mask = polygon2mask(shape, contour)
   1731    coords = np.column_stack(np.where(mask)).astype(int)
   1732    values = interp2d.ev(coords[:, 0], coords[:, 1])
   1733    return np.sum(values), np.max(values), np.mean(values)
   1734    # Calculate the total phase intensity within the contour mask
   1735
   1736
   1737def get_contour_intensity(contour, interp2d):
   1738    data = interp2d.ev(contour.T[0], contour.T[1])
   1739    return data
   1740
   1741
   1742def find_signal_peaks(signal, maximum):
   1743    peaks, _ = find_peaks(signal, prominence=0.5, height=(maximum * 0.5))
   1744    return len(peaks)
   1745
   1746
   1747def measure_contour_variability(signal):
   1748    extra_contour_intensities = np.concatenate([signal, signal[:10]])
   1749    contour_intensities_variability = np.array(
   1750        [
   1751            np.std(extra_contour_intensities[ss : ss + 10])
   1752            for ss in range(1, len(extra_contour_intensities) - 10)
   1753        ]
   1754    )
   1755    return contour_intensities_variability
   1756
   1757
   1758def get_kurtosis(signal):
   1759    return kurtosis(signal)
   1760
   1761
   1762def get_skew(signal):
   1763    return skew(signal)
   1764
   1765
   1766def sinuosity(midline):
   1767    end_to_end_dist = distance(midline[0], midline[-1])
   1768    if end_to_end_dist <= 1e-9:
   1769        raise ValueError("Unable to calculate sinuosity for a closed contour.")
   1770
   1771    length = line_length(midline)
   1772    sinuosity = round(length / end_to_end_dist, 3)
   1773
   1774    if sinuosity < 1:
   1775        sinuosity = 1.0
   1776
   1777    return sinuosity
   1778
   1779
   1780def get_area_asymmetry(step_length, width_no):
   1781    result = [step * width for step, width in zip(step_length, width_no)]
   1782
   1783    # Calculate the two areas
   1784    area1 = sum(result[: len(result) // 2])
   1785    area2 = sum(result[len(result) // 2 :])
   1786
   1787    if area1 == area2:
   1788        return 1.0
   1789    # Determine which area is smaller and which is larger
   1790    smaller_area = min(area1, area2)
   1791    larger_area = max(area1, area2)
   1792
   1793    # Calculate the ratio of the smaller area over the larger area
   1794    asymmetry_ratio = smaller_area / larger_area
   1795
   1796    return asymmetry_ratio
   1797
   1798
   1799#  ------------------------------------------------------------------------------------------------
   1800# %% Object features (multiple contours)
   1801def get_object_params(object_contours, px):
   1802    length_list, width_list, box_list, x_y_list = [], [], [], []
   1803
   1804    for object_contour in object_contours:
   1805        object_contour = object_contour.astype(np.float32)
   1806        rect = cv.minAreaRect(object_contour)
   1807        box = cv.boxPoints(rect).tolist()
   1808        x_y, (width, length), _ = rect
   1809        x_y_list.append(x_y)
   1810        box_list.append(box)
   1811        length_list.append(max(width, length) * px)
   1812        width_list.append(min(width, length) * px)
   1813
   1814    total_length = sum(length_list)
   1815    average_width = sum(width_list) / len(width_list)
   1816
   1817    return total_length, average_width
   1818
   1819
   1820def get_object_width_variability(widths_no):
   1821
   1822    width_var_list = []
   1823
   1824    for width_no in widths_no:
   1825        width_var = get_cell_width_variability(width_no)
   1826        width_var_list.append(width_var)
   1827
   1828    avg_width_var = [np.nan] if not width_var_list else np.mean(width_var_list)
   1829
   1830    return width_var_list, avg_width_var
   1831
   1832
   1833def get_mesh_coordinates(mesh):
   1834    x1 = mesh[:, 0]
   1835    y1 = mesh[:, 1]
   1836    x2 = mesh[:, 2]
   1837    y2 = mesh[:, 3]
   1838    return x1, y1, x2, y2
   1839
   1840
   1841def get_object_rect_length(nucleoid_contours, px):
   1842    rect_length_list = []
   1843    rect_width_list = []
   1844
   1845    for nucleoid_contour in nucleoid_contours:
   1846        nucleoid_contour = nucleoid_contour.astype(np.float32)
   1847
   1848        x_y, width_height, angle_of_rotation = cv.minAreaRect(nucleoid_contour)
   1849
   1850        length = max(width_height) * px
   1851        rect_length_list.append(length)
   1852        width = min(width_height) * px
   1853        rect_width_list.append(width)
   1854
   1855    total_rect_length = sum(rect_length_list)
   1856    # Calculate average width, handling division by zero
   1857    average_rect_width = [np.nan] if not rect_width_list else np.mean(rect_width_list)
   1858
   1859    return rect_length_list, rect_width_list, total_rect_length, average_rect_width
   1860
   1861
   1862def _max_feret_diameter_px(coords):
   1863    """Exact point-set diameter with bounded pairwise-distance memory."""
   1864    if len(coords) < 2:
   1865        return 0.0
   1866    if len(coords) > 2:
   1867        try:
   1868            vertices = coords[spatial.ConvexHull(coords).vertices]
   1869        except spatial.QhullError:
   1870            # A one-pixel-wide, collinear mask has no 2D convex hull.
   1871            axis = int(np.argmax(np.ptp(coords, axis=0)))
   1872            vertices = coords[
   1873                [np.argmin(coords[:, axis]), np.argmax(coords[:, axis])]
   1874            ]
   1875    else:
   1876        vertices = coords
   1877
   1878    maximum = 0.0
   1879    block_size = 256
   1880    for start in range(0, len(vertices), block_size):
   1881        first = vertices[start : start + block_size]
   1882        for second_start in range(start, len(vertices), block_size):
   1883            second = vertices[second_start : second_start + block_size]
   1884            maximum = max(
   1885                maximum,
   1886                float(np.max(spatial.distance.cdist(first, second))),
   1887            )
   1888    return maximum
   1889
   1890
   1891def get_additional_regionprops_features(mask, px):
   1892    mask = mask.astype(bool)
   1893    props = regionprops(mask.astype(int))[0]
   1894    # Calculate Feret diameters
   1895    coords = np.column_stack(np.where(mask))
   1896    max_feret_diameter = _max_feret_diameter_px(coords) * px
   1897
   1898    # Calculate radii
   1899    center = props.centroid
   1900    y, x = np.nonzero(mask)
   1901    distances2 = np.sqrt((x - center[1]) ** 2 + (y - center[0]) ** 2) * px
   1902    max_radius = np.max(distances2)
   1903    mean_radius = np.mean(distances2)
   1904    median_radius = np.median(distances2)
   1905
   1906    # Create features dictionary
   1907    features = {
   1908        "F_MAX_FERET_DIAMETER": max_feret_diameter,
   1909        "F_EQUIVALENT_DIAMETER": props.equivalent_diameter * px,
   1910        "F_MAXIMUM_RADIUS": max_radius,
   1911        "F_MEAN_RADIUS": mean_radius,
   1912        "F_MEDIAN_RADIUS": median_radius,
   1913        "F_EXTENT": props.extent,
   1914    }
   1915
   1916    # Add moments and other properties
   1917    moment_properties = [
   1918        "moments",
   1919        "moments_central",
   1920        "moments_normalized",
   1921        "moments_hu",
   1922        "inertia_tensor",
   1923        "inertia_tensor_eigvals",
   1924    ]
   1925
   1926    for prop in moment_properties:
   1927        if prop in ["moments_hu", "inertia_tensor_eigvals"]:
   1928            for i, value in enumerate(getattr(props, prop)):
   1929                features[f"{prop}_{i}"] = value
   1930        else:
   1931            moments = getattr(props, prop)
   1932            for i in range(moments.shape[0]):
   1933                for j in range(moments.shape[1]):
   1934                    # Skip moments_normalized for indices (0,1), (0,0), and (1,0)
   1935                    if prop == "moments_normalized" and (
   1936                        i == 0 and (j == 1 or j == 0) or i == 1 and j == 0
   1937                    ):
   1938                        continue
   1939                    features[f"{prop}_{i}_{j}"] = moments[i, j]
   1940
   1941    return features
   1942
   1943
   1944def get_object_avg_width(nucleoid_masks, px):
   1945
   1946    nucleoid_avg_width_list = []
   1947    width_not_ordered_list = []
   1948
   1949    for nucleoid_mask in nucleoid_masks:
   1950
   1951        # Convert nucleoid_mask to 8-bit unsigned integer type
   1952        nucleoid_mask = nucleoid_mask.astype(np.uint8)
   1953        # Ensure both arrays have the same data type
   1954        distance = cv.distanceTransform(
   1955            nucleoid_mask, distanceType=cv.DIST_L2, maskSize=5
   1956        ).astype(np.float64)
   1957        # Extract values at skeleton locations
   1958        skeleton = skeletonize(nucleoid_mask)
   1959
   1960        width = distance[skeleton] * 2 * px
   1961        width_not_ordered_list.append(width)
   1962        nucleoid_avg_width_list.append(np.mean(width))
   1963
   1964    return np.mean(nucleoid_avg_width_list), width_not_ordered_list
   1965
   1966
   1967def get_object_aspect_ratio(rect_length_lists, rect_width_lists):
   1968
   1969    aspect_ratio_list = []
   1970    # wrong: has to be length_list/width_list
   1971    for length_list, width_list in zip(rect_length_lists, rect_width_lists):
   1972        aspect_ratio = length_list / width_list
   1973        aspect_ratio_list.append(aspect_ratio)
   1974    avg_aspect_ratio = [np.nan] if not aspect_ratio_list else np.mean(aspect_ratio_list)
   1975
   1976    return aspect_ratio_list, avg_aspect_ratio
   1977
   1978
   1979def get_object_area(contours, px):
   1980
   1981    area_list = []
   1982    for contour in contours:
   1983
   1984        poly = Polygon(contour)
   1985        area = poly.area
   1986        area_list.append(area * px * px)
   1987
   1988    total_area = np.sum(area_list)
   1989    return area_list, total_area
   1990
   1991
   1992def get_object_curvature_characteristics(contours, px):
   1993    curvatures = []
   1994    max_c_list = []
   1995    min_c_list = []
   1996    mean_c_list = []
   1997    std_c_list = []
   1998
   1999    for contour in contours:
   2000        # Calculate gradients of x and y coordinates
   2001        curvature, max_c, min_c, mean_c, std_c = get_curvature_characteristics(
   2002            contour, px
   2003        )
   2004        curvatures.append(curvature)
   2005        max_c_list.append(max_c)
   2006        min_c_list.append(min_c)
   2007        mean_c_list.append(mean_c)
   2008        std_c_list.append(std_c)
   2009
   2010    avg_mean_c = [np.nan] if not mean_c_list else np.mean(mean_c_list)
   2011    avg_std_c = [np.nan] if not std_c_list else np.mean(std_c_list)
   2012
   2013    return (
   2014        curvatures,
   2015        max_c_list,
   2016        min_c_list,
   2017        mean_c_list,
   2018        std_c_list,
   2019        avg_mean_c,
   2020        avg_std_c,
   2021    )
   2022
   2023
   2024def get_object_perimeter_measurements(contours, areas, px):
   2025    perimeter_list = []
   2026    circularity_list = []
   2027    compactness_list = []
   2028    sphericity_list = []
   2029
   2030    for contour, area in zip(contours, areas):
   2031        perimeter, circularity, compactness, sphericity = (
   2032            get_cell_perimeter_measurements(contour, area, px)
   2033        )
   2034
   2035        perimeter_list.append(perimeter)
   2036        circularity_list.append(circularity)
   2037        compactness_list.append(compactness)
   2038        sphericity_list.append(sphericity)
   2039
   2040    avg_perimeter = [np.nan] if not perimeter_list else np.mean(perimeter_list)
   2041    avg_circularity = [np.nan] if not circularity_list else np.mean(circularity_list)
   2042    avg_compactness = [np.nan] if not compactness_list else np.mean(compactness_list)
   2043    avg_sphericity = [np.nan] if not sphericity_list else np.mean(sphericity_list)
   2044
   2045    return (
   2046        perimeter_list,
   2047        circularity_list,
   2048        compactness_list,
   2049        sphericity_list,
   2050        avg_perimeter,
   2051        avg_circularity,
   2052        avg_compactness,
   2053        avg_sphericity,
   2054    )
   2055
   2056
   2057def get_cell_convexity(mask):
   2058    labels = label(mask)
   2059
   2060    properties = regionprops(labels)
   2061    prop = properties[0]
   2062
   2063    image_perimeter = prop.perimeter
   2064    image_eccentricity = prop.eccentricity
   2065    image_solidity = prop.solidity
   2066    convex_image_perimeter = regionprops(prop.convex_image * 1)[0].perimeter
   2067    convexity = round(convex_image_perimeter / image_perimeter, 3)
   2068
   2069    return image_eccentricity, image_solidity, convexity
   2070
   2071
   2072def get_object_convexity(cropped_masks):
   2073
   2074    convexity_list = []
   2075    eccentricity_list = []
   2076    solidity_list = []
   2077
   2078    for nucleoid_mask in cropped_masks:
   2079
   2080        ecc, sol, conv = get_cell_convexity(nucleoid_mask)
   2081
   2082        convexity_list.append(conv)
   2083        eccentricity_list.append(ecc)
   2084        solidity_list.append(sol)
   2085
   2086        avg_convexity = [np.nan] if not convexity_list else np.mean(convexity_list)
   2087        avg_eccentricity = (
   2088            [np.nan] if not eccentricity_list else np.mean(eccentricity_list)
   2089        )
   2090        avg_solidity = [np.nan] if not solidity_list else np.mean(solidity_list)
   2091
   2092    return (
   2093        convexity_list,
   2094        eccentricity_list,
   2095        solidity_list,
   2096        avg_convexity,
   2097        avg_eccentricity,
   2098        avg_solidity,
   2099    )
   2100
   2101
   2102def get_approx_nucleoid_volume(cell_volume, nucleoid_areas, cell_area):
   2103
   2104    nucleoid_volume_list = []
   2105
   2106    for nucleoid_area in nucleoid_areas:
   2107        # The power 3/2 was used to convert the estimated nucleoid area fraction into a volume fraction.
   2108        nucleoid_volume = cell_volume * (nucleoid_area / cell_area) ** (3 / 2)
   2109        # https://doi.org/10.1016/j.cell.2021.05.037
   2110        nucleoid_volume_list.append(nucleoid_volume)
   2111        total_volume = np.sum(nucleoid_volume_list)
   2112
   2113    return nucleoid_volume_list, total_volume
   2114
   2115
   2116def get_nucleoid_width_variability(widths):
   2117    # Because the CV is unitless and usually expressed as a percentage, it is used instead of the SD to compare the spread of data
   2118    # sets that have different units of measurements or have the same units of measurements but differs greatly in magnitude.
   2119    # Ratio between the standard deviation and the mean ( which in some case could be negative)
   2120
   2121    width_var_list = []
   2122    if widths:
   2123
   2124        for width in widths:
   2125            width_var = np.std(width) / np.mean(width)
   2126            width_var_list.append(width_var)
   2127
   2128    avg_width_var = [np.nan] if not width_var_list else np.mean(width_var_list)
   2129
   2130    return width_var_list, avg_width_var
   2131
   2132
   2133def get_obj_mesh_coords(nucleoid_meshes):
   2134    x1_list, y1_list, x2_list, y2_list = [], [], [], []
   2135
   2136    for nucleoid_mesh in nucleoid_meshes:
   2137        if np.any(nucleoid_mesh):
   2138            x1 = nucleoid_mesh[:, 0]
   2139            y1 = nucleoid_mesh[:, 1]
   2140            x2 = nucleoid_mesh[:, 2]
   2141            y2 = nucleoid_mesh[:, 3]
   2142
   2143            x1_list.append(x1)
   2144            y1_list.append(y1)
   2145            x2_list.append(x2)
   2146            y2_list.append(y2)
   2147        else:
   2148            x1_list.append([])
   2149            y1_list.append([])
   2150            x2_list.append([])
   2151            y2_list.append([])
   2152
   2153    return x1_list, y1_list, x2_list, y2_list
   2154
   2155
   2156def measure_along_objects_midline_interp2d(
   2157    midlines, im_interp2d, width=7, subpixel=0.5
   2158):
   2159
   2160    prf_list = []
   2161
   2162    for midline in midlines:
   2163        if np.any(midline):
   2164            # Calculate a unit vector perpendicular to the midline
   2165            unit_dxy = unit_perpendicular_vector(midline, closed=False)
   2166
   2167            # Calculate a normalized vector for the width and subpixel resolution
   2168            width_normalized_dxy = unit_dxy * subpixel
   2169            # It prepares the vectors for use in generating profiles at varying distances from the midline.
   2170
   2171            # Interpolate image values along the midline
   2172            data = im_interp2d.ev(midline.T[0], midline.T[1])
   2173
   2174            # Iterate over width steps and calculate profile values
   2175            for i in range(1, 1 + int(width * 0.5 / subpixel)):
   2176
   2177                dxy = width_normalized_dxy * i
   2178                # This vector represents a displacement in the direction perpendicular to the midline
   2179
   2180                v1 = midline + dxy
   2181                v2 = midline - dxy
   2182                # calculate two points, v1 and v2, located at distances i units away from the midline in opposite directions.
   2183
   2184                p1 = im_interp2d.ev(v1.T[0], v1.T[1])
   2185                ##calculate a value p1 by evaluating the image data at the coordinates v1
   2186
   2187                p2 = im_interp2d.ev(v2.T[0], v2.T[1])
   2188
   2189                data = np.vstack([p1, data, p2])
   2190                # Stack arrays in sequence vertically (row wise).
   2191
   2192            # Calculate the average profile
   2193            prf = np.average(data, axis=0)
   2194            sigma = 2  # standard deviation of Gaussian filter
   2195
   2196            prf = gaussian_filter1d(prf, sigma)
   2197            # A Gaussian Filter is a low pass filter used for reducing noise (high frequency components)
   2198            # and blurring regions of an image.
   2199            prf_list.append(prf)
   2200        else:
   2201            prf_list.append([])
   2202    return prf_list
   2203
   2204
   2205# calculates the constriction degree in relative and absolute value,
   2206# and the relative position along the length of the nucleoid
   2207# relative constr degree needs to divided by the height (hgt) of the local minima
   2208def constr_degree(intensities, new_step_lengths, width_not_ordered):
   2209    constrDegree_list = []
   2210    relPos_list = []
   2211    constrDegree_abs_list = []
   2212    ctpos_list = []
   2213
   2214    for intensity, width_no, new_step_length in zip(
   2215        intensities, width_not_ordered, new_step_lengths
   2216    ):
   2217        if not np.any(intensity):
   2218            fill_nan_values(
   2219                constrDegree_list, constrDegree_abs_list, ctpos_list, relPos_list
   2220            )
   2221        else:
   2222            weighted_intensity = intensity * width_no
   2223            minima = get_minima(weighted_intensity)
   2224            if all(not x for x in minima) or np.sum(weighted_intensity.T) == 0:
   2225                minsize, minsizeabs, ctpos, ctposr = 0, 0, np.nan, np.nan
   2226            else:
   2227                index, dh, dhi, hgt = identify_minima(weighted_intensity, minima)
   2228                fix = np.argmax(dh)
   2229                minsizeabs = dhi[fix]
   2230                minsize = minsizeabs / hgt[fix]
   2231                ctpos = index[fix]
   2232                ctposr = np.cumsum(new_step_length)[ctpos] / np.sum(new_step_length)
   2233
   2234            assign_values(
   2235                constrDegree_list,
   2236                constrDegree_abs_list,
   2237                ctpos_list,
   2238                relPos_list,
   2239                minsize,
   2240                minsizeabs,
   2241                ctpos,
   2242                ctposr,
   2243            )
   2244
   2245    avg_consDegree = [np.nan] if not constrDegree_list else np.mean(constrDegree_list)
   2246    avg_abs_consDegree = (
   2247        [np.nan] if not constrDegree_abs_list else np.mean(constrDegree_abs_list)
   2248    )
   2249    return (
   2250        constrDegree_list,
   2251        relPos_list,
   2252        constrDegree_abs_list,
   2253        ctpos_list,
   2254        avg_consDegree,
   2255        avg_abs_consDegree,
   2256    )
   2257
   2258
   2259def fill_nan_values(constrDegree_list, constrDegree_abs_list, ctpos_list, relPos_list):
   2260    for _ in range(4):
   2261        constrDegree_list.append(np.nan)
   2262    ctpos_list.append(np.nan)
   2263    relPos_list.append(np.nan)
   2264
   2265
   2266def get_minima(weighted_intensity):
   2267    return np.concatenate(
   2268        (
   2269            [False],
   2270            (weighted_intensity.T[1:-1] < weighted_intensity.T[:-2])
   2271            & (weighted_intensity.T[1:-1] < weighted_intensity.T[2:]),
   2272            [False],
   2273        )
   2274    )
   2275
   2276
   2277def identify_minima(intensity, minima):
   2278    index = np.where(minima)[0]
   2279    dh, dhi, hgt = np.zeros(index.shape), np.zeros(index.shape), np.zeros(index.shape)
   2280    for i, k in enumerate(index):
   2281        half1 = intensity.T[: k - 1]
   2282        half2 = intensity.T[k + 1 :]
   2283        try:
   2284            dh1 = np.max(half1) - intensity.T[k]
   2285            dh2 = np.max(half2) - intensity.T[k]
   2286            dh[i] = np.min([dh1, dh2])
   2287            dhi[i] = np.mean([dh1, dh2])
   2288            hgt[i] = intensity.T[k] + dhi[i]
   2289        except ValueError:
   2290            pass
   2291    return index, dh, dhi, hgt
   2292
   2293
   2294def assign_values(
   2295    constrDegree_list,
   2296    constrDegree_abs_list,
   2297    ctpos_list,
   2298    relPos_list,
   2299    minsize,
   2300    minsizeabs,
   2301    ctpos,
   2302    ctposr,
   2303):
   2304    constrDegree_list.append(minsize)
   2305    constrDegree_abs_list.append(minsizeabs)
   2306    ctpos_list.append(ctpos)
   2307    relPos_list.append(ctposr)
   2308
   2309
   2310# calculate the average width of segments defined by coordinates (x1, y1) and (x2, y2)
   2311def get_obj_avg_width(x1s, y1s, x2s, y2s, px):
   2312    width_not_ordered_list = []
   2313    width_list = []
   2314
   2315    for x1, y1, x2, y2 in zip(x1s, y1s, x2s, y2s):
   2316        if np.any(x1):
   2317            width_not_ordered = np.sqrt((x1 - x2) ** 2 + (y1 - y2) ** 2) * px
   2318            width_not_ordered_list.append(width_not_ordered)
   2319
   2320            sorted_width = sorted(width_not_ordered, reverse=True)
   2321            width = sum(sorted_width[: math.floor(len(sorted_width) / 3)]) / math.floor(
   2322                len(sorted_width) / 3
   2323            )
   2324            width_list.append(width)
   2325        else:
   2326            width_not_ordered_list.append(np.nan)
   2327            width_list.append(np.nan)  # Assign np.nan for empty lists
   2328
   2329    average_width = (
   2330        np.nan if not width_list else np.nanmean(width_list)
   2331    )  # Calculate average handling NaNs
   2332    return width_list, width_not_ordered_list, average_width
   2333
   2334
   2335def get_nucleoid_volume(x1s, y1s, x2s, y2s, step_lengths, px):
   2336
   2337    nucleoid_volume_list = []
   2338
   2339    for x1, y1, x2, y2, step_length in zip(x1s, y1s, x2s, y2s, step_lengths):
   2340        if np.any(x1):
   2341            d = np.sqrt((x1 - x2) ** 2 + (y1 - y2) ** 2)
   2342            volume = np.trapz((np.pi * (d / 2) ** 2), dx=step_length)
   2343            nucleoid_volume_list.append(volume * px * px)
   2344
   2345    total_volume = np.nan if not nucleoid_volume_list else np.sum(nucleoid_volume_list)
   2346    return nucleoid_volume_list, total_volume
   2347
   2348
   2349def get_object_step_length(x1s, y1s, x2s, y2s, px):
   2350
   2351    nuc_step_length_list = []
   2352    nuc_length_list = []
   2353    for x1, y1, x2, y2 in zip(x1s, y1s, x2s, y2s):
   2354        if np.any(x1):
   2355
   2356            # calculated as the difference between corresponding x-coordinates of two line segments
   2357            dx = x1[1:] + x2[1:] - x1[:-1] - x2[:-1]
   2358            # x1[1:] => start from index 1
   2359
   2360            dy = y1[1:] + y2[1:] - y1[:-1] - y2[:-1]
   2361
   2362            # Calculate the step length for each segment and multiply by pixel size
   2363            step_length = (np.sqrt(dx**2 + dy**2) / 2) * px
   2364            nuc_step_length_list.append(step_length)
   2365            nuc_length_list.append(np.sum(step_length))
   2366        else:
   2367            nuc_step_length_list.append([])
   2368            nuc_length_list.append([])
   2369    return nuc_step_length_list, nuc_length_list
   2370
   2371
   2372def get_object_surface_area(width_nos, step_lengths):
   2373
   2374    sum_surface_area_list = []
   2375
   2376    for width_no, step_length in zip(width_nos, step_lengths):
   2377        if not np.any(np.isnan(width_no)):
   2378
   2379            widths = width_no[1:]
   2380
   2381            # Calculate the surface areas for each segment
   2382            surface_areas = 2 * np.pi * (widths / 2) * step_length
   2383
   2384            # Calculate the total surface area by summing up individual surface areas
   2385            sum_surface_area = np.sum(surface_areas)
   2386            sum_surface_area_list.append(sum_surface_area)
   2387        else:
   2388            sum_surface_area_list.append(np.nan)
   2389        total_surface_area = (
   2390            [np.nan] if not sum_surface_area_list else np.mean(sum_surface_area_list)
   2391        )
   2392
   2393    return sum_surface_area_list, total_surface_area
   2394
   2395
   2396def get_object_surface_area_over_volume(surface_areas, volumes):
   2397
   2398    ratio_list = []
   2399
   2400    for surface_area, volume in zip(surface_areas, volumes):
   2401
   2402        ratio = surface_area / volume
   2403        ratio_list.append(ratio)
   2404
   2405    avg_ratio = [np.nan] if not ratio_list else np.mean(ratio_list)
   2406
   2407    return ratio_list, avg_ratio
   2408
   2409
   2410def object_bending_energy(data, px):
   2411    """
   2412    Simplification of the formula from https://doi.org/10.1016/j.comgeo.2013.09.001
   2413    """
   2414    bending_energy_list = []
   2415    curvatures, _, _, _, _, _, _ = get_object_curvature_characteristics(data, px)
   2416    for line, curvature in zip(data, curvatures):
   2417        BE = bending_energy(line, curvature, px)
   2418        bending_energy_list.append(np.round(BE, 5))
   2419
   2420    avg_bending_energy = (
   2421        [np.nan] if not bending_energy_list else np.mean(bending_energy_list)
   2422    )
   2423    return bending_energy_list, avg_bending_energy
   2424
   2425
   2426def bending_energy(line, curvature, px):
   2427    length = line_length(line) * px
   2428    sum_squared_curvature = np.sum(curvature**2)
   2429    bending_energy = sum_squared_curvature / length
   2430    return np.round(bending_energy, 5)
   2431
   2432
   2433def midline_midpoint(midline):
   2434    middle_index = len(midline) // 2
   2435    middle_point = midline[middle_index]
   2436    return middle_point
   2437
   2438
   2439def get_object_midpoint(midlines):
   2440    centroid_coordinates = []
   2441    for midline in midlines:
   2442        middle_point = midline_midpoint(midline)
   2443        centroid_coordinates.append(middle_point)
   2444    return centroid_coordinates
   2445
   2446
   2447def get_centroid(cropped_masks):
   2448
   2449    labels = label(cropped_masks)
   2450    props = regionprops_table(labels, properties=("label", "centroid"))
   2451    centroids = [props["centroid-0"], props["centroid-1"]]
   2452
   2453    return [centroids[0][0], centroids[1][0]]
   2454
   2455
   2456def get_object_centroid(cropped_masks):
   2457
   2458    centroid_list = []
   2459
   2460    for nucleoid_mask in cropped_masks:
   2461        centroid = get_centroid(nucleoid_mask)
   2462        centroid_list.append(centroid)
   2463
   2464    return centroid_list
   2465
   2466
   2467def get_interobject_distance(centroid_coordinates, px):
   2468    num_pairs = len(list(combinations(centroid_coordinates, 2)))
   2469
   2470    if num_pairs == 0:
   2471        return np.nan, np.nan  # No pairs found
   2472
   2473    distances = [distance(c1, c2) for c1, c2 in combinations(centroid_coordinates, 2)]
   2474    total_distance = sum(distances) * px
   2475    average_distance = total_distance / num_pairs
   2476
   2477    return average_distance, distances
   2478
   2479
   2480def get_object_sinuosity(midlines):
   2481
   2482    sinuosity_list = []
   2483
   2484    for midline in midlines:
   2485        temp_sinuosity = sinuosity(midline)
   2486        sinuosity_list.append(temp_sinuosity)
   2487    avg_sinuosity = [np.nan] if not sinuosity_list else np.mean(sinuosity_list)
   2488
   2489    return sinuosity_list, avg_sinuosity
   2490
   2491
   2492def get_mesh_area(width_no, step_length):
   2493    result = [step * width for step, width in zip(step_length, width_no)]
   2494    mesh_area = sum(result)
   2495    return mesh_area
   2496
   2497
   2498def ratio_mesh_over_total_area(step_length, width_no, areas):
   2499    ratio_list = []
   2500    for step_length, width_no, area in zip(step_length, width_no, areas):
   2501        mesh_area = get_mesh_area(width_no, step_length)
   2502        ratio = mesh_area / area
   2503        ratio_list.append(ratio)
   2504    return ratio_list
   2505
   2506
   2507def get_shortest_pole_object_distances(pole_coordinates, centroid_coordinates, px):
   2508    if len(centroid_coordinates) == 0:
   2509        return np.nan, np.nan  # No centroid coordinates found
   2510
   2511    distances = [
   2512        distance(pole, c) for pole in pole_coordinates for c in centroid_coordinates
   2513    ]
   2514    average_distance = np.mean(distances) * px if distances else np.nan
   2515
   2516    return average_distance, distances
   2517
   2518
   2519def get_pole_coordinates(mesh):
   2520    pole1 = [mesh[0, 0], mesh[0, 1]]
   2521    pole2 = [mesh[-1, 0], mesh[-1, 1]]
   2522    return pole1, pole2
   2523
   2524
   2525def get_pole_object_distances(mesh, centroid_coordinates, px):
   2526    if len(centroid_coordinates) == 0:
   2527        return np.nan, np.nan  # No centroid coordinates found
   2528    pole_coordinates = get_pole_coordinates(mesh)
   2529    distances = [
   2530        distance(pole, c) for pole in pole_coordinates for c in centroid_coordinates
   2531    ]
   2532    avg_distance = np.mean(distances) * px if distances else np.nan
   2533    # Sort the distances
   2534    sorted_distances = np.sort(distances)
   2535    # Take the two smallest distances
   2536    min_2_distances = sorted_distances[:2] * px
   2537    return avg_distance, min_2_distances, distances
   2538
   2539
   2540def angle_between_vector(v1, v2):
   2541    d1 = np.sqrt(np.sum(v1**2))
   2542    d2 = np.sqrt(np.sum(v2**2))
   2543    return np.arccos(np.dot(v1, v2) / (d1 * d2))
   2544
   2545
   2546def projection_on_midline(centroids, midline, width, px):
   2547
   2548    L = line_length(midline)
   2549    sorted_width = sorted(width, reverse=True)
   2550    avg_width = sum(sorted_width[: math.floor(len(sorted_width) / 3)]) / math.floor(
   2551        len(sorted_width) / 3
   2552    )
   2553    l_list = []
   2554    dl_list = []
   2555    norm_l_list = []
   2556    norm_dl_list = []
   2557    for pp in centroids:
   2558        v_p = np.array(pp)[np.newaxis, :] - midline
   2559        v_p_length = np.sqrt(np.sum(v_p**2, axis=1))
   2560
   2561        nearest_loc = np.argmin(v_p_length)
   2562        if nearest_loc == len(midline) - 1:
   2563            v_m = midline[-1] - midline[-2]
   2564        else:
   2565            v_m = midline[nearest_loc + 1] - midline[nearest_loc]
   2566
   2567        dl = np.dot(v_p[nearest_loc], np.array([-v_m[1], v_m[0]])) / np.linalg.norm(v_m)
   2568
   2569        # Determine the sign of dl based on the orientation of the midline
   2570        if np.cross(v_m, v_p[nearest_loc]) < 0:
   2571            dl = -dl
   2572
   2573        # Normalize signed distance
   2574        norm_dl = dl / avg_width
   2575
   2576        v_x_proj = np.array([1, 0])
   2577        ang1 = angle_between_vector(v_m, v_x_proj)
   2578        ang2 = angle_between_vector(v_p[nearest_loc], v_x_proj)
   2579        ang12 = ang1 - ang2
   2580        dx = np.cos(ang12) * dl
   2581        dy = np.sin(ang12) * dl
   2582
   2583        if nearest_loc == 0:
   2584            l = 0
   2585        else:
   2586            l = line_length(midline[: nearest_loc + 1])
   2587
   2588        w = width[nearest_loc]
   2589
   2590        if w == 0:
   2591            norm_dy = dy / (2 * np.max(width))
   2592            norm_dx = dx / L
   2593        else:
   2594            norm_dx, norm_dy = dx / L, dy / w
   2595        norm_l = l / L
   2596
   2597        l_list.append(np.round(l, 5) * px)
   2598        dl_list.append(np.round(dl, 5) * px)
   2599        norm_l_list.append(np.round(norm_l, 5))
   2600        norm_dl_list.append(np.round(norm_dl, 5))
   2601
   2602    return l_list, dl_list, norm_l_list, norm_dl_list
   2603
   2604
   2605def get_avg_distance_from_center(l_norm_list):
   2606    average_pos = np.mean(l_norm_list)
   2607    avg_dist = average_pos - 0.5
   2608    return avg_dist
   2609
   2610
   2611def get_object_intensities(object_masks, cropped_signal):
   2612    individual_intensity_values = []
   2613
   2614    # Loop through each mask to get the intensity values
   2615    for mask in object_masks:
   2616        intensity_values = cropped_signal[mask != 0]
   2617        individual_intensity_values.append(intensity_values)
   2618
   2619    # Flatten the list of individual intensities to get the combined intensities
   2620    combined_intensity_values = np.concatenate(individual_intensity_values)
   2621
   2622    return combined_intensity_values, individual_intensity_values
   2623
   2624
   2625def get_objects_mean(intensity_values_list):
   2626    means = [np.mean(intensity_values) for intensity_values in intensity_values_list]
   2627    return means
   2628
   2629
   2630def get_objects_std(intensity_values_list):
   2631    stds = [np.std(intensity_values) for intensity_values in intensity_values_list]
   2632    return stds
   2633
   2634
   2635def get_objects_kurtosis(intensity_values_list):
   2636    kurtoses = [
   2637        kurtosis(intensity_values) for intensity_values in intensity_values_list
   2638    ]
   2639    return kurtoses
   2640
   2641
   2642def get_objects_skewness(intensity_values_list):
   2643    skewnesses = [skew(intensity_values) for intensity_values in intensity_values_list]
   2644    return skewnesses
   2645
   2646
   2647def get_variance(data):
   2648    return np.var(data)
   2649
   2650
   2651# %% Profiling
   2652def measure_smoothened_intensity(midline, im_interp2d, width=7, subpixel=0.5):
   2653
   2654    data = measure_intensity_interp2d(
   2655        midline, im_interp2d, width=width, subpixel=subpixel
   2656    )
   2657    prf = gaussian_smoothing(data)
   2658    return prf
   2659
   2660
   2661def gaussian_smoothing(data):
   2662    prf = np.average(data, axis=0)
   2663    sigma = 2  # standard deviation of Gaussian filter
   2664    prf = gaussian_filter1d(prf, sigma)
   2665    return prf.flatten()
   2666
   2667
   2668def measure_intensity_interp2d(midline, im_interp2d, width=7, subpixel=0.5):
   2669
   2670    unit_dxy = unit_perpendicular_vector(midline, closed=False)
   2671    width_normalized_dxy = unit_dxy * subpixel
   2672
   2673    data = im_interp2d.ev(midline.T[0], midline.T[1])
   2674    for i in range(1, 1 + int(width * 0.5 / subpixel)):
   2675        dxy = width_normalized_dxy * i
   2676        v1 = midline + dxy
   2677        v2 = midline - dxy
   2678        p1 = im_interp2d.ev(v1.T[0], v1.T[1])
   2679        p2 = im_interp2d.ev(v2.T[0], v2.T[1])
   2680        data = np.vstack([p1, data, p2])
   2681    return data
   2682
   2683
   2684def normalize_per_cell(profile):
   2685    """
   2686    Normalize the values in a 1D profile using Min-Max scaling.
   2687
   2688    Parameters:
   2689        profile (array-like): A 1D array-like object representing the profile.
   2690
   2691    Returns:
   2692        normalized_profile (ndarray): Normalized profile.
   2693    """
   2694    # Convert the input to a numpy array for ease of manipulation
   2695    profile = np.array(profile)
   2696
   2697    # Calculate the minimum and maximum values in the profile
   2698    min_val = np.min(profile)
   2699    max_val = np.max(profile)
   2700
   2701    # Check if min_val equals max_val to avoid division by zero
   2702    if min_val == max_val:
   2703        return np.zeros_like(
   2704            profile
   2705        )  # Return an array of zeros if all values are the same
   2706
   2707    # Normalize each value to the range [0, 1]
   2708    normalized_profile = (profile - min_val) / (max_val - min_val)
   2709
   2710    return normalized_profile
   2711
   2712
   2713def get_cumulative_step_length_midline(midline, px):
   2714    """
   2715    Calculate cumulative step lengths for each segment in the midline.
   2716
   2717    Parameters:
   2718        midline (ndarray): An array containing X and Y coordinates of the midline.
   2719        px (float): The pixel size.
   2720
   2721    Returns:
   2722        cumulative_step_lengths (ndarray): Array containing cumulative step lengths for each segment in the midline.
   2723    """
   2724    # Extract X and Y coordinates from the midline
   2725    x = midline[:, 0]
   2726    y = midline[:, 1]
   2727
   2728    # Compute the step lengths
   2729    dx = x[1:] - x[:-1]
   2730    dy = y[1:] - y[:-1]
   2731    step_lengths = np.sqrt(dx**2 + dy**2) * px
   2732
   2733    # Calculate cumulative step lengths
   2734    cumulative_step_lengths = np.zeros_like(x)
   2735    cumulative_step_lengths[1:] = np.cumsum(step_lengths)
   2736
   2737    return cumulative_step_lengths
   2738
   2739
   2740def linear_interpolation(array, num_values):
   2741    # Original x-coordinates
   2742    x_original = np.arange(len(array))
   2743    # New x-coordinates for interpolation
   2744    x_new = np.linspace(0, len(array) - 1, num_values)
   2745    # Perform linear interpolation
   2746    interpolated_values = np.interp(x_new, x_original, array)
   2747    return interpolated_values
   2748
   2749
   2750# %% Contour shift functions
   2751def get_phase_mask(phase_cropped_img, phase_log_sigma=0.5):
   2752
   2753    inverted_phase_cropped_image = cv.bitwise_not(phase_cropped_img)
   2754
   2755    log_filtered_img_phase = gaussian_laplace(
   2756        inverted_phase_cropped_image, sigma=phase_log_sigma
   2757    )
   2758    threshold_value = threshold_otsu(log_filtered_img_phase)
   2759    cropped_phase_mask = log_filtered_img_phase > threshold_value
   2760    cropped_phase_mask = cropped_phase_mask.astype(np.uint8)
   2761
   2762    return cropped_phase_mask
   2763
   2764
   2765def closing_mask(cropped_phase_mask, closing_level=2):
   2766
   2767    mask_closed = closing(
   2768        cropped_phase_mask, footprint=np.ones((closing_level, closing_level))
   2769    )
   2770
   2771    return mask_closed
   2772
   2773
   2774def find_centroid_in_largest_mask(phase_mask_closed, phase_cropped_contour):
   2775    # Find the largest connected component in the phase mask
   2776    labeled_mask, num_labels = label(phase_mask_closed)
   2777
   2778    if num_labels > 0:
   2779        label_sizes = ndimage.sum(
   2780            phase_mask_closed, labeled_mask, range(1, num_labels + 1)
   2781        )
   2782        largest_label = np.argmax(label_sizes) + 1  # +1 to adjust for 0-based indexing
   2783        largest_mask = labeled_mask == largest_label
   2784
   2785        # Get the centroid of the largest mask
   2786        phase_centroid = get_centroid(largest_mask)
   2787
   2788        # Check if the centroid lies within the cell contour
   2789        phase_cropped_contour = phase_cropped_contour.astype(np.float32)
   2790        if cv.pointPolygonTest(phase_cropped_contour, phase_centroid, False) >= 0:
   2791            return [phase_centroid]
   2792
   2793    return []
   2794
   2795
   2796def align(signal_centroid, phase_centroid, cropped_contour):
   2797
   2798    # Calculate the translation vector
   2799    translation_vector = np.array(signal_centroid) - np.array(phase_centroid)
   2800
   2801    # Ensure translation_vector is a 1D array
   2802    translation_vector = translation_vector.flatten()
   2803    aligned_contour = cropped_contour.copy()
   2804
   2805    # colume of contour or midline is 2
   2806    if aligned_contour.shape == (len(aligned_contour), 2):
   2807        aligned_contour = aligned_contour + translation_vector
   2808
   2809    # colume of mesh is 4
   2810    if aligned_contour.shape == (len(aligned_contour), 4):
   2811        # Add translation vector to the first and second columns
   2812        aligned_contour[:, :2] += translation_vector
   2813
   2814        # Add translation vector to the third and fourth columns
   2815        aligned_contour[:, 2:] += translation_vector
   2816
   2817    return aligned_contour
   2818
   2819
   2820def align_objects(
   2821    signal_centroid, phase_centroid, cropped_mesh, cropped_contour, cropped_midline
   2822):
   2823    # Aligned contour
   2824    aligned_contour = align(signal_centroid, phase_centroid, cropped_contour)
   2825
   2826    aligned_mesh = align(signal_centroid, phase_centroid, cropped_mesh)
   2827
   2828    aligned_midline = align(signal_centroid, phase_centroid, cropped_midline)
   2829
   2830    return aligned_contour, aligned_midline, aligned_mesh
   2831
   2832
   2833def get_signal_mask(signal, contour, log_sigma, min_overlap_ratio, max_external_ratio):
   2834
   2835    cropped_img, _, cropped_contour, x_offset, y_offset = crop_image(
   2836        image=signal, contour=contour, mask_to_crop=None
   2837    )
   2838    cropped_mask = draw_mask(cropped_contour, cropped_img)
   2839    bgr = mean_background(cropped_img, cropped_mask)
   2840    cropped_img_bgr = np.maximum(cropped_img - bgr, 0)
   2841    # Apply Laplacian of Gaussian (LoG) filter to the cropped_img
   2842    log_filtered_img = gaussian_laplace(cropped_img_bgr, sigma=log_sigma)
   2843
   2844    threshold_value = threshold_otsu(log_filtered_img)
   2845
   2846    object_mask = log_filtered_img < threshold_value
   2847
   2848    object_mask = keep_significant_masks(
   2849        object_mask, cropped_mask, min_overlap_ratio, max_external_ratio
   2850    )
   2851
   2852    return object_mask, cropped_img, x_offset, y_offset
   2853
   2854
   2855def shift_contour(
   2856    phase_image,
   2857    signal_image,
   2858    contour,
   2859    mesh,
   2860    midline,
   2861    log_sigma=1.5,
   2862    kernel_width=3,
   2863    min_overlap_ratio=0.3,
   2864    max_external_ratio=0.8,
   2865    phase_log_sigma=0.5,
   2866    phase_closing_level=2,
   2867    signal_closing_level=12,
   2868    max_shift_correction=10,
   2869):
   2870    """
   2871    Aligns the contour, mesh, and midline based on the centroids of the phase and signal masks.
   2872
   2873    Parameters:
   2874    - phase_image (array): The phase image.
   2875    - signal_image (array): The signal image.
   2876    - contour (array): The contour data.
   2877    - mesh (array): The mesh data.
   2878    - midline (array): The midline data.
   2879    - log_sigma (float): Sigma for the log filter for the signal mask.
   2880    - kernel_width (int): Width of the dilation kernel.
   2881    - min_overlap_ratio (float): Minimum overlap ratio for the signal mask.
   2882    - max_external_ratio (float): Maximum external ratio for the signal mask.
   2883    - phase_log_sigma (float): Sigma for the log filter for the phase mask.
   2884    - phase_closing_level (int): Closing level for the phase mask.
   2885    - signal_closing_level (int): Closing level for the signal mask.
   2886    - max_shift_correction (float): Maximum allowable distance between centroids for correction.
   2887
   2888    Returns:
   2889    - aligned_contour (array): The aligned contour data.
   2890    - aligned_mesh (array): The aligned mesh data.
   2891    - aligned_midline (array): The aligned midline data.
   2892    """
   2893    # Phase image
   2894    phase_cropped_img, _, phase_cropped_contour, x, y = crop_image(phase_image, contour)
   2895
   2896    # FM4-64 image
   2897    signal_mask, cropped_signal, _, _ = get_signal_mask(
   2898        signal_image, contour, log_sigma, min_overlap_ratio, max_external_ratio
   2899    )
   2900
   2901    cropped_contour, cropped_midline, cropped_mesh = get_cropped_cell_data(
   2902        contour, midline, mesh, x, y
   2903    )
   2904
   2905    # Get phase mask
   2906    cropped_phase_mask = get_phase_mask(
   2907        phase_cropped_img, phase_log_sigma=phase_log_sigma
   2908    )
   2909
   2910    # Fill the hole within mask
   2911    phase_mask_closed = closing_mask(
   2912        cropped_phase_mask, closing_level=phase_closing_level
   2913    )
   2914    signal_mask_closed = closing_mask(signal_mask, closing_level=signal_closing_level)
   2915
   2916    kernel = np.ones(
   2917        (kernel_width, kernel_width), np.uint8
   2918    )  # You can adjust the kernel size as needed
   2919
   2920    phase_mask_closed = cv.dilate(
   2921        phase_mask_closed.astype(np.uint8), kernel, iterations=2
   2922    )
   2923    signal_mask_closed = cv.dilate(
   2924        signal_mask_closed.astype(np.uint8), kernel, iterations=2
   2925    )
   2926
   2927    # Find largest mask
   2928    phase_centroid = find_centroid_in_largest_mask(
   2929        phase_mask_closed, phase_cropped_contour
   2930    )
   2931    signal_centroid = get_centroid(signal_mask_closed)
   2932
   2933    # Calculate the distance between the centroids
   2934    distance = np.linalg.norm(np.array(phase_centroid) - np.array(signal_centroid))
   2935
   2936    # Check if the distance exceeds the maximum allowable shift correction
   2937    if distance > max_shift_correction:
   2938        return None, None, None
   2939
   2940    cropped_aligned_contour, cropped_aligned_midline, cropped_aligned_mesh = (
   2941        align_objects(
   2942            signal_centroid,
   2943            phase_centroid,
   2944            cropped_mesh,
   2945            cropped_contour,
   2946            cropped_midline,
   2947        )
   2948    )
   2949    aligned_contour, aligned_midline, aligned_mesh = get_uncropped_cell_data(
   2950        cropped_aligned_contour, cropped_aligned_midline, cropped_aligned_mesh, x, y
   2951    )
   2952
   2953    return aligned_contour, aligned_mesh, aligned_midline
   2954
   2955
   2956def crop_signal_with_optional_shift(image_obj, channel, contour, shift_signal):
   2957    phase_img = image_obj.image if shift_signal else None
   2958    cropped_signal, cropped_mask, cropped_contour, x_offset, y_offset = crop_image(
   2959        image_obj.bg_channels[channel], contour, phase_img=phase_img
   2960    )
   2961    return cropped_signal, cropped_mask, cropped_contour, x_offset, y_offset
   2962
   2963
   2964# %% Membrane features
   2965
   2966
   2967def change_contour_array_length(cell_perimeter, contour, scale=0.02, smooth=1):
   2968    num = int(cell_perimeter / scale)
   2969    adjust_eroded_contour = spline_approximation(
   2970        contour, n=num, smooth_factor=smooth, closed=True
   2971    )
   2972    return adjust_eroded_contour
   2973
   2974def normalized_intensity_by_image(profile, cropped_image):
   2975    # Convert the input to a numpy array for ease of manipulation
   2976    prof = np.array(profile)
   2977
   2978    # Flatten the arrays to avoid issues with multidimensional input
   2979    prof_flat = prof.flatten()
   2980    cropped_image_flat = np.array(cropped_image).flatten()
   2981
   2982    # Calculate the mean of the lowest 5% values in the profile
   2983    min_percentile_value = np.percentile(cropped_image_flat, 5)
   2984    max_percentile_value = np.percentile(prof_flat, 95)
   2985
   2986    # Calculate the mean of the lowest 5% values in the cropped image
   2987    min_intensity = np.mean(cropped_image_flat[cropped_image_flat <= min_percentile_value])
   2988    # Calculate the mean of the highest 5% values in the profile
   2989    max_intensity = np.mean(prof_flat[prof_flat >= max_percentile_value])
   2990
   2991    # Check if min_intensity equals max_intensity to avoid division by zero
   2992    if np.isclose(max_intensity, min_intensity):
   2993        return np.zeros_like(prof)  # Return an array of zeros if all values are the same
   2994
   2995    # Normalize each value to the range [0, 1] using the calculated min and max intensities
   2996    normalized_profile = (prof - min_intensity) / (max_intensity - min_intensity)
   2997
   2998    return normalized_profile
   2999
   3000def get_homogeneity(normalized_intensity):
   3001    # for normalized_int in normalized_intensity:
   3002    gradient = np.gradient(normalized_intensity)
   3003    homogeneity = abs(np.median(abs(gradient)) - np.mean(abs(gradient)))
   3004    return homogeneity
   3005
   3006
   3007def set_bg_zero(contour, cropped_img):
   3008    cell_cropped_mask = draw_mask(contour, cropped_img)
   3009    bg_zero_img = cropped_img * cell_cropped_mask
   3010    return bg_zero_img
   3011
   3012
   3013def rescale_image_255(cropped_image, max_value=255):
   3014
   3015    # Find the minimum and maximum intensity values
   3016    min_intensity = np.min(cropped_image)
   3017    max_intensity = np.max(cropped_image)
   3018
   3019    # Scale the intensity values to the range [0, 255]
   3020    scaled_image = (
   3021        (cropped_image - min_intensity) / (max_intensity - min_intensity) * max_value
   3022    )
   3023
   3024    # Float images are not supported by graycomatrix. Convert the image to an unsigned integer type.
   3025    scaled_image = scaled_image.astype(np.uint16)
   3026
   3027    return scaled_image
   3028
   3029
   3030def create_co_occurrence_matrix(
   3031    cropped_img,
   3032    distances=[1],
   3033    angles=[0, np.pi / 4, np.pi / 2, 3 * np.pi / 4],
   3034    levels=256,
   3035    symmetric=True,
   3036    normed=True,
   3037):
   3038    """
   3039    Create a co-occurrence matrix from a cropped image.
   3040
   3041    Parameters:
   3042    - cropped_image: The cropped grayscale image.
   3043    - distances: List of pixel pair distances.
   3044    - angles: List of angles in radians for pixel pairs.
   3045    - levels: Number of gray levels (typically 256 for 8-bit images).
   3046    - symmetric: If True, the output matrix P[:, :, d, theta] is symmetric.
   3047                 This is accomplished by ignoring the order of value pairs, so both (i, j) and (j, i) are accumulated when (i, j) is encountered for a given offset.
   3048    - normed: If True, normalize each matrix P[:, :, d, theta] by dividing by the total number of accumulated co-occurrences for the given offset.
   3049
   3050    Returns:
   3051    - co_occurrence_matrix: The co-occurrence matrix.
   3052    """
   3053    # Compute the co-occurrence matrix
   3054    co_occurrence_matrix = graycomatrix(
   3055        cropped_img,
   3056        distances=distances,
   3057        angles=angles,
   3058        levels=levels,
   3059        symmetric=True,
   3060        normed=True,
   3061    )
   3062    co_occurrence_matrix[:, :, :, :] *= np.where(
   3063        (np.arange(levels) == 0)[:, None, None, None]
   3064        | (np.arange(levels) == 0)[None, :, None, None],
   3065        0,
   3066        1,
   3067    )
   3068
   3069    return co_occurrence_matrix
   3070
   3071
   3072def get_glcm_feature_means(glcm):
   3073    """
   3074    Compute the mean Haralick properties for a GLCM in one pass.
   3075
   3076    The formulas mirror skimage.feature.graycoprops so callers keep the same
   3077    numerical behavior while avoiding repeated scans of the same matrix.
   3078    """
   3079    if glcm.ndim != 4:
   3080        raise ValueError("GLCM must be a 4D array.")
   3081
   3082    num_level, num_level2, num_dist, num_angle = glcm.shape
   3083    if num_level != num_level2:
   3084        raise ValueError("num_level and num_level2 must be equal.")
   3085    if num_dist <= 0:
   3086        raise ValueError("num_dist must be positive.")
   3087    if num_angle <= 0:
   3088        raise ValueError("num_angle must be positive.")
   3089
   3090    normalized_glcm = glcm.astype(np.float64)
   3091    glcm_sums = np.sum(normalized_glcm, axis=(0, 1), keepdims=True)
   3092    glcm_sums[glcm_sums == 0] = 1
   3093    normalized_glcm /= glcm_sums
   3094
   3095    I = np.arange(num_level, dtype=np.float64).reshape((num_level, 1, 1, 1))
   3096    J = np.arange(num_level, dtype=np.float64).reshape((1, num_level, 1, 1))
   3097    diff = I - J
   3098    squared_diff = diff**2
   3099
   3100    dissimilarity = np.sum(normalized_glcm * np.abs(diff), axis=(0, 1)).mean()
   3101    contrast = np.sum(normalized_glcm * squared_diff, axis=(0, 1)).mean()
   3102    energy = np.sqrt(np.sum(normalized_glcm**2, axis=(0, 1))).mean()
   3103    homogeneity = graycoprops(glcm, "homogeneity").mean()
   3104
   3105    diff_i = I - np.sum(I * normalized_glcm, axis=(0, 1))
   3106    diff_j = J - np.sum(J * normalized_glcm, axis=(0, 1))
   3107    std_i = np.sqrt(np.sum(normalized_glcm * (diff_i**2), axis=(0, 1)))
   3108    std_j = np.sqrt(np.sum(normalized_glcm * (diff_j**2), axis=(0, 1)))
   3109    cov = np.sum(normalized_glcm * (diff_i * diff_j), axis=(0, 1))
   3110
   3111    correlation = np.ones((num_dist, num_angle), dtype=np.float64)
   3112    non_zero_mask = (std_i >= 1e-15) & (std_j >= 1e-15)
   3113    correlation[non_zero_mask] = (
   3114        cov[non_zero_mask] / (std_i[non_zero_mask] * std_j[non_zero_mask])
   3115    )
   3116    entropy = shannon_entropy(glcm, base=2)
   3117
   3118    return (
   3119        dissimilarity,
   3120        correlation.mean(),
   3121        homogeneity,
   3122        energy,
   3123        contrast,
   3124        entropy,
   3125    )
   3126
   3127
   3128def calculate_haralick_features(
   3129    contour,
   3130    cropped_img,
   3131    distances=[2],
   3132    angles=[0, np.pi / 4, np.pi / 2, 3 * np.pi / 4],
   3133    levels=256,
   3134):
   3135    bg_zero_img = set_bg_zero(contour, cropped_img)
   3136    bg_zero_img_rescaled = rescale_image_255(bg_zero_img, max_value=255)
   3137    glcm = create_co_occurrence_matrix(
   3138        bg_zero_img_rescaled,
   3139        distances=distances,
   3140        angles=[0, np.pi / 4, np.pi / 2, 3 * np.pi / 4],
   3141        levels=256,
   3142        symmetric=True,
   3143        normed=True,
   3144    )
   3145    return get_glcm_feature_means(glcm)
   3146
   3147def align_contour_start(eroded_contour,first_point):
   3148    """
   3149    Rotates the eroded_contour so that its point closest to the first point
   3150    of p2p_contour comes first.
   3151
   3152    Parameters
   3153    ----------
   3154    eroded_contour : np.ndarray
   3155        Array of shape (N, 2), representing the contour after erosion.
   3156    p2p_contour : np.ndarray
   3157        Original contour array of shape (N, 2) before erosion.
   3158
   3159    Returns
   3160    -------
   3161    np.ndarray
   3162        Rotated eroded_contour starting from the closest point to the original start.
   3163    """
   3164
   3165    distances = np.linalg.norm(eroded_contour - first_point, axis=1)
   3166    closest_idx = np.argmin(distances)
   3167
   3168    # Rotate so closest point is first
   3169    aligned_contour = np.roll(eroded_contour, -closest_idx, axis=0)
   3170    return aligned_contour
   3171
   3172def calculate_object_haralick_features(
   3173    object_masks,
   3174    cropped_img,
   3175    distances=[2],
   3176    angles=[0, np.pi / 4, np.pi / 2, 3 * np.pi / 4],
   3177    levels=256,
   3178):
   3179    combined_mask = np.sum(object_masks, axis=0)
   3180    bg_zero_img = cropped_img * combined_mask
   3181    bg_zero_img_rescaled = rescale_image_255(bg_zero_img, max_value=255)
   3182    glcm = create_co_occurrence_matrix(
   3183        bg_zero_img_rescaled,
   3184        distances=distances,
   3185        angles=[0, np.pi / 4, np.pi / 2, 3 * np.pi / 4],
   3186        levels=256,
   3187        symmetric=True,
   3188        normed=True,
   3189    )
   3190    return get_glcm_feature_means(glcm)
   3191
   3192
   3193def calculate_contour_haralick_features(
   3194    contour_intensities,
   3195    distances=[2],
   3196    angles=[0, np.pi / 4, np.pi / 2, 3 * np.pi / 4],
   3197    levels=256,
   3198):
   3199    rescaled_contour_intensities = rescale_image_255(contour_intensities, max_value=255)
   3200    glcm = create_co_occurrence_matrix(
   3201        rescaled_contour_intensities,
   3202        distances=distances,
   3203        angles=[0, np.pi / 4, np.pi / 2, 3 * np.pi / 4],
   3204        levels=256,
   3205        symmetric=True,
   3206        normed=True,
   3207    )
   3208    return get_glcm_feature_means(glcm)
   3209
   3210
   3211# Function to split the row's array into two halves
   3212def split_array(row):
   3213    midpoint = len(row) // 2
   3214    first_half = row[:midpoint]
   3215    second_half = row[midpoint:]
   3216    return first_half, second_half
   3217
   3218
   3219def get_complemented_contour_intensity(row):
   3220
   3221    first_half, second_half = split_array(row)
   3222
   3223    # Invert the 'second_half' array
   3224    inverted_second_half = second_half[::-1]
   3225
   3226    # Add corresponding elements from 'first_half' and inverted 'second_half' arrays
   3227    complemented_contour_intensity_profile = np.array(
   3228        [(a + b) / 2 for a, b in zip(first_half, inverted_second_half)]
   3229    )
   3230
   3231    return complemented_contour_intensity_profile
   3232
   3233
   3234# %% Shift functions
   3235def shift_image(img, shift):
   3236    """
   3237    Correct XY drift between phase contrast image and fluorescent image(s).
   3238    :param img: Input image
   3239    :param shift: Subpixel XY drift
   3240    :return: Drift corrected image
   3241    """
   3242    offset_image = fourier_shift(np.fft.fftn(img), shift)
   3243    offset_image = np.fft.ifftn(offset_image)
   3244    offset_image = np.round(offset_image.real)
   3245    offset_image[offset_image <= 0] = 10
   3246
   3247    return offset_image.astype(np.uint16)